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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8795 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2012-09-18 22:55:16 +00:00
parent 6367c8b601
commit 1bc399ec36
5 changed files with 364 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files

65
src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
View File

@ -13,6 +13,7 @@
/* ----------------------------------------------------------------------
Contributing authors: Amalie Frischknecht and Ahmed Ismail (SNL)
simpler force assignment added by Rolf Isele-Holder (Aachen University)
------------------------------------------------------------------------- */
#include "math.h"
@ -75,7 +76,7 @@ PairLJCutCoulLongTIP4P::~PairLJCutCoulLongTIP4P()
void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype,itable;
int i,j,ii,jj,inum,jnum,itype,jtype,itable,key;
int n,vlist[6];
int iH1,iH2,jH1,jH2;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
@ -277,6 +278,7 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
// vlist stores 2,4,6 atoms whose forces contribute to virial
n = 0;
key = 0;
if (itype != typeO) {
f[i][0] += delx * cforce;
@ -290,33 +292,23 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
v[3] = x[i][0] * dely * cforce;
v[4] = x[i][0] * delz * cforce;
v[5] = x[i][1] * delz * cforce;
vlist[n++] = i;
}
vlist[n++] = i;
} else {
key++;
fd[0] = delx*cforce;
fd[1] = dely*cforce;
fd[2] = delz*cforce;
delxOM = x[i][0] - x1[0];
delyOM = x[i][1] - x1[1];
delzOM = x[i][2] - x1[2];
fO[0] = fd[0]*(1 - alpha);
fO[1] = fd[1]*(1 - alpha);
fO[2] = fd[2]*(1 - alpha);
ddotf = (delxOM * fd[0] + delyOM * fd[1] + delzOM * fd[2]) /
(qdist*qdist);
f1[0] = ddotf * delxOM;
f1[1] = ddotf * delyOM;
f1[2] = ddotf * delzOM;
fO[0] = fd[0] - alpha * (fd[0] - f1[0]);
fO[1] = fd[1] - alpha * (fd[1] - f1[1]);
fO[2] = fd[2] - alpha * (fd[2] - f1[2]);
fH[0] = 0.5 * alpha * (fd[0] - f1[0]);
fH[1] = 0.5 * alpha * (fd[1] - f1[1]);
fH[2] = 0.5 * alpha * (fd[2] - f1[2]);
fH[0] = 0.5 * alpha * fd[0];
fH[1] = 0.5 * alpha * fd[1];
fH[2] = 0.5 * alpha * fd[2];
f[i][0] += fO[0];
f[i][1] += fO[1];
@ -340,12 +332,11 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
}
vlist[n++] = i;
vlist[n++] = iH1;
vlist[n++] = iH2;
}
}
if (jtype != typeO) {
f[j][0] -= delx * cforce;
@ -359,33 +350,23 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
v[3] -= x[j][0] * dely * cforce;
v[4] -= x[j][0] * delz * cforce;
v[5] -= x[j][1] * delz * cforce;
vlist[n++] = j;
}
vlist[n++] = j;
} else {
key += 2;
fd[0] = -delx*cforce;
fd[1] = -dely*cforce;
fd[2] = -delz*cforce;
delxOM = x[j][0] - x2[0];
delyOM = x[j][1] - x2[1];
delzOM = x[j][2] - x2[2];
fO[0] = fd[0]*(1 - alpha);
fO[1] = fd[1]*(1 - alpha);
fO[2] = fd[2]*(1 - alpha);
ddotf = (delxOM * fd[0] + delyOM * fd[1] + delzOM * fd[2]) /
(qdist*qdist);
f1[0] = ddotf * delxOM;
f1[1] = ddotf * delyOM;
f1[2] = ddotf * delzOM;
fO[0] = fd[0] - alpha * (fd[0] - f1[0]);
fO[1] = fd[1] - alpha * (fd[1] - f1[1]);
fO[2] = fd[2] - alpha * (fd[2] - f1[2]);
fH[0] = 0.5 * alpha * (fd[0] - f1[0]);
fH[1] = 0.5 * alpha * (fd[1] - f1[1]);
fH[2] = 0.5 * alpha * (fd[2] - f1[2]);
fH[0] = 0.5 * alpha * fd[0];
fH[1] = 0.5 * alpha * fd[1];
fH[2] = 0.5 * alpha * fd[2];
f[j][0] += fO[0];
f[j][1] += fO[1];
@ -409,12 +390,11 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
}
vlist[n++] = j;
vlist[n++] = jH1;
vlist[n++] = jH2;
}
}
if (eflag) {
if (!ncoultablebits || rsq <= tabinnersq)
@ -426,7 +406,7 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
if (evflag) ev_tally_list(n,vlist,ecoul,v);
if (evflag) ev_tally_list(ecoul,vlist,v,alpha,key);
}
}
}
@ -614,3 +594,4 @@ double PairLJCutCoulLongTIP4P::memory_usage()
bytes += 2 * nmax * sizeof(double);
return bytes;
}

268
src/KSPACE/pppm_tip4p.cpp
View File

@ -22,8 +22,10 @@
#include "force.h"
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define OFFSET 16384
@ -38,7 +40,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PPPMTIP4P::PPPMTIP4P(LAMMPS *lmp, int narg, char **arg) :
PPPMOld(lmp, narg, arg) {}
PPPM(lmp, narg, arg) {}
/* ---------------------------------------------------------------------- */
@ -49,7 +51,7 @@ void PPPMTIP4P::init()
if (force->newton == 0)
error->all(FLERR,"Kspace style pppm/tip4p requires newton on");
PPPMOld::init();
PPPM::init();
}
/* ----------------------------------------------------------------------
@ -162,6 +164,16 @@ void PPPMTIP4P::make_rho()
------------------------------------------------------------------------- */
void PPPMTIP4P::fieldforce()
{
if (differentiation_flag == 1) fieldforce_ad();
else fieldforce_ik();
}
/* ----------------------------------------------------------------------
interpolate from grid to get electric field & force on my particles for ik
------------------------------------------------------------------------- */
void PPPMTIP4P::fieldforce_ik()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz,x0,y0,z0;
@ -170,14 +182,12 @@ void PPPMTIP4P::fieldforce()
int iH1,iH2;
double xM[3];
double fx,fy,fz;
double ddotf, rOMx, rOMy, rOMz, f1x, f1y, f1z;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
double *q = atom->q;
double **x = atom->x;
double **f = atom->f;
@ -231,27 +241,245 @@ void PPPMTIP4P::fieldforce()
fz = qfactor * ekz;
find_M(i,iH1,iH2,xM);
rOMx = xM[0] - x[i][0];
rOMy = xM[1] - x[i][1];
rOMz = xM[2] - x[i][2];
f[i][0] += fx*(1 - alpha);
f[i][1] += fy*(1 - alpha);
f[i][2] += fz*(1 - alpha);
ddotf = (rOMx * fx + rOMy * fy + rOMz * fz) / (qdist * qdist);
f[iH1][0] += 0.5*alpha*fx;
f[iH1][1] += 0.5*alpha*fy;
f[iH1][2] += 0.5*alpha*fz;
f1x = ddotf * rOMx;
f1y = ddotf * rOMy;
f1z = ddotf * rOMz;
f[iH2][0] += 0.5*alpha*fx;
f[iH2][1] += 0.5*alpha*fy;
f[iH2][2] += 0.5*alpha*fz;
}
}
}
f[i][0] += fx - alpha * (fx - f1x);
f[i][1] += fy - alpha * (fy - f1y);
f[i][2] += fz - alpha * (fz - f1z);
/* ----------------------------------------------------------------------
interpolate from grid to get electric field & force on my particles for ad
------------------------------------------------------------------------- */
f[iH1][0] += 0.5*alpha*(fx - f1x);
f[iH1][1] += 0.5*alpha*(fy - f1y);
f[iH1][2] += 0.5*alpha*(fz - f1z);
void PPPMTIP4P::fieldforce_ad()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz,x0,y0,z0,dx0,dy0,dz0;
FFT_SCALAR ekx,eky,ekz;
double *xi;
int iH1,iH2;
double xM[3];
double s1,s2,s3;
double sf;
double *prd;
double fx,fy,fz;
f[iH2][0] += 0.5*alpha*(fx - f1x);
f[iH2][1] += 0.5*alpha*(fy - f1y);
f[iH2][2] += 0.5*alpha*(fz - f1z);
if (triclinic == 0) prd = domain->prd;
else prd = domain->prd_lamda;
double xprd = prd[0];
double yprd = prd[1];
double zprd = prd[2];
double zprd_slab = zprd*slab_volfactor;
double hx_inv = nx_pppm/xprd;
double hy_inv = ny_pppm/yprd;
double hz_inv = nz_pppm/zprd;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
double *q = atom->q;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (type[i] == typeO) {
find_M(i,iH1,iH2,xM);
xi = xM;
} else xi = x[i];
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
compute_rho1d(dx,dy,dz);
compute_drho1d(dx,dy,dz);
ekx = eky = ekz = ZEROF;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
for (m = nlower; m <= nupper; m++) {
my = m+ny;
for (l = nlower; l <= nupper; l++) {
mx = l+nx;
ekx += drho1d[0][l]*rho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
eky += rho1d[0][l]*drho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
ekz += rho1d[0][l]*rho1d[1][m]*drho1d[2][n]*u_brick[mz][my][mx];
}
}
}
ekx *= hx_inv;
eky *= hy_inv;
ekz *= hz_inv;
// convert E-field to force and substract self forces
const double qfactor = force->qqrd2e * scale;
s1 = x[i][0]*hx_inv;
s2 = x[i][1]*hy_inv;
s3 = x[i][2]*hz_inv;
sf = sf_coeff[0]*sin(2*MY_PI*s1);
sf += sf_coeff[1]*sin(4*MY_PI*s1);
sf *= 2*q[i]*q[i];
fx = qfactor*(ekx*q[i] - sf);
sf = sf_coeff[2]*sin(2*MY_PI*s2);
sf += sf_coeff[3]*sin(4*MY_PI*s2);
sf *= 2*q[i]*q[i];
fy = qfactor*(eky*q[i] - sf);
sf = sf_coeff[4]*sin(2*MY_PI*s3);
sf += sf_coeff[5]*sin(4*MY_PI*s3);
sf *= 2*q[i]*q[i];
fz = qfactor*(ekz*q[i] - sf);
if (type[i] != typeO) {
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
} else {
find_M(i,iH1,iH2,xM);
f[i][0] += fx*(1 - alpha);
f[i][1] += fy*(1 - alpha);
f[i][2] += fz*(1 - alpha);
f[iH1][0] += 0.5*alpha*fx;
f[iH1][1] += 0.5*alpha*fy;
f[iH1][2] += 0.5*alpha*fz;
f[iH2][0] += 0.5*alpha*fx;
f[iH2][1] += 0.5*alpha*fy;
f[iH2][2] += 0.5*alpha*fz;
}
}
}
/* ----------------------------------------------------------------------
interpolate from grid to get electric field & force on my particles
------------------------------------------------------------------------- */
void PPPMTIP4P::fieldforce_peratom()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz,x0,y0,z0;
double *xi;
int iH1,iH2;
double xM[3];
FFT_SCALAR u_pa,v0,v1,v2,v3,v4,v5;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
double *q = atom->q;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (type[i] == typeO) {
find_M(i,iH1,iH2,xM);
xi = xM;
} else xi = x[i];
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
dx = nx+shiftone - (xi[0]-boxlo[0])*delxinv;
dy = ny+shiftone - (xi[1]-boxlo[1])*delyinv;
dz = nz+shiftone - (xi[2]-boxlo[2])*delzinv;
compute_rho1d(dx,dy,dz);
u_pa = v0 = v1 = v2 = v3 = v4 = v5 = ZEROF;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
z0 = rho1d[2][n];
for (m = nlower; m <= nupper; m++) {
my = m+ny;
y0 = z0*rho1d[1][m];
for (l = nlower; l <= nupper; l++) {
mx = l+nx;
x0 = y0*rho1d[0][l];
if (eflag_atom) u_pa += x0*u_brick[mz][my][mx];
if (vflag_atom) {
v0 += x0*v0_brick[mz][my][mx];
v1 += x0*v1_brick[mz][my][mx];
v2 += x0*v2_brick[mz][my][mx];
v3 += x0*v3_brick[mz][my][mx];
v4 += x0*v4_brick[mz][my][mx];
v5 += x0*v5_brick[mz][my][mx];
}
}
}
}
if (eflag_atom) {
if (type[i] != typeO) {
eatom[i] += q[i]*u_pa;
} else {
eatom[i] += q[i]*u_pa*(1-alpha);
eatom[iH1] += q[i]*u_pa*alpha*0.5;
eatom[iH2] += q[i]*u_pa*alpha*0.5;
}
}
if (vflag_atom) {
if (type[i] != typeO) {
vatom[i][0] += v0*q[i];
vatom[i][1] += v1*q[i];
vatom[i][2] += v2*q[i];
vatom[i][3] += v3*q[i];
vatom[i][4] += v4*q[i];
vatom[i][5] += v5*q[i];
} else {
vatom[i][0] += v0*(1-alpha)*q[i];
vatom[i][1] += v1*(1-alpha)*q[i];
vatom[i][2] += v2*(1-alpha)*q[i];
vatom[i][3] += v3*(1-alpha)*q[i];
vatom[i][4] += v4*(1-alpha)*q[i];
vatom[i][5] += v5*(1-alpha)*q[i];
vatom[iH1][0] += v0*alpha*0.5*q[i];
vatom[iH1][1] += v1*alpha*0.5*q[i];
vatom[iH1][2] += v2*alpha*0.5*q[i];
vatom[iH1][3] += v3*alpha*0.5*q[i];
vatom[iH1][4] += v4*alpha*0.5*q[i];
vatom[iH1][5] += v5*alpha*0.5*q[i];
vatom[iH2][0] += v0*alpha*0.5*q[i];
vatom[iH2][1] += v1*alpha*0.5*q[i];
vatom[iH2][2] += v2*alpha*0.5*q[i];
vatom[iH2][3] += v3*alpha*0.5*q[i];
vatom[iH2][4] += v4*alpha*0.5*q[i];
vatom[iH2][5] += v5*alpha*0.5*q[i];
}
}
}
}

7
src/KSPACE/pppm_tip4p.h
View File

@ -20,11 +20,11 @@ KSpaceStyle(pppm/tip4p,PPPMTIP4P)
#ifndef LMP_PPPM_TIP4P_H
#define LMP_PPPM_TIP4P_H
#include "pppm_old.h"
#include "pppm.h"
namespace LAMMPS_NS {
class PPPMTIP4P : public PPPMOld {
class PPPMTIP4P : public PPPM {
public:
PPPMTIP4P(class LAMMPS *, int, char **);
virtual ~PPPMTIP4P () {};
@ -34,6 +34,9 @@ class PPPMTIP4P : public PPPMOld {
virtual void particle_map();
virtual void make_rho();
virtual void fieldforce();
void fieldforce_ik();
void fieldforce_ad();
virtual void fieldforce_peratom();
private:
void find_M(int, int &, int &, double *);

83
src/pair.cpp
View File

@ -755,6 +755,8 @@ void Pair::ev_tally4(int i, int j, int k, int m, double evdwl,
}
/* ----------------------------------------------------------------------
NOTE: this is older version of ev tally for TIP4P
NOTE: delete this when all pair TIP4P derivatives have been updated
tally ecoul and virial into each of n atoms in list
called by TIP4P potential, newton_pair is always on
changes v values by dividing by n
@ -802,6 +804,87 @@ void Pair::ev_tally_list(int n, int *list, double ecoul, double *v)
}
}
/* ----------------------------------------------------------------------
tally ecoul and virial into each of atoms in list, list length set by key
called by TIP4P potential, newton_pair is always on
------------------------------------------------------------------------- */
void Pair::ev_tally_list(double ecoul, int *list, double *v,
double alpha, int key)
{
int i,j;
if (eflag_either) {
if (eflag_global) eng_coul += ecoul;
if (eflag_atom) {
if (key == 0) {
eatom[list[0]] += 0.5*ecoul;
eatom[list[1]] += 0.5*ecoul;
} else if (key == 1) {
eatom[list[0]] += 0.5*ecoul*(1-alpha);
eatom[list[1]] += 0.25*ecoul*alpha;
eatom[list[2]] += 0.25*ecoul*alpha;
eatom[list[3]] += 0.5*ecoul;
} else if (key == 2) {
eatom[list[0]] += 0.5*ecoul;
eatom[list[1]] += 0.5*ecoul*(1-alpha);
eatom[list[2]] += 0.25*ecoul*alpha;
eatom[list[3]] += 0.25*ecoul*alpha;
} else {
eatom[list[0]] += 0.5*ecoul*(1-alpha);
eatom[list[1]] += 0.25*ecoul*alpha;
eatom[list[2]] += 0.25*ecoul*alpha;
eatom[list[3]] += 0.5*ecoul*(1-alpha);
eatom[list[4]] += 0.25*ecoul*alpha;
eatom[list[5]] += 0.25*ecoul*alpha;
}
}
}
if (vflag_either) {
if (vflag_global) {
virial[0] += v[0];
virial[1] += v[1];
virial[2] += v[2];
virial[3] += v[3];
virial[4] += v[4];
virial[5] += v[5];
}
if (vflag_atom) {
if (key == 0) {
for (i = 0; i <= 5; i++) {
vatom[list[0]][i] += 0.5*v[i];
vatom[list[1]][i] += 0.5*v[i];
}
} else if (key == 1) {
for (i = 0; i <= 5; i++) {
vatom[list[0]][i] += 0.5*v[i]*(1-alpha);
vatom[list[1]][i] += 0.25*v[i]*alpha;
vatom[list[2]][i] += 0.25*v[i]*alpha;
vatom[list[3]][i] += 0.5*v[i];
}
} else if (key == 2) {
for (i = 0; i <= 5; i++) {
vatom[list[0]][i] += 0.5*v[i];
vatom[list[1]][i] += 0.5*v[i]*(1-alpha);
vatom[list[2]][i] += 0.25*v[i]*alpha;
vatom[list[3]][i] += 0.25*v[i]*alpha;
}
} else {
for (i = 0; i <= 5; i++) {
vatom[list[0]][i] += 0.5*v[i]*(1-alpha);
vatom[list[1]][i] += 0.25*v[i]*alpha;
vatom[list[2]][i] += 0.25*v[i]*alpha;
vatom[list[3]][i] += 0.5*v[i]*(1-alpha);
vatom[list[4]][i] += 0.25*v[i]*alpha;
vatom[list[5]][i] += 0.25*v[i]*alpha;
}
}
}
}
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by REAX/C potential, newton_pair is always on

1
src/pair.h
View File

@ -173,6 +173,7 @@ class Pair : protected Pointers {
void ev_tally4(int, int, int, int, double,
double *, double *, double *, double *, double *, double *);
void ev_tally_list(int, int *, double, double *);
void ev_tally_list(double, int *, double *, double, int);
void v_tally2(int, int, double, double *);
void v_tally_tensor(int, int, int, int,
double, double, double, double, double, double);