Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85792412
pppm_tip4p_cg.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Oct 2, 03:55
Size
20 KB
Mime Type
text/x-c
Expires
Fri, Oct 4, 03:55 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
21270966
Attached To
rLAMMPS lammps
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
Log In to Comment