diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index 821c2191b..9d4dab929 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -1,833 +1,847 @@ /* ---------------------------------------------------------------------- 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: Roy Pollock (LLNL), Paul Crozier (SNL) ------------------------------------------------------------------------- */ #include "mpi.h" #include "stdlib.h" #include "stdio.h" #include "math.h" #include "ewald.h" #include "atom.h" #include "comm.h" #include "force.h" #include "pair.h" #include "domain.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; +#define SMALL 0.00001 + #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) /* ---------------------------------------------------------------------- */ Ewald::Ewald(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) { if (narg != 1) error->all("Illegal kspace_style ewald command"); precision = atof(arg[0]); PI = 4.0*atan(1.0); kmax = 0; kxvecs = kyvecs = kzvecs = NULL; ug = NULL; eg = vg = NULL; sfacrl = sfacim = sfacrl_all = sfacim_all = NULL; nmax = 0; ek = NULL; cs = sn = NULL; kcount = 0; } /* ---------------------------------------------------------------------- free all memory ------------------------------------------------------------------------- */ Ewald::~Ewald() { deallocate(); memory->destroy_2d_double_array(ek); memory->destroy_3d_double_array(cs,-kmax_created); memory->destroy_3d_double_array(sn,-kmax_created); } /* ---------------------------------------------------------------------- */ void Ewald::init() { if (comm->me == 0) { if (screen) fprintf(screen,"Ewald initialization ...\n"); if (logfile) fprintf(logfile,"Ewald initialization ...\n"); } // error check + if (domain->triclinic) error->all("Cannot use Ewald with triclinic box"); if (force->dimension == 2) error->all("Cannot use Ewald with 2d simulation"); + if (atom->q == NULL) + error->all("Must use charged atom style with kspace style"); + if (slabflag == 0 && domain->nonperiodic > 0) error->all("Cannot use nonperiodic boundaries with Ewald"); if (slabflag == 1) { if (domain->xperiodic != 1 || domain->yperiodic != 1 || domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) error->all("Incorrect boundaries with slab Ewald"); } // insure use of long (or table) pair_style with long-range Coulombics // set cutoff to Pair's short-range Coulombic cutoff qqrd2e = force->qqrd2e; Pair *pair = force->pair_match("long"); if (pair == NULL) pair = force->pair_match("table"); if (pair == NULL) error->all("KSpace style is incompatible with Pair style"); double cutoff; pair->extract_long(&cutoff); - // compute qsum & qsqsum + // compute qsum & qsqsum and warn if not charge-neutral - double tmp; + qsum = qsqsum = 0.0; + for (int i = 0; i < atom->nlocal; i++) { + qsum += atom->q[i]; + qsqsum += atom->q[i]*atom->q[i]; + } - qsum = 0.0; - for (int i = 0; i < atom->nlocal; i++) qsum += atom->q[i]; + double tmp; MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); qsum = tmp; - - qsqsum = 0.0; - for (int i = 0; i < atom->nlocal; i++) qsqsum += atom->q[i]*atom->q[i]; MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); qsqsum = tmp; + if (qsqsum == 0.0) + error->all("Cannot use kspace solver on system with no charge"); + if (fabs(qsum) > SMALL && comm->me == 0) { + char str[128]; + sprintf(str,"System is not charge neutral, net charge = %g",qsum); + error->warning(str); + } + // setup K-space resolution g_ewald = (1.35 - 0.15*log(precision))/cutoff; gsqmx = -4.0*g_ewald*g_ewald*log(precision); if (comm->me == 0) { if (screen) fprintf(screen," G vector = %g\n",g_ewald); if (logfile) fprintf(logfile," G vector = %g\n",g_ewald); } } /* ---------------------------------------------------------------------- adjust Ewald coeffs, called initially and whenever volume has changed ------------------------------------------------------------------------- */ void Ewald::setup() { // volume-dependent factors double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; // adjustment of z dimension for 2d slab Ewald // 3d Ewald just uses zprd since slab_volfactor = 1.0 double zprd_slab = zprd*slab_volfactor; volume = xprd * yprd * zprd_slab; unitk[0] = 2.0*PI/xprd; unitk[1] = 2.0*PI/yprd; unitk[2] = 2.0*PI/zprd_slab; // determine kmax // function of current box size, precision, G_ewald (short-range cutoff) int nkxmx = static_cast ((g_ewald*xprd/PI) * sqrt(-log(precision))); int nkymx = static_cast ((g_ewald*yprd/PI) * sqrt(-log(precision))); int nkzmx = static_cast ((g_ewald*zprd_slab/PI) * sqrt(-log(precision))); int kmax_old = kmax; kmax = MAX(nkxmx,nkymx); kmax = MAX(kmax,nkzmx); kmax3d = 4*kmax*kmax*kmax + 6*kmax*kmax + 3*kmax; // if size has grown, reallocate k-dependent and nlocal-dependent arrays if (kmax > kmax_old) { deallocate(); allocate(); memory->destroy_2d_double_array(ek); memory->destroy_3d_double_array(cs,-kmax_created); memory->destroy_3d_double_array(sn,-kmax_created); nmax = atom->nmax; ek = memory->create_2d_double_array(nmax,3,"ewald:ek"); cs = memory->create_3d_double_array(-kmax,kmax,3,nmax,"ewald:cs"); sn = memory->create_3d_double_array(-kmax,kmax,3,nmax,"ewald:sn"); kmax_created = kmax; } // pre-compute Ewald coefficients int kcount_old = kcount; coeffs(); // if array sizes changed, print out new sizes if (kmax != kmax_old || kcount != kcount_old) { if (comm->me == 0) { if (screen) fprintf(screen," vectors: actual 1d max = %d %d %d\n", kcount,kmax,kmax3d); if (logfile) fprintf(logfile," vectors: actual 1d max = %d %d %d\n", kcount,kmax,kmax3d); } } } /* ---------------------------------------------------------------------- compute the Ewald long-range force, energy, virial ------------------------------------------------------------------------- */ void Ewald::compute(int eflag, int vflag) { int i,k,n; energy = 0.0; if (vflag) for (n = 0; n < 6; n++) virial[n] = 0.0; // extend size of nlocal-dependent arrays if necessary int nlocal = atom->nlocal; if (nlocal > nmax) { memory->destroy_2d_double_array(ek); memory->destroy_3d_double_array(cs,-kmax_created); memory->destroy_3d_double_array(sn,-kmax_created); nmax = atom->nmax; ek = memory->create_2d_double_array(nmax,3,"ewald:ek"); cs = memory->create_3d_double_array(-kmax,kmax,3,nmax,"ewald:cs"); sn = memory->create_3d_double_array(-kmax,kmax,3,nmax,"ewald:sn"); kmax_created = kmax; } // partial structure factors on each processor // total structure factor by summing over procs eik_dot_r(); MPI_Allreduce(sfacrl,sfacrl_all,kcount,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(sfacim,sfacim_all,kcount,MPI_DOUBLE,MPI_SUM,world); // K-space portion of electric field // double loop over K-vectors and local atoms double **f = atom->f; double *q = atom->q; int kx,ky,kz; double cypz,sypz,exprl,expim,partial; for (i = 0; i < nlocal; i++) { ek[i][0] = 0.0; ek[i][1] = 0.0; ek[i][2] = 0.0; } for (k = 0; k < kcount; k++) { kx = kxvecs[k]; ky = kyvecs[k]; kz = kzvecs[k]; for (i = 0; i < nlocal; i++) { cypz = cs[ky][1][i]*cs[kz][2][i] - sn[ky][1][i]*sn[kz][2][i]; sypz = sn[ky][1][i]*cs[kz][2][i] + cs[ky][1][i]*sn[kz][2][i]; exprl = cs[kx][0][i]*cypz - sn[kx][0][i]*sypz; expim = sn[kx][0][i]*cypz + cs[kx][0][i]*sypz; partial = expim*sfacrl_all[k] - exprl*sfacim_all[k]; ek[i][0] += partial*eg[k][0]; ek[i][1] += partial*eg[k][1]; ek[i][2] += partial*eg[k][2]; } } // convert E-field to force for (i = 0; i < nlocal; i++) { f[i][0] += qqrd2e*q[i]*ek[i][0]; f[i][1] += qqrd2e*q[i]*ek[i][1]; f[i][2] += qqrd2e*q[i]*ek[i][2]; } // energy if requested if (eflag) { for (k = 0; k < kcount; k++) energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]); PI = 4.0*atan(1.0); energy -= g_ewald*qsqsum/1.772453851 + 0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume); energy *= qqrd2e; } // virial if requested if (vflag) { double uk; for (k = 0; k < kcount; k++) { uk = ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]); for (n = 0; n < 6; n++) virial[n] += uk*vg[k][n]; } for (n = 0; n < 6; n++) virial[n] *= qqrd2e; } if (slabflag) slabcorr(eflag); } /* ---------------------------------------------------------------------- */ void Ewald::eik_dot_r() { int i,k,l,m,n,ic; double cstr1,sstr1,cstr2,sstr2,cstr3,sstr3,cstr4,sstr4; double sqk,clpm,slpm; double **x = atom->x; double *q = atom->q; int nlocal = atom->nlocal; n = 0; // (k,0,0), (0,l,0), (0,0,m) for (ic = 0; ic < 3; ic++) { sqk = unitk[ic]*unitk[ic]; if (sqk <= gsqmx) { cstr1 = 0.0; sstr1 = 0.0; for (i = 0; i < nlocal; i++) { cs[0][ic][i] = 1.0; sn[0][ic][i] = 0.0; cs[1][ic][i] = cos(unitk[ic]*x[i][ic]); sn[1][ic][i] = sin(unitk[ic]*x[i][ic]); cs[-1][ic][i] = cs[1][ic][i]; sn[-1][ic][i] = -sn[1][ic][i]; cstr1 += q[i]*cs[1][ic][i]; sstr1 += q[i]*sn[1][ic][i]; } sfacrl[n] = cstr1; sfacim[n++] = sstr1; } } for (m = 2; m <= kmax; m++) { for (ic = 0; ic < 3; ic++) { sqk = m*unitk[ic] * m*unitk[ic]; if (sqk <= gsqmx) { cstr1 = 0.0; sstr1 = 0.0; for (i = 0; i < nlocal; i++) { cs[m][ic][i] = cs[m-1][ic][i]*cs[1][ic][i] - sn[m-1][ic][i]*sn[1][ic][i]; sn[m][ic][i] = sn[m-1][ic][i]*cs[1][ic][i] + cs[m-1][ic][i]*sn[1][ic][i]; cs[-m][ic][i] = cs[m][ic][i]; sn[-m][ic][i] = -sn[m][ic][i]; cstr1 += q[i]*cs[m][ic][i]; sstr1 += q[i]*sn[m][ic][i]; } sfacrl[n] = cstr1; sfacim[n++] = sstr1; } } } // 1 = (k,l,0), 2 = (k,-l,0) for (k = 1; k <= kmax; k++) { for (l = 1; l <= kmax; l++) { sqk = (k*unitk[0] * k*unitk[0]) + (l*unitk[1] * l*unitk[1]); if (sqk <= gsqmx) { cstr1 = 0.0; sstr1 = 0.0; cstr2 = 0.0; sstr2 = 0.0; for (i = 0; i < nlocal; i++) { cstr1 += q[i]*(cs[k][0][i]*cs[l][1][i] - sn[k][0][i]*sn[l][1][i]); sstr1 += q[i]*(sn[k][0][i]*cs[l][1][i] + cs[k][0][i]*sn[l][1][i]); cstr2 += q[i]*(cs[k][0][i]*cs[l][1][i] + sn[k][0][i]*sn[l][1][i]); sstr2 += q[i]*(sn[k][0][i]*cs[l][1][i] - cs[k][0][i]*sn[l][1][i]); } sfacrl[n] = cstr1; sfacim[n++] = sstr1; sfacrl[n] = cstr2; sfacim[n++] = sstr2; } } } // 1 = (0,l,m), 2 = (0,l,-m) for (l = 1; l <= kmax; l++) { for (m = 1; m <= kmax; m++) { sqk = (l*unitk[1] * l*unitk[1]) + (m*unitk[2] * m*unitk[2]); if (sqk <= gsqmx) { cstr1 = 0.0; sstr1 = 0.0; cstr2 = 0.0; sstr2 = 0.0; for (i = 0; i < nlocal; i++) { cstr1 += q[i]*(cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i]); sstr1 += q[i]*(sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i]); cstr2 += q[i]*(cs[l][1][i]*cs[m][2][i] + sn[l][1][i]*sn[m][2][i]); sstr2 += q[i]*(sn[l][1][i]*cs[m][2][i] - cs[l][1][i]*sn[m][2][i]); } sfacrl[n] = cstr1; sfacim[n++] = sstr1; sfacrl[n] = cstr2; sfacim[n++] = sstr2; } } } // 1 = (k,0,m), 2 = (k,0,-m) for (k = 1; k <= kmax; k++) { for (m = 1; m <= kmax; m++) { sqk = (k*unitk[0] * k*unitk[0]) + (m*unitk[2] * m*unitk[2]); if (sqk <= gsqmx) { cstr1 = 0.0; sstr1 = 0.0; cstr2 = 0.0; sstr2 = 0.0; for (i = 0; i < nlocal; i++) { cstr1 += q[i]*(cs[k][0][i]*cs[m][2][i] - sn[k][0][i]*sn[m][2][i]); sstr1 += q[i]*(sn[k][0][i]*cs[m][2][i] + cs[k][0][i]*sn[m][2][i]); cstr2 += q[i]*(cs[k][0][i]*cs[m][2][i] + sn[k][0][i]*sn[m][2][i]); sstr2 += q[i]*(sn[k][0][i]*cs[m][2][i] - cs[k][0][i]*sn[m][2][i]); } sfacrl[n] = cstr1; sfacim[n++] = sstr1; sfacrl[n] = cstr2; sfacim[n++] = sstr2; } } } // 1 = (k,l,m), 2 = (k,-l,m), 3 = (k,l,-m), 4 = (k,-l,-m) for (k = 1; k <= kmax; k++) { for (l = 1; l <= kmax; l++) { for (m = 1; m <= kmax; m++) { sqk = (k*unitk[0] * k*unitk[0]) + (l*unitk[1] * l*unitk[1]) + (m*unitk[2] * m*unitk[2]); if (sqk <= gsqmx) { cstr1 = 0.0; sstr1 = 0.0; cstr2 = 0.0; sstr2 = 0.0; cstr3 = 0.0; sstr3 = 0.0; cstr4 = 0.0; sstr4 = 0.0; for (i = 0; i < nlocal; i++) { clpm = cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i]; slpm = sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i]; cstr1 += q[i]*(cs[k][0][i]*clpm - sn[k][0][i]*slpm); sstr1 += q[i]*(sn[k][0][i]*clpm + cs[k][0][i]*slpm); clpm = cs[l][1][i]*cs[m][2][i] + sn[l][1][i]*sn[m][2][i]; slpm = -sn[l][1][i]*cs[m][2][i] + cs[l][1][i]*sn[m][2][i]; cstr2 += q[i]*(cs[k][0][i]*clpm - sn[k][0][i]*slpm); sstr2 += q[i]*(sn[k][0][i]*clpm + cs[k][0][i]*slpm); clpm = cs[l][1][i]*cs[m][2][i] + sn[l][1][i]*sn[m][2][i]; slpm = sn[l][1][i]*cs[m][2][i] - cs[l][1][i]*sn[m][2][i]; cstr3 += q[i]*(cs[k][0][i]*clpm - sn[k][0][i]*slpm); sstr3 += q[i]*(sn[k][0][i]*clpm + cs[k][0][i]*slpm); clpm = cs[l][1][i]*cs[m][2][i] - sn[l][1][i]*sn[m][2][i]; slpm = -sn[l][1][i]*cs[m][2][i] - cs[l][1][i]*sn[m][2][i]; cstr4 += q[i]*(cs[k][0][i]*clpm - sn[k][0][i]*slpm); sstr4 += q[i]*(sn[k][0][i]*clpm + cs[k][0][i]*slpm); } sfacrl[n] = cstr1; sfacim[n++] = sstr1; sfacrl[n] = cstr2; sfacim[n++] = sstr2; sfacrl[n] = cstr3; sfacim[n++] = sstr3; sfacrl[n] = cstr4; sfacim[n++] = sstr4; } } } } } /* ---------------------------------------------------------------------- pre-compute coefficients for each Ewald K-vector ------------------------------------------------------------------------- */ void Ewald::coeffs() { int k,l,m; double sqk,vterm; double unitkx = unitk[0]; double unitky = unitk[1]; double unitkz = unitk[2]; double g_ewald_sq_inv = 1.0 / (g_ewald*g_ewald); double preu = 4.0*PI/volume; kcount = 0; // (k,0,0), (0,l,0), (0,0,m) for (m = 1; m <= kmax; m++) { sqk = (m*unitkx) * (m*unitkx); if (sqk <= gsqmx) { kxvecs[kcount] = m; kyvecs[kcount] = 0; kzvecs[kcount] = 0; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*m*ug[kcount]; eg[kcount][1] = 0.0; eg[kcount][2] = 0.0; vterm = -2.0*(1.0/sqk + 0.25*g_ewald_sq_inv); vg[kcount][0] = 1.0 + vterm*(unitkx*m)*(unitkx*m); vg[kcount][1] = 1.0; vg[kcount][2] = 1.0; vg[kcount][3] = 0.0; vg[kcount][4] = 0.0; vg[kcount][5] = 0.0; kcount++; } sqk = (m*unitky) * (m*unitky); if (sqk <= gsqmx) { kxvecs[kcount] = 0; kyvecs[kcount] = m; kzvecs[kcount] = 0; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 0.0; eg[kcount][1] = 2.0*unitky*m*ug[kcount]; eg[kcount][2] = 0.0; vterm = -2.0*(1.0/sqk + 0.25*g_ewald_sq_inv); vg[kcount][0] = 1.0; vg[kcount][1] = 1.0 + vterm*(unitky*m)*(unitky*m); vg[kcount][2] = 1.0; vg[kcount][3] = 0.0; vg[kcount][4] = 0.0; vg[kcount][5] = 0.0; kcount++; } sqk = (m*unitkz) * (m*unitkz); if (sqk <= gsqmx) { kxvecs[kcount] = 0; kyvecs[kcount] = 0; kzvecs[kcount] = m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 0.0; eg[kcount][1] = 0.0; eg[kcount][2] = 2.0*unitkz*m*ug[kcount]; vterm = -2.0*(1.0/sqk + 0.25*g_ewald_sq_inv); vg[kcount][0] = 1.0; vg[kcount][1] = 1.0; vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = 0.0; vg[kcount][4] = 0.0; vg[kcount][5] = 0.0; kcount++; } } // 1 = (k,l,0), 2 = (k,-l,0) for (k = 1; k <= kmax; k++) { for (l = 1; l <= kmax; l++) { sqk = (unitkx*k) * (unitkx*k) + (unitky*l) * (unitky*l); if (sqk <= gsqmx) { kxvecs[kcount] = k; kyvecs[kcount] = l; kzvecs[kcount] = 0; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*k*ug[kcount]; eg[kcount][1] = 2.0*unitky*l*ug[kcount]; eg[kcount][2] = 0.0; vterm = -2.0*(1.0/sqk + 0.25*g_ewald_sq_inv); vg[kcount][0] = 1.0 + vterm*(unitkx*k)*(unitkx*k); vg[kcount][1] = 1.0 + vterm*(unitky*l)*(unitky*l); vg[kcount][2] = 1.0; vg[kcount][3] = vterm*unitkx*k*unitky*l; vg[kcount][4] = 0.0; vg[kcount][5] = 0.0; kcount++; kxvecs[kcount] = k; kyvecs[kcount] = -l; kzvecs[kcount] = 0; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*k*ug[kcount]; eg[kcount][1] = -2.0*unitky*l*ug[kcount]; eg[kcount][2] = 0.0; vg[kcount][0] = 1.0 + vterm*(unitkx*k)*(unitkx*k); vg[kcount][1] = 1.0 + vterm*(unitky*l)*(unitky*l); vg[kcount][2] = 1.0; vg[kcount][3] = -vterm*unitkx*k*unitky*l; vg[kcount][4] = 0.0; vg[kcount][5] = 0.0; kcount++;; } } } // 1 = (0,l,m), 2 = (0,l,-m) for (l = 1; l <= kmax; l++) { for (m = 1; m <= kmax; m++) { sqk = (unitky*l) * (unitky*l) + (unitkz*m) * (unitkz*m); if (sqk <= gsqmx) { kxvecs[kcount] = 0; kyvecs[kcount] = l; kzvecs[kcount] = m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 0.0; eg[kcount][1] = 2.0*unitky*l*ug[kcount]; eg[kcount][2] = 2.0*unitkz*m*ug[kcount]; vterm = -2.0*(1.0/sqk + 0.25*g_ewald_sq_inv); vg[kcount][0] = 1.0; vg[kcount][1] = 1.0 + vterm*(unitky*l)*(unitky*l); vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = 0.0; vg[kcount][4] = 0.0; vg[kcount][5] = vterm*unitky*l*unitkz*m; kcount++; kxvecs[kcount] = 0; kyvecs[kcount] = l; kzvecs[kcount] = -m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 0.0; eg[kcount][1] = 2.0*unitky*l*ug[kcount]; eg[kcount][2] = -2.0*unitkz*m*ug[kcount]; vg[kcount][0] = 1.0; vg[kcount][1] = 1.0 + vterm*(unitky*l)*(unitky*l); vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = 0.0; vg[kcount][4] = 0.0; vg[kcount][5] = -vterm*unitky*l*unitkz*m; kcount++; } } } // 1 = (k,0,m), 2 = (k,0,-m) for (k = 1; k <= kmax; k++) { for (m = 1; m <= kmax; m++) { sqk = (unitkx*k) * (unitkx*k) + (unitkz*m) * (unitkz*m); if (sqk <= gsqmx) { kxvecs[kcount] = k; kyvecs[kcount] = 0; kzvecs[kcount] = m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*k*ug[kcount]; eg[kcount][1] = 0.0; eg[kcount][2] = 2.0*unitkz*m*ug[kcount]; vterm = -2.0*(1.0/sqk + 0.25*g_ewald_sq_inv); vg[kcount][0] = 1.0 + vterm*(unitkx*k)*(unitkx*k); vg[kcount][1] = 1.0; vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = 0.0; vg[kcount][4] = vterm*unitkx*k*unitkz*m; vg[kcount][5] = 0.0; kcount++; kxvecs[kcount] = k; kyvecs[kcount] = 0; kzvecs[kcount] = -m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*k*ug[kcount]; eg[kcount][1] = 0.0; eg[kcount][2] = -2.0*unitkz*m*ug[kcount]; vg[kcount][0] = 1.0 + vterm*(unitkx*k)*(unitkx*k); vg[kcount][1] = 1.0; vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = 0.0; vg[kcount][4] = -vterm*unitkx*k*unitkz*m; vg[kcount][5] = 0.0; kcount++; } } } // 1 = (k,l,m), 2 = (k,-l,m), 3 = (k,l,-m), 4 = (k,-l,-m) for (k = 1; k <= kmax; k++) { for (l = 1; l <= kmax; l++) { for (m = 1; m <= kmax; m++) { sqk = (unitkx*k) * (unitkx*k) + (unitky*l) * (unitky*l) + (unitkz*m) * (unitkz*m); if (sqk <= gsqmx) { kxvecs[kcount] = k; kyvecs[kcount] = l; kzvecs[kcount] = m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*k*ug[kcount]; eg[kcount][1] = 2.0*unitky*l*ug[kcount]; eg[kcount][2] = 2.0*unitkz*m*ug[kcount]; vterm = -2.0*(1.0/sqk + 0.25*g_ewald_sq_inv); vg[kcount][0] = 1.0 + vterm*(unitkx*k)*(unitkx*k); vg[kcount][1] = 1.0 + vterm*(unitky*l)*(unitky*l); vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = vterm*unitkx*k*unitky*l; vg[kcount][4] = vterm*unitkx*k*unitkz*m; vg[kcount][5] = vterm*unitky*l*unitkz*m; kcount++; kxvecs[kcount] = k; kyvecs[kcount] = -l; kzvecs[kcount] = m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*k*ug[kcount]; eg[kcount][1] = -2.0*unitky*l*ug[kcount]; eg[kcount][2] = 2.0*unitkz*m*ug[kcount]; vg[kcount][0] = 1.0 + vterm*(unitkx*k)*(unitkx*k); vg[kcount][1] = 1.0 + vterm*(unitky*l)*(unitky*l); vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = -vterm*unitkx*k*unitky*l; vg[kcount][4] = vterm*unitkx*k*unitkz*m; vg[kcount][5] = -vterm*unitky*l*unitkz*m; kcount++; kxvecs[kcount] = k; kyvecs[kcount] = l; kzvecs[kcount] = -m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*k*ug[kcount]; eg[kcount][1] = 2.0*unitky*l*ug[kcount]; eg[kcount][2] = -2.0*unitkz*m*ug[kcount]; vg[kcount][0] = 1.0 + vterm*(unitkx*k)*(unitkx*k); vg[kcount][1] = 1.0 + vterm*(unitky*l)*(unitky*l); vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = vterm*unitkx*k*unitky*l; vg[kcount][4] = -vterm*unitkx*k*unitkz*m; vg[kcount][5] = -vterm*unitky*l*unitkz*m; kcount++; kxvecs[kcount] = k; kyvecs[kcount] = -l; kzvecs[kcount] = -m; ug[kcount] = preu*exp(-0.25*sqk*g_ewald_sq_inv)/sqk; eg[kcount][0] = 2.0*unitkx*k*ug[kcount]; eg[kcount][1] = -2.0*unitky*l*ug[kcount]; eg[kcount][2] = -2.0*unitkz*m*ug[kcount]; vg[kcount][0] = 1.0 + vterm*(unitkx*k)*(unitkx*k); vg[kcount][1] = 1.0 + vterm*(unitky*l)*(unitky*l); vg[kcount][2] = 1.0 + vterm*(unitkz*m)*(unitkz*m); vg[kcount][3] = -vterm*unitkx*k*unitky*l; vg[kcount][4] = -vterm*unitkx*k*unitkz*m; vg[kcount][5] = vterm*unitky*l*unitkz*m; kcount++;; } } } } } /* ---------------------------------------------------------------------- allocate memory that depends on # of K-vectors ------------------------------------------------------------------------- */ void Ewald::allocate() { kxvecs = new int[kmax3d]; kyvecs = new int[kmax3d]; kzvecs = new int[kmax3d]; ug = new double[kmax3d]; eg = memory->create_2d_double_array(kmax3d,3,"ewald:eg"); vg = memory->create_2d_double_array(kmax3d,6,"ewald:vg"); sfacrl = new double[kmax3d]; sfacim = new double[kmax3d]; sfacrl_all = new double[kmax3d]; sfacim_all = new double[kmax3d]; } /* ---------------------------------------------------------------------- deallocate memory that depends on # of K-vectors ------------------------------------------------------------------------- */ void Ewald::deallocate() { delete [] kxvecs; delete [] kyvecs; delete [] kzvecs; delete [] ug; memory->destroy_2d_double_array(eg); memory->destroy_2d_double_array(vg); delete [] sfacrl; delete [] sfacim; delete [] sfacrl_all; delete [] sfacim_all; } /* ---------------------------------------------------------------------- Slab-geometry correction term to dampen inter-slab interactions between periodically repeating slabs. Yields good approximation to 2-D Ewald if adequate empty space is left between repeating slabs (J. Chem. Phys. 111, 3155). Slabs defined here to be parallel to the xy plane. ------------------------------------------------------------------------- */ void Ewald::slabcorr(int eflag) { // compute local contribution to global dipole moment double *q = atom->q; double **x = atom->x; int nlocal = atom->nlocal; double dipole = 0.0; for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2]; // sum local contributions to get global dipole moment double dipole_all; MPI_Allreduce(&dipole,&dipole_all,1,MPI_DOUBLE,MPI_SUM,world); // compute corrections double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume; if (eflag) energy += qqrd2e*e_slabcorr; // add on force corrections double ffact = -4.0*PI*dipole_all/volume; double **f = atom->f; for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*q[i]*ffact; } /* ---------------------------------------------------------------------- memory usage of local arrays ------------------------------------------------------------------------- */ int Ewald::memory_usage() { int bytes = 3 * kmax3d * sizeof(int); bytes += (1 + 3 + 6) * kmax3d * sizeof(double); bytes += 4 * kmax3d * sizeof(double); bytes += nmax*3 * sizeof(double); bytes += 2 * (2*kmax+1)*3*nmax * sizeof(double); return bytes; } diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp index 3c37218a9..540f9d5d6 100644 --- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp @@ -1,544 +1,544 @@ /* ---------------------------------------------------------------------- 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) ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_lj_cut_coul_long_tip4p.h" #include "angle.h" #include "atom.h" #include "bond.h" #include "comm.h" #include "domain.h" #include "force.h" #include "kspace.h" #include "update.h" #include "respa.h" #include "memory.h" #include "neighbor.h" #include "error.h" using namespace LAMMPS_NS; #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define EWALD_F 1.12837917 #define EWALD_P 0.3275911 #define A1 0.254829592 #define A2 -0.284496736 #define A3 1.421413741 #define A4 -1.453152027 #define A5 1.061405429 /* ---------------------------------------------------------------------- */ PairLJCutCoulLongTIP4P::PairLJCutCoulLongTIP4P(LAMMPS *lmp) : PairLJCutCoulLong(lmp) { single_enable = 0; } /* ---------------------------------------------------------------------- */ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) { int i,j,k,numneigh,itype,jtype,itable; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fraction,table; double delx1,dely1,delz1,delx2,dely2,delz2,delx3,dely3,delz3; double r,r2inv,r6inv,forcecoul,forcelj,cforce,negforce; double factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; double phicoul,philj; int iH1,iH2,jH1,jH2; double xiM[3],xjM[3]; double *x1,*x2; double fO[3],fH[3]; int *neighs; double **f; float rsq; int *int_rsq = (int *) &rsq; eng_vdwl = eng_coul = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = tvirial[i] = 0.0; if (vflag == 2) { f = update->f_pair; tf = atom->f; } else f = atom->f; double **x = atom->x; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; double qqrd2e = force->qqrd2e; // loop over neighbors of my atoms for (i = 0; i < nlocal; i++) { qtmp = q[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; if (itype == typeO) { find_M(i,iH1,iH2,xiM); x1 = xiM; } else x1 = x[i]; neighs = neighbor->firstneigh[i]; numneigh = neighbor->numneigh[i]; for (k = 0; k < numneigh; k++) { j = neighs[k]; if (j < nall) factor_coul = factor_lj = 1.0; else { factor_coul = special_coul[j/nall]; factor_lj = special_lj[j/nall]; j %= nall; } 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]) { r2inv = 1.0/rsq; if (rsq < cut_ljsq[itype][jtype]) { r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); forcelj *= factor_lj * r2inv; f[i][0] += delx*forcelj; f[i][1] += dely*forcelj; f[i][2] += delz*forcelj; f[j][0] -= delx*forcelj; f[j][1] -= dely*forcelj; f[j][2] -= delz*forcelj; if (eflag) { philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; eng_vdwl += factor_lj*philj; } } // adjust rsq for off-site O charge(s) if (itype == typeO || jtype == typeO) { if (jtype == typeO) { find_M(j,jH1,jH2,xjM); x2 = xjM; } else x2 = x[j]; delx = x1[0] - x2[0]; dely = x1[1] - x2[1]; delz = x1[2] - x2[2]; rsq = delx*delx + dely*dely + delz*delz; } // test current rsq against cutoff and compute Coulombic force if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrtf(rsq); r2inv = 1 / rsq; grij = g_ewald * r; expm2 = exp(-grij*grij); t = 1.0 / (1.0 + EWALD_P*grij); erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; prefactor = qqrd2e * qtmp*q[j]/r; forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); if (factor_coul < 1.0) { forcecoul -= (1.0-factor_coul)*prefactor; } } else { r2inv = 1 / rsq; itable = *int_rsq & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq - rtable[itable]) * drtable[itable]; table = ftable[itable] + fraction*dftable[itable]; forcecoul = qtmp*q[j] * table; if (factor_coul < 1.0) { table = ctable[itable] + fraction*dctable[itable]; prefactor = qtmp*q[j] * table; forcecoul -= (1.0-factor_coul)*prefactor; } } cforce = forcecoul * r2inv; // if i,j are not O atoms, force is applied directly // if i or j are O atoms, force is on fictitious atoms // spread force to all 3 atoms in water molecule // formulas due to Feenstra et al, J Comp Chem, 20, 786 (1999) if (itype != typeO) { if (vflag == 0) { f[i][0] += delx * cforce; f[i][1] += dely * cforce; f[i][2] += delz * cforce; } else { tf[i][0] += delx * cforce; tf[i][1] += dely * cforce; tf[i][2] += delz * cforce; tvirial[0] += 0.5 * delx * delx * cforce; tvirial[1] += 0.5 * dely * dely * cforce; tvirial[2] += 0.5 * delz * delz * cforce; tvirial[3] += 0.5 * dely * delx * cforce; tvirial[4] += 0.5 * delz * delx * cforce; tvirial[5] += 0.5 * delz * dely * cforce; } } else { fO[0] = delx*cforce*(1.0-2.0*alpha); fO[1] = dely*cforce*(1.0-2.0*alpha); fO[2] = delz*cforce*(1.0-2.0*alpha); fH[0] = alpha * (delx*cforce); fH[1] = alpha * (dely*cforce); fH[2] = alpha * (delz*cforce); if (vflag == 0) { f[i][0] += fO[0]; f[i][1] += fO[1]; f[i][2] += fO[2]; f[iH1][0] += fH[0]; f[iH1][1] += fH[1]; f[iH1][2] += fH[2]; f[iH2][0] += fH[0]; f[iH2][1] += fH[1]; f[iH2][2] += fH[2]; } else { tf[i][0] += fO[0]; tf[i][1] += fO[1]; tf[i][2] += fO[2]; tf[iH1][0] += fH[0]; tf[iH1][1] += fH[1]; tf[iH1][2] += fH[2]; tf[iH2][0] += fH[0]; tf[iH2][1] += fH[1]; tf[iH2][2] += fH[2]; delx1 = x[i][0] - x2[0]; dely1 = x[i][1] - x2[1]; delz1 = x[i][2] - x2[2]; - domain->minimum_image(&delx1,&dely1,&delz1); + domain->minimum_image(delx1,dely1,delz1); delx2 = x[iH1][0] - x2[0]; dely2 = x[iH1][1] - x2[1]; delz2 = x[iH1][2] - x2[2]; - domain->minimum_image(&delx2,&dely2,&delz2); + domain->minimum_image(delx2,dely2,delz2); delx3 = x[iH2][0] - x2[0]; dely3 = x[iH2][1] - x2[1]; delz3 = x[iH2][2] - x2[2]; - domain->minimum_image(&delx3,&dely3,&delz3); + domain->minimum_image(delx3,dely3,delz3); tvirial[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); tvirial[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); tvirial[2] += 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); tvirial[3] += 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); tvirial[4] += 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); tvirial[5] += 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); } } if (jtype != typeO) { if (vflag == 0) { f[j][0] -= delx * cforce; f[j][1] -= dely * cforce; f[j][2] -= delz * cforce; } else { tf[j][0] -= delx * cforce; tf[j][1] -= dely * cforce; tf[j][2] -= delz * cforce; tvirial[0] += 0.5 * (delx * delx * cforce); tvirial[1] += 0.5 * (dely * dely * cforce); tvirial[2] += 0.5 * (delz * delz * cforce); tvirial[3] += 0.5 * (dely * delx * cforce); tvirial[4] += 0.5 * (delz * delx * cforce); tvirial[5] += 0.5 * (delz * dely * cforce); } } else { negforce = -cforce; fO[0] = delx*negforce*(1.0-2.0*alpha); fO[1] = dely*negforce*(1.0-2.0*alpha); fO[2] = delz*negforce*(1.0-2.0*alpha); fH[0] = alpha * (delx*negforce); fH[1] = alpha * (dely*negforce); fH[2] = alpha * (delz*negforce); if (vflag != 2) { f[j][0] += fO[0]; f[j][1] += fO[1]; f[j][2] += fO[2]; f[jH1][0] += fH[0]; f[jH1][1] += fH[1]; f[jH1][2] += fH[2]; f[jH2][0] += fH[0]; f[jH2][1] += fH[1]; f[jH2][2] += fH[2]; } else { tf[j][0] += fO[0]; tf[j][1] += fO[1]; tf[j][2] += fO[2]; tf[jH1][0] += fH[0]; tf[jH1][1] += fH[1]; tf[jH1][2] += fH[2]; tf[jH2][0] += fH[0]; tf[jH2][1] += fH[1]; tf[jH2][2] += fH[2]; delx1 = x[j][0] - x1[0]; dely1 = x[j][1] - x1[1]; delz1 = x[j][2] - x1[2]; - domain->minimum_image(&delx1,&dely1,&delz1); + domain->minimum_image(delx1,dely1,delz1); delx2 = x[jH1][0] - x1[0]; dely2 = x[jH1][1] - x1[1]; delz2 = x[jH1][2] - x1[2]; - domain->minimum_image(&delx2,&dely2,&delz2); + domain->minimum_image(delx2,dely2,delz2); delx3 = x[jH2][0] - x1[0]; dely3 = x[jH2][1] - x1[1]; delz3 = x[jH2][2] - x1[2]; - domain->minimum_image(&delx3,&dely3,&delz3); + domain->minimum_image(delx3,dely3,delz3); tvirial[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); tvirial[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); tvirial[2] += 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); tvirial[3] += 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); tvirial[4] += 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); tvirial[5] += 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); } } if (eflag) { if (!ncoultablebits || rsq <= tabinnersq) phicoul = prefactor*erfc; else { table = etable[itable] + fraction*detable[itable]; phicoul = qtmp*q[j] * table; } if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; eng_coul += phicoul; } } } } } if (vflag == 2) { virial_compute(); for (int i = 0; i < 6; i++) virial[i] += tvirial[i]; } } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairLJCutCoulLongTIP4P::settings(int narg, char **arg) { if (narg < 6 || narg > 7) error->all("Illegal pair_style command"); typeO = atoi(arg[0]); typeH = atoi(arg[1]); typeB = atoi(arg[2]); typeA = atoi(arg[3]); qdist = atof(arg[4]); cut_lj_global = atof(arg[5]); if (narg == 6) cut_coul = cut_lj_global; else cut_coul = atof(arg[6]); // reset cutoffs that have been explicitly set if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; } } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairLJCutCoulLongTIP4P::init_style() { int i,j; if (atom->tag_enable == 0) error->all("Pair style lj/cut/coul/long/tip4p requires atom IDs"); if (!force->newton_pair) error->all("Pair style lj/cut/coul/long/tip4p requires newton pair on"); if (atom->q == NULL) error->all("Must use charged atom style with this pair style"); cut_coulsq = cut_coul * cut_coul; // set & error check interior rRESPA cutoffs if (strcmp(update->integrate_style,"respa") == 0) { if (((Respa *) update->integrate)->level_inner >= 0) { cut_respa = ((Respa *) update->integrate)->cutoff; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) error->all("Pair cutoff < Respa interior cutoff"); } } else cut_respa = NULL; // insure use of correct KSpace long-range solver, set g_ewald if (force->kspace == NULL) error->all("Pair style is incompatible with KSpace style"); if (strcmp(force->kspace_style,"pppm/tip4p") == 0) g_ewald = force->kspace->g_ewald; else error->all("Pair style is incompatible with KSpace style"); // setup force tables if (ncoultablebits) init_tables(); // set alpha parameter double theta = force->angle->equilibrium_angle(typeA); double blen = force->bond->equilibrium_distance(typeB); alpha = qdist / (2.0 * cos(0.5*theta) * blen); } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJCutCoulLongTIP4P::write_restart_settings(FILE *fp) { fwrite(&typeO,sizeof(int),1,fp); fwrite(&typeH,sizeof(int),1,fp); fwrite(&typeB,sizeof(int),1,fp); fwrite(&typeA,sizeof(int),1,fp); fwrite(&qdist,sizeof(double),1,fp); fwrite(&cut_lj_global,sizeof(double),1,fp); fwrite(&cut_coul,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJCutCoulLongTIP4P::read_restart_settings(FILE *fp) { if (comm->me == 0) { fread(&typeO,sizeof(int),1,fp); fread(&typeH,sizeof(int),1,fp); fread(&typeB,sizeof(int),1,fp); fread(&typeA,sizeof(int),1,fp); fread(&qdist,sizeof(double),1,fp); fread(&cut_lj_global,sizeof(double),1,fp); fread(&cut_coul,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&typeO,1,MPI_INT,0,world); MPI_Bcast(&typeH,1,MPI_INT,0,world); MPI_Bcast(&typeB,1,MPI_INT,0,world); MPI_Bcast(&typeA,1,MPI_INT,0,world); MPI_Bcast(&qdist,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- 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 PairLJCutCoulLongTIP4P::find_M(int i, int &iH1, int &iH2, double *xM) { // test that O is correctly bonded to 2 succesive H atoms iH1 = atom->map(atom->tag[i] + 1); iH2 = atom->map(atom->tag[i] + 2); if (iH1 == -1 || iH2 == -1) error->one("TIP4P hydrogen is missing"); if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) error->one("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); + 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); + domain->minimum_image(delx2,dely2,delz2); xM[0] = x[i][0] + alpha * (delx1 + delx2); xM[1] = x[i][1] + alpha * (dely1 + dely2); xM[2] = x[i][2] + alpha * (delz1 + delz2); } /* ---------------------------------------------------------------------- */ void PairLJCutCoulLongTIP4P::extract_tip4p(double *p_qdist, int *p_typeO, int *p_typeH, int *p_typeA, int *p_typeB) { *p_qdist = qdist; *p_typeO = typeO; *p_typeH = typeH; *p_typeA = typeA; *p_typeB = typeB; } diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 3bd785738..79dec5543 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -1,1859 +1,1895 @@ /* ---------------------------------------------------------------------- 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: Roy Pollock (LLNL), Paul Crozier (SNL) ------------------------------------------------------------------------- */ #include "mpi.h" #include "string.h" #include "stdio.h" #include "stdlib.h" #include "math.h" #include "pppm.h" #include "atom.h" #include "comm.h" #include "neighbor.h" #include "force.h" #include "pair.h" #include "bond.h" #include "angle.h" #include "domain.h" #include "fft3d_wrap.h" #include "remap_wrap.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - #define MAXORDER 7 #define OFFSET 4096 #define SMALL 0.00001 #define LARGE 10000.0 #define EPS_HOC 1.0e-7 +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + /* ---------------------------------------------------------------------- */ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) { if (narg != 1) error->all("Illegal kspace_style pppm command"); precision = atof(arg[0]); PI = 4.0*atan(1.0); nfactors = 3; factors = new int[nfactors]; factors[0] = 2; factors[1] = 3; factors[2] = 5; MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); density_brick = vdx_brick = vdy_brick = vdz_brick = NULL; density_fft = NULL; greensfn = NULL; work1 = work2 = NULL; vg = NULL; fkx = fky = fkz = NULL; buf1 = buf2 = NULL; gf_b = NULL; rho1d = rho_coeff = NULL; fft1 = fft2 = NULL; remap = NULL; nmax = 0; part2grid = NULL; } /* ---------------------------------------------------------------------- free all memory ------------------------------------------------------------------------- */ PPPM::~PPPM() { delete [] factors; deallocate(); memory->destroy_2d_int_array(part2grid); } /* ---------------------------------------------------------------------- called once before run ------------------------------------------------------------------------- */ void PPPM::init() { if (me == 0) { if (screen) fprintf(screen,"PPPM initialization ...\n"); if (logfile) fprintf(logfile,"PPPM initialization ...\n"); } // error check + if (domain->triclinic) + error->all("Cannot (yet) use PPPM with triclinic box"); if (force->dimension == 2) error->all("Cannot use PPPM with 2d simulation"); + if (atom->q == NULL) + error->all("Must use charged atom style with kspace style"); + if (slabflag == 0 && domain->nonperiodic > 0) error->all("Cannot use nonperiodic boundaries with PPPM"); if (slabflag == 1) { if (domain->xperiodic != 1 || domain->yperiodic != 1 || domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) error->all("Incorrect boundaries with slab PPPM"); } if (order > MAXORDER) { char str[128]; sprintf(str,"PPPM order cannot be greater than %d",MAXORDER); error->all(str); } // free all arrays previously allocated deallocate(); // insure use of long (or table) pair_style with long-range Coulombics // set cutoff to Pair's short-range Coulombic cutoff qqrd2e = force->qqrd2e; Pair *pair = force->pair_match("long"); if (pair == NULL) pair = force->pair_match("table"); if (pair == NULL) error->all("KSpace style is incompatible with Pair style"); pair->extract_long(&cutoff); // insure use of TIP4P pair_style with TIP4P long-range Coulombics // set TIP4P params from Pair's params qdist = 0.0; pair = force->pair_match("tip4p"); if (strcmp(force->kspace_style,"pppm/tip4p") != 0 && pair != NULL) error->all("KSpace style is incompatible with Pair style"); if (strcmp(force->kspace_style,"pppm/tip4p") == 0 && pair == NULL) error->all("KSpace style is incompatible with Pair style"); if (pair) { int typeA,typeB; pair->extract_tip4p(&qdist,&typeO,&typeH,&typeA,&typeB); if (force->angle == NULL || force->bond == NULL) error->all("Bond and angle potentials must be defined for TIP4P"); double theta = force->angle->equilibrium_angle(typeA); double blen = force->bond->equilibrium_distance(typeB); alpha = qdist / (2.0 * cos(0.5*theta) * blen); } // compute qsum & qsqsum and warn if not charge-neutral qsum = qsqsum = 0.0; for (int i = 0; i < atom->nlocal; i++) { qsum += atom->q[i]; qsqsum += atom->q[i]*atom->q[i]; } double tmp; MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); qsum = tmp; MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); qsqsum = tmp; + if (qsqsum == 0.0) + error->all("Cannot use kspace solver on system with no charge"); if (fabs(qsum) > SMALL && me == 0) { char str[128]; sprintf(str,"System is not charge neutral, net charge = %g",qsum); error->warning(str); } // setup FFT grid resolution and g_ewald set_grid(); if (nx_pppm >= OFFSET || ny_pppm >= OFFSET || nz_pppm >= OFFSET) error->all("PPPM grid is too large"); // global indices of PPPM grid range from 0 to N-1 // nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of // global PPPM grid that I own without ghost cells // for slab PPPM, assign z grid as if it were not extended nxlo_in = comm->myloc[0]*nx_pppm / comm->procgrid[0]; nxhi_in = (comm->myloc[0]+1)*nx_pppm / comm->procgrid[0] - 1; nylo_in = comm->myloc[1]*ny_pppm / comm->procgrid[1]; nyhi_in = (comm->myloc[1]+1)*ny_pppm / comm->procgrid[1] - 1; nzlo_in = comm->myloc[2] * (static_cast (nz_pppm/slab_volfactor)) / comm->procgrid[2]; nzhi_in = (comm->myloc[2]+1) * (static_cast (nz_pppm/slab_volfactor)) / comm->procgrid[2] - 1; // nlower,nupper = stencil size for mapping particles to PPPM grid nlower = -(order-1)/2; nupper = order/2; // shift values for particle <-> grid mapping // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1 if (order % 2) shift = OFFSET + 0.5; else shift = OFFSET; if (order % 2) shiftone = 0.0; else shiftone = 0.5; // nlo_out,nhi_out = lower/upper limits of the 3d sub-brick of // global PPPM grid that my particles can contribute charge to // effectively nlo_in,nhi_in + ghost cells // nlo,nhi = global coords of grid pt to "lower left" of smallest/largest // position a particle in my box can be at - // particle position bound = subbox + skin/2.0 + qdist + // dist[3] = particle position bound = subbox + skin/2.0 + qdist // qdist = offset due to TIP4P fictitious charge + // convert to triclinic if necessary // nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping // for slab PPPM, assign z grid as if it were not extended - int nlo,nhi; - double cuthalf = 0.5 * neighbor->skin + qdist; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + triclinic = domain->triclinic; + double *prd,*sublo,*subhi; + + if (triclinic == 0) { + prd = domain->prd; + boxlo = domain->boxlo; + sublo = domain->sublo; + subhi = domain->subhi; + } else { + prd = domain->prd_lamda; + boxlo = domain->boxlo_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; double zprd_slab = zprd*slab_volfactor; - nlo = static_cast ((domain->subxlo-cuthalf-domain->boxxlo) * + double dist[3]; + double cuthalf = 0.5 * neighbor->skin + qdist; + if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf; + else { + dist[0] = cuthalf/domain->prd[0]; + dist[1] = cuthalf/domain->prd[1]; + dist[2] = cuthalf/domain->prd[2]; + } + + int nlo,nhi; + + nlo = static_cast ((sublo[0]-dist[0]-boxlo[0]) * nx_pppm/xprd + shift) - OFFSET; - nhi = static_cast ((domain->subxhi+cuthalf-domain->boxxlo) * + nhi = static_cast ((subhi[0]+dist[0]-boxlo[0]) * nx_pppm/xprd + shift) - OFFSET; nxlo_out = nlo + nlower; nxhi_out = nhi + nupper; - nlo = static_cast ((domain->subylo-cuthalf-domain->boxylo) * + nlo = static_cast ((sublo[1]-dist[1]-boxlo[1]) * ny_pppm/yprd + shift) - OFFSET; - nhi = static_cast ((domain->subyhi+cuthalf-domain->boxylo) * + nhi = static_cast ((subhi[1]+dist[1]-boxlo[1]) * ny_pppm/yprd + shift) - OFFSET; nylo_out = nlo + nlower; nyhi_out = nhi + nupper; - nlo = static_cast ((domain->subzlo-cuthalf-domain->boxzlo) * + nlo = static_cast ((sublo[2]-dist[2]-boxlo[2]) * nz_pppm/zprd_slab + shift) - OFFSET; - nhi = static_cast ((domain->subzhi+cuthalf-domain->boxzlo) * + nhi = static_cast ((subhi[2]+dist[2]-boxlo[2]) * nz_pppm/zprd_slab + shift) - OFFSET; nzlo_out = nlo + nlower; nzhi_out = nhi + nupper; // for slab PPPM, change the grid boundary for processors at +z end // to include the empty volume between periodically repeating slabs // for slab PPPM, want charge data communicated from -z proc to +z proc, // but not vice versa, also want field data communicated from +z proc to // -z proc, but not vice versa // this is accomplished by nzhi_in = nzhi_out on +z end (no ghost cells) if (slabflag && ((comm->myloc[2]+1) == (comm->procgrid[2]))) { nzhi_in = nz_pppm - 1; nzhi_out = nz_pppm - 1; } // nlo_ghost,nhi_ghost = # of planes I will recv from 6 directions // that overlay domain I own // proc in that direction tells me via sendrecv() // if no neighbor proc, value comes from self since I have ghosts regardless int nplanes; MPI_Status status; nplanes = nxlo_in - nxlo_out; if (comm->procneigh[0][0] != me) MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[0][0],0, &nxhi_ghost,1,MPI_INT,comm->procneigh[0][1],0,world,&status); else nxhi_ghost = nplanes; nplanes = nxhi_out - nxhi_in; if (comm->procneigh[0][1] != me) MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[0][1],0, &nxlo_ghost,1,MPI_INT,comm->procneigh[0][0],0,world,&status); else nxlo_ghost = nplanes; nplanes = nylo_in - nylo_out; if (comm->procneigh[1][0] != me) MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[1][0],0, &nyhi_ghost,1,MPI_INT,comm->procneigh[1][1],0,world,&status); else nyhi_ghost = nplanes; nplanes = nyhi_out - nyhi_in; if (comm->procneigh[1][1] != me) MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[1][1],0, &nylo_ghost,1,MPI_INT,comm->procneigh[1][0],0,world,&status); else nylo_ghost = nplanes; nplanes = nzlo_in - nzlo_out; if (comm->procneigh[2][0] != me) MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[2][0],0, &nzhi_ghost,1,MPI_INT,comm->procneigh[2][1],0,world,&status); else nzhi_ghost = nplanes; nplanes = nzhi_out - nzhi_in; if (comm->procneigh[2][1] != me) MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[2][1],0, &nzlo_ghost,1,MPI_INT,comm->procneigh[2][0],0,world,&status); else nzlo_ghost = nplanes; // test that ghost overlap is not bigger than my sub-domain int flag = 0; if (nxlo_ghost > nxhi_in-nxlo_in+1) flag = 1; if (nxhi_ghost > nxhi_in-nxlo_in+1) flag = 1; if (nylo_ghost > nyhi_in-nylo_in+1) flag = 1; if (nyhi_ghost > nyhi_in-nylo_in+1) flag = 1; if (nzlo_ghost > nzhi_in-nzlo_in+1) flag = 1; if (nzhi_ghost > nzhi_in-nzlo_in+1) flag = 1; int flag_all; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all("PPPM stencil extends too far, reduce PPPM order"); // decomposition of FFT mesh // global indices range from 0 to N-1 // proc owns entire x-dimension, clump of columns in y,z dimensions // npey_fft,npez_fft = # of procs in y,z dims // if nprocs is small enough, proc can own 1 or more entire xy planes, // else proc owns 2d sub-blocks of yz plane // me_y,me_z = which proc (0-npe_fft-1) I am in y,z dimensions // nlo_fft,nhi_fft = lower/upper limit of the section // of the global FFT mesh that I own int npey_fft,npez_fft; if (nz_pppm >= nprocs) { npey_fft = 1; npez_fft = nprocs; } else procs2grid2d(nprocs,ny_pppm,nz_pppm,&npey_fft,&npez_fft); int me_y = me % npey_fft; int me_z = me / npey_fft; nxlo_fft = 0; nxhi_fft = nx_pppm - 1; nylo_fft = me_y*ny_pppm/npey_fft; nyhi_fft = (me_y+1)*ny_pppm/npey_fft - 1; nzlo_fft = me_z*nz_pppm/npez_fft; nzhi_fft = (me_z+1)*nz_pppm/npez_fft - 1; // PPPM grid for this proc, including ghosts ngrid = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) * (nzhi_out-nzlo_out+1); // FFT arrays on this proc, without ghosts // nfft = FFT points in FFT decomposition on this proc // nfft_brick = FFT points in 3d brick-decomposition on this proc // nfft_both = greater of 2 values nfft = (nxhi_fft-nxlo_fft+1) * (nyhi_fft-nylo_fft+1) * (nzhi_fft-nzlo_fft+1); int nfft_brick = (nxhi_in-nxlo_in+1) * (nyhi_in-nylo_in+1) * (nzhi_in-nzlo_in+1); nfft_both = MAX(nfft,nfft_brick); // buffer space for use in brick2fft and fillbrick // idel = max # of ghost planes to send or recv in +/- dir of each dim // nx,ny,nz = owned planes (including ghosts) in each dim // nxx,nyy,nzz = max # of grid cells to send in each dim // nbuf = max in any dim, augment by 3x for components of vd_xyz in fillbrick int idelx,idely,idelz,nx,ny,nz,nxx,nyy,nzz; idelx = MAX(nxlo_ghost,nxhi_ghost); idelx = MAX(idelx,nxhi_out-nxhi_in); idelx = MAX(idelx,nxlo_in-nxlo_out); idely = MAX(nylo_ghost,nyhi_ghost); idely = MAX(idely,nyhi_out-nyhi_in); idely = MAX(idely,nylo_in-nylo_out); idelz = MAX(nzlo_ghost,nzhi_ghost); idelz = MAX(idelz,nzhi_out-nzhi_in); idelz = MAX(idelz,nzlo_in-nzlo_out); nx = nxhi_out - nxlo_out + 1; ny = nyhi_out - nylo_out + 1; nz = nzhi_out - nzlo_out + 1; nxx = idelx * ny * nz; nyy = idely * nx * nz; nzz = idelz * nx * ny; nbuf = MAX(nxx,nyy); nbuf = MAX(nbuf,nzz); nbuf *= 3; // print stats int ngrid_max,nfft_both_max,nbuf_max; MPI_Allreduce(&ngrid,&ngrid_max,1,MPI_INT,MPI_MAX,world); MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world); MPI_Allreduce(&nbuf,&nbuf_max,1,MPI_INT,MPI_MAX,world); if (me == 0) { if (screen) fprintf(screen," brick FFT buffer size/proc = %d %d %d\n", ngrid_max,nfft_both_max,nbuf_max); if (logfile) fprintf(logfile," brick FFT buffer size/proc = %d %d %d\n", ngrid_max,nfft_both_max,nbuf_max); } // allocate K-space dependent memory allocate(); // pre-compute Green's function denomiator expansion // pre-compute 1d charge distribution coefficients compute_gf_denom(); compute_rho_coeff(); } /* ---------------------------------------------------------------------- adjust PPPM coeffs, called initially and whenever volume has changed ------------------------------------------------------------------------- */ void PPPM::setup() { int i,j,k,l,m,n; + double *prd; // volume-dependent factors + // adjust z dimension for 2d slab PPPM + // z dimension for 3d PPPM is zprd since slab_volfactor = 1.0 - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - - // adjustment of z dimension for 2d slab PPPM - // 3d PPPM just uses zprd since slab_volfactor = 1.0 + 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; - volume = xprd * yprd * zprd_slab; delxinv = nx_pppm/xprd; delyinv = ny_pppm/yprd; delzinv = nz_pppm/zprd_slab; delvolinv = delxinv*delyinv*delzinv; double unitkx = (2.0*PI/xprd); double unitky = (2.0*PI/yprd); double unitkz = (2.0*PI/zprd_slab); // fkx,fky,fkz for my FFT grid pts double per; for (i = nxlo_fft; i <= nxhi_fft; i++) { per = i - nx_pppm*(2*i/nx_pppm); fkx[i] = unitkx*per; } for (i = nylo_fft; i <= nyhi_fft; i++) { per = i - ny_pppm*(2*i/ny_pppm); fky[i] = unitky*per; } for (i = nzlo_fft; i <= nzhi_fft; i++) { per = i - nz_pppm*(2*i/nz_pppm); fkz[i] = unitkz*per; } // virial coefficients double sqk,vterm; n = 0; for (k = nzlo_fft; k <= nzhi_fft; k++) { for (j = nylo_fft; j <= nyhi_fft; j++) { for (i = nxlo_fft; i <= nxhi_fft; i++) { sqk = fkx[i]*fkx[i] + fky[j]*fky[j] + fkz[k]*fkz[k]; if (sqk == 0.0) { vg[n][0] = 0.0; vg[n][1] = 0.0; vg[n][2] = 0.0; vg[n][3] = 0.0; vg[n][4] = 0.0; vg[n][5] = 0.0; } else { vterm = -2.0 * (1.0/sqk + 0.25/(g_ewald*g_ewald)); vg[n][0] = 1.0 + vterm*fkx[i]*fkx[i]; vg[n][1] = 1.0 + vterm*fky[j]*fky[j]; vg[n][2] = 1.0 + vterm*fkz[k]*fkz[k]; vg[n][3] = vterm*fkx[i]*fky[j]; vg[n][4] = vterm*fkx[i]*fkz[k]; vg[n][5] = vterm*fky[j]*fkz[k]; } n++; } } } // modified (Hockney-Eastwood) Coulomb Green's function int nx,ny,nz,kper,lper,mper; double snx,sny,snz,snx2,sny2,snz2; double argx,argy,argz,wx,wy,wz,sx,sy,sz,qx,qy,qz; double sum1,dot1,dot2; double numerator,denominator; int nbx = static_cast ((g_ewald*xprd/(PI*nx_pppm)) * pow(-log(EPS_HOC),0.25)); int nby = static_cast ((g_ewald*yprd/(PI*ny_pppm)) * pow(-log(EPS_HOC),0.25)); int nbz = static_cast ((g_ewald*zprd_slab/(PI*nz_pppm)) * pow(-log(EPS_HOC),0.25)); double form = 1.0; n = 0; for (m = nzlo_fft; m <= nzhi_fft; m++) { mper = m - nz_pppm*(2*m/nz_pppm); snz = sin(0.5*unitkz*mper*zprd_slab/nz_pppm); snz2 = snz*snz; for (l = nylo_fft; l <= nyhi_fft; l++) { lper = l - ny_pppm*(2*l/ny_pppm); sny = sin(0.5*unitky*lper*yprd/ny_pppm); sny2 = sny*sny; for (k = nxlo_fft; k <= nxhi_fft; k++) { kper = k - nx_pppm*(2*k/nx_pppm); snx = sin(0.5*unitkx*kper*xprd/nx_pppm); snx2 = snx*snx; sqk = pow(unitkx*kper,2.0) + pow(unitky*lper,2.0) + pow(unitkz*mper,2.0); if (sqk != 0.0) { numerator = form*12.5663706/sqk; denominator = gf_denom(snx2,sny2,snz2); sum1 = 0.0; for (nx = -nbx; nx <= nbx; nx++) { qx = unitkx*(kper+nx_pppm*nx); sx = exp(-.25*pow(qx/g_ewald,2.0)); wx = 1.0; argx = 0.5*qx*xprd/nx_pppm; if (argx != 0.0) wx = pow(sin(argx)/argx,order); for (ny = -nby; ny <= nby; ny++) { qy = unitky*(lper+ny_pppm*ny); sy = exp(-.25*pow(qy/g_ewald,2.0)); wy = 1.0; argy = 0.5*qy*yprd/ny_pppm; if (argy != 0.0) wy = pow(sin(argy)/argy,order); for (nz = -nbz; nz <= nbz; nz++) { qz = unitkz*(mper+nz_pppm*nz); sz = exp(-.25*pow(qz/g_ewald,2.0)); wz = 1.0; argz = 0.5*qz*zprd_slab/nz_pppm; if (argz != 0.0) wz = pow(sin(argz)/argz,order); dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; dot2 = qx*qx+qy*qy+qz*qz; sum1 += (dot1/dot2) * sx*sy*sz * pow(wx*wy*wz,2.0); } } } greensfn[n++] = numerator*sum1/denominator; } else greensfn[n++] = 0.0; } } } } /* ---------------------------------------------------------------------- compute the PPPM long-range force, energy, virial ------------------------------------------------------------------------- */ void PPPM::compute(int eflag, int vflag) { int i; + // 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 nlocal-dependent arrays if necessary if (atom->nlocal > nmax) { memory->destroy_2d_int_array(part2grid); nmax = atom->nmax; part2grid = memory->create_2d_int_array(nmax,3,"pppm:part2grid"); } energy = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; // 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 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 poisson(eflag,vflag); // all procs communicate E-field values to fill ghost cells // surrounding their 3d bricks fillbrick(); // calculate the force on my particles fieldforce(); // sum energy across procs and add in volume-dependent term if (eflag) { 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/1.772453851 + 0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume); energy *= qqrd2e; } // sum virial across procs if (vflag) { double virial_all[6]; MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*volume*virial_all[i]; } // 2d slab correction if (slabflag) slabcorr(eflag); + + // convert atoms back from lamda to box coords + + if (triclinic) domain->lamda2x(atom->nlocal); } /* ---------------------------------------------------------------------- allocate memory that depends on # of K-vectors and order ------------------------------------------------------------------------- */ void PPPM::allocate() { density_brick = memory->create_3d_double_array(nzlo_out,nzhi_out,nylo_out,nyhi_out, nxlo_out,nxhi_out,"pppm:density_brick"); vdx_brick = memory->create_3d_double_array(nzlo_out,nzhi_out,nylo_out,nyhi_out, nxlo_out,nxhi_out,"pppm:vdx_brick"); vdy_brick = memory->create_3d_double_array(nzlo_out,nzhi_out,nylo_out,nyhi_out, nxlo_out,nxhi_out,"pppm:vdy_brick"); vdz_brick = memory->create_3d_double_array(nzlo_out,nzhi_out,nylo_out,nyhi_out, nxlo_out,nxhi_out,"pppm:vdz_brick"); density_fft = new double[nfft_both]; greensfn = new double[nfft_both]; work1 = new double[2*nfft_both]; work2 = new double[2*nfft_both]; vg = memory->create_2d_double_array(nfft_both,6,"pppm:vg"); fkx = memory->create_1d_double_array(nxlo_fft,nxhi_fft,"pppm:fkx"); fky = memory->create_1d_double_array(nylo_fft,nyhi_fft,"pppm:fky"); fkz = memory->create_1d_double_array(nzlo_fft,nzhi_fft,"pppm:fkz"); buf1 = new double[nbuf]; buf2 = new double[nbuf]; // summation coeffs gf_b = new double[order]; rho1d = memory->create_2d_double_array(3,-order/2,order/2,"pppm:rho1d"); rho_coeff = memory->create_2d_double_array(order,(1-order)/2,order/2, "pppm:rho_coeff"); // create 2 FFTs and a Remap // 1st FFT keeps data in FFT decompostion // 2nd FFT returns data in 3d brick decomposition // remap takes data from 3d brick to FFT decomposition int tmp; fft1 = new FFT3d(lmp,world,nx_pppm,ny_pppm,nz_pppm, nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft, nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft, 0,0,&tmp); fft2 = new FFT3d(lmp,world,nx_pppm,ny_pppm,nz_pppm, nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft, nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, 0,0,&tmp); remap = new Remap(lmp,world, nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft, 1,0,0,2); } /* ---------------------------------------------------------------------- deallocate memory that depends on # of K-vectors and order ------------------------------------------------------------------------- */ void PPPM::deallocate() { memory->destroy_3d_double_array(density_brick,nzlo_out,nylo_out,nxlo_out); memory->destroy_3d_double_array(vdx_brick,nzlo_out,nylo_out,nxlo_out); memory->destroy_3d_double_array(vdy_brick,nzlo_out,nylo_out,nxlo_out); memory->destroy_3d_double_array(vdz_brick,nzlo_out,nylo_out,nxlo_out); delete [] density_fft; delete [] greensfn; delete [] work1; delete [] work2; memory->destroy_2d_double_array(vg); memory->destroy_1d_double_array(fkx,nxlo_fft); memory->destroy_1d_double_array(fky,nylo_fft); memory->destroy_1d_double_array(fkz,nzlo_fft); delete [] buf1; delete [] buf2; delete [] gf_b; memory->destroy_2d_double_array(rho1d,-order/2); memory->destroy_2d_double_array(rho_coeff,(1-order)/2); delete fft1; delete fft2; delete remap; } /* ---------------------------------------------------------------------- set size of FFT grid (nx,ny,nz_pppm) and g_ewald ------------------------------------------------------------------------- */ void PPPM::set_grid() { // see JCP 109, pg. 7698 for derivation of coefficients // higher order coefficients may be computed if needed double **acons = memory->create_2d_double_array(8,7,"pppm:acons"); acons[1][0] = 2.0 / 3.0; acons[2][0] = 1.0 / 50.0; acons[2][1] = 5.0 / 294.0; acons[3][0] = 1.0 / 588.0; acons[3][1] = 7.0 / 1440.0; acons[3][2] = 21.0 / 3872.0; acons[4][0] = 1.0 / 4320.0; acons[4][1] = 3.0 / 1936.0; acons[4][2] = 7601.0 / 2271360.0; acons[4][3] = 143.0 / 28800.0; acons[5][0] = 1.0 / 23232.0; acons[5][1] = 7601.0 / 13628160.0; acons[5][2] = 143.0 / 69120.0; acons[5][3] = 517231.0 / 106536960.0; acons[5][4] = 106640677.0 / 11737571328.0; acons[6][0] = 691.0 / 68140800.0; acons[6][1] = 13.0 / 57600.0; acons[6][2] = 47021.0 / 35512320.0; acons[6][3] = 9694607.0 / 2095994880.0; acons[6][4] = 733191589.0 / 59609088000.0; acons[6][5] = 326190917.0 / 11700633600.0; acons[7][0] = 1.0 / 345600.0; acons[7][1] = 3617.0 / 35512320.0; acons[7][2] = 745739.0 / 838397952.0; acons[7][3] = 56399353.0 / 12773376000.0; acons[7][4] = 25091609.0 / 1560084480.0; acons[7][5] = 1755948832039.0 / 36229939200000.0; acons[7][6] = 4887769399.0 / 37838389248.0; double q2 = qsqsum / force->dielectric; double natoms = atom->natoms; - // adjustment of z dimension for 2d slab PPPM + // use xprd,yprd,zprd even if triclinic so grid size is the same + // adjust z dimension for 2d slab PPPM // 3d PPPM just uses zprd since slab_volfactor = 1.0 double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; double zprd_slab = zprd*slab_volfactor; // make initial g_ewald estimate // based on desired error and real space cutoff // fluid-occupied volume used to estimate real-space error // zprd used rather than zprd_slab if (!gewaldflag) g_ewald = sqrt(-log(precision*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*q2))) / cutoff; // set optimal nx_pppm,ny_pppm,nz_pppm based on order and precision // nz_pppm uses extended zprd_slab instead of zprd double h,h1,h2,err,er1,er2,lpr; int ncount; if (!gridflag) { h = 1.0; h1 = 2.0; ncount = 0; err = LARGE; er1 = 0.0; while (fabs(err) > SMALL) { lpr = rms(h,xprd,natoms,q2,acons); err = log(lpr) - log(precision); er2 = er1; er1 = err; h2 = h1; h1 = h; if ((er1 - er2) == 0.0) h = h1 + er1; else h = h1 + er1*(h2 - h1)/(er1 - er2); ncount++; if (ncount > LARGE) error->all("Cannot compute PPPM X grid spacing"); } nx_pppm = static_cast (xprd/h + 1); ncount = 0; err = LARGE; er1 = 0.0; while (fabs(err) > SMALL) { lpr = rms(h,yprd,natoms,q2,acons); err = log(lpr) - log(precision); er2 = er1; er1 = err; h2 = h1; h1 = h; if ((er1 - er2) == 0.0) h = h1 + er1; else h = h1 + er1*(h2 - h1)/(er1 - er2); ncount++; if (ncount > LARGE) error->all("Cannot compute PPPM Y grid spacing"); } ny_pppm = static_cast (yprd/h + 1); ncount = 0; err = LARGE; er1 = 0.0; while (fabs(err) > SMALL) { lpr = rms(h,zprd_slab,natoms,q2,acons); err = log(lpr) - log(precision); er2 = er1; er1 = err; h2 = h1; h1 = h; if ((er1 - er2) == 0.0) h = h1 + er1; else h = h1 + er1*(h2 - h1)/(er1 - er2); ncount++; if (ncount > LARGE) error->all("Cannot compute PPPM Z grid spacing"); } nz_pppm = static_cast (zprd_slab/h + 1); } // convert grid into sizes that are factorable while (!factorable(nx_pppm)) nx_pppm++; while (!factorable(ny_pppm)) ny_pppm++; while (!factorable(nz_pppm)) nz_pppm++; // adjust g_ewald for new grid double dx = xprd/nx_pppm; double dy = yprd/ny_pppm; double dz = zprd_slab/nz_pppm; if (!gewaldflag) { double gew1,gew2,lprx,lpry,lprz,spr; gew1 = g_ewald + 0.01; ncount = 0; err = LARGE; er1 = 0.0; while (fabs(err) > SMALL) { lprx = rms(dx,xprd,natoms,q2,acons); lpry = rms(dy,yprd,natoms,q2,acons); lprz = rms(dz,zprd_slab,natoms,q2,acons); lpr = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0); spr = 2.0*q2 * exp(-g_ewald*g_ewald*cutoff*cutoff) / sqrt(natoms*cutoff*xprd*yprd*zprd_slab); err = log(lpr) - log(spr); er2 = er1; er1 = err; gew2 = gew1; gew1 = g_ewald; if ((er1 - er2) == 0.0) g_ewald = gew1 + er1; else g_ewald = gew1 + er1*(gew2 - gew1)/(er1 - er2); ncount++; if (ncount > LARGE) error->all("Cannot compute PPPM G"); } } // compute final RMS precision double lprx = rms(dx,xprd,natoms,q2,acons); double lpry = rms(dy,yprd,natoms,q2,acons); double lprz = rms(dz,zprd_slab,natoms,q2,acons); lpr = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0); double spr = 2.0*q2 * exp(-g_ewald*g_ewald*cutoff*cutoff) / sqrt(natoms*cutoff*xprd*yprd*zprd_slab); // free local memory memory->destroy_2d_double_array(acons); // print info if (me == 0) { if (screen) { fprintf(screen," G vector = %g\n",g_ewald); fprintf(screen," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm); fprintf(screen," RMS precision = %g\n",MAX(lpr,spr)); } if (logfile) { fprintf(logfile," G vector = %g\n",g_ewald); fprintf(logfile," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm); fprintf(logfile," RMS precision = %g\n",MAX(lpr,spr)); } } } /* ---------------------------------------------------------------------- check if all factors of n are in list of factors return 1 if yes, 0 if no ------------------------------------------------------------------------- */ int PPPM::factorable(int n) { int i; while (n > 1) { for (i = 0; i < nfactors; i++) { if (n % factors[i] == 0) { n /= factors[i]; break; } } if (i == nfactors) return 0; } return 1; } /* ---------------------------------------------------------------------- compute RMS precision for a dimension ------------------------------------------------------------------------- */ double PPPM::rms(double h, double prd, double natoms, double q2, double **acons) { double sum = 0.0; for (int m = 0; m < order; m++) sum += acons[order][m] * pow(h*g_ewald,2.0*m); double value = q2 * pow(h*g_ewald,order) * sqrt(g_ewald*prd*sqrt(2.0*PI)*sum/natoms) / (prd*prd); return value; } /* ---------------------------------------------------------------------- denominator for Hockney-Eastwood Green's function of x,y,z = sin(kx*deltax/2), etc inf n-1 S(n,k) = Sum W(k+pi*j)**2 = Sum b(l)*(z*z)**l j=-inf l=0 = -(z*z)**n /(2n-1)! * (d/dx)**(2n-1) cot(x) at z = sin(x) gf_b = denominator expansion coeffs ------------------------------------------------------------------------- */ double PPPM::gf_denom(double x, double y, double z) { double sx,sy,sz; sz = sy = sx = 0.0; for (int l = order-1; l >= 0; l--) { sx = gf_b[l] + sx*x; sy = gf_b[l] + sy*y; sz = gf_b[l] + sz*z; } double s = sx*sy*sz; return s*s; } /* ---------------------------------------------------------------------- pre-compute Green's function denominator expansion coeffs, Gamma(2n) ------------------------------------------------------------------------- */ void PPPM::compute_gf_denom() { int k,l,m; for (l = 1; l < order; l++) gf_b[l] = 0.0; gf_b[0] = 1.0; for (m = 1; m < order; m++) { for (l = m; l > 0; l--) gf_b[l] = 4.0 * (gf_b[l]*(l-m)*(l-m-0.5)-gf_b[l-1]*(l-m-1)*(l-m-1)); gf_b[0] = 4.0 * (gf_b[0]*(l-m)*(l-m-0.5)); } int ifact = 1; for (k = 1; k < 2*order; k++) ifact *= k; double gaminv = 1.0/ifact; for (l = 0; l < order; l++) gf_b[l] *= gaminv; } /* ---------------------------------------------------------------------- ghost-swap to accumulate full density in brick decomposition remap density from 3d brick decomposition to FFT decomposition ------------------------------------------------------------------------- */ void PPPM::brick2fft() { int i,n,ix,iy,iz; MPI_Request request; MPI_Status status; // pack my ghosts for +x processor // pass data to self or +x processor // unpack and sum recv data into my real cells n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy <= nyhi_out; iy++) for (ix = nxhi_in+1; ix <= nxhi_out; ix++) buf1[n++] = density_brick[iz][iy][ix]; if (comm->procneigh[0][1] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[0][0],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[0][1],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy <= nyhi_out; iy++) for (ix = nxlo_in; ix < nxlo_in+nxlo_ghost; ix++) density_brick[iz][iy][ix] += buf2[n++]; // pack my ghosts for -x processor // pass data to self or -x processor // unpack and sum recv data into my real cells n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy <= nyhi_out; iy++) for (ix = nxlo_out; ix < nxlo_in; ix++) buf1[n++] = density_brick[iz][iy][ix]; if (comm->procneigh[0][0] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[0][1],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[0][0],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy <= nyhi_out; iy++) for (ix = nxhi_in-nxhi_ghost+1; ix <= nxhi_in; ix++) density_brick[iz][iy][ix] += buf2[n++]; // pack my ghosts for +y processor // pass data to self or +y processor // unpack and sum recv data into my real cells n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nyhi_in+1; iy <= nyhi_out; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) buf1[n++] = density_brick[iz][iy][ix]; if (comm->procneigh[1][1] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[1][0],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[1][1],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_in; iy < nylo_in+nylo_ghost; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) density_brick[iz][iy][ix] += buf2[n++]; // pack my ghosts for -y processor // pass data to self or -y processor // unpack and sum recv data into my real cells n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy < nylo_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) buf1[n++] = density_brick[iz][iy][ix]; if (comm->procneigh[1][0] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[1][1],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[1][0],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nyhi_in-nyhi_ghost+1; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) density_brick[iz][iy][ix] += buf2[n++]; // pack my ghosts for +z processor // pass data to self or +z processor // unpack and sum recv data into my real cells n = 0; for (iz = nzhi_in+1; iz <= nzhi_out; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) buf1[n++] = density_brick[iz][iy][ix]; if (comm->procneigh[2][1] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[2][0],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[2][1],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_in; iz < nzlo_in+nzlo_ghost; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) density_brick[iz][iy][ix] += buf2[n++]; // pack my ghosts for -z processor // pass data to self or -z processor // unpack and sum recv data into my real cells n = 0; for (iz = nzlo_out; iz < nzlo_in; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) buf1[n++] = density_brick[iz][iy][ix]; if (comm->procneigh[2][0] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[2][1],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[2][0],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzhi_in-nzhi_ghost+1; iz <= nzhi_in; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) density_brick[iz][iy][ix] += buf2[n++]; // remap from 3d brick decomposition to FFT decomposition // copy grabs inner portion of density from 3d brick // remap could be done as pre-stage of FFT, // but this works optimally on only double values, not complex values n = 0; for (iz = nzlo_in; iz <= nzhi_in; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) density_fft[n++] = density_brick[iz][iy][ix]; remap->perform(density_fft,density_fft,work1); } /* ---------------------------------------------------------------------- ghost-swap to fill ghost cells of my brick with field values ------------------------------------------------------------------------- */ void PPPM::fillbrick() { int i,n,ix,iy,iz; MPI_Request request; MPI_Status status; // pack my real cells for +z processor // pass data to self or +z processor // unpack and sum recv data into my ghost cells n = 0; for (iz = nzhi_in-nzhi_ghost+1; iz <= nzhi_in; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { buf1[n++] = vdx_brick[iz][iy][ix]; buf1[n++] = vdy_brick[iz][iy][ix]; buf1[n++] = vdz_brick[iz][iy][ix]; } if (comm->procneigh[2][1] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[2][0],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[2][1],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz < nzlo_in; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { vdx_brick[iz][iy][ix] = buf2[n++]; vdy_brick[iz][iy][ix] = buf2[n++]; vdz_brick[iz][iy][ix] = buf2[n++]; } // pack my real cells for -z processor // pass data to self or -z processor // unpack and sum recv data into my ghost cells n = 0; for (iz = nzlo_in; iz < nzlo_in+nzlo_ghost; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { buf1[n++] = vdx_brick[iz][iy][ix]; buf1[n++] = vdy_brick[iz][iy][ix]; buf1[n++] = vdz_brick[iz][iy][ix]; } if (comm->procneigh[2][0] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[2][1],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[2][0],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzhi_in+1; iz <= nzhi_out; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { vdx_brick[iz][iy][ix] = buf2[n++]; vdy_brick[iz][iy][ix] = buf2[n++]; vdz_brick[iz][iy][ix] = buf2[n++]; } // pack my real cells for +y processor // pass data to self or +y processor // unpack and sum recv data into my ghost cells n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nyhi_in-nyhi_ghost+1; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { buf1[n++] = vdx_brick[iz][iy][ix]; buf1[n++] = vdy_brick[iz][iy][ix]; buf1[n++] = vdz_brick[iz][iy][ix]; } if (comm->procneigh[1][1] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[1][0],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[1][1],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy < nylo_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { vdx_brick[iz][iy][ix] = buf2[n++]; vdy_brick[iz][iy][ix] = buf2[n++]; vdz_brick[iz][iy][ix] = buf2[n++]; } // pack my real cells for -y processor // pass data to self or -y processor // unpack and sum recv data into my ghost cells n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_in; iy < nylo_in+nylo_ghost; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { buf1[n++] = vdx_brick[iz][iy][ix]; buf1[n++] = vdy_brick[iz][iy][ix]; buf1[n++] = vdz_brick[iz][iy][ix]; } if (comm->procneigh[1][0] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[1][1],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[1][0],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nyhi_in+1; iy <= nyhi_out; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { vdx_brick[iz][iy][ix] = buf2[n++]; vdy_brick[iz][iy][ix] = buf2[n++]; vdz_brick[iz][iy][ix] = buf2[n++]; } // pack my real cells for +x processor // pass data to self or +x processor // unpack and sum recv data into my ghost cells n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy <= nyhi_out; iy++) for (ix = nxhi_in-nxhi_ghost+1; ix <= nxhi_in; ix++) { buf1[n++] = vdx_brick[iz][iy][ix]; buf1[n++] = vdy_brick[iz][iy][ix]; buf1[n++] = vdz_brick[iz][iy][ix]; } if (comm->procneigh[0][1] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[0][0],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[0][1],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy <= nyhi_out; iy++) for (ix = nxlo_out; ix < nxlo_in; ix++) { vdx_brick[iz][iy][ix] = buf2[n++]; vdy_brick[iz][iy][ix] = buf2[n++]; vdz_brick[iz][iy][ix] = buf2[n++]; } // pack my real cells for -x processor // pass data to self or -x processor // unpack and sum recv data into my ghost cells n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy <= nyhi_out; iy++) for (ix = nxlo_in; ix < nxlo_in+nxlo_ghost; ix++) { buf1[n++] = vdx_brick[iz][iy][ix]; buf1[n++] = vdy_brick[iz][iy][ix]; buf1[n++] = vdz_brick[iz][iy][ix]; } if (comm->procneigh[0][0] == me) for (i = 0; i < n; i++) buf2[i] = buf1[i]; else { MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[0][1],0,world,&request); MPI_Send(buf1,n,MPI_DOUBLE,comm->procneigh[0][0],0,world); MPI_Wait(&request,&status); } n = 0; for (iz = nzlo_out; iz <= nzhi_out; iz++) for (iy = nylo_out; iy <= nyhi_out; iy++) for (ix = nxhi_in+1; ix <= nxhi_out; ix++) { vdx_brick[iz][iy][ix] = buf2[n++]; vdy_brick[iz][iy][ix] = buf2[n++]; vdz_brick[iz][iy][ix] = buf2[n++]; } } /* ---------------------------------------------------------------------- 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 PPPM::particle_map() { int nx,ny,nz; double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; int flag = 0; for (int i = 0; i < nlocal; 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 ((x[i][0]-boxxlo)*delxinv+shift) - OFFSET; - ny = static_cast ((x[i][1]-boxylo)*delyinv+shift) - OFFSET; - nz = static_cast ((x[i][2]-boxzlo)*delzinv+shift) - OFFSET; + nx = static_cast ((x[i][0]-boxlo[0])*delxinv+shift) - OFFSET; + ny = static_cast ((x[i][1]-boxlo[1])*delyinv+shift) - OFFSET; + nz = static_cast ((x[i][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("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 PPPM::make_rho() { int i,l,m,n,nx,ny,nz,mx,my,mz; double dx,dy,dz,x0,y0,z0; // clear 3d density array double *vec = &density_brick[nzlo_out][nylo_out][nxlo_out]; for (i = 0; i < ngrid; i++) vec[i] = 0.0; // 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 double *q = atom->q; double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; for (int i = 0; i < nlocal; i++) { nx = part2grid[i][0]; ny = part2grid[i][1]; nz = part2grid[i][2]; - dx = nx+shiftone - (x[i][0]-boxxlo)*delxinv; - dy = ny+shiftone - (x[i][1]-boxylo)*delyinv; - dz = nz+shiftone - (x[i][2]-boxzlo)*delzinv; + 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); 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]; } } } } } /* ---------------------------------------------------------------------- FFT-based Poisson solver ------------------------------------------------------------------------- */ void PPPM::poisson(int eflag, int vflag) { int i,j,k,n; double eng; // transform charge density (r -> k) n = 0; for (i = 0; i < nfft; i++) { work1[n++] = density_fft[i]; work1[n++] = 0.0; } fft1->compute(work1,work1,1); // if requested, compute energy and virial contribution double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm); double s2 = scaleinv*scaleinv; if (eflag || vflag) { if (vflag) { n = 0; for (i = 0; i < nfft; i++) { eng = s2 * greensfn[i] * (work1[n]*work1[n] + work1[n+1]*work1[n+1]); for (j = 0; j < 6; j++) virial[j] += eng*vg[i][j]; energy += eng; n += 2; } } else { n = 0; for (i = 0; i < nfft; i++) { energy += s2 * greensfn[i] * (work1[n]*work1[n] + work1[n+1]*work1[n+1]); n += 2; } } } // scale by 1/total-grid-pts to get rho(k) // multiply by Green's function to get V(k) n = 0; for (i = 0; i < nfft; i++) { work1[n++] *= scaleinv * greensfn[i]; work1[n++] *= scaleinv * greensfn[i]; } // compute gradients of V(r) in each of 3 dims by transformimg -ik*V(k) // FFT leaves data in 3d brick decomposition // copy it into inner portion of vdx,vdy,vdz arrays // x direction gradient n = 0; for (k = nzlo_fft; k <= nzhi_fft; k++) for (j = nylo_fft; j <= nyhi_fft; j++) for (i = nxlo_fft; i <= nxhi_fft; i++) { work2[n] = fkx[i]*work1[n+1]; work2[n+1] = -fkx[i]*work1[n]; n += 2; } fft2->compute(work2,work2,-1); n = 0; for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { vdx_brick[k][j][i] = work2[n]; n += 2; } // y direction gradient n = 0; for (k = nzlo_fft; k <= nzhi_fft; k++) for (j = nylo_fft; j <= nyhi_fft; j++) for (i = nxlo_fft; i <= nxhi_fft; i++) { work2[n] = fky[j]*work1[n+1]; work2[n+1] = -fky[j]*work1[n]; n += 2; } fft2->compute(work2,work2,-1); n = 0; for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { vdy_brick[k][j][i] = work2[n]; n += 2; } // z direction gradient n = 0; for (k = nzlo_fft; k <= nzhi_fft; k++) for (j = nylo_fft; j <= nyhi_fft; j++) for (i = nxlo_fft; i <= nxhi_fft; i++) { work2[n] = fkz[k]*work1[n+1]; work2[n+1] = -fkz[k]*work1[n]; n += 2; } fft2->compute(work2,work2,-1); n = 0; for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { vdz_brick[k][j][i] = work2[n]; n += 2; } } /* ---------------------------------------------------------------------- interpolate from grid to get electric field & force on my particles ------------------------------------------------------------------------- */ void PPPM::fieldforce() { int i,l,m,n,nx,ny,nz,mx,my,mz; double dx,dy,dz,x0,y0,z0; double ek[3]; // 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 nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; for (i = 0; i < nlocal; i++) { nx = part2grid[i][0]; ny = part2grid[i][1]; nz = part2grid[i][2]; - dx = nx+shiftone - (x[i][0]-boxxlo)*delxinv; - dy = ny+shiftone - (x[i][1]-boxylo)*delyinv; - dz = nz+shiftone - (x[i][2]-boxzlo)*delzinv; + 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); ek[0] = ek[1] = ek[2] = 0.0; 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]; ek[0] -= x0*vdx_brick[mz][my][mx];; ek[1] -= x0*vdy_brick[mz][my][mx];; ek[2] -= x0*vdz_brick[mz][my][mx];; } } } // convert E-field to force f[i][0] += qqrd2e*q[i]*ek[0]; f[i][1] += qqrd2e*q[i]*ek[1]; f[i][2] += qqrd2e*q[i]*ek[2]; } } /* ---------------------------------------------------------------------- map nprocs to NX by NY grid as PX by PY procs - return optimal px,py ------------------------------------------------------------------------- */ void PPPM::procs2grid2d(int nprocs, int nx, int ny, int *px, int *py) { // loop thru all possible factorizations of nprocs // surf = surface area of largest proc sub-domain // innermost if test minimizes surface area and surface/volume ratio int bestsurf = 2 * (nx + ny); int bestboxx = 0; int bestboxy = 0; int boxx,boxy,surf,ipx,ipy; ipx = 1; while (ipx <= nprocs) { if (nprocs % ipx == 0) { ipy = nprocs/ipx; boxx = nx/ipx; if (nx % ipx) boxx++; boxy = ny/ipy; if (ny % ipy) boxy++; surf = boxx + boxy; if (surf < bestsurf || (surf == bestsurf && boxx*boxy > bestboxx*bestboxy)) { bestsurf = surf; bestboxx = boxx; bestboxy = boxy; *px = ipx; *py = ipy; } } ipx++; } } /* ---------------------------------------------------------------------- charge assignment into rho1d dx,dy,dz = distance of particle from "lower left" grid point ------------------------------------------------------------------------- */ void PPPM::compute_rho1d(double dx, double dy, double dz) { int k,l; for (k = (1-order)/2; k <= order/2; k++) { rho1d[0][k] = 0.0; rho1d[1][k] = 0.0; rho1d[2][k] = 0.0; for (l = order-1; l >= 0; l--) { rho1d[0][k] = rho_coeff[l][k] + rho1d[0][k]*dx; rho1d[1][k] = rho_coeff[l][k] + rho1d[1][k]*dy; rho1d[2][k] = rho_coeff[l][k] + rho1d[2][k]*dz; } } } /* ---------------------------------------------------------------------- generate coeffients for the weight function of order n (n-1) Wn(x) = Sum wn(k,x) , Sum is over every other integer k=-(n-1) For k=-(n-1),-(n-1)+2, ....., (n-1)-2,n-1 k is odd integers if n is even and even integers if n is odd --- | n-1 | Sum a(l,j)*(x-k/2)**l if abs(x-k/2) < 1/2 wn(k,x) = < l=0 | | 0 otherwise --- a coeffients are packed into the array rho_coeff to eliminate zeros rho_coeff(l,((k+mod(n+1,2))/2) = a(l,k) ------------------------------------------------------------------------- */ void PPPM::compute_rho_coeff() { int j,k,l,m; double s; double **a = memory->create_2d_double_array(order,-order,order,"pppm:a"); for (k = -order; k <= order; k++) for (l = 0; l < order; l++) a[l][k] = 0.0; a[0][0] = 1.0; for (j = 1; j < order; j++) { for (k = -j; k <= j; k += 2) { s = 0.0; for (l = 0; l < j; l++) { a[l+1][k] = (a[l][k+1]-a[l][k-1]) / (l+1); s += pow(0.5,(double) l+1) * (a[l][k-1] + pow(-1.0,(double) l) * a[l][k+1]) / (l+1); } a[0][k] = s; } } m = (1-order)/2; for (k = -(order-1); k < order; k += 2) { for (l = 0; l < order; l++) rho_coeff[l][m] = a[l][k]; m++; } memory->destroy_2d_double_array(a,-order); } /* ---------------------------------------------------------------------- Slab-geometry correction term to dampen inter-slab interactions between periodically repeating slabs. Yields good approximation to 2D Ewald if adequate empty space is left between repeating slabs (J. Chem. Phys. 111, 3155). Slabs defined here to be parallel to the xy plane. ------------------------------------------------------------------------- */ void PPPM::slabcorr(int eflag) { // compute local contribution to global dipole moment double *q = atom->q; double **x = atom->x; int nlocal = atom->nlocal; double dipole = 0.0; for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2]; // sum local contributions to get global dipole moment double dipole_all; MPI_Allreduce(&dipole,&dipole_all,1,MPI_DOUBLE,MPI_SUM,world); // compute corrections double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume; if (eflag) energy += qqrd2e*e_slabcorr; // add on force corrections double ffact = -4.0*PI*dipole_all/volume; double **f = atom->f; for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*q[i]*ffact; } /* ---------------------------------------------------------------------- perform and time the 4 FFTs required for N timesteps ------------------------------------------------------------------------- */ void PPPM::timing(int n, double &time3d, double &time1d) { double time1,time2; for (int i = 0; i < 2*nfft_both; i++) work1[i] = 0.0; MPI_Barrier(world); time1 = MPI_Wtime(); for (int i = 0; i < n; i++) { fft1->compute(work1,work1,1); fft2->compute(work1,work1,-1); fft2->compute(work1,work1,-1); fft2->compute(work1,work1,-1); } MPI_Barrier(world); time2 = MPI_Wtime(); time3d = time2 - time1; MPI_Barrier(world); time1 = MPI_Wtime(); for (int i = 0; i < n; i++) { fft1->timing1d(work1,nfft_both,1); fft2->timing1d(work1,nfft_both,-1); fft2->timing1d(work1,nfft_both,-1); fft2->timing1d(work1,nfft_both,-1); } MPI_Barrier(world); time2 = MPI_Wtime(); time1d = time2 - time1; } /* ---------------------------------------------------------------------- memory usage of local arrays ------------------------------------------------------------------------- */ int PPPM::memory_usage() { int bytes = nmax*3 * sizeof(double); int nbrick = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) * (nzhi_out-nzlo_out+1); bytes += 4 * nbrick * sizeof(double); bytes += 6 * nfft_both * sizeof(double); bytes += nfft_both*6 * sizeof(double); bytes += 2 * nbuf * sizeof(double); return bytes; } diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index 2acaa9f30..021d6dfc4 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -1,95 +1,97 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef PPPM_H #define PPPM_H #include "kspace.h" namespace LAMMPS_NS { class PPPM : public KSpace { public: PPPM(class LAMMPS *, int, char **); ~PPPM(); void init(); void setup(); void compute(int, int); void timing(int, double &, double &); int memory_usage(); protected: int me,nprocs; double PI; double precision; int nfactors; int *factors; double qsum,qsqsum; double qqrd2e; double cutoff; double volume; double delxinv,delyinv,delzinv,delvolinv; double shift,shiftone; int nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in; int nxlo_out,nylo_out,nzlo_out,nxhi_out,nyhi_out,nzhi_out; int nxlo_ghost,nxhi_ghost,nylo_ghost,nyhi_ghost,nzlo_ghost,nzhi_ghost; int nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft; int nlower,nupper; int ngrid,nfft,nbuf,nfft_both; double ***density_brick; double ***vdx_brick,***vdy_brick,***vdz_brick; double *greensfn; double **vg; double *fkx,*fky,*fkz; double *density_fft; double *work1,*work2; double *buf1,*buf2; double *gf_b; double **rho1d,**rho_coeff; class FFT3d *fft1,*fft2; class Remap *remap; int **part2grid; // storage for particle -> grid mapping int nmax; + int triclinic; // domain settings, orthog or triclinic + double *boxlo; // TIP4P settings int typeH,typeO; // atom types of TIP4P water H and O atoms double qdist; // distance from O site to negative charge double alpha; // geometric factor void set_grid(); void allocate(); void deallocate(); int factorable(int); double rms(double, double, double, double, double **); void compute_gf_denom(); double gf_denom(double, double, double); virtual void particle_map(); virtual void make_rho(); void brick2fft(); void fillbrick(); void poisson(int, int); virtual void fieldforce(); void procs2grid2d(int,int,int,int *, int*); void compute_rho1d(double, double, double); void compute_rho_coeff(); void slabcorr(int); }; } #endif diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp index 3486e9528..dfab00794 100644 --- a/src/KSPACE/pppm_tip4p.cpp +++ b/src/KSPACE/pppm_tip4p.cpp @@ -1,264 +1,255 @@ /* ---------------------------------------------------------------------- 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) ------------------------------------------------------------------------- */ #include "math.h" #include "pppm_tip4p.h" #include "atom.h" #include "domain.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define OFFSET 4096 /* ---------------------------------------------------------------------- */ PPPMTIP4P::PPPMTIP4P(LAMMPS *lmp, int narg, char **arg) : PPPM(lmp, narg, arg) {} /* ---------------------------------------------------------------------- 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 PPPMTIP4P::particle_map() { int nx,ny,nz,iH1,iH2; double *xi,xM[3]; int *type = atom->type; double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; 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]-boxxlo)*delxinv+shift) - OFFSET; - ny = static_cast ((xi[1]-boxylo)*delyinv+shift) - OFFSET; - nz = static_cast ((xi[2]-boxzlo)*delzinv+shift) - OFFSET; + nx = static_cast ((xi[0]-boxlo[0])*delxinv+shift) - OFFSET; + ny = static_cast ((xi[1]-boxlo[1])*delyinv+shift) - OFFSET; + nz = static_cast ((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("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 PPPMTIP4P::make_rho() { int i,l,m,n,nx,ny,nz,mx,my,mz,iH1,iH2; double dx,dy,dz,x0,y0,z0; double *xi,xM[3]; // clear 3d density array double *vec = &density_brick[nzlo_out][nylo_out][nxlo_out]; for (i = 0; i < ngrid; i++) vec[i] = 0.0; // 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; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; 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]-boxxlo)*delxinv; - dy = ny+shiftone - (xi[1]-boxylo)*delyinv; - dz = nz+shiftone - (xi[2]-boxzlo)*delzinv; + 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 ------------------------------------------------------------------------- */ void PPPMTIP4P::fieldforce() { int i,l,m,n,nx,ny,nz,mx,my,mz; double dx,dy,dz,x0,y0,z0; double ek[3]; 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; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; 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]-boxxlo)*delxinv; - dy = ny+shiftone - (xi[1]-boxylo)*delyinv; - dz = nz+shiftone - (xi[2]-boxzlo)*delzinv; + 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); ek[0] = ek[1] = ek[2] = 0.0; 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]; ek[0] -= x0*vdx_brick[mz][my][mx]; ek[1] -= x0*vdy_brick[mz][my][mx]; ek[2] -= x0*vdz_brick[mz][my][mx]; } } } // convert E-field to force if (type[i] != typeO) { f[i][0] += qqrd2e*q[i]*ek[0]; f[i][1] += qqrd2e*q[i]*ek[1]; f[i][2] += qqrd2e*q[i]*ek[2]; } else { fx = qqrd2e * q[i] * ek[0]; fy = qqrd2e * q[i] * ek[1]; fz = qqrd2e * q[i] * ek[2]; find_M(i,iH1,iH2,xM); f[i][0] += fx*(1.0-2.0*alpha); f[i][1] += fy*(1.0-2.0*alpha); f[i][2] += fz*(1.0-2.0*alpha); f[iH1][0] += alpha*(fx); f[iH1][1] += alpha*(fy); f[iH1][2] += alpha*(fz); f[iH2][0] += alpha*(fx); f[iH2][1] += alpha*(fy); f[iH2][2] += alpha*(fz); } } } /* ---------------------------------------------------------------------- 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 PPPMTIP4P::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("TIP4P hydrogen is missing"); if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) error->one("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); + 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); + domain->minimum_image(delx2,dely2,delz2); xM[0] = x[i][0] + alpha * (delx1 + delx2); xM[1] = x[i][1] + alpha * (dely1 + dely2); xM[2] = x[i][2] + alpha * (delz1 + delz2); }