Forcefield development

(updated )

Comments

Torsion conventions

  • Use MATLAB/Python to plot the torsional potential as presented in the paper, using the listed parameters and the given functional form.
  • Check the torsion energy against physical intuition to ascertain the dihedral angle convention (i.e., is trans 180 or 0 degree?). Typically the OPLS functional form uses the protein convention while the Ryckaert-Bellemans functional form uses the polymer convention, but papers may not follow these conventions, which is why we need to perform this step.
  • Refer to the GROMACS/LAMMPS/Topmon documentation to make sure the same convention is used. We may need to modify the torsion parameters to account for the different conventions.
  • For fitting new torsions, it is probably good to be consistent with a recent example, which uses torsion type 6:
  1. Cosine power series (generalized OPLS form)
  2. One phase shift
  3. Protein convention, where cis == 0 degrees
  • Note that torsion energy profile must have a period of 2π and satisfy: U(φ) = U(-φ) and U(φ) = U(2π - φ). The latter two identities are easy to see from the cosine series we use to fit them, but we can also think about it physically. Picture a C-C-C-C dihedral torsion at 0° and rotate it to 30° and -30° — we can easily see that the two conformations are mirror images of each other and must therefore have the same energy. A torsion scan from -20° to 160° is therefore incomplete, because [-20°, 0] and [0, 20°] contain identical information. An easy way to avoid this issue is to always start with the scanned dihedral at either 0 or 180°, and go in one direction to cover a 180° range.

Intramolecular 1-4 interactions

  • Different forcefields treat 1-4 interactions (non-bonded interactions for atoms within the same molecule that are separated by three bonds) differently in their parametrization. Early OPLS-UA excludes 1-4 Lennard-Jones (LJ) interactions while OPLS-AA includes both 1-4 LJ and Coulomb interactions but scales them down by a factor of 0.5. Early AMBER also uses 1/2 LJ + 1/2 Coulomb 1-4 interactions but more recent AMBER uses 1/2 LJ + 5/6 Coulomb 1-4 interactions. CHARMM does not include 1-4 interactions at all. This convention along with the torsion potentials together define the forcefield.
  • It was found that to use AMBER torsion potentials with LJ parameters from OPLS-UA, it was necessary to scale the LJ -14 interactions by 1/8 instead of 1/2, suggesting that the LJ parameters in OPLS-UA are more repulsive than in AMBER. Since early TraPPE papers often borrow torsion potentials from OPLS and Amber and LJ parameters in TraPPE are closer to OPLS-UA than AMBER, most TraPPE papers uses zero LJ + 1/2 Coulomb 1-4 interactions. See below for notes on a messy history.
  • For partially rigid molecules such as substituted aromatics, OPLS-UA and TraPPE exclude both LJ and Coulomb interactions between the atoms at the β position (separated from aromatic rings by two bonds) and all rigid atoms. The rigid atoms include the aromatic rings and the atoms at the α position. The result is that, for an isolated molecule, the conformation of the β atom is completely controlled by the bond length (α-β), angle (aromatic-α-β), and torsion (aromatic2-aromatic-α-β).
  • However, these are rules of thumb. It is critical to always double check the original references where forcefield parameters are taken.
  • In topmon, including a dihedral torsion will by default turn off 1-4 interactions. Setting lq14scale=T and qscale=0.5d0 will add back Coulomb interactions with the specified scale factor. INTRAMOLECULAR_SPECIAL can be used to add more exotic exceptions to these rules. Note that 1-4 pairs for which dihedral torsions are not defined will have both LJ and Coulomb interactions. As a result, when fitting torsion potentials, it is necessary to use a dummy torsion potential (e.g., with all parameters set to zero).
  • Topmon also automatically excludes the interactions of the β atoms connected to the first growpoint. If the molecule has multiple side chains, β atoms from the other side chains need to be excluded in INTRAMOLECULAR_SPECIAL.
  • In GROMACS, 1) setting nrexcl=3 will remove all NORMAL 1-2, 1-3, 1-4 nonbonded interactions; 2) including all 1-4 pairs under [ pairs ] for each molecule’s .itp file will add back 1-4 interactions; and 3) setting gen-pairs=yes, fudgeLJ=0, fudgeQQ=0.5 will set the correct interaction parameters. The [ exclusions ] section can be used to exclude the interactions of β atoms in a partially rigid molecule, see the documentation.
  • In LAMMPS, special_bonds lj 0.0 0.0 0.0 coul 0.0 0.0 0.5 angle yes dihedral yes will exclude LJ 1-2, 1-3, 1-4 and Coulomb 1-2, 1-3 interactions, and scale 1-4 Coulomb by 0.5. To exclude the interactions of β atoms in a partially rigid molecule, a hack is to add bonds between the β atoms and the rigid atoms with a zero force constant. The angle yes dihedral yes keywords is important with this hack. They specify that the 1-3 and 1-4 interactions are scaled or excluded only if the corresponding angle or dihedral is defined. With angle no dihedral no, the fictitious bonds can lead to extra 1-3 and 1-4 neighbors that would be incorrectly excluded.
  • Always check the actual exclusion tables as seen by the program when setting up a new calculation. For topmon, these are listed under LJ INCLUSION TABLE, CHARGE INCLUSION TABLE, 1-4 LJ SCALING FACTORS, and 1-4 CHARGE SCALING FACTORS in the output file. For GROMACS, check excls:, LJ-14:, and Coulomb-14: in the output of gmx_mpi_d dump -s NAME.tpr|less. For LAMMPS, check the output log and/or screen file for the Finding 1-2 1-3 1-4 neighbors section where the special bonds scaling factors should be printed.
  • More notes on 1-4 conventions in the TraPPE forcefield: 1-4 LJ interactions are almost never included, but TraPPE hasn’t been very consistent with 1-4 Coulomb interactions. Earlier TraPPE papers dealt with systems where at least one of the two atoms in the 1-4 pairs has zero partial charge, so there are no 1-4 Coulomb interactions to speak of. These papers then simply said that 1-4 nonbonded interactions are not included. However, since TraPPE borrowed many torsions from early OPLS-UA/AMBER, zero LJ + 1/2 Coulomb 1-4 interactions should therefore be used; see above and, for example, this paper, where 1-4 Coulomb is an important consideration and explicitly discussed. In other cases, including the TraPPE website, I suspect they simply quoted the statements in earlier TraPPE papers. When we fit missing parameters for a new molecule, even if we exclude 1-4 Coulomb (incorrectly), we MIGHT still be able to correct for this missing energy in the newly fitted potential. This is certainly not ideal, but may be what has happened over the years. In the MCCCS code, some of the input files for the test suites have 1-4 Coulomb turned on, while (fewer) others have it off.

Calculation of CM5 charges (Gaussian/cm5pac)

Stretching/bending/torsion scans

Comments

Comments powered by Disqus