tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add e3b water model

Browse Source
This commit is contained in:
Steven Strong
2019-04-16 11:32:47 -05:00
parent 898860328b
commit 69e7a2a237
16 changed files with 1300 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_compute.txt
View File

@ -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,

1
doc/src/Commands_pair.txt
View File

@ -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,

15
doc/src/Eqs/e3b.tex Normal file
View File

@ -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 & r<R_s \\
\displaystyle\frac{(R_f-r)^2(R_f-3R_s+2r)}{(R_f-R_s)^3} & R_s\leq r\leq R_f \\
0 & r>R_f\\
\end{cases}
\end{align*}
\end{document}

1
doc/src/compute.txt
View File

@ -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 -

60
doc/src/compute_pe_e3b.txt Normal file
View File

@ -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

1
doc/src/computes.txt
View File

@ -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

138
doc/src/pair_e3b.txt Normal file
View File

@ -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)

1
doc/src/pair_style.txt
View File

@ -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

1
doc/src/pairs.txt
View File

@ -32,6 +32,7 @@ Pair Styles :h1
pair_dpd
pair_dpd_fdt
pair_dsmc
pair_e3b
pair_eam
pair_edip
pair_eff

10
examples/USER/e3b/README Normal file
View File

@ -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.

115
examples/USER/e3b/in.lammps Normal file
View File

@ -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

61
examples/USER/e3b/tip4p2005.mol Normal file
View File

@ -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

86
src/USER-MISC/compute_pe_e3b.cpp Normal file
View File

@ -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;
}

45
src/USER-MISC/compute_pe_e3b.h Normal file
View File

@ -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

690
src/USER-MISC/pair_e3b.cpp Normal file
View File

@ -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 <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
//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 rsq<rc2sq
//accumulate info about each OH pair for later 3body stuff
//test OO distance with augmented cutoff to account for dangling Hs
if (rsq < rc3deltaSq) {
//pairO and pairH are set here even if no Hs are within the cutoff
//in that case, npair is not incremented and they will be overwritten
pairO[npair][0] = i;
pairO[npair][1] = j;
addedH = false;
for (kk=0; kk<NUMO; kk++) {
k = pairO[npair][kk];
otherO = pairO[npair][(kk+1)%2];
for (hh=0; hh<NUMH; hh++) {
h=atom->map(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; kk<NUMO; kk++) {
i = pairO[ii][kk];
otherO = (kk+1)%2;
partB = eb*( sumExp[tag[pairO[ii][otherO] ]-1]
+ sumExp[tag[pairH[ii][otherO][0]]-1]
+ sumExp[tag[pairH[ii][otherO][1]]-1]
- 2*(exps[ii][otherO][0] + exps[ii][otherO][1]));
partC = ec*(sumExp[tag[i]-1] - exps[ii][kk][0] - exps[ii][kk][1]);
for (hh=0; hh<NUMH; hh++) {
j = pairH[ii][kk][hh];
//type A
otherH = (hh+1)%2;
j2 = pairH[ii][kk][otherH];
partA = ea*(sumExp[tag[j2]-1] - exps[ii][kk][otherH]); //not full energy yet
fpair = partA*fpair3[ii][kk][hh];
fxtmp = fpair*del3[ii][kk][hh][0];
fytmp = fpair*del3[ii][kk][hh][1];
fztmp = fpair*del3[ii][kk][hh][2];
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
f[j][0] -= fxtmp;
f[j][1] -= fytmp;
f[j][2] -= fztmp;
if (evflag) {
evdwl = partA*exps[ii][kk][hh]*0.5; //mult by exp on this H
ev_tally(i,j,nlocal,newton_pair,evdwl, 0.0,fpair,
del3[ii][kk][hh][0],del3[ii][kk][hh][1],del3[ii][kk][hh][2]);
etot[0] += evdwl;
}
//type B
fpair = partB*fpair3[ii][kk][hh];
fxtmp = fpair*del3[ii][kk][hh][0];
fytmp = fpair*del3[ii][kk][hh][1];
fztmp = fpair*del3[ii][kk][hh][2];
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
f[j][0] -= fxtmp;
f[j][1] -= fytmp;
f[j][2] -= fztmp;
if (evflag) {
evdwl = partB*exps[ii][kk][hh]*0.5; //mult by exp on this H
ev_tally(i,j,nlocal,newton_pair,evdwl, 0.0,fpair,
del3[ii][kk][hh][0],del3[ii][kk][hh][1],del3[ii][kk][hh][2]);
etot[1] += evdwl;
}
//type C
fpair = partC*fpair3[ii][kk][hh];
fxtmp = fpair*del3[ii][kk][hh][0];
fytmp = fpair*del3[ii][kk][hh][1];
fztmp = fpair*del3[ii][kk][hh][2];
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
f[j][0] -= fxtmp;
f[j][1] -= fytmp;
f[j][2] -= fztmp;
if (evflag) {
evdwl = partC*exps[ii][kk][hh]*0.5; //mult by exp on this H
ev_tally(i,j,nlocal,newton_pair,evdwl, 0.0,fpair,
del3[ii][kk][hh][0],del3[ii][kk][hh][1],del3[ii][kk][hh][2]);
etot[2] += evdwl;
}
} //end for hh in NUMH
} //end for kk in NUMO
} //end for ii in npairs
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairE3B::allocate()
{
allocated = 1;
int n = atom->ntypes;
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(nRem<nVal) {
char str[256];
snprintf(str,256,"Too few arguments to \"%s\" keyword.",test);
error->all(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");
}

74
src/USER-MISC/pair_e3b.h Normal file
View File

@ -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