diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp new file mode 100644 index 0000000000..cd55f0fbb2 --- /dev/null +++ b/src/KSPACE/msm.cpp @@ -0,0 +1,1693 @@ +/* ---------------------------------------------------------------------- + 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: Paul Crozier, Stan Moore, Stephen Bond, (all SNL) +------------------------------------------------------------------------- */ + +#include "lmptype.h" +#include "mpi.h" +#include "string.h" +#include "stdio.h" +#include "stdlib.h" +#include "math.h" +#include "msm.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 "memory.h" +#include "error.h" +#include "modify.h" +#include "fix.h" + +#include "math_const.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define MAXORDER 7 +#define MAX_LEVELS 10 +#define OFFSET 16384 +#define SMALL 0.00001 +#define LARGE 10000.0 +#define EPS_HOC 1.0e-7 + +/* ---------------------------------------------------------------------- */ + +MSM::MSM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) +{ + if (narg < 1) error->all(FLERR,"Illegal kspace_style msm command"); + + accuracy_relative = atof(arg[0]); + + nfactors = 1; + factors = new int[nfactors]; + factors[0] = 2; + + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + + buf1 = buf2 = NULL; + + phi1d = dphi1d = NULL; + + nmax = 0; + part2grid = NULL; + + g_direct = NULL; + + levels = 0; + +} + +/* ---------------------------------------------------------------------- + free all memory +------------------------------------------------------------------------- */ + +MSM::~MSM() +{ + delete [] factors; + deallocate(); + memory->destroy(part2grid); + memory->destroy(g_direct); + deallocate_levels(); +} + +/* ---------------------------------------------------------------------- + called once before run +------------------------------------------------------------------------- */ + +void MSM::init() +{ + if (me == 0) { + if (screen) fprintf(screen,"MSM initialization ...\n"); + if (logfile) fprintf(logfile,"MSM initialization ...\n"); + } + + // error check + + if (domain->triclinic) + error->all(FLERR,"Cannot (yet) use MSM with triclinic box"); + + if (domain->dimension == 2) error->all(FLERR,"Cannot (yet) use MSM with 2d simulation"); + + if (!atom->q_flag) error->all(FLERR,"Kspace style requires atom attribute q"); + + if (slabflag == 1) + error->all(FLERR,"Cannot use slab correction with MSM"); + + if (domain->nonperiodic > 0) + error->all(FLERR,"Cannot (yet) use nonperiodic boundaries with MSM"); + + for (int i = 0; i < modify->nfix; i++) { + if ((strcmp(modify->fix[i]->style,"npt") == 0) || + (strcmp(modify->fix[i]->style,"nph") == 0)) { + error->all(FLERR,"Cannot (yet) use MSM with npt/nph simulation"); + } + } + + order = 4; // 4th order interpolation scheme has been implemented for MSM + + if (order > MAXORDER) { + char str[128]; + sprintf(str,"MSM order cannot be greater than %d",MAXORDER); + error->all(FLERR,str); + } + + // free all arrays previously allocated + + deallocate(); + + // extract short-range Coulombic cutoff from pair style + + qqrd2e = force->qqrd2e; + scale = 1.0; + + if (force->pair == NULL) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + int itmp; + double *p_cutoff = (double *) force->pair->extract("cut_msm",itmp); + if (p_cutoff == NULL) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + cutoff = *p_cutoff; + + if ((strcmp(force->kspace_style,"pppm/tip4p") == 0) || + (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0)) { + error->all(FLERR,"KSpace style is incompatible with Pair style"); + } + + // compute qsum & qsqsum and give error 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; + q2 = qsqsum * force->qqrd2e / force->dielectric; + + if (qsqsum == 0.0) + error->all(FLERR,"Cannot use kspace solver on system with no charge"); + if (fabs(qsum) > SMALL) { + char str[128]; + sprintf(str,"System is not charge neutral, net charge = %g",qsum); + error->all(FLERR,str); // Not yet sure of the correction needed for non-neutral systems + } + + // 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 grid resolution + + set_grid(); + + int flag_global = 0; + + // loop over grid levels + + for (int n=0; n= OFFSET || ny_msm[n] >= OFFSET || nz_msm[n] >= OFFSET) + error->all(FLERR,"MSM grid is too large"); + + // global indices of MSM grid range from 0 to N-1 + // nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of + // global MSM grid that I own without ghost cells + + if (n == 0) { + + nxlo_in_d[n] = static_cast (comm->xsplit[comm->myloc[0]] * nx_msm[n]); + nxhi_in_d[n] = static_cast (comm->xsplit[comm->myloc[0]+1] * nx_msm[n]) - 1; + + nylo_in_d[n] = static_cast (comm->ysplit[comm->myloc[1]] * ny_msm[n]); + nyhi_in_d[n] = static_cast (comm->ysplit[comm->myloc[1]+1] * ny_msm[n]) - 1; + + nzlo_in_d[n] = static_cast (comm->zsplit[comm->myloc[2]] * nz_msm[n]); + nzhi_in_d[n] = static_cast (comm->zsplit[comm->myloc[2]+1] * nz_msm[n]) - 1; + + } else { + + nxlo_in_d[n] = 0; + nxhi_in_d[n] = nx_msm[n] - 1; + + nylo_in_d[n] = 0; + nyhi_in_d[n] = ny_msm[n] - 1; + + nzlo_in_d[n] = 0; + nzhi_in_d[n] = nz_msm[n] - 1; + } + + + // Use simple method of parallel communication for now + + nxlo_in[n] = 0; + nxhi_in[n] = nx_msm[n] - 1; + + nylo_in[n] = 0; + nyhi_in[n] = ny_msm[n] - 1; + + nzlo_in[n] = 0; + nzhi_in[n] = nz_msm[n] - 1; + + // nlower,nupper = stencil size for mapping particles to MSM 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 + + // nlo_out,nhi_out = lower/upper limits of the 3d sub-brick of + // global MSM 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 + // dist[3] = particle position bound = subbox + skin/2.0 + // nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping + + double *prd,*sublo,*subhi; + + //prd = domain->prd; + //boxlo = domain->boxlo; + //sublo = domain->sublo; + //subhi = domain->subhi; + + // Use only one partition for now + + prd = domain->prd; + boxlo = domain->boxlo; + sublo = boxlo; + subhi = domain->boxhi; + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; + + double dist[3]; + double cuthalf = 0.0; + if (n == 0) cuthalf = 0.5*neighbor->skin; // Only applies to finest grid + dist[0] = dist[1] = dist[2] = cuthalf; + + int nlo,nhi; + + nlo = static_cast ((sublo[0]-dist[0]-boxlo[0]) * + nx_msm[n]/xprd + OFFSET) - OFFSET; + nhi = static_cast ((subhi[0]+dist[0]-boxlo[0]) * + nx_msm[n]/xprd + OFFSET) - OFFSET; + nxlo_out[n] = nlo + nlower; + nxhi_out[n] = nhi + nupper; + + nlo = static_cast ((sublo[1]-dist[1]-boxlo[1]) * + ny_msm[n]/yprd + OFFSET) - OFFSET; + nhi = static_cast ((subhi[1]+dist[1]-boxlo[1]) * + ny_msm[n]/yprd + OFFSET) - OFFSET; + nylo_out[n] = nlo + nlower; + nyhi_out[n] = nhi + nupper; + + nlo = static_cast ((sublo[2]-dist[2]-boxlo[2]) * + nz_msm[n]/zprd + OFFSET) - OFFSET; + nhi = static_cast ((subhi[2]+dist[2]-boxlo[2]) * + nz_msm[n]/zprd + OFFSET) - OFFSET; + nzlo_out[n] = nlo + nlower; + nzhi_out[n] = nhi + nupper; + + // 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 is from self since I have ghosts regardless + + int nplanes; + MPI_Status status; + + nplanes = nxlo_in[n] - nxlo_out[n]; + if (comm->procneigh[0][0] != me) + MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[0][0],0, + &nxhi_ghost[n],1,MPI_INT,comm->procneigh[0][1],0, + world,&status); + else nxhi_ghost[n] = nplanes; + + nplanes = nxhi_out[n] - nxhi_in[n]; + if (comm->procneigh[0][1] != me) + MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[0][1],0, + &nxlo_ghost[n],1,MPI_INT,comm->procneigh[0][0], + 0,world,&status); + else nxlo_ghost[n] = nplanes; + + nplanes = nylo_in[n] - nylo_out[n]; + if (comm->procneigh[1][0] != me) + MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[1][0],0, + &nyhi_ghost[n],1,MPI_INT,comm->procneigh[1][1],0, + world,&status); + else nyhi_ghost[n] = nplanes; + + nplanes = nyhi_out[n] - nyhi_in[n]; + if (comm->procneigh[1][1] != me) + MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[1][1],0, + &nylo_ghost[n],1,MPI_INT,comm->procneigh[1][0],0, + world,&status); + else nylo_ghost[n] = nplanes; + + nplanes = nzlo_in[n] - nzlo_out[n]; + if (comm->procneigh[2][0] != me) + MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[2][0],0, + &nzhi_ghost[n],1,MPI_INT,comm->procneigh[2][1],0, + world,&status); + else nzhi_ghost[n] = nplanes; + + nplanes = nzhi_out[n] - nzhi_in[n]; + if (comm->procneigh[2][1] != me) + MPI_Sendrecv(&nplanes,1,MPI_INT,comm->procneigh[2][1],0, + &nzlo_ghost[n],1,MPI_INT,comm->procneigh[2][0],0, + world,&status); + else nzlo_ghost[n] = nplanes; + + int flag = 0; + + if (n == 0) { + + // test that ghost overlap is not bigger than my sub-domain + + if (nxlo_ghost[n] > nxhi_in[n]-nxlo_in[n]+1) flag = 1; + if (nxhi_ghost[n] > nxhi_in[n]-nxlo_in[n]+1) flag = 1; + if (nylo_ghost[n] > nyhi_in[n]-nylo_in[n]+1) flag = 1; + if (nyhi_ghost[n] > nyhi_in[n]-nylo_in[n]+1) flag = 1; + if (nzlo_ghost[n] > nzhi_in[n]-nzlo_in[n]+1) flag = 1; + if (nzhi_ghost[n] > nzhi_in[n]-nzlo_in[n]+1) flag = 1; + } + + int flag_all; + MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); + + if (flag_all != 0) { + char str[128]; + sprintf(str,"MSM parallel communication error, try reducing number of procs"); + error->all(FLERR,str); + } + } + + // Largest MSM grid for this proc, including ghosts + + ngrid = (nxhi_out[0]-nxlo_out[0]+1) * (nyhi_out[0]-nylo_out[0]+1) * + (nzhi_out[0]-nzlo_out[0]+1); + + + // buffer space for use in ghost_swap 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[0],nxhi_ghost[0]); + idelx = MAX(idelx,nxhi_out[0]-nxhi_in[0]); + idelx = MAX(idelx,nxlo_in[0]-nxlo_out[0]); + + idely = MAX(nylo_ghost[0],nyhi_ghost[0]); + idely = MAX(idely,nyhi_out[0]-nyhi_in[0]); + idely = MAX(idely,nylo_in[0]-nylo_out[0]); + + idelz = MAX(nzlo_ghost[0],nzhi_ghost[0]); + idelz = MAX(idelz,nzhi_out[0]-nzhi_in[0]); + idelz = MAX(idelz,nzlo_in[0]-nzlo_out[0]); + + nx = nxhi_out[0] - nxlo_out[0] + 1; + ny = nyhi_out[0] - nylo_out[0] + 1; + nz = nzhi_out[0] - nzlo_out[0] + 1; + + nxx = idelx * ny * nz; + nyy = idely * nx * nz; + nzz = idelz * nx * ny; + + nbuf = MAX(nxx,nyy); + nbuf = MAX(nbuf,nzz); + nbuf *= 3; + + double estimated_error = estimate_total_error(); + + // print stats + + int ngrid_max,nbuf_max; + + // All processors have a copy of the complete grid at each level + + nbuf_max = nbuf; + ngrid_max = ngrid; + + if (me == 0) { + if (screen) { + fprintf(screen," brick buffer size/proc = %d %d\n", + ngrid_max,nbuf_max); + fprintf(screen," estimated absolute RMS force accuracy = %g\n", + estimated_error); + fprintf(screen," estimated relative force accuracy = %g\n", + estimated_error/two_charge_force); + } + if (logfile) { + fprintf(logfile," brick buffer size/proc = %d %d\n", + ngrid_max,nbuf_max); + fprintf(logfile," estimated absolute RMS force accuracy = %g\n", + estimated_error); + fprintf(logfile," estimated relative force accuracy = %g\n", + estimated_error/two_charge_force); + } + } + + // allocate K-space dependent memory + + allocate(); +} + +/* ---------------------------------------------------------------------- + estimate 1d grid RMS force error for MSM from Sec 3.1 of Hardy's thesis +------------------------------------------------------------------------- */ + +double MSM::estimate_1d_error(double h, double prd) +{ + double a = cutoff; + int p = order - 1; + double error_1d = pow(h,(p-1))/pow(a,(p+1)); + error_1d *= 24.0*q2/sqrt(atom->natoms*a*prd*prd); + return error_1d; +} + +/* ---------------------------------------------------------------------- + estimate 3d grid RMS force error +------------------------------------------------------------------------- */ + +double MSM::estimate_3d_error() +{ + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double error_x = estimate_1d_error(xprd/nx_msm[0],xprd); + double error_y = estimate_1d_error(yprd/ny_msm[0],yprd); + double error_z = estimate_1d_error(zprd/nz_msm[0],zprd); + double error_3d = + sqrt(error_x*error_x + error_y*error_y + error_z*error_z) / sqrt(3.0); + return error_3d; +} + +/* ---------------------------------------------------------------------- + estimate total RMS force error +------------------------------------------------------------------------- */ + +double MSM::estimate_total_error() +{ + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + bigint natoms = atom->natoms; + + double grid_error = estimate_3d_error(); + double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd); + double short_range_error = 0.0; + double table_error = + estimate_table_accuracy(q2_over_sqrt,short_range_error); + double estimated_total_error = sqrt(grid_error*grid_error + + short_range_error*short_range_error + table_error*table_error); + + return estimated_total_error; +} + +/* ---------------------------------------------------------------------- + adjust MSM coeffs, called initially and whenever volume has changed +------------------------------------------------------------------------- */ + +void MSM::setup() +{ + int i,j,k,l,m,n; + double *prd; + + double a = cutoff; + + // volume-dependent factors + + prd = domain->prd; + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; + volume = xprd * yprd * zprd; + + // loop over grid levels + + for (int n=0; n (2.0*a*delxinv[0]); + nxlo_direct = -nxhi_direct; + nyhi_direct = static_cast (2.0*a*delyinv[0]); + nylo_direct = -nyhi_direct; + nzhi_direct = static_cast (2.0*a*delzinv[0]); + nzlo_direct = -nzhi_direct; + + nmax_direct = 8*(nxhi_direct+1)*(nyhi_direct+1)*(nzhi_direct+1); + + get_g_direct(); + + boxlo = domain->boxlo; +} + +/* ---------------------------------------------------------------------- + compute the MSM long-range force, energy, virial +------------------------------------------------------------------------- */ + +void MSM::compute(int eflag, int vflag) +{ + int i; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = evflag_atom = eflag_global = vflag_global = + eflag_atom = vflag_atom = 0; + + // extend size of per-atom arrays if necessary + + if (atom->nlocal > nmax) { + memory->destroy(part2grid); + nmax = atom->nmax; + memory->create(part2grid,nmax,3,"msm:part2grid"); + } + + energy = 0.0; + + // find grid points for all my particles + // map my particle charge onto my local 3d density grid (aninterpolation) + + particle_map(); + make_rho(); + + + // all procs communicate density values from their ghost cells + // to fully sum contribution in their 3d bricks + + ghost_swap(0); + + charge_swap(0); + + // Direct sum on finest grid level is parallel + + direct(eflag_global,vflag_global,0); + + potential_swap(0); + + restrict(eflag_global,vflag_global,0); + + // compute potential gradient on my MSM grid and + // portion of e_long on this proc's MSM grid + // return gradients (electric fields) in 3d brick decomposition + + for (int n=1; n=0; n--) + prolongate(eflag_global,vflag_global,n); + + // all procs communicate E-field values + // to fill ghost cells surrounding their 3d bricks + + fillbrick(0); + + // calculate the force on my particles (interpolation) + + fieldforce(); + + // sum energy across procs and add in volume-dependent term + + if (eflag_global) { + double e_self = qsqsum*gamma(0.0)/cutoff; // Self-energy term + energy -= e_self; + double energy_all; + energy *= 0.5*qqrd2e*scale; + } + +} + +/* ---------------------------------------------------------------------- + allocate memory that depends on # of K-vectors and order +------------------------------------------------------------------------- */ + +void MSM::allocate() +{ + memory->create(buf1,nbuf,"msm:buf1"); + memory->create(buf2,nbuf,"msm:buf2"); + + // summation coeffs + + memory->create2d_offset(phi1d,3,nlower-2,nupper+2,"msm:phi1d"); + memory->create2d_offset(dphi1d,3,nlower-2,nupper+2,"msm:dphi1d"); + + // allocate grid levels + + for (int n=0; ncreate3d_offset(qgrid[n],nzlo_out[n],nzhi_out[n], + nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:qgrid"); + memory->create3d_offset(egrid[n],nzlo_out[n],nzhi_out[n], + nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:egrid"); + } + +} + +/* ---------------------------------------------------------------------- + deallocate memory that depends on # of K-vectors and order +------------------------------------------------------------------------- */ + +void MSM::deallocate() +{ + memory->destroy(buf1); + memory->destroy(buf2); + + memory->destroy2d_offset(phi1d,nlower-2); + memory->destroy2d_offset(dphi1d,nlower-2); + + // deallocate grid levels + + for (int n=0; ndestroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); + if (egrid[n]) + memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]); + } +} + +/* ---------------------------------------------------------------------- + allocate memory that depends on # of grid levels +------------------------------------------------------------------------- */ + +void MSM::allocate_levels() +{ + nx_msm = new int[levels]; + ny_msm = new int[levels]; + nz_msm = new int[levels]; + + nxlo_in = new int[levels]; + nylo_in = new int[levels]; + nzlo_in = new int[levels]; + + nxhi_in = new int[levels]; + nyhi_in = new int[levels]; + nzhi_in = new int[levels]; + + nxlo_in_d = new int[levels]; + nylo_in_d = new int[levels]; + nzlo_in_d = new int[levels]; + + nxhi_in_d = new int[levels]; + nyhi_in_d = new int[levels]; + nzhi_in_d = new int[levels]; + + nxlo_out = new int[levels]; + nylo_out = new int[levels]; + nzlo_out = new int[levels]; + + nxhi_out = new int[levels]; + nyhi_out = new int[levels]; + nzhi_out = new int[levels]; + + nxlo_ghost = new int[levels]; + nylo_ghost = new int[levels]; + nzlo_ghost = new int[levels]; + + nxhi_ghost = new int[levels]; + nyhi_ghost = new int[levels]; + nzhi_ghost = new int[levels]; + + delxinv = new double[levels]; + delyinv = new double[levels]; + delzinv = new double[levels]; + delvolinv = new double[levels]; + + qgrid = new double***[levels]; + egrid = new double***[levels]; + +} + +/* ---------------------------------------------------------------------- + deallocate memory that depends on # of grid levels +------------------------------------------------------------------------- */ + +void MSM::deallocate_levels() +{ + + delete [] nx_msm; + delete [] ny_msm; + delete [] nz_msm; + + delete [] nxlo_in; + delete [] nylo_in; + delete [] nzlo_in; + + delete [] nxhi_in; + delete [] nyhi_in; + delete [] nzhi_in; + + delete [] nxlo_in_d; + delete [] nylo_in_d; + delete [] nzlo_in_d; + + delete [] nxhi_in_d; + delete [] nyhi_in_d; + delete [] nzhi_in_d; + + delete [] nxlo_out; + delete [] nylo_out; + delete [] nzlo_out; + + delete [] nxhi_out; + delete [] nyhi_out; + delete [] nzhi_out; + + delete [] nxlo_ghost; + delete [] nylo_ghost; + delete [] nzlo_ghost; + + delete [] nxhi_ghost; + delete [] nyhi_ghost; + delete [] nzhi_ghost; + + delete [] delxinv; + delete [] delyinv; + delete [] delzinv; + delete [] delvolinv; + + delete [] qgrid; + delete [] egrid; + +} + +/* ---------------------------------------------------------------------- + set size of MSM grids +------------------------------------------------------------------------- */ + +void MSM::set_grid() +{ + if (accuracy_relative <= 0.0) + error->all(FLERR,"KSpace accuracy must be > 0"); + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + int nx_max,ny_max,nz_max; + + if (!gridflag) { + nx_max = ny_max = nz_max = 2; + double hx = xprd/nx_max; + double hy = yprd/ny_max; + double hz = zprd/nz_max; + + double x_error = 2.0*accuracy; + double y_error = 2.0*accuracy; + double z_error = 2.0*accuracy; + + while (x_error > accuracy) { + nx_max *= 2; + hx = xprd/nx_max; + x_error = estimate_1d_error(hx,xprd); + } + + while (y_error > accuracy) { + ny_max *= 2; + hy = yprd/ny_max; + y_error = estimate_1d_error(hy,yprd); + } + + while (z_error > accuracy) { + nz_max *= 2; + hz = zprd/nz_max; + z_error = estimate_1d_error(hz,zprd); + } + } else { + nx_max = nx_msm_max; + ny_max = ny_msm_max; + nz_max = nz_msm_max; + } + + // boost grid size until it is factorable + + int flag = 0; + int xlevels,ylevels,zlevels; + + while (!factorable(nx_max,flag,xlevels)) nx_max++; + while (!factorable(ny_max,flag,ylevels)) ny_max++; + while (!factorable(nz_max,flag,zlevels)) nz_max++; + + if (flag) + error->warning(FLERR,"Number of MSM mesh points increased to be a multiple of 2"); + + // Find maximum number of levels + + levels = MAX(xlevels,ylevels); + levels = MAX(levels,zlevels); + + if (levels > MAX_LEVELS) + error->all(FLERR,"Too many MSM grid levels"); + + allocate_levels(); + + for (int n = 0; n < levels; n++) { + + if (xlevels-n-1 > 0) + nx_msm[n] = static_cast (pow(2.0,xlevels-n-1)); + else + nx_msm[n] = 1; + + if (ylevels-n-1 > 0) + ny_msm[n] = static_cast (pow(2.0,ylevels-n-1)); + else + ny_msm[n] = 1; + + if (zlevels-n-1 > 0) + nz_msm[n] = static_cast (pow(2.0,zlevels-n-1)); + else + nz_msm[n] = 1; + } + + // Output grid stats + + if (me == 0) { + if (screen) { + fprintf(screen," grid = %d %d %d\n",nx_msm[0],ny_msm[0],nz_msm[0]); + fprintf(screen," stencil order = %d\n",order); + } + if (logfile) { + fprintf(logfile," grid = %d %d %d\n",nx_msm[0],ny_msm[0],nz_msm[0]); + fprintf(logfile," stencil order = %d\n",order); + } + } + +} + +/* ---------------------------------------------------------------------- + check if all factors of n are in list of factors + return 1 if yes, 0 if no +------------------------------------------------------------------------- */ + +int MSM::factorable(int n, int &flag, int &levels) +{ + int i,norig; + norig = n; + levels = 1; + + while (n > 1) { + for (i = 0; i < nfactors; i++) { + if (n % factors[i] == 0) { + n /= factors[i]; + levels++; + break; + } + } + if (i == nfactors) { + flag = 1; + return 0; + } + } + + return 1; +} + +/* ---------------------------------------------------------------------- + MPI-Reduce so each processor has all the info it needs +------------------------------------------------------------------------- */ +void MSM::charge_swap(int n) +{ + double ***qgridn = qgrid[n]; + double ***qgridn_all; + memory->create3d_offset(qgridn_all,nzlo_out[n],nzhi_out[n],nylo_out[n],nyhi_out[n], + nxlo_out[n],nxhi_out[n],"msm:qgrid_all"); + + memset(&(qgridn_all[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid*sizeof(double)); + + MPI_Allreduce(&(qgridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]), + &(qgridn_all[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]), + ngrid,MPI_DOUBLE,MPI_SUM,world); + + // Swap pointers between qgridn and qgridn_all to avoid need of copy operation + + double ***tmp; + tmp = qgridn; + qgrid[n] = qgridn_all; + qgridn_all = tmp; + + memory->destroy3d_offset(qgridn_all,nzlo_out[n],nylo_out[n],nxlo_out[n]); + +} + +/* ---------------------------------------------------------------------- + MPI-Reduce so each processor has all the info it needs +------------------------------------------------------------------------- */ +void MSM::potential_swap(int n) +{ + double ***egridn = egrid[n]; + double ***egridn_all; + memory->create3d_offset(egridn_all,nzlo_out[n],nzhi_out[n],nylo_out[n],nyhi_out[n], + nxlo_out[n],nxhi_out[n],"msm:qgrid_all"); + + memset(&(egridn_all[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid*sizeof(double)); + + MPI_Allreduce(&(egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]), + &(egridn_all[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]), + ngrid,MPI_DOUBLE,MPI_SUM,world); + + // Swap pointers between egridn and egridn_all to avoid need of copy operation + + double ***tmp; + tmp = egridn; + egrid[n] = egridn_all; + egridn_all = tmp; + + memory->destroy3d_offset(egridn_all,nzlo_out[n],nylo_out[n],nxlo_out[n]); + +} + +/* ---------------------------------------------------------------------- + ghost-swap to accumulate full density in brick decomposition +------------------------------------------------------------------------- */ + +void MSM::ghost_swap(int n) +{ + double ***qgridn = qgrid[n]; + + int i,k,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 + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy <= nyhi_out[n]; iy++) + for (ix = nxhi_in[n]+1; ix <= nxhi_out[n]; ix++) + buf1[k++] = qgridn[iz][iy][ix]; + + if (comm->procneigh[0][1] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[0][0],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[0][1],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy <= nyhi_out[n]; iy++) + for (ix = nxlo_in[n]; ix < nxlo_in[n]+nxlo_ghost[n]; ix++) + qgridn[iz][iy][ix] += buf2[k++]; + + // pack my ghosts for -x processor + // pass data to self or -x processor + // unpack and sum recv data into my real cells + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy <= nyhi_out[n]; iy++) + for (ix = nxlo_out[n]; ix < nxlo_in[n]; ix++) + buf1[k++] = qgridn[iz][iy][ix]; + + if (comm->procneigh[0][0] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[0][1],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[0][0],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy <= nyhi_out[n]; iy++) + for (ix = nxhi_in[n]-nxhi_ghost[n]+1; ix <= nxhi_in[n]; ix++) + qgridn[iz][iy][ix] += buf2[k++]; + + // pack my ghosts for +y processor + // pass data to self or +y processor + // unpack and sum recv data into my real cells + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nyhi_in[n]+1; iy <= nyhi_out[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) + buf1[k++] = qgridn[iz][iy][ix]; + + if (comm->procneigh[1][1] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[1][0],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[1][1],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_in[n]; iy < nylo_in[n]+nylo_ghost[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) + qgridn[iz][iy][ix] += buf2[k++]; + + // pack my ghosts for -y processor + // pass data to self or -y processor + // unpack and sum recv data into my real cells + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy < nylo_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) + buf1[k++] = qgridn[iz][iy][ix]; + + if (comm->procneigh[1][0] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[1][1],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[1][0],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nyhi_in[n]-nyhi_ghost[n]+1; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) + qgridn[iz][iy][ix] += buf2[k++]; + + // pack my ghosts for +z processor + // pass data to self or +z processor + // unpack and sum recv data into my real cells + + k = 0; + for (iz = nzhi_in[n]+1; iz <= nzhi_out[n]; iz++) + for (iy = nylo_in[n]; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) + buf1[k++] = qgridn[iz][iy][ix]; + + if (comm->procneigh[2][1] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[2][0],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[2][1],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_in[n]; iz < nzlo_in[n]+nzlo_ghost[n]; iz++) + for (iy = nylo_in[n]; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) + qgridn[iz][iy][ix] += buf2[k++]; + + // pack my ghosts for -z processor + // pass data to self or -z processor + // unpack and sum recv data into my real cells + + k = 0; + for (iz = nzlo_out[n]; iz < nzlo_in[n]; iz++) + for (iy = nylo_in[n]; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) + buf1[k++] = qgridn[iz][iy][ix]; + + if (comm->procneigh[2][0] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[2][1],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[2][0],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzhi_in[n]-nzhi_ghost[n]+1; iz <= nzhi_in[n]; iz++) + for (iy = nylo_in[n]; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) + qgridn[iz][iy][ix] += buf2[k++]; + +} + +/* ---------------------------------------------------------------------- + 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 MSM::particle_map() +{ + int nx,ny,nz; + + double **x = atom->x; + int nlocal = atom->nlocal; + + 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]-boxlo[0])*delxinv[0]+OFFSET) - OFFSET; + ny = static_cast ((x[i][1]-boxlo[1])*delyinv[0]+OFFSET) - OFFSET; + nz = static_cast ((x[i][2]-boxlo[2])*delzinv[0]+OFFSET) - 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[0] || nx+nupper > nxhi_out[0] || + ny+nlower < nylo_out[0] || ny+nupper > nyhi_out[0] || + nz+nlower < nzlo_out[0] || nz+nupper > nzhi_out[0]) flag = 1; + } + + if (flag) error->one(FLERR,"Out of range atoms - cannot compute MSM"); +} + +/* ---------------------------------------------------------------------- + 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 MSM::make_rho() +{ + //fprintf(screen,"MSM aninterpolation\n\n"); + + int i,l,m,n,nn,nx,ny,nz,mx,my,mz; + double dx,dy,dz,x0,y0,z0; + + // clear 3d density array + + double ***qgridn = qgrid[0]; + + memset(&(qgridn[nzlo_out[0]][nylo_out[0]][nxlo_out[0]]),0,ngrid*sizeof(double)); + + // 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; + + for (i = 0; i < nlocal; i++) { + + nx = part2grid[i][0]; + ny = part2grid[i][1]; + nz = part2grid[i][2]; + dx = nx - (x[i][0]-boxlo[0])*delxinv[0]; + dy = ny - (x[i][1]-boxlo[1])*delyinv[0]; + dz = nz - (x[i][2]-boxlo[2])*delzinv[0]; + + compute_phis_and_dphis(dx,dy,dz); + + z0 = q[i]; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + y0 = z0*phi1d[2][n]; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + x0 = y0*phi1d[1][m]; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + qgridn[mz][my][mx] += x0*phi1d[0][l]; + } + } + } + } + +} + +/* ---------------------------------------------------------------------- + MSM direct part procedure for intermediate grid levels +------------------------------------------------------------------------- */ + +void MSM::direct(int eflag, int vflag, int n) +{ + //fprintf(screen,"Direct contribution on level %i\n\n",n); + + double ***egridn = egrid[n]; + double ***qgridn = qgrid[n]; + + // bitmask for PBCs (only works for power of 2 numbers) + + int PBCx,PBCy,PBCz; + + PBCx = nx_msm[n]-1; + PBCy = ny_msm[n]-1; + PBCz = nz_msm[n]-1; + + // zero out electric field brick + + for (int icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) + for (int icy = nylo_in[n]; icy <= nyhi_in[n]; icy++) + for (int icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++) + egridn[icz][icy][icx] = 0.0; + + // Simple parallelization of direct sum + + for (int icz = nzlo_in_d[n]; icz <= nzhi_in_d[n]; icz++) { + for (int icy = nylo_in_d[n]; icy <= nyhi_in_d[n]; icy++) { + for (int icx = nxlo_in_d[n]; icx <= nxhi_in_d[n]; icx++) { + + + // do double loop over points on the intermediate grid level + // for now, assume I own all points on the intermediate grid level + + int k = 0; + for (int iz = nzlo_direct; iz <= nzhi_direct; iz++) { + for (int iy = nylo_direct; iy <= nyhi_direct; iy++) { + for (int ix = nxlo_direct; ix <= nxhi_direct; ix++) { + egridn[icz][icy][icx] += g_direct[n][k++] + * qgridn[(icz+iz) & PBCz][(icy+iy) & PBCy][(icx+ix) & PBCx]; + } + } + } + + } + } + } + +} + +/* ---------------------------------------------------------------------- + MSM restrict procedure for intermediate grid levels +------------------------------------------------------------------------- */ + +void MSM::restrict(int eflag, int vflag,int n) +{ + //fprintf(screen,"Restricting from level %i to %i\n\n",n,n+1); + + int p = order-1; + + double ***qgrid1 = qgrid[n]; + double ***qgrid2 = qgrid[n+1]; + + // bitmask for PBCs (only works for power of 2 numbers) + + int PBCx,PBCy,PBCz; + + PBCx = nx_msm[n]-1; + PBCy = ny_msm[n]-1; + PBCz = nz_msm[n]-1; + + //restrict grid (going from grid n to grid n+1, i.e. to a coarser grid) + + for (int nu=-p; nu<=p; nu++) { + phi1d[0][nu] = compute_phi(nu*delxinv[n+1]/delxinv[n]); + phi1d[1][nu] = compute_phi(nu*delyinv[n+1]/delyinv[n]); + phi1d[2][nu] = compute_phi(nu*delzinv[n+1]/delzinv[n]); + } + + for (int kp = nzlo_in[n+1]; kp <= nzhi_in[n+1]; kp++) + for (int jp = nylo_in[n+1]; jp <= nyhi_in[n+1]; jp++) + for (int ip = nxlo_in[n+1]; ip <= nxhi_in[n+1]; ip++) { + qgrid2[kp][jp][ip] = 0.0; + + int ic = static_cast (ip*delxinv[n]/delxinv[n+1]); + int jc = static_cast (jp*delyinv[n]/delyinv[n+1]); + int kc = static_cast (kp*delzinv[n]/delzinv[n+1]); + for (int k=-p; k<=p; k++) // Could make this faster by eliminating zeros + for (int j=-p; j<=p; j++) + for (int i=-p; i<=p; i++) + qgrid2[kp][jp][ip] += + qgrid1[(kc+k)&PBCz][(jc+j)&PBCy][(ic+i)&PBCx] * + phi1d[0][i]*phi1d[1][j]*phi1d[2][k]; + } + + +} + +/* ---------------------------------------------------------------------- + MSM prolongate procedure for intermediate grid levels +------------------------------------------------------------------------- */ + +void MSM::prolongate(int eflag, int vflag,int n) +{ + //fprintf(screen,"Prolongating from level %i to %i\n\n",n+1,n); + + int p = order-1; + + double ***egrid1 = egrid[n]; + double ***egrid2 = egrid[n+1]; + + // bitmask for PBCs (only works for power of 2 numbers) + + int PBCx,PBCy,PBCz; + + PBCx = nx_msm[n]-1; + PBCy = ny_msm[n]-1; + PBCz = nz_msm[n]-1; + + //prolongate grid (going from grid n to grid n-1, i.e. to a finer grid) + + for (int nu=-p; nu<=p; nu++) { + phi1d[0][nu] = compute_phi(nu*delxinv[n+1]/delxinv[n]); + phi1d[1][nu] = compute_phi(nu*delyinv[n+1]/delyinv[n]); + phi1d[2][nu] = compute_phi(nu*delzinv[n+1]/delzinv[n]); + } + + for (int kp = nzlo_in[n+1]; kp <= nzhi_in[n+1]; kp++) + for (int jp = nylo_in[n+1]; jp <= nyhi_in[n+1]; jp++) + for (int ip = nxlo_in[n+1]; ip <= nxhi_in[n+1]; ip++) { + + int ic = static_cast (ip*delxinv[n]/delxinv[n+1]); + int jc = static_cast (jp*delyinv[n]/delyinv[n+1]); + int kc = static_cast (kp*delzinv[n]/delzinv[n+1]); + for (int k=-p; k<=p; k++) // Could make this faster by eliminating zeros + for (int j=-p; j<=p; j++) + for (int i=-p; i<=p; i++) + egrid1[(kc+k)&PBCz][(jc+j)&PBCy][(ic+i)&PBCx] += + egrid2[kp][jp][ip] * phi1d[0][i]*phi1d[1][j]*phi1d[2][k]; + } + +} + +/* ---------------------------------------------------------------------- + ghost-swap to fill ghost cells of my brick with field values +------------------------------------------------------------------------- */ + +void MSM::fillbrick(int n) +{ + double ***egridn = egrid[n]; + + int i,k,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 + + k = 0; + for (iz = nzhi_in[n]-nzhi_ghost[n]+1; iz <= nzhi_in[n]; iz++) + for (iy = nylo_in[n]; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) { + buf1[k++] = egridn[iz][iy][ix]; + } + + if (comm->procneigh[2][1] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[2][0],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[2][1],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz < nzlo_in[n]; iz++) + for (iy = nylo_in[n]; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) { + egridn[iz][iy][ix] = buf2[k++]; + } + + // pack my real cells for -z processor + // pass data to self or -z processor + // unpack and sum recv data into my ghost cells + + k = 0; + for (iz = nzlo_in[n]; iz < nzlo_in[n]+nzlo_ghost[n]; iz++) + for (iy = nylo_in[n]; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) { + buf1[k++] = egridn[iz][iy][ix]; + } + + if (comm->procneigh[2][0] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[2][1],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[2][0],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzhi_in[n]+1; iz <= nzhi_out[n]; iz++) + for (iy = nylo_in[n]; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) { + egridn[iz][iy][ix] = buf2[k++]; + } + + // pack my real cells for +y processor + // pass data to self or +y processor + // unpack and sum recv data into my ghost cells + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nyhi_in[n]-nyhi_ghost[n]+1; iy <= nyhi_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) { + buf1[k++] = egridn[iz][iy][ix]; + } + + if (comm->procneigh[1][1] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[1][0],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[1][1],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy < nylo_in[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) { + egridn[iz][iy][ix] = buf2[k++]; + } + + // pack my real cells for -y processor + // pass data to self or -y processor + // unpack and sum recv data into my ghost cells + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_in[n]; iy < nylo_in[n]+nylo_ghost[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) { + buf1[k++] = egridn[iz][iy][ix]; + } + + if (comm->procneigh[1][0] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[1][1],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[1][0],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nyhi_in[n]+1; iy <= nyhi_out[n]; iy++) + for (ix = nxlo_in[n]; ix <= nxhi_in[n]; ix++) { + egridn[iz][iy][ix] = buf2[k++]; + } + + // pack my real cells for +x processor + // pass data to self or +x processor + // unpack and sum recv data into my ghost cells + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy <= nyhi_out[n]; iy++) + for (ix = nxhi_in[n]-nxhi_ghost[n]+1; ix <= nxhi_in[n]; ix++) { + buf1[k++] = egridn[iz][iy][ix]; + } + + if (comm->procneigh[0][1] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[0][0],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[0][1],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy <= nyhi_out[n]; iy++) + for (ix = nzlo_out[n]; ix < nxlo_in[n]; ix++) { + egridn[iz][iy][ix] = buf2[k++]; + } + + // pack my real cells for -x processor + // pass data to self or -x processor + // unpack and sum recv data into my ghost cells + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy <= nyhi_out[n]; iy++) + for (ix = nxlo_in[n]; ix < nxlo_in[n]+nxlo_ghost[n]; ix++) { + buf1[k++] = egridn[iz][iy][ix]; + } + + if (comm->procneigh[0][0] == me) + for (i = 0; i < k; i++) buf2[i] = buf1[i]; + else { + MPI_Irecv(buf2,nbuf,MPI_DOUBLE,comm->procneigh[0][1],0,world,&request); + MPI_Send(buf1,k,MPI_DOUBLE,comm->procneigh[0][0],0,world); + MPI_Wait(&request,&status); + } + + k = 0; + for (iz = nzlo_out[n]; iz <= nzhi_out[n]; iz++) + for (iy = nylo_out[n]; iy <= nyhi_out[n]; iy++) + for (ix = nxhi_in[n]+1; ix <= nxhi_out[n]; ix++) { + egridn[iz][iy][ix] = buf2[k++]; + } +} + +/* ---------------------------------------------------------------------- + interpolate from grid to get force on my particles +------------------------------------------------------------------------- */ + +void MSM::fieldforce() +{ + //fprintf(screen,"MSM interpolation\n\n"); + + double ***egridn = egrid[0]; + double ***qgridn = qgrid[0]; + + int i,l,m,n,nx,ny,nz,mx,my,mz; + double dx,dy,dz; + double phi_x,phi_y,phi_z; + double dphi_x,dphi_y,dphi_z; + double ekx,eky,ekz; + + + // 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; + + for (i = 0; i < nlocal; i++) { + nx = part2grid[i][0]; + ny = part2grid[i][1]; + nz = part2grid[i][2]; + dx = nx - (x[i][0]-boxlo[0])*delxinv[0]; + dy = ny - (x[i][1]-boxlo[1])*delyinv[0]; + dz = nz - (x[i][2]-boxlo[2])*delzinv[0]; + + compute_phis_and_dphis(dx,dy,dz); + + ekx = eky = ekz = 0.0; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + phi_z = phi1d[2][n]; + dphi_z = dphi1d[2][n]; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + phi_y = phi1d[1][m]; + dphi_y = dphi1d[1][m]; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + phi_x = phi1d[0][l]; + dphi_x = dphi1d[0][l]; + ekx += dphi_x*phi_y*phi_z*egridn[mz][my][mx]; + eky += phi_x*dphi_y*phi_z*egridn[mz][my][mx]; + ekz += phi_x*phi_y*dphi_z*egridn[mz][my][mx]; + } + } + } + + ekx *= delxinv[0]; + eky *= delyinv[0]; + ekz *= delzinv[0]; + + // convert E-field to force + const double qfactor = force->qqrd2e*scale*q[i]; + f[i][0] += qfactor*ekx; + f[i][1] += qfactor*eky; + f[i][2] += qfactor*ekz; + + } + + // Sum total long-range energy + + for (int kp=0; kpdestroy(g_direct); + memory->create(g_direct,levels,nmax_direct,"msm:g_direct"); + + double a = cutoff; + + for (int n=0; n grid mapping + int nmax; + + int triclinic; // domain settings, orthog or triclinic + double *boxlo; + + void set_grid(); + double estimate_1d_error(double,double); + double estimate_3d_error(); + double estimate_total_error(); + void allocate(); + void deallocate(); + void allocate_levels(); + void deallocate_levels(); + int factorable(int,int&,int&); + void compute_gf_denom(); + double gf_denom(double, double, double); + void particle_map(); + void make_rho(); + void ghost_swap(int); + void charge_swap(int); + void potential_swap(int); + void direct(int, int, int); + void restrict(int,int,int); + void prolongate(int,int,int); + void fillbrick(int); + void fieldforce(); + void procs2grid2d(int,int,int,int*,int*); + void compute_phis_and_dphis(const double &, const double &, const double &); + double compute_phi(const double &); + double compute_dphi(const double &); + void get_g_direct(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Cannot (yet) use MSM with triclinic box + +This feature is not yet supported. + +E: Cannot (yet) use MSM with 2d simulation + +This feature is not yet supported. + +E: Kspace style requires atom attribute q + +The atom style defined does not have these attributes. + +E: Cannot (yet) use nonperiodic boundaries with MSM + +This feature is not yet supported. + +E: Cannot use slab correction with MSM + +Slab correction can only be used with Ewald and PPPM, not MSM + +E: MSM order cannot be < 2 or > than %d + +This is a limitation of the MSM implementation in LAMMPS. +Currently the only available order is 4. + +E: KSpace style is incompatible with Pair style + +Setting a kspace style requires that a pair style with a long-range +Coulombic component be selected that is compatible with MSM. Note +that TIP4P is not (yet) supported by MSM + +E: Cannot use kspace solver on system with no charge + +No atoms in system have a non-zero charge. + +E: System is not charge neutral, net charge = %g + +The total charge on all atoms on the system is not 0.0, which +is not valid for MSM. + +W: MSM parallel communication error, try reducing accuracy or number of procs + +Currently only nearest neighbor communication between processors is implemented in MSM. +If charge from an atom spans more than one processor domain this error will result. + +E: MSM grid is too large + +The global MSM grid is larger than OFFSET in one or more dimensions. +OFFSET is currently set to 16384. You likely need to decrease the +requested accuracy. + +E: KSpace accuracy must be > 0 + +The kspace accuracy designated in the input must be greater than zero. + +E: Out of range atoms - cannot compute MSM + +One or more atoms are attempting to map their charge to a MSM grid +point that is not owned by a processor. This is likely for one of two +reasons, both of them bad. First, it may mean that an atom near the +boundary of a processor's sub-domain has moved more than 1/2 the +"neighbor skin distance"_neighbor.html without neighbor lists being +rebuilt and atoms being migrated to new processors. This also means +you may be missing pairwise interactions that need to be computed. +The solution is to change the re-neighboring criteria via the +"neigh_modify"_neigh_modify command. The safest settings are "delay 0 +every 1 check yes". Second, it may mean that an atom has moved far +outside a processor's sub-domain or even the entire simulation box. +This indicates bad physics, e.g. due to highly overlapping atoms, too +large a timestep, etc. + +*/ diff --git a/src/KSPACE/pair_born_coul_long.h b/src/KSPACE/pair_born_coul_long.h index 5f06a43c97..dc60fb13f6 100644 --- a/src/KSPACE/pair_born_coul_long.h +++ b/src/KSPACE/pair_born_coul_long.h @@ -37,8 +37,8 @@ class PairBornCoulLong : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); - double single(int, int, int, int, double, double, double, double &); - void *extract(const char *, int &); + virtual double single(int, int, int, int, double, double, double, double &); + virtual void *extract(const char *, int &); protected: double cut_lj_global; diff --git a/src/KSPACE/pair_born_coul_msm.cpp b/src/KSPACE/pair_born_coul_msm.cpp new file mode 100644 index 0000000000..85d960145d --- /dev/null +++ b/src/KSPACE/pair_born_coul_msm.cpp @@ -0,0 +1,194 @@ +/* ---------------------------------------------------------------------- + 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: Stan Moore (SNL), Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_born_coul_msm.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairBornCoulMSM::PairBornCoulMSM(LAMMPS *lmp) : PairBornCoulLong(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void PairBornCoulMSM::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj; + double egamma,fgamma,prefactor; + double r,rexp; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + 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_coulsq) { + r = sqrt(rsq); + prefactor = qqrd2e * qtmp*q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]); + forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv + + born3[itype][jtype]*r2inv*r6inv; + } else forceborn = 0.0; + + fpair = (forcecoul + factor_lj*forceborn) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + ecoul = prefactor*egamma; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv + + d[itype][jtype]*r6inv*r2inv - offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +double PairBornCoulMSM::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,rexp,egamma,fgamma,prefactor; + double forcecoul,forceborn,phicoul,phiborn; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + r = sqrt(rsq); + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + r = sqrt(rsq); + rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]); + forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv + + born3[itype][jtype]*r2inv*r6inv; + } else forceborn = 0.0; + fforce = (forcecoul + factor_lj*forceborn) * r2inv; + + double eng = 0.0; + if (rsq < cut_coulsq) { + phicoul = prefactor*egamma; + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + if (rsq < cut_ljsq[itype][jtype]) { + phiborn = a[itype][jtype]*rexp - c[itype][jtype]*r6inv + + d[itype][jtype]*r2inv*r6inv - offset[itype][jtype]; + eng += factor_lj*phiborn; + } + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairBornCoulMSM::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_msm") == 0) return (void *) &cut_coul; + return NULL; +} diff --git a/src/KSPACE/pair_born_coul_msm.h b/src/KSPACE/pair_born_coul_msm.h new file mode 100644 index 0000000000..b4d8368d1f --- /dev/null +++ b/src/KSPACE/pair_born_coul_msm.h @@ -0,0 +1,68 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_CLASS + +PairStyle(born/coul/msm,PairBornCoulMSM) + +#else + +#ifndef LMP_PAIR_BORN_COUL_MSM_H +#define LMP_PAIR_BORN_COUL_MSM_H + +#include "pair_born_coul_long.h" + +namespace LAMMPS_NS { + +class PairBornCoulMSM : public PairBornCoulLong { + public: + PairBornCoulMSM(class LAMMPS *); + virtual ~PairBornCoulMSM(){}; + virtual void compute(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + virtual void *extract(const char *, int &); + +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: All pair coeffs are not set + +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. + +E: Pair style born/coul/long requires atom attribute q + +An atom style that defines this attribute must be used. + +E: Pair style is incompatible with KSpace style + +If a pair style with a long-range Coulombic component is selected, +then a kspace style must also be used. + +*/ diff --git a/src/KSPACE/pair_buck_coul_long.h b/src/KSPACE/pair_buck_coul_long.h index 17cd53056e..05903b6870 100644 --- a/src/KSPACE/pair_buck_coul_long.h +++ b/src/KSPACE/pair_buck_coul_long.h @@ -28,7 +28,7 @@ class PairBuckCoulLong : public Pair { public: PairBuckCoulLong(class LAMMPS *); virtual ~PairBuckCoulLong(); - void compute(int, int); + virtual void compute(int, int); void settings(int, char **); void coeff(int, char **); void init_style(); @@ -37,8 +37,8 @@ class PairBuckCoulLong : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); - double single(int, int, int, int, double, double, double, double &); - void *extract(const char *, int &); + virtual double single(int, int, int, int, double, double, double, double &); + virtual void *extract(const char *, int &); protected: double cut_lj_global; diff --git a/src/KSPACE/pair_buck_coul_msm.cpp b/src/KSPACE/pair_buck_coul_msm.cpp new file mode 100644 index 0000000000..ffaf111992 --- /dev/null +++ b/src/KSPACE/pair_buck_coul_msm.cpp @@ -0,0 +1,189 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_buck_coul_msm.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + + +/* ---------------------------------------------------------------------- */ + +PairBuckCoulMSM::PairBuckCoulMSM(LAMMPS *lmp) : PairBuckCoulLong(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void PairBuckCoulMSM::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double rsq,r2inv,r6inv,forcecoul,forcebuck,factor_coul,factor_lj; + double egamma,fgamma,prefactor; + double r,rexp; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + 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_coulsq) { + r = sqrt(rsq); + prefactor = qqrd2e * qtmp*q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + rexp = exp(-r*rhoinv[itype][jtype]); + forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; + } else forcebuck = 0.0; + + fpair = (forcecoul + factor_lj*forcebuck) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + ecoul = prefactor*egamma; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv - + offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +double PairBuckCoulMSM::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,rexp,egamma,fgamma,prefactor; + double forcecoul,forcebuck,phicoul,phibuck; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + r = sqrt(rsq); + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + r = sqrt(rsq); + rexp = exp(-r*rhoinv[itype][jtype]); + forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; + } else forcebuck = 0.0; + fforce = (forcecoul + factor_lj*forcebuck) * r2inv; + + double eng = 0.0; + if (rsq < cut_coulsq) { + phicoul = prefactor*egamma; + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + if (rsq < cut_ljsq[itype][jtype]) { + phibuck = a[itype][jtype]*rexp - c[itype][jtype]*r6inv - + offset[itype][jtype]; + eng += factor_lj*phibuck; + } + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairBuckCoulMSM::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_msm") == 0) return (void *) &cut_coul; + return NULL; +} diff --git a/src/KSPACE/pair_buck_coul_msm.h b/src/KSPACE/pair_buck_coul_msm.h new file mode 100644 index 0000000000..95f63268a8 --- /dev/null +++ b/src/KSPACE/pair_buck_coul_msm.h @@ -0,0 +1,68 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_CLASS + +PairStyle(buck/coul/msm,PairBuckCoulMSM) + +#else + +#ifndef LMP_PAIR_BUCK_COUL_MSM_H +#define LMP_PAIR_BUCK_COUL_MSM_H + +#include "pair_buck_coul_long.h" + +namespace LAMMPS_NS { + +class PairBuckCoulMSM : public PairBuckCoulLong { + public: + PairBuckCoulMSM(class LAMMPS *); + virtual ~PairBuckCoulMSM(){}; + virtual void compute(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + virtual void *extract(const char *, int &); + +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: All pair coeffs are not set + +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. + +E: Pair style buck/coul/long requires atom attribute q + +The atom style defined does not have these attributes. + +E: Pair style is incompatible with KSpace style + +If a pair style with a long-range Coulombic component is selected, +then a kspace style must also be used. + +*/ diff --git a/src/KSPACE/pair_coul_long.h b/src/KSPACE/pair_coul_long.h index ec7398f5a0..70d97c1897 100644 --- a/src/KSPACE/pair_coul_long.h +++ b/src/KSPACE/pair_coul_long.h @@ -38,7 +38,7 @@ class PairCoulLong : public Pair { virtual void write_restart_settings(FILE *); virtual void read_restart_settings(FILE *); virtual double single(int, int, int, int, double, double, double, double &); - void *extract(const char *, int &); + virtual void *extract(const char *, int &); protected: double cut_coul,cut_coulsq; @@ -52,7 +52,7 @@ class PairCoulLong : public Pair { int ncoulshiftbits,ncoulmask; void allocate(); - void init_tables(); + virtual void init_tables(); void free_tables(); }; diff --git a/src/KSPACE/pair_coul_msm.cpp b/src/KSPACE/pair_coul_msm.cpp new file mode 100644 index 0000000000..46742252f1 --- /dev/null +++ b/src/KSPACE/pair_coul_msm.cpp @@ -0,0 +1,379 @@ +/* ---------------------------------------------------------------------- + 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: Stan Moore (SNL), Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_coul_msm.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairCoulMSM::PairCoulMSM(LAMMPS *lmp) : PairCoulLong(lmp) +{ + +} + +/* ---------------------------------------------------------------------- */ + +void PairCoulMSM::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itable,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; + double fraction,table; + double r,r2inv,forcecoul,factor_coul; + double egamma,fgamma,prefactor; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + 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 < cut_coulsq) { + r2inv = 1.0/rsq; + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = scale[itype][jtype] * qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = scale[itype][jtype] * qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + + fpair = forcecoul * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*egamma; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = scale[itype][jtype] * qtmp*q[j] * table; + } + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + 0.0,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + + setup force tables used in compute routines +------------------------------------------------------------------------- */ + +void PairCoulMSM::init_tables() +{ + int masklo,maskhi; + double r,egamma,fgamma,rsw; + double qqrd2e = force->qqrd2e; + + tabinnersq = tabinner*tabinner; + init_bitmap(tabinner,cut_coul,ncoultablebits, + masklo,maskhi,ncoulmask,ncoulshiftbits); + + int ntable = 1; + for (int i = 0; i < ncoultablebits; i++) ntable *= 2; + + // linear lookup tables of length N = 2^ncoultablebits + // stored value = value at lower edge of bin + // d values = delta from lower edge to upper edge of bin + + if (ftable) free_tables(); + + memory->create(rtable,ntable,"pair:rtable"); + memory->create(ftable,ntable,"pair:ftable"); + memory->create(ctable,ntable,"pair:ctable"); + memory->create(etable,ntable,"pair:etable"); + memory->create(drtable,ntable,"pair:drtable"); + memory->create(dftable,ntable,"pair:dftable"); + memory->create(dctable,ntable,"pair:dctable"); + memory->create(detable,ntable,"pair:detable"); + + if (cut_respa == NULL) { + vtable = ptable = dvtable = dptable = NULL; + } else { + memory->create(vtable,ntable,"pair:vtable"); + memory->create(ptable,ntable,"pair:ptable"); + memory->create(dvtable,ntable,"pair:dvtable"); + memory->create(dptable,ntable,"pair:dptable"); + } + + union_int_float_t rsq_lookup; + union_int_float_t minrsq_lookup; + int itablemin; + minrsq_lookup.i = 0 << ncoulshiftbits; + minrsq_lookup.i |= maskhi; + + for (int i = 0; i < ntable; i++) { + rsq_lookup.i = i << ncoulshiftbits; + rsq_lookup.i |= masklo; + if (rsq_lookup.f < tabinnersq) { + rsq_lookup.i = i << ncoulshiftbits; + rsq_lookup.i |= maskhi; + } + r = sqrtf(rsq_lookup.f); + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + if (cut_respa == NULL) { + rtable[i] = rsq_lookup.f; + ftable[i] = qqrd2e/r * fgamma; + ctable[i] = qqrd2e/r; + etable[i] = qqrd2e/r * egamma; + } else { + rtable[i] = rsq_lookup.f; + ftable[i] = qqrd2e/r * (fgamma - 1.0); + ctable[i] = 0.0; + etable[i] = qqrd2e/r * egamma; + ptable[i] = qqrd2e/r; + vtable[i] = qqrd2e/r * fgamma; + if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { + if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + ftable[i] = qqrd2e/r * fgamma; + ctable[i] = qqrd2e/r; + } + } + } + minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); + } + + tabinnersq = minrsq_lookup.f; + + int ntablem1 = ntable - 1; + + for (int i = 0; i < ntablem1; i++) { + drtable[i] = 1.0/(rtable[i+1] - rtable[i]); + dftable[i] = ftable[i+1] - ftable[i]; + dctable[i] = ctable[i+1] - ctable[i]; + detable[i] = etable[i+1] - etable[i]; + } + + if (cut_respa) { + for (int i = 0; i < ntablem1; i++) { + dvtable[i] = vtable[i+1] - vtable[i]; + dptable[i] = ptable[i+1] - ptable[i]; + } + } + + // get the delta values for the last table entries + // tables are connected periodically between 0 and ntablem1 + + drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); + dftable[ntablem1] = ftable[0] - ftable[ntablem1]; + dctable[ntablem1] = ctable[0] - ctable[ntablem1]; + detable[ntablem1] = etable[0] - etable[ntablem1]; + if (cut_respa) { + dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; + dptable[ntablem1] = ptable[0] - ptable[ntablem1]; + } + + // get the correct delta values at itablemax + // smallest r is in bin itablemin + // largest r is in bin itablemax, which is itablemin-1, + // or ntablem1 if itablemin=0 + // deltas at itablemax only needed if corresponding rsq < cut*cut + // if so, compute deltas between rsq and cut*cut + + double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp; + itablemin = minrsq_lookup.i & ncoulmask; + itablemin >>= ncoulshiftbits; + int itablemax = itablemin - 1; + if (itablemin == 0) itablemax = ntablem1; + rsq_lookup.i = itablemax << ncoulshiftbits; + rsq_lookup.i |= maskhi; + + if (rsq_lookup.f < cut_coulsq) { + rsq_lookup.f = cut_coulsq; + r = sqrtf(rsq_lookup.f); + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + + if (cut_respa == NULL) { + f_tmp = qqrd2e/r * fgamma; + c_tmp = qqrd2e/r; + e_tmp = qqrd2e/r * egamma; + } else { + f_tmp = qqrd2e/r * (fgamma - 1.0); + c_tmp = 0.0; + e_tmp = qqrd2e/r * egamma; + p_tmp = qqrd2e/r; + v_tmp = qqrd2e/r * fgamma; + if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { + if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + f_tmp = qqrd2e/r * fgamma; + c_tmp = qqrd2e/r; + } + } + } + + drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]); + dftable[itablemax] = f_tmp - ftable[itablemax]; + dctable[itablemax] = c_tmp - ctable[itablemax]; + detable[itablemax] = e_tmp - etable[itablemax]; + if (cut_respa) { + dvtable[itablemax] = v_tmp - vtable[itablemax]; + dptable[itablemax] = p_tmp - ptable[itablemax]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +double PairCoulMSM::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r,egamma,fgamma,prefactor; + double fraction,table,forcecoul,phicoul; + int itable; + + r2inv = 1.0/rsq; + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * fgamma; + + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = atom->q[i]*atom->q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + fforce = forcecoul * r2inv; + + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*egamma; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + + return phicoul; +} + +/* ---------------------------------------------------------------------- */ + +void *PairCoulMSM::extract(const char *str, int &dim) +{ + if (strcmp(str,"cut_msm") == 0) { + dim = 0; + return (void *) &cut_coul; + } + if (strcmp(str,"scale") == 0) { + dim = 2; + return (void *) scale; + } + return NULL; +} diff --git a/src/KSPACE/pair_coul_msm.h b/src/KSPACE/pair_coul_msm.h new file mode 100644 index 0000000000..2a9bc340fd --- /dev/null +++ b/src/KSPACE/pair_coul_msm.h @@ -0,0 +1,70 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_CLASS + +PairStyle(coul/msm,PairCoulMSM) + +#else + +#ifndef LMP_PAIR_COUL_MSM_H +#define LMP_PAIR_COUL_MSM_H + +#include "pair_coul_long.h" + +namespace LAMMPS_NS { + +class PairCoulMSM : public PairCoulLong { + public: + PairCoulMSM(class LAMMPS *); + virtual ~PairCoulMSM(){}; + virtual void compute(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + virtual void *extract(const char *, int &); + + protected: + virtual void init_tables(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style lj/cut/coul/msm requires atom attribute q + +The atom style defined does not have this attribute. + +E: Pair cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +E: Pair style is incompatible with KSpace style + +If a pair style with a long-range Coulombic component is selected, +then a kspace style must also be used. + +*/ diff --git a/src/KSPACE/pair_lj_charmm_coul_long.h b/src/KSPACE/pair_lj_charmm_coul_long.h index 76c022651c..4c28126bf2 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.h +++ b/src/KSPACE/pair_lj_charmm_coul_long.h @@ -39,12 +39,12 @@ class PairLJCharmmCoulLong : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); - double single(int, int, int, int, double, double, double, double &); + virtual double single(int, int, int, int, double, double, double, double &); void compute_inner(); void compute_middle(); - void compute_outer(int, int); - void *extract(const char *, int &); + virtual void compute_outer(int, int); + virtual void *extract(const char *, int &); protected: int implicit; @@ -65,7 +65,7 @@ class PairLJCharmmCoulLong : public Pair { int ncoulshiftbits,ncoulmask; void allocate(); - void init_tables(); + virtual void init_tables(); void free_tables(); }; diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.cpp b/src/KSPACE/pair_lj_charmm_coul_msm.cpp new file mode 100644 index 0000000000..3539ded08a --- /dev/null +++ b/src/KSPACE/pair_lj_charmm_coul_msm.cpp @@ -0,0 +1,647 @@ +/* ---------------------------------------------------------------------- + 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: Paul Crozier (SNL), Stan Moore (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_charmm_coul_msm.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairLJCharmmCoulMSM::PairLJCharmmCoulMSM(LAMMPS *lmp) : PairLJCharmmCoulLong(lmp) +{ + +} + +/* ---------------------------------------------------------------------- */ +void PairLJCharmmCoulMSM::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double egamma,fgamma,prefactor; + double philj,switch1,switch2; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cut_bothsq) { + r2inv = 1.0/rsq; + + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + prefactor = qqrd2e * qtmp*q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - 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; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq) { + r6inv = r2inv*r2inv*r2inv; + jtype = type[j]; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + switch2 = 12.0*rsq * (cut_ljsq-rsq) * + (rsq-cut_lj_innersq) / denom_lj; + philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); + forcelj = forcelj*switch1 + philj*switch2; + } + } else forcelj = 0.0; + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*egamma; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + } + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + + if (rsq < cut_ljsq) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + evdwl *= switch1; + } + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double egamma,fgamma,prefactor; + double philj,switch1,switch2; + double rsw; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = listouter->inum; + ilist = listouter->ilist; + numneigh = listouter->numneigh; + firstneigh = listouter->firstneigh; + + double cut_in_off = cut_respa[2]; + double cut_in_on = cut_respa[3]; + + double cut_in_diff = cut_in_on - cut_in_off; + double cut_in_off_sq = cut_in_off*cut_in_off; + double cut_in_on_sq = cut_in_on*cut_in_on; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + 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 < cut_bothsq) { + r2inv = 1.0/rsq; + + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + prefactor = qqrd2e * qtmp*q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * (fgamma - 1.0); + if (rsq > cut_in_off_sq) { + if (rsq < cut_in_on_sq) { + rsw = (r - cut_in_off)/cut_in_diff; + forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw); + if (factor_coul < 1.0) + forcecoul -= + (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw); + } else { + forcecoul += prefactor; + if (factor_coul < 1.0) + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - 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; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq && rsq > cut_in_off_sq) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + switch2 = 12.0*rsq * (cut_ljsq-rsq) * + (rsq-cut_lj_innersq) / denom_lj; + philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); + forcelj = forcelj*switch1 + philj*switch2; + } + if (rsq < cut_in_on_sq) { + rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + forcelj *= rsw*rsw*(3.0 - 2.0*rsw); + } + } else forcelj = 0.0; + + fpair = (forcecoul + forcelj) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + ecoul = prefactor*egamma; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ptable[itable] + fraction*dptable[itable]; + prefactor = qtmp*q[j] * table; + ecoul -= (1.0-factor_coul)*prefactor; + } + } + } else ecoul = 0.0; + + if (rsq < cut_ljsq) { + r6inv = r2inv*r2inv*r2inv; + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + evdwl *= switch1; + } + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (vflag) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + table = vtable[itable] + fraction*dvtable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ptable[itable] + fraction*dptable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq <= cut_in_off_sq) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + switch2 = 12.0*rsq * (cut_ljsq-rsq) * + (rsq-cut_lj_innersq) / denom_lj; + philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); + forcelj = forcelj*switch1 + philj*switch2; + } + } else if (rsq <= cut_in_on_sq) { + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + switch2 = 12.0*rsq * (cut_ljsq-rsq) * + (rsq-cut_lj_innersq) / denom_lj; + philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); + forcelj = forcelj*switch1 + philj*switch2; + } + } + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } +} + + +/* ---------------------------------------------------------------------- + setup force tables used in compute routines +------------------------------------------------------------------------- */ + +void PairLJCharmmCoulMSM::init_tables() +{ + int masklo,maskhi; + double r,egamma,fgamma,rsw; + double qqrd2e = force->qqrd2e; + + tabinnersq = tabinner*tabinner; + init_bitmap(tabinner,cut_coul,ncoultablebits, + masklo,maskhi,ncoulmask,ncoulshiftbits); + + int ntable = 1; + for (int i = 0; i < ncoultablebits; i++) ntable *= 2; + + // linear lookup tables of length N = 2^ncoultablebits + // stored value = value at lower edge of bin + // d values = delta from lower edge to upper edge of bin + + if (ftable) free_tables(); + + memory->create(rtable,ntable,"pair:rtable"); + memory->create(ftable,ntable,"pair:ftable"); + memory->create(ctable,ntable,"pair:ctable"); + memory->create(etable,ntable,"pair:etable"); + memory->create(drtable,ntable,"pair:drtable"); + memory->create(dftable,ntable,"pair:dftable"); + memory->create(dctable,ntable,"pair:dctable"); + memory->create(detable,ntable,"pair:detable"); + + if (cut_respa == NULL) { + vtable = ptable = dvtable = dptable = NULL; + } else { + memory->create(vtable,ntable,"pair:vtable"); + memory->create(ptable,ntable,"pair:ptable"); + memory->create(dvtable,ntable,"pair:dvtable"); + memory->create(dptable,ntable,"pair:dptable"); + } + + union_int_float_t rsq_lookup; + union_int_float_t minrsq_lookup; + int itablemin; + minrsq_lookup.i = 0 << ncoulshiftbits; + minrsq_lookup.i |= maskhi; + + for (int i = 0; i < ntable; i++) { + rsq_lookup.i = i << ncoulshiftbits; + rsq_lookup.i |= masklo; + if (rsq_lookup.f < tabinnersq) { + rsq_lookup.i = i << ncoulshiftbits; + rsq_lookup.i |= maskhi; + } + r = sqrtf(rsq_lookup.f); + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + if (cut_respa == NULL) { + rtable[i] = rsq_lookup.f; + ftable[i] = qqrd2e/r * fgamma; + ctable[i] = qqrd2e/r; + etable[i] = qqrd2e/r * egamma; + } else { + rtable[i] = rsq_lookup.f; + ftable[i] = qqrd2e/r * (fgamma - 1.0); + ctable[i] = 0.0; + etable[i] = qqrd2e/r * egamma; + ptable[i] = qqrd2e/r; + vtable[i] = qqrd2e/r * fgamma; + if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { + if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + ftable[i] = qqrd2e/r * fgamma; + ctable[i] = qqrd2e/r; + } + } + } + minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); + } + + tabinnersq = minrsq_lookup.f; + + int ntablem1 = ntable - 1; + + for (int i = 0; i < ntablem1; i++) { + drtable[i] = 1.0/(rtable[i+1] - rtable[i]); + dftable[i] = ftable[i+1] - ftable[i]; + dctable[i] = ctable[i+1] - ctable[i]; + detable[i] = etable[i+1] - etable[i]; + } + + if (cut_respa) { + for (int i = 0; i < ntablem1; i++) { + dvtable[i] = vtable[i+1] - vtable[i]; + dptable[i] = ptable[i+1] - ptable[i]; + } + } + + // get the delta values for the last table entries + // tables are connected periodically between 0 and ntablem1 + + drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); + dftable[ntablem1] = ftable[0] - ftable[ntablem1]; + dctable[ntablem1] = ctable[0] - ctable[ntablem1]; + detable[ntablem1] = etable[0] - etable[ntablem1]; + if (cut_respa) { + dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; + dptable[ntablem1] = ptable[0] - ptable[ntablem1]; + } + + // get the correct delta values at itablemax + // smallest r is in bin itablemin + // largest r is in bin itablemax, which is itablemin-1, + // or ntablem1 if itablemin=0 + // deltas at itablemax only needed if corresponding rsq < cut*cut + // if so, compute deltas between rsq and cut*cut + + double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp; + itablemin = minrsq_lookup.i & ncoulmask; + itablemin >>= ncoulshiftbits; + int itablemax = itablemin - 1; + if (itablemin == 0) itablemax = ntablem1; + rsq_lookup.i = itablemax << ncoulshiftbits; + rsq_lookup.i |= maskhi; + + if (rsq_lookup.f < cut_coulsq) { + rsq_lookup.f = cut_coulsq; + r = sqrtf(rsq_lookup.f); + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + + if (cut_respa == NULL) { + f_tmp = qqrd2e/r * fgamma; + c_tmp = qqrd2e/r; + e_tmp = qqrd2e/r * egamma; + } else { + f_tmp = qqrd2e/r * (fgamma - 1.0); + c_tmp = 0.0; + e_tmp = qqrd2e/r * egamma; + p_tmp = qqrd2e/r; + v_tmp = qqrd2e/r * fgamma; + if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { + if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + f_tmp = qqrd2e/r * fgamma; + c_tmp = qqrd2e/r; + } + } + } + + drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]); + dftable[itablemax] = f_tmp - ftable[itablemax]; + dctable[itablemax] = c_tmp - ctable[itablemax]; + detable[itablemax] = e_tmp - etable[itablemax]; + if (cut_respa) { + dvtable[itablemax] = v_tmp - vtable[itablemax]; + dptable[itablemax] = p_tmp - ptable[itablemax]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,egamma,fgamma,prefactor; + double switch1,switch2,fraction,table,forcecoul,forcelj,phicoul,philj; + int itable; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = atom->q[i]*atom->q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + if (rsq < cut_ljsq) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + switch2 = 12.0*rsq * (cut_ljsq-rsq) * + (rsq-cut_lj_innersq) / denom_lj; + philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); + forcelj = forcelj*switch1 + philj*switch2; + } + } else forcelj = 0.0; + fforce = (forcecoul + factor_lj*forcelj) * r2inv; + + double eng = 0.0; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*egamma; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + + if (rsq < cut_ljsq) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + philj *= switch1; + } + eng += factor_lj*philj; + } + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJCharmmCoulMSM::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"lj14_1") == 0) return (void *) lj14_1; + if (strcmp(str,"lj14_2") == 0) return (void *) lj14_2; + if (strcmp(str,"lj14_3") == 0) return (void *) lj14_3; + if (strcmp(str,"lj14_4") == 0) return (void *) lj14_4; + + dim = 0; + if (strcmp(str,"implicit") == 0) return (void *) &implicit; + if (strcmp(str,"cut_msm") == 0) return (void *) &cut_coul; + + return NULL; +} diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.h b/src/KSPACE/pair_lj_charmm_coul_msm.h new file mode 100644 index 0000000000..aa437ef51e --- /dev/null +++ b/src/KSPACE/pair_lj_charmm_coul_msm.h @@ -0,0 +1,80 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_CLASS + +PairStyle(lj/charmm/coul/msm,PairLJCharmmCoulMSM) + +#else + +#ifndef LMP_PAIR_LJ_CHARMM_COUL_MSM_H +#define LMP_PAIR_LJ_CHARMM_COUL_MSM_H + +#include "pair_lj_charmm_coul_long.h" + +namespace LAMMPS_NS { + +class PairLJCharmmCoulMSM : public PairLJCharmmCoulLong { + public: + PairLJCharmmCoulMSM(class LAMMPS *); + virtual ~PairLJCharmmCoulMSM(){}; + virtual void compute(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + virtual void compute_outer(int, int); + virtual void *extract(const char *, int &); + + protected: + virtual void init_tables(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style lj/charmm/coul/long requires atom attribute q + +The atom style defined does not have these attributes. + +E: Pair inner cutoff >= Pair outer cutoff + +The specified cutoffs for the pair style are inconsistent. + +E: Pair cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +E: Pair inner cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +E: Pair style is incompatible with KSpace style + +If a pair style with a long-range Coulombic component is selected, +then a kspace style must also be used. + +*/ diff --git a/src/KSPACE/pair_lj_cut_coul_long.h b/src/KSPACE/pair_lj_cut_coul_long.h index fa8acedcda..6f65573a50 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.h +++ b/src/KSPACE/pair_lj_cut_coul_long.h @@ -25,6 +25,7 @@ PairStyle(lj/cut/coul/long,PairLJCutCoulLong) namespace LAMMPS_NS { class PairLJCutCoulLong : public Pair { + public: PairLJCutCoulLong(class LAMMPS *); virtual ~PairLJCutCoulLong(); @@ -42,8 +43,8 @@ class PairLJCutCoulLong : public Pair { void compute_inner(); void compute_middle(); - void compute_outer(int, int); - void *extract(const char *, int &); + virtual void compute_outer(int, int); + virtual void *extract(const char *, int &); protected: double cut_lj_global; @@ -61,7 +62,7 @@ class PairLJCutCoulLong : public Pair { int ncoulshiftbits,ncoulmask; void allocate(); - void init_tables(); + virtual void init_tables(); void free_tables(); }; diff --git a/src/KSPACE/pair_lj_cut_coul_msm.cpp b/src/KSPACE/pair_lj_cut_coul_msm.cpp new file mode 100644 index 0000000000..087a90d384 --- /dev/null +++ b/src/KSPACE/pair_lj_cut_coul_msm.cpp @@ -0,0 +1,590 @@ +/* ---------------------------------------------------------------------- + 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: Stan Moore (SNL), Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_cut_coul_msm.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulMSM::PairLJCutCoulMSM(LAMMPS *lmp) : PairLJCutCoulLong(lmp) +{ + +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulMSM::compute(int eflag, int vflag) +{ + int i,ii,j,jj,inum,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double egamma,fgamma,prefactor; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + 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_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + prefactor = qqrd2e * qtmp*q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - 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; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else forcelj = 0.0; + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*egamma; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + } + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulMSM::compute_outer(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double egamma,fgamma,prefactor; + double rsw; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = listouter->inum; + ilist = listouter->ilist; + numneigh = listouter->numneigh; + firstneigh = listouter->firstneigh; + + double cut_in_off = cut_respa[2]; + double cut_in_on = cut_respa[3]; + + double cut_in_diff = cut_in_on - cut_in_off; + double cut_in_off_sq = cut_in_off*cut_in_off; + double cut_in_on_sq = cut_in_on*cut_in_on; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + 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_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + prefactor = qqrd2e * qtmp*q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * (fgamma - 1.0); + if (rsq > cut_in_off_sq) { + if (rsq < cut_in_on_sq) { + rsw = (r - cut_in_off)/cut_in_diff; + forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw); + if (factor_coul < 1.0) + forcecoul -= + (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw); + } else { + forcecoul += prefactor; + if (factor_coul < 1.0) + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - 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; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq < cut_in_on_sq) { + rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + forcelj *= rsw*rsw*(3.0 - 2.0*rsw); + } + } else forcelj = 0.0; + + fpair = (forcecoul + forcelj) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + ecoul = prefactor*egamma; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ptable[itable] + fraction*dptable[itable]; + prefactor = qtmp*q[j] * table; + ecoul -= (1.0-factor_coul)*prefactor; + } + } + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (vflag) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + table = vtable[itable] + fraction*dvtable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ptable[itable] + fraction*dptable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq <= cut_in_off_sq) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else if (rsq <= cut_in_on_sq) + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } +} +/* ---------------------------------------------------------------------- + setup force tables used in compute routines +------------------------------------------------------------------------- */ + +void PairLJCutCoulMSM::init_tables() +{ + int masklo,maskhi; + double r,egamma,fgamma,rsw; + double qqrd2e = force->qqrd2e; + + tabinnersq = tabinner*tabinner; + init_bitmap(tabinner,cut_coul,ncoultablebits, + masklo,maskhi,ncoulmask,ncoulshiftbits); + + int ntable = 1; + for (int i = 0; i < ncoultablebits; i++) ntable *= 2; + + // linear lookup tables of length N = 2^ncoultablebits + // stored value = value at lower edge of bin + // d values = delta from lower edge to upper edge of bin + + if (ftable) free_tables(); + + memory->create(rtable,ntable,"pair:rtable"); + memory->create(ftable,ntable,"pair:ftable"); + memory->create(ctable,ntable,"pair:ctable"); + memory->create(etable,ntable,"pair:etable"); + memory->create(drtable,ntable,"pair:drtable"); + memory->create(dftable,ntable,"pair:dftable"); + memory->create(dctable,ntable,"pair:dctable"); + memory->create(detable,ntable,"pair:detable"); + + if (cut_respa == NULL) { + vtable = ptable = dvtable = dptable = NULL; + } else { + memory->create(vtable,ntable*sizeof(double),"pair:vtable"); + memory->create(ptable,ntable*sizeof(double),"pair:ptable"); + memory->create(dvtable,ntable*sizeof(double),"pair:dvtable"); + memory->create(dptable,ntable*sizeof(double),"pair:dptable"); + } + + union_int_float_t rsq_lookup; + union_int_float_t minrsq_lookup; + int itablemin; + minrsq_lookup.i = 0 << ncoulshiftbits; + minrsq_lookup.i |= maskhi; + + for (int i = 0; i < ntable; i++) { + rsq_lookup.i = i << ncoulshiftbits; + rsq_lookup.i |= masklo; + if (rsq_lookup.f < tabinnersq) { + rsq_lookup.i = i << ncoulshiftbits; + rsq_lookup.i |= maskhi; + } + r = sqrtf(rsq_lookup.f); + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + + if (cut_respa == NULL) { + rtable[i] = rsq_lookup.f; + ftable[i] = qqrd2e/r * fgamma; + ctable[i] = qqrd2e/r; + etable[i] = qqrd2e/r * egamma; + } else { + rtable[i] = rsq_lookup.f; + ftable[i] = qqrd2e/r * (fgamma - 1.0); + ctable[i] = 0.0; + etable[i] = qqrd2e/r * egamma; + ptable[i] = qqrd2e/r; + vtable[i] = qqrd2e/r * fgamma; + if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { + if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + ftable[i] = qqrd2e/r * fgamma; + ctable[i] = qqrd2e/r; + } + } + } + minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); + } + + tabinnersq = minrsq_lookup.f; + + int ntablem1 = ntable - 1; + + for (int i = 0; i < ntablem1; i++) { + drtable[i] = 1.0/(rtable[i+1] - rtable[i]); + dftable[i] = ftable[i+1] - ftable[i]; + dctable[i] = ctable[i+1] - ctable[i]; + detable[i] = etable[i+1] - etable[i]; + } + + if (cut_respa) { + for (int i = 0; i < ntablem1; i++) { + dvtable[i] = vtable[i+1] - vtable[i]; + dptable[i] = ptable[i+1] - ptable[i]; + } + } + + // get the delta values for the last table entries + // tables are connected periodically between 0 and ntablem1 + + drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); + dftable[ntablem1] = ftable[0] - ftable[ntablem1]; + dctable[ntablem1] = ctable[0] - ctable[ntablem1]; + detable[ntablem1] = etable[0] - etable[ntablem1]; + if (cut_respa) { + dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; + dptable[ntablem1] = ptable[0] - ptable[ntablem1]; + } + + // get the correct delta values at itablemax + // smallest r is in bin itablemin + // largest r is in bin itablemax, which is itablemin-1, + // or ntablem1 if itablemin=0 + // deltas at itablemax only needed if corresponding rsq < cut*cut + // if so, compute deltas between rsq and cut*cut + + double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp; + itablemin = minrsq_lookup.i & ncoulmask; + itablemin >>= ncoulshiftbits; + int itablemax = itablemin - 1; + if (itablemin == 0) itablemax = ntablem1; + rsq_lookup.i = itablemax << ncoulshiftbits; + rsq_lookup.i |= maskhi; + + if (rsq_lookup.f < cut_coulsq) { + rsq_lookup.f = cut_coulsq; + r = sqrtf(rsq_lookup.f); + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + + if (cut_respa == NULL) { + f_tmp = qqrd2e/r * fgamma; + c_tmp = qqrd2e/r; + e_tmp = qqrd2e/r * egamma; + } else { + f_tmp = qqrd2e/r * (fgamma - 1.0); + c_tmp = 0.0; + e_tmp = qqrd2e/r * egamma; + p_tmp = qqrd2e/r; + v_tmp = qqrd2e/r * fgamma; + if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { + if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + f_tmp = qqrd2e/r * fgamma; + c_tmp = qqrd2e/r; + } + } + } + + drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]); + dftable[itablemax] = f_tmp - ftable[itablemax]; + dctable[itablemax] = c_tmp - ctable[itablemax]; + detable[itablemax] = e_tmp - etable[itablemax]; + if (cut_respa) { + dvtable[itablemax] = v_tmp - vtable[itablemax]; + dptable[itablemax] = p_tmp - ptable[itablemax]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCutCoulMSM::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,egamma,fgamma,prefactor; + double fraction,table,forcecoul,forcelj,phicoul,philj; + int itable; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); + fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); + forcecoul = prefactor * fgamma; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup_single; + rsq_lookup_single.f = rsq; + itable = rsq_lookup_single.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = atom->q[i]*atom->q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else forcelj = 0.0; + + fforce = (forcecoul + factor_lj*forcelj) * r2inv; + + double eng = 0.0; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*egamma; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + + if (rsq < cut_ljsq[itype][jtype]) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + eng += factor_lj*philj; + } + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJCutCoulMSM::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_msm") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + return NULL; +} diff --git a/src/KSPACE/pair_lj_cut_coul_msm.h b/src/KSPACE/pair_lj_cut_coul_msm.h new file mode 100644 index 0000000000..6e69c3f727 --- /dev/null +++ b/src/KSPACE/pair_lj_cut_coul_msm.h @@ -0,0 +1,71 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_CLASS + +PairStyle(lj/cut/coul/msm,PairLJCutCoulMSM) + +#else + +#ifndef LMP_PAIR_LJ_CUT_COUL_MSM_H +#define LMP_PAIR_LJ_CUT_COUL_MSM_H + +#include "pair_lj_cut_coul_long.h" + +namespace LAMMPS_NS { + +class PairLJCutCoulMSM : public PairLJCutCoulLong { + public: + PairLJCutCoulMSM(class LAMMPS *); + virtual ~PairLJCutCoulMSM(){}; + virtual void compute(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + virtual void compute_outer(int, int); + virtual void *extract(const char *, int &); + + protected: + virtual void init_tables(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style lj/cut/coul/msm requires atom attribute q + +The atom style defined does not have this attribute. + +E: Pair style is incompatible with KSpace style + +If a pair style with a long-range Coulombic component is selected, +then a kspace style must also be used. + +E: Pair cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +*/