ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

USER-DPD: whitespace and indentation fixes

Browse Source
This commit is contained in:
Tim Mattox
2016-10-06 16:01:32 -04:00
parent 9789f047d7
commit 9507a786f0
3 changed files with 270 additions and 267 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/USER-DPD/fix_dpd_energy.cpp
View File

@ -35,9 +35,9 @@ 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");
}

209
src/USER-DPD/pair_dpd_fdt.cpp
View File

@ -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];
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;
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];
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;
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;
// 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;
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 (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;
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];
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);
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;
// 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;
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;
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();
}

322
src/USER-DPD/pair_dpd_fdt_energy.cpp
View File

@ -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];
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;
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];
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;
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;
// 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;
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];
// 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;
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];
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;
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();
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;
// 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);
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;
// 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;
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;
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*(1.0/dpdTheta[i] - 1.0/dpdTheta[j])*wd;
uTmp += randPair;
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 (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
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;
}
// 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;
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 (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");
}