Page MenuHomec4science
Files F90663867

pppm_tip4p_cg.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Nov 3, 16:45

pppm_tip4p_cg.cpp
View Options

/* ----------------------------------------------------------------------
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, Ahmed Ismail and Axel Kohlmeyer
------------------------------------------------------------------------- */
#include "math.h"
#include "pppm_tip4p_cg.h"
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "neighbor.h"
#include "memory.h"
#include "error.h"
#include "commgrid.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define OFFSET 16384
#define SMALLQ 0.00001
enum{REVERSE_RHO};
enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM};
#ifdef FFT_SINGLE
#define ZEROF 0.0f
#define ONEF 1.0f
#else
#define ZEROF 0.0
#define ONEF 1.0
#endif
/* ---------------------------------------------------------------------- */
PPPMTIP4PCG::PPPMTIP4PCG(LAMMPS *lmp, int narg, char **arg) :
PPPMTIP4P(lmp, narg, arg)
{
if ((narg < 1) || (narg > 2))
error->all(FLERR,"Illegal kspace_style pppm/tip4p/cg command");
if (narg == 2) smallq = fabs(force->numeric(FLERR,arg[1]));
else smallq = SMALLQ;
num_charged = -1;
is_charged = NULL;
}
/* ----------------------------------------------------------------------
free all memory
------------------------------------------------------------------------- */
PPPMTIP4PCG::~PPPMTIP4PCG()
{
memory->destroy(is_charged);
}
/* ----------------------------------------------------------------------
compute the PPPM long-range force, energy, virial
------------------------------------------------------------------------- */
void PPPMTIP4PCG::compute(int eflag, int vflag)
{
// set energy/virial flags
// invoke allocate_peratom() if needed for first time
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = evflag_atom = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
if (evflag_atom && !peratom_allocate_flag) {
allocate_peratom();
cg_peratom->ghost_notify();
cg_peratom->setup();
}
// convert atoms from box to lamda coords
if (triclinic == 0) boxlo = domain->boxlo;
else {
boxlo = domain->boxlo_lamda;
domain->x2lamda(atom->nlocal);
}
// extend size of per-atom arrays if necessary
if (atom->nlocal > nmax) {
memory->destroy(part2grid);
memory->destroy(is_charged);
nmax = atom->nmax;
memory->create(part2grid,nmax,3,"pppm:part2grid");
memory->create(is_charged,nmax,"pppm/cg:is_charged");
}
// one time setup message
if (num_charged < 0) {
bigint charged_all, charged_num;
double charged_frac, charged_fmax, charged_fmin;
num_charged=0;
for (int i=0; i < atom->nlocal; ++i)
if (fabs(atom->q[i]) > smallq)
++num_charged;
// get fraction of charged particles per domain
if (atom->nlocal > 0)
charged_frac = static_cast<double>(num_charged) * 100.0
/ static_cast<double>(atom->nlocal);
else
charged_frac = 0.0;
MPI_Reduce(&charged_frac,&charged_fmax,1,MPI_DOUBLE,MPI_MAX,0,world);
MPI_Reduce(&charged_frac,&charged_fmin,1,MPI_DOUBLE,MPI_MIN,0,world);
// get fraction of charged particles overall
charged_num = num_charged;
MPI_Reduce(&charged_num,&charged_all,1,MPI_LMP_BIGINT,MPI_SUM,0,world);
charged_frac = static_cast<double>(charged_all) * 100.0
/ static_cast<double>(atom->natoms);
if (me == 0) {
if (screen)
fprintf(screen,
" PPPM/cg optimization cutoff: %g\n"
" Total charged atoms: %.1f%%\n"
" Min/max charged atoms/proc: %.1f%% %.1f%%\n",
smallq,charged_frac,charged_fmin,charged_fmax);
if (logfile)
fprintf(logfile,
" PPPM/cg optimization cutoff: %g\n"
" Total charged atoms: %.1f%%\n"
" Min/max charged atoms/proc: %.1f%% %.1f%%\n",
smallq,charged_frac,charged_fmin,charged_fmax);
}
}
// only need to rebuild this list after a neighbor list update
if (neighbor->ago == 0) {
num_charged = 0;
for (int i = 0; i < atom->nlocal; ++i) {
if (fabs(atom->q[i]) > smallq) {
is_charged[num_charged] = i;
++num_charged;
}
}
}
// find grid points for all my particles
// map my particle charge onto my local 3d density grid
particle_map();
make_rho();
// all procs communicate density values from their ghost cells
// to fully sum contribution in their 3d bricks
// remap from 3d decomposition to FFT decomposition
cg->reverse_comm(this,REVERSE_RHO);
brick2fft();
// compute potential gradient on my FFT grid and
// portion of e_long on this proc's FFT grid
// return gradients (electric fields) in 3d brick decomposition
// also performs per-atom calculations via poisson_peratom()
poisson();
// all procs communicate E-field values
// to fill ghost cells surrounding their 3d bricks
if (differentiation_flag == 1) cg->forward_comm(this,FORWARD_AD);
else cg->forward_comm(this,FORWARD_IK);
// extra per-atom energy/virial communication
if (evflag_atom) {
if (differentiation_flag == 1 && vflag_atom)
cg_peratom->forward_comm(this,FORWARD_AD_PERATOM);
else if (differentiation_flag == 0)
cg_peratom->forward_comm(this,FORWARD_IK_PERATOM);
}
// calculate the force on my particles
fieldforce();
// extra per-atom energy/virial communication
if (evflag_atom) fieldforce_peratom();
// sum global energy across procs and add in volume-dependent term
const double qscale = force->qqrd2e * scale;
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy = energy_all;
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qscale;
}
// sum global virial across procs
if (vflag_global) {
double virial_all[6];
MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i];
}
// per-atom energy/virial
// energy includes self-energy correction
if (evflag_atom) {
const double * const q = atom->q;
if (eflag_atom) {
for (int j = 0; j < num_charged; j++) {
const int i = is_charged[j];
eatom[i] *= 0.5;
eatom[i] -= g_ewald*q[i]*q[i]/MY_PIS + MY_PI2*q[i]*qsum /
(g_ewald*g_ewald*volume);
eatom[i] *= qscale;
}
}
if (vflag_atom) {
for (int n = 0; n < num_charged; n++) {
const int i = is_charged[n];
for (int j = 0; j < 6; j++) vatom[i][j] *= 0.5*qscale;
}
}
}
// 2d slab correction
if (slabflag == 1) slabcorr();
// convert atoms back from lamda to box coords
if (triclinic) domain->lamda2x(atom->nlocal);
}
/* ----------------------------------------------------------------------
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 PPPMTIP4PCG::particle_map()
{
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 j = 0; j < num_charged; j++) {
int i = is_charged[j];
if (type[i] == typeO) {
find_M_cg(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<int> ((xi[0]-boxlo[0])*delxinv+shift) - OFFSET;
ny = static_cast<int> ((xi[1]-boxlo[1])*delyinv+shift) - OFFSET;
nz = static_cast<int> ((xi[2]-boxlo[2])*delzinv+shift) - OFFSET;
part2grid[i][0] = nx;
part2grid[i][1] = ny;
part2grid[i][2] = nz;
// check that entire stencil around nx,ny,nz will fit in my 3d brick
if (nx+nlower < nxlo_out || nx+nupper > nxhi_out ||
ny+nlower < nylo_out || ny+nupper > nyhi_out ||
nz+nlower < nzlo_out || nz+nupper > nzhi_out) flag++;
}
int flag_all;
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
if (flag_all) error->all(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 PPPMTIP4PCG::make_rho()
{
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 j = 0; j < num_charged; j++) {
int i = is_charged[j];
if (type[i] == typeO) {
find_M_cg(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);
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
------------------------------------------------------------------------- */
void PPPMTIP4PCG::fieldforce_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 (int j = 0; j < num_charged; j++) {
int i = is_charged[j];
if (type[i] == typeO) {
find_M_cg(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);
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_cg(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
------------------------------------------------------------------------- */
void PPPMTIP4PCG::fieldforce_ad()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz;
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;
if (triclinic == 0) prd = domain->prd;
else prd = domain->prd_lamda;
double xprd = prd[0];
double yprd = prd[1];
double zprd = prd[2];
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 (int j = 0; j < num_charged; j++) {
int i = is_charged[j];
if (type[i] == typeO) {
find_M_cg(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);
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.0*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.0*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.0*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_cg(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 PPPMTIP4PCG::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;
int *type = atom->type;
int nlocal = atom->nlocal;
for (int j = 0; j < num_charged; j++) {
int i = is_charged[j];
if (type[i] == typeO) {
find_M_cg(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];
}
}
}
}
/* ----------------------------------------------------------------------
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 PPPMTIP4PCG::find_M_cg(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);
}

Event Timeline

c4science · Help