jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

new CHARMM pair styles with force swithing/shifting

This commit is contained in:
Steve Plimpton 2017-03-24 13:53:07 -06:00
parent e6fcaefe95
commit 394e9b42b0
24 changed files with 3302 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/Section_commands.txt
View File

@ -940,6 +940,8 @@ KOKKOS, o = USER-OMP, t = OPT.
"lj/charmm/coul/charmm/implicit (ko)"_pair_charmm.html,
"lj/charmm/coul/long (giko)"_pair_charmm.html,
"lj/charmm/coul/msm"_pair_charmm.html,
"lj/charmmfsw/coul/charmmfsh"_pair_charmm.html,
"lj/charmmfsw/coul/long"_pair_charmm.html,
"lj/class2 (gko)"_pair_class2.html,
"lj/class2/coul/cut (ko)"_pair_class2.html,
"lj/class2/coul/long (gko)"_pair_class2.html,
@ -1148,6 +1150,7 @@ USER-OMP, t = OPT.
"zero"_dihedral_zero.html,
"hybrid"_dihedral_hybrid.html,
"charmm (ko)"_dihedral_charmm.html,
"charmmfsh"_dihedral_charmm.html,
"class2 (ko)"_dihedral_class2.html,
"harmonic (io)"_dihedral_harmonic.html,
"helix (o)"_dihedral_helix.html,

8
doc/src/Section_example.txt
View File

@ -51,9 +51,11 @@ Lowercase directories :h4
accelerate: run with various acceleration options (OpenMP, GPU, Phi)
balance: dynamic load balancing, 2d system
body: body particles, 2d system
cmap: CMAP 5-body contributions to CHARMM force field
colloid: big colloid particles in a small particle solvent, 2d system
comb: models using the COMB potential
coreshell: core/shell model using CORESHELL package
controller: use of fix controller as a thermostat
crack: crack propagation in a 2d solid
deposit: deposit atoms and molecules on a surface
dipole: point dipolar particles, 2d system
@ -62,6 +64,8 @@ eim: NaCl using the EIM potential
ellipse: ellipsoidal particles in spherical solvent, 2d system
flow: Couette and Poiseuille flow in a 2d channel
friction: frictional contact of spherical asperities between 2d surfaces
gcmc: Grand Canonical Monte Carlo (GCMC) via the fix gcmc command
granregion: use of fix wall/region/gran as boundary on granular particles
hugoniostat: Hugoniostat shock dynamics
indent: spherical indenter into a 2d solid
kim: use of potentials in Knowledge Base for Interatomic Models (KIM)
@ -69,6 +73,7 @@ meam: MEAM test for SiC and shear (same as shear examples)
melt: rapid melt of 3d LJ system
micelle: self-assembly of small lipid-like molecules into 2d bilayers
min: energy minimization of 2d LJ melt
mscg: parameterize a multi-scale coarse-graining (MSCG) model
msst: MSST shock dynamics
nb3b: use of nonbonded 3-body harmonic pair style
neb: nudged elastic band (NEB) calculation for barrier finding
@ -87,7 +92,8 @@ snap: NVE dynamics for BCC tantalum crystal using SNAP potential
srd: stochastic rotation dynamics (SRD) particles as solvent
streitz: use of Streitz/Mintmire potential with charge equilibration
tad: temperature-accelerated dynamics of vacancy diffusion in bulk Si
vashishta: use of the Vashishta potential :tb(s=:)
vashishta: use of the Vashishta potential
voronoi: Voronoi tesselation via compute voronoi/atom command :tb(s=:)
Here is how you can run and visualize one of the sample problems:

9
doc/src/Section_howto.txt
View File

@ -204,7 +204,10 @@ documentation for the formula it computes.
"bond_style"_bond_harmonic.html harmonic
"angle_style"_angle_charmm.html charmm
"dihedral_style"_dihedral_charmm.html charmmfsh
"dihedral_style"_dihedral_charmm.html charmm
"pair_style"_pair_charmm.html lj/charmmfsw/coul/charmmfsh
"pair_style"_pair_charmm.html lj/charmmfsw/coul/long
"pair_style"_pair_charmm.html lj/charmm/coul/charmm
"pair_style"_pair_charmm.html lj/charmm/coul/charmm/implicit
"pair_style"_pair_charmm.html lj/charmm/coul/long :ul
@ -212,6 +215,12 @@ documentation for the formula it computes.
"special_bonds"_special_bonds.html charmm
"special_bonds"_special_bonds.html amber :ul
NOTE: For CHARMM, the newer {charmmfsw} or {charmmfsh} styles were
released in March 2017. We recommend they be used instead of the
older {charmm} styles. See discussion of the differences on the "pair
charmm"_pair_charmm.html and "dihedral charmm"_dihedral_charmm.html
doc pages.
DREIDING is a generic force field developed by the "Goddard
group"_http://www.wag.caltech.edu at Caltech and is useful for
predicting structures and dynamics of organic, biological and

19
doc/src/compute_modify.txt
View File

@ -19,7 +19,7 @@ keyword = {extra/dof} or {extra} or {dynamic/dof} or {dynamic} :l
N = # of extra degrees of freedom to subtract
{extra} syntax is identical to {extra/dof}, will be disabled at some point
{dynamic/dof} value = {yes} or {no}
yes/no = do or do not recompute the number of atoms contributing to the temperature
yes/no = do or do not recompute the number of degrees of freedom (DOF) contributing to the temperature
{dynamic} syntax is identical to {dynamic/dof}, will be disabled at some point :pre
:ule
@ -46,13 +46,16 @@ degrees-of-freedom. See the "compute
temp/asphere"_compute_temp_asphere.html command for an example.
The {dynamic/dof} or {dynamic} keyword determines whether the number
of atoms N in the compute group is re-computed each time a temperature
is computed. Only compute styles that calculate a temperature use
this option. By default, N is assumed to be constant. If you are
adding atoms to the system (see the "fix pour"_fix_pour.html, "fix
deposit"_fix_deposit.html and "fix gcmc"_fix_gcmc.html commands) or
expect atoms to be lost (e.g. due to evaporation), then this option
should be used to insure the temperature is correctly normalized.
of atoms N in the compute group and their associated degrees of
freedom are re-computed each time a temperature is computed. Only
compute styles that calculate a temperature use this option. By
default, N and their DOF are assumed to be constant. If you are
adding atoms or molecules to the system (see the "fix
pour"_fix_pour.html, "fix deposit"_fix_deposit.html, and "fix
gcmc"_fix_gcmc.html commands) or expect atoms or molecules to be lost
(e.g. due to exiting the simulation box or via "fix
evaporation"_fix_evaporation.html), then this option should be used to
insure the temperature is correctly normalized.
NOTE: The {extra} and {dynamic} keywords should not be used as they
are deprecated (March 2017) and will eventually be disabled. Instead,

37
doc/src/dihedral_charmm.txt
View File

@ -10,21 +10,25 @@ dihedral_style charmm command :h3
dihedral_style charmm/intel command :h3
dihedral_style charmm/kk command :h3
dihedral_style charmm/omp command :h3
dihedral_style charmmfsh command :h3
[Syntax:]
dihedral_style charmm :pre
dihedral_style style :pre
style = {charmm} or {charmmfsh} :ul
[Examples:]
dihedral_style charmm
dihedral_style charmmfsh
dihedral_coeff 1 0.2 1 180 1.0
dihedral_coeff 2 1.8 1 0 1.0
dihedral_coeff 1 3.1 2 180 0.5 :pre
[Description:]
The {charmm} dihedral style uses the potential
The {charmm} and {charmmfsh} dihedral styles use the potential
:c,image(Eqs/dihedral_charmm.jpg)
@ -34,6 +38,11 @@ field (see comment on weighting factors below). See
"(Cornell)"_#dihedral-Cornell for a description of the AMBER force
field.
NOTE: The newer {charmmfsh} style was released in March 2017. We
recommend it be used instead of the older {charmm} style when running
a simulation with the CHARMM force field. See the discussion below
and more details on the "pair_style charmm"_pair_charmm.html doc page.
The following coefficients must be defined for each dihedral type via the
"dihedral_coeff"_dihedral_coeff.html command as in the example above, or in
the data file or restart files read by the "read_data"_read_data.html
@ -73,13 +82,23 @@ special_bonds 1-4 scaling factor to 0.0 (which is the
default). Otherwise 1-4 non-bonded interactions in dihedrals will be
computed twice.
Also note that for AMBER force fields, which use pair styles with
"lj/cut", the special_bonds 1-4 scaling factor should be set to the
AMBER defaults (1/2 and 5/6) and all the dihedral weighting factors
(4th coeff above) must be set to 0.0. In this case, you can use any
pair style you wish, since the dihedral does not need any
Lennard-Jones parameter information and will not compute any 1-4
non-bonded interactions.
For simulations using the CHARMM force field, the difference between
the {charmm} and {charmmfsh} styles is in the computation of the 1-4
non-bond interactions, if the distance between the two atoms is within
the switching distance of the pairwise potential defined by the
corresponding CHARMM pair style, i.e. between the inner and outer
cutoffs specified for the pair style. See the discussion on the
"CHARMM pair_style"_pair_charmm.html doc page for details.
Note that for AMBER force fields, which use pair styles with "lj/cut",
the special_bonds 1-4 scaling factor should be set to the AMBER
defaults (1/2 and 5/6) and all the dihedral weighting factors (4th
coeff above) must be set to 0.0. In this case, you can use any pair
style you wish, since the dihedral does not need any Lennard-Jones
parameter information and will not compute any 1-4 non-bonded
interactions. Likewise the {charmm} or {charmmfsh} styles are
identical in this case since no 1-4 non-bonded interactions are
computed.
:line

12
doc/src/kspace_style.txt
View File

@ -182,13 +182,13 @@ currently support the -DFFT_SINGLE compiler switch.
:line
The {msm} style invokes a multi-level summation method MSM solver,
"(Hardy)"_#Hardy2006 or "(Hardy2)"_#Hardy2009, which maps atom charge to a 3d
mesh, and uses a multi-level hierarchy of coarser and coarser meshes
on which direct coulomb solves are done. This method does not use
FFTs and scales as N. It may therefore be faster than the other
"(Hardy)"_#Hardy2006 or "(Hardy2)"_#Hardy2009, which maps atom charge
to a 3d mesh, and uses a multi-level hierarchy of coarser and coarser
meshes on which direct coulomb solves are done. This method does not
use FFTs and scales as N. It may therefore be faster than the other
K-space solvers for relatively large problems when running on large
core counts. MSM can also be used for non-periodic boundary conditions and
for mixed periodic and non-periodic boundaries.
core counts. MSM can also be used for non-periodic boundary conditions
and for mixed periodic and non-periodic boundaries.
MSM is most competitive versus Ewald and PPPM when only relatively
low accuracy forces, about 1e-4 relative error or less accurate,

138
doc/src/pair_charmm.txt
View File

@ -17,12 +17,14 @@ pair_style lj/charmm/coul/long/opt command :h3
pair_style lj/charmm/coul/long/omp command :h3
pair_style lj/charmm/coul/msm command :h3
pair_style lj/charmm/coul/msm/omp command :h3
pair_style lj/charmmfsw/coul/charmmfsh command :h3
pair_style lj/charmmfsw/coul/long command :h3
[Syntax:]
pair_style style args :pre
style = {lj/charmm/coul/charmm} or {lj/charmm/coul/charmm/implicit} or {lj/charmm/coul/long} or {lj/charmm/coul/msm}
style = {lj/charmm/coul/charmm} or {lj/charmm/coul/charmm/implicit} or {lj/charmm/coul/long} or {lj/charmm/coul/msm} or {lj/charmmfsw/coul/charmmfsh} or {lj/charmmfsw/coul/long}
args = list of arguments for a particular style :ul
{lj/charmm/coul/charmm} args = inner outer (inner2) (outer2)
inner, outer = global switching cutoffs for Lennard Jones (and Coulombic if only 2 args)
@ -35,12 +37,20 @@ args = list of arguments for a particular style :ul
cutoff = global cutoff for Coulombic (optional, outer is Coulombic cutoff if only 2 args)
{lj/charmm/coul/msm} args = inner outer (cutoff)
inner, outer = global switching cutoffs for LJ (and Coulombic if only 2 args)
cutoff = global cutoff for Coulombic (optional, outer is Coulombic cutoff if only 2 args)
{lj/charmmfsw/coul/charmmfsh} args = inner outer (cutoff)
inner, outer = global cutoffs for LJ (and Coulombic if only 2 args)
cutoff = global cutoff for Coulombic (optional, outer is Coulombic cutoff if only 2 args)
{lj/charmmfsw/coul/long} args = inner outer (cutoff)
inner, outer = global cutoffs for LJ (and Coulombic if only 2 args)
cutoff = global cutoff for Coulombic (optional, outer is Coulombic cutoff if only 2 args) :pre
[Examples:]
pair_style lj/charmm/coul/charmm 8.0 10.0
pair_style lj/charmm/coul/charmm 8.0 10.0 7.0 9.0
pair_style lj/charmmfsw/coul/charmmfsh 8.0 10.0
pair_style lj/charmmfsw/coul/charmmfsh 8.0 10.0 7.0 9.0
pair_coeff * * 100.0 2.0
pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
@ -51,6 +61,8 @@ pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
pair_style lj/charmm/coul/long 8.0 10.0
pair_style lj/charmm/coul/long 8.0 10.0 9.0
pair_style lj/charmmfsw/coul/long 8.0 10.0
pair_style lj/charmmfsw/coul/long 8.0 10.0 9.0
pair_coeff * * 100.0 2.0
pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
@ -61,21 +73,55 @@ pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
[Description:]
The {lj/charmm} styles compute LJ and Coulombic interactions with an
additional switching function S(r) that ramps the energy and force
smoothly to zero between an inner and outer cutoff. It is a widely
used potential in the "CHARMM"_http://www.scripps.edu/brooks MD code.
See "(MacKerell)"_#pair-MacKerell for a description of the CHARMM force
field.
These pair styles compute Lennard Jones (LJ) and Coulombic
interactions with additional switching or shifting functions that ramp
the energy and/or force smoothly to zero between an inner and outer
cutoff. They are implementations of the widely used CHARMM force
field used in the "CHARMM"_http://www.scripps.edu/brooks MD code (and
others). See "(MacKerell)"_#pair-MacKerell for a description of the
CHARMM force field.
The styles with {charmm} (not {charmmfsw} or {charmmfsh}) in their
name are the older, original LAMMPS implemetations. They compute the
LJ and Coulombic interactions with an energy switching function (esw,
a cubic polynomial), which ramps the energy smoothly to zero between
the inner and outer cutoff. This can cause irregularities in
pair-wise forces (due to the discontinuous 2nd derivative of energy at
the boundaries of the switching region), which in some cases can
result in substantial artifacts in an MD simulation.
The newer styles with {charmmfsw} or {charmmfsh} in their name replace
the energy switching with force switching (fsw) and force shifting
(fsh) functions, for LJ and Coulombic interactions respectively.
These follow the formulas and description given in
"(Steinbach)"_#Steinbach and "(Brooks)"_#Brooks to minimize these
artifacts.
NOTE: The newer {charmmfsw} or {charmmfsh} styles were released in
March 2017. We recommend they be used instead of the older {charmm}
styles. Eventually code from the new styles will propagate into the
related pair styles (e.g. implicit, accelerator, free energy
variants).
The general CHARMM formulas are as follows
:c,image(Eqs/pair_charmm.jpg)
Both the LJ and Coulombic terms require an inner and outer cutoff.
They can be the same for both formulas or different depending on
whether 2 or 4 arguments are used in the pair_style command. In each
case, the inner cutoff distance must be less than the outer cutoff.
It it typical to make the difference between the 2 cutoffs about 1.0
Angstrom.
where the specific form of the switching or shifting function S(r) is
given in the "(Steinbach)"_#Steinbach paper.
When using the {lj/charmm/coul/charmm styles}, both the LJ and
Coulombic terms require an inner and outer cutoff. They can be the
same for both formulas or different depending on whether 2 or 4
arguments are used in the pair_style command. For the
{lj/charmmfsw/coul/charmmfsh} style, the LJ term requires both an
inner and outer cutoff, while the Coulombic term requires only one
cutoff. If the Coulomb cutoff is not specified (2 instead of 3
arguments), the LJ outer cutoff is used for the Coulombic cutoff. In
all cases where an inner and outer cutoff are specified, the inner
cutoff distance must be less than the outer cutoff. It is typical to
make the difference between the inner and outer cutoffs about 2.0
Angstroms.
Style {lj/charmm/coul/charmm/implicit} computes the same formulas as
style {lj/charmm/coul/charmm} except that an additional 1/r term is
@ -86,12 +132,18 @@ screening. It is designed for use in a simulation of an unsolvated
biomolecule (no explicit water molecules).
Styles {lj/charmm/coul/long} and {lj/charmm/coul/msm} compute the same
formulas as style {lj/charmm/coul/charmm} except that an additional
damping factor is applied to the Coulombic term, as described for the
"lj/cut"_pair_lj.html pair styles. Only one Coulombic cutoff is
specified for {lj/charmm/coul/long} and {lj/charmm/coul/msm}; if only
2 arguments are used in the pair_style command, then the outer LJ
cutoff is used as the single Coulombic cutoff.
formulas as style {lj/charmm/coul/charmm} and style
{lj/charmmfsw/coul/long} computes the same formulas as style
{lj/charmmfsw/coul/charmmfsh}, except that an additional damping
factor is applied to the Coulombic term, so it can be used in
conjunction with the "kspace_style"_kspace_style.html command and its
{ewald} or {pppm} or {msm} option. Only one Coulombic cutoff is
specified for these styles; if only 2 arguments are used in the
pair_style command, then the outer LJ cutoff is used as the single
Coulombic cutoff. The Coulombic cutoff specified for these styles
means that pairwise interactions within this distance are computed
directly; interactions outside that distance are computed in
reciprocal space.
The following coefficients must be defined for each pair of atoms
types via the "pair_coeff"_pair_coeff.html command as in the examples
@ -111,7 +163,7 @@ sigma.
The latter 2 coefficients are optional. If they are specified, they
are used in the LJ formula between 2 atoms of these types which are
also first and fourth atoms in any dihedral. No cutoffs are specified
because this CHARMM force field does not allow varying cutoffs for
because the CHARMM force field does not allow varying cutoffs for
individual atom pairs; all pairs use the global cutoff(s) specified in
the pair_style command.
@ -148,38 +200,41 @@ mixed. The default mix value is {arithmetic} to coincide with the
usual settings for the CHARMM force field. See the "pair_modify"
command for details.
None of the lj/charmm pair styles support the
None of the {lj/charmm} or {lj/charmmfsw} pair styles support the
"pair_modify"_pair_modify.html shift option, since the Lennard-Jones
portion of the pair interaction is smoothed to 0.0 at the cutoff.
The {lj/charmm/coul/long} style supports the
"pair_modify"_pair_modify.html table option since it can tabulate the
short-range portion of the long-range Coulombic interaction.
The {lj/charmm/coul/long} and {lj/charmmfsw/coul/long} styles support
the "pair_modify"_pair_modify.html table option since they can
tabulate the short-range portion of the long-range Coulombic
interaction.
None of the lj/charmm pair styles support the
None of the {lj/charmm} or {lj/charmmfsw} pair styles support the
"pair_modify"_pair_modify.html tail option for adding long-range tail
corrections to energy and pressure, since the Lennard-Jones portion of
the pair interaction is smoothed to 0.0 at the cutoff.
All of the lj/charmm pair styles write their information to "binary
restart files"_restart.html, so pair_style and pair_coeff commands do
not need to be specified in an input script that reads a restart file.
All of the {lj/charmm} and {lj/charmmfsw} pair styles write their
information to "binary restart files"_restart.html, so pair_style and
pair_coeff commands do not need to be specified in an input script
that reads a restart file.
The lj/charmm/coul/long pair style supports the use of the {inner},
{middle}, and {outer} keywords of the "run_style respa"_run_style.html
command, meaning the pairwise forces can be partitioned by distance at
different levels of the rRESPA hierarchy. The other styles only
support the {pair} keyword of run_style respa. See the
"run_style"_run_style.html command for details.
The {lj/charmm/coul/long} and {lj/charmmfsw/coul/long} pair styles
support the use of the {inner}, {middle}, and {outer} keywords of the
"run_style respa"_run_style.html command, meaning the pairwise forces
can be partitioned by distance at different levels of the rRESPA
hierarchy. The other styles only support the {pair} keyword of
run_style respa. See the "run_style"_run_style.html command for
details.
:line
[Restrictions:]
The {lj/charmm/coul/charmm} and {lj/charmm/coul/charmm/implicit}
styles are part of the MOLECULE package. The {lj/charmm/coul/long}
style is part of the KSPACE package. They are only enabled if LAMMPS
was built with those packages. See the "Making
All the styles with {coul/charmm} or {coul/charmmfsh} styles are part
of the MOLECULE package. All the styles with {coul/long} style are
part of the KSPACE package. They are only enabled if LAMMPS was built
with those packages. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info. Note that
the MOLECULE and KSPACE packages are installed by default.
@ -191,6 +246,13 @@ the MOLECULE and KSPACE packages are installed by default.
:line
:link(Brooks)
[(Brooks)] Brooks, et al, J Comput Chem, 30, 1545 (2009).
:link(pair-MacKerell)
[(MacKerell)] MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
:link(Stenibach)
[(Steinbach)] Steinbach, Brooks, J Comput Chem, 15, 667 (1994).

19
doc/src/pair_momb.txt
View File

@ -25,18 +25,19 @@ pair_coeff 1 2 momb 0.0 1.0 1.0 10.2847 2.361 :pre
[Description:]
Style {momb} computes pairwise van der Waals (vdW) and short-range interactions using the
Morse potential and "(Grimme)"_#Grimme method implemented in the Many-Body
Metal-Organic (MOMB) force field described comprehensively in "(Fichthorn)"_#Fichthorn and
"(Zhou)"_#Zhou4. Grimme's method is widely used to correct for dispersion
in density functional theory calculations.
Style {momb} computes pairwise van der Waals (vdW) and short-range
interactions using the Morse potential and "(Grimme)"_#Grimme method
implemented in the Many-Body Metal-Organic (MOMB) force field
described comprehensively in "(Fichthorn)"_#Fichthorn and
"(Zhou)"_#Zhou4. Grimme's method is widely used to correct for
dispersion in density functional theory calculations.
:c,image(Eqs/pair_momb.jpg)
For the {momb} pair style, the following coefficients must be defined for each
pair of atoms types via the "pair_coeff"_pair_coeff.html command as in the
examples above, or in the data file or restart files read by the
"read_data"_read_data.html as described below:
For the {momb} pair style, the following coefficients must be defined
for each pair of atoms types via the "pair_coeff"_pair_coeff.html
command as in the examples above, or in the data file or restart files
read by the "read_data"_read_data.html as described below:
D0 (energy units)
alpha (1/distance units)

5
examples/README
View File

@ -63,8 +63,8 @@ balance: dynamic load balancing, 2d system
body: body particles, 2d system
cmap: CMAP 5-body contributions to CHARMM force field
colloid: big colloid particles in a small particle solvent, 2d system
coreshell: adiabatic core/shell model
comb: models using the COMB potential
coreshell: adiabatic core/shell model
controller: use of fix controller as a thermostat
crack: crack propagation in a 2d solid
deposit: deposition of atoms and molecules onto a 3d substrate
@ -74,6 +74,7 @@ eim: NaCl using the EIM potential
ellipse: ellipsoidal particles in spherical solvent, 2d system
flow: Couette and Poiseuille flow in a 2d channel
friction: frictional contact of spherical asperities between 2d surfaces
gcmc: Grand Canonical Monte Carlo (GCMC) via the fix gcmc command
granregion: use of fix wall/region/gran as boundary on granular particles
hugoniostat: Hugoniostat shock dynamics
indent: spherical indenter into a 2d solid
@ -103,7 +104,7 @@ snap: NVE dynamics for BCC tantalum crystal using SNAP potential
streitz: Streitz-Mintmire potential for Al2O3
tad: temperature-accelerated dynamics of vacancy diffusion in bulk Si
vashishta: models using the Vashishta potential
voronoi: test of Voronoi tesselation in compute voronoi/atom
voronoi: Voronoi tesselation via compute voronoi/atom command
Here is a src/Make.py command which will perform a parallel build of a
LAMMPS executable "lmp_mpi" with all the packages needed by all the

44
examples/gcmc/CO2.txt Normal file
View File

@ -0,0 +1,44 @@
# CO2 molecule file. TraPPE model.
3 atoms
2 bonds
1 angles
Coords
1 0.0 0.0 0.0
2 -1.16 0.0 0.0
3 1.16 0.0 0.0
Types
1 1
2 2
3 2
Charges
1 0.7
2 -0.35
3 -0.35
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

83
examples/gcmc/in.gcmc.co2 Normal file
View File

@ -0,0 +1,83 @@
# GCMC for CO2 molecular fluid, rigid/small/nvt dynamics
# Rigid CO2 TraPPE model
# [Potoff and J.I. Siepmann, Vapor-liquid equilibria of
# mixtures containing alkanes, carbon dioxide and
# nitrogen AIChE J., 47,1676-1682 (2001)].
# variables available on command line
variable mu index -8.1
variable disp index 0.5
variable temp index 338.0
variable lbox index 10.0
variable spacing index 5.0
# global model settings
units real
atom_style full
boundary p p p
pair_style lj/cut/coul/long 14
pair_modify mix arithmetic tail yes
kspace_style ewald 0.0001
bond_style harmonic
angle_style harmonic
# box, start molecules on simple cubic lattice
lattice sc ${spacing}
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} units box
create_box 2 box &
bond/types 1 &
angle/types 1 &
extra/bond/per/atom 2 &
extra/angle/per/atom 1 &
extra/special/per/atom 2
molecule co2mol CO2.txt
create_atoms 0 box mol co2mol 464563 units box
# rigid CO2 TraPPE model
pair_coeff 1 1 0.053649 2.8
pair_coeff 2 2 0.156973 3.05
bond_coeff 1 0 1.16
angle_coeff 1 0 180
# masses
mass 1 12.0107
mass 2 15.9994
# MD settings
group co2 type 1 2
neighbor 2.0 bin
neigh_modify every 1 delay 10 check yes
velocity all create ${temp} 54654
timestep 1.0
# rigid constraints with thermostat
fix myrigidnvt all rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
fix_modify myrigidnvt dynamic/dof no
# gcmc
variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix mygcmc all gcmc 100 100 100 0 54341 ${temp} ${mu} ${disp} mol &
co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
compute_modify thermo_temp dynamic/dof yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc
thermo 1000
# run
run 20000

42
examples/gcmc/in.gcmc.lj Normal file
View File

@ -0,0 +1,42 @@
# GCMC for LJ simple fluid, no dynamics
# variables available on command line
variable mu index -21.0
variable disp index 1.0
variable temp index 2.0
variable lbox index 10.0
# global model settings
units lj
atom_style atomic
pair_style lj/cut 3.0
pair_modify tail yes
# box
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
create_box 1 box
# lj parameters
pair_coeff * * 1.0 1.0
mass * 1.0
# gcmc
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
compute_modify thermo_temp dynamic yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
thermo 100
# run
run 1000

177
examples/gcmc/log.24Mar17.gcmc.co2.g++.1 Normal file
View File

@ -0,0 +1,177 @@
LAMMPS (17 Mar 2017)
# GCMC for CO2 molecular fluid, rigid/small/nvt dynamics
# Rigid CO2 TraPPE model
# [Potoff and J.I. Siepmann, Vapor-liquid equilibria of
# mixtures containing alkanes, carbon dioxide and
# nitrogen AIChE J., 47,1676-1682 (2001)].
# variables available on command line
variable mu index -8.1
variable disp index 0.5
variable temp index 338.0
variable lbox index 10.0
variable spacing index 5.0
# global model settings
units real
atom_style full
boundary p p p
pair_style lj/cut/coul/long 14
pair_modify mix arithmetic tail yes
kspace_style ewald 0.0001
bond_style harmonic
angle_style harmonic
# box, start molecules on simple cubic lattice
lattice sc ${spacing}
lattice sc 5.0
Lattice spacing in x,y,z = 5 5 5
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} units box
region box block 0 10.0 0 ${lbox} 0 ${lbox} units box
region box block 0 10.0 0 10.0 0 ${lbox} units box
region box block 0 10.0 0 10.0 0 10.0 units box
create_box 2 box bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
Created orthogonal box = (0 0 0) to (10 10 10)
1 by 1 by 1 MPI processor grid
molecule co2mol CO2.txt
Read molecule co2mol:
3 atoms with 2 types
2 bonds with 1 types
1 angles with 1 types
0 dihedrals with 0 types
0 impropers with 0 types
create_atoms 0 box mol co2mol 464563 units box
Created 24 atoms
# rigid CO2 TraPPE model
pair_coeff 1 1 0.053649 2.8
pair_coeff 2 2 0.156973 3.05
bond_coeff 1 0 1.16
angle_coeff 1 0 180
# masses
mass 1 12.0107
mass 2 15.9994
# MD settings
group co2 type 1 2
24 atoms in group co2
neighbor 2.0 bin
neigh_modify every 1 delay 10 check yes
velocity all create ${temp} 54654
velocity all create 338.0 54654
timestep 1.0
# rigid constraints with thermostat
fix myrigidnvt all rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
fix myrigidnvt all rigid/nvt/small molecule temp 338.0 ${temp} 100 mol co2mol
fix myrigidnvt all rigid/nvt/small molecule temp 338.0 338.0 100 mol co2mol
8 rigid bodies with 24 atoms
1.16 = max distance from body owner to body atom
fix_modify myrigidnvt dynamic/dof no
# gcmc
variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix mygcmc all gcmc 100 100 100 0 54341 ${temp} ${mu} ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
fix mygcmc all gcmc 100 100 100 0 54341 338.0 ${mu} ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 0.5 mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 0.5 mol co2mol tfac_insert 1.66666666666667 group co2 rigid myrigidnvt
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
compute_modify thermo_temp dynamic/dof yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc
thermo 1000
# run
run 20000
Ewald initialization ...
WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.164636
estimated absolute RMS force accuracy = 0.0332064
estimated relative force accuracy = 0.0001
KSpace vectors: actual max1d max3d = 16 2 62
kxmax kymax kzmax = 2 2 2
WARNING: Fix gcmc using full_energy option (../fix_gcmc.cpp:439)
0 atoms in group FixGCMC:gcmc_exclusion_group:mygcmc
0 atoms in group FixGCMC:rotation_gas_atoms:mygcmc
WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (../neighbor.cpp:472)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 16
ghost atom cutoff = 16
binsize = 8, bins = 2 2 2
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 15.61 | 15.61 | 15.61 Mbytes
Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_racc
0 364.27579 4238.8631 -9.6809388 13.391989 0.5846359 24 0 0 0 0
1000 311.39835 -327.93481 -8.6795381 9.9010062 0.51155641 21 0.13302848 0.12331626 0.6894397 0.90997852
WARNING: Using kspace solver on system with no charge (../kspace.cpp:289)
2000 905.66812 319.43347 -0.50350961 6.2991241 0.14615898 6 0.20952183 0.20430213 0.71797992 0.92626683
3000 275.57393 -719.89718 -26.534978 14.238181 0.80387436 33 0.21291069 0.20460696 0.72899202 0.9133259
4000 254.70771 -245.01902 -20.981537 13.160079 0.80387436 33 0.17245726 0.16974613 0.70145764 0.90542759
5000 96.073601 -517.98124 -34.019065 5.441166 0.87695385 36 0.14174575 0.13607057 0.6776754 0.90155771
6000 397.57265 148.92645 -7.2012893 10.665797 0.43847693 18 0.12299956 0.1202471 0.66165464 0.90274793
7000 455.4271 -347.44181 -5.9244703 12.217875 0.43847693 18 0.15182038 0.14791307 0.67904236 0.90560829
8000 301.03124 -627.45324 -13.251012 11.066909 0.5846359 24 0.16687346 0.16315516 0.6936719 0.91129375
9000 256.5747 -565.67983 -17.814128 11.981874 0.73079488 30 0.15458482 0.15131825 0.68966283 0.90993975
10000 443.60076 89.586912 -6.077863 11.900606 0.43847693 18 0.16092552 0.16020353 0.69882461 0.91422145
11000 436.43777 64.412921 -6.7128469 11.708443 0.43847693 18 0.17453966 0.17480683 0.70679243 0.91369445
12000 594.42207 849.07743 -3.3708621 10.040536 0.29231795 12 0.17461606 0.17568622 0.71175869 0.91333367
13000 426.85849 -1093.1334 -17.524618 17.813377 0.65771539 27 0.17742896 0.17792831 0.71363306 0.91450124
14000 317.75995 336.31107 -10.46774 11.681912 0.5846359 24 0.18331181 0.18427921 0.71715557 0.91652256
15000 272.65129 317.50536 -26.428336 14.087176 0.80387436 33 0.17449167 0.175957 0.71122398 0.91528038
16000 344.28567 -577.91079 -18.177927 16.077919 0.73079488 30 0.1661682 0.16781514 0.70485136 0.91508882
17000 134.55928 -193.5668 -30.297136 7.6208177 0.87695385 36 0.15965609 0.1605036 0.69658104 0.9140445
18000 231.87302 -446.07671 -14.875027 9.6763722 0.65771539 27 0.15270985 0.15351831 0.69002918 0.91372795
19000 328.6835 -280.22365 -20.001303 16.982214 0.80387436 33 0.15201017 0.15272181 0.69023195 0.91272534
20000 0 -20.39554 -0.14872889 -0 0 0 0.15600204 0.15750795 0.69503275 0.9138765
Loop time of 30.9008 on 1 procs for 20000 steps with 0 atoms
Performance: 55.921 ns/day, 0.429 hours/ns, 647.233 timesteps/s
99.8% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.1985 | 2.1985 | 2.1985 | 0.0 | 7.11
Bond | 0.029596 | 0.029596 | 0.029596 | 0.0 | 0.10
Kspace | 0.23123 | 0.23123 | 0.23123 | 0.0 | 0.75
Neigh | 0.16141 | 0.16141 | 0.16141 | 0.0 | 0.52
Comm | 0.20628 | 0.20628 | 0.20628 | 0.0 | 0.67
Output | 0.00068831 | 0.00068831 | 0.00068831 | 0.0 | 0.00
Modify | 28.022 | 28.022 | 28.022 | 0.0 | 90.69
Other | | 0.05058 | | | 0.16
Nlocal: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 0
Neighbor list builds = 40367
Dangerous builds = 118
Total wall time: 0:00:30

179
examples/gcmc/log.24Mar17.gcmc.co2.g++.4 Normal file
View File

@ -0,0 +1,179 @@
LAMMPS (17 Mar 2017)
# GCMC for CO2 molecular fluid, rigid/small/nvt dynamics
# Rigid CO2 TraPPE model
# [Potoff and J.I. Siepmann, Vapor-liquid equilibria of
# mixtures containing alkanes, carbon dioxide and
# nitrogen AIChE J., 47,1676-1682 (2001)].
# variables available on command line
variable mu index -8.1
variable disp index 0.5
variable temp index 338.0
variable lbox index 10.0
variable spacing index 5.0
# global model settings
units real
atom_style full
boundary p p p
pair_style lj/cut/coul/long 14
pair_modify mix arithmetic tail yes
kspace_style ewald 0.0001
bond_style harmonic
angle_style harmonic
# box, start molecules on simple cubic lattice
lattice sc ${spacing}
lattice sc 5.0
Lattice spacing in x,y,z = 5 5 5
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} units box
region box block 0 10.0 0 ${lbox} 0 ${lbox} units box
region box block 0 10.0 0 10.0 0 ${lbox} units box
region box block 0 10.0 0 10.0 0 10.0 units box
create_box 2 box bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
Created orthogonal box = (0 0 0) to (10 10 10)
1 by 2 by 2 MPI processor grid
molecule co2mol CO2.txt
Read molecule co2mol:
3 atoms with 2 types
2 bonds with 1 types
1 angles with 1 types
0 dihedrals with 0 types
0 impropers with 0 types
create_atoms 0 box mol co2mol 464563 units box
Created 24 atoms
# rigid CO2 TraPPE model
pair_coeff 1 1 0.053649 2.8
pair_coeff 2 2 0.156973 3.05
bond_coeff 1 0 1.16
angle_coeff 1 0 180
# masses
mass 1 12.0107
mass 2 15.9994
# MD settings
group co2 type 1 2
24 atoms in group co2
neighbor 2.0 bin
neigh_modify every 1 delay 10 check yes
velocity all create ${temp} 54654
velocity all create 338.0 54654
timestep 1.0
# rigid constraints with thermostat
fix myrigidnvt all rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
fix myrigidnvt all rigid/nvt/small molecule temp 338.0 ${temp} 100 mol co2mol
fix myrigidnvt all rigid/nvt/small molecule temp 338.0 338.0 100 mol co2mol
8 rigid bodies with 24 atoms
1.16 = max distance from body owner to body atom
fix_modify myrigidnvt dynamic/dof no
# gcmc
variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix mygcmc all gcmc 100 100 100 0 54341 ${temp} ${mu} ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
fix mygcmc all gcmc 100 100 100 0 54341 338.0 ${mu} ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 0.5 mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 0.5 mol co2mol tfac_insert 1.66666666666667 group co2 rigid myrigidnvt
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
compute_modify thermo_temp dynamic/dof yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc
thermo 1000
# run
run 20000
Ewald initialization ...
WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.164636
estimated absolute RMS force accuracy = 0.0332064
estimated relative force accuracy = 0.0001
KSpace vectors: actual max1d max3d = 16 2 62
kxmax kymax kzmax = 2 2 2
WARNING: Fix gcmc using full_energy option (../fix_gcmc.cpp:439)
0 atoms in group FixGCMC:gcmc_exclusion_group:mygcmc
0 atoms in group FixGCMC:rotation_gas_atoms:mygcmc
WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (../neighbor.cpp:472)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 16
ghost atom cutoff = 16
binsize = 8, bins = 2 2 2
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 15.4 | 15.4 | 15.4 Mbytes
Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_racc
0 386.52184 23582.465 -3.2433417 14.209828 0.5846359 24 0 0 0 0
WARNING: Using kspace solver on system with no charge (../kspace.cpp:289)
1000 760.80877 -39.270882 -3.5239626 12.851016 0.29231795 12 0.24161633 0.22984103 0.71087092 0.85283311
2000 308.0739 -255.061 -20.411926 14.386853 0.73079488 30 0.26075352 0.24898725 0.73128383 0.88590474
3000 432.34358 -1361.3278 -12.644057 15.894387 0.5846359 24 0.21121583 0.21051229 0.70036003 0.86735027
4000 631.524 -63.488785 -3.6517158 13.804656 0.36539744 15 0.22486443 0.22886173 0.72358173 0.87172606
5000 730.61244 -1029.284 -6.2144028 19.600352 0.43847693 18 0.23017182 0.22740779 0.72281887 0.87820845
6000 752.43412 503.4547 -3.7053679 16.447663 0.36539744 15 0.22943971 0.226183 0.71450085 0.87447436
7000 660.68448 828.51735 -10.592278 21.006666 0.51155641 21 0.24702096 0.24218506 0.71815602 0.8740222
8000 331.58822 -621.22187 -5.3705759 7.2482776 0.36539744 15 0.23211903 0.22906813 0.70281376 0.86269411
9000 413.91538 869.51669 -11.28701 15.216905 0.5846359 24 0.23246466 0.22923961 0.70832684 0.86244176
10000 242.20861 -808.23311 -5.4533937 5.2945044 0.36539744 15 0.22024676 0.22031775 0.70785097 0.85712561
11000 348.20046 -372.16895 -3.4663358 7.6114092 0.36539744 15 0.2252033 0.22688969 0.71513402 0.86123263
12000 251.99682 303.30092 -18.58289 11.768089 0.73079488 30 0.20916844 0.21068047 0.694787 0.84635875
13000 306.83592 -1582.0137 -20.810287 14.329041 0.73079488 30 0.19494837 0.196527 0.67554784 0.83056119
14000 476.57411 268.94927 -14.185859 19.888076 0.65771539 27 0.19779631 0.20016859 0.67957528 0.83375167
15000 267.03534 730.71183 -9.3348616 9.8171066 0.5846359 24 0.19468305 0.19814971 0.68032974 0.83810439
16000 639.83235 2190.3244 -9.6666503 26.701062 0.65771539 27 0.19520687 0.19848997 0.68514387 0.84100361
17000 2237.1203 -222.59057 -0.18248834 4.4456205 0.073079488 3 0.20412446 0.20757814 0.69175318 0.8434939
18000 754.44841 205.54694 -10.501943 27.736031 0.5846359 24 0.2129422 0.21508015 0.69665031 0.84758331
19000 1610.1148 1293.6003 -0.20849836 3.1996309 0.073079488 3 0.22061668 0.22356929 0.69949369 0.84851405
20000 231.61458 -39.696514 -4.6452226 5.0629266 0.36539744 15 0.21984893 0.22246517 0.69914635 0.85058457
Loop time of 21.1019 on 4 procs for 20000 steps with 15 atoms
Performance: 81.888 ns/day, 0.293 hours/ns, 947.781 timesteps/s
98.9% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.31897 | 0.41973 | 0.49748 | 10.1 | 1.99
Bond | 0.014808 | 0.015063 | 0.015289 | 0.2 | 0.07
Kspace | 0.3813 | 0.46228 | 0.56585 | 9.8 | 2.19
Neigh | 0.049173 | 0.050043 | 0.050868 | 0.3 | 0.24
Comm | 0.9755 | 0.99686 | 1.0205 | 1.9 | 4.72
Output | 0.0014546 | 0.0015051 | 0.0016098 | 0.2 | 0.01
Modify | 19.043 | 19.062 | 19.085 | 0.4 | 90.33
Other | | 0.09438 | | | 0.45
Nlocal: 3.75 ave 6 max 3 min
Histogram: 3 0 0 0 0 0 0 0 0 1
Nghost: 876.5 ave 937 max 818 min
Histogram: 1 1 0 0 0 0 0 0 1 1
Neighs: 490.5 ave 647 max 363 min
Histogram: 1 0 1 0 0 1 0 0 0 1
Total # of neighbors = 1962
Ave neighs/atom = 130.8
Ave special neighs/atom = 2
Neighbor list builds = 40070
Dangerous builds = 115
Total wall time: 0:00:21

103
examples/gcmc/log.24Mar17.gcmc.lj.g++.1 Normal file
View File

@ -0,0 +1,103 @@
LAMMPS (17 Mar 2017)
# GCMC for LJ simple fluid, no dynamics
# variables available on command line
variable mu index -21.0
variable disp index 1.0
variable temp index 2.0
variable lbox index 10.0
# global model settings
units lj
atom_style atomic
pair_style lj/cut 3.0
pair_modify tail yes
# box
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
region box block 0 10.0 0 ${lbox} 0 ${lbox}
region box block 0 10.0 0 10.0 0 ${lbox}
region box block 0 10.0 0 10.0 0 10.0
create_box 1 box
Created orthogonal box = (0 0 0) to (10 10 10)
1 by 1 by 1 MPI processor grid
# lj parameters
pair_coeff * * 1.0 1.0
mass * 1.0
# gcmc
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
compute_modify thermo_temp dynamic yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
thermo 100
# run
run 1000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 0.4369 | 0.4369 | 0.4369 Mbytes
Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc
0 0 0 0 -0 0 0 0 0 0
100 1.9042848 0.39026453 -1.7692765 2.8466449 0.292 292 0.3619855 0.30247792 0.40278761
200 1.8651924 0.47815517 -1.8494955 2.7886155 0.305 305 0.34021109 0.30357196 0.37759189
300 2.0626994 0.52068504 -1.8197295 3.0834166 0.291 291 0.32055605 0.3003043 0.36103862
400 2.0394818 0.53751435 -1.7636699 3.0482184 0.278 278 0.31698808 0.29995864 0.35441275
500 1.9628066 0.54594742 -1.7145336 2.9339513 0.287 287 0.31211861 0.29724228 0.35161407
600 1.9845913 0.40846162 -1.8199325 2.9669308 0.299 299 0.30976643 0.29612711 0.34933559
700 1.8582606 0.53445462 -1.7869306 2.777974 0.296 296 0.30642103 0.29446478 0.34633665
800 2.0340641 0.66057698 -1.7075279 3.0403148 0.283 283 0.30730979 0.29746793 0.34768045
900 2.0830765 0.63731971 -1.894775 3.114911 0.322 322 0.30636338 0.29737705 0.34737644
1000 1.9688933 0.50024802 -1.7013944 2.9428299 0.281 281 0.3053174 0.29772245 0.34788254
Loop time of 3.98286 on 1 procs for 1000 steps with 281 atoms
Performance: 108464.750 tau/day, 251.076 timesteps/s
99.9% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.10563 | 0.10563 | 0.10563 | 0.0 | 2.65
Neigh | 0.33428 | 0.33428 | 0.33428 | 0.0 | 8.39
Comm | 0.027969 | 0.027969 | 0.027969 | 0.0 | 0.70
Output | 0.00017285 | 0.00017285 | 0.00017285 | 0.0 | 0.00
Modify | 3.5096 | 3.5096 | 3.5096 | 0.0 | 88.12
Other | | 0.005197 | | | 0.13
Nlocal: 281 ave 281 max 281 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 977 ave 977 max 977 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 5902 ave 5902 max 5902 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 5902
Ave neighs/atom = 21.0036
Neighbor list builds = 1000
Dangerous builds = 0
Total wall time: 0:00:03

103
examples/gcmc/log.24Mar17.gcmc.lj.g++.4 Normal file
View File

@ -0,0 +1,103 @@
LAMMPS (17 Mar 2017)
# GCMC for LJ simple fluid, no dynamics
# variables available on command line
variable mu index -21.0
variable disp index 1.0
variable temp index 2.0
variable lbox index 10.0
# global model settings
units lj
atom_style atomic
pair_style lj/cut 3.0
pair_modify tail yes
# box
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
region box block 0 10.0 0 ${lbox} 0 ${lbox}
region box block 0 10.0 0 10.0 0 ${lbox}
region box block 0 10.0 0 10.0 0 10.0
create_box 1 box
Created orthogonal box = (0 0 0) to (10 10 10)
1 by 2 by 2 MPI processor grid
# lj parameters
pair_coeff * * 1.0 1.0
mass * 1.0
# gcmc
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
compute_modify thermo_temp dynamic yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
thermo 100
# run
run 1000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 0.434 | 0.434 | 0.434 Mbytes
Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc
0 0 0 0 -0 0 0 0 0 0
100 2.0328045 0.58661762 -1.6812724 3.0385824 0.287 287 0.35917318 0.30067507 0.38663622
200 1.9594279 0.50682399 -1.7308396 2.9287927 0.284 284 0.33788365 0.30337335 0.37300293
300 2.0602937 0.7028247 -1.9278541 3.0806296 0.315 315 0.31882007 0.29697498 0.36167185
400 1.995183 0.4328246 -1.8715454 2.983026 0.307 307 0.31527654 0.29681901 0.35673374
500 2.1390101 0.48232215 -1.554138 3.1960306 0.257 257 0.31372975 0.30003067 0.35558858
600 2.0584244 0.4929049 -1.6995569 3.0767263 0.283 283 0.31114213 0.29801665 0.35160109
700 1.9155066 0.49654243 -1.5770611 2.8624174 0.265 265 0.31056419 0.29944173 0.35157337
800 2.0883562 0.52731947 -1.8261112 3.1220925 0.3 300 0.30730979 0.29704354 0.34898892
900 2.0470677 0.5605993 -2.0130053 3.0610656 0.322 322 0.30484441 0.29586719 0.34678883
1000 2.004135 0.50642204 -1.6956257 2.9955798 0.283 283 0.30396929 0.29634309 0.34770304
Loop time of 3.688 on 4 procs for 1000 steps with 283 atoms
Performance: 117136.751 tau/day, 271.150 timesteps/s
99.2% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.024644 | 0.026027 | 0.027483 | 0.6 | 0.71
Neigh | 0.085449 | 0.088998 | 0.092893 | 0.9 | 2.41
Comm | 0.045756 | 0.051296 | 0.056578 | 1.7 | 1.39
Output | 0.00028491 | 0.00030857 | 0.00035262 | 0.0 | 0.01
Modify | 3.5189 | 3.5191 | 3.5194 | 0.0 | 95.42
Other | | 0.002221 | | | 0.06
Nlocal: 70.75 ave 77 max 68 min
Histogram: 1 2 0 0 0 0 0 0 0 1
Nghost: 514.25 ave 520 max 507 min
Histogram: 1 0 0 0 1 0 0 1 0 1
Neighs: 1483.5 ave 1715 max 1359 min
Histogram: 2 0 0 1 0 0 0 0 0 1
Total # of neighbors = 5934
Ave neighs/atom = 20.9682
Neighbor list builds = 1000
Dangerous builds = 0
Total wall time: 0:00:03

1078
src/KSPACE/pair_lj_charmmfsw_coul_long.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

110
src/KSPACE/pair_lj_charmmfsw_coul_long.h Normal file
View File

@ -0,0 +1,110 @@
/* -*- 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/charmmfsw/coul/long,PairLJCharmmfswCoulLong)
#else
#ifndef LMP_PAIR_LJ_CHARMMFSW_COUL_LONG_H
#define LMP_PAIR_LJ_CHARMMFSW_COUL_LONG_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJCharmmfswCoulLong : public Pair {
public:
PairLJCharmmfswCoulLong(class LAMMPS *);
virtual ~PairLJCharmmfswCoulLong();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
virtual void init_style();
void init_list(int, class NeighList *);
virtual double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
void compute_inner();
void compute_middle();
virtual void compute_outer(int, int);
virtual void *extract(const char *, int &);
protected:
int implicit;
int dihedflag;
double cut_lj_inner,cut_lj,cut_ljinv,cut_lj_innerinv;
double cut_lj_innersq,cut_ljsq;
double cut_lj3inv,cut_lj_inner3inv,cut_lj3,cut_lj_inner3;
double cut_lj6inv,cut_lj_inner6inv,cut_lj6,cut_lj_inner6;
double cut_coul,cut_coulsq;
double cut_bothsq;
double denom_lj,denom_lj12,denom_lj6;
double **epsilon,**sigma,**eps14,**sigma14;
double **lj1,**lj2,**lj3,**lj4,**offset;
double **lj14_1,**lj14_2,**lj14_3,**lj14_4;
double *cut_respa;
double g_ewald;
virtual void allocate();
};
}
#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/charmmfsw/coul/long requires atom attribute q
The atom style defined does not have these attributes.
E: Pair inner cutoff >= Pair outer cutoff
The specified cutoffs for the pair style are inconsistent.
E: Pair cutoff < Respa interior cutoff
One or more pairwise cutoffs are too short to use with the specified
rRESPA cutoffs.
E: Pair inner cutoff < Respa interior cutoff
One or more pairwise cutoffs are too short to use with the specified
rRESPA cutoffs.
E: Pair style requires a KSpace style
No kspace style is defined.
*/

6
src/MANYBODY/pair_airebo.cpp
View File

@ -3183,7 +3183,9 @@ double PairAIREBO::piRCSpline(double Nij, double Nji, double Nijconj,
for (i=0; i<64; i++) coeffs[i]=0.0;
if (typei==0 && typej==0) {
//if the inputs are out of bounds set them back to a point in bounds
// if the inputs are out of bounds set them back to a point in bounds
if (Nij<piCCdom[0][0]) Nij=piCCdom[0][0];
if (Nij>piCCdom[0][1]) Nij=piCCdom[0][1];
if (Nji<piCCdom[1][0]) Nji=piCCdom[1][0];
@ -3213,10 +3215,10 @@ double PairAIREBO::piRCSpline(double Nij, double Nji, double Nijconj,
}
}
// CH interaction
if ((typei==0 && typej==1) || (typei==1 && typej==0)) {
// if the inputs are out of bounds set them back to a point in bounds
if (Nij<piCHdom[0][0] || Nij>piCHdom[0][1] ||

482
src/MOLECULE/dihedral_charmmfsh.cpp Normal file
View File

@ -0,0 +1,482 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Paul Crozier (SNL)
The force-shifted sections were provided by Robert Meissner
and Lucio Colombi Ciacchi of Bremen University, Bremen, Germany,
with additional assistance from Robert A. Latour, Clemson University
------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include <stdlib.h>
#include "dihedral_charmmfsh.h"
#include "atom.h"
#include "comm.h"
#include "neighbor.h"
#include "domain.h"
#include "force.h"
#include "pair.h"
#include "update.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 0.05
/* ---------------------------------------------------------------------- */
DihedralCharmmfsh::DihedralCharmmfsh(LAMMPS *lmp) : Dihedral(lmp)
{
weightflag = 0;
writedata = 1;
}
/* ---------------------------------------------------------------------- */
DihedralCharmmfsh::~DihedralCharmmfsh()
{
if (allocated && !copymode) {
memory->destroy(setflag);
memory->destroy(k);
memory->destroy(multiplicity);
memory->destroy(shift);
memory->destroy(cos_shift);
memory->destroy(sin_shift);
memory->destroy(weight);
}
}
/* ---------------------------------------------------------------------- */
void DihedralCharmmfsh::compute(int eflag, int vflag)
{
int i1,i2,i3,i4,i,m,n,type;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
double edihedral,f1[3],f2[3],f3[3],f4[3];
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
double df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb;
double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz;
double c,s,p,sx2,sy2,sz2;
int itype,jtype;
double delx,dely,delz,rsq,r2inv,r6inv,r;
double forcecoul,forcelj,fpair,ecoul,evdwl;
edihedral = evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = 0;
// insure pair->ev_tally() will use 1-4 virial contribution
if (weightflag && vflag_global == 2)
force->pair->vflag_either = force->pair->vflag_global = 1;
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *atomtype = atom->type;
int **dihedrallist = neighbor->dihedrallist;
int ndihedrallist = neighbor->ndihedrallist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
double qqrd2e = force->qqrd2e;
for (n = 0; n < ndihedrallist; n++) {
i1 = dihedrallist[n][0];
i2 = dihedrallist[n][1];
i3 = dihedrallist[n][2];
i4 = dihedrallist[n][3];
type = dihedrallist[n][4];
// 1st bond
vb1x = x[i1][0] - x[i2][0];
vb1y = x[i1][1] - x[i2][1];
vb1z = x[i1][2] - x[i2][2];
// 2nd bond
vb2x = x[i3][0] - x[i2][0];
vb2y = x[i3][1] - x[i2][1];
vb2z = x[i3][2] - x[i2][2];
vb2xm = -vb2x;
vb2ym = -vb2y;
vb2zm = -vb2z;
// 3rd bond
vb3x = x[i4][0] - x[i3][0];
vb3y = x[i4][1] - x[i3][1];
vb3z = x[i4][2] - x[i3][2];
ax = vb1y*vb2zm - vb1z*vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm;
rasq = ax*ax + ay*ay + az*az;
rbsq = bx*bx + by*by + bz*bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
rg = sqrt(rgsq);
rginv = ra2inv = rb2inv = 0.0;
if (rg > 0) rginv = 1.0/rg;
if (rasq > 0) ra2inv = 1.0/rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq;
rabinv = sqrt(ra2inv*rb2inv);
c = (ax*bx + ay*by + az*bz)*rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
m = multiplicity[type];
p = 1.0;
ddf1 = df1 = 0.0;
for (i = 0; i < m; i++) {
ddf1 = p*c - df1*s;
df1 = p*s + df1*c;
p = ddf1;
}
p = p*cos_shift[type] + df1*sin_shift[type];
df1 = df1*cos_shift[type] - ddf1*sin_shift[type];
df1 *= -m;
p += 1.0;
if (m == 0) {
p = 1.0 + cos_shift[type];
df1 = 0.0;
}
if (eflag) edihedral = k[type] * p;
fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm;
hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm;
fga = fg*ra2inv*rginv;
hgb = hg*rb2inv*rginv;
gaa = -ra2inv*rg;
gbb = rb2inv*rg;
dtfx = gaa*ax;
dtfy = gaa*ay;
dtfz = gaa*az;
dtgx = fga*ax - hgb*bx;
dtgy = fga*ay - hgb*by;
dtgz = fga*az - hgb*bz;
dthx = gbb*bx;
dthy = gbb*by;
dthz = gbb*bz;
df = -k[type] * df1;
sx2 = df*dtgx;
sy2 = df*dtgy;
sz2 = df*dtgz;
f1[0] = df*dtfx;
f1[1] = df*dtfy;
f1[2] = df*dtfz;
f2[0] = sx2 - f1[0];
f2[1] = sy2 - f1[1];
f2[2] = sz2 - f1[2];
f4[0] = df*dthx;
f4[1] = df*dthy;
f4[2] = df*dthz;
f3[0] = -sx2 - f4[0];
f3[1] = -sy2 - f4[1];
f3[2] = -sz2 - f4[2];
// apply force to each of 4 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += f1[0];
f[i1][1] += f1[1];
f[i1][2] += f1[2];
}
if (newton_bond || i2 < nlocal) {
f[i2][0] += f2[0];
f[i2][1] += f2[1];
f[i2][2] += f2[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f3[0];
f[i3][1] += f3[1];
f[i3][2] += f3[2];
}
if (newton_bond || i4 < nlocal) {
f[i4][0] += f4[0];
f[i4][1] += f4[1];
f[i4][2] += f4[2];
}
if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
// 1-4 LJ and Coulomb interactions
// tally energy/virial in pair, using newton_bond as newton flag
if (weight[type] > 0.0) {
itype = atomtype[i1];
jtype = atomtype[i4];
delx = x[i1][0] - x[i4][0];
dely = x[i1][1] - x[i4][1];
delz = x[i1][2] - x[i4][2];
rsq = delx*delx + dely*dely + delz*delz;
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
// modifying coul and LJ force and energies to apply
// force_shift and force_switch as in CHARMM pairwise
// LJ interactions between 1-4 atoms should usually be
// for r < cut_inner, so switching not applied
r = sqrt(rsq);
if (implicit) forcecoul = qqrd2e * q[i1]*q[i4]*r2inv;
else if (dihedflag) forcecoul = qqrd2e * q[i1]*q[i4]*sqrt(r2inv);
else forcecoul = qqrd2e * q[i1]*q[i4]*(sqrt(r2inv) -
r*cut_coulinv14*cut_coulinv14);
forcelj = r6inv * (lj14_1[itype][jtype]*r6inv - lj14_2[itype][jtype]);
fpair = weight[type] * (forcelj+forcecoul)*r2inv;
if (eflag) {
if (dihedflag) ecoul = weight[type] * forcecoul;
else ecoul = weight[type] * qqrd2e * q[i1]*q[i4] *
(sqrt(r2inv) + r*cut_coulinv14*cut_coulinv14 -
2.0*cut_coulinv14);
evdwl14_12 = r6inv*lj14_3[itype][jtype]*r6inv -
lj14_3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
evdwl14_6 = -lj14_4[itype][jtype]*r6inv +
lj14_4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
evdwl = evdwl14_12 + evdwl14_6;
evdwl *= weight[type];
}
if (newton_bond || i1 < nlocal) {
f[i1][0] += delx*fpair;
f[i1][1] += dely*fpair;
f[i1][2] += delz*fpair;
}
if (newton_bond || i4 < nlocal) {
f[i4][0] -= delx*fpair;
f[i4][1] -= dely*fpair;
f[i4][2] -= delz*fpair;
}
if (evflag) force->pair->ev_tally(i1,i4,nlocal,newton_bond,
evdwl,ecoul,fpair,delx,dely,delz);
}
}
}
/* ---------------------------------------------------------------------- */
void DihedralCharmmfsh::allocate()
{
allocated = 1;
int n = atom->ndihedraltypes;
memory->create(k,n+1,"dihedral:k");
memory->create(multiplicity,n+1,"dihedral:k");
memory->create(shift,n+1,"dihedral:shift");
memory->create(cos_shift,n+1,"dihedral:cos_shift");
memory->create(sin_shift,n+1,"dihedral:sin_shift");
memory->create(weight,n+1,"dihedral:weight");
memory->create(setflag,n+1,"dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one type
------------------------------------------------------------------------- */
void DihedralCharmmfsh::coeff(int narg, char **arg)
{
if (narg != 5) error->all(FLERR,"Incorrect args for dihedral coefficients");
if (!allocated) allocate();
int ilo,ihi;
force->bounds(FLERR,arg[0],atom->ndihedraltypes,ilo,ihi);
// require integer values of shift for backwards compatibility
// arbitrary phase angle shift could be allowed, but would break
// backwards compatibility and is probably not needed
double k_one = force->numeric(FLERR,arg[1]);
int multiplicity_one = force->inumeric(FLERR,arg[2]);
int shift_one = force->inumeric(FLERR,arg[3]);
double weight_one = force->numeric(FLERR,arg[4]);
if (multiplicity_one < 0)
error->all(FLERR,"Incorrect multiplicity arg for dihedral coefficients");
if (weight_one < 0.0 || weight_one > 1.0)
error->all(FLERR,"Incorrect weight arg for dihedral coefficients");
if (weight_one > 0.0) weightflag=1;
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k[i] = k_one;
shift[i] = shift_one;
cos_shift[i] = cos(MY_PI*shift_one/180.0);
sin_shift[i] = sin(MY_PI*shift_one/180.0);
multiplicity[i] = multiplicity_one;
weight[i] = weight_one;
setflag[i] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients");
}
/* ----------------------------------------------------------------------
error check and initialize all values needed for force computation
------------------------------------------------------------------------- */
void DihedralCharmmfsh::init_style()
{
// insure use of CHARMM pair_style if any weight factors are non-zero
// set local ptrs to LJ 14 arrays setup by Pair
if (weightflag) {
int itmp;
if (force->pair == NULL)
error->all(FLERR,"Dihedral charmmfsh is incompatible with Pair style");
lj14_1 = (double **) force->pair->extract("lj14_1",itmp);
lj14_2 = (double **) force->pair->extract("lj14_2",itmp);
lj14_3 = (double **) force->pair->extract("lj14_3",itmp);
lj14_4 = (double **) force->pair->extract("lj14_4",itmp);
int *ptr = (int *) force->pair->extract("implicit",itmp);
if (!lj14_1 || !lj14_2 || !lj14_3 || !lj14_4 || !ptr)
error->all(FLERR,"Dihedral charmmfsh is incompatible with Pair style");
implicit = *ptr;
}
// constants for applying force switch (LJ) and force_shift (coul)
// to 1/4 dihedral atoms to match CHARMM pairwise interactions
int itmp;
int *p_dihedflag = (int *) force->pair->extract("dihedflag",itmp);
double *p_cutljinner = (double *) force->pair->extract("cut_lj_inner",itmp);
double *p_cutlj = (double *) force->pair->extract("cut_lj",itmp);
double *p_cutcoul = (double *) force->pair->extract("cut_coul",itmp);
if (p_cutcoul == NULL || p_cutljinner == NULL ||
p_cutlj == NULL || p_dihedflag == NULL)
error->all(FLERR,"Dihedral charmmfsh is incompatible with Pair style");
dihedflag = *p_dihedflag;
cut_coul14 = *p_cutcoul;
cut_lj_inner14 = *p_cutljinner;
cut_lj14 = *p_cutlj;
cut_coulinv14 = 1/cut_coul14;
cut_lj_inner3inv = (1/cut_lj_inner14) * (1/cut_lj_inner14) *
(1/cut_lj_inner14);
cut_lj_inner6inv = cut_lj_inner3inv * cut_lj_inner3inv;
cut_lj3inv = (1/cut_lj14) * (1/cut_lj14) * (1/cut_lj14);
cut_lj6inv = cut_lj3inv * cut_lj3inv;
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void DihedralCharmmfsh::write_restart(FILE *fp)
{
fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp);
fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp);
fwrite(&shift[1],sizeof(int),atom->ndihedraltypes,fp);
fwrite(&weight[1],sizeof(double),atom->ndihedraltypes,fp);
fwrite(&weightflag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void DihedralCharmmfsh::read_restart(FILE *fp)
{
allocate();
if (comm->me == 0) {
fread(&k[1],sizeof(double),atom->ndihedraltypes,fp);
fread(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp);
fread(&shift[1],sizeof(int),atom->ndihedraltypes,fp);
fread(&weight[1],sizeof(double),atom->ndihedraltypes,fp);
fread(&weightflag,sizeof(int),1,fp);
}
MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
MPI_Bcast(&multiplicity[1],atom->ndihedraltypes,MPI_INT,0,world);
MPI_Bcast(&shift[1],atom->ndihedraltypes,MPI_INT,0,world);
MPI_Bcast(&weight[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
MPI_Bcast(&weightflag,1,MPI_INT,0,world);
for (int i = 1; i <= atom->ndihedraltypes; i++) {
setflag[i] = 1;
cos_shift[i] = cos(MY_PI*shift[i]/180.0);
sin_shift[i] = sin(MY_PI*shift[i]/180.0);
}
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void DihedralCharmmfsh::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ndihedraltypes; i++)
fprintf(fp,"%d %g %d %d %g\n",i,k[i],multiplicity[i],shift[i],weight[i]);
}

81
src/MOLECULE/dihedral_charmmfsh.h Normal file
View File

@ -0,0 +1,81 @@
/* -*- 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 DIHEDRAL_CLASS
DihedralStyle(charmmfsh,DihedralCharmmfsh)
#else
#ifndef LMP_DIHEDRAL_CHARMMFSH_H
#define LMP_DIHEDRAL_CHARMMFSH_H
#include <stdio.h>
#include "dihedral.h"
namespace LAMMPS_NS {
class DihedralCharmmfsh : public Dihedral {
public:
DihedralCharmmfsh(class LAMMPS *);
virtual ~DihedralCharmmfsh();
virtual void compute(int, int);
virtual void coeff(int, char **);
virtual void init_style();
void write_restart(FILE *);
void read_restart(FILE *);
void write_data(FILE *);
protected:
int implicit,weightflag,dihedflag;
double cut_lj_inner14,cut_lj14,cut_coul14;
double evdwl14_12,evdwl14_6,cut_coulinv14;
double cut_lj_inner3inv,cut_lj_inner6inv,cut_lj3inv,cut_lj6inv;
double *k,*weight,*cos_shift,*sin_shift;
int *multiplicity,*shift;
double **lj14_1,**lj14_2,**lj14_3,**lj14_4;
virtual void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
W: Dihedral problem: %d %ld %d %d %d %d
Conformation of the 4 listed dihedral atoms is extreme; you may want
to check your simulation geometry.
E: Incorrect args for dihedral coefficients
Self-explanatory. Check the input script or data file.
E: Incorrect multiplicity arg for dihedral coefficients
Self-explanatory. Check the input script or data file.
E: Incorrect weight arg for dihedral coefficients
Self-explanatory. Check the input script or data file.
E: Dihedral charmmfsh is incompatible with Pair style
Dihedral style charmmfsh must be used with a pair style charmm
in order for the 1-4 epsilon/sigma parameters to be defined.
*/

546
src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp Normal file
View File

@ -0,0 +1,546 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Paul Crozier (SNL)
The lj-fsw/coul-fsh (force-switched and force-shifted) sections
were provided by Robert Meissner
and Lucio Colombi Ciacchi of Bremen University, Bremen, Germany,
with additional assistance from Robert A. Latour, Clemson University
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_lj_charmmfsw_coul_charmmfsh.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairLJCharmmfswCoulCharmmfsh::PairLJCharmmfswCoulCharmmfsh(LAMMPS *lmp) :
Pair(lmp)
{
implicit = 0;
mix_flag = ARITHMETIC;
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairLJCharmmfswCoulCharmmfsh::~PairLJCharmmfswCoulCharmmfsh()
{
if (!copymode) {
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(eps14);
memory->destroy(sigma14);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(lj14_1);
memory->destroy(lj14_2);
memory->destroy(lj14_3);
memory->destroy(lj14_4);
}
}
}
/* ---------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwl12,evdwl6,ecoul,fpair;
double r,rinv,r3inv,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
double switch1;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cut_bothsq) {
r2inv = 1.0/rsq;
r = sqrt(rsq);
if (rsq < cut_coulsq) {
forcecoul = qqrd2e * qtmp*q[j]*
(sqrt(r2inv) - r*cut_coulinv*cut_coulinv);
} else forcecoul = 0.0;
if (rsq < cut_ljsq) {
r6inv = r2inv*r2inv*r2inv;
jtype = type[j];
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
forcelj = forcelj*switch1;
}
} else forcelj = 0.0;
fpair = (factor_coul*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) {
ecoul = qqrd2e * qtmp*q[j]*
(sqrt(r2inv) + cut_coulinv*cut_coulinv*r - 2.0*cut_coulinv);
ecoul *= factor_coul;
} else ecoul = 0.0;
if (rsq < cut_ljsq) {
if (rsq > cut_lj_innersq) {
rinv = 1.0/r;
r3inv = rinv*rinv*rinv;
evdwl12 = lj3[itype][jtype]*cut_lj6*denom_lj12 *
(r6inv - cut_lj6inv)*(r6inv - cut_lj6inv);
evdwl6 = -lj4[itype][jtype]*cut_lj3*denom_lj6 *
(r3inv - cut_lj3inv)*(r3inv - cut_lj3inv);;
evdwl = evdwl12 + evdwl6;
} else {
evdwl12 = r6inv*lj3[itype][jtype]*r6inv -
lj3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
evdwl6 = -lj4[itype][jtype]*r6inv +
lj4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
evdwl = evdwl12 + evdwl6;
}
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();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(eps14,n+1,n+1,"pair:eps14");
memory->create(sigma14,n+1,n+1,"pair:sigma14");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(lj14_1,n+1,n+1,"pair:lj14_1");
memory->create(lj14_2,n+1,n+1,"pair:lj14_2");
memory->create(lj14_3,n+1,n+1,"pair:lj14_3");
memory->create(lj14_4,n+1,n+1,"pair:lj14_4");
}
/* ----------------------------------------------------------------------
global settings
unlike other pair styles,
there are no individual pair settings that these override
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::settings(int narg, char **arg)
{
if (narg != 2 && narg != 3)
error->all(FLERR,"Illegal pair_style command");
cut_lj_inner = force->numeric(FLERR,arg[0]);
cut_lj = force->numeric(FLERR,arg[1]);
if (narg == 2) {
cut_coul = cut_lj;
} else {
cut_coul = force->numeric(FLERR,arg[2]);
}
// indicates pair_style being used for dihedral_charmm
dihedflag = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::coeff(int narg, char **arg)
{
if (narg != 4 && narg != 6)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
double epsilon_one = force->numeric(FLERR,arg[2]);
double sigma_one = force->numeric(FLERR,arg[3]);
double eps14_one = epsilon_one;
double sigma14_one = sigma_one;
if (narg == 6) {
eps14_one = force->numeric(FLERR,arg[4]);
sigma14_one = force->numeric(FLERR,arg[5]);
}
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
eps14[i][j] = eps14_one;
sigma14[i][j] = sigma14_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style lj/charmmfsw/coul/charmmfsh "
"requires atom attribute q");
neighbor->request(this,instance_me);
// require cut_lj_inner < cut_lj
if (cut_lj_inner >= cut_lj)
error->all(FLERR,"Pair inner lj cutoff >= Pair outer lj cutoff");
cut_lj_innersq = cut_lj_inner * cut_lj_inner;
cut_ljsq = cut_lj * cut_lj;
cut_ljinv = 1.0/cut_lj;
cut_lj_innerinv = 1.0/cut_lj_inner;
cut_lj3 = cut_lj * cut_lj * cut_lj;
cut_lj3inv = cut_ljinv * cut_ljinv * cut_ljinv;
cut_lj_inner3inv = cut_lj_innerinv * cut_lj_innerinv * cut_lj_innerinv;
cut_lj_inner3 = cut_lj_inner * cut_lj_inner * cut_lj_inner;
cut_lj6 = cut_ljsq * cut_ljsq * cut_ljsq;
cut_lj6inv = cut_lj3inv * cut_lj3inv;
cut_lj_inner6inv = cut_lj_inner3inv * cut_lj_inner3inv;
cut_lj_inner6 = cut_lj_innersq * cut_lj_innersq * cut_lj_innersq;
cut_coulsq = cut_coul * cut_coul;
cut_coulinv = 1.0/cut_coul;
cut_bothsq = MAX(cut_ljsq,cut_coulsq);
denom_lj = (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) *
(cut_ljsq-cut_lj_innersq);
denom_lj12 = 1.0/(cut_lj6 - cut_lj_inner6);
denom_lj6 = 1.0/(cut_lj3 - cut_lj_inner3);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJCharmmfswCoulCharmmfsh::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
sigma[i][i],sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
eps14[i][j] = mix_energy(eps14[i][i],eps14[j][j],
sigma14[i][i],sigma14[j][j]);
sigma14[i][j] = mix_distance(sigma14[i][i],sigma14[j][j]);
}
double cut = MAX(cut_lj,cut_coul);
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj14_1[i][j] = 48.0 * eps14[i][j] * pow(sigma14[i][j],12.0);
lj14_2[i][j] = 24.0 * eps14[i][j] * pow(sigma14[i][j],6.0);
lj14_3[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],12.0);
lj14_4[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],6.0);
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
lj14_1[j][i] = lj14_1[i][j];
lj14_2[j][i] = lj14_2[i][j];
lj14_3[j][i] = lj14_3[i][j];
lj14_4[j][i] = lj14_4[i][j];
return cut;
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d %g %g %g %g\n",
i,epsilon[i][i],sigma[i][i],eps14[i][i],sigma14[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g %g %g\n",i,j,
epsilon[i][j],sigma[i][j],eps14[i][j],sigma14[i][j]);
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&eps14[i][j],sizeof(double),1,fp);
fwrite(&sigma14[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&eps14[i][j],sizeof(double),1,fp);
fread(&sigma14[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&eps14[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma14[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_inner,sizeof(double),1,fp);
fwrite(&cut_lj,sizeof(double),1,fp);
fwrite(&cut_coul,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJCharmmfswCoulCharmmfsh::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_lj_inner,sizeof(double),1,fp);
fread(&cut_lj,sizeof(double),1,fp);
fread(&cut_coul,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_lj_inner,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_lj,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairLJCharmmfswCoulCharmmfsh::
single(int i, int j, int itype, int jtype,
double rsq, double factor_coul, double factor_lj, double &fforce)
{
double r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj;
double phicoul,philj,philj12,philj6;
double switch1;
r2inv = 1.0/rsq;
r = sqrt(rsq);
rinv = 1.0/r;
if (rsq < cut_coulsq) {
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j] *
(sqrt(r2inv) - r*cut_coulinv*cut_coulinv);
} else forcecoul = 0.0;
if (rsq < cut_ljsq) {
r6inv = r2inv*r2inv*r2inv;
r3inv = rinv*rinv*rinv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
forcelj = forcelj*switch1;
}
} else forcelj = 0.0;
fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
double eng = 0.0;
if (rsq < cut_coulsq) {
phicoul = force->qqrd2e * atom->q[i]*atom->q[j] *
(sqrt(r2inv) + cut_coulinv*cut_coulinv*r - 2.0*cut_coulinv);
eng += factor_coul*phicoul;
}
if (rsq < cut_ljsq) {
if (rsq > cut_lj_innersq) {
philj12 = lj3[itype][jtype]*cut_lj6*denom_lj12 *
(r6inv - cut_lj6inv)*(r6inv - cut_lj6inv);
philj6 = -lj4[itype][jtype]*cut_lj3*denom_lj6 *
(r3inv - cut_lj3inv)*(r3inv - cut_lj3inv);;
philj = philj12 + philj6;
} else {
philj12 = r6inv*lj3[itype][jtype]*r6inv -
lj3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
philj6 = -lj4[itype][jtype]*r6inv +
lj4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
philj = philj12 + philj6;
}
eng += factor_lj*philj;
}
return eng;
}
/* ---------------------------------------------------------------------- */
void *PairLJCharmmfswCoulCharmmfsh::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"lj14_1") == 0) return (void *) lj14_1;
if (strcmp(str,"lj14_2") == 0) return (void *) lj14_2;
if (strcmp(str,"lj14_3") == 0) return (void *) lj14_3;
if (strcmp(str,"lj14_4") == 0) return (void *) lj14_4;
dim = 0;
if (strcmp(str,"implicit") == 0) return (void *) &implicit;
// info extracted by dihedral_charmmf
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
if (strcmp(str,"cut_lj_inner") == 0) return (void *) &cut_lj_inner;
if (strcmp(str,"cut_lj") == 0) return (void *) &cut_lj;
if (strcmp(str,"dihedflag") == 0) return (void *) &dihedflag;
return NULL;
}

88
src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.h Normal file
View File

@ -0,0 +1,88 @@
/* -*- 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/charmmfsw/coul/charmmfsh,PairLJCharmmfswCoulCharmmfsh)
#else
#ifndef LMP_PAIR_LJ_CHARMMFSW_COUL_CHARMMFSH_H
#define LMP_PAIR_LJ_CHARMMFSW_COUL_CHARMMFSH_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJCharmmfswCoulCharmmfsh : public Pair {
public:
PairLJCharmmfswCoulCharmmfsh(class LAMMPS *);
virtual ~PairLJCharmmfswCoulCharmmfsh();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
virtual void init_style();
virtual double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
virtual void *extract(const char *, int &);
protected:
int implicit;
int dihedflag;
double cut_lj_inner,cut_lj,cut_coul,cut_coulinv,cut_ljinv,cut_lj_innerinv;
double cut_lj_innersq,cut_ljsq,cut_coulsq,cut_bothsq;
double cut_lj3inv,cut_lj_inner3inv,cut_lj3,cut_lj_inner3;
double cut_lj6inv,cut_lj_inner6inv,cut_lj6,cut_lj_inner6;
double denom_lj,denom_lj12,denom_lj6;
double **epsilon,**sigma,**eps14,**sigma14;
double **lj1,**lj2,**lj3,**lj4;
double **lj14_1,**lj14_2,**lj14_3,**lj14_4;
virtual void allocate();
};
}
#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/charmmfsw/coul/charmmfsh requires atom attribute q
The atom style defined does not have these attributes.
E: Pair inner cutoff >= Pair outer cutoff
The specified cutoffs for the pair style are inconsistent.
*/

7
src/USER-MISC/pair_momb.cpp
View File

@ -9,8 +9,11 @@
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------
Contributed by Kristen Fichthorn, Tonnam Balankura, Ya Zhou @ Penn State University
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Kristen Fichthorn, Tonnam Balankura,
Ya Zhou (Penn State University)
------------------------------------------------------------------------- */
#include <math.h>