diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 2dcea05a06..5bacae48be 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -700,10 +700,13 @@ void PPPM::compute(int eflag, int vflag) // per-atom energy/virial // energy includes self-energy correction + // notal accounts for TIP4P tallying eatom/vatom for ghost atoms if (evflag_atom) { double *q = atom->q; int nlocal = atom->nlocal; + int ntotal = nlocal; + if (tip4pflag) ntotal += atom->nghost; if (eflag_atom) { for (i = 0; i < nlocal; i++) { @@ -712,10 +715,11 @@ void PPPM::compute(int eflag, int vflag) (g_ewald*g_ewald*volume); eatom[i] *= qscale; } + for (i = nlocal; i < ntotal; i++) eatom[i] *= 0.5*qscale; } if (vflag_atom) { - for (i = 0; i < nlocal; i++) + for (i = 0; i < ntotal; i++) for (j = 0; j < 6; j++) vatom[i][j] *= 0.5*qscale; } } diff --git a/src/KSPACE/pppm_disp_tip4p.cpp b/src/KSPACE/pppm_disp_tip4p.cpp index 1de82c867e..72d2020e00 100755 --- a/src/KSPACE/pppm_disp_tip4p.cpp +++ b/src/KSPACE/pppm_disp_tip4p.cpp @@ -1,518 +1,519 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Amalie Frischknecht and Ahmed Ismail (SNL) - Rolf Isele-Holder (Aachen University) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "pppm_disp_tip4p.h" -#include "pppm_disp.h" -#include "atom.h" -#include "domain.h" -#include "force.h" -#include "memory.h" -#include "error.h" -#include "math_const.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define OFFSET 16384 - -#ifdef FFT_SINGLE -#define ZEROF 0.0f -#define ONEF 1.0f -#else -#define ZEROF 0.0 -#define ONEF 1.0 -#endif - -/* ---------------------------------------------------------------------- */ - -PPPMDispTIP4P::PPPMDispTIP4P(LAMMPS *lmp, int narg, char **arg) : - PPPMDisp(lmp, narg, arg) -{ - triclinic_support = 0; - tip4pflag = 1; -} - -/* ---------------------------------------------------------------------- */ - -void PPPMDispTIP4P::init() -{ - // TIP4P PPPM requires newton on, b/c it computes forces on ghost atoms - - if (force->newton == 0) - error->all(FLERR,"Kspace style pppm/disp/tip4p requires newton on"); - - PPPMDisp::init(); -} - -/* ---------------------------------------------------------------------- - find center grid pt for each of my particles - check that full stencil for the particle will fit in my 3d brick - store central grid pt indices in part2grid array -------------------------------------------------------------------------- */ - -void PPPMDispTIP4P::particle_map_c(double delx, double dely, double delz, - double sft, int** p2g, int nup, int nlow, - int nxlo, int nylo, int nzlo, - int nxhi, int nyhi, int nzhi) -{ - int nx,ny,nz,iH1,iH2; - double *xi,xM[3]; - - int *type = atom->type; - double **x = atom->x; - int nlocal = atom->nlocal; - - int flag = 0; - for (int i = 0; i < nlocal; i++) { - if (type[i] == typeO) { - find_M(i,iH1,iH2,xM); - xi = xM; - } else xi = x[i]; - - // (nx,ny,nz) = global coords of grid pt to "lower left" of charge - // current particle coord can be outside global and local box - // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1 - - nx = static_cast ((xi[0]-boxlo[0])*delx+sft) - OFFSET; - ny = static_cast ((xi[1]-boxlo[1])*dely+sft) - OFFSET; - nz = static_cast ((xi[2]-boxlo[2])*delz+sft) - OFFSET; - - p2g[i][0] = nx; - p2g[i][1] = ny; - p2g[i][2] = nz; - - - // check that entire stencil around nx,ny,nz will fit in my 3d brick - - if (nx+nlow < nxlo || nx+nup > nxhi || - ny+nlow < nylo || ny+nup > nyhi || - nz+nlow < nzlo || nz+nup > nzhi) - flag = 1; - } - - if (flag) error->one(FLERR,"Out of range atoms - cannot compute PPPM"); -} - -/* ---------------------------------------------------------------------- - create discretized "density" on section of global grid due to my particles - density(x,y,z) = charge "density" at grid points of my 3d brick - (nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts) - in global grid -------------------------------------------------------------------------- */ - -void PPPMDispTIP4P::make_rho_c() -{ - int i,l,m,n,nx,ny,nz,mx,my,mz,iH1,iH2; - FFT_SCALAR dx,dy,dz,x0,y0,z0; - double *xi,xM[3]; - - // clear 3d density array - - FFT_SCALAR *vec = &density_brick[nzlo_out][nylo_out][nxlo_out]; - for (i = 0; i < ngrid; i++) vec[i] = ZEROF; - - // loop over my charges, add their contribution to 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 - - int *type = atom->type; - double *q = atom->q; - double **x = atom->x; - int nlocal = atom->nlocal; - - for (int 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, order, rho_coeff, rho1d); - - z0 = delvolinv * q[i]; - for (n = nlower; n <= nupper; n++) { - mz = n+nz; - y0 = z0*rho1d[2][n]; - for (m = nlower; m <= nupper; m++) { - my = m+ny; - x0 = y0*rho1d[1][m]; - for (l = nlower; l <= nupper; l++) { - mx = l+nx; - density_brick[mz][my][mx] += x0*rho1d[0][l]; - } - } - } - } -} - -/* ---------------------------------------------------------------------- - interpolate from grid to get electric field & force on my particles - for ik differentiation -------------------------------------------------------------------------- */ - -void PPPMDispTIP4P::fieldforce_c_ik() -{ - int i,l,m,n,nx,ny,nz,mx,my,mz; - FFT_SCALAR dx,dy,dz,x0,y0,z0; - FFT_SCALAR ekx,eky,ekz; - double *xi; - int iH1,iH2; - double xM[3]; - double fx,fy,fz; - - // 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, order, rho_coeff, rho1d); - - ekx = eky = ekz = 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]; - ekx -= x0*vdx_brick[mz][my][mx]; - eky -= x0*vdy_brick[mz][my][mx]; - ekz -= x0*vdz_brick[mz][my][mx]; - } - } - } - - // convert E-field to force - - const double qfactor = force->qqrd2e * scale * q[i]; - if (type[i] != typeO) { - f[i][0] += qfactor*ekx; - f[i][1] += qfactor*eky; - if (slabflag != 2) f[i][2] += qfactor*ekz; - - } else { - fx = qfactor * ekx; - fy = qfactor * eky; - fz = qfactor * ekz; - find_M(i,iH1,iH2,xM); - - f[i][0] += fx*(1 - alpha); - f[i][1] += fy*(1 - alpha); - if (slabflag != 2) f[i][2] += fz*(1 - alpha); - - f[iH1][0] += 0.5*alpha*fx; - f[iH1][1] += 0.5*alpha*fy; - if (slabflag != 2) f[iH1][2] += 0.5*alpha*fz; - - f[iH2][0] += 0.5*alpha*fx; - f[iH2][1] += 0.5*alpha*fy; - if (slabflag != 2) f[iH2][2] += 0.5*alpha*fz; - } - } -} - -/* ---------------------------------------------------------------------- - interpolate from grid to get electric field & force on my particles - for ad scheme -------------------------------------------------------------------------- */ - -void PPPMDispTIP4P::fieldforce_c_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 *prd; - double fx,fy,fz; - double sf; - - 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_slab; - - - - // 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, order, rho_coeff, rho1d); - compute_drho1d(dx,dy,dz, order, drho_coeff, drho1d); - - 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; - if (slabflag != 2) f[i][2] += fz; - - } else { - find_M(i,iH1,iH2,xM); - - f[i][0] += fx*(1 - alpha); - f[i][1] += fy*(1 - alpha); - if (slabflag != 2) f[i][2] += fz*(1 - alpha); - - f[iH1][0] += 0.5*alpha*fx; - f[iH1][1] += 0.5*alpha*fy; - if (slabflag != 2) f[iH1][2] += 0.5*alpha*fz; - - f[iH2][0] += 0.5*alpha*fx; - f[iH2][1] += 0.5*alpha*fy; - if (slabflag != 2) f[iH2][2] += 0.5*alpha*fz; - } - } -} - - -/* ---------------------------------------------------------------------- - interpolate from grid to get electric field & force on my particles -------------------------------------------------------------------------- */ - -void PPPMDispTIP4P::fieldforce_c_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, order, rho_coeff, rho1d); - - 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]; - } - } - } -} - -/* ---------------------------------------------------------------------- - find 2 H atoms bonded to O atom i - compute position xM of fictitious charge site for O atom - also return local indices iH1,iH2 of H atoms -------------------------------------------------------------------------- */ - -void PPPMDispTIP4P::find_M(int i, int &iH1, int &iH2, double *xM) -{ - iH1 = atom->map(atom->tag[i] + 1); - iH2 = atom->map(atom->tag[i] + 2); - - if (iH1 == -1 || iH2 == -1) error->one(FLERR,"TIP4P hydrogen is missing"); - if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) - error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); - - double **x = atom->x; - - double delx1 = x[iH1][0] - x[i][0]; - double dely1 = x[iH1][1] - x[i][1]; - double delz1 = x[iH1][2] - x[i][2]; - domain->minimum_image(delx1,dely1,delz1); - - double delx2 = x[iH2][0] - x[i][0]; - double dely2 = x[iH2][1] - x[i][1]; - double delz2 = x[iH2][2] - x[i][2]; - domain->minimum_image(delx2,dely2,delz2); - - xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2); - xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2); - xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2); -} +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Amalie Frischknecht and Ahmed Ismail (SNL) + Rolf Isele-Holder (Aachen University) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "pppm_disp_tip4p.h" +#include "pppm_disp.h" +#include "atom.h" +#include "domain.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "math_const.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define OFFSET 16384 + +#ifdef FFT_SINGLE +#define ZEROF 0.0f +#define ONEF 1.0f +#else +#define ZEROF 0.0 +#define ONEF 1.0 +#endif + +/* ---------------------------------------------------------------------- */ + +PPPMDispTIP4P::PPPMDispTIP4P(LAMMPS *lmp, int narg, char **arg) : + PPPMDisp(lmp, narg, arg) +{ + triclinic_support = 0; + tip4pflag = 1; +} + +/* ---------------------------------------------------------------------- */ + +void PPPMDispTIP4P::init() +{ + // TIP4P PPPM requires newton on, b/c it computes forces on ghost atoms + + if (force->newton == 0) + error->all(FLERR,"Kspace style pppm/disp/tip4p requires newton on"); + + PPPMDisp::init(); +} + +/* ---------------------------------------------------------------------- + find center grid pt for each of my particles + check that full stencil for the particle will fit in my 3d brick + store central grid pt indices in part2grid array +------------------------------------------------------------------------- */ + +void PPPMDispTIP4P::particle_map_c(double delx, double dely, double delz, + double sft, int** p2g, int nup, int nlow, + int nxlo, int nylo, int nzlo, + int nxhi, int nyhi, int nzhi) +{ + int nx,ny,nz,iH1,iH2; + double *xi,xM[3]; + + int *type = atom->type; + double **x = atom->x; + int nlocal = atom->nlocal; + + int flag = 0; + for (int i = 0; i < nlocal; i++) { + if (type[i] == typeO) { + find_M(i,iH1,iH2,xM); + xi = xM; + } else xi = x[i]; + + // (nx,ny,nz) = global coords of grid pt to "lower left" of charge + // current particle coord can be outside global and local box + // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1 + + nx = static_cast ((xi[0]-boxlo[0])*delx+sft) - OFFSET; + ny = static_cast ((xi[1]-boxlo[1])*dely+sft) - OFFSET; + nz = static_cast ((xi[2]-boxlo[2])*delz+sft) - OFFSET; + + p2g[i][0] = nx; + p2g[i][1] = ny; + p2g[i][2] = nz; + + + // check that entire stencil around nx,ny,nz will fit in my 3d brick + + if (nx+nlow < nxlo || nx+nup > nxhi || + ny+nlow < nylo || ny+nup > nyhi || + nz+nlow < nzlo || nz+nup > nzhi) + flag = 1; + } + + if (flag) error->one(FLERR,"Out of range atoms - cannot compute PPPM"); +} + +/* ---------------------------------------------------------------------- + create discretized "density" on section of global grid due to my particles + density(x,y,z) = charge "density" at grid points of my 3d brick + (nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts) + in global grid +------------------------------------------------------------------------- */ + +void PPPMDispTIP4P::make_rho_c() +{ + int i,l,m,n,nx,ny,nz,mx,my,mz,iH1,iH2; + FFT_SCALAR dx,dy,dz,x0,y0,z0; + double *xi,xM[3]; + + // clear 3d density array + + FFT_SCALAR *vec = &density_brick[nzlo_out][nylo_out][nxlo_out]; + for (i = 0; i < ngrid; i++) vec[i] = ZEROF; + + // loop over my charges, add their contribution to 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 + + int *type = atom->type; + double *q = atom->q; + double **x = atom->x; + int nlocal = atom->nlocal; + + for (int 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, order, rho_coeff, rho1d); + + z0 = delvolinv * q[i]; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + y0 = z0*rho1d[2][n]; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + x0 = y0*rho1d[1][m]; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + density_brick[mz][my][mx] += x0*rho1d[0][l]; + } + } + } + } +} + +/* ---------------------------------------------------------------------- + interpolate from grid to get electric field & force on my particles + for ik differentiation +------------------------------------------------------------------------- */ + +void PPPMDispTIP4P::fieldforce_c_ik() +{ + int i,l,m,n,nx,ny,nz,mx,my,mz; + FFT_SCALAR dx,dy,dz,x0,y0,z0; + FFT_SCALAR ekx,eky,ekz; + double *xi; + int iH1,iH2; + double xM[3]; + double fx,fy,fz; + + // 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, order, rho_coeff, rho1d); + + ekx = eky = ekz = 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]; + ekx -= x0*vdx_brick[mz][my][mx]; + eky -= x0*vdy_brick[mz][my][mx]; + ekz -= x0*vdz_brick[mz][my][mx]; + } + } + } + + // convert E-field to force + + const double qfactor = force->qqrd2e * scale * q[i]; + if (type[i] != typeO) { + f[i][0] += qfactor*ekx; + f[i][1] += qfactor*eky; + if (slabflag != 2) f[i][2] += qfactor*ekz; + + } else { + fx = qfactor * ekx; + fy = qfactor * eky; + fz = qfactor * ekz; + find_M(i,iH1,iH2,xM); + + f[i][0] += fx*(1 - alpha); + f[i][1] += fy*(1 - alpha); + if (slabflag != 2) f[i][2] += fz*(1 - alpha); + + f[iH1][0] += 0.5*alpha*fx; + f[iH1][1] += 0.5*alpha*fy; + if (slabflag != 2) f[iH1][2] += 0.5*alpha*fz; + + f[iH2][0] += 0.5*alpha*fx; + f[iH2][1] += 0.5*alpha*fy; + if (slabflag != 2) f[iH2][2] += 0.5*alpha*fz; + } + } +} + +/* ---------------------------------------------------------------------- + interpolate from grid to get electric field & force on my particles + for ad scheme +------------------------------------------------------------------------- */ + +void PPPMDispTIP4P::fieldforce_c_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 *prd; + double fx,fy,fz; + double sf; + + 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_slab; + + + + // 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, order, rho_coeff, rho1d); + compute_drho1d(dx,dy,dz, order, drho_coeff, drho1d); + + 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; + if (slabflag != 2) f[i][2] += fz; + + } else { + find_M(i,iH1,iH2,xM); + + f[i][0] += fx*(1 - alpha); + f[i][1] += fy*(1 - alpha); + if (slabflag != 2) f[i][2] += fz*(1 - alpha); + + f[iH1][0] += 0.5*alpha*fx; + f[iH1][1] += 0.5*alpha*fy; + if (slabflag != 2) f[iH1][2] += 0.5*alpha*fz; + + f[iH2][0] += 0.5*alpha*fx; + f[iH2][1] += 0.5*alpha*fy; + if (slabflag != 2) f[iH2][2] += 0.5*alpha*fz; + } + } +} + + +/* ---------------------------------------------------------------------- + interpolate from grid to get electric field & force on my particles +------------------------------------------------------------------------- */ + +void PPPMDispTIP4P::fieldforce_c_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, order, rho_coeff, rho1d); + + 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]; + } + } + } + } + + const double qfactor = 0.5*force->qqrd2e * scale * q[i]; + + if (eflag_atom) { + if (type[i] != typeO) { + eatom[i] += qfactor*u_pa; + } else { + eatom[i] += qfactor*u_pa*(1-alpha); + eatom[iH1] += qfactor*u_pa*alpha*0.5; + eatom[iH2] += qfactor*u_pa*alpha*0.5; + } + } + if (vflag_atom) { + if (type[i] != typeO) { + vatom[i][0] += v0*qfactor; + vatom[i][1] += v1*qfactor; + vatom[i][2] += v2*qfactor; + vatom[i][3] += v3*qfactor; + vatom[i][4] += v4*qfactor; + vatom[i][5] += v5*qfactor; + } else { + vatom[i][0] += v0*(1-alpha)*qfactor; + vatom[i][1] += v1*(1-alpha)*qfactor; + vatom[i][2] += v2*(1-alpha)*qfactor; + vatom[i][3] += v3*(1-alpha)*qfactor; + vatom[i][4] += v4*(1-alpha)*qfactor; + vatom[i][5] += v5*(1-alpha)*qfactor; + vatom[iH1][0] += v0*alpha*0.5*qfactor; + vatom[iH1][1] += v1*alpha*0.5*qfactor; + vatom[iH1][2] += v2*alpha*0.5*qfactor; + vatom[iH1][3] += v3*alpha*0.5*qfactor; + vatom[iH1][4] += v4*alpha*0.5*qfactor; + vatom[iH1][5] += v5*alpha*0.5*qfactor; + vatom[iH2][0] += v0*alpha*0.5*qfactor; + vatom[iH2][1] += v1*alpha*0.5*qfactor; + vatom[iH2][2] += v2*alpha*0.5*qfactor; + vatom[iH2][3] += v3*alpha*0.5*qfactor; + vatom[iH2][4] += v4*alpha*0.5*qfactor; + vatom[iH2][5] += v5*alpha*0.5*qfactor; + } + } + } +} + +/* ---------------------------------------------------------------------- + find 2 H atoms bonded to O atom i + compute position xM of fictitious charge site for O atom + also return local indices iH1,iH2 of H atoms +------------------------------------------------------------------------- */ + +void PPPMDispTIP4P::find_M(int i, int &iH1, int &iH2, double *xM) +{ + iH1 = atom->map(atom->tag[i] + 1); + iH2 = atom->map(atom->tag[i] + 2); + + if (iH1 == -1 || iH2 == -1) error->one(FLERR,"TIP4P hydrogen is missing"); + if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) + error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); + + double **x = atom->x; + + double delx1 = x[iH1][0] - x[i][0]; + double dely1 = x[iH1][1] - x[i][1]; + double delz1 = x[iH1][2] - x[i][2]; + domain->minimum_image(delx1,dely1,delz1); + + double delx2 = x[iH2][0] - x[i][0]; + double dely2 = x[iH2][1] - x[i][1]; + double delz2 = x[iH2][2] - x[i][2]; + domain->minimum_image(delx2,dely2,delz2); + + xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2); + xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2); + xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2); +} diff --git a/src/compute_pe_atom.cpp b/src/compute_pe_atom.cpp index 91bdb9120c..7d4797b3ec 100755 --- a/src/compute_pe_atom.cpp +++ b/src/compute_pe_atom.cpp @@ -97,14 +97,17 @@ void ComputePEAtom::compute_peratom() // b/c some bonds/dihedrals call pair::ev_tally with pairwise info // nbond includes ghosts if newton_bond is set // ntotal includes ghosts if either newton flag is set + // KSpace includes ghosts if tip4pflag is set int nlocal = atom->nlocal; int npair = nlocal; int nbond = nlocal; int ntotal = nlocal; + int nkspace = nlocal; if (force->newton) npair += atom->nghost; if (force->newton_bond) nbond += atom->nghost; if (force->newton) ntotal += atom->nghost; + if (force->kspace->tip4pflag) nkspace += atom->nghost; // clear local energy array @@ -137,17 +140,16 @@ void ComputePEAtom::compute_peratom() for (i = 0; i < nbond; i++) energy[i] += eatom[i]; } - // communicate ghost energy between neighbor procs - - if (force->newton) comm->reverse_comm_compute(this); - - // KSpace contribution is already per local atom - if (kspaceflag && force->kspace) { double *eatom = force->kspace->eatom; - for (i = 0; i < nlocal; i++) energy[i] += eatom[i]; + for (i = 0; i < nkspace; i++) energy[i] += eatom[i]; } + // communicate ghost energy between neighbor procs + + if (force->newton || force->kspace->tip4pflag) + comm->reverse_comm_compute(this); + // zero energy of atoms not in group // only do this after comm since ghost contributions must be included diff --git a/src/compute_stress_atom.cpp b/src/compute_stress_atom.cpp index 19a5c31a8c..5eb3f0f94d 100644 --- a/src/compute_stress_atom.cpp +++ b/src/compute_stress_atom.cpp @@ -111,14 +111,17 @@ void ComputeStressAtom::compute_peratom() // b/c some bonds/dihedrals call pair::ev_tally with pairwise info // nbond includes ghosts if newton_bond is set // ntotal includes ghosts if either newton flag is set + // KSpace includes ghosts if tip4pflag is set int nlocal = atom->nlocal; int npair = nlocal; int nbond = nlocal; int ntotal = nlocal; + int nkspace = nlocal; if (force->newton) npair += atom->nghost; if (force->newton_bond) nbond += atom->nghost; if (force->newton) ntotal += atom->nghost; + if (force->kspace->tip4pflag) nkspace += atom->nghost; // clear local stress array @@ -163,6 +166,13 @@ void ComputeStressAtom::compute_peratom() stress[i][j] += vatom[i][j]; } + if (kspaceflag && force->kspace) { + double **vatom = force->kspace->vatom; + for (i = 0; i < nkspace; i++) + for (j = 0; j < 6; j++) + stress[i][j] += vatom[i][j]; + } + // add in per-atom contributions from relevant fixes // skip if vatom = NULL // possible during setup phase if fix has not initialized its vatom yet @@ -180,18 +190,10 @@ void ComputeStressAtom::compute_peratom() } } - // communicate ghost atom virials between neighbor procs + // communicate ghost virials between neighbor procs - if (force->newton) comm->reverse_comm_compute(this); - - // KSpace contribution is already per local atom - - if (kspaceflag && force->kspace) { - double **vatom = force->kspace->vatom; - for (i = 0; i < nlocal; i++) - for (j = 0; j < 6; j++) - stress[i][j] += vatom[i][j]; - } + if (force->newton || force->kspace->tip4pflag) + comm->reverse_comm_compute(this); // zero virial of atoms not in group // only do this after comm since ghost contributions must be included diff --git a/src/force.cpp b/src/force.cpp index 8170c874f8..43ef2cf8e9 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -669,7 +669,6 @@ void Force::bounds(char *str, int nmax, int &nlo, int &nhi) /* ---------------------------------------------------------------------- read a floating point value from a string generate an error if not a legitimate floating point value - called by force fields to check validity of their arguments ------------------------------------------------------------------------- */ double Force::numeric(char *str) @@ -689,7 +688,6 @@ double Force::numeric(char *str) /* ---------------------------------------------------------------------- read an integer value from a string generate an error if not a legitimate integer value - called by force fields to check validity of their arguments ------------------------------------------------------------------------- */ int Force::inumeric(char *str) diff --git a/src/kspace.cpp b/src/kspace.cpp index f7db79568c..c29d9a0795 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -188,12 +188,12 @@ void KSpace::ev_setup(int eflag, int vflag) // reallocate per-atom arrays if necessary - if (eflag_atom && atom->nlocal > maxeatom) { + if (eflag_atom && atom->nmax > maxeatom) { maxeatom = atom->nmax; memory->destroy(eatom); memory->create(eatom,maxeatom,"kspace:eatom"); } - if (vflag_atom && atom->nlocal > maxvatom) { + if (vflag_atom && atom->nmax > maxvatom) { maxvatom = atom->nmax; memory->destroy(vatom); memory->create(vatom,maxvatom,6,"kspace:vatom"); @@ -205,10 +205,12 @@ void KSpace::ev_setup(int eflag, int vflag) if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; if (eflag_atom) { n = atom->nlocal; + if (tip4pflag) n += atom->nghost; for (i = 0; i < n; i++) eatom[i] = 0.0; } if (vflag_atom) { n = atom->nlocal; + if (tip4pflag) n += atom->nghost; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; vatom[i][1] = 0.0; diff --git a/src/lattice.cpp b/src/lattice.cpp index 195cde8f0b..4fcba7cbf9 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -18,6 +18,7 @@ #include "update.h" #include "domain.h" #include "comm.h" +#include "force.h" #include "memory.h" #include "error.h" @@ -161,9 +162,9 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) if (dim == 0) orient = orientx; else if (dim == 1) orient = orienty; else if (dim == 2) orient = orientz; - orient[0] = atoi(arg[iarg+2]); - orient[1] = atoi(arg[iarg+3]); - orient[2] = atoi(arg[iarg+4]); + orient[0] = force->inumeric(arg[iarg+2]); + orient[1] = force->inumeric(arg[iarg+3]); + orient[2] = force->inumeric(arg[iarg+4]); iarg += 5; } else if (strcmp(arg[iarg],"spacing") == 0) {