This file outlines the force-field formulas used in LAMMPS. Read this file in conjunction with the data_format and units files.

The sections of this page are as follows:

- Nonbond Coulomb
- Nonbond van der Waals
- Mixing Rules for Lennard-Jones
- Bonds
- Angles
- Dihedrals
- Impropers
- Class 2 Force Field

E = C q1 q2 / (epsilon * r) r = distance (computed by LAMMPS) C = hardwired constant to convert to energy units q1,q2 = charge of each atom in electron units (proton = +1), specified in "Atoms" entry in data file epsilon = dielectric constant (vacuum = 1.0), set by user in input command file

E = Do [ (Ro/r)^12 - 2 (Ro/r)^6 ] standard Lennard Jones potential r = distance (computed by LAMMPS) coeff1 = Do (energy) coeff2 = Ro (distance) 2 coeffs are listed in data file or set in input script 1 cutoff is set in input script

E = Do [ (Ro/r)^12 - 2(Ro/r)^6 ] for r < r_inner = spline fit for r_inner < r < cutoff = 0 for r > cutoff switching function (spline fit) is applied to standard LJ within a switching region (from r_inner to cutoff) so that energy and force go smoothly to zero spline coefficients are computed by LAMMPS so that at inner cutoff (r_inner) the potential, force, and 1st-derivative of force are all continuous, and at outer cutoff (cutoff) the potential and force both go to zero r = distance (computed by LAMMPS) coeff1 = Do (energy) coeff2 = Ro (distance) 2 coeffs are listed in data file or set in input script 2 cutoffs (r_inner and cutoff) are set in input script

E = Do [ (Ro/(r - delta))^12 - 2(Ro/(r - delta))^6 ] same as lj/cutoff except that r is shifted by delta r = distance (computed by LAMMPS) coeff1 = Do (energy) coeff2 = Ro (distance) coeff3 = delta (distance) 3 coeffs are listed in data file or set in input script 1 cutoff is set in input script

E = A * [ 1 + cos( pi * r / cutoff ) ] useful for pushing apart overlapping atoms by ramping A over time r = distance (computed by LAMMPS) coeff1 = prefactor A at start of run (energy) coeff2 = prefactor A at end of run (energy) 2 coeffs are listed in data file or set in input script 1 cutoff is set in input script

E = epsilon [ 2 (sigma/r)^9 - 3 (sigma/r)^6 ] used with class2 bonded force field r = distance (computed by LAMMPS) coeff1 = epsilon (energy) coeff2 = sigma (distance) 2 coeffs are listed in data file or set in input script 1 cutoff is set in input script

E = Do [ (Ro/r)^12 - 2(Ro/r)^6 ] for r < r_inner = switch * E for r_inner < r < cutoff = 0 for r > cutoff where switch = [(cutoff^2 - r^2)^2 * (cutoff^2 + 2*r^2 - 3*r_inner)] / [(cutoff^2 - r_inner^2)^3] switching function is applied to standard LJ within a switching region (from r_inner to cutoff) so that energy and force go smoothly to zero switching function causes that at inner cutoff (r_inner) the potential and force are continuous, and at outer cutoff (cutoff) the potential and force both go to zero r = distance (computed by LAMMPS) coeff1 = Do (energy) coeff2 = Ro (distance) coeff3 = Do for 1-4 interactions (energy) coeff4 = Ro for 1-4 interactions (distance) 4 coeffs are listed in data file or set in input script 2 cutoffs (r_inner and cutoff) are set in input script

E = Do [(6/6-y) exp (y(1-r/Ro) - (y/y-6)(Ro/r)^6 ] for r < cutoff = 0 for r > cutoff Input parameters are Do, Ro, and y. The internal representation of expo_6 is

E = A*exp(-CR) - B/r^6

with combination rules calculated with A,B,C as follows:

Aij = sqrt(Aii Ajj) Bij = sqrt(Bii Bjj) Cij = 1/2(Cii + Cjj)

E = Do [(6/6-y) exp (y(1-r/Ro) - (y/y-6)(Ro/r)^6 ] for r < r_inner = switch * E for r_inner < r < cutoff = 0 for r > cutoff where switch = [(r_outer^2 - r^2)^2 * (r_outer + 2*r - 3*r_inner)] / [(r_outer^2 - r_outer^2)^3] switching function ensures continuity of energy and forces at r_inner. Forces are not monotonically vanishing to zero. A ratio r_outer/r_inner >= 5/4 is recommended.

E = Do [(6/6-y) exp (y(1-r/Ro) - (y/y-6)(Ro/r)^6 ] for r < cutoff = switch for cutoff < r < r_outer = 0 for r > r_outer where switch = a*(r-cutoff)^3 + d1*(r-cutoff) + e1

where

a = -4 d1^3/27 e1^2

d1 = dE/dr at r=cutoff

e1 = E(r=cutoff)

r_outer = cutoff - ( 3e1/2d1 ) switching function ensures continuity of energy and forces at r=cutoff. Forces and energy monotonically vanish to zero at r_outer. The ratio r_outer/r_inner is approximately equal to 5/4. The info flag 1 prints out the values of the coefficients and cutoffs.

Aij = sqrt(Aii * Ajj)

Bij = sqrt(Bii * Bjj)

Cij = 1/2(Cii * Cjj)

where the exponential-6 E = Do [(6/6-y) exp (y(1-r/Ro) - (y/y-6)(Ro/r)^6 ] is internally handled as

E(r)=A*exp(-Cr) - B/r^6

geometric: Do_ij = sqrt(Do_i * Do_j) Ro_ij = sqrt(Ro_i * Ro_j) arithmetic: Do_ij = sqrt(Do_i * Do_j) Ro_ij = (Ro_i + Ro_j) / 2 sixthpower: Do_ij = (2 * sqrt(Do_i*Do_j) * Ro_i^3 * Ro_j^3) / (Ro_i^6 + Ro_j^6) Ro_ij= ((Ro_i**6 + Ro_j**6) / 2) ^ (1/6)The default mixing rule for nonbond styles lj/cutoff, lj/switch, lj/shift, and soft is "geometric". The default for nonbond style class2/cutoff is "sixthpower".

The default can be overridden using the "mixing style" command. Three exceptions to this are (1) the exponential-6 form, (2) for the nonbond style soft, for which only an Do prefactor is input. This is always mixed geometrically. Also, for nonbond style lj/shift, the delta coefficient is always mixed using the rule

- delta_ij = (delta_i + delta_j) / 2

E = K (r - r0)^2 standard harmonic spring r = distance (computed by LAMMPS) coeff1 = K (energy/distance^2) (the usual 1/2 is included in the K) coeff2 = r0 (distance) 2 coeffs are listed in data file or set in input script

E = -0.5 K R0^2 * ln[1 - (r/R0)^2] + 4 epsilon [(sigma/r)^12 - (sigma/r)^6] + epsilon finite extensible nonlinear elastic (FENE) potential for polymer bead-spring models see Kremer, Grest, J Chem Phys, 92, p 5057 (1990) r = distance (computed by LAMMPS) coeff1 = K (energy/distance^2) coeff2 = R0 (distance) coeff3 = epsilon (energy) coeff4 = sigma (distance) 1st term is attraction, 2nd term is repulsion (shifted LJ) 1st term extends to R0 2nd term only extends to the minimum of the LJ potential, a cutoff distance computed by LAMMPS (2^(1/6) * sigma) 4 coeffs are listed in data file or set in input script

E = -0.5 K R0^2 * ln[1 - ((r - delta)/R0)^2] + 4 epsilon [(sigma/(r - delta))^12 - (sigma/(r - delta))^6] + epsilon same as FENE/standard expect that r is shifted by delta r = distance (computed by LAMMPS) coeff1 = K (energy/distance^2) coeff2 = R0 (distance) coeff3 = epsilon (energy) coeff4 = sigma (distance) coeff5 = delta (distance) 1st term is attraction, 2nd term is repulsion (shifted LJ) 1st term extends to R0 2nd term only extends to the minimum of the LJ potential, a cutoff distance computed by LAMMPS (2^(1/6) * sigma + delta) 5 coeffs are listed in data file or set in input script

E = epsilon (r - r0)^2 / [ lamda^2 - (r - r0)^2 ] non-harmonic spring of equilibrium length r0 with finite extension of lamda see Rector, Van Swol, Henderson, Molecular Physics, 82, p 1009 (1994) r = distance (computed by LAMMPS) coeff1 = epsilon (energy) coeff2 = r0 (distance) coeff3 = lamda (distance) 3 coeffs are listed in data file or set in input script

E = K2 (r - r0)^2 + K3 (r - r0)^3 + K4 (r - r0)^4 r = distance (computed by LAMMPS) coeff1 = r0 (distance) coeff2 = K2 (energy/distance^2) coeff3 = K3 (energy/distance^3) coeff4 = K4 (energy/distance^4) 4 coeffs are listed in data file - cannot be set in input script

E = 1/2 K (theta - theta0)^2 theta = radians (computed by LAMMPS) coeff1 = K (energy/radian^2) (the usual 1/2 is NOT included in the K) coeff2 = theta0 (degrees) (converted to radians within LAMMPS) 2 coeffs are listed in data file or set in input script

E = K2 (theta - theta0)^2 + K3 (theta - theta0)^3 + K4 (theta - theta0)^4 theta = radians (computed by LAMMPS) coeff1 = theta0 (degrees) (converted to radians within LAMMPS) coeff2 = K2 (energy/radian^2) coeff3 = K3 (energy/radian^3) coeff4 = K4 (energy/radian^4) 4 coeffs are listed in data file - cannot be set in input script

(harmonic + Urey-Bradley) E = K (theta - theta0)^2 + K_UB (r_13 - r_UB)^2 theta = radians (computed by LAMMPS) r_13 = distance (computed by LAMMPS) coeff1 = K (energy/radian^2) (the usual 1/2 is included in the K) coeff2 = theta0 (degrees) (converted to radians within LAMMPS) coeff3 = K_UB (energy/distance^2) coeff4 = r_UB (distance) 4 coeffs are listed in data file or set in input script

E = K [1 + d * cos (n*phi) ] phi = radians (computed by LAMMPS) coeff1 = K (energy) coeff2 = d (+1 or -1) coeff3 = n (1,2,3,4,6) Cautions when comparing to other force fields: some force fields reverse the sign convention on d so that E = K [1 - d * cos(n*phi)] some force fields divide/multiply K by the number of multiple torsions that contain the j-k bond in an i-j-k-l torsion some force fields let n be positive or negative which corresponds to d = 1,-1 in the LAMMPS force field, the trans position = 180 degrees, while in some force fields trans = 0 degrees 3 coeffs are listed in data file or set in input script

E = SUM(n=1,3) { K_n [ 1 - cos( n*Phi - Phi0_n ) ] } phi = radians (computed by LAMMPS) coeff1 = K_1 (energy) coeff2 = Phi0_1 (degrees) (converted to radians within LAMMPS) coeff3 = K_2 (energy) coeff4 = Phi0_2 (degrees) (converted to radians within LAMMPS) coeff5 = K_3 (energy) coeff6 = Phi0_3 (degrees) (converted to radians within LAMMPS) 6 coeffs are listed in data file - cannot be set in input script

E = SUM(n=1,5) { A_n * cos(Phi)^n } phi = radians (computed by LAMMPS) coeff1 = A_1 coeff2 = A_2 coeff3 = A_3 coeff4 = A_4

coeff5 = A_5 5 coeffs are listed in data file or set in input script

(harmonic + 1-4 interactions) E = K [1 + cos (n*phi + d) ] phi = radians (computed by LAMMPS) coeff1 = K (energy) coeff2 = n (1,2,3,4,6) coeff3 = d (0 or 180 degrees) (converted to radians within LAMMPS) coeff4 = weighting factor to turn on/off 1-4 neighbor nonbond interactions coeff4 weight values are from 0.0 to 1.0 and are used to multiply the energy and force interaction (both Coulombic and LJ) between the 2 atoms weight of 0.0 means no interaction weight of 1.0 means full interaction must be used with the special bonds charmm command (special bonds 0 0 0) which shuts off the uniform special bonds and allows pair-specific special bonds for the 1-4 interactions to be defined in the data file LAMMPS assumes that all 1-4 interaction distances, which are generally less than 6 angstroms, are less than the smallest of the inner LJ and Coulombic cutoffs, which are generally at least 8 angstroms. 4 coeffs are listed in data file or set in input script

E = SUM 1/2 Kn [1 - d * cos (n*phi) ] n phi = radians (computed by LAMMPS), 0 for the cis conformation, 180 trans phi is positive angle clockwise when looking form J toward K. coeff1 = Kn (energy) coeff2 = d (+1 or -1) coeff3 = n (1,2,3,4,6) Cautions when comparing to other force fields: notice the 1/2 in front of the input force constant (data file) as well as the negative sign in front of d. The value of Kn in dreiding is divided by the number of multiple torsions that contain the j-k bond in an i-j-k-l torsion. Kn represents 1/2 the rotational barrier, n the periodicity. d=+1, the cis conformation is the minimum, while for d=-1 the cis conformation is a maximum. In addition the value of K is rescaled for exocyclic "aromatic rings" according to the force field used in the translation of the file that creates the data file. All fourier terms (up to six) can be specified although one is sufficient often. Use the same type for all terms in consecutive lines in the data or input files. 3 coeffs are listed in data file or set in input script for each term in each dihedral type dreiding.

E = K (chi - chi0)^2 chi = radians (computed by LAMMPS) coeff1 = K (energy/radian^2) (the usual 1/2 is included in the K) coeff2 = chi0 (degrees) (converted to radians within LAMMPS) 2 coeffs are listed in data file or set in input script

E = K [1 + d * cos (n*chi) ] chi = radians (computed by LAMMPS) coeff1 = K (energy) coeff2 = d (+1 or -1) coeff3 = n (0,1,2,3,4,6) 3 coeffs are listed in data file or set in input script

same formula, coeffs, and meaning as "harmonic" except that LAMMPS averages all 3 angle-contributions to chi in class 2 this is called a Wilson out-of-plane interaction 2 coeffs are listed in data file - cannot be set in input script

E = K (r - r0) * (r' - r0') r,r' = distance (computed by LAMMPS) coeff1 = K (energy/distance^2) coeff2 = r0 (distance) coeff3 = r0' (distance) 3 coeffs are input in data file

E = K_n (r - r0_n) * (theta - theta0) r = distance (computed by LAMMPS) theta = radians (computed by LAMMPS) coeff1 = K_1 (energy/distance-radians) coeff2 = K_2 (energy/distance-radians) coeff3 = r0_1 (distance) coeff4 = r0_2 (distance) Note: theta0 is known from angle coeffs so don't need it specified here 4 coeffs are listed in data file

E = (r - r0) * [ F1*cos(phi) + F2*cos(2*phi) + F3*cos(3*phi) ] r = distance (computed by LAMMPS) phi = radians (computed by LAMMPS) coeff1 = F1 (energy/distance) coeff2 = F2 (energy/distance) coeff3 = F3 (energy/distance) coeff4 = r0 (distance) 4 coeffs are listed in data file

E = (r - r0_n) * [ F1_n*cos(phi) + F2_n*cos(2*phi) + F3_n*cos(3*phi) ] r = distance (computed by LAMMPS) phi = radians (computed by LAMMPS) coeff1 = F1_1 (energy/distance) coeff2 = F2_1 (energy/distance) coeff3 = F3_1 (energy/distance) coeff4 = F1_2 (energy/distance) coeff5 = F2_3 (energy/distance) coeff6 = F3_3 (energy/distance) coeff7 = r0_1 (distance) coeff8 = r0_2 (distance) 8 coeffs are listed in data file

E = (theta - theta0) * [ F1_n*cos(phi) + F2_n*cos(2*phi) + F3_n*cos(3*phi) ] theta = radians (computed by LAMMPS) phi = radians (computed by LAMMPS) coeff1 = F1_1 (energy/radians) coeff2 = F2_1 (energy/radians) coeff3 = F3_1 (energy/radians) coeff4 = F1_2 (energy/radians) coeff5 = F2_3 (energy/radians) coeff6 = F3_3 (energy/radians) coeff7 = theta0_1 (degrees) (converted to radians within LAMMPS) coeff8 = theta0_2 (degrees) (converted to radians within LAMMPS) 8 coeffs are listed in data file

E = K (theta - theta0) * (theta' - theta0') * (phi - phi0) theta,theta' = radians (computed by LAMMPS) phi = radians (computed by LAMMPS) coeff1 = K (energy/radians^3) coeff2 = theta0 (degrees) (converted to radians within LAMMPS) coeff3 = theta0' (degrees) (converted to radians within LAMMPS) Note: phi0 is known from dihedral coeffs so don't need it specified here 3 coeffs are listed in data file

E = K * (r1 - r10)*(r3 - r30) r1,r3 = bond lengths of bonds 1 and 3 (computed by LAMMPS) coeff1 = K (energy/distance^2) coeff2 = r10 (distance) = equilibrium bond length for bond 1 coeff3 = r30 (distance) = equilibrium bond length for bond 3 K is only non-zero for aromatic rings 3 coeffs are listed in data file

E = K_n (theta - theta0_n) * (theta' - theta0_n') theta,theta' = radians (computed by LAMMPS) coeff1 = K_1 (energy/radians^2) coeff2 = K_2 (energy/radians^2) coeff3 = K_3 (energy/radians^2) coeff4 = theta0_1 (degrees) (converted to radians within LAMMPS) coeff5 = theta0_2 (degrees) (converted to radians within LAMMPS) coeff6 = theta0_3 (degrees) (converted to radians within LAMMPS) 6 coeffs are listed in data file