From 9507a786f013709b8fc4600198c8dbdda25da94c Mon Sep 17 00:00:00 2001 From: Tim Mattox <timothy.mattox@engilitycorp.com> Date: Thu, 6 Oct 2016 16:01:32 -0400 Subject: [PATCH] USER-DPD: whitespace and indentation fixes --- src/USER-DPD/fix_dpd_energy.cpp | 6 +- src/USER-DPD/pair_dpd_fdt.cpp | 231 +++++++++--------- src/USER-DPD/pair_dpd_fdt_energy.cpp | 342 ++++++++++++++------------- 3 files changed, 291 insertions(+), 288 deletions(-) diff --git a/src/USER-DPD/fix_dpd_energy.cpp b/src/USER-DPD/fix_dpd_energy.cpp index df01b12ff8..e6b7a9e83e 100644 --- a/src/USER-DPD/fix_dpd_energy.cpp +++ b/src/USER-DPD/fix_dpd_energy.cpp @@ -34,10 +34,10 @@ FixDPDenergy::FixDPDenergy(LAMMPS *lmp, int narg, char **arg) : pairDPDE = NULL; pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy",1); - - if(pairDPDE == NULL) + + if (pairDPDE == NULL) error->all(FLERR,"Must use pair_style dpd/fdt/energy with fix dpd/energy"); - if(!(atom->dpd_flag)) + if (!(atom->dpd_flag)) error->all(FLERR,"Must use atom_style dpd/fdt/energy with fix dpd/energy"); } diff --git a/src/USER-DPD/pair_dpd_fdt.cpp b/src/USER-DPD/pair_dpd_fdt.cpp index 0660cc4b64..d151b77c85 100644 --- a/src/USER-DPD/pair_dpd_fdt.cpp +++ b/src/USER-DPD/pair_dpd_fdt.cpp @@ -93,129 +93,130 @@ void PairDPDfdt::compute(int eflag, int vflag) // loop over neighbors of my atoms - if(splitFDT_flag){ - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_dpd = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r = sqrt(rsq); - if (r < EPSILON) continue; // r can be 0.0 in DPD systems - rinv = 1.0/r; - wr = 1.0 - r/cut[itype][jtype]; - wd = wr*wr; - - // conservative force = a0 * wr - fpair = a0[itype][jtype]*wr; - fpair *= factor_dpd*rinv; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + if (splitFDT_flag) { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_dpd = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + if (r < EPSILON) continue; // r can be 0.0 in DPD systems + rinv = 1.0/r; + wr = 1.0 - r/cut[itype][jtype]; + wd = wr*wr; + + // conservative force = a0 * wr + fpair = a0[itype][jtype]*wr; + fpair *= factor_dpd*rinv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + // unshifted eng of conservative term: + // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); + // eng shifted to 0.0 at cutoff + evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd; + evdwl *= factor_dpd; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); } - - if (eflag) { - // unshifted eng of conservative term: - // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); - // eng shifted to 0.0 at cutoff - evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd; - evdwl *= factor_dpd; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); } } - } } else { - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - vxtmp = v[i][0]; - vytmp = v[i][1]; - vztmp = v[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_dpd = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r = sqrt(rsq); - if (r < EPSILON) continue; // r can be 0.0 in DPD systems - rinv = 1.0/r; - delvx = vxtmp - v[j][0]; - delvy = vytmp - v[j][1]; - delvz = vztmp - v[j][2]; - dot = delx*delvx + dely*delvy + delz*delvz; - wr = 1.0 - r/cut[itype][jtype]; - wd = wr*wr; - randnum = random->gaussian(); - gamma_ij = sigma[itype][jtype]*sigma[itype][jtype]/(2.0*force->boltz*temperature); - - // conservative force = a0 * wd - // drag force = -gamma * wd^2 * (delx dot delv) / r - // random force = sigma * wd * rnd * dtinvsqrt; - - fpair = a0[itype][jtype]*wr; - fpair -= gamma_ij*wd*dot*rinv; - fpair += sigma[itype][jtype]*wr*randnum*dtinvsqrt; - fpair *= factor_dpd*rinv; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_dpd = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + if (r < EPSILON) continue; // r can be 0.0 in DPD systems + rinv = 1.0/r; + delvx = vxtmp - v[j][0]; + delvy = vytmp - v[j][1]; + delvz = vztmp - v[j][2]; + dot = delx*delvx + dely*delvy + delz*delvz; + wr = 1.0 - r/cut[itype][jtype]; + wd = wr*wr; + randnum = random->gaussian(); + gamma_ij = sigma[itype][jtype]*sigma[itype][jtype] + / (2.0*force->boltz*temperature); + + // conservative force = a0 * wd + // drag force = -gamma * wd^2 * (delx dot delv) / r + // random force = sigma * wd * rnd * dtinvsqrt; + + fpair = a0[itype][jtype]*wr; + fpair -= gamma_ij*wd*dot*rinv; + fpair += sigma[itype][jtype]*wr*randnum*dtinvsqrt; + fpair *= factor_dpd*rinv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + // unshifted eng of conservative term: + // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); + // eng shifted to 0.0 at cutoff + evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd; + evdwl *= factor_dpd; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); } - - if (eflag) { - // unshifted eng of conservative term: - // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); - // eng shifted to 0.0 at cutoff - evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd; - evdwl *= factor_dpd; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); } } } - } if (vflag_fdotr) virial_fdotr_compute(); } diff --git a/src/USER-DPD/pair_dpd_fdt_energy.cpp b/src/USER-DPD/pair_dpd_fdt_energy.cpp index 544eee4367..902ca00525 100644 --- a/src/USER-DPD/pair_dpd_fdt_energy.cpp +++ b/src/USER-DPD/pair_dpd_fdt_energy.cpp @@ -59,7 +59,7 @@ PairDPDfdtEnergy::~PairDPDfdtEnergy() memory->destroy(a0); memory->destroy(sigma); memory->destroy(kappa); - if(!splitFDT_flag){ + if (!splitFDT_flag) { memory->destroy(duCond); memory->destroy(duMech); } @@ -108,187 +108,189 @@ void PairDPDfdtEnergy::compute(int eflag, int vflag) // loop over neighbors of my atoms - if(splitFDT_flag){ - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_dpd = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r = sqrt(rsq); - if (r < EPSILON) continue; // r can be 0.0 in DPD systems - rinv = 1.0/r; - wr = 1.0 - r/cut[itype][jtype]; - wd = wr*wr; - - // conservative force = a0 * wr - fpair = a0[itype][jtype]*wr; - fpair *= factor_dpd*rinv; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + if (splitFDT_flag) { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_dpd = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + if (r < EPSILON) continue; // r can be 0.0 in DPD systems + rinv = 1.0/r; + wr = 1.0 - r/cut[itype][jtype]; + wd = wr*wr; + + // conservative force = a0 * wr + fpair = a0[itype][jtype]*wr; + fpair *= factor_dpd*rinv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + // unshifted eng of conservative term: + // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); + // eng shifted to 0.0 at cutoff + evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd; + evdwl *= factor_dpd; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); } - - if (eflag) { - // unshifted eng of conservative term: - // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); - // eng shifted to 0.0 at cutoff - evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd; - evdwl *= factor_dpd; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); } } - } } else { - // Allocate memory for duCond and duMech - if (allocated) { - memory->destroy(duCond); - memory->destroy(duMech); - } - memory->create(duCond,nlocal+nghost,"pair:duCond"); - memory->create(duMech,nlocal+nghost,"pair:duMech"); - for (int ii = 0; ii < nlocal+nghost; ii++) { - duCond[ii] = double(0.0); - duMech[ii] = double(0.0); - } + // Allocate memory for duCond and duMech + if (allocated) { + memory->destroy(duCond); + memory->destroy(duMech); + } + memory->create(duCond,nlocal+nghost,"pair:duCond"); + memory->create(duMech,nlocal+nghost,"pair:duMech"); + for (int ii = 0; ii < nlocal+nghost; ii++) { + duCond[ii] = 0.0; + duMech[ii] = 0.0; + } - // loop over neighbors of my atoms - for (int ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - vxtmp = v[i][0]; - vytmp = v[i][1]; - vztmp = v[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_dpd = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r = sqrt(rsq); - if (r < EPSILON) continue; // r can be 0.0 in DPD systems - rinv = 1.0/r; - wr = 1.0 - r/cut[itype][jtype]; - wd = wr*wr; - - delvx = vxtmp - v[j][0]; - delvy = vytmp - v[j][1]; - delvz = vztmp - v[j][2]; - dot = delx*delvx + dely*delvy + delz*delvz; - randnum = random->gaussian(); - - // Compute the current temperature - theta_ij = double(0.5)*(double(1.0)/dpdTheta[i] + double(1.0)/dpdTheta[j]); - theta_ij = double(1.0)/theta_ij; + // loop over neighbors of my atoms + for (int ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_dpd = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + if (r < EPSILON) continue; // r can be 0.0 in DPD systems + rinv = 1.0/r; + wr = 1.0 - r/cut[itype][jtype]; + wd = wr*wr; + + delvx = vxtmp - v[j][0]; + delvy = vytmp - v[j][1]; + delvz = vztmp - v[j][2]; + dot = delx*delvx + dely*delvy + delz*delvz; + randnum = random->gaussian(); + + // Compute the current temperature + theta_ij = 0.5*(1.0/dpdTheta[i] + 1.0/dpdTheta[j]); + theta_ij = 1.0/theta_ij; - gamma_ij = sigma[itype][jtype]*sigma[itype][jtype]/(2.0*force->boltz*theta_ij); - - // conservative force = a0 * wr - // drag force = -gamma * wr^2 * (delx dot delv) / r - // random force = sigma * wr * rnd * dtinvsqrt; - - fpair = a0[itype][jtype]*wr; - fpair -= gamma_ij*wd*dot*rinv; - fpair += sigma[itype][jtype]*wr*randnum*dtinvsqrt; - fpair *= factor_dpd*rinv; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; - } - - if (rmass) { - mass_i = rmass[i]; - mass_j = rmass[j]; - } else { - mass_i = mass[itype]; - mass_j = mass[jtype]; - } - massinv_i = double(1.0) / mass_i; - massinv_j = double(1.0) / mass_j; - - // Compute the mechanical and conductive energy, uMech and uCond - mu_ij = massinv_i + massinv_j; - mu_ij *= force->ftm2v; - - uTmp = gamma_ij*wd*rinv*rinv*dot*dot - double(0.5)*sigma[itype][jtype]*sigma[itype][jtype]*mu_ij*wd; - uTmp -= sigma[itype][jtype]*wr*rinv*dot*randnum*dtinvsqrt; - uTmp *= double(0.5); - - duMech[i] += uTmp; - if (newton_pair || j < nlocal) { - duMech[j] += uTmp; - } + gamma_ij = sigma[itype][jtype]*sigma[itype][jtype] + / (2.0*force->boltz*theta_ij); + + // conservative force = a0 * wr + // drag force = -gamma * wr^2 * (delx dot delv) / r + // random force = sigma * wr * rnd * dtinvsqrt; + + fpair = a0[itype][jtype]*wr; + fpair -= gamma_ij*wd*dot*rinv; + fpair += sigma[itype][jtype]*wr*randnum*dtinvsqrt; + fpair *= factor_dpd*rinv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (rmass) { + mass_i = rmass[i]; + mass_j = rmass[j]; + } else { + mass_i = mass[itype]; + mass_j = mass[jtype]; + } + massinv_i = 1.0 / mass_i; + massinv_j = 1.0 / mass_j; + + // Compute the mechanical and conductive energy, uMech and uCond + mu_ij = massinv_i + massinv_j; + mu_ij *= force->ftm2v; + + uTmp = gamma_ij*wd*rinv*rinv*dot*dot + - 0.5*sigma[itype][jtype]*sigma[itype][jtype]*mu_ij*wd; + uTmp -= sigma[itype][jtype]*wr*rinv*dot*randnum*dtinvsqrt; + uTmp *= 0.5; + + duMech[i] += uTmp; + if (newton_pair || j < nlocal) { + duMech[j] += uTmp; + } - // Compute uCond - randnum = random->gaussian(); - kappa_ij = kappa[itype][jtype]; - alpha_ij = sqrt(2.0*force->boltz*kappa_ij); - randPair = alpha_ij*wr*randnum*dtinvsqrt; - - uTmp = kappa_ij*(double(1.0)/dpdTheta[i] - double(1.0)/dpdTheta[j])*wd; - uTmp += randPair; + // Compute uCond + randnum = random->gaussian(); + kappa_ij = kappa[itype][jtype]; + alpha_ij = sqrt(2.0*force->boltz*kappa_ij); + randPair = alpha_ij*wr*randnum*dtinvsqrt; + + uTmp = kappa_ij*(1.0/dpdTheta[i] - 1.0/dpdTheta[j])*wd; + uTmp += randPair; - duCond[i] += uTmp; - if (newton_pair || j < nlocal) { - duCond[j] -= uTmp; - } + duCond[i] += uTmp; + if (newton_pair || j < nlocal) { + duCond[j] -= uTmp; + } - if (eflag) { - // unshifted eng of conservative term: - // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); - // eng shifted to 0.0 at cutoff - evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd; - evdwl *= factor_dpd; + if (eflag) { + // unshifted eng of conservative term: + // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); + // eng shifted to 0.0 at cutoff + evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd; + evdwl *= factor_dpd; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); } } - } - // Communicate the ghost delta energies to the locally owned atoms - comm->reverse_comm_pair(this); + // Communicate the ghost delta energies to the locally owned atoms + comm->reverse_comm_pair(this); } if (vflag_fdotr) virial_fdotr_compute(); @@ -316,7 +318,7 @@ void PairDPDfdtEnergy::allocate() memory->create(a0,n+1,n+1,"pair:a0"); memory->create(sigma,n+1,n+1,"pair:sigma"); memory->create(kappa,n+1,n+1,"pair:kappa"); - if(!splitFDT_flag){ + if (!splitFDT_flag) { memory->create(duCond,nlocal+nghost+1,"pair:duCond"); memory->create(duMech,nlocal+nghost+1,"pair:duMech"); } -- GitLab