diff --git a/doc/src/Commands_bond.rst b/doc/src/Commands_bond.rst index 1a0876e88f..0f6cac6a5b 100644 --- a/doc/src/Commands_bond.rst +++ b/doc/src/Commands_bond.rst @@ -71,6 +71,7 @@ OPT. * * * + * :doc:`amoeba ` * :doc:`charmm (iko) ` * :doc:`class2 (ko) ` * :doc:`class2/p6 ` @@ -149,6 +150,7 @@ OPT. * * * + * :doc:`amoeba ` * :doc:`class2 (ko) ` * :doc:`cossq (o) ` * :doc:`cvff (io) ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 45a75ff394..01beef5ea6 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -27,6 +27,8 @@ OPT. * :doc:`adapt/fep ` * :doc:`addforce ` * :doc:`addtorque ` + * :doc:`amoeba/bitorsion ` + * :doc:`amoeba/pitorsion ` * :doc:`append/atoms ` * :doc:`atc ` * :doc:`atom/swap ` diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 7cf4e7635b..3765937e1a 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -38,6 +38,7 @@ OPT. * :doc:`agni (o) ` * :doc:`airebo (io) ` * :doc:`airebo/morse (io) ` + * :doc:`amoeba ` * :doc:`atm ` * :doc:`awpmd/cut ` * :doc:`beck (go) ` @@ -122,6 +123,7 @@ OPT. * :doc:`hbond/dreiding/lj (o) ` * :doc:`hbond/dreiding/morse (o) ` * :doc:`hdnnp ` + * :doc:`hippo ` * :doc:`ilp/graphene/hbn ` * :doc:`kolmogorov/crespi/full ` * :doc:`kolmogorov/crespi/z ` diff --git a/doc/src/Howto_amoeba.rst b/doc/src/Howto_amoeba.rst index 9460d54200..6615e319a6 100644 --- a/doc/src/Howto_amoeba.rst +++ b/doc/src/Howto_amoeba.rst @@ -2,36 +2,87 @@ AMOEBA and HIPPO force fields ============================= The AMOEBA and HIPPO polarizable force fields were developed by Jay -Ponder's group at U Washington at St Louis. Their implementation in -LAMMPS was done using F90 code provided by the Ponder group from their -`Tinker MD code `_ +Ponder's group at the U Washington at St Louis. Their implementation +in LAMMPS was done using F90 code provided by the Ponder group from +their `Tinker MD code `_. -NOTE: what version of AMOEBA and HIPPO does LAMMPS implement? +NOTE: What version of AMOEBA and HIPPO does LAMMPS implement? These force fields can be used when polarization effects are desired in simulations of water, organic molecules, and biomolecules including -proteins, provided that parameterizations (force field files) are -available for the systems you are interested in. Files in the LAMMPS -potentials directory with a "amoeba" or "hippo" suffix can be used. -The Tinker distribution and website may have other force field files. +proteins, provided that parameterizations (Tinker PRM force field +files) are available for the systems you are interested in. Files in +the LAMMPS potentials directory with a "amoeba" or "hippo" suffix can +be used. The Tinker distribution and website have additional force +field files as well. + +NOTE: Can we include a few Tinker FF files in the LAMMPS distro, +e.g. ubiquitin.prm as potentials/ubiquitin.prm.amoeba ? Note that currently, HIPPO can only be used for water systems, but HIPPO files for a variety of small organic and biomolecules are in preparation by the Ponder group. Those force field files will be included in the LAMMPS distribution when available. -The :doc:`pair_style amoeba ` doc page gives a brief -description of the AMOEBA and HIPPO force fields. Further details for -AMOEBA are in these papers: :ref:`(Ren) `, :ref:`(Shi) -`. Further details for HIPPO are in this paper: -:ref:`(Rackers) `. +---------- + +The AMOEBA and HIPPO force fields contain the following terms in their +energy (U) computation. Further details for AMOEBA are in these +papers: :ref:`(Ren) `, :ref:`(Shi) `. Further +details for HIPPO are in this paper: :ref:`(Rackers) +`. + +.. math:: + + U & = U_{intermolecular} + U_{intramolecular) \\ + U_{intermolecular} & = U_{hal} + U_{repulsion} + U_{dispersion} + U_{multipole} + U_{polar} + U_{qxfer} \\ + U_{intramolecular} & = U_{bond} + U_{angle} + U_{torsion} + U_{oop} + U_{b\theta} + U_{Urey-Bradley} + U_{pitorsion} + U_{bitorsion} + +For intermolecular terms, the AMOEBA force field includes only the +:math:`U_{hal}`, :math:`U_{multipole}`, :math:`U_{polar}` terms. The +HIPPO force field includes all but the :math:`U_{hal}` term. In +LAMMPS, these are all computed by the :doc:`pair_style ` +command. + +For intramolecular terms, the :math:`U_{bond}`, :math:`U_{angle}`, +:math:`U_{torsion}`, :math:`U_{oop}` terms are computed by the +:doc:`bond_style class2 ` :doc:`angle_style amoeba +`, :doc:`dihedral_style fourier `, and +:doc:`improper_style amoeba ` commands respectively. +The :doc:`angle_style amoeba ` command includes the +:math:`U_{b\theta}` bond-angle cross term, and the +:math:`U_{Urey-Bradley}` term for a bond contribution between the 1-3 +atoms in the angle. + +The :math:`U_{pitorsion}` term is computed by the :doc:`fix +amoeba/pitorsion ` command. It computes 6-body +interaction between a pair of bonded atoms which each have 2 +additional bond partners. + +The :math:`U_{bitorsion}` term is computed by the :doc:`fix +amoeba/bitorsion ` command. It computes 5-body +interaction between two 4-body torsions (dihedrals) which overlap, +having 3 atoms in common. + +These command doc pages have additional details on the terms they +compute: + +* :doc:`pair_style amoeba or hippo ` +* :doc:`bond_style class2 ` +* :doc:`angle_style amoeba ` +* :doc:`dihedral_style fourier ` +* :doc:`improper_style amoeba ` +* :doc:`fix amoeba/pitorsion ` +* :doc:`fix amoeba/bitorsion ` ---------- -To use the AMOEBA force field in LAMMPS you should use command like -these appropriately in your input script. The only change needed for -a HIPPO simulation is for the pair_style and pair_coeff commands. See -examples/amoeba for example input scripts for both AMOEBA and HIPPO. +To use the AMOEBA or HIPPO force fields in LAMMPS, use commands like +the following appropriately in your input script. The only change +needed for AMOEBA vs HIPPO simulation is for the :doc:`pair_style +` and :doc:`pair_coeff ` commands, as shown +below. See examples/amoeba for example input scripts for both AMOEBA +and HIPPO. .. code-block:: LAMMPS @@ -66,73 +117,183 @@ examples/amoeba for example input scripts for both AMOEBA and HIPPO. special_bonds lj/coul 0.5 0.5 0.5 one/five yes # 1-5 neighbors +The data file read by the :doc:`read_data ` command should +be created by the tools/tinker/tinker2lmp.py conversion program +described below. It will create a section in the data file with the +header "Tinker Types". A :doc:`fix property/atom ` +command for the data must be specified befroe the read_data command. +In the example above the fix ID is amtype. -NOTE: that some options not needed for simpler systems +Similarly, if the system you are simulating defines AMOEBA/HIPPO +pitorsion or bitorsion interactions, there will be entries in the data +file for those interactions. They require a :doc:`fix +amoeba/pitortion ` and :doc:`fix +amoeba/bitorsion ` command be defined. In the +example above, the IDs for these two fixes are "pit" and "bit". -NOTE: tinker2lmp.py tool +Of course, if the system being modeled does not have one or more of +the following -- bond, angle, dihedral, improper, pitorision, +bitorsion interactions -- then the corresponding style and fix +commands above do not need to be used. See the example scripts in +examples/amoeba for water systems as examples; they are simpler than +what is listed above. -NOTE: are some commands not needed if system is simple and not ubi2 ? +The two :doc:`fix property/atom ` commands with IDs +"extra" and "polaxe" are also needed to define internal per-atom +quantities used by the AMOEBA and HIPPO force fields. -Both AMOEBA and HIPPO use amoeba for bond/angle/dihedral styles, -assuming the molecular system has bonds, angles, or dihedrals. These -style commands must be be used before the data file is read, as the -data file defines the coefficients for the various bond/angle/dihedral -types. The pair_style command can come after the read_data command, -as no pair coefficients are defined in the data file. +The :doc:`pair_coeff ` command used for either the AMOEBA +or HIPPO force field takes two arguments for Tinker force field files, +namely a PRM and KEY file. The keyfile can be specified as NULL and +default values for a various settings will be used. Note that these 2 +files are meant to allow use of native Tinker files as-is. However +LAMMPS does not support all the options which can be included +in a Tinker PRM or KEY file. See specifis below. -The :doc:`fix property/atom ` commands are required -to create per-atom data used by the AMOEBA/HIPPO force fields, some of -which is defined in the LAMMPS data file. The LAMMPS data file should -be produced as a pre-processing step by the tinker2lmp.py Python -script in tools/amoeba with a command like one of the following. -The tools/amoeba/README file has more details on using this tool. +A :doc:`special_bonds ` command with the "one/five" +option is required, since the AMOEBA/HIPPO force fields define +weighting factors for not only 1-2, 1-3, 1-4 interactions, but also +1-5 interactions. This command will trigger a per-atom list of 1-5 +neighbors to be generated. The AMOEBA and HIPPO force fields define +their own custom weighting factors for all the 1-2, 1-3, 1-4, 1-5 +terms which in the Tinker PRM and KEY files; they can be different for +different terms in the force field. -.. code-block:: bash +In addition to the list above, these command doc pages have additional +details: - % python tinker2lmp.py -xyz tinker.xyz -amoeba protein.prm.amoeba -data lmp.data # AMOEBA for non-periodic systems - % python tinker2lmp.py -xyz tinker.xyz -amoeba protein.prm.amoeba -data lmp.data -pbc 10.0 20.0 15.0 # AMOEBA for periodic systems - % python tinker2lmp.py -xyz tinker.xyz -hippo water.prm.hippo -data lmp.data # HIPPO for non-periodic systems - % python tinker2lmp.py -xyz tinker.xyz -hippo water.prm.hippo -data lmp.data -pbc 10.0 20.0 15.0 # HIPPO for periodic systems - -Note that two input files are needed and a LAMMPS data file (lmp.data) -is produced. The data file will have information on Tinker atom types -and AMOEBA/HIPPO force field parameters for bonds, angles, and -dihedrals. - -The first input is an XYZ file listing all the atoms in the -system. - -NOTE: is this a Tinker-augmented-XYZ format or standard? In either -case, how do we suggest LAMMPS users come up with these files? - -The second input is an AMOEBA or HIPPO PRM (force field) file. The -format of these files is defined by Tinker. A few such files are -provided in the LAMMPS potentials directory. Others may be available -in the Tinker distribution or from the Ponder group. - -The pair_coeff command should specify the same PRM file, and -optionally a Tinker-format KEY file. See the :doc:`pair_style amoeba -` doc page for more information about Tinker PRM and KEY -files. - -Finally, the :doc:`special_bonds ` command is used to -set all LJ and Coulombic 1-2, 1-3, 1-4 weighting factors to non-zero -and non-unity values, and to generate a per-atom list of 1-5 neighbors -as well. This is to insure all bond-topology neighbors are included -in the neighbor lists used by AMOEBA/HIPPO. These force fields apply -their own custom weighting factors to all these terms, including the -1-5 neighbors. +* :doc:`atom_style amoeba ` +* :doc:`fix property/atom ` +* :doc:`special_bonds ` ---------- -These command doc pages have additional details: +Tinker PRM and KEY files -* :doc:`pair_style amoeba or hippo ` -* :doc:`bond_style amoeba ` -* :doc:`angle_style amoeba ` -* :doc:`dihedral_style amoeba ` -* :doc:`fix property/atom ` -* :doc:`special_bonds ` +A Tinker PRM file is composed of sections, each of which has multiple +lines. This is the list of sections LAMMPS knows how to parse and +use. Any other sections are skipped: + +* Angle Bending Parameters +* Atom Type Definitions +* Atomic Multipole Parameters +* Bond Stretching Parameters +* Charge Penetration Parameters +* Charge Transfer Parameters +* Dipole Polarizability Parameters +* Dispersion Parameters +* Force Field Definition +* Literature References +* Out-of-Plane Bend Parameters +* Pauli Repulsion Parameters +* Pi-Torsion Parameters +* Stretch-Bend Parameters +* Torsion-Torsion Parameters +* Torsional Parameters +* Urey-Bradley Parameters +* Van der Waals Pair Parameters +* Van der Waals Parameters + +A Tinker KEY file is composed of lines, each of which has a keyword, +which can be followed by zero or more parameters. This is the list of +keywords LAMMPS knows how to parse and use in the same manner Tinker +does. Any other keywords are skipped. The value following the equal +sign is the default value for the keyword if it is not specified, or +if the keyfile in the :doc:`pair_coeff ` command is +specified as NULL: + +* a-axis = 0.0 +* b-axis = 0.0 +* c-axis = 0.0 +* ctrn-cutoff = 6.0 +* ctrn-taper = 0.9 * ctrn-cutoff +* cutoff +* delta-halgren = 0.07 +* dewald = no use of long-range dispersion unless specified +* dewald-alpha = 0.4 +* dewald-cutoff = 7.0 +* dispersion-cutoff = 9.0 +* dispersion-taper = 9.0 * dispersion-cutoff +* dpme-grid +* dpme-order = 4 +* ewald = no use of long-range electrostatics unless specified +* ewald-alpha = 0.4 +* ewald-cutoff = 7.0 +* gamma-halgren = 0.12 +* mpole-cutoff = 9.0 +* mpole-taper 0.65 * mpole-cutoff +* pcg-guess = enabled (default) +* pcg-noguess = disable pcg-guess +* pcg-noprecond = disable pcg-precond +* pcg-peek = 1.0 +* pcg-precond = enabled (default) +* pewald-alpha = 0.4 +* pme-grid +* pme-order = 5 +* polar-eps = 1.0e-6 +* polar-iter = 100 +* polar-predict = no prediction operation unless specified +* ppme-order = 5 +* repulsion-cutoff = 6.0 +* repulsion-taper = 0.9 * repulsion-cutoff +* taper +* usolve-cutoff = 4.5 +* usolve-diag = 2.0 +* vdw-cutoff = 9.0 +* vdw-taper = 0.9 * vdw-cutoff + +---------- + +Tinker2lmp.py tool + +This conversion tool is found in the tools/tinker directory. +As listed in examples/amoeba/README, these commands produce +the data files found in examples/amoeba, and also illustrate +all the options available to use with the tinker2lmp.py script: + +.. code-block:: bash + + % python tinker2lmp.py -xyz water_dimer.xyz -amoeba amoeba_water.prm -data data.water_dimer.amoeba # AMOEBA non-periodic system + % python tinker2lmp.py -xyz water_dimer.xyz -hippo hippo_water.prm -data data.water_dimer.hippo # HIPPO non-periodic system + % python tinker2lmp.py -xyz water_box.xyz -amoeba amoeba_water.prm -data data.water_box.amoeba -pbc 18.643 18.643 18.643 # AMOEBA periodic system + % python tinker2lmp.py -xyz water_box.xyz -hippo hippo_water.prm -data data.water_box.hippo -pbc 18.643 18.643 18.643 # HIPPO periodic system + % python tinker2lmp.py -xyz ubiquitin.xyz -amoeba amoeba_ubiquitin.prm -data data.ubiquitin.new -pbc 54.99 41.91 41.91 -bitorsion bitorsion.ubiquitin.data.new # system with bitorsions + +Switches and their arguments may be specified in any order. + +The -xyz switch is required and specifies an input XYZ file as an +argument. The format of this file is an extended XYZ format used by +Tinker for its input. Example *.xyz files are in the examples/amoeba +directory. The file lists the atoms in the system. Each atom has the +following information: Tinker species name (ignored by LAMMPS), xyz +coordinates, Tinker numeric type, and a list of atom IDs the atom is +bonded to. + +NOTE: is this a Tinker-unique augmented XYZ format or standard? Where +can a LAMMPS user get or generate this file for a system they want +to simulate? + +The -amoeba or -hippo switch is required. It specifies an input +AMOEBA or HIPPO PRM force field file as an argument. This should be +the same file used by the :doc:`pair_style ` command in +the input script. + +The -data switch is required. It specifies an output file name for +the LAMMPS data file that will be produced. + +For periodic systems, the -pbc switch is required. It specifies the +periodic box size for each dimension (x,y,z). For a Tinker simulation +these are specified in the KEY file. + +NOTE: What about a system with a free surface. What about a triclinic +box. + +The -bitorsion switch is only needed if the system contains Tinker +bitorsion interactions. The data for each type of bitorsion +interaction will be written to the specified file, and read by the +:doc:`fix amoeba/bitorsion ` command. The data +includes 2d arrays of values to which splines are fit, and thus is not +compatible with the LAMMPS data file format. ---------- diff --git a/doc/src/angle_amoeba.rst b/doc/src/angle_amoeba.rst new file mode 100644 index 0000000000..9a20299f45 --- /dev/null +++ b/doc/src/angle_amoeba.rst @@ -0,0 +1,137 @@ +.. index:: angle_style amoeba + +angle_style amoeba command +========================== + + +Syntax +"""""" + +.. code-block:: LAMMPS + + angle_style amoeba + +Examples +"""""""" + +.. code-block:: LAMMPS + + angle_style amoeba + angle_coeff * 75.0 -25.0 1.0 0.3 0.02 0.003 + angle_coeff * ba 3.6551 24.895 1.0119 1.5228 + angle_coeff * ub -7.6 1.5537 + +Description +""""""""""" + +The *amoeba* angle style uses the potential + +.. math:: + + E & = E_a + E_{ba} + E_{ub} \\ + E_a = K_2\left(\theta - \theta_0\right)^2 + K_3\left(\theta - \theta_0\right)^3 + K_4\left(\theta - \theta_0\right)^4 + K_5\left(\theta - \theta_0\right)^5 + K_6\left(\theta - \theta_0\right)^6 + E_{ba} & = N_1 (r_{ij} - r_1) (\theta - \theta_0) + N_2(r_{jk} - r_2)(\theta - \theta_0) + E_{ub} & = K_{ub} (r_{ik} - r_{ub})^2) + +where :math:`E_a` is the angle term, :math:`E_{ba}` is a bond-angle +term, ::math:`E_{ub}` is a Urey-Bradley bond term, :math:`\theta_0` is +the equilibrium angle, :math:`r_1` and :math:`r_2` are the equilibrium +bond lengths, and :math:`r_{ub}` is the equilibrium Urey-Bradley bond +length between atoms I and K in the angle IJK. + +These formulas match how the Tinker MD code does its angle +calculations for the AMOEBA and HIPPO force fields. See the +doc:`Howto amoeba ` doc page for more information about +the implemention of AMOEBA and HIPPO in LAMMPS. + +Note that the :math:`E_a` and :math:`E_{ba}` formulas are identical to +those used for the :doc:`angle_style class2/p6 ` +command, however there is no bond-bond cross term formula for +:math:`E_{bb}`. Additionally, there is a :math:`E_{ub}` term for a +Urey-Bradley bond. It is effectively a harmonic bond between the I +and K atoms of angle IJK, even though that bond is not enumerated in +the "Bonds" section of the data file. + +There are also two ways that Tinker computes the angle :math:`\theta` +in the :math:`E_a` formula. The first is the standard way of treating +IJK as an "in-plane" angle. The second is an "out-of-plane" method +which Tinker may use if the center atom J in the angle is bonded to +one additional atom in addition to I and K. In this case, all 4 atoms +are used to compute the :math:`E_a` formula, resulting in forces on +all 4 atoms. In the Tinker PRM file, these 2 options are denoted by +"angle" versus "anglep" entries in the "Angle Bending Parameters" +section of the PRM force field file. The *pflag* coefficient +described below selects between the 2 options. + +---------- + + +Coefficients for the :math:`E_a`, :math:`E_{bb}`, and :math:`E_{ub}` +formulas must be defined for each angle type via the :doc:`angle_coeff +` command as in the example above, or in the data file or +restart files read by the :doc:`read_data ` or +:doc:`read_restart ` commands. + +These are the 8 coefficients for the :math:`E_a` formula: + +* pflag = 0 or 1 +* ublag = 0 or 1 +* :math:`\theta_0` (degrees) +* :math:`K_2` (energy) +* :math:`K_3` (energy) +* :math:`K_4` (energy) +* :math:`K_5` (energy) +* :math:`K_6` (energy) + +A pflag value of 0 vs 1 selects between the "in-plane" and +"out-of-plane" options described above. Ubflag is 1 if there is a +Urey-Bradley term associated with this angle type, else it is 0. +:math:`\theta_0` is specified in degrees, but LAMMPS converts it to +radians internally; hence the various :math:`K` are effectively energy +per radian\^2 or radian\^3 or radian\^4 or radian\^5 or radian\^6. + +For the :math:`E_{ba}` formula, each line in a :doc:`angle_coeff +` command in the input script lists 5 coefficients, the +first of which is "ba" to indicate they are BondAngle coefficients. +In a data file, these coefficients should be listed under a "BondAngle +Coeffs" heading and you must leave out the "ba", i.e. only list 4 +coefficients after the angle type. + +* ba +* :math:`N_1` (energy/distance\^2) +* :math:`N_2` (energy/distance\^2) +* :math:`r_1` (distance) +* :math:`r_2` (distance) + +The :math:`\theta_0` value in the :math:`E_{ba}` formula is not specified, +since it is the same value from the :math:`E_a` formula. + +For the :math:`E_{ub}` formula, each line in a :doc:`angle_coeff +` command in the input script lists 3 coefficients, the +first of which is "ub" to indicate they are UreyBradley coefficients. +In a data file, these coefficients should be listed under a +"UreyBradley Coeffs" heading and you must leave out the "ub", +i.e. only list 2 coefficients after the angle type. + +* ub +* :math:`K_{ub}` (energy/distance\^2) +* :math:`r_{ub}` (distance) + +---------- + +Restrictions +"""""""""""" + +This angle style can only be used if LAMMPS was built with the AMOEBA +package. See the :doc:`Build package ` doc page for +more info. + +Related commands +"""""""""""""""" + +:doc:`angle_coeff ` + +Default +""""""" + +none diff --git a/doc/src/angle_class2.rst b/doc/src/angle_class2.rst index f257d96dc3..418df242aa 100644 --- a/doc/src/angle_class2.rst +++ b/doc/src/angle_class2.rst @@ -24,7 +24,7 @@ Examples .. code-block:: LAMMPS angle_style class2 - angle_coeff * 75.0 + angle_coeff * 75.0 25.0 0.3 0.002 angle_coeff 1 bb 10.5872 1.0119 1.5228 angle_coeff * ba 3.6551 24.895 1.0119 1.5228 diff --git a/doc/src/angle_style.rst b/doc/src/angle_style.rst index 0715625d45..88b135918a 100644 --- a/doc/src/angle_style.rst +++ b/doc/src/angle_style.rst @@ -73,6 +73,7 @@ of (g,i,k,o,t) to indicate which accelerated styles exist. * :doc:`zero ` - topology but no interactions * :doc:`hybrid ` - define multiple styles of angle interactions +* :doc:`amoeba ` - AMOEBA angle * :doc:`charmm ` - CHARMM angle * :doc:`class2 ` - COMPASS (class 2) angle * :doc:`class2/p6 ` - COMPASS (class 2) angle expanded to 6th order diff --git a/doc/src/atom_style.rst b/doc/src/atom_style.rst index bade8c2f79..55b661555b 100644 --- a/doc/src/atom_style.rst +++ b/doc/src/atom_style.rst @@ -10,7 +10,7 @@ Syntax atom_style style args -* style = *angle* or *atomic* or *body* or *bond* or *charge* or *dipole* or *dpd* or *edpd* or *electron* or *ellipsoid* or *full* or *line* or *mdpd* or *molecular* or *oxdna* or *peri* or *smd* or *sph* or *sphere* or *spin* or *tdpd* or *tri* or *template* or *hybrid* +* style = *amoeba* or *angle* or *atomic* or *body* or *bond* or *charge* or *dipole* or *dpd* or *edpd* or *electron* or *ellipsoid* or *full* or *line* or *mdpd* or *molecular* or *oxdna* or *peri* or *smd* or *sph* or *sphere* or *spin* or *tdpd* or *tri* or *template* or *hybrid* .. parsed-literal:: @@ -77,6 +77,8 @@ coordinates, velocities, atom IDs and types. See the :doc:`set ` commands for info on how to set these various quantities. ++--------------+-----------------------------------------------------+--------------------------------------+ +| *amoeba* | molecular + charge + 1/5 neighbors | AMOEBA/HIPPO polarized force fields | +--------------+-----------------------------------------------------+--------------------------------------+ | *angle* | bonds and angles | bead-spring polymers with stiffness | +--------------+-----------------------------------------------------+--------------------------------------+ @@ -134,15 +136,18 @@ quantities. .. note:: It is possible to add some attributes, such as a molecule ID, to - atom styles that do not have them via the :doc:`fix property/atom ` command. This command also - allows new custom attributes consisting of extra integer or - floating-point values to be added to atoms. See the :doc:`fix property/atom ` page for examples of cases - where this is useful and details on how to initialize, access, and - output the custom values. + atom styles that do not have them via the :doc:`fix property/atom + ` command. This command also allows new custom + attributes consisting of extra integer or floating-point values to + be added to atoms. See the :doc:`fix property/atom + ` page for examples of cases where this is + useful and details on how to initialize, access, and output the + custom values. All of the above styles define point particles, except the *sphere*, *ellipsoid*, *electron*, *peri*, *wavepacket*, *line*, *tri*, and -*body* styles, which define finite-size particles. See the :doc:`Howto spherical ` page for an overview of using +*body* styles, which define finite-size particles. See the +:doc:`Howto spherical ` page for an overview of using finite-size particle models with LAMMPS. All of the point-particle styles assign mass to particles on a @@ -266,16 +271,17 @@ showing the use of the *template* atom style versus *molecular*. .. note:: - When using the *template* style with a :doc:`molecule template ` that contains multiple molecules, you should - insure the atom types, bond types, angle_types, etc in all the - molecules are consistent. E.g. if one molecule represents H2O and - another CO2, then you probably do not want each molecule file to - define 2 atom types and a single bond type, because they will conflict - with each other when a mixture system of H2O and CO2 molecules is - defined, e.g. by the :doc:`read_data ` command. Rather the - H2O molecule should define atom types 1 and 2, and bond type 1. And - the CO2 molecule should define atom types 3 and 4 (or atom types 3 and - 2 if a single oxygen type is desired), and bond type 2. + When using the *template* style with a :doc:`molecule template + ` that contains multiple molecules, you should insure the + atom types, bond types, angle_types, etc in all the molecules are + consistent. E.g. if one molecule represents H2O and another CO2, + then you probably do not want each molecule file to define 2 atom + types and a single bond type, because they will conflict with each + other when a mixture system of H2O and CO2 molecules is defined, + e.g. by the :doc:`read_data ` command. Rather the H2O + molecule should define atom types 1 and 2, and bond type 1. And + the CO2 molecule should define atom types 3 and 4 (or atom types 3 + and 2 if a single oxygen type is desired), and bond type 2. For the *body* style, the particles are arbitrary bodies with internal attributes defined by the "style" of the bodies, which is specified by @@ -352,6 +358,8 @@ Many of the styles listed above are only enabled if LAMMPS was built with a specific package, as listed below. See the :doc:`Build package ` page for more info. +The *amoeba* style is part of the AMOEBA package. + The *angle*, *bond*, *full*, *molecular*, and *template* styles are part of the MOLECULE package. @@ -363,9 +371,11 @@ The *dipole* style is part of the DIPOLE package. The *peri* style is part of the PERI package for Peridynamics. -The *oxdna* style is part of the CG-DNA package for coarse-grained simulation of DNA and RNA. +The *oxdna* style is part of the CG-DNA package for coarse-grained +simulation of DNA and RNA. -The *electron* style is part of the EFF package for :doc:`electronic force fields `. +The *electron* style is part of the EFF package for :doc:`electronic +force fields `. The *dpd* style is part of the DPD-REACT package for dissipative particle dynamics (DPD). @@ -376,7 +386,8 @@ dissipative particle dynamics (mDPD), and transport dissipative particle dynamics (tDPD), respectively. The *sph* style is part of the SPH package for smoothed particle -hydrodynamics (SPH). See `this PDF guide `_ to using SPH in LAMMPS. +hydrodynamics (SPH). See `this PDF guide +`_ to using SPH in LAMMPS. The *mesont* style is part of the MESONT package. diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 81c0c87320..9e9089b312 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -170,6 +170,8 @@ accelerated styles exist. * :doc:`adapt/fep ` - enhanced version of fix adapt * :doc:`addforce ` - add a force to each atom * :doc:`addtorque ` - add a torque to a group of atoms +* :doc:`amoeba/bitorsion ` - torsion/torsion terms in AMOEBA force field +* :doc:`amoeba/pitorsion ` - 6-body terms in AMOEBA force field * :doc:`append/atoms ` - append atoms to a running simulation * :doc:`atc ` - initiates a coupled MD/FE simulation * :doc:`atom/swap ` - Monte Carlo atom type swapping @@ -194,7 +196,7 @@ accelerated styles exist. * :doc:`box/relax ` - relax box size during energy minimization * :doc:`charge/regulation ` - Monte Carlo sampling of charge regulation * :doc:`client/md ` - MD client for client/server simulations -* :doc:`cmap ` - enables CMAP cross-terms of the CHARMM force field +* :doc:`cmap ` - CMAP torsion/torsion terms in CHARMM force field * :doc:`colvars ` - interface to the collective variables "Colvars" library * :doc:`controller ` - apply control loop feedback mechanism * :doc:`deform ` - change the simulation box size/shape diff --git a/doc/src/fix_amoeba_bitorsion.rst b/doc/src/fix_amoeba_bitorsion.rst new file mode 100644 index 0000000000..b05c2128d6 --- /dev/null +++ b/doc/src/fix_amoeba_bitorsion.rst @@ -0,0 +1,167 @@ +.. index:: fix amoeba/bitorsion + +fix amoeba/bitorsion command +================ + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID ameoba/bitorsion filename + +* ID, group-ID are documented in :doc:`fix ` command +* amoeba/bitorsion = style name of this fix command +* filename = force-field file with AMOEBA bitorsion coefficients + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix bit all amoeba/bitorsion bitorsion.ubiquitin.data + read_data proteinX.data fix bit bitorsions BiTorsions + fix_modify bit energy yes + +Description +""""""""""" + +This command enables 5-body torsion/torsion interactions to be added +to simulations which use the AMOEBA and HIPPO force fields. It +matches how the Tinker MD code computes its torsion/torsion +interactions for the AMOEBA and HIPPO force fields. See the +doc:`Howto amoeba ` doc page for more information about +the implemention of AMOEBA and HIPPO in LAMMPS. + +Bitorsion interactions add additional potential energy contributions +to pairs of overlapping phi-psi dihedrals of amino-acids, which are +important to properly represent their conformational behavior. + +The examples/amoeba directory has a sample input script and data file +for ubiquitin, which illustrates use of the fix amoeba/bitorsion +command. + +As in the example above, this fix should be used before reading a data +file that contains a listing of bitorsion interactions. The +*filename* specified should contain the bitorsion parameters for the +AMOEBA or HIPPO force field. + +The data file read by the "read_data" must contain the topology of all +the bitorsion interactions, similar to the topology data for bonds, +angles, dihedrals, etc. Specially it should have a line like this in +its header section: + +.. parsed-literal:: + + N bitorsions + +where N is the number of bitorsion 5-body interactions. It should +also have a section in the body of the data file like this with N +lines: + +.. parsed-literal:: + + BiTorsions + + 1 1 8 10 12 18 20 + 2 5 18 20 22 25 27 + [...] + N 3 314 315 317 318 330 + +The first column is an index from 1 to N to enumerate the bitorsion +5-atom tuples; it is ignored by LAMMPS. The second column is the +"type" of the interaction; it is an index into the bitorsion force +field file. The remaining 5 columns are the atom IDs of the atoms in +the two 4-atom dihedrals that overlap to create the bitorsion 5-body +interaction. Note that the "bitorsions" and "BiTorsions" keywords for +the header and body sections match those specified in the read_data +command following the data file name; see the :doc:`read_data +` page for more details. + +The data file should be generated by using the +tools/tinker/tinker2lmp.py conversion script which creates a LAMMPS +data file from Tinker input files, including its PRM file which +contains the parameters necessary for computing bitorsion +interactions. The script must be invoked with the optional +"-bitorsion" flag to do this; see the example for the ubiquitin system +in the tools/tinker/README file. The same conversion script also +creates the file of bitorsion coefficient data which is read by this +command. + +The potential energy associated with bitorsion interactions can be +output as described below. It can also be included in the total +potential energy of the system, as output by the :doc:`thermo_style +` command, if the :doc:`fix_modify energy ` +command is used, as in the example above. See the note below about +how to include the bitorsion energy when performing an :doc:`energy +minimization `. + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +This fix writes the list of bitorsion interactions to :doc:`binary +restart files `. See the :doc:`read_restart ` +command for info on how to re-specify a fix in an input script that +reads a restart file, so that the operation of the fix continues in an +uninterrupted fashion. + +The :doc:`fix_modify ` *energy* option is supported by +this fix to add the potential energy of the bitorsion interactions to +both the global potential energy and peratom potential energies of the +system as part of :doc:`thermodynamic output ` or output +by the :doc:`compute pe/atom ` command. The default +setting for this fix is :doc:`fix_modify energy yes `. + +The :doc:`fix_modify ` *virial* option is supported by +this fix to add the contribution due to the bitorsion interactions to +both the global pressure and per-atom stress of the system via the +:doc:`compute pressure ` and :doc:`compute +stress/atom ` commands. The former can be +accessed by :doc:`thermodynamic output `. The default +setting for this fix is :doc:`fix_modify virial yes `. + +This fix computes a global scalar which can be accessed by various +:doc:`output commands `. The scalar is the potential +energy discussed above. The scalar value calculated by this fix is +"extensive". + +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. + +The forces due to this fix are imposed during an energy minimization, +invoked by the :doc:`minimize ` command. + +The :doc:`fix_modify ` *respa* option is supported by this +fix. This allows to set at which level of the :doc:`r-RESPA +` integrator the fix is adding its forces. Default is the +outermost level. + +.. note:: + + For energy minimization, if you want the potential energy + associated with the bitorsion terms forces to be included in the + total potential energy of the system (the quantity being + minimized), you MUST not disable the :doc:`fix_modify ` + *energy* option for this fix. + +Restrictions +"""""""""""" + +To function as expected this fix command must be issued *before* a +:doc:`read_data ` command but *after* a :doc:`read_restart +` command. + +This fix can only be used if LAMMPS was built with the AMOEBA package. +See the :doc:`Build package ` page for more info. + +Related commands +"""""""""""""""" + +:doc:`fix_modify `, :doc:`read_data ` + +Default +""""""" + +none diff --git a/doc/src/fix_amoeba_pitorsion.rst b/doc/src/fix_amoeba_pitorsion.rst new file mode 100644 index 0000000000..b80ca3a285 --- /dev/null +++ b/doc/src/fix_amoeba_pitorsion.rst @@ -0,0 +1,186 @@ +.. index:: fix amoeba/pitorsion + +fix amoeba/pitorsion command +================ + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID ameoba/pitorsion + +* ID, group-ID are documented in :doc:`fix ` command +* amoeba/pitorsion = style name of this fix command + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix pit all amoeba/pitorsion + read_data proteinX.data fix pit "pitorsion types" "PiTorsion Coeffs" & + fix pit pitorsions PiTorsions + fix_modify pit energy yes + +Description +""""""""""" + +This command enables 6-body pitorsion interactions to be added to +simulations which use the AMOEBA and HIPPO force fields. It matches +how the Tinker MD code computes its pitorsion interactions for the +AMOEBA and HIPPO force fields. See the doc:`Howto amoeba +` doc page for more information about the implemention +of AMOEBA and HIPPO in LAMMPS. + +Pitorsion interactions add additional potential energy contributions +to 6-tuples of atoms IJKLMN which have a bond between atoms K and L, +where both K and L are additionally bonded to exactly two other atoms. +Namely K is also bonded to I and J. And L is also bonded to M and N. + +The examples/amoeba directory has a sample input script and data file +for ubiquitin, which illustrates use of the fix amoeba/pitorsion +command. + +As in the example above, this fix should be used before reading a data +file that contains a listing of pitorsion interactions. The +*filename* specified should contain the pitorsion parameters for the +AMOEBA or HIPPO force field. + +The data file read by the "read_data" must contain the topology of all +the pitorsion interactions, similar to the topology data for bonds, +angles, dihedrals, etc. Specially it should have two lines like these +in its header section: + +.. parsed-literal:: + + M pitorsion types + N pitorsions + +where N is the number of pitorsion 5-body interactions and M is the +number of pitorsion types. It should also have two sections in the body +of the data file like these with M and N lines: + +.. parsed-literal:: + + PiTorsion Coeffs + + 1 6.85 + 2 10.2 + [...] + M 6.85 + +.. parsed-literal:: + + PiTorsions + + 1 1 8 10 12 18 20 + 2 5 18 20 22 25 27 + [...] + N 3 314 315 317 318 330 + +For PiTorsion Coeffs, the first column is an index from 1 to M to +enumerate the pitorsion types. The second column is the single +coefficient needed for each type. + +For PiTorsions, the first column is an index from 1 to N to enumerate +the pitorsion 5-atom tuples; it is ignored by LAMMPS. The second +column is the "type" of the interaction; it is an index into the +PiTorsion Coeffs. The remaining 5 columns are the atom IDs of the +atoms in the two 4-atom dihedrals that overlap to create the pitorsion +5-body interaction. + +Note that the "pitorsion types" and "pitorsions" and "PiTorsion +Coeffs" and "PiTorsions" keywords for the header and body sections +match those specified in the read_data command following the data file +name; see the :doc:`read_data ` page for more details. + +A data file containing pitorsion 6-body interactions can be generated +from using the tinker2lmp.py script in the tools/tinker directory of +the LAMMPS distribution. See the examples for its use for the +ubiquitin system in the tools/tinker/README file. + +The data file should be generated by using the +tools/tinker/tinker2lmp.py conversion script which creates a LAMMPS +data file from Tinker input files, including its PRM file which +contains the parameters necessary for computing pitorsion +interactions. See the example for the ubiquitin system in the +tools/tinker/README file. + +The potential energy associated with pitorsion interactions can be +output as described below. It can also be included in the total +potential energy of the system, as output by the :doc:`thermo_style +` command, if the :doc:`fix_modify energy ` +command is used, as in the example above. See the note below about +how to include the pitorsion energy when performing an :doc:`energy +minimization `. + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +This fix writes the list of pitorsion interactions to :doc:`binary +restart files `. See the :doc:`read_restart ` +command for info on how to re-specify a fix in an input script that +reads a restart file, so that the operation of the fix continues in an +uninterrupted fashion. + +The :doc:`fix_modify ` *energy* option is supported by +this fix to add the potential energy of the pitorsion interactions to +both the global potential energy and peratom potential energies of the +system as part of :doc:`thermodynamic output ` or output +by the :doc:`compute pe/atom ` command. The default +setting for this fix is :doc:`fix_modify energy yes `. + +The :doc:`fix_modify ` *virial* option is supported by +this fix to add the contribution due to the pitorsion interactions to +both the global pressure and per-atom stress of the system via the +:doc:`compute pressure ` and :doc:`compute +stress/atom ` commands. The former can be +accessed by :doc:`thermodynamic output `. The default +setting for this fix is :doc:`fix_modify virial yes `. + +This fix computes a global scalar which can be accessed by various +:doc:`output commands `. The scalar is the potential +energy discussed above. The scalar value calculated by this fix is +"extensive". + +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. + +The forces due to this fix are imposed during an energy minimization, +invoked by the :doc:`minimize ` command. + +The :doc:`fix_modify ` *respa* option is supported by this +fix. This allows to set at which level of the :doc:`r-RESPA +` integrator the fix is adding its forces. Default is the +outermost level. + +.. note:: + + For energy minimization, if you want the potential energy + associated with the pitorsion terms forces to be included in the + total potential energy of the system (the quantity being + minimized), you MUST not disable the :doc:`fix_modify ` + *energy* option for this fix. + +Restrictions +"""""""""""" + +To function as expected this fix command must be issued *before* a +:doc:`read_data ` command but *after* a :doc:`read_restart +` command. + +This fix can only be used if LAMMPS was built with the AMOEBA package. +See the :doc:`Build package ` page for more info. + +Related commands +"""""""""""""""" + +:doc:`fix_modify `, :doc:`read_data ` + +Default +""""""" + +none diff --git a/doc/src/fix_cmap.rst b/doc/src/fix_cmap.rst index 6475e538a2..19d8cd2d52 100644 --- a/doc/src/fix_cmap.rst +++ b/doc/src/fix_cmap.rst @@ -26,13 +26,13 @@ Examples Description """"""""""" -This command enables CMAP cross-terms to be added to simulations which -use the CHARMM force field. These are relevant for any CHARMM model -of a peptide or protein sequences that is 3 or more amino-acid -residues long; see :ref:`(Buck) ` and :ref:`(Brooks) ` -for details, including the analytic energy expressions for CMAP -interactions. The CMAP cross-terms add additional potential energy -contributions to pairs of overlapping phi-psi dihedrals of +This command enables CMAP 5-body interactions to be added to +simulations which use the CHARMM force field. These are relevant for +any CHARMM model of a peptide or protein sequences that is 3 or more +amino-acid residues long; see :ref:`(Buck) ` and :ref:`(Brooks) +` for details, including the analytic energy expressions for +CMAP interactions. The CMAP 5-body terms add additional potential +energy contributions to pairs of overlapping phi-psi dihedrals of amino-acids, which are important to properly represent their conformational behavior. @@ -47,15 +47,15 @@ lammps/potentials directory: charmm22.cmap and charmm36.cmap. The data file read by the "read_data" must contain the topology of all the CMAP interactions, similar to the topology data for bonds, angles, -dihedrals, etc. Specially it should have a line like this -in its header section: +dihedrals, etc. Specially it should have a line like this in its +header section: .. parsed-literal:: N crossterms -where N is the number of CMAP cross-terms. It should also have a section -in the body of the data file like this with N lines: +where N is the number of CMAP 5-body interactions. It should also +have a section in the body of the data file like this with N lines: .. parsed-literal:: @@ -66,28 +66,29 @@ in the body of the data file like this with N lines: [...] N 3 314 315 317 318 330 -The first column is an index from 1 to N to enumerate the CMAP terms; -it is ignored by LAMMPS. The second column is the "type" of the -interaction; it is an index into the CMAP force field file. The +The first column is an index from 1 to N to enumerate the CMAP 5-atom +tuples; it is ignored by LAMMPS. The second column is the "type" of +the interaction; it is an index into the CMAP force field file. The remaining 5 columns are the atom IDs of the atoms in the two 4-atom -dihedrals that overlap to create the CMAP 5-body interaction. Note -that the "crossterm" and "CMAP" keywords for the header and body -sections match those specified in the read_data command following the -data file name; see the :doc:`read_data ` page for -more details. +dihedrals that overlap to create the CMAP interaction. Note that the +"crossterm" and "CMAP" keywords for the header and body sections match +those specified in the read_data command following the data file name; +see the :doc:`read_data ` page for more details. -A data file containing CMAP cross-terms can be generated from a PDB -file using the charmm2lammps.pl script in the tools/ch2lmp directory -of the LAMMPS distribution. The script must be invoked with the -optional "-cmap" flag to do this; see the tools/ch2lmp/README file for -more information. +A data file containing CMAP 5-body interactions can be generated from +a PDB file using the charmm2lammps.pl script in the tools/ch2lmp +directory of the LAMMPS distribution. The script must be invoked with +the optional "-cmap" flag to do this; see the tools/ch2lmp/README file +for more information. The same conversion script also creates the +file of CMAP coefficient data which is read by this command. The potential energy associated with CMAP interactions can be output as described below. It can also be included in the total potential -energy of the system, as output by the -:doc:`thermo_style ` command, if the :doc:`fix_modify energy ` command is used, as in the example above. See -the note below about how to include the CMAP energy when performing an -:doc:`energy minimization `. +energy of the system, as output by the :doc:`thermo_style +` command, if the :doc:`fix_modify energy ` +command is used, as in the example above. See the note below about +how to include the CMAP energy when performing an :doc:`energy +minimization `. ---------- @@ -134,10 +135,11 @@ outermost level. .. note:: - If you want the potential energy associated with the CMAP terms - forces to be included in the total potential energy of the system - (the quantity being minimized), you MUST not disable the - :doc:`fix_modify ` *energy* option for this fix. + For energy minimization, if you want the potential energy + associated with the CMAP terms forces to be included in the total + potential energy of the system (the quantity being minimized), you + MUST not disable the :doc:`fix_modify ` *energy* option + for this fix. Restrictions """""""""""" diff --git a/doc/src/improper_amoeba.rst b/doc/src/improper_amoeba.rst new file mode 100644 index 0000000000..75f32abaf7 --- /dev/null +++ b/doc/src/improper_amoeba.rst @@ -0,0 +1,77 @@ +.. index:: improper_style amoeba + +improper_style harmonic command +=============================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + improper_style amoeba + +Examples +"""""""" + +.. code-block:: LAMMPS + + improper_style amoeba + improper_coeff 1 49.6 + +Description +""""""""""" + +The *amoeba* improper style uses the potential + +.. math:: + + E = K (\chi)^2 + +where :math:`\chi` is the improper angle and :math:`K` is a prefactor. +Note that the usual 1/2 factor is included in :math:`K`. + +This formula seems like a simplified version of the formula for the +:doc:`improper_style harmonic ` command with +:math:`\chi_0` = 0.0. However the computation of the angle +:math:`\chi` is done differently to match how the Tinker MD code +computes its out-of-plane improper for the AMOEBA and HIPPO force +fields. See the doc:`Howto amoeba ` doc page for more +information about the implemention of AMOEBA and HIPPO in LAMMPS. + +If the 4 atoms in an improper quadruplet (listed in the data file read +by the :doc:`read_data ` command) are ordered I,J,K,L then +atoms I,K,L form a plane and atom J is out-of-place. The angle +:math:`\chi_0` is computed as the Allinger angle which is defined as +the angle between the plane of I,K,L, and the vector from atom I to +atom J. + +The following coefficient must be defined for each improper type via +the :doc:`improper_coeff ` command as in the example +above, or in the data file or restart files read by the +:doc:`read_data ` or :doc:`read_restart ` +commands: + +* :math:`K` (energy) + +Note that the angle :math:`\chi` is computed in radians; hence +:math:`K` is effectively energy per radian\^2. + +---------- + +Restrictions +"""""""""""" + +This improper style can only be used if LAMMPS was built with the +AMOEBA package. See the :doc:`Build package ` doc page +for more info. + +Related commands +"""""""""""""""" + +:doc:`improper_coeff `, `improper_harmonic +:doc:` + +Default +""""""" + +none diff --git a/doc/src/improper_style.rst b/doc/src/improper_style.rst index ca7c31a923..3fd7eeb1e8 100644 --- a/doc/src/improper_style.rst +++ b/doc/src/improper_style.rst @@ -77,6 +77,7 @@ more of (g,i,k,o,t) to indicate which accelerated styles exist. * :doc:`zero ` - topology but no interactions * :doc:`hybrid ` - define multiple styles of improper interactions +* :doc:`amoeba ` - AMOEBA out-of-plane improper * :doc:`class2 ` - COMPASS (class 2) improper * :doc:`cossq ` - improper with a cosine squared term * :doc:`cvff ` - CVFF improper diff --git a/doc/src/pair_amoeba.rst b/doc/src/pair_amoeba.rst index 67c8e02557..c301caa332 100644 --- a/doc/src/pair_amoeba.rst +++ b/doc/src/pair_amoeba.rst @@ -33,100 +33,73 @@ Examples Additional info """"""""""""""" -doc:`Howto amoeba ` -doc:`bond_style amoeba ` -doc:`angle_style amoeba ` -doc:`dihedral_style amoeba ` -examples/amoeba -tools/amoeba -potentials/\*.prm.ameoba -potentials/\*.key.ameoba -potentials/\*.prm.hippo -potentials/\*.key.hippo +* doc:`Howto amoeba ` +* examples/amoeba +* tools/amoeba +* potentials/*.amoeba +* potentials/*.hippo Description """"""""""" -NOTE: edit this paragraph however you wish including adding or changing -citations (see bottom of this file) - The *amoeba* style computes the AMOEBA polarizeable field formulated -by Jay Ponder's group at U Washington at St Louis :ref:`(Ren) -`, :ref:`(Shi) `. The *hippo* style -computes the HIPPO polarizeable force field , an extension to AMOEBA, -formulated by Josh Rackers and collaborators in the Ponder group -:ref:`(Rackers) `. +by Jay Ponder's group at the U Washington at St Louis :ref:`(Ren) +`, :ref:`(Shi) `. The *hippo* style computes +the HIPPO polarizeable force field, an extension to AMOEBA, formulated +by Josh Rackers and collaborators in the Ponder group :ref:`(Rackers) +`. These force fields can be used when polarization effects are desired in simulations of water, organic molecules, and biomolecules including -proteins, provided that parameterizations (force field files) are -available for the systems you are interested in. Files in the LAMMPS -potentials with a "amoeba" or "hippo" suffix can be used. The Tinker -distribution and website may have other such files. +proteins, provided that parameterizations (Tinker PRM force field +files) are available for the systems you are interested in. Files in +the LAMMPS potentials directory with a "amoeba" or "hippo" suffix can +be used. The Tinker distribution and website have additional force +field files as well. -NOTE: replace these LaTeX formulas with a list of AMOEBA and HIPPO -terms in some simple notation. If desired, a more detailed -mathematical expression for each term could also be included below -the initial 2 formulas. - -The AMOEBA force field can be written as a collection of terms +As discussed on the doc:`Howto amoeba ` doc page, the +intermolecular (non-bonded) portion of the AMOEBA force field contains +these terms: .. math:: - E = 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - - \left(\frac{\sigma}{r}\right)^6 \right] - \qquad r < r_c + U_{amoeba} = U_{hal} + U_{multipole} + U_{polar} -:math:`r_c` is the cutoff, blah is the dipole, etc - -The HIPPO force field is similar but alters some of the terms +while the HIPPO force field contains these terms: .. math:: - E = 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - - \left(\frac{\sigma}{r}\right)^6 \right] - \qquad r < r_c + U_{hippo} = U_{repulsion} + U_{dispersion} + U_{multipole} + U_{polar} + U_{qxfer} -:math:`r_c` is the cutoff, blah is the dipole, etc +Conceptually, these terms compute the following interactions: + +* :math:`U_{hal}` = buffered 14-7 van der Waals with offsets applied to hydrogen atoms +* :math:`U_{repulsion}` = Pauli repulsion due to rearrangement of electron density +* :math:`U_{dispersion}` = dispersion between correlated, instantaneous induced dipole moments +* :math:`U_{multipole}` = electrostatics between permanent point charges, dipoles, and quadrupoles +* :math:`U_{polar}` = electronic polarization bewteen induced point dipoles +* :math:`U_{qxfer}` = charge transfer effects -NOTE: Add a sentence for each term to explain what physical effects -the FFs are encoding. E.g. The polar term iteratively computes an -induced dipole for each atom, then calculates dipole-dipole -interactions ... - -See the AMOEBA and HIPPO papers for further details. +Note that the AMOEBA versus HIPPO force fields typically compute the +same term differently using their own formulas. The references on +this doc page give full details for both force fields. .. note:: The AMOEBA and HIPPO force fields compute long-range charge, dipole, - and quadrupole interactions (NOTE: also long-range dispersion?). - However, unlike other models with long-range interactions in LAMMPS, - this does not require use of a KSpace style via the + and quadrupole interactions as well as long-range dispersion + effects. However, unlike other models with long-range interactions + in LAMMPS, this does not require use of a KSpace style via the :doc:`kspace_style ` command. That is because for AMOEBA and HIPPO the long-range computations are intertwined with - the pairwise compuations. So these pair style include both short- + the pairwise computations. So these pair style include both short- and long-range computations. This means the energy and virial computed by the pair style as well as the "Pair" timing reported by - LAMMPS will include the long-range components. - -.. note:: - - As explained on the :doc:`Howto amoeba ` doc page, use - of these pair styles to run a simulation with the AMOEBA or HIPPO - force fields requires your input script to use the :doc:`atom_style - hybrid full amoeba ` atom style, AMOEBA versions of - bond/angle/dihedral styles, the :doc:`special_bonds one/five - ` option, and the :doc:`fix property/atom one/five - ` command to define several additional per-atom - properties. The latter requires a "Tinker Types" section be - included in the LAMMPS data file. This can be auto-generated using - the tools/amoeba/tinker2lmp.py Python script. See the :doc:`Howto - amoeba ` doc page and tools/amoeba/README file for - details on using that tool. + LAMMPS will include the long-range calculations. The implementation of the AMOEBA and HIPPO force fields in LAMMPS was -done using code provided by the Ponder group from their `Tinker MD -code `_ written in F90. +done using F90 code provided by the Ponder group from their `Tinker MD +code `_. NOTE: what version of AMOEBA and HIPPO does LAMMPS implement? @@ -140,22 +113,22 @@ Only a single pair_coeff command is used with either the *amoeba* and pair_coeff * * ../potentials/protein.prm.amoeba ../potentials/protein.key.amoeba pair_coeff * * ../potentials/water.prm.hippo ../potentials/water.key.hippo -See examples of these files in the potentials directory. +Examples of the PRM files are in the potentials directory with an +*.amoeba or *.hippo suffix. The examples/amoeba directory has +examples of both PRM and KEY files. -The format of a PRM file is a collection of sections, each with -multiple lines. These are the sections which LAMMPS recognizes: +A Tinker PRM file composed of sections, each of which has multiple +lines. A Tinker KEY file is composed of lines, each of which has a keyword, +which can be followed by zero or more parameters. -NOTE: need to list these, possibly there are some LAMMPS skips - or which need to be removed for LAMMPS to use the PRM file? +The list of PRM sections and KEY keywords which LAMMPS recognizes are +listed on the doc:`Howto amoeba ` doc page. If not +recognized, the section or keyword is skipped. -The format of a KEY file is a series of lines, with one parameter and -its value per line. These are the parameters which LAMMPS recognizes: - -NOTE: need to list these, possibly there are some keywords LAMMPS skips - or which need to be removed for LAMMPS to use the KEY file? - -NOTE: Any other info about PRM and KEY files we should explain -for LAMMPS users here? +Note that if the KEY file is specified as NULL, then no file is +required; default values for various AMOEBA/HIPPO settings are used. +The doc:`Howto amoeba ` doc page also gives the default +settings. ---------- @@ -184,14 +157,38 @@ enabled if LAMMPS was built with that package. See the :doc:`Build package ` doc page for more info. The AMOEBA and HIPPO potential (PRM) and KEY files provided with -LAMMPS in the potentials directory are Tinker files parameterized for -Tinker units. Their numeric parameters are converted by LAMMPS to its -real units :doc:`units `. You can only use these pair styles -with real units. +LAMMPS in the potentials and examples/amoeva directories are Tinker +files parameterized for Tinker units. Their numeric parameters are +converted by LAMMPS to its real units :doc:`units `. Thus uou +can only use these pair styles with real units. These potentials do not yet calculate per-atom energy or virial contributions. +As explained on the :doc:`Howto amoeba ` doc page, use +of these pair styles to run a simulation with the AMOEBA or HIPPO +force fields requires several things. + +The first is a data file generated by the tools/tinker/tinker2lmp.py +conversion script which uses Tinker file force field file input to +create a data file compatible with LAMMPS. + +The second is use of these commands: + +* :doc:`atom_style amoeba ` +* :doc:`fix property/atom ` +* :doc:`special_bonds one/five ` + +And third, depending on the model being simulated, these +commands for intramolecular interactions may also be required: + +* :doc:`bond_style class2 ` +* :doc:`angle_style amoeba ` +* :doc:`dihedral_style fourier ` +* :doc:`improper_style amoeba ` +* :doc:`fix amoeba/pitorsion ` +* :doc:`fix amoeba/bitorsion ` + ---------- Related commands diff --git a/doc/src/special_bonds.rst b/doc/src/special_bonds.rst index 4b197fbf3a..8cb81c4cdb 100644 --- a/doc/src/special_bonds.rst +++ b/doc/src/special_bonds.rst @@ -11,7 +11,7 @@ Syntax special_bonds keyword values ... * one or more keyword/value pairs may be appended -* keyword = *amber* or *charmm* or *dreiding* or *fene* or *lj/coul* or *lj* or *coul* or *angle* or *dihedral* +* keyword = *amber* or *charmm* or *dreiding* or *fene* or *lj/coul* or *lj* or *coul* or *angle* or *dihedral* or *one/five* .. parsed-literal:: @@ -27,6 +27,7 @@ Syntax w1,w2,w3 = weights (0.0 to 1.0) on pairwise Coulombic interactions *angle* value = *yes* or *no* *dihedral* value = *yes* or *no* + *one/five* value = *yes* or *no* Examples """""""" @@ -45,10 +46,10 @@ Description Set weighting coefficients for pairwise energy and force contributions between pairs of atoms that are also permanently bonded to each other, either directly or via one or two intermediate bonds. These weighting -factors are used by nearly all :doc:`pair styles ` in LAMMPS -that compute simple pairwise interactions. Permanent bonds between -atoms are specified by defining the bond topology in the data file -read by the :doc:`read_data ` command. Typically a +factors are used by nearly all :doc:`pair styles ` in +LAMMPS that compute simple pairwise interactions. Permanent bonds +between atoms are specified by defining the bond topology in the data +file read by the :doc:`read_data ` command. Typically a :doc:`bond_style ` command is also used to define a bond potential. The rationale for using these weighting factors is that the interaction between a pair of bonded atoms is all (or mostly) @@ -58,31 +59,34 @@ atoms should be excluded (or reduced by a weighting factor). .. note:: - These weighting factors are NOT used by :doc:`pair styles ` that compute many-body interactions, since the - "bonds" that result from such interactions are not permanent, but are - created and broken dynamically as atom conformations change. Examples - of pair styles in this category are EAM, MEAM, Stillinger-Weber, - Tersoff, COMB, AIREBO, and ReaxFF. In fact, it generally makes no - sense to define permanent bonds between atoms that interact via these - potentials, though such bonds may exist elsewhere in your system, - e.g. when using the :doc:`pair_style hybrid ` command. - Thus LAMMPS ignores special_bonds settings when many-body potentials - are calculated. Please note, that the existence of explicit bonds - for atoms that are described by a many-body potential will alter the - neighbor list and thus can render the computation of those interactions - invalid, since those pairs are not only used to determine direct - pairwise interactions but also neighbors of neighbors and more. - The recommended course of action is to remove such bonds, or - if - that is not possible - use a special bonds setting of 1.0 1.0 1.0. + These weighting factors are NOT used by :doc:`pair styles + ` that compute many-body interactions, since the + "bonds" that result from such interactions are not permanent, but + are created and broken dynamically as atom conformations change. + Examples of pair styles in this category are EAM, MEAM, + Stillinger-Weber, Tersoff, COMB, AIREBO, and ReaxFF. In fact, it + generally makes no sense to define permanent bonds between atoms + that interact via these potentials, though such bonds may exist + elsewhere in your system, e.g. when using the :doc:`pair_style + hybrid ` command. Thus LAMMPS ignores special_bonds + settings when many-body potentials are calculated. Please note, + that the existence of explicit bonds for atoms that are described + by a many-body potential will alter the neighbor list and thus can + render the computation of those interactions invalid, since those + pairs are not only used to determine direct pairwise interactions + but also neighbors of neighbors and more. The recommended course + of action is to remove such bonds, or - if that is not possible - + use a special bonds setting of 1.0 1.0 1.0. .. note:: Unlike some commands in LAMMPS, you cannot use this command multiple times in an incremental fashion: e.g. to first set the LJ - settings and then the Coulombic ones. Each time you use this command - it sets all the coefficients to default values and only overrides the - one you specify, so you should set all the options you need each time - you use it. See more details at the bottom of this page. + settings and then the Coulombic ones. Each time you use this + command it sets all the coefficients to default values and only + overrides the one you specify, so you should set all the options + you need each time you use it. See more details at the bottom of + this page. The Coulomb factors are applied to any Coulomb (charge interaction) term that the potential calculates. The LJ factors are applied to the @@ -94,14 +98,14 @@ exclude it completely. The first of the 3 coefficients (LJ or Coulombic) is the weighting factor on 1-2 atom pairs, which are pairs of atoms directly bonded to -each other. The second coefficient is the weighting factor on 1-3 atom -pairs which are those separated by 2 bonds (e.g. the two H atoms in a -water molecule). The third coefficient is the weighting factor on 1-4 -atom pairs which are those separated by 3 bonds (e.g. the first and fourth -atoms in a dihedral interaction). Thus if the 1-2 coefficient is set -to 0.0, then the pairwise interaction is effectively turned off for -all pairs of atoms bonded to each other. If it is set to 1.0, then -that interaction will be at full strength. +each other. The second coefficient is the weighting factor on 1-3 +atom pairs which are those separated by 2 bonds (e.g. the two H atoms +in a water molecule). The third coefficient is the weighting factor +on 1-4 atom pairs which are those separated by 3 bonds (e.g. the first +and fourth atoms in a dihedral interaction). Thus if the 1-2 +coefficient is set to 0.0, then the pairwise interaction is +effectively turned off for all pairs of atoms bonded to each other. +If it is set to 1.0, then that interaction will be at full strength. .. note:: @@ -184,21 +188,27 @@ interaction between atoms 2 and 5 will be unaffected (full weighting of 1.0). If the *dihedral* keyword is specified as *no* which is the default, then the 2,5 interaction will also be weighted by 0.5. +The *one/five* keyword enable calculation and storage of a list of 1-5 +neighbors in the molecular topology for each atom. It is required by +some pair styles, such as :doc:`pair_style amoeba ` and +:doc:`pair_style hippo `. + ---------- .. note:: LAMMPS stores and maintains a data structure with a list of the - first, second, and third neighbors of each atom (within the bond topology of - the system). If new bonds are created (or molecules added containing - atoms with more special neighbors), the size of this list needs to - grow. Note that adding a single bond always adds a new first neighbor - but may also induce \*many\* new second and third neighbors, depending on the - molecular topology of your system. Using the *extra/special/per/atom* - keyword to either :doc:`read_data ` or :doc:`create_box ` - reserves empty space in the list for this N additional first, second, or third - neighbors to be added. If you do not do this, you may get an error - when bonds (or molecules) are added. + first, second, and third neighbors of each atom (within the bond + topology of the system). If new bonds are created (or molecules + added containing atoms with more special neighbors), the size of + this list needs to grow. Note that adding a single bond always + adds a new first neighbor but may also induce \*many\* new second + and third neighbors, depending on the molecular topology of your + system. Using the *extra/special/per/atom* keyword to either + :doc:`read_data ` or :doc:`create_box ` + reserves empty space in the list for this N additional first, + second, or third neighbors to be added. If you do not do this, you + may get an error when bonds (or molecules) are added. ---------- @@ -226,16 +236,16 @@ In the first case you end up with (after the second command): LJ: 0.0 0.0 0.0 Coul: 0.0 0.0 1.0 -while only in the second case, you get the desired settings of: +while only in the second case do you get the desired settings of: .. parsed-literal:: LJ: 0.0 1.0 1.0 Coul: 0.0 0.0 1.0 -This happens because the LJ (and Coul) settings are reset to -their default values before modifying them, each time the -*special_bonds* command is issued. +This happens because the LJ (and Coul) settings are reset to their +default values before modifying them, each time the *special_bonds* +command is issued. Restrictions """""""""""" diff --git a/src/AMOEBA/amoeba_file.cpp b/src/AMOEBA/amoeba_file.cpp index 2469c08f43..f67eacd879 100644 --- a/src/AMOEBA/amoeba_file.cpp +++ b/src/AMOEBA/amoeba_file.cpp @@ -493,9 +493,6 @@ void PairAmoeba::read_keyfile(char *filename) } else if (strcmp(keyword,"pcg-peek") == 0) { if (nwords != 2) error->all(FLERR,"AMOEBA keyfile line is invalid"); pcgpeek = utils::numeric(FLERR,words[1],true,lmp); - } else if (strcmp(keyword,"usolve-diag") == 0) { - if (nwords != 2) error->all(FLERR,"AMOEBA keyfile line is invalid"); - udiag = utils::numeric(FLERR,words[1],true,lmp); } else {} diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp index 4288de6e36..ebf4a19bc6 100644 --- a/src/AMOEBA/pair_amoeba.cpp +++ b/src/AMOEBA/pair_amoeba.cpp @@ -620,14 +620,6 @@ void PairAmoeba::coeff(int narg, char **arg) if (narg == 3) read_keyfile(NULL); else read_keyfile(arg[3]); - // required for now, until 1-5 settings are implemented in LAMMPS - - //if (special_hal[4] != 1.0 || special_repel[4] != 1.0 || - // special_disp[4] != 1.0 || special_mpole[4] != 1.0 || - // special_polar_pscale[4] != 1.0 || special_polar_piscale[4] != 1.0 || - // special_polar_wscale[4] != 1.0) - // error->all(FLERR,"AMOEBA 1-5 weights must be 1.0 for now"); - // compute Vdwl mixing rules, only for AMOEBA if (amoeba) {