diff --git a/src/CORESHELL/Install.sh b/src/CORESHELL/Install.sh index e6650adcf2..b89b586865 100644 --- a/src/CORESHELL/Install.sh +++ b/src/CORESHELL/Install.sh @@ -30,3 +30,5 @@ action pair_born_coul_long_cs.cpp pair_born_coul_long.cpp action pair_buck_coul_long_cs.cpp pair_buck_coul_long.cpp action pair_born_coul_long_cs.h pair_born_coul_long.h action pair_buck_coul_long_cs.h pair_buck_coul_long.h +action pair_coul_long_cs.cpp pair_coul_long.cpp +action pair_coul_long_cs.h pair_coul_long.h diff --git a/src/CORESHELL/pair_coul_long_cs.cpp b/src/CORESHELL/pair_coul_long_cs.cpp new file mode 100644 index 0000000000..fc561ac219 --- /dev/null +++ b/src/CORESHELL/pair_coul_long_cs.cpp @@ -0,0 +1,159 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Added epsilon for use with CORESHELL and USER-DRUDE +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_coul_long_cs.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; + +#define EWALD_F 1.12837917 +#define EWALD_P 9.95473818e-1 +#define B0 -0.1335096380159268 +#define B1 -2.57839507e-1 +#define B2 -1.37203639e-1 +#define B3 -8.88822059e-3 +#define B4 -5.80844129e-3 +#define B5 1.14652755e-1 + +/* ---------------------------------------------------------------------- */ + +PairCoulLongCS::PairCoulLongCS(LAMMPS *lmp) : PairCoulLong(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void PairCoulLongCS::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 grij,expm2,prefactor,t,erfc,u; + 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; + r = sqrt(rsq); + if (!ncoultablebits || rsq <= tabinnersq) { + prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]; + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + u = 1. - t; + erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2; + prefactor /= r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + 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*erfc; + 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(); +} diff --git a/src/CORESHELL/pair_coul_long_cs.h b/src/CORESHELL/pair_coul_long_cs.h new file mode 100644 index 0000000000..a6dcaf4db8 --- /dev/null +++ b/src/CORESHELL/pair_coul_long_cs.h @@ -0,0 +1,58 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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/long/cs,PairCoulLongCS) + +#else + +#ifndef LMP_PAIR_COUL_LONG_CS_H +#define LMP_PAIR_COUL_LONG_CS_H + +#include "pair_coul_long.h" + +namespace LAMMPS_NS { + +class PairCoulLongCS : public PairCoulLong { + public: + PairCoulLongCS(class LAMMPS *); + virtual void compute(int, 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: Pair style lj/cut/coul/long/drude requires atom attribute q + +The atom style defined does not have this attribute. + +E: Pair style requires a KSpace style + +No kspace style is defined. + +*/ diff --git a/src/Makefile b/src/Makefile index aa692caa71..88b7b09b7f 100755 --- a/src/Makefile +++ b/src/Makefile @@ -46,9 +46,9 @@ PACKAGE = asphere body class2 colloid coreshell dipole fld gpu granular kim \ python qeq reax replica rigid shock snap srd voronoi xtc PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars user-cuda \ - user-diffraction user-eff user-fep user-intel user-lb user-misc \ - user-molfile user-omp user-phonon user-qmmm user-qtb user-quip \ - user-reaxc user-sph + user-diffraction user-drude user-eff user-fep user-intel user-lb \ + user-misc user-molfile user-omp user-phonon user-qmmm user-qtb \ + user-quip user-reaxc user-sph PACKLIB = gpu kim kokkos meam poems python reax voronoi \ user-atc user-awpmd user-colvars user-cuda user-molfile \ diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp index 62dcea023d..1dbd21726f 100644 --- a/src/compute_chunk_atom.cpp +++ b/src/compute_chunk_atom.cpp @@ -722,12 +722,17 @@ int ComputeChunkAtom::setup_chunks() // assign chunk IDs to atoms // will exclude atoms not in group or in optional region // for binning styles, need to setup bins and their volume first + // else chunk_volume_scalar = box volume // IDs are needed to scan for max ID and for compress() if (binflag) { nchunk = setup_bins(); bin_volumes(); + } else { + chunk_volume_scalar = domain->xprd * domain->yprd; + if (domain->dimension == 3) chunk_volume_scalar *= domain->zprd; } + assign_chunk_ids(); // set nchunk for chunk styles other than binning diff --git a/src/fix.cpp b/src/fix.cpp index 64b84c38d2..ac7e014393 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -66,6 +66,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) wd_header = wd_section = 0; dynamic_group_allow = 0; dof_flag = 0; + special_alter_flag = 0; cudable_comm = 0; scalar_flag = vector_flag = array_flag = 0; diff --git a/src/fix.h b/src/fix.h index 6a11e18844..f0836de021 100644 --- a/src/fix.h +++ b/src/fix.h @@ -50,6 +50,7 @@ class Fix : protected Pointers { int wd_section; // # of sections fix writes to data file int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 int dof_flag; // 1 if has dof() method (not min_dof()) + int special_alter_flag; // 1 if has special_alter() meth for spec lists int cudable_comm; // 1 if fix has CUDA-enabled communication int scalar_flag; // 0/1 if compute_scalar() function exists @@ -194,6 +195,8 @@ class Fix : protected Pointers { virtual void zero_momentum() {} virtual void zero_rotation() {} + virtual void rebuild_special() {} + virtual int modify_param(int, char **) {return 0;} virtual void *extract(const char *, int &) {return NULL;} diff --git a/src/neighbor.cpp b/src/neighbor.cpp index ade7dc956c..5fc91aee04 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -310,6 +310,7 @@ void Neighbor::init() // flag = 1 if both LJ/Coulomb special values are 1.0 // flag = 2 otherwise or if KSpace solver is enabled // pairwise portion of KSpace solver uses all 1-2,1-3,1-4 neighbors + // or selected Coulomb-approixmation pair styles require it if (force->special_lj[1] == 0.0 && force->special_coul[1] == 0.0) special_flag[1] = 0; @@ -329,8 +330,8 @@ void Neighbor::init() special_flag[3] = 1; else special_flag[3] = 2; - if (force->kspace || force->pair_match("coul/wolf",0) - || force->pair_match("coul/dsf",0)) + if (force->kspace || force->pair_match("coul/wolf",0) || + force->pair_match("coul/dsf",0) || force->pair_match("thole",0)) special_flag[1] = special_flag[2] = special_flag[3] = 2; // maxwt = max multiplicative factor on atom indices stored in neigh list diff --git a/src/special.cpp b/src/special.cpp index e685a54a27..ebabbff281 100644 --- a/src/special.cpp +++ b/src/special.cpp @@ -18,6 +18,8 @@ #include "atom_vec.h" #include "force.h" #include "comm.h" +#include "modify.h" +#include "fix.h" #include "accelerator_kokkos.h" #include "memory.h" #include "error.h" @@ -191,6 +193,7 @@ void Special::build() force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0) { dedup(); combine(); + fix_alteration(); return; } @@ -307,6 +310,7 @@ void Special::build() dedup(); if (force->special_angle) angle_trim(); combine(); + fix_alteration(); return; } @@ -419,6 +423,7 @@ void Special::build() if (force->special_angle) angle_trim(); if (force->special_dihedral) dihedral_trim(); combine(); + fix_alteration(); } /* ---------------------------------------------------------------------- @@ -1127,3 +1132,17 @@ void Special::ring_eight(int ndatum, char *cbuf) i += 2; } } + +/* ---------------------------------------------------------------------- + allow fixes to alter special list + currently, only fix drude does this + so that both the Drude core and electron are same level of neighbor +------------------------------------------------------------------------- */ + +void Special::fix_alteration() +{ + for (int ifix = 0; ifix < modify->nfix; ifix++) + if (modify->fix[ifix]->special_alter_flag) + modify->fix[ifix]->rebuild_special(); +} + diff --git a/src/special.h b/src/special.h index 987ca490ea..fbaad5ea34 100644 --- a/src/special.h +++ b/src/special.h @@ -37,6 +37,7 @@ class Special : protected Pointers { void angle_trim(); void dihedral_trim(); void combine(); + void fix_alteration(); // static variable for ring communication callback to access class data // callback functions for ring communication