diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 99fe7204ef..1f94f35110 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -60,8 +60,8 @@ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) { if (narg < 1) error->all(FLERR,"Illegal kspace_style pppm command"); - precision = atof(arg[0]); - + accuracy_relative = atof(arg[0]); + nfactors = 3; factors = new int[nfactors]; factors[0] = 2; @@ -222,6 +222,11 @@ void PPPM::init() error->warning(FLERR,str); } + // set accuracy (force units) from accuracy_relative or accuracy_absolute + + if (accuracy_absolute >= 0.0) accuracy = accuracy_absolute; + else accuracy = accuracy_relative * two_charge_force; + // setup FFT grid resolution and g_ewald // normally one iteration thru while loop is all that is required // if grid stencil extends beyond neighbor proc, reduce order and try again @@ -229,7 +234,6 @@ void PPPM::init() int iteration = 0; while (order > 0) { - if (iteration && me == 0) error->warning(FLERR,"Reducing PPPM order b/c stencil extends " "beyond neighbor processor"); @@ -970,9 +974,7 @@ void PPPM::set_grid() acons[7][5] = 1755948832039.0 / 36229939200000.0; acons[7][6] = 4887769399.0 / 37838389248.0; - //double q2 = qsqsum * force->qqrd2e / force->dielectric; - double q2 = qsqsum / force->dielectric; - bigint natoms = atom->natoms; + double q2 = qsqsum * force->qqrd2e / force->dielectric; // use xprd,yprd,zprd even if triclinic so grid size is the same // adjust z dimension for 2d slab PPPM @@ -984,20 +986,21 @@ void PPPM::set_grid() double zprd_slab = zprd*slab_volfactor; // make initial g_ewald estimate - // based on desired error and real space cutoff + // based on desired accuracy and real space cutoff // fluid-occupied volume used to estimate real-space error // zprd used rather than zprd_slab double h_x,h_y,h_z; + bigint natoms = atom->natoms; if (!gewaldflag) - g_ewald = sqrt(-log(precision*sqrt(natoms*cutoff*xprd*yprd*zprd) / + g_ewald = sqrt(-log(accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*q2))) / cutoff; - // set optimal nx_pppm,ny_pppm,nz_pppm based on order and precision + // set optimal nx_pppm,ny_pppm,nz_pppm based on order and accuracy // nz_pppm uses extended zprd_slab instead of zprd // h = 1/g_ewald is upper bound on h such that h*g_ewald <= 1 - // reduce it until precision target is met + // reduce it until accuracy target is met if (!gridflag) { double err; @@ -1008,21 +1011,21 @@ void PPPM::set_grid() nz_pppm = static_cast (zprd_slab/h_z + 1); err = rms(h_x,xprd,natoms,q2,acons); - while (err > precision) { + while (err > accuracy) { err = rms(h_x,xprd,natoms,q2,acons); nx_pppm++; h_x = xprd/nx_pppm; } err = rms(h_y,yprd,natoms,q2,acons); - while (err > precision) { + while (err > accuracy) { err = rms(h_y,yprd,natoms,q2,acons); ny_pppm++; h_y = yprd/ny_pppm; } err = rms(h_z,zprd_slab,natoms,q2,acons); - while (err > precision) { + while (err > accuracy) { err = rms(h_z,zprd_slab,natoms,q2,acons); nz_pppm++; h_z = zprd_slab/nz_pppm; @@ -1067,7 +1070,7 @@ void PPPM::set_grid() } } - // final RMS precision + // final RMS accuracy double lprx = rms(h_x,xprd,natoms,q2,acons); double lpry = rms(h_y,yprd,natoms,q2,acons); @@ -1092,14 +1095,18 @@ void PPPM::set_grid() fprintf(screen," G vector (1/distance)= %g\n",g_ewald); fprintf(screen," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm); fprintf(screen," stencil order = %d\n",order); - fprintf(screen," RMS precision = %g\n",MAX(lpr,spr)); + fprintf(screen," absolute RMS force accuracy = %g\n",MAX(lpr,spr)); + fprintf(screen," relative force accuracy = %g\n", + MAX(lpr,spr)/two_charge_force); fprintf(screen," using %s precision FFTs\n",fft_prec); } if (logfile) { fprintf(logfile," G vector (1/distance) = %g\n",g_ewald); fprintf(logfile," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm); fprintf(logfile," stencil order = %d\n",order); - fprintf(logfile," RMS precision = %g\n",MAX(lpr,spr)); + fprintf(logfile," absolute RMS force accuracy = %g\n",MAX(lpr,spr)); + fprintf(logfile," relative force accuracy = %g\n", + MAX(lpr,spr)/two_charge_force); fprintf(logfile," using %s precision FFTs\n",fft_prec); } } @@ -1128,7 +1135,7 @@ int PPPM::factorable(int n) } /* ---------------------------------------------------------------------- - compute RMS precision for a dimension + compute RMS accuracy for a dimension ------------------------------------------------------------------------- */ double PPPM::rms(double h, double prd, bigint natoms, @@ -1143,7 +1150,7 @@ double PPPM::rms(double h, double prd, bigint natoms, } /* ---------------------------------------------------------------------- - compute difference in real-space and kspace RMS precision + compute difference in real-space and KSpace RMS accuracy ------------------------------------------------------------------------- */ double PPPM::diffpr(double h_x, double h_y, double h_z, double q2, diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index dd327fb170..852f303d5d 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -47,7 +47,6 @@ class PPPM : public KSpace { protected: int me,nprocs; - double precision; int nfactors; int *factors; double qsum,qsqsum; @@ -221,7 +220,7 @@ E: PPPM grid is too large The global PPPM grid is larger than OFFSET in one or more dimensions. OFFSET is currently set to 4096. You likely need to decrease the -requested precision. +requested accuracy. E: PPPM order has been reduced to 0 diff --git a/src/USER-EWALDN/ewald_n.cpp b/src/USER-EWALDN/ewald_n.cpp index fd08295fd4..3f7503db06 100644 --- a/src/USER-EWALDN/ewald_n.cpp +++ b/src/USER-EWALDN/ewald_n.cpp @@ -1,710 +1,1023 @@ -/* ---------------------------------------------------------------------- - 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 author: Pieter J. in 't Veld (SNL) -------------------------------------------------------------------------- */ - -#include "mpi.h" -#include "string.h" -#include "stdio.h" -#include "stdlib.h" -#include "math.h" -#include "ewald_n.h" -#include "math_vector.h" -#include "math_const.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; -using namespace MathConst; - -#define KSPACE_ILLEGAL "Illegal kspace_style ewald/n command" -#define KSPACE_ORDER "Unsupported order in kspace_style ewald/n for" -#define KSPACE_MIX "Unsupported mixing rule in kspace_style ewald/n for" - -enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // same as in pair.cpp - -//#define DEBUG - -/* ---------------------------------------------------------------------- */ - -EwaldN::EwaldN(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) -{ - if (narg!=1) error->all(FLERR,KSPACE_ILLEGAL); - precision = fabs(atof(arg[0])); - memset(function, 0, EWALD_NORDER*sizeof(int)); - kenergy = kvirial = NULL; - cek_local = cek_global = NULL; - ekr_local = NULL; - hvec = NULL; - kvec = NULL; - B = NULL; - first_output = 0; -} - -EwaldN::~EwaldN() -{ - deallocate(); - delete [] ekr_local; - delete [] B; -} - -/* --------------------------------------------------------------------- */ - -void EwaldN::init() -{ - nkvec = nkvec_max = nevec = nevec_max = 0; - nfunctions = nsums = sums = 0; - nbox = -1; - bytes = 0.0; - - if (!comm->me) { // output message - if (screen) fprintf(screen,"EwaldN initialization ...\n"); - if (logfile) fprintf(logfile,"EwaldN initialization ...\n"); - } - - if (domain->dimension == 2) // check for errors - error->all(FLERR,"Cannot use EwaldN with 2d simulation"); - if (slabflag == 0 && domain->nonperiodic > 0) - error->all(FLERR,"Cannot use nonperiodic boundaries with EwaldN"); - if (slabflag == 1) { - if (domain->xperiodic != 1 || domain->yperiodic != 1 || - domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) - error->all(FLERR,"Incorrect boundaries with slab EwaldN"); - } - - scale = 1.0; - //mumurd2e = force->mumurd2e; - //dielectric = force->dielectric; - mumurd2e = dielectric = 1.0; - - int tmp; - Pair *pair = force->pair; - int *ptr = pair ? (int *) pair->extract("ewald_order",tmp) : NULL; - double *cutoff = pair ? (double *) pair->extract("cut_coul",tmp) : NULL; - if (!(ptr||cutoff)) - error->all(FLERR,"KSpace style is incompatible with Pair style"); - int ewald_order = ptr ? *((int *) ptr) : 1<<1; - int ewald_mix = ptr ? *((int *) pair->extract("ewald_mix",tmp)) : GEOMETRIC; - memset(function, 0, EWALD_NFUNCS*sizeof(int)); - for (int i=0; i<=EWALD_NORDER; ++i) // transcribe order - if (ewald_order&(1<pair_style); - error->all(FLERR,str); - default: - sprintf(str, "%s pair_style %s", KSPACE_ORDER, force->pair_style); - error->all(FLERR,str); - } - nfunctions += function[k] = 1; - nsums += n[k]; - } - - g_ewald = (1.35 - 0.15*log(precision))/ *cutoff; // determine resolution - g2_max = -4.0*g_ewald*g_ewald*log(precision); - - if (!comm->me) { // output results - if (screen) fprintf(screen, " G vector = %g\n", g_ewald); - if (logfile) fprintf(logfile, " G vector = %g\n", g_ewald); - } -} - - -/* ---------------------------------------------------------------------- - adjust EwaldN coeffs, called initially and whenever volume has changed -------------------------------------------------------------------------- */ - -void EwaldN::setup() -{ - volume = shape_det(domain->h)*slab_volfactor; // cell volume - memcpy(unit, domain->h_inv, sizeof(shape)); // wave vector units - shape_scalar_mult(unit, 2.0*MY_PI); - unit[2] /= slab_volfactor; - - //int nbox_old = nbox, nkvec_old = nkvec; - if (precision>=1) nbox = 0; - else { - vector n = {1.0, 1.0, 1.0}; // based on cutoff - vec_scalar_mult(n, g_ewald*sqrt(-log(precision))/MY_PI); - shape_vec_dot(n, n, domain->h); - n[2] *= slab_volfactor; - nbox = (int) n[0]; - if (nbox<(int) n[1]) nbox = (int) n[1]; - if (nbox<(int) n[2]) nbox = (int) n[2]; - } - reallocate(); - - coefficients(); // compute coeffs - init_coeffs(); - init_coeff_sums(); - init_self(); - - if (!(first_output||comm->me)) { // output on first - first_output = 1; - if (screen) fprintf(screen, - " vectors: nbox = %d, nkvec = %d\n", nbox, nkvec); - if (logfile) fprintf(logfile, - " vectors: nbox = %d, nkvec = %d\n", nbox, nkvec); - } -} - - -void EwaldN::reallocate() // allocate memory -{ - int ix, iy, iz; - int nkvec_max = nkvec; - vector h; - - nkvec = 0; // determine size(kvec) - int *kflag = new int[(nbox+1)*(2*nbox+1)*(2*nbox+1)]; - int *flag = kflag; - - for (ix=0; ix<=nbox; ++ix) - for (iy=-nbox; iy<=nbox; ++iy) - for (iz=-nbox; iz<=nbox; ++iz) - if (!(ix||iy||iz)) *(flag++) = 0; - else if ((!ix)&&(iy<0)) *(flag++) = 0; - else if ((!(ix||iy))&&(iz<0)) *(flag++) = 0; // use symmetry - else { - h[0] = unit[0]*ix; - h[1] = unit[5]*ix+unit[1]*iy; - h[2] = unit[4]*ix+unit[3]*iy+unit[2]*iz; - if ((*(flag++) = h[0]*h[0]+h[1]*h[1]+h[2]*h[2]<=g2_max)) ++nkvec; - } - - if (nkvec>nkvec_max) { - deallocate(); // free memory - hvec = new hvector[nkvec]; // hvec - bytes += (nkvec-nkvec_max)*sizeof(hvector); - kvec = new kvector[nkvec]; // kvec - bytes += (nkvec-nkvec_max)*sizeof(kvector); - kenergy = new double[nkvec*nfunctions]; // kenergy - bytes += (nkvec-nkvec_max)*nfunctions*sizeof(double); - kvirial = new double[6*nkvec*nfunctions]; // kvirial - bytes += 6*(nkvec-nkvec_max)*nfunctions*sizeof(double); - cek_local = new complex[nkvec*nsums]; // cek_local - bytes += (nkvec-nkvec_max)*nsums*sizeof(complex); - cek_global = new complex[nkvec*nsums]; // cek_global - bytes += (nkvec-nkvec_max)*nsums*sizeof(complex); - nkvec_max = nkvec; - } - - flag = kflag; // create index and - kvector *k = kvec; // wave vectors - hvector *hi = hvec; - for (ix=0; ix<=nbox; ++ix) - for (iy=-nbox; iy<=nbox; ++iy) - for (iz=-nbox; iz<=nbox; ++iz) - if (*(flag++)) { - hi->x = unit[0]*ix; - hi->y = unit[5]*ix+unit[1]*iy; - (hi++)->z = unit[4]*ix+unit[3]*iy+unit[2]*iz; - k->x = ix+nbox; k->y = iy+nbox; (k++)->z = iz+nbox; } - - delete [] kflag; -} - -void EwaldN::reallocate_atoms() -{ - if ((nevec = atom->nmax*(2*nbox+1))<=nevec_max) return; - delete [] ekr_local; - ekr_local = new cvector[nevec]; - bytes += (nevec-nevec_max)*sizeof(cvector); - nevec_max = nevec; -} - - -void EwaldN::deallocate() // free memory -{ - delete [] hvec; hvec = NULL; - delete [] kvec; kvec = NULL; - delete [] kenergy; kenergy = NULL; - delete [] kvirial; kvirial = NULL; - delete [] cek_local; cek_local = NULL; - delete [] cek_global; cek_global = NULL; -} - - -void EwaldN::coefficients() // set up pre-factors -{ - vector h; - hvector *hi = hvec, *nh; - double eta2 = 0.25/(g_ewald*g_ewald); - double b1, b2, expb2, h1, h2, c1, c2; - double *ke = kenergy, *kv = kvirial; - int func0 = function[0], func12 = function[1]||function[2], - func3 = function[3]; - - for (nh = (hi = hvec)+nkvec; hintypes; - - if (function[1]) { // geometric 1/r^6 - double **b = (double **) force->pair->extract("B",tmp); - delete [] B; - B = new double[n+1]; - bytes += (n+1)*sizeof(double); - for (int i=0; i<=n; ++i) B[i] = sqrt(fabs(b[i][i])); - } - if (function[2]) { // arithmetic 1/r^6 - double **epsilon = (double **) force->pair->extract("epsilon",tmp); - double **sigma = (double **) force->pair->extract("sigma",tmp); - if (!(epsilon&&sigma)) - error->all(FLERR,"epsilon or sigma reference not set by pair style in ewald/n"); - double eps_i, sigma_i, sigma_n, *bi = B = new double[7*n+7]; - double c[7] = { - 1.0, sqrt(6.0), sqrt(15.0), sqrt(20.0), sqrt(15.0), sqrt(6.0), 1.0}; - for (int i=0; i<=n; ++i) { - eps_i = sqrt(epsilon[i][i]); - sigma_i = sigma[i][i]; - sigma_n = 1.0; - for (int j=0; j<7; ++j) { - *(bi++) = sigma_n*eps_i*c[j]; sigma_n *= sigma_i; - } - } - } -} - - -void EwaldN::init_coeff_sums() // sums based on atoms -{ - if (sums) return; // calculated only once - sums = 1; - - Sum sum_local[EWALD_MAX_NSUMS]; - - memset(sum_local, 0, EWALD_MAX_NSUMS*sizeof(Sum)); - if (function[0]) { // 1/r - double *q = atom->q, *qn = q+atom->nlocal; - for (double *i=q; itype, *ntype = type+atom->nlocal; - for (int *i=type; itype, *ntype = type+atom->nlocal; - for (int *i=type; imu; - int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) - sum_local[9].x2 += mu[i][3]*mu[i][3]; - } - MPI_Allreduce(sum_local, sum, 2*EWALD_MAX_NSUMS, MPI_DOUBLE, MPI_SUM, world); -} - - -void EwaldN::init_self() -{ - double g1 = g_ewald, g2 = g1*g1, g3 = g1*g2; - - memset(energy_self, 0, EWALD_NFUNCS*sizeof(double)); // self energy - memset(virial_self, 0, EWALD_NFUNCS*sizeof(double)); - const double qscale = force->qqrd2e * scale; - - if (function[0]) { // 1/r - virial_self[0] = -0.5*MY_PI*qscale/(g2*volume)*sum[0].x*sum[0].x; - energy_self[0] = sum[0].x2*qscale*g1/sqrt(MY_PI)-virial_self[0]; - } - if (function[1]) { // geometric 1/r^6 - virial_self[1] = MY_PI*sqrt(MY_PI)*g3/(6.0*volume)*sum[1].x*sum[1].x; - energy_self[1] = -sum[1].x2*g3*g3/12.0+virial_self[1]; - } - if (function[2]) { // arithmetic 1/r^6 - virial_self[2] = MY_PI*sqrt(MY_PI)*g3/(48.0*volume)*(sum[2].x*sum[8].x+ - sum[3].x*sum[7].x+sum[4].x*sum[6].x+0.5*sum[5].x*sum[5].x); - energy_self[2] = -sum[2].x2*g3*g3/3.0+virial_self[2]; - } - if (function[3]) { // dipole - virial_self[3] = 0; // in surface - energy_self[3] = sum[9].x2*mumurd2e*2.0*g3/3.0/sqrt(MY_PI)-virial_self[3]; - } -} - -/* ---------------------------------------------------------------------- - compute the EwaldN long-range force, energy, virial -------------------------------------------------------------------------- */ - -void EwaldN::compute(int eflag, int vflag) -{ - if (!nbox) return; - reallocate_atoms(); - compute_ek(); - compute_force(); - compute_surface(); - compute_energy(eflag); - compute_virial(vflag); -} - - -void EwaldN::compute_ek() -{ - cvector *ekr = ekr_local; - int lbytes = (2*nbox+1)*sizeof(cvector); - hvector *h; - kvector *k, *nk = kvec+nkvec; - cvector *z = new cvector[2*nbox+1]; - cvector z1, *zx, *zy, *zz, *zn = z+2*nbox; - complex *cek, zxyz, zxy, cx; - vector mui; - double *x = atom->x[0], *xn = x+3*atom->nlocal, *q = atom->q, qi, bi, ci[7]; - double *mu = atom->mu ? atom->mu[0] : NULL; - int i, kx, ky, n = nkvec*nsums, *type = atom->type, tri = domain->triclinic; - int func[EWALD_NFUNCS]; - - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - memset(cek_local, 0, n*sizeof(complex)); // reset sums - while (xx, 1, 0); C_SET(zz->y, 1, 0); C_SET(zz->z, 1, 0); // z[0] - if (tri) { // triclinic z[1] - C_ANGLE(z1.x, unit[0]*x[0]+unit[5]*x[1]+unit[4]*x[2]); - C_ANGLE(z1.y, unit[1]*x[1]+unit[3]*x[2]); - C_ANGLE(z1.z, x[2]*unit[2]); x += 3; - } - else { // orthogonal z[1] - C_ANGLE(z1.x, *(x++)*unit[0]); - C_ANGLE(z1.y, *(x++)*unit[1]); - C_ANGLE(z1.z, *(x++)*unit[2]); - } - for (; zzx, zz->x, z1.x); // 3D k-vector - C_RMULT(zy->y, zz->y, z1.y); C_CONJ(zx->y, zy->y); - C_RMULT(zy->z, zz->z, z1.z); C_CONJ(zx->z, zy->z); - } - kx = ky = -1; - cek = cek_local; - if (func[0]) qi = *(q++); - if (func[1]) bi = B[*type]; - if (func[2]) memcpy(ci, B+7*type[0], 7*sizeof(double)); - if (func[3]) { - memcpy(mui, mu, sizeof(vector)); mu += 3; - vec_scalar_mult(mui, mu[3]); - mu += 4; - h = hvec; - } - for (k=kvec; ky) { // based on order in - if (kx!=k->x) cx = z[kx = k->x].x; // reallocate - C_RMULT(zxy, z[ky = k->y].y, cx); - } - C_RMULT(zxyz, z[k->z].z, zxy); - if (func[0]) { - cek->re += zxyz.re*qi; (cek++)->im += zxyz.im*qi; - } - if (func[1]) { - cek->re += zxyz.re*bi; (cek++)->im += zxyz.im*bi; - } - if (func[2]) for (i=0; i<7; ++i) { - cek->re += zxyz.re*ci[i]; (cek++)->im += zxyz.im*ci[i]; - } - if (func[3]) { - register double muk = mui[0]*h->x+mui[1]*h->y+mui[2]*h->z; ++h; - cek->re += zxyz.re*muk; (cek++)->im += zxyz.im*muk; - } - } - ekr = (cvector *) ((char *) memcpy(ekr, z, lbytes)+lbytes); - ++type; - } - MPI_Allreduce(cek_local, cek_global, 2*n, MPI_DOUBLE, MPI_SUM, world); - - delete [] z; -} - -void EwaldN::compute_force() -{ - kvector *k; - hvector *h, *nh; - cvector *z = ekr_local; - vector sum[EWALD_MAX_NSUMS], mui; - complex *cek, zc, zx, zxy; - double *f = atom->f[0], *fn = f+3*atom->nlocal, *q = atom->q, *t = NULL; - double *mu = atom->mu ? atom->mu[0] : NULL; - const double qscale = force->qqrd2e * scale; - double *ke, c[EWALD_NFUNCS] = { - 8.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(12.0*volume), - 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 8.0*MY_PI*mumurd2e/volume}; - double kt = 4.0*pow(g_ewald, 3.0)/3.0/sqrt(MY_PI)/c[3]; - int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type; - int func[EWALD_NFUNCS]; - - if (atom->torque) t = atom->torque[0]; - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector)); // fj = -dE/dr = - for (; fy) { // based on order in - if (kx!=k->x) zx = z[kx = k->x].x; // reallocate - C_RMULT(zxy, z[ky = k->y].y, zx); - } - C_CRMULT(zc, z[k->z].z, zxy); - if (func[0]) { // 1/r - register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek; - sum[0][0] += h->x*im; sum[0][1] += h->y*im; sum[0][2] += h->z*im; - } - if (func[1]) { // geometric 1/r^6 - register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek; - sum[1][0] += h->x*im; sum[1][1] += h->y*im; sum[1][2] += h->z*im; - } - if (func[2]) { // arithmetic 1/r^6 - register double im, c = *(ke++); - for (i=2; i<9; ++i) { - im = c*(zc.im*cek->re+cek->im*zc.re); ++cek; - sum[i][0] += h->x*im; sum[i][1] += h->y*im; sum[i][2] += h->z*im; - } - } - if (func[3]) { // dipole - register double im = *(ke++)*(zc.im*cek->re+ - cek->im*zc.re)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z); ++cek; - sum[9][0] += h->x*im; sum[9][1] += h->y*im; sum[9][2] += h->z*im; - } - } - if (func[0]) { // 1/r - register double qi = *(q++)*c[0]; - f[0] -= sum[0][0]*qi; f[1] -= sum[0][1]*qi; f[2] -= sum[0][2]*qi; - } - if (func[1]) { // geometric 1/r^6 - register double bi = B[*type]*c[1]; - f[0] -= sum[1][0]*bi; f[1] -= sum[1][1]*bi; f[2] -= sum[1][2]*bi; - } - if (func[2]) { // arithmetic 1/r^6 - register double *bi = B+7*type[0]+7; - for (i=2; i<9; ++i) { - register double c2 = (--bi)[0]*c[2]; - f[0] -= sum[i][0]*c2; f[1] -= sum[i][1]*c2; f[2] -= sum[i][2]*c2; - } - } - if (func[3]) { // dipole - f[0] -= sum[9][0]; f[1] -= sum[9][1]; f[2] -= sum[9][2]; - *(t++) -= mui[1]*sum[0][2]+mui[2]*sum[0][1]-mui[0]*kt; // torque - *(t++) -= mui[2]*sum[0][0]+mui[0]*sum[0][2]-mui[1]*kt; - *(t++) -= mui[0]*sum[0][1]+mui[1]*sum[0][0]-mui[2]*kt; - } - z = (cvector *) ((char *) z+lbytes); - ++type; - } -} - -void EwaldN::compute_surface() -{ - if (!function[3]) return; - - vector sum_local = VECTOR_NULL, sum_total; - double **mu = atom->mu; - int nlocal = atom->nlocal; - - for (int i=0; i < nlocal; i++) { - register double di = mu[i][3]; - sum_local[0] += di*mu[i][0]; - sum_local[1] += di*mu[i][1]; - sum_local[2] += di*mu[i][2]; - } - MPI_Allreduce(sum_local, sum_total, 3, MPI_DOUBLE, MPI_SUM, world); - - energy_self[3] += virial_self[3]; - virial_self[3] = - mumurd2e*(2.0*MY_PI*vec_dot(sum_total,sum_total)/(2.0*dielectric+1)/volume); - energy_self[3] -= virial_self[3]; -} - -void EwaldN::compute_energy(int eflag) -{ - energy = 0.0; - if (!eflag) return; - - complex *cek = cek_global; - double *ke = kenergy; - const double qscale = force->qqrd2e * scale; - double c[EWALD_NFUNCS] = { - 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), - 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; - double sum[EWALD_NFUNCS]; - int func[EWALD_NFUNCS]; - - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - memset(sum, 0, EWALD_NFUNCS*sizeof(double)); // reset sums - for (int k=0; kre*cek->re+cek->im*cek->im); ++cek; } - if (func[1]) { // geometric 1/r^6 - sum[1] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; } - if (func[2]) { // arithmetic 1/r^6 - register double r = - (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+ - (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+ - (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+ - 0.5*(cek[3].re*cek[3].re+cek[3].im*cek[3].im); cek += 7; - sum[2] += *(ke++)*r; - } - if (func[3]) { // dipole - sum[3] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; } - } - for (int k=0; kqqrd2e * scale; - double c[EWALD_NFUNCS] = { - 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), - 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; - shape sum[EWALD_NFUNCS]; - int func[EWALD_NFUNCS]; - - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - memset(sum, 0, EWALD_NFUNCS*sizeof(shape)); - for (int k=0; kre*cek->re+cek->im*cek->im; ++cek; - sum[0][0] += *(kv++)*r; sum[0][1] += *(kv++)*r; sum[0][2] += *(kv++)*r; - sum[0][3] += *(kv++)*r; sum[0][4] += *(kv++)*r; sum[0][5] += *(kv++)*r; - } - if (func[1]) { // geometric 1/r^6 - register double r = cek->re*cek->re+cek->im*cek->im; ++cek; - sum[1][0] += *(kv++)*r; sum[1][1] += *(kv++)*r; sum[1][2] += *(kv++)*r; - sum[1][3] += *(kv++)*r; sum[1][4] += *(kv++)*r; sum[1][5] += *(kv++)*r; - } - if (func[2]) { // arithmetic 1/r^6 - register double r = - (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+ - (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+ - (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+ - 0.5*(cek[3].re*cek[3].re+cek[3].im*cek[3].im); cek += 7; - sum[2][0] += *(kv++)*r; sum[2][1] += *(kv++)*r; sum[2][2] += *(kv++)*r; - sum[2][3] += *(kv++)*r; sum[2][4] += *(kv++)*r; sum[2][5] += *(kv++)*r; - } - if (func[3]) { - register double r = cek->re*cek->re+cek->im*cek->im; ++cek; - sum[3][0] += *(kv++)*r; sum[3][1] += *(kv++)*r; sum[3][2] += *(kv++)*r; - sum[3][3] += *(kv++)*r; sum[3][4] += *(kv++)*r; sum[3][5] += *(kv++)*r; - } - } - for (int k=0; kq; - double *x = atom->x[0]-1, *xn = x+3*atom->nlocal-1; - double dipole = 0.0, dipole_all; - - while ((x+=3)qqrd2e * scale; - - double ffact = -4.0*MY_PI*qscale*dipole_all/volume; - double *f = atom->f[0]-1, *fn = f+3*atom->nlocal-3; - - q = atom->q; - while ((f+=3)all(FLERR,KSPACE_ILLEGAL); + precision = fabs(atof(arg[0])); + memset(function, 0, EWALD_NORDER*sizeof(int)); + kenergy = kvirial = NULL; + cek_local = cek_global = NULL; + ekr_local = NULL; + hvec = NULL; + kvec = NULL; + B = NULL; + first_output = 0; + energy_self_peratom = NULL; + virial_self_peratom = NULL; +} + +EwaldN::~EwaldN() +{ + deallocate(); + deallocate_peratom(); + delete [] ekr_local; + delete [] B; +} + +/* --------------------------------------------------------------------- */ + +void EwaldN::init() +{ + nkvec = nkvec_max = nevec = nevec_max = 0; + nfunctions = nsums = sums = 0; + nbox = -1; + bytes = 0.0; + + if (!comm->me) { // output message + if (screen) fprintf(screen,"EwaldN initialization ...\n"); + if (logfile) fprintf(logfile,"EwaldN initialization ...\n"); + } + + if (domain->dimension == 2) // check for errors + error->all(FLERR,"Cannot use EwaldN with 2d simulation"); + if (slabflag == 0 && domain->nonperiodic > 0) + error->all(FLERR,"Cannot use nonperiodic boundaries with EwaldN"); + if (slabflag == 1) { + if (domain->xperiodic != 1 || domain->yperiodic != 1 || + domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) + error->all(FLERR,"Incorrect boundaries with slab EwaldN"); + } + + scale = 1.0; + //mumurd2e = force->mumurd2e; + //dielectric = force->dielectric; + mumurd2e = dielectric = 1.0; + + int tmp; + Pair *pair = force->pair; + int *ptr = pair ? (int *) pair->extract("ewald_order",tmp) : NULL; + double *cutoff = pair ? (double *) pair->extract("cut_coul",tmp) : NULL; + if (!(ptr||cutoff)) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + int ewald_order = ptr ? *((int *) ptr) : 1<<1; + int ewald_mix = ptr ? *((int *) pair->extract("ewald_mix",tmp)) : GEOMETRIC; + memset(function, 0, EWALD_NFUNCS*sizeof(int)); + for (int i=0; i<=EWALD_NORDER; ++i) // transcribe order + if (ewald_order&(1<pair_style); + error->all(FLERR,str); + default: + sprintf(str, "%s pair_style %s", KSPACE_ORDER, force->pair_style); + error->all(FLERR,str); + } + nfunctions += function[k] = 1; + nsums += n[k]; + } + + g_ewald = (1.35 - 0.15*log(precision))/ *cutoff; // determine resolution + g2_max = -4.0*g_ewald*g_ewald*log(precision); + + if (!comm->me) { // output results + if (screen) fprintf(screen, " G vector = %g\n", g_ewald); + if (logfile) fprintf(logfile, " G vector = %g\n", g_ewald); + } + + deallocate_peratom(); + peratom_allocate_flag = 0; +} + + +/* ---------------------------------------------------------------------- + adjust EwaldN coeffs, called initially and whenever volume has changed +------------------------------------------------------------------------- */ + +void EwaldN::setup() +{ + volume = shape_det(domain->h)*slab_volfactor; // cell volume + memcpy(unit, domain->h_inv, sizeof(shape)); // wave vector units + shape_scalar_mult(unit, 2.0*MY_PI); + unit[2] /= slab_volfactor; + + //int nbox_old = nbox, nkvec_old = nkvec; + if (precision>=1) nbox = 0; + else { + vector n = {1.0, 1.0, 1.0}; // based on cutoff + vec_scalar_mult(n, g_ewald*sqrt(-log(precision))/MY_PI); + shape_vec_dot(n, n, domain->h); + n[2] *= slab_volfactor; + nbox = (int) n[0]; + if (nbox<(int) n[1]) nbox = (int) n[1]; + if (nbox<(int) n[2]) nbox = (int) n[2]; + } + reallocate(); + + coefficients(); // compute coeffs + init_coeffs(); + init_coeff_sums(); + init_self(); + + if (!(first_output||comm->me)) { // output on first + first_output = 1; + if (screen) fprintf(screen, + " vectors: nbox = %d, nkvec = %d\n", nbox, nkvec); + if (logfile) fprintf(logfile, + " vectors: nbox = %d, nkvec = %d\n", nbox, nkvec); + } +} + + +void EwaldN::reallocate() // allocate memory +{ + int ix, iy, iz; + int nkvec_max = nkvec; + vector h; + + nkvec = 0; // determine size(kvec) + int *kflag = new int[(nbox+1)*(2*nbox+1)*(2*nbox+1)]; + int *flag = kflag; + + for (ix=0; ix<=nbox; ++ix) + for (iy=-nbox; iy<=nbox; ++iy) + for (iz=-nbox; iz<=nbox; ++iz) + if (!(ix||iy||iz)) *(flag++) = 0; + else if ((!ix)&&(iy<0)) *(flag++) = 0; + else if ((!(ix||iy))&&(iz<0)) *(flag++) = 0; // use symmetry + else { + h[0] = unit[0]*ix; + h[1] = unit[5]*ix+unit[1]*iy; + h[2] = unit[4]*ix+unit[3]*iy+unit[2]*iz; + if ((*(flag++) = h[0]*h[0]+h[1]*h[1]+h[2]*h[2]<=g2_max)) ++nkvec; + } + + if (nkvec>nkvec_max) { + deallocate(); // free memory + hvec = new hvector[nkvec]; // hvec + bytes += (nkvec-nkvec_max)*sizeof(hvector); + kvec = new kvector[nkvec]; // kvec + bytes += (nkvec-nkvec_max)*sizeof(kvector); + kenergy = new double[nkvec*nfunctions]; // kenergy + bytes += (nkvec-nkvec_max)*nfunctions*sizeof(double); + kvirial = new double[6*nkvec*nfunctions]; // kvirial + bytes += 6*(nkvec-nkvec_max)*nfunctions*sizeof(double); + cek_local = new complex[nkvec*nsums]; // cek_local + bytes += (nkvec-nkvec_max)*nsums*sizeof(complex); + cek_global = new complex[nkvec*nsums]; // cek_global + bytes += (nkvec-nkvec_max)*nsums*sizeof(complex); + nkvec_max = nkvec; + } + + flag = kflag; // create index and + kvector *k = kvec; // wave vectors + hvector *hi = hvec; + for (ix=0; ix<=nbox; ++ix) + for (iy=-nbox; iy<=nbox; ++iy) + for (iz=-nbox; iz<=nbox; ++iz) + if (*(flag++)) { + hi->x = unit[0]*ix; + hi->y = unit[5]*ix+unit[1]*iy; + (hi++)->z = unit[4]*ix+unit[3]*iy+unit[2]*iz; + k->x = ix+nbox; k->y = iy+nbox; (k++)->z = iz+nbox; } + + delete [] kflag; +} + +void EwaldN::reallocate_atoms() +{ + if (eflag_atom || vflag_atom) + if (atom->nlocal > atom->nmax) { + deallocate_peratom(); + allocate_peratom(); + } + + if ((nevec = atom->nmax*(2*nbox+1))<=nevec_max) return; + delete [] ekr_local; + ekr_local = new cvector[nevec]; + bytes += (nevec-nevec_max)*sizeof(cvector); + nevec_max = nevec; +} + +void EwaldN::allocate_peratom() +{ +memory->create(energy_self_peratom,atom->nmax,EWALD_NFUNCS,"ewald/n:energy_self_peratom"); +memory->create(virial_self_peratom,atom->nmax,EWALD_NFUNCS,"ewald/n:virial_self_peratom"); +} + + +void EwaldN::deallocate_peratom() // free memory +{ +memory->destroy(energy_self_peratom); +memory->destroy(virial_self_peratom); +} + + +void EwaldN::deallocate() // free memory +{ + delete [] hvec; hvec = NULL; + delete [] kvec; kvec = NULL; + delete [] kenergy; kenergy = NULL; + delete [] kvirial; kvirial = NULL; + delete [] cek_local; cek_local = NULL; + delete [] cek_global; cek_global = NULL; +} + + +void EwaldN::coefficients() // set up pre-factors +{ + vector h; + hvector *hi = hvec, *nh; + double eta2 = 0.25/(g_ewald*g_ewald); + double b1, b2, expb2, h1, h2, c1, c2; + double *ke = kenergy, *kv = kvirial; + int func0 = function[0], func12 = function[1]||function[2], + func3 = function[3]; + + for (nh = (hi = hvec)+nkvec; hintypes; + + if (function[1]) { // geometric 1/r^6 + double **b = (double **) force->pair->extract("B",tmp); + delete [] B; + B = new double[n+1]; + bytes += (n+1)*sizeof(double); + for (int i=0; i<=n; ++i) B[i] = sqrt(fabs(b[i][i])); + } + if (function[2]) { // arithmetic 1/r^6 + double **epsilon = (double **) force->pair->extract("epsilon",tmp); + double **sigma = (double **) force->pair->extract("sigma",tmp); + if (!(epsilon&&sigma)) + error->all(FLERR,"epsilon or sigma reference not set by pair style in ewald/n"); + double eps_i, sigma_i, sigma_n, *bi = B = new double[7*n+7]; + double c[7] = { + 1.0, sqrt(6.0), sqrt(15.0), sqrt(20.0), sqrt(15.0), sqrt(6.0), 1.0}; + for (int i=0; i<=n; ++i) { + eps_i = sqrt(epsilon[i][i]); + sigma_i = sigma[i][i]; + sigma_n = 1.0; + for (int j=0; j<7; ++j) { + *(bi++) = sigma_n*eps_i*c[j]; sigma_n *= sigma_i; + } + } + } +} + + +void EwaldN::init_coeff_sums() // sums based on atoms +{ + if (sums) return; // calculated only once + sums = 1; + + Sum sum_local[EWALD_MAX_NSUMS]; + + memset(sum_local, 0, EWALD_MAX_NSUMS*sizeof(Sum)); + if (function[0]) { // 1/r + double *q = atom->q, *qn = q+atom->nlocal; + for (double *i=q; itype, *ntype = type+atom->nlocal; + for (int *i=type; itype, *ntype = type+atom->nlocal; + for (int *i=type; imu; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) + sum_local[9].x2 += mu[i][3]*mu[i][3]; + } + MPI_Allreduce(sum_local, sum, 2*EWALD_MAX_NSUMS, MPI_DOUBLE, MPI_SUM, world); +} + + +void EwaldN::init_self() +{ + double g1 = g_ewald, g2 = g1*g1, g3 = g1*g2; + + memset(energy_self, 0, EWALD_NFUNCS*sizeof(double)); // self energy + memset(virial_self, 0, EWALD_NFUNCS*sizeof(double)); + const double qscale = force->qqrd2e * scale; + + if (function[0]) { // 1/r + virial_self[0] = -0.5*MY_PI*qscale/(g2*volume)*sum[0].x*sum[0].x; + energy_self[0] = sum[0].x2*qscale*g1/sqrt(MY_PI)-virial_self[0]; + } + if (function[1]) { // geometric 1/r^6 + virial_self[1] = MY_PI*sqrt(MY_PI)*g3/(6.0*volume)*sum[1].x*sum[1].x; + energy_self[1] = -sum[1].x2*g3*g3/12.0+virial_self[1]; + } + if (function[2]) { // arithmetic 1/r^6 + virial_self[2] = MY_PI*sqrt(MY_PI)*g3/(48.0*volume)*(sum[2].x*sum[8].x+ + sum[3].x*sum[7].x+sum[4].x*sum[6].x+0.5*sum[5].x*sum[5].x); + energy_self[2] = -sum[2].x2*g3*g3/3.0+virial_self[2]; + } + if (function[3]) { // dipole + virial_self[3] = 0; // in surface + energy_self[3] = sum[9].x2*mumurd2e*2.0*g3/3.0/sqrt(MY_PI)-virial_self[3]; + } +} + + +void EwaldN::init_self_peratom() +{ + if (!(vflag_atom || eflag_atom)) return; + + double g1 = g_ewald, g2 = g1*g1, g3 = g1*g2; + const double qscale = force->qqrd2e * scale; + int nlocal = atom->nlocal; + + for (int k = 0; k < 3; k++) + for (int i = 0; i < nlocal; i++) { + energy_self_peratom[i][k] = 0; + virial_self_peratom[i][k] = 0; + } + + if (function[0]) { // 1/r + double *q = atom->q; + for (int i = 0; i < nlocal; i++) { + virial_self_peratom[i][0] = -0.5*MY_PI*qscale/(g2*volume)*q[i]*sum[0].x; + energy_self_peratom[i][0] = q[i]*q[i]*qscale*g1/sqrt(MY_PI)-virial_self_peratom[i][0]; + } + } + if (function[1]) { // geometric 1/r^6 + int *type = atom->type; + for (int i = 0; i < nlocal; i++) { + virial_self_peratom[i][1] = MY_PI*sqrt(MY_PI)*g3/(6.0*volume)*B[type[i]]*sum[1].x; + energy_self_peratom[i][1] = -B[type[i]]*B[type[i]]*g3*g3/12.0+virial_self_peratom[i][1]; + } + } + if (function[2]) { // arithmetic 1/r^6 + double *bi; + int *type = atom->type; + for (int i = 0; i < nlocal; i++) { + bi = B+7*type[i]+7; + for (int k=2; k<9; ++k) { + virial_self_peratom[i][2] += 0.5*MY_PI*sqrt(MY_PI)*g3/(48.0*volume)*sum[k].x*(--bi)[0]; + } + energy_self_peratom[i][2] = -bi[0]*bi[6]*g3*g3/3.0+virial_self_peratom[i][2]; + } + } + if (function[3]) { // dipole + double **mu = atom->mu; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) { + virial_self_peratom[i][3] = 0; // in surface + energy_self_peratom[i][3] = mu[i][3]*mu[i][3]*mumurd2e*2.0*g3/3.0/sqrt(MY_PI)-virial_self_peratom[i][3]; + } + } +} + +/* ---------------------------------------------------------------------- + compute the EwaldN long-range force, energy, virial +------------------------------------------------------------------------- */ + +void EwaldN::compute(int eflag, int vflag) +{ + if (!nbox) return; + + // set energy/virial flags + // invoke allocate_peratom() if needed for first time + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = eflag_global = vflag_global = eflag_atom = vflag_atom = 0; + + if (!peratom_allocate_flag && (eflag_atom || vflag_atom)) { + allocate_peratom(); + peratom_allocate_flag = 1; + } + + reallocate_atoms(); + init_self_peratom(); + compute_ek(); + compute_force(); + compute_surface(); + compute_surface_peratom(); + compute_energy(); + compute_energy_peratom(); + compute_virial(); + compute_virial_peratom(); +} + + +void EwaldN::compute_ek() +{ + cvector *ekr = ekr_local; + int lbytes = (2*nbox+1)*sizeof(cvector); + hvector *h; + kvector *k, *nk = kvec+nkvec; + cvector *z = new cvector[2*nbox+1]; + cvector z1, *zx, *zy, *zz, *zn = z+2*nbox; + complex *cek, zxyz, zxy, cx; + vector mui; + double *x = atom->x[0], *xn = x+3*atom->nlocal, *q = atom->q, qi, bi, ci[7]; + double *mu = atom->mu ? atom->mu[0] : NULL; + int i, kx, ky, n = nkvec*nsums, *type = atom->type, tri = domain->triclinic; + int func[EWALD_NFUNCS]; + + memcpy(func, function, EWALD_NFUNCS*sizeof(int)); + memset(cek_local, 0, n*sizeof(complex)); // reset sums + while (xx, 1, 0); C_SET(zz->y, 1, 0); C_SET(zz->z, 1, 0); // z[0] + if (tri) { // triclinic z[1] + C_ANGLE(z1.x, unit[0]*x[0]+unit[5]*x[1]+unit[4]*x[2]); + C_ANGLE(z1.y, unit[1]*x[1]+unit[3]*x[2]); + C_ANGLE(z1.z, x[2]*unit[2]); x += 3; + } + else { // orthogonal z[1] + C_ANGLE(z1.x, *(x++)*unit[0]); + C_ANGLE(z1.y, *(x++)*unit[1]); + C_ANGLE(z1.z, *(x++)*unit[2]); + } + for (; zzx, zz->x, z1.x); // 3D k-vector + C_RMULT(zy->y, zz->y, z1.y); C_CONJ(zx->y, zy->y); + C_RMULT(zy->z, zz->z, z1.z); C_CONJ(zx->z, zy->z); + } + kx = ky = -1; + cek = cek_local; + if (func[0]) qi = *(q++); + if (func[1]) bi = B[*type]; + if (func[2]) memcpy(ci, B+7*type[0], 7*sizeof(double)); + if (func[3]) { + memcpy(mui, mu, sizeof(vector)); mu += 3; + vec_scalar_mult(mui, mu[3]); + mu += 4; + h = hvec; + } + for (k=kvec; ky) { // based on order in + if (kx!=k->x) cx = z[kx = k->x].x; // reallocate + C_RMULT(zxy, z[ky = k->y].y, cx); + } + C_RMULT(zxyz, z[k->z].z, zxy); + if (func[0]) { + cek->re += zxyz.re*qi; (cek++)->im += zxyz.im*qi; + } + if (func[1]) { + cek->re += zxyz.re*bi; (cek++)->im += zxyz.im*bi; + } + if (func[2]) for (i=0; i<7; ++i) { + cek->re += zxyz.re*ci[i]; (cek++)->im += zxyz.im*ci[i]; + } + if (func[3]) { + register double muk = mui[0]*h->x+mui[1]*h->y+mui[2]*h->z; ++h; + cek->re += zxyz.re*muk; (cek++)->im += zxyz.im*muk; + } + } + ekr = (cvector *) ((char *) memcpy(ekr, z, lbytes)+lbytes); + ++type; + } + MPI_Allreduce(cek_local, cek_global, 2*n, MPI_DOUBLE, MPI_SUM, world); + + delete [] z; +} + +void EwaldN::compute_force() +{ + kvector *k; + hvector *h, *nh; + cvector *z = ekr_local; + vector sum[EWALD_MAX_NSUMS], mui; + complex *cek, zc, zx, zxy; + double *f = atom->f[0], *fn = f+3*atom->nlocal, *q = atom->q, *t = NULL; + double *mu = atom->mu ? atom->mu[0] : NULL; + const double qscale = force->qqrd2e * scale; + double *ke, c[EWALD_NFUNCS] = { + 8.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(12.0*volume), + 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 8.0*MY_PI*mumurd2e/volume}; + double kt = 4.0*pow(g_ewald, 3.0)/3.0/sqrt(MY_PI)/c[3]; + int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type; + int func[EWALD_NFUNCS]; + + if (atom->torque) t = atom->torque[0]; + memcpy(func, function, EWALD_NFUNCS*sizeof(int)); + memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector)); // fj = -dE/dr = + for (; fy) { // based on order in + if (kx!=k->x) zx = z[kx = k->x].x; // reallocate + C_RMULT(zxy, z[ky = k->y].y, zx); + } + C_CRMULT(zc, z[k->z].z, zxy); + if (func[0]) { // 1/r + register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek; + sum[0][0] += h->x*im; sum[0][1] += h->y*im; sum[0][2] += h->z*im; + } + if (func[1]) { // geometric 1/r^6 + register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek; + sum[1][0] += h->x*im; sum[1][1] += h->y*im; sum[1][2] += h->z*im; + } + if (func[2]) { // arithmetic 1/r^6 + register double im, c = *(ke++); + for (i=2; i<9; ++i) { + im = c*(zc.im*cek->re+cek->im*zc.re); ++cek; + sum[i][0] += h->x*im; sum[i][1] += h->y*im; sum[i][2] += h->z*im; + } + } + if (func[3]) { // dipole + register double im = *(ke++)*(zc.im*cek->re+ + cek->im*zc.re)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z); ++cek; + sum[9][0] += h->x*im; sum[9][1] += h->y*im; sum[9][2] += h->z*im; + } + } + if (func[0]) { // 1/r + register double qi = *(q++)*c[0]; + f[0] -= sum[0][0]*qi; f[1] -= sum[0][1]*qi; f[2] -= sum[0][2]*qi; + } + if (func[1]) { // geometric 1/r^6 + register double bi = B[*type]*c[1]; + f[0] -= sum[1][0]*bi; f[1] -= sum[1][1]*bi; f[2] -= sum[1][2]*bi; + } + if (func[2]) { // arithmetic 1/r^6 + register double *bi = B+7*type[0]+7; + for (i=2; i<9; ++i) { + register double c2 = (--bi)[0]*c[2]; + f[0] -= sum[i][0]*c2; f[1] -= sum[i][1]*c2; f[2] -= sum[i][2]*c2; + } + } + if (func[3]) { // dipole + f[0] -= sum[9][0]; f[1] -= sum[9][1]; f[2] -= sum[9][2]; + *(t++) -= mui[1]*sum[0][2]+mui[2]*sum[0][1]-mui[0]*kt; // torque + *(t++) -= mui[2]*sum[0][0]+mui[0]*sum[0][2]-mui[1]*kt; + *(t++) -= mui[0]*sum[0][1]+mui[1]*sum[0][0]-mui[2]*kt; + } + z = (cvector *) ((char *) z+lbytes); + ++type; + } +} + +void EwaldN::compute_surface() +{ + if (!function[3]) return; + + vector sum_local = VECTOR_NULL, sum_total; + double **mu = atom->mu; + int nlocal = atom->nlocal; + + for (int i=0; i < nlocal; i++) { + register double di = mu[i][3]; + sum_local[0] += di*mu[i][0]; + sum_local[1] += di*mu[i][1]; + sum_local[2] += di*mu[i][2]; + } + MPI_Allreduce(sum_local, sum_total, 3, MPI_DOUBLE, MPI_SUM, world); + + energy_self[3] += virial_self[3]; + virial_self[3] = + mumurd2e*(2.0*MY_PI*vec_dot(sum_total,sum_total)/(2.0*dielectric+1)/volume); + energy_self[3] -= virial_self[3]; +} + +void EwaldN::compute_surface_peratom() +{ + if (!(vflag_atom || eflag_atom)) return; + if (!function[3]) return; + + vector sum_local = VECTOR_NULL, sum_total; + double **mu = atom->mu; + int nlocal = atom->nlocal; + + for (int i=0; i < nlocal; i++) { + register double di = mu[i][3]; + sum_local[0] += di*mu[i][0]; + sum_local[1] += di*mu[i][1]; + sum_local[2] += di*mu[i][2]; + } + + MPI_Allreduce(sum_local, sum_total, 3, MPI_DOUBLE, MPI_SUM, world); + + + for (int i = 0; i < nlocal; i++) { + energy_self_peratom[i][3] += virial_self_peratom[i][3]; + register double di = mu[i][3]; + virial_self_peratom[i][3] = + mumurd2e*(2.0*MY_PI*(di*mu[i][0]*sum_total[0] + di*mu[i][1]*sum_total[1] + + di*mu[i][2]*sum_total[2])/(2.0*dielectric+1)/volume); + energy_self_peratom[i][3] -= virial_self_peratom[i][3]; + } + +} + +void EwaldN::compute_energy() +{ + energy = 0.0; + if (!eflag_global) return; + + complex *cek = cek_global; + double *ke = kenergy; + const double qscale = force->qqrd2e * scale; + double c[EWALD_NFUNCS] = { + 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), + 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; + double sum[EWALD_NFUNCS]; + int func[EWALD_NFUNCS]; + + memcpy(func, function, EWALD_NFUNCS*sizeof(int)); + memset(sum, 0, EWALD_NFUNCS*sizeof(double)); // reset sums + for (int k=0; kre*cek->re+cek->im*cek->im); ++cek; } + if (func[1]) { // geometric 1/r^6 + sum[1] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; } + if (func[2]) { // arithmetic 1/r^6 + register double r = + (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+ + (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+ + (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+ + 0.5*(cek[3].re*cek[3].re+cek[3].im*cek[3].im); cek += 7; + sum[2] += *(ke++)*r; + } + if (func[3]) { // dipole + sum[3] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; } + } + for (int k=0; kq; + double *mu = atom->mu ? atom->mu[0] : NULL; + const double qscale = force->qqrd2e * scale; + double *ke = kenergy; + double c[EWALD_NFUNCS] = { + 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), + 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; + int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type; + int func[EWALD_NFUNCS]; + + memcpy(func, function, EWALD_NFUNCS*sizeof(int)); + for (int j = 0; j < atom->nlocal; j++) { + k = kvec; + kx = ky = -1; + ke = kenergy; + cek = cek_global; + memset(sum, 0, EWALD_MAX_NSUMS*sizeof(double)); + if (func[3]) { + register double di = mu[3] * c[3]; + mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[1]; mui[2] = di*(mu++)[2]; + mu++; + } + for (nh = (h = hvec)+nkvec; hy) { // based on order in + if (kx!=k->x) zx = z[kx = k->x].x; // reallocate + C_RMULT(zxy, z[ky = k->y].y, zx); + } + C_CRMULT(zc, z[k->z].z, zxy); + if (func[0]) { // 1/r + sum[0] += *(ke++)*(cek->re*zc.re - cek->im*zc.im); ++cek; } + if (func[1]) { // geometric 1/r^6 + sum[1] += *(ke++)*(cek->re*zc.re - cek->im*zc.im); ++cek; } + if (func[2]) { // arithmetic 1/r^6 + register double im, c = *(ke++); + for (i=2; i<9; ++i) { + im = c*(cek->re*zc.re - cek->im*zc.im); ++cek; + sum[i] += im; + } + } + if (func[3]) { // dipole + sum[9] += *(ke++)*(cek->re*zc.re - + cek->im*zc.im)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z); ++cek; + } + } + + if (func[0]) { // 1/r + register double qj = *(q++)*c[0]; + eatom[j] += sum[0]*qj - energy_self_peratom[j][0]; + } + if (func[1]) { // geometric 1/r^6 + register double bj = B[*type]*c[1]; + eatom[j] += sum[1]*bj - energy_self_peratom[j][1]; + } + if (func[2]) { // arithmetic 1/r^6 + register double *bj = B+7*type[0]+7; + for (i=2; i<9; ++i) { + register double c2 = (--bj)[0]*c[2]; + eatom[j] += 0.5*sum[i]*c2; + } + eatom[j] -= energy_self_peratom[j][2]; + } + if (func[3]) { // dipole + eatom[j] += sum[9] - energy_self_peratom[j][3]; + } + z = (cvector *) ((char *) z+lbytes); + ++type; + } +} + +#define swap(a, b) { register double t = a; a= b; b = t; } + +void EwaldN::compute_virial() +{ + memset(virial, 0, sizeof(shape)); + if (!vflag_global) return; + + complex *cek = cek_global; + double *kv = kvirial; + const double qscale = force->qqrd2e * scale; + double c[EWALD_NFUNCS] = { + 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), + 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; + shape sum[EWALD_NFUNCS]; + int func[EWALD_NFUNCS]; + + memcpy(func, function, EWALD_NFUNCS*sizeof(int)); + memset(sum, 0, EWALD_NFUNCS*sizeof(shape)); + for (int k=0; kre*cek->re+cek->im*cek->im; ++cek; + sum[0][0] += *(kv++)*r; sum[0][1] += *(kv++)*r; sum[0][2] += *(kv++)*r; + sum[0][3] += *(kv++)*r; sum[0][4] += *(kv++)*r; sum[0][5] += *(kv++)*r; + } + if (func[1]) { // geometric 1/r^6 + register double r = cek->re*cek->re+cek->im*cek->im; ++cek; + sum[1][0] += *(kv++)*r; sum[1][1] += *(kv++)*r; sum[1][2] += *(kv++)*r; + sum[1][3] += *(kv++)*r; sum[1][4] += *(kv++)*r; sum[1][5] += *(kv++)*r; + } + if (func[2]) { // arithmetic 1/r^6 + register double r = + (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+ + (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+ + (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+ + 0.5*(cek[3].re*cek[3].re+cek[3].im*cek[3].im); cek += 7; + sum[2][0] += *(kv++)*r; sum[2][1] += *(kv++)*r; sum[2][2] += *(kv++)*r; + sum[2][3] += *(kv++)*r; sum[2][4] += *(kv++)*r; sum[2][5] += *(kv++)*r; + } + if (func[3]) { + register double r = cek->re*cek->re+cek->im*cek->im; ++cek; + sum[3][0] += *(kv++)*r; sum[3][1] += *(kv++)*r; sum[3][2] += *(kv++)*r; + sum[3][3] += *(kv++)*r; sum[3][4] += *(kv++)*r; sum[3][5] += *(kv++)*r; + } + } + for (int k=0; kq; + double *mu = atom->mu ? atom->mu[0] : NULL; + const double qscale = force->qqrd2e * scale; + double c[EWALD_NFUNCS] = { + 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), + 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; + shape sum[EWALD_MAX_NSUMS]; + int func[EWALD_NFUNCS]; + + memcpy(func, function, EWALD_NFUNCS*sizeof(int)); + int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type; + for (int j = 0; j < atom->nlocal; j++) { + k = kvec; + kx = ky = -1; + kv = kvirial; + cek = cek_global; + memset(sum, 0, EWALD_MAX_NSUMS*sizeof(shape)); + if (func[3]) { + register double di = mu[3] * c[3]; + mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[1]; mui[2] = di*(mu++)[2]; + mu++; + } + for (nh = (h = hvec)+nkvec; hy) { // based on order in + if (kx!=k->x) zx = z[kx = k->x].x; // reallocate + C_RMULT(zxy, z[ky = k->y].y, zx); + } + C_CRMULT(zc, z[k->z].z, zxy); + if (func[0]) { // 1/r + register double r = cek->re*zc.re - cek->im*zc.im; ++cek; + sum[0][0] += *(kv++)*r; sum[0][1] += *(kv++)*r; sum[0][2] += *(kv++)*r; + sum[0][3] += *(kv++)*r; sum[0][4] += *(kv++)*r; sum[0][5] += *(kv++)*r; + } + if (func[1]) { // geometric 1/r^6 + register double r = cek->re*zc.re - cek->im*zc.im; ++cek; + sum[1][0] += *(kv++)*r; sum[1][1] += *(kv++)*r; sum[1][2] += *(kv++)*r; + sum[1][3] += *(kv++)*r; sum[1][4] += *(kv++)*r; sum[1][5] += *(kv++)*r; + } + if (func[2]) { // arithmetic 1/r^6 + register double r; + for (i=2; i<9; ++i) { + r = cek->re*zc.re - cek->im*zc.im; ++cek; + sum[i][0] += *(kv++)*r; sum[i][1] += *(kv++)*r; sum[i][2] += *(kv++)*r; + sum[i][3] += *(kv++)*r; sum[i][4] += *(kv++)*r; sum[i][5] += *(kv++)*r; + kv -= 6; + } + kv += 6; + } + if (func[3]) { // dipole + register double r = (cek->re*zc.re - + cek->im*zc.im)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z); ++cek; + sum[9][0] += *(kv++)*r; sum[9][1] += *(kv++)*r; sum[9][2] += *(kv++)*r; + sum[9][3] += *(kv++)*r; sum[9][4] += *(kv++)*r; sum[9][5] += *(kv++)*r; + } + } + + if (func[0]) { // 1/r + register double qi = *(q++)*c[0]; + for (int n = 0; n < 6; n++) vatom[j][n] += sum[0][n]*qi; + } + if (func[1]) { // geometric 1/r^6 + register double bi = B[*type]*c[1]; + for (int n = 0; n < 6; n++) vatom[j][n] += sum[1][n]*bi; + double dum = sum[1][2]; + } + if (func[2]) { // arithmetic 1/r^6 + register double *bj = B+7*type[0]+7; + for (i=2; i<9; ++i) { + register double c2 = (--bj)[0]*c[2]; + for (int n = 0; n < 6; n++) { + vatom[j][n] += 0.5*sum[i][n]*c2; + } + } + } + if (func[3]) { // dipole + for (int n = 0; n < 6; n++) vatom[j][n] += sum[9][n]; + } + + for (int k=0; kq; + 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 + + const double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; + const double qscale = force->qqrd2e * scale; + + if (eflag_global) energy += qscale * e_slabcorr; + + // per-atom energy + + if (eflag_atom) { + double efact = 2.0*MY_PI*dipole_all/volume; + for (int i = 0; i < nlocal; i++) eatom[i] += qscale * q[i]*x[i][2]*efact; + } + + // add on force corrections + + double ffact = -4.0*MY_PI*dipole_all/volume; + double **f = atom->f; + + for (int i = 0; i < nlocal; i++) f[i][2] += qscale * q[i]*ffact; +} + diff --git a/src/USER-EWALDN/ewald_n.h b/src/USER-EWALDN/ewald_n.h index 2bec485323..f7ba0a692d 100644 --- a/src/USER-EWALDN/ewald_n.h +++ b/src/USER-EWALDN/ewald_n.h @@ -1,82 +1,91 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#ifdef KSPACE_CLASS - -KSpaceStyle(ewald/n,EwaldN) - -#else - -#ifndef LMP_EWALD_N_H -#define LMP_EWALD_N_H - -#include "kspace.h" -#include "math_complex.h" - -namespace LAMMPS_NS { - -#define EWALD_NORDER 6 -#define EWALD_NFUNCS 4 -#define EWALD_MAX_NSUMS 10 -#define EWALD_NSUMS {1, 1, 7, 1} - -typedef struct cvector { complex x, y, z; } cvector; -typedef struct hvector { double x, y, z; } hvector; -typedef struct kvector { long x, y, z; } kvector; - -class EwaldN : public KSpace { - public: - EwaldN(class LAMMPS *, int, char **); - ~EwaldN(); - void init(); - void setup(); - void compute(int, int); - double memory_usage() {return bytes;} - - private: - double unit[6]; - int function[EWALD_NFUNCS], first_output; - - int nkvec, nkvec_max, nevec, nevec_max, - nbox, nfunctions, nsums, sums; - double bytes; - double precision, g2_max; - double *kenergy, energy_self[EWALD_NFUNCS]; - double *kvirial, virial_self[EWALD_NFUNCS]; - cvector *ekr_local; - hvector *hvec; - kvector *kvec; - - double mumurd2e, dielectric, *B, volume; - struct Sum { double x, x2; } sum[EWALD_MAX_NSUMS]; - complex *cek_local, *cek_global; - - void reallocate(); - void reallocate_atoms(); - void deallocate(); - void coefficients(); - void init_coeffs(); - void init_coeff_sums(); - void init_self(); - void compute_ek(); - void compute_force(); - void compute_surface(); - void compute_energy(int); - void compute_virial(int); - void compute_slabcorr(int); -}; - -} - -#endif -#endif +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef KSPACE_CLASS + +KSpaceStyle(ewald/n,EwaldN) + +#else + +#ifndef LMP_EWALD_N_H +#define LMP_EWALD_N_H + +#include "kspace.h" +#include "math_complex.h" + +namespace LAMMPS_NS { + +#define EWALD_NORDER 6 +#define EWALD_NFUNCS 4 +#define EWALD_MAX_NSUMS 10 +#define EWALD_NSUMS {1, 1, 7, 1} + +typedef struct cvector { complex x, y, z; } cvector; +typedef struct hvector { double x, y, z; } hvector; +typedef struct kvector { long x, y, z; } kvector; + +class EwaldN : public KSpace { + public: + EwaldN(class LAMMPS *, int, char **); + ~EwaldN(); + void init(); + void setup(); + void compute(int, int); + double memory_usage() {return bytes;} + + private: + double unit[6]; + int function[EWALD_NFUNCS], first_output; + + int nkvec, nkvec_max, nevec, nevec_max, + nbox, nfunctions, nsums, sums; + int peratom_allocate_flag; + double bytes; + double precision, g2_max; + double *kenergy, energy_self[EWALD_NFUNCS]; + double *kvirial, virial_self[EWALD_NFUNCS]; + double **energy_self_peratom; + double **virial_self_peratom; + cvector *ekr_local; + hvector *hvec; + kvector *kvec; + + double mumurd2e, dielectric, *B, volume; + struct Sum { double x, x2; } sum[EWALD_MAX_NSUMS]; + complex *cek_local, *cek_global; + + void reallocate(); + void allocate_peratom(); + void reallocate_atoms(); + void deallocate(); + void deallocate_peratom(); + void coefficients(); + void init_coeffs(); + void init_coeff_sums(); + void init_self(); + void init_self_peratom(); + void compute_ek(); + void compute_force(); + void compute_surface(); + void compute_surface_peratom(); + void compute_energy(); + void compute_energy_peratom(); + void compute_virial(); + void compute_virial_peratom(); + void compute_slabcorr(); +}; + +} + +#endif +#endif diff --git a/src/fix_box_relax.cpp b/src/fix_box_relax.cpp index a372841361..24891ee950 100644 --- a/src/fix_box_relax.cpp +++ b/src/fix_box_relax.cpp @@ -324,7 +324,8 @@ void FixBoxRelax::init() temperature = modify->compute[icompute]; icompute = modify->find_compute(id_press); - if (icompute < 0) error->all(FLERR,"Pressure ID for fix box/relax does not exist"); + if (icompute < 0) + error->all(FLERR,"Pressure ID for fix box/relax does not exist"); pressure = modify->compute[icompute]; pv2e = 1.0 / force->nktv2p; diff --git a/src/kspace.cpp b/src/kspace.cpp index 7db770c033..6ea1499b2f 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -16,6 +16,7 @@ #include "kspace.h" #include "atom.h" #include "comm.h" +#include "force.h" #include "memory.h" #include "error.h" @@ -32,6 +33,11 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) slabflag = 0; slab_volfactor = 1; + accuracy_absolute = -1.0; + two_charge_force = force->qqr2e * + (force->qelectron * force->qelectron) / + (force->angstrom * force->angstrom); + maxeatom = maxvatom = 0; eatom = NULL; vatom = NULL; @@ -121,6 +127,10 @@ void KSpace::modify_params(int narg, char **arg) if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); order = atoi(arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"force") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); + accuracy_absolute = atof(arg[iarg+1]); + iarg += 2; } else if (strcmp(arg[iarg],"gewald") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); g_ewald = atof(arg[iarg+1]); diff --git a/src/kspace.h b/src/kspace.h index 83333da746..3549a35cfa 100644 --- a/src/kspace.h +++ b/src/kspace.h @@ -45,6 +45,13 @@ class KSpace : protected Pointers { int slabflag; double scale; double slab_volfactor; + + double accuracy; // accuracy of KSpace solver (force units) + double accuracy_absolute; // user-specifed accuracy in force units + double accuracy_relative; // user-specified dimensionless accuracy + // accurary = acc_rel * two_charge_force + double two_charge_force; // force in user units of two point + // charges separated by 1 Angstrom int evflag,evflag_atom; int eflag_either,eflag_global,eflag_atom;