Merge branch 'master' into next_lammps_version

This commit is contained in:
Axel Kohlmeyer 2020-04-15 05:08:53 -04:00
commit 3deece168e
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
8 changed files with 692 additions and 57 deletions

View File

@ -107,7 +107,7 @@ install(TARGETS lmp EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR})
option(CMAKE_VERBOSE_MAKEFILE "Generate verbose Makefiles" OFF)
set(DEFAULT_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS CORESHELL DIPOLE
set(STANDARD_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE
GRANULAR KSPACE LATTE MANYBODY MC MESSAGE MISC MOLECULE PERI POEMS QEQ
REPLICA RIGID SHOCK SPIN SNAP SRD KIM PYTHON MSCG MPIIO VORONOI
USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-MESODPD USER-CGSDK USER-COLVARS
@ -116,8 +116,8 @@ set(DEFAULT_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS CORESHELL DIPOLE
USER-NETCDF USER-PHONON USER-PLUMED USER-PTM USER-QTB USER-REACTION
USER-REAXC USER-SCAFACOS USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY
USER-UEF USER-VTK USER-QUIP USER-QMMM USER-YAFF USER-ADIOS)
set(ACCEL_PACKAGES USER-OMP KOKKOS OPT USER-INTEL GPU)
foreach(PKG ${DEFAULT_PACKAGES} ${ACCEL_PACKAGES})
set(SUFFIX_PACKAGES CORESHELL USER-OMP KOKKOS OPT USER-INTEL GPU)
foreach(PKG ${STANDARD_PACKAGES} ${SUFFIX_PACKAGES})
option(PKG_${PKG} "Build ${PKG} Package" OFF)
endforeach()
@ -362,7 +362,7 @@ RegisterStyles(${LAMMPS_SOURCE_DIR})
##############################################
# add sources of enabled packages
############################################
foreach(PKG ${DEFAULT_PACKAGES})
foreach(PKG ${STANDARD_PACKAGES})
set(${PKG}_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/${PKG})
file(GLOB ${PKG}_SOURCES ${${PKG}_SOURCES_DIR}/[^.]*.cpp)
@ -390,7 +390,7 @@ foreach(PKG MPIIO)
endforeach()
# dedicated check for entire contents of accelerator packages
foreach(PKG ${ACCEL_PACKAGES})
foreach(PKG ${SUFFIX_PACKAGES})
set(${PKG}_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/${PKG})
file(GLOB ${PKG}_SOURCES ${${PKG}_SOURCES_DIR}/[^.]*.cpp)
@ -481,7 +481,7 @@ target_include_directories(lammps PRIVATE ${LAMMPS_STYLE_HEADERS_DIR})
######################################
set(temp "#ifndef LMP_INSTALLED_PKGS_H\n#define LMP_INSTALLED_PKGS_H\n")
set(temp "${temp}const char * LAMMPS_NS::LAMMPS::installed_packages[] = {\n")
set(temp_PKG_LIST ${DEFAULT_PACKAGES} ${ACCEL_PACKAGES})
set(temp_PKG_LIST ${STANDARD_PACKAGES} ${SUFFIX_PACKAGES})
list(SORT temp_PKG_LIST)
foreach(PKG ${temp_PKG_LIST})
if(PKG_${PKG})
@ -655,7 +655,7 @@ include(CodeCoverage)
###############################################################################
# Print package summary
###############################################################################
foreach(PKG ${DEFAULT_PACKAGES} ${ACCEL_PACKAGES})
foreach(PKG ${STANDARD_PACKAGES} ${SUFFIX_PACKAGES})
if(PKG_${PKG})
message(STATUS "Building package: ${PKG}")
endif()

View File

@ -131,6 +131,7 @@ OPT.
* :doc:`lj/class2/coul/cut (ko) <pair_class2>`
* :doc:`lj/class2/coul/cut/soft <pair_fep_soft>`
* :doc:`lj/class2/coul/long (gko) <pair_class2>`
* :doc:`lj/class2/coul/long/cs <pair_cs>`
* :doc:`lj/class2/coul/long/soft <pair_fep_soft>`
* :doc:`lj/class2/soft <pair_fep_soft>`
* :doc:`lj/cubic (go) <pair_lj_cubic>`

View File

@ -30,6 +30,9 @@ pair_style coul/wolf/cs command
pair_style lj/cut/coul/long/cs command
=======================================
pair_style lj/class2/coul/long/cs command
==========================================
Syntax
""""""
@ -37,7 +40,7 @@ Syntax
pair_style style args
* style = *born/coul/dsf/cs* or *born/coul/long/cs* or *born/coul/wolf/cs* or *buck/coul/long/cs* or *coul/long/cs* or *coul/wolf/cs* or *lj/cut/coul/long/cs*
* style = *born/coul/dsf/cs* or *born/coul/long/cs* or *born/coul/wolf/cs* or *buck/coul/long/cs* or *coul/long/cs* or *coul/wolf/cs* or *lj/cut/coul/long/cs* or *lj/class2/coul/long/cs*
* args = list of arguments for a particular style
.. parsed-literal::
@ -64,6 +67,9 @@ Syntax
*lj/cut/coul/long/cs* args = cutoff (cutoff2)
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
*lj/class2/coul/long/cs* args = cutoff (cutoff2)
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
Examples
""""""""
@ -115,6 +121,7 @@ the "/cs" in the name:
* :doc:`pair_style coul/long <pair_coul>`
* :doc:`pair_style coul/wolf <pair_coul>`
* :doc:`pair_style lj/cut/coul/long <pair_lj>`
* :doc:`pair_style lj/class2/coul/long <pair_class2>`
except that they correctly treat the special case where the distance
between two charged core and shell atoms in the same core/shell pair

View File

@ -111,41 +111,41 @@ accelerated styles exist.
* :doc:`born <pair_born>` - Born-Mayer-Huggins potential
* :doc:`born/coul/dsf <pair_born>` - Born with damped-shifted-force model
* :doc:`born/coul/dsf/cs <pair_cs>` - Born with damped-shifted-force and core/shell model
* :doc:`born/coul/long <pair_born>` - Born with long-range Coulombics
* :doc:`born/coul/long/cs <pair_cs>` - Born with long-range Coulombics and core/shell
* :doc:`born/coul/msm <pair_born>` - Born with long-range MSM Coulombics
* :doc:`born/coul/wolf <pair_born>` - Born with Wolf potential for Coulombics
* :doc:`born/coul/wolf/cs <pair_cs>` - Born with Wolf potential for Coulombics and core/shell model
* :doc:`born/coul/long <pair_born>` - Born with long-range Coulomb
* :doc:`born/coul/long/cs <pair_cs>` - Born with long-range Coulomb and core/shell
* :doc:`born/coul/msm <pair_born>` - Born with long-range MSM Coulomb
* :doc:`born/coul/wolf <pair_born>` - Born with Wolf potential for Coulomb
* :doc:`born/coul/wolf/cs <pair_cs>` - Born with Wolf potential for Coulomb and core/shell model
* :doc:`brownian <pair_brownian>` - Brownian potential for Fast Lubrication Dynamics
* :doc:`brownian/poly <pair_brownian>` - Brownian potential for Fast Lubrication Dynamics with polydispersity
* :doc:`buck <pair_buck>` - Buckingham potential
* :doc:`buck/coul/cut <pair_buck>` - Buckingham with cutoff Coulomb
* :doc:`buck/coul/long <pair_buck>` - Buckingham with long-range Coulombics
* :doc:`buck/coul/long/cs <pair_cs>` - Buckingham with long-range Coulombics and core/shell
* :doc:`buck/coul/msm <pair_buck>` - Buckingham with long-range MSM Coulombics
* :doc:`buck/long/coul/long <pair_buck_long>` - long-range Buckingham with long-range Coulombics
* :doc:`buck/coul/long <pair_buck>` - Buckingham with long-range Coulomb
* :doc:`buck/coul/long/cs <pair_cs>` - Buckingham with long-range Coulomb and core/shell
* :doc:`buck/coul/msm <pair_buck>` - Buckingham with long-range MSM Coulomb
* :doc:`buck/long/coul/long <pair_buck_long>` - long-range Buckingham with long-range Coulomb
* :doc:`buck/mdf <pair_mdf>` - Buckingham with a taper function
* :doc:`buck6d/coul/gauss/dsf <pair_buck6d_coul_gauss>` - dispersion-damped Buckingham with damped-shift-force model
* :doc:`buck6d/coul/gauss/long <pair_buck6d_coul_gauss>` - dispersion-damped Buckingham with long-range Coulombics
* :doc:`buck6d/coul/gauss/long <pair_buck6d_coul_gauss>` - dispersion-damped Buckingham with long-range Coulomb
* :doc:`colloid <pair_colloid>` - integrated colloidal potential
* :doc:`comb <pair_comb>` - charge-optimized many-body (COMB) potential
* :doc:`comb3 <pair_comb>` - charge-optimized many-body (COMB3) potential
* :doc:`cosine/squared <pair_cosine_squared>` - Cooke-Kremer-Deserno membrane model potential
* :doc:`coul/cut <pair_coul>` - cutoff Coulombic potential
* :doc:`coul/cut/soft <pair_fep_soft>` - Coulombic potential with a soft core
* :doc:`coul/debye <pair_coul>` - cutoff Coulombic potential with Debye screening
* :doc:`coul/cut <pair_coul>` - cutoff Coulomb potential
* :doc:`coul/cut/soft <pair_fep_soft>` - Coulomb potential with a soft core
* :doc:`coul/debye <pair_coul>` - cutoff Coulomb potential with Debye screening
* :doc:`coul/diel <pair_coul_diel>` - Coulomb potential with dielectric permittivity
* :doc:`coul/dsf <pair_coul>` - Coulombics with damped-shifted-force model
* :doc:`coul/long <pair_coul>` - long-range Coulombic potential
* :doc:`coul/long/cs <pair_cs>` - long-range Coulombic potential and core/shell
* :doc:`coul/long/soft <pair_fep_soft>` - long-range Coulombic potential with a soft core
* :doc:`coul/msm <pair_coul>` - long-range MSM Coulombics
* :doc:`coul/slater/cut <pair_coul>` - smeared out Coulombics
* :doc:`coul/slater/long <pair_coul>` - long-range smeared out Coulombics
* :doc:`coul/shield <pair_coul_shield>` - Coulombics for boron nitride for use with :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>` potential
* :doc:`coul/streitz <pair_coul>` - Coulombics via Streitz/Mintmire Slater orbitals
* :doc:`coul/wolf <pair_coul>` - Coulombics via Wolf potential
* :doc:`coul/wolf/cs <pair_cs>` - ditto with core/shell adjustments
* :doc:`coul/dsf <pair_coul>` - Coulomb with damped-shifted-force model
* :doc:`coul/long <pair_coul>` - long-range Coulomb potential
* :doc:`coul/long/cs <pair_cs>` - long-range Coulomb potential and core/shell
* :doc:`coul/long/soft <pair_fep_soft>` - long-range Coulomb potential with a soft core
* :doc:`coul/msm <pair_coul>` - long-range MSM Coulomb
* :doc:`coul/slater/cut <pair_coul>` - smeared out Coulomb
* :doc:`coul/slater/long <pair_coul>` - long-range smeared out Coulomb
* :doc:`coul/shield <pair_coul_shield>` - Coulomb for boron nitride for use with :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>` potential
* :doc:`coul/streitz <pair_coul>` - Coulomb via Streitz/Mintmire Slater orbitals
* :doc:`coul/wolf <pair_coul>` - Coulomb via Wolf potential
* :doc:`coul/wolf/cs <pair_cs>` - Coulomb via Wolf potential with core/shell adjustments
* :doc:`dpd <pair_dpd>` - dissipative particle dynamics (DPD)
* :doc:`dpd/fdt <pair_dpd_fdt>` - DPD for constant temperature and pressure
* :doc:`dpd/fdt/energy <pair_dpd_fdt>` - DPD for constant energy and enthalpy
@ -189,44 +189,45 @@ accelerated styles exist.
* :doc:`lj/charmm/coul/charmm/implicit <pair_charmm>` - CHARMM for implicit solvent
* :doc:`lj/charmm/coul/long <pair_charmm>` - CHARMM with long-range Coulomb
* :doc:`lj/charmm/coul/long/soft <pair_fep_soft>` - CHARMM with long-range Coulomb and a soft core
* :doc:`lj/charmm/coul/msm <pair_charmm>` - CHARMM with long-range MSM Coulombics
* :doc:`lj/charmm/coul/msm <pair_charmm>` - CHARMM with long-range MSM Coulomb
* :doc:`lj/charmmfsw/coul/charmmfsh <pair_charmm>` - CHARMM with force switching and shifting
* :doc:`lj/charmmfsw/coul/long <pair_charmm>` - CHARMM with force switching and long-rnage Coulombics
* :doc:`lj/class2 <pair_class2>` - COMPASS (class 2) force field with no Coulomb
* :doc:`lj/charmmfsw/coul/long <pair_charmm>` - CHARMM with force switching and long-rnage Coulomb
* :doc:`lj/class2 <pair_class2>` - COMPASS (class 2) force field without Coulomb
* :doc:`lj/class2/coul/cut <pair_class2>` - COMPASS with cutoff Coulomb
* :doc:`lj/class2/coul/cut/soft <pair_fep_soft>` - COMPASS with cutoff Coulomb with a soft core
* :doc:`lj/class2/coul/long <pair_class2>` - COMPASS with long-range Coulomb
* :doc:`lj/class2/coul/long/cs <pair_cs>` - COMPASS with long-range Coulomb with core/shell adjustments
* :doc:`lj/class2/coul/long/soft <pair_fep_soft>` - COMPASS with long-range Coulomb with a soft core
* :doc:`lj/class2/soft <pair_fep_soft>` - COMPASS (class 2) force field with no Coulomb with a soft core
* :doc:`lj/cubic <pair_lj_cubic>` - LJ with cubic after inflection point
* :doc:`lj/cut <pair_lj>` - cutoff Lennard-Jones potential with no Coulomb
* :doc:`lj/cut <pair_lj>` - cutoff Lennard-Jones potential without Coulomb
* :doc:`lj/cut/coul/cut <pair_lj>` - LJ with cutoff Coulomb
* :doc:`lj/cut/coul/cut/soft <pair_fep_soft>` - LJ with cutoff Coulomb with a soft core
* :doc:`lj/cut/coul/debye <pair_lj>` - LJ with Debye screening added to Coulomb
* :doc:`lj/cut/coul/dsf <pair_lj>` - LJ with Coulombics via damped shifted forces
* :doc:`lj/cut/coul/long <pair_lj>` - LJ with long-range Coulombics
* :doc:`lj/cut/coul/long/cs <pair_cs>` - ditto with core/shell adjustments
* :doc:`lj/cut/coul/long/soft <pair_fep_soft>` - LJ with long-range Coulombics with a soft core
* :doc:`lj/cut/coul/msm <pair_lj>` - LJ with long-range MSM Coulombics
* :doc:`lj/cut/coul/wolf <pair_lj>` - LJ with Coulombics via Wolf potential
* :doc:`lj/cut/coul/dsf <pair_lj>` - LJ with Coulomb via damped shifted forces
* :doc:`lj/cut/coul/long <pair_lj>` - LJ with long-range Coulomb
* :doc:`lj/cut/coul/long/cs <pair_cs>` - LJ with long-range Coulomb with core/shell adjustments
* :doc:`lj/cut/coul/long/soft <pair_fep_soft>` - LJ with long-range Coulomb with a soft core
* :doc:`lj/cut/coul/msm <pair_lj>` - LJ with long-range MSM Coulomb
* :doc:`lj/cut/coul/wolf <pair_lj>` - LJ with Coulomb via Wolf potential
* :doc:`lj/cut/dipole/cut <pair_dipole>` - point dipoles with cutoff
* :doc:`lj/cut/dipole/long <pair_dipole>` - point dipoles with long-range Ewald
* :doc:`lj/cut/soft <pair_fep_soft>` - LJ with a soft core
* :doc:`lj/cut/thole/long <pair_thole>` - LJ with Coulombics with thole damping
* :doc:`lj/cut/thole/long <pair_thole>` - LJ with Coulomb with thole damping
* :doc:`lj/cut/tip4p/cut <pair_lj>` - LJ with cutoff Coulomb for TIP4P water
* :doc:`lj/cut/tip4p/long <pair_lj>` - LJ with long-range Coulomb for TIP4P water
* :doc:`lj/cut/tip4p/long/soft <pair_fep_soft>` - LJ with cutoff Coulomb for TIP4P water with a soft core
* :doc:`lj/expand <pair_lj_expand>` - Lennard-Jones for variable size particles
* :doc:`lj/expand/coul/long <pair_lj_expand>` - Lennard-Jones for variable size particles with long-range Coulombics
* :doc:`lj/expand/coul/long <pair_lj_expand>` - Lennard-Jones for variable size particles with long-range Coulomb
* :doc:`lj/gromacs <pair_gromacs>` - GROMACS-style Lennard-Jones potential
* :doc:`lj/gromacs/coul/gromacs <pair_gromacs>` - GROMACS-style LJ and Coulombic potential
* :doc:`lj/long/coul/long <pair_lj_long>` - long-range LJ and long-range Coulombics
* :doc:`lj/gromacs/coul/gromacs <pair_gromacs>` - GROMACS-style LJ and Coulomb potential
* :doc:`lj/long/coul/long <pair_lj_long>` - long-range LJ and long-range Coulomb
* :doc:`lj/long/dipole/long <pair_dipole>` - long-range LJ and long-range point dipoles
* :doc:`lj/long/tip4p/long <pair_lj_long>` - long-range LJ and long-range Coulombics for TIP4P water
* :doc:`lj/long/tip4p/long <pair_lj_long>` - long-range LJ and long-range Coulomb for TIP4P water
* :doc:`lj/mdf <pair_mdf>` - LJ potential with a taper function
* :doc:`lj/sdk <pair_sdk>` - LJ for SDK coarse-graining
* :doc:`lj/sdk/coul/long <pair_sdk>` - LJ for SDK coarse-graining with long-range Coulombics
* :doc:`lj/sdk/coul/msm <pair_sdk>` - LJ for SDK coarse-graining with long-range Coulombics via MSM
* :doc:`lj/sdk/coul/long <pair_sdk>` - LJ for SDK coarse-graining with long-range Coulomb
* :doc:`lj/sdk/coul/msm <pair_sdk>` - LJ for SDK coarse-graining with long-range Coulomb via MSM
* :doc:`lj/sf/dipole/sf <pair_dipole>` - LJ with dipole interaction with shifted forces
* :doc:`lj/smooth <pair_lj_smooth>` - smoothed Lennard-Jones potential
* :doc:`lj/smooth/linear <pair_lj_smooth_linear>` - linear smoothed LJ potential
@ -255,7 +256,7 @@ accelerated styles exist.
* :doc:`nb3b/harmonic <pair_nb3b_harmonic>` - non-bonded 3-body harmonic potential
* :doc:`nm/cut <pair_nm>` - N-M potential
* :doc:`nm/cut/coul/cut <pair_nm>` - N-M potential with cutoff Coulomb
* :doc:`nm/cut/coul/long <pair_nm>` - N-M potential with long-range Coulombics
* :doc:`nm/cut/coul/long <pair_nm>` - N-M potential with long-range Coulomb
* :doc:`oxdna/coaxstk <pair_oxdna>` -
* :doc:`oxdna/excv <pair_oxdna>` -
* :doc:`oxdna/hbond <pair_oxdna>` -
@ -315,7 +316,7 @@ accelerated styles exist.
* :doc:`tersoff/zbl <pair_tersoff_zbl>` - Tersoff/ZBL 3-body potential
* :doc:`thole <pair_thole>` - Coulomb interactions with thole damping
* :doc:`tip4p/cut <pair_coul>` - Coulomb for TIP4P water w/out LJ
* :doc:`tip4p/long <pair_coul>` - long-range Coulombics for TIP4P water w/out LJ
* :doc:`tip4p/long <pair_coul>` - long-range Coulomb for TIP4P water w/out LJ
* :doc:`tip4p/long/soft <pair_fep_soft>` -
* :doc:`tri/lj <pair_tri_lj>` - LJ potential between triangles
* :doc:`ufm <pair_ufm>` -

View File

@ -41,6 +41,8 @@ action pair_coul_long_cs.cpp pair_coul_long.cpp
action pair_coul_long_cs.h pair_coul_long.h
action pair_lj_cut_coul_long_cs.cpp pair_lj_cut_coul_long.cpp
action pair_lj_cut_coul_long_cs.h pair_lj_cut_coul_long.h
action pair_lj_class2_coul_long_cs.cpp pair_lj_class2_coul_long.cpp
action pair_lj_class2_coul_long_cs.h pair_lj_class2_coul_long.h
action pair_born_coul_wolf_cs.cpp pair_born_coul_wolf.cpp
action pair_born_coul_wolf_cs.h pair_born_coul_wolf.h

View File

@ -0,0 +1,561 @@
/* ----------------------------------------------------------------------
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 "pair_lj_class2_coul_long_cs.h"
#include <cmath>
#include "atom.h"
#include "force.h"
#include "neigh_list.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
#define EPSILON 1.0e-20
#define EPS_EWALD 1.0e-6
#define EPS_EWALD_SQR 1.0e-12
/* ---------------------------------------------------------------------- */
PairLJClass2CoulLongCS::PairLJClass2CoulLongCS(LAMMPS *lmp) : PairLJClass2CoulLong(lmp)
{
ewaldflag = pppmflag = 1;
respa_enable = 0; // TODO: r-RESPA handling is inconsistent and thus disabled until fixed
writedata = 1;
ftable = NULL;
}
/* ---------------------------------------------------------------------- */
void PairLJClass2CoulLongCS::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itable,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double fraction,table;
double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj;
double grij,expm2,prefactor,t,erfc,u;
double factor_coul,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
ev_init(eflag,vflag);
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]) {
rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond;
r2inv = 1.0/rsq;
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
prefactor = qqrd2e * qtmp*q[j];
if (factor_coul < 1.0) {
// When bonded parts are being calculated a minimal distance (EPS_EWALD)
// has to be added to the prefactor and erfc in order to make the
// used approximation functions for the Ewald correction valid
grij = g_ewald * (r+EPS_EWALD);
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
u = 1.0 - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
prefactor /= (r+EPS_EWALD);
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - (1.0-factor_coul));
// Additionally r2inv needs to be accordingly modified since the later
// scaling of the overall force shall be consistent
r2inv = 1.0/(rsq + EPS_EWALD_SQR);
} else {
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
u = 1.0 - 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);
}
} 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]) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - 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*erfc;
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]*r3inv-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 PairLJClass2CoulLongCS::compute_inner()
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
double rsw;
int *ilist,*jlist,*numneigh,**firstneigh;
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_inner;
ilist = list->ilist_inner;
numneigh = list->numneigh_inner;
firstneigh = list->firstneigh_inner;
double cut_out_on = cut_respa[0];
double cut_out_off = cut_respa[1];
double cut_out_diff = cut_out_off - cut_out_on;
double cut_out_on_sq = cut_out_on*cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off;
// 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_out_off_sq) {
rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond;
r2inv = 1.0/rsq;
forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
jtype = type[j];
if (rsq < cut_ljsq[itype][jtype]) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
} else forcelj = 0.0;
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
}
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;
}
}
}
}
}
/* ---------------------------------------------------------------------- */
void PairLJClass2CoulLongCS::compute_middle()
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
double rsw;
int *ilist,*jlist,*numneigh,**firstneigh;
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_middle;
ilist = list->ilist_middle;
numneigh = list->numneigh_middle;
firstneigh = list->firstneigh_middle;
double cut_in_off = cut_respa[0];
double cut_in_on = cut_respa[1];
double cut_out_on = cut_respa[2];
double cut_out_off = cut_respa[3];
double cut_in_diff = cut_in_on - cut_in_off;
double cut_out_diff = cut_out_off - cut_out_on;
double cut_in_off_sq = cut_in_off*cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on;
double cut_out_on_sq = cut_out_on*cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off;
// 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_out_off_sq && rsq > cut_in_off_sq) {
r2inv = 1.0/rsq;
forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
jtype = type[j];
if (rsq < cut_ljsq[itype][jtype]) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
} else forcelj = 0.0;
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
}
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
}
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;
}
}
}
}
}
/* ---------------------------------------------------------------------- */
void PairLJClass2CoulLongCS::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,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc,u;
double rsw;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
evdwl = ecoul = 0.0;
ev_init(eflag,vflag);
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;
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);
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 = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 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) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - 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*erfc;
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]) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
evdwl = r6inv*(lj3[itype][jtype]*r3inv-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 * (erfc + EWALD_F*grij*expm2);
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) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
} else if (rsq <= cut_in_on_sq) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
}
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,ecoul,fpair,delx,dely,delz);
}
}
}
}

View File

@ -0,0 +1,67 @@
/* -*- 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(lj/class2/coul/long/cs,PairLJClass2CoulLongCS)
#else
#ifndef LMP_PAIR_LJ_CLASS2_COUL_LONG_CS_H
#define LMP_PAIR_LJ_CLASS2_COUL_LONG_CS_H
#include "pair_lj_class2_coul_long.h"
namespace LAMMPS_NS {
class PairLJClass2CoulLongCS : public PairLJClass2CoulLong {
public:
PairLJClass2CoulLongCS(class LAMMPS *);
virtual void compute(int, int);
void compute_inner();
void compute_middle();
void compute_outer(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/class2/coul/long 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.
E: Pair cutoff < Respa interior cutoff
One or more pairwise cutoffs are too short to use with the specified
rRESPA cutoffs.
*/

View File

@ -41,7 +41,7 @@ using namespace LAMMPS_NS;
PairLJCutCoulLongCS::PairLJCutCoulLongCS(LAMMPS *lmp) : PairLJCutCoulLong(lmp)
{
ewaldflag = pppmflag = 1;
respa_enable = 1;
respa_enable = 0; // TODO: r-RESPA handling is inconsistent and thus disabled until fixed
writedata = 1;
ftable = NULL;
qdist = 0.0;
@ -173,7 +173,6 @@ void PairLJCutCoulLongCS::compute(int eflag, int vflag)
}
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];
@ -433,7 +432,6 @@ void PairLJCutCoulLongCS::compute_outer(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond;
r2inv = 1.0/rsq;
if (rsq < cut_coulsq) {
@ -442,9 +440,8 @@ void PairLJCutCoulLongCS::compute_outer(int eflag, int vflag)
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
u = 1. - t;
u = 1. - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
//erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
if (rsq > cut_in_off_sq) {
@ -543,7 +540,6 @@ void PairLJCutCoulLongCS::compute_outer(int eflag, int vflag)
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
}
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
}