diff --git a/doc/src/Commands_compute.txt b/doc/src/Commands_compute.txt index f566702609..e36234a306 100644 --- a/doc/src/Commands_compute.txt +++ b/doc/src/Commands_compute.txt @@ -91,6 +91,7 @@ KOKKOS, o = USER-OMP, t = OPT. "pe/atom"_compute_pe_atom.html, "pe/mol/tally"_compute_tally.html, "pe/tally"_compute_tally.html, +"pe/e3b"_compute_pe_e3b.html, "plasticity/atom"_compute_plasticity_atom.html, "pressure"_compute_pressure.html, "pressure/cylinder"_compute_pressure_cylinder.html, diff --git a/doc/src/Commands_pair.txt b/doc/src/Commands_pair.txt index e887f0178a..f56c991322 100644 --- a/doc/src/Commands_pair.txt +++ b/doc/src/Commands_pair.txt @@ -80,6 +80,7 @@ OPT. "dpd/fdt/energy (k)"_pair_dpd_fdt.html, "dpd/tstat (go)"_pair_dpd.html, "dsmc"_pair_dsmc.html, +"e3b"_pair_e3b.html, "eam (gikot)"_pair_eam.html, "eam/alloy (gikot)"_pair_eam.html, "eam/cd (o)"_pair_eam.html, diff --git a/doc/src/Eqs/e3b.tex b/doc/src/Eqs/e3b.tex new file mode 100644 index 0000000000..550538bf35 --- /dev/null +++ b/doc/src/Eqs/e3b.tex @@ -0,0 +1,15 @@ +\documentclass[12pt]{article} +\usepackage{amsmath} +\begin{document} + +\begin{align*} +E =& E_2 \sum_{i,j}e^{-k_2 r_{ij}} + E_A \sum_{\substack{i,j,k,\ell \\\in \textrm{type A}}} f(r_{ij})f(r_{k\ell}) + E_B \sum_{\substack{i,j,k,\ell \\\in \textrm{type B}}} f(r_{ij})f(r_{k\ell}) + E_C \sum_{\substack{i,j,k,\ell \\\in \textrm{type C}}} f(r_{ij})f(r_{k\ell}) \\ +f(r) =& e^{-k_3 r}s(r) \\ +s(r) =& \begin{cases} + 1 & rR_f\\ +\end{cases} +\end{align*} + +\end{document} diff --git a/doc/src/compute.txt b/doc/src/compute.txt index 87dbee57d6..047d838b2d 100644 --- a/doc/src/compute.txt +++ b/doc/src/compute.txt @@ -240,6 +240,7 @@ compute"_Commands_compute.html doc page are followed by one or more of "pe/atom"_compute_pe_atom.html - potential energy for each atom "pe/mol/tally"_compute_tally.html - "pe/tally"_compute_tally.html - +"pe/e3b"_compute_pe_e3b.html - potential energy from pair_style e3b "plasticity/atom"_compute_plasticity_atom.html - Peridynamic plasticity for each atom "pressure"_compute_pressure.html - total pressure and pressure tensor "pressure/cylinder"_compute_pressure_cylinder.html - diff --git a/doc/src/compute_pe_e3b.txt b/doc/src/compute_pe_e3b.txt new file mode 100644 index 0000000000..a47d8beaa9 --- /dev/null +++ b/doc/src/compute_pe_e3b.txt @@ -0,0 +1,60 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +compute pe/e3b command :h3 + +[Syntax:] + +compute ID group-ID pe/e3b :pre + +ID, group-ID are documented in "compute"_compute.html command +pe/e3b = style name of this compute command :ul + +[Examples:] + +compute 1 all pe/e3b :pre + +[Description:] + +Define a computation that calculates the contribution of "pair_style e3b"_pair_e3b.html to the potential energy. +The specified group must be "all". +See the "compute pe/atom"_compute_pe_atom.html command if you want per-atom +energies. +These per-atom values could be summed for a group of atoms via the "compute reduce"_compute_reduce.html command. + +The "pair_style e3b"_pair_e3b.html potential is composed of 4 terms. +This compute calculates the total e3b contribution to the energy as well as each of the four terms. +The four terms are stored as a 4-element vector in the order pe_Ea, pe_Eb, pe_Ec, pe_E2. +See "pair_style e3b"_pair_e3b.html for more details. + +:line + +[Output info:] + +This compute calculates a global scalar (the total e3b energy) and a global +vector of length 4 (the four energy terms), which can be accessed by indices +1-4. These values can be used by any command that uses global scalar +or vector values from a compute as input. See the "Howto +output"_Howto_output.html doc page for an overview of LAMMPS output +options. + +The scalar and vector values calculated by this compute are +"extensive" and in energy +"units"_units.html. + +[Restrictions:] + +This compute must be used with "pair_style e3b"_pair_e3b.html. + +[Related commands:] + +"pair_style e3b"_pair_e3b.html, +"compute pe"_compute_pe.html, +"compute pe/atom"_compute_pe_atom.html + +[Default:] none diff --git a/doc/src/computes.txt b/doc/src/computes.txt index 926b8da222..eefbe5116c 100644 --- a/doc/src/computes.txt +++ b/doc/src/computes.txt @@ -66,6 +66,7 @@ Computes :h1 compute_pair_local compute_pe compute_pe_atom + compute_pe_e3b compute_plasticity_atom compute_pressure compute_pressure_cylinder diff --git a/doc/src/pair_e3b.txt b/doc/src/pair_e3b.txt new file mode 100644 index 0000000000..1121cd4034 --- /dev/null +++ b/doc/src/pair_e3b.txt @@ -0,0 +1,138 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +pair_style e3b command :h3 + +[Syntax:] + +pair_style e3b Otype :pre +Otype = atom type for oxygen :l + +pair_coeff * * keyword :pre +one or more keyword/value pairs must be appended. :l +keyword = {preset} or {Ea} or {Eb} or {Ec} or {E2} or {K3} or {K2} or {Rs} or {Rc3} or {Rc2} or {bondL} or {neigh} :l +If the {preset} keyword is given, no others are needed. +Otherwise, all are mandatory except for {neigh}. +The {neigh} keyword is always optional. :l + {preset} arg = {2011} or {2015} = which set of predefined parameters to use + 2011 = use the potential parameters from "(Tainter 2011)"_#Tainter2011 + 2015 = use the potential parameters from "(Tainter 2015)"_#Tainter2015 + {Ea} arg = three-body energy for type A hydrogen bonding interactions (energy units) + {Eb} arg = three-body energy for type B hydrogen bonding interactions (energy units) + {Ec} arg = three-body energy for type C hydrogen bonding interactions (energy units) + {E2} arg = two-body energy correction (energy units) + {K3} arg = three-body exponential constant (inverse distance units) + {K2} arg = two-body exponential constant (inverse distance units) + {Rc3} arg = three-body cutoff (distance units) + {Rc2} arg = two-body cutoff (distance units) + {Rs} arg = three-body switching function cutoff (distance units) + {bondL} arg = intramolecular OH bond length (distance units) + {neigh} arg = approximate integer number of molecules within Rc3 of an oxygen atom :pre + +[Examples:] + +pair_style e3b 1 +pair_coeff * * Ea 35.85 Eb -240.2 Ec 449.3 E2 108269.9 K3 1.907 K2 4.872 Rc3 5.2 Rc2 5.2 Rs 5.0 bondL 0.9572 :pre + +pair_style hybrid/overlay e3b 1 lj/cut/tip4p/long 1 2 1 1 0.15 8.5 +pair_coeff * * e3b preset 2011 :pre + +[Description:] + +The {e3b} style computes an \"explicit three-body\" (E3B) potential for water "(Kumar 2008)"_#Kumar. + +:c,image(Eqs/e3b.jpg) + +This potential was developed as a water model that includes the three-body cooperativity of hydrogen bonding explicitly. +To use it in this way, it must be applied in conjunction with a conventional two-body water model, through {pair_style hybrid/overlay}. +The three body interactions are split into three types: A, B, and C. +Type A corresponds to anti-cooperative double hydrogen bond donor interactions. +Type B corresponds to the cooperative interaction of molecules that both donate and accept a hydrogen bond. +Type C corresponds to anti-cooperative double hydrogen bond acceptor interactions. +The three-body interactions are smoothly cutoff by the switching function s(r) between Rs and Rc3. +The two-body interactions are designed to correct for the effective many-body interactions implicitly included in the conventional two-body potential. +The two-body interactions are cut off sharply at Rc2, because K3 is typically significantly smaller than K2. +See "(Kumar 2008)"_#Kumar for more details. + +Only a single {pair_coeff} command is used with the {e3b} style. +The 1st two arguments must be * *. +The oxygen atom type for the pair style is passed as the only argument to the {pair_style} command, not in the {pair_coeff} command. +The hydrogen atom type is inferred by the ordering of the atoms. + +NOTE: Every atom of type Otype must be part of a water molecule. +Each water molecule must have consecutive IDs with the oxygen first. +This pair style does not test that this criteria is met. + +The {pair_coeff} command must have at least one keyword/value pair, as described above. +The {preset} keyword sets the potential parameters to the values used in "(Tainter 2011)"_#Tainter2011 or "(Tainter 2015)"_#Tainter2015. +To use the water models defined in those references, the {e3b} style should always be used in conjunction with an {lj/cut/tip4p/long} style through {pair_style hybrid/overlay}, as demonstrated in the second example above. +The {preset 2011} option should be used with the "TIP4P water model"_Howto_tip4p.html. +The {preset 2015} option should be used with the "TIP4P/2005 water model"_Howto_tip4p.html. +If the {preset} keyword is used, no other keyword is needed. +Changes to the preset parameters can be made by specifying the {preset} keyword followed by the specific parameter to change, like {Ea}. +Note that the other keywords must come after {preset} in the pair_style command. +The {e3b} style can also be used to implement any three-body potential of the same form by specifying all the keywords except {neigh}: {Ea}, {Eb}, {Ec}, {E2}, {K3}, {K2}, {Rc3}, {Rc2}, {Rs}, and {bondL}. +The keyword {bondL} specifies the intramolecular OH bond length of the water model being used. +This is needed to include H atoms that are within the cutoff even when the attached oxygen atom is not. + +This pair style allocates arrays sized according to the number of pairwise interactions within Rc3. +To do this it needs an estimate for the number of water molecules within Rc3 of an oxygen atom. +This estimate defaults to 10 and can be changed using the {neigh} keyword, which takes an integer as an argument. +If the neigh setting is too small, the simulation will fail with the error "neigh is too small". +If the neigh setting is too large, the pair style will use more memory than necessary. + +This pair style makes 4 different contributions to the potential energy from the E2, Ea, Eb, and Ec terms above. +The value of each of these terms can be computed using "compute pe/e3b"_compute_pe_e3b.html. + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +This pair style does not support the "pair_modify"_pair_modify.html +shift, table, and tail options. + +This pair style does not write its information to "binary restart +files"_restart.html, since it is stored in potential files. Thus, you +need to re-specify the pair_style and pair_coeff commands in an input +script that reads a restart file. + +This pair style is incompatible with "respa"_run_style.html. + +:line + +[Restrictions:] + +This pair style is part of the USER-MISC package. It is only enabled +if LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +This pair style requires the "newton"_newton.html setting to be "on" +for pair interactions. + +This pair style requires a fixed number of atoms in the simulation, so it is incompatible with fixes like "fix deposit"_fix_deposit.html. +If the number of atoms changes between runs, this pair style must be re-initialized by calling the {pair_style} and {pair_coeffs} commands. +This is not a fundamental limitation of the pair style, but the code currently does not support a variable number of atoms. + +The {preset} keyword currently only works with real, metal, si, and cgs "units"_units.html. + +[Related commands:] + +"pair_coeff"_pair_coeff.html, "compute pe/e3b"_compute_pe_e3b.html + +[Default:] + +The option default for the {neigh} keyword is 10. + +:line + +:link(Kumar) +[(Kumar)] Kumar and Skinner, J. Phys. Chem. B, 112, 8311 (2008) +:link(Tainter2011) +[(Tainter 2011)] Tainter, Pieniazek, Lin, and Skinner, J. Chem. Phys., 134, 184501 (2011) +:link(Tainter2015) +[(Tainter 2015)] Tainter, Shi, and Skinner, 11, 2268 (2015) diff --git a/doc/src/pair_style.txt b/doc/src/pair_style.txt index f6fcd110d8..27cc168236 100644 --- a/doc/src/pair_style.txt +++ b/doc/src/pair_style.txt @@ -147,6 +147,7 @@ accelerated styles exist. "dpd/fdt/energy"_pair_dpd_fdt.html - DPD for constant energy and enthalpy "dpd/tstat"_pair_dpd.html - pair-wise DPD thermostatting "dsmc"_pair_dsmc.html - Direct Simulation Monte Carlo (DSMC) +"e3b"_pair_e3b.html - Explicit-three body (E3B) water model "eam"_pair_eam.html - embedded atom method (EAM) "eam/alloy"_pair_eam.html - alloy EAM "eam/cd"_pair_eam.html - concentration-dependent EAM diff --git a/doc/src/pairs.txt b/doc/src/pairs.txt index 30dcc8fd4b..119771cd0d 100644 --- a/doc/src/pairs.txt +++ b/doc/src/pairs.txt @@ -32,6 +32,7 @@ Pair Styles :h1 pair_dpd pair_dpd_fdt pair_dsmc + pair_e3b pair_eam pair_edip pair_eff diff --git a/examples/USER/e3b/README b/examples/USER/e3b/README new file mode 100644 index 0000000000..7f13c60386 --- /dev/null +++ b/examples/USER/e3b/README @@ -0,0 +1,10 @@ +The input script in.lammps simulates bulk water using the 2015 E3B potential. +It can be modified to use the 2011 E3B potential by + 1) using a tip4p molecule file instead of tip4p2005.mol + 2) changing the qdist parameter for lj/cut/tip4p/long in the pair_style + hybrid/overlay command + 3) using the 2011 "preset" command in the e3b pair_coeff command +This script also demonstrates the use of compute pe/e3b to calculate the +potential energy contribution of the e3b pair style. These potential energy +contributions can be found in the output file e3b.txt. See the LAMMPS +documentation for more details. \ No newline at end of file diff --git a/examples/USER/e3b/in.lammps b/examples/USER/e3b/in.lammps new file mode 100644 index 0000000000..3960c9ad5e --- /dev/null +++ b/examples/USER/e3b/in.lammps @@ -0,0 +1,115 @@ +#LAMMPS input file +#to simulate bulk E3B3 water model + +##################################################################### +clear + +variable samp_rate equal 10 +variable thermo_rate equal 10 +variable Wlat equal 3.10744 #for water density 0.997g/mL +variable L equal 3 #(L*2)^3 = Nmolec, L=3 -> N=216 + +variable equil equal 100 +variable run equal 100 + +variable ts equal 2.0 +variable Tdamp equal 100*${ts} +variable Pdamp equal 1000*${ts} +variable myT equal 298.0 +variable myP equal 1.0 + +units real +atom_style full + +dimension 3 + +boundary p p p + +############################################################################# +#setup box + +lattice sc ${Wlat} +region simbox block -$L $L -$L $L -$L $L units lattice + +############################################################################# +#set up potential + +create_box 2 simbox bond/types 1 angle/types 1 extra/bond/per/atom 2 & + extra/special/per/atom 2 extra/angle/per/atom 1 + +molecule h2o tip4p2005.mol +mass 1 15.9994 #oxygen +mass 2 1.008 #hydrogen + +pair_style hybrid/overlay e3b 1 lj/cut/tip4p/long 1 2 1 1 0.1546 8.5 +pair_modify table 0 table/disp 0 shift yes + +bond_style harmonic +angle_style harmonic + +kspace_style pppm/tip4p 1.0e-6 + +pair_coeff * * lj/cut/tip4p/long 0.0 0.0 +pair_coeff 1 1 lj/cut/tip4p/long 0.1852 3.1589 +pair_coeff * * e3b preset 2015 + +#intramolecular bond/angle coeffs very stiff for minimization +bond_coeff 1 100000 0.9572 +angle_coeff 1 100000 104.52 + +############################################################################# +#setup for run +thermo ${thermo_rate} +thermo_style custom step vol temp epair pe etotal press density + +timestep ${ts} +run_style verlet + +neighbor 2.0 bin +neigh_modify every 1 delay 3 check yes + +############################################################################# +#make atoms and rough minimze +create_atoms 0 box mol h2o 15856 + +#dump positions only in first batch run +dump 7 all custom ${samp_rate} dump.lammpstrj id x y z +dump_modify 7 sort id + +min_style cg +minimize 1.0e-4 1.0e-4 10000 100000 + +#potential coeffs aren't too important since will be rigid anyways +bond_coeff 1 554.13 0.9572 +angle_coeff 1 45.769 104.52 + + +############################################################################# +#initialize velocity and rigid constraint + +fix rigid all shake 1.0e-8 100 0 b 1 a 1 t 1 2 mol h2o +velocity all create ${myT} 15856 dist gaussian rot yes mom yes + +#scale velocity +run 0 +velocity all scale ${myT} + +compute e3b all pe/e3b +fix e3b all ave/time 1 1 ${thermo_rate} c_e3b c_e3b[*] & + file e3b.txt title2 "step pe_e3b pe_ea pe_eb pe_ec pe_e2" + +############################################################################# +#equilibrate bulk water at NVT + +fix 1 all nvt temp ${myT} ${myT} ${Tdamp} +run ${equil} + +############################################################################# +#run at NVT + +dump 1 all custom ${samp_rate} dump.lammpstrj id x y z type +dump_modify 1 sort id + +run ${run} + +write_restart lammps.restart diff --git a/examples/USER/e3b/tip4p2005.mol b/examples/USER/e3b/tip4p2005.mol new file mode 100644 index 0000000000..8dff52bebe --- /dev/null +++ b/examples/USER/e3b/tip4p2005.mol @@ -0,0 +1,61 @@ +#TIP4P/2005 H20 water molecule, parameters from section 6.9 of LAMMPS manual +3 atoms +2 bonds +1 angles + +Coords + +1 0.0 0.0 0.0 +2 0.58588228 0.75695033 0.0 +3 0.58588228 -0.75695033 0.0 + +Types + +1 1 +2 2 +3 2 + +Charges + +1 -1.1128 +2 0.5564 +3 0.5564 + +Bonds + +1 1 1 2 +2 1 1 3 + +Angles + +1 1 2 1 3 + +Special Bond Counts + +1 2 0 0 +2 1 1 0 +3 1 1 0 + +Special Bonds + +1 2 3 +2 1 3 +3 1 2 + +Shake Flags + +1 1 +2 1 +3 1 + +Shake Atoms + +1 1 2 3 +2 1 2 3 +3 1 2 3 + +Shake Bond Types + +1 1 1 1 +2 1 1 1 +3 1 1 1 diff --git a/src/USER-MISC/compute_pe_e3b.cpp b/src/USER-MISC/compute_pe_e3b.cpp new file mode 100644 index 0000000000..9301883ac2 --- /dev/null +++ b/src/USER-MISC/compute_pe_e3b.cpp @@ -0,0 +1,86 @@ +/* ---------------------------------------------------------------------- + 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. + +Adapted from USER-TALLY/compute_pe_tally.cpp by Steven E Strong +Also used code from compute_pe.cpp +------------------------------------------------------------------------- */ + +#include "compute_pe_e3b.h" +#include "pair_e3b.h" + +#include "pair.h" +#include "update.h" +#include "error.h" +#include "force.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputePEE3B::ComputePEE3B(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), e3b(NULL) +{ + // 0 1 2 + //compute ID grp pe/e3b + if (narg != 3) error->all(FLERR,"Illegal compute pe/e3b command"); + + scalar_flag = 1; + vector_flag = 1; + size_vector = 4; //etotA,etotB,etotC,etot2 + extvector = extscalar = 1; + timeflag = 1; + + peflag = 1; // we need Pair::ev_tally() to be run + + invoked_vector = invoked_scalar = -1; + vector = new double[size_vector]; +} + +/* ---------------------------------------------------------------------- */ + +ComputePEE3B::~ComputePEE3B() +{ + delete[] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputePEE3B::init() { + Pair *pair = force->pair_match("e3b",false,0); + if (pair==NULL) + error->all(FLERR,"This compute must be used with pair_style e3b"); + + e3b = (PairE3B *) pair; + if (e3b==NULL) + error->all(FLERR,"something went wrong"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputePEE3B::compute_vector() +{ + invoked_vector = update->ntimestep; + if (update->eflag_global != invoked_scalar) + error->all(FLERR,"Energy was not tallied on needed timestep"); + + // sum energies across procs + MPI_Allreduce(e3b->etot,vector,4,MPI_DOUBLE,MPI_SUM,world); +} + +double ComputePEE3B::compute_scalar() { + invoked_scalar = update->ntimestep; + if (invoked_scalar != invoked_vector) + compute_vector(); + + scalar = vector[0]+vector[1]+vector[2]+vector[3]; + return scalar; +} diff --git a/src/USER-MISC/compute_pe_e3b.h b/src/USER-MISC/compute_pe_e3b.h new file mode 100644 index 0000000000..a23c48b38d --- /dev/null +++ b/src/USER-MISC/compute_pe_e3b.h @@ -0,0 +1,45 @@ +/* -*- 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 COMPUTE_CLASS + +ComputeStyle(pe/e3b,ComputePEE3B) + +#else + +#ifndef LMP_COMPUTE_PEE3B_H +#define LMP_COMPUTE_PEE3B_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputePEE3B : public Compute { + + public: + ComputePEE3B(class LAMMPS *, int, char **); + virtual ~ComputePEE3B(); + + void init(); + + double compute_scalar(); + void compute_vector(); + + private: + class PairE3B *e3b; +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/pair_e3b.cpp b/src/USER-MISC/pair_e3b.cpp new file mode 100644 index 0000000000..84a3400787 --- /dev/null +++ b/src/USER-MISC/pair_e3b.cpp @@ -0,0 +1,690 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Written by Steven E Strong and Nicholas J Hestand + Adapted from MANYBODY/pair_sw.cpp and an + implementation of E3B in GROMACS by Craig Tainter (?) +------------------------------------------------------------------------- */ + +#include "pair_e3b.h" + +#include "atom.h" +#include "neighbor.h" +#include "neigh_request.h" +#include "neigh_list.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "error.h" +#include "update.h" +#include "domain.h" +#include "citeme.h" + +#include +#include +#include +#include + +//these are defined here to avoid confusing hardcoded indicies, but +//they do not allow flexibility. If they are changed the code will break +#define DIM 3 +#define NUMH 2 //number of hydrogen atoms per water molecule +#define NUMO 2 //number of oxygen atoms per pair of water molecules +#define BOND_DELTA 1.01 //buffer for OH bonds that aren't perfectly rigid + +using namespace LAMMPS_NS; +/* ---------------------------------------------------------------------- */ + +PairE3B::PairE3B(LAMMPS *lmp) : Pair(lmp),pairPerAtom(10) +{ + single_enable = 0; + restartinfo = 0; + one_coeff = 1; + manybody_flag = 1; + + allocatedE3B = false; + pairO = NULL; + pairH = NULL; + exps = NULL; + del3 = NULL; + fpair3 = NULL; + sumExp = NULL; +} + +/* ---------------------------------------------------------------------- + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairE3B::~PairE3B() +{ + if (copymode) return; + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + } + if (allocatedE3B) { + memory->destroy(pairO); + memory->destroy(pairH); + memory->destroy(exps); + memory->destroy(del3); + memory->destroy(fpair3); + memory->destroy(sumExp); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairE3B::compute(int eflag, int vflag) +{ + int i,j,k,h,ii,jj,hh,kk,inum,jnum,otherO; + tagint itag,jtag; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,rsq,tmpexp; + double fxtmp,fytmp,fztmp,fix,fiy,fiz; + double delxh,delyh,delzh,rsqh,tmpr; + double scFact1,scFact2,scEng,scDer; + int *ilist,*jlist,*numneigh,**firstneigh; + bool addedH; + + if (natoms != atom->natoms) + error->all(FLERR,"pair E3B requires a fixed number of atoms"); + + //clear sumExp array + memset(sumExp,0.0,nbytes); + + evdwl = 0.0; + etot[0]=etot[1]=etot[2]=etot[3]=0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + tagint *tag = atom->tag; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + int npair = 0; + // loop over half neighbor list of my atoms + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (type[i]!=typeO) + continue; + + itag = tag[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + fix = fiy = fiz = 0.0; + + // two-body interactions + jlist = firstneigh[i]; + jnum = numneigh[i]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + //skip unless O-O interaction + if (type[j]!=typeO) + continue; + + jtag = tag[j]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; //OO distance + + //two body interaction + //not shifted b/c k2=4.87/A, so at cutoff (5.2A) e^(-kr) = 1e-11 + if (rsq < rc2sq) { + tmpr = sqrt(rsq); + tmpexp = e2 * exp(-k2*tmpr); + fpair = k2 * tmpexp / tmpr; + + fxtmp = delx*fpair; + fytmp = dely*fpair; + fztmp = delz*fpair; + fix += fxtmp; + fiy += fytmp; + fiz += fztmp; + f[j][0] -= fxtmp; + f[j][1] -= fytmp; + f[j][2] -= fztmp; + + if (evflag) { + ev_tally(i,j,nlocal,newton_pair,tmpexp,0.0,fpair,delx,dely,delz); + etot[3] += tmpexp; + } + } //end if rsqmap(tag[otherO]+hh+1); + //if hydrogen atom is missing, bond potential or shake will + //catch this, so don't need to check here + //if (h<0) + // error->one(FLERR,"hydrogen atom missing"); + h = domain->closest_image(otherO,h); + pairH[npair][kk][hh] = h; + + delxh = x[k][0] - x[h][0]; + delyh = x[k][1] - x[h][1]; + delzh = x[k][2] - x[h][2]; + rsqh = delxh*delxh + delyh*delyh + delzh*delzh; + + if (rsqh < rc3sq) { + + tmpr = sqrt(rsqh); + tmpexp = exp(-k3*tmpr); + if (tmpr > rs) { + scFact1 = rc3-tmpr; + scFact2 = sc_num + 2*tmpr; + scEng = scFact1*scFact1*scFact2*sc_denom; + scDer = k3*scEng - 6*scFact1*(rs-tmpr)*sc_denom; + } else { + scDer = k3; + scEng = 1.0; + } + + //need to keep fpair3 separate from del3 for virial + fpair3[npair][kk][hh] = scDer*tmpexp/tmpr; + tmpexp *= scEng; + exps[npair][kk][hh] = tmpexp; + del3[npair][kk][hh][0] = delxh; + del3[npair][kk][hh][1] = delyh; + del3[npair][kk][hh][2] = delzh; + + //accumulate global vector of sum(e^kr) + //tags start at 1, so subtract one to index sumExp + sumExp[tag[k]-1] += tmpexp; + sumExp[tag[h]-1] += tmpexp; + + addedH = true; + } else { + exps [npair][kk][hh] = 0.0; + fpair3[npair][kk][hh] = 0.0; + } //if < rc3sq + } //end loop through 2 Hs + } //end for kk in NUMO + //if added a pair, check if array is too small and grow + if (addedH) { + npair++; + if (npair >= pairmax) + error->one(FLERR,"neigh is too small"); + } + } //end if < rc3deltaSq + } //end for jj neigh + + //add 2-body forces on i + f[i][0] += fix; + f[i][1] += fiy; + f[i][2] += fiz; + } //end for ii + + //communicate sumExp array + //tested that no change in speed with MPI_IN_PLACE + MPI_Allreduce(MPI_IN_PLACE,sumExp,maxID,MPI_DOUBLE,MPI_SUM,world); + + //now loop through list of pairs, calculating 3body forces + int j2,otherH; + double partA,partB,partC; + for (ii = 0; ii < npair; ii++) { + + for (kk=0; kkntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + memory->create(cutsq,n+1,n+1,"pair:cutsq"); +} + +void PairE3B::allocateE3B() +{ + allocatedE3B = true; + + //TODO: get memory->grow working for 4d arrays + pairmax = atom->nlocal*pairPerAtom; //initial guess for size of pair lists + memory->create(pairO ,pairmax,NUMO ,"pair:pairO"); + memory->create(pairH ,pairmax,NUMO,NUMH ,"pair:pairH"); + memory->create(exps ,pairmax,NUMO,NUMH ,"pair:exps"); + memory->create(fpair3,pairmax,NUMO,NUMH ,"pair:fpair3"); + memory->create(del3 ,pairmax,NUMO,NUMH,DIM,"pair:del3"); + + natoms=atom->natoms; + maxID=find_maxID(); + if (!natoms) + error->all(FLERR,"No atoms found"); + memory->create(sumExp,maxID,"pair:sumExp"); + nbytes = sizeof(double) * maxID; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairE3B::settings(int narg, char **arg) +{ + if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + typeO=force->inumeric(FLERR,arg[0]); + if (typeO<1 || typeO>atom->ntypes) + error->all(FLERR,"Invalid Otype: out of bounds"); +} + +/* ---------------------------------------------------------------------- + coeffs must be * * keyword/value + keyword/values set the potential parameters +------------------------------------------------------------------------- */ +void PairE3B::coeff(int narg, char **arg) +{ + if (!allocated) allocate(); + + //1=* 2=* 3/4=1st keyword/value + if (narg < 4) + error->all(FLERR,"There must be at least one keyword given to pair_coeff"); + + // ensure I,J args are * * + if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) + error->all(FLERR,"Incorrect args for pair coefficients"); + + // clear setflag since coeff() called once with I,J = * * + int n = atom->ntypes; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + setflag[typeO][typeO]=1; + + //parse keyword/value pairs + double bondL=0.0; //OH bond length + bool repeatFlag=false; + int presetFlag; + + //clear parameters + e2=ea=eb=ec=k3=k2=NAN; + rs=rc3=rc2=0.0; + + int iarg = 2; //beginning of keyword/value pairs + while(iarg < narg) { + char *keyword = arg[iarg++]; + if (checkKeyword(keyword,"Ea",1,narg-iarg)) + ea=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"Eb",1,narg-iarg)) + eb=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"Ec",1,narg-iarg)) + ec=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"K3",1,narg-iarg)) + k3=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"Rs",1,narg-iarg)) + rs=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"Rc3",1,narg-iarg)) + rc3=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"Rc2",1,narg-iarg)) + rc2=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"bondL",1,narg-iarg)) + bondL=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"E2",1,narg-iarg)) + e2=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"K2",1,narg-iarg)) + k2=force->numeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"neigh",1,narg-iarg)) + pairPerAtom=force->inumeric(FLERR,arg[iarg++]); + else if (checkKeyword(keyword,"preset",1,narg-iarg)) { + presetFlag=force->inumeric(FLERR,arg[iarg++]); + presetParam(presetFlag,repeatFlag,bondL); + } else { + char str[256]; + snprintf(str,256,"Keyword %s is unknown",keyword); + error->all(FLERR,str); + } + } + + checkInputs(bondL); + + //cutmax for neighbor listing + cutmax = std::max(rc2,rc3); + rc2sq = rc2*rc2; + rc3sq = rc3*rc3; + rc3deltaSq = (rc3+bondL)*(rc3+bondL); + + double tmpfact=1.0/(rc3-rs); + sc_denom=tmpfact*tmpfact*tmpfact; + sc_num=rc3-3*rs; +} + +/* ---------------------------------------------------------------------- + init specific to this pair styles +------------------------------------------------------------------------- */ + +void PairE3B::init_style() +{ + if (atom->tag_enable == 0) + error->all(FLERR,"Pair style E3B requires atom IDs"); + if (force->newton_pair == 0) + error->all(FLERR,"Pair style E3B requires newton pair on"); + + // need a half neighbor list + int irequest = neighbor->request(this,instance_me); + //don't need this, half is default + //neighbor->requests[irequest]->half = 0; + + if (!force->pair_match("tip4p",false,0)) + if (comm->me==0) error->warning(FLERR,"E3B pair_style is designed for use with hybrid/overlay tip4p style"); + + if (!allocatedE3B) allocateE3B(); +} + + +static const char cite_E3B1[] = + "Explicit Three-Body (E3B) potential for water:\n\n" + "@article{kumar_water_2008,\n" + "title = {Water Simulation Model with Explicit Three-Molecule Interactions},\n" + "volume = {112},\n" + "doi = {10.1021/jp8009468},\n" + "number = {28},\n" + "journal = {J Phys. Chem. B},\n" + "author = {Kumar, R. and Skinner, J. L.},\n" + "year = {2008},\n" + "pages = {8311--8318}\n" + "}\n\n"; + +static const char cite_E3B2[] = + "Explicit Three-Body (E3B) potential for water:\n\n" + "@article{tainter_robust_2011,\n" + "title = {Robust three-body water simulation model},\n" + "volume = {134},\n" + "doi = {10.1063/1.3587053},\n" + "number = {18},\n" + "journal = {J. Chem. Phys},\n" + "author = {Tainter, C. J. and Pieniazek, P. A. and Lin, Y.-S. and Skinner, J. L.},\n" + "year = {2011},\n" + "pages = {184501}\n" + "}\n\n"; + +static const char cite_E3B3[] = + "Explicit Three-Body (E3B) potential for water:\n\n" + "@article{tainter_reparametrized_2015,\n" + "title = {Reparametrized {E3B} (Explicit Three-Body) Water Model Using the {TIP4P/2005} Model as a Reference},\n" + "volume = {11},\n" + "doi = {10.1021/acs.jctc.5b00117},\n" + "number = {5},\n" + "journal = {J. Chem. Theory Comput.},\n" + "author = {Tainter, Craig J. and Shi, Liang and Skinner, James L.},\n" + "year = {2015},\n" + "pages = {2268--2277}\n" + "}\n\n"; + +void PairE3B::presetParam(const int flag,bool &repeatFlag,double &bondL) { + if (repeatFlag) { + error->all(FLERR, + "Cannot request two different sets of preset parameters"); + } + repeatFlag=true; + + if (!isnan(ea) || !isnan(eb) || !isnan(ec) || !isnan(e2) || bondL!=0.0 || + !isnan(k3) || !isnan(k2) || rs!=0.0 || rc3!=0.0 || rc2!=0.0 ) + error->all(FLERR,"Preset keyword will overwrite another keyword setting"); + + double econv,lconv; + if (strcmp(update->unit_style,"real") == 0) { + econv=1.0/4.184; + lconv=1.0; + } else if (strcmp(update->unit_style,"metal") == 0) { + econv=0.103653271; + lconv=1.0; + } else if (strcmp(update->unit_style,"si") == 0) { + econv=1.660578e-21; + lconv=1e-10; + } else if (strcmp(update->unit_style,"cgs") == 0) { + econv=1.660578e-14; + lconv=1e-8; + } else { + char str[256]; + snprintf(str,256, + "Pre-defined E3B parameters have not been set for %s units.", + update->unit_style); + error->all(FLERR,str); + } + + //here parameters are defined in kJ/mol and A + //they will be converted to the lammps units after + if (flag==2008) { + error->all(FLERR,"\"preset 2008\" is not yet supported, because this would require distinct k3 coefficients, use \"preset 2011\" or \"preset 2015\""); + if (lmp->citeme) lmp->citeme->add(cite_E3B1); + ea = 4699.6; + eb =-2152.9; + ec = 1312.7; + //ka = 1.0/1.88; + //kb = 1.0/1.71; + //kc = 1.0/1.56; + e2 = 1.925e6; + k2 = 4.67; + rs = 5.0; + rc3 = 5.2; + rc2 = 5.2; + bondL = 0.9572; + } else if (flag==2011) { + if (lmp->citeme) lmp->citeme->add(cite_E3B2); + ea = 1745.7; + eb =-4565.0; + ec = 7606.8; + k3 = 1.907; + e2 = 2.349e6; + k2 = 4.872; + rs = 5.0; + rc3 = 5.2; + rc2 = 5.2; + bondL = 0.9572; + } else if (flag==2015) { + if (lmp->citeme) lmp->citeme->add(cite_E3B3); + ea = 150.0; + eb =-1005.0; + ec = 1880.0; + k3 = 1.907; + e2 = 0.453e6; + k2 = 4.872; + rs = 5.0; + rc3 = 5.2; + rc2 = 5.2; + bondL = 0.9572; + } else + error->all(FLERR,"Unknown argument: preset only takes 2011 or 2015 as arguments"); + + //convert units + ea *= econv; + eb *= econv; + ec *= econv; + e2 *= econv; + k3 /= lconv; + k2 /= lconv; + rs *= lconv; + rc2 *= lconv; + rc3 *= lconv; + bondL *= lconv*BOND_DELTA; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ +//pair.cpp::init uses this to set cutsq array, used for neighboring, etc +double PairE3B::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + return cutmax; +} + +bool PairE3B::checkKeyword(const char *thiskey,const char *test,const int nVal, const int nRem) { + if (strcmp(thiskey,test) == 0) { + if(nRemall(FLERR,str); + } + return true; + } + return false; +} + +//find max atom ID for all atoms +//from fix_deposit.cpp +tagint PairE3B::find_maxID() +{ + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + + tagint max = 0; + tagint maxID; + for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); + MPI_Allreduce(&max,&maxID,1,MPI_LMP_TAGINT,MPI_MAX,world); + + return maxID; +} + +void PairE3B::checkInputs(const double &bondL) { + //first check that all necessary values were set + if (rc2==0.0) + error->all(FLERR,"rc2 keyword missing"); + if (rs==0.0) + error->all(FLERR,"Rs keyword missing"); + if (rc3==0.0) + error->all(FLERR,"Rc3 keyword missing"); + if (bondL==0.0) + error->all(FLERR,"bondL keyword missing"); + if (isnan(ea)) + error->all(FLERR,"Ea keyword missing"); + if (isnan(eb)) + error->all(FLERR,"Eb keyword missing"); + if (isnan(ec)) + error->all(FLERR,"Ec keyword missing"); + if (isnan(k3)) + error->all(FLERR,"K3 keyword missing"); + if (isnan(e2)) + error->all(FLERR,"E2 keyword missing"); + if (isnan(k2)) + error->all(FLERR,"K2 keyword missing"); + + //now test that values are within acceptable ranges + if (k2 < 0.0 or k3 < 0.0) + error->all(FLERR,"exponential decay is negative"); + if (bondL<0.0) + error->all(FLERR,"OH bond length is negative"); + if (rc2 < 0.0 || rc3 < 0.0 || rs < 0.0) + error->all(FLERR,"potential cutoff is negative"); + if (rs > rc3) + error->all(FLERR,"potential switching distance is larger than cutoff"); + if (rs==rc3) + error->warning(FLERR,"potential switching distance is equal to cutoff: this is untested and not conserve energy"); + if (pairPerAtom<0) + error->all(FLERR,"neigh is negative"); +} diff --git a/src/USER-MISC/pair_e3b.h b/src/USER-MISC/pair_e3b.h new file mode 100644 index 0000000000..f247b51b7a --- /dev/null +++ b/src/USER-MISC/pair_e3b.h @@ -0,0 +1,74 @@ +/* -*- 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(e3b,PairE3B) + +#else + +#ifndef LMP_PAIR_E3B_H +#define LMP_PAIR_E3B_H + +#include "pair.h" +#include "compute_pe_e3b.h" + +namespace LAMMPS_NS { + +class PairE3B : public Pair { + public: + PairE3B(class LAMMPS *); + virtual ~PairE3B(); + virtual void compute(int, int); + void settings(int, char **); + virtual void coeff(int, char **); + virtual double init_one(int, int); + virtual void init_style(); + + //allow compute pe/e3b to access etot vector + friend void ComputePEE3B::compute_vector(); + +protected: + //potential parameters + int typeO; + double ea,eb,ec; //three body energies + double k3; //three body exponential decay (units inverse length) + double rs,rc3,rc2; //rs: switching cutuff, rc3: cutoff for 3-body + double e2,k2; //2-body energy and exp decay + double cutmax; //max cutoff of all interactions + double rc2sq,rc3sq,rc3deltaSq; + double sc_denom,sc_num; + + //list of indexes of Os and Hs in each pair + int pairmax,pairPerAtom; // size of pair list + int **pairO,***pairH; // pair lists + double ***exps,****del3,***fpair3,*sumExp; + int maxID; //size of global sumExp array + size_t nbytes; //size of sumExp array in bytes + int natoms; //to make sure number of atoms is constant + + double etot[4]; //etotA,etotB,etotC,etot2 + + virtual void allocate(); + void allocateE3B(); + bool allocatedE3B; + //for reading settings from pair_style input + bool checkKeyword(const char *,const char *,const int, const int); + void checkInputs(const double &bondL); + void presetParam(const int flag,bool &repeatFlag,double &bondL); + tagint find_maxID(); +}; +} + +#endif +#endif