Amber masthead
Filler image AmberTools24 Amber24 Manuals Tutorials Force Fields Contacts History
Filler image

Useful links:

Amber Home
Download Amber
Installation
Amber Citations
GPU Support
Updates
Mailing Lists
For Educators
File Formats
Contributors
Workshops
Developing Nonstandard Parameters

Placing Customized Extra Points in an Amber Topology File

The mdgx &parmedit module serves a purpose with overlap in parmed and tleap, and while those programs could, in time, take over this functionality it was expedient to have the feature in mdgx as there is no file format precedent for it. Extra points, or "virtual sites" are a powerful way to improve the monopole distribution of a molecule. The mdgx program has the ability to fit partial charges for these extra points once they are part of the topology. Prior to AmberTools21, the extra points had to be added as "post-translational modifications" to the topology, only mdgx could read them, and mdgx could only add them at runtime (no topology, and hence no coordinate set containing the new points, could be produced). With the addition of topology writing capabilities and a new format for encoding the extra points, these features of the system can start to become more permanent. The current master branch of pmemd and pmemd.cuda support the upgraded topology format, and is available in Amber22. Users can experiment with adding extra points, fitting partial charges or bonded parameters using the associated mdgx modules, and run short simulations with the new models using the mdgx CPU PME dynamics capability.

Topologies are still made with tleap, but can now be modified using a new file format that makes use of the mdgx &rule namelist. Like all other mdgx namelists, the list of keywords for &rule and &parmedit can be obtained by running:

mdgx -RULE
mdgx -PARMEDIT

on the command line. The &parmedit module merely takes directives for the names of the new topology and coordiante files along with an optional topology title and requests to check the total charge and connectivity. The &rule namelist is much more detailed, allowing users to set all non-bonded aspects of existing as well as new, massless particles. The ability to modify existing particles (including atoms with mass) is crucial to conserve charges on a system if a new extra point is added. Seven other, unique "frame" types may be specified, denoting how the extra point will be positioned according to the locations of atoms that do have mass.

For the first example, we will convert TIP3P water into TIP4P-Ew water. Because the TIP3P water molecule is rigid, there are several frame types that would produce equivalent results, and we will choose the simplest of them. The following &rule file will guide the conversion:

&rule
  ExtraPoint  'EPW',
  ResidueName 'WAT',
  FrameType   'Flex-3',
  ParentAtom  'O',
  FrameAtom2  'H1',
  FrameAtom3  'H2',
  Vector12    0.106641378,
  Vector13    0.106641378,
  Charge     -1.04844
&end

&rule
  ResidueName 'WAT',
  FrameType   'None',
  AtomName    'O',
  Charge      0.0,
  Sigma       3.164349,
  Epsilon     0.162750,
&end

&rule
  ResidueName 'WAT',
  FrameType   'None',
  AtomName    'H1',
  Charge      0.52422,
&end

&rule
  ResidueName 'WAT',
  FrameType   'None',
  AtomName    'H2',
  Charge      0.52422,
&end

This input, taking the &rule namelist file to be "tip4p.rules", will then convert any topology containing TIP3P waters:

&files
  -p     tip3p.prmtop
  -c     tip3p.inpcrd
  -xpt   tip4p.rules
&end

&parmedit
  NewPrmtop   tip4pew.prmtop,
  NewInpcrd   tip4pew.inpcrd,
&end

Running the above input as "convert.mdin" with this topology file and these input coordinates produces a new topology file and coordinates wherein all of the waters have been converted to TIP4P-Ew parameters. Note that, while this is a handy trick, it is for demonstrative purposes. There is no advantage to changing the water model in this way rather than loading TIP4P-Ew from the start. Furthermore, there is no way to modify the names of atoms or residues using &rule namelists, if these are important for later analysis. The new topology contains some new fields towards the bottom:

%FLAG VIRTUAL_SITE_FRAMES
%FORMAT(10I8)
       4       1       2       3      -1       5       8       5       6       7
      -1       5      12       9      10      11      -1       5      16      13
 (...)
    2444    2441    2442    2443      -1       5    2448    2445    2446    2447
      -1       5
%FLAG VIRTUAL_SITE_FRAME_DETAILS
%FORMAT(5E16.8)
  1.06641378E-01  1.06641378E-01  0.00000000E+00  1.06641378E-01  1.06641378E-01
  0.00000000E+00  1.06641378E-01  1.06641378E-01  0.00000000E+00  1.06641378E-01
 (...)
  1.06641378E-01  1.06641378E-01  0.00000000E+00  1.06641378E-01  1.06641378E-01
  0.00000000E+00

These new sections will direct pmemd, pmemd.cuda, and mdgx to see the relevant connections between the new, massless atoms and their frames. The VIRTUAL_SITE_FRAMES section contains one tuple of six integers per extra point. They describe the index of the new particle, the indices of up to four frame atoms (marked -1 if there is no second, third, or fourth frame atom), and the frame type index (an internal numbering system dictated by the string found after the obligatory "FrameType" keyword in each &rule). The VIRTUAL_SITE_FRAME_DETAILS dection contains tuples of three reals for every frame, marked with zeros if any of the three numbers is not relevant.

As mentioned before, were the TIP4P-Ew water model to be loaded from the beginning, these new topology segments would not be present and all of the Amber programs would default to their original behavior of inferring the extra point types by seeing atoms with no mass and then looking at connectivities recorded as actual bonds between each extra point and one other atom. Also, note that tip4pew.prmtop contains Lennard-Jones coefficients for three atom types whereas tip3p.prmtop contains coefficients only for two. There are Lennard-Jones ACOEF and BCOEF values for the interactions of TIP3P and TIP4P-Ew oxygen atoms, even though such a mixed system would probably never be built. mdgx has given the re-parameterized oxygen atom a new type index (3), and since all of the oxygens used to be of type 1 there simply are no longer any atoms of this type. The new extra point, like the hydrogens, has no Lennard-Jones properties and takes type index 2, just like them. The two topologies do not work the same way, and the custom extra point format will not work in sander as of yet, but in our faster engines the two topologies will create water with identical properties.

It is also helpful to check the resulting mdout. The topology checking features are engaged by default, and they report that all of the residues in TIP4P-Ew water still have neutral net charges. A total of 612 water molecules received an extra point.

It is more common to want to add extra points to a specific ligand or side chain. The next example will place extra points on a bromobenzene molecule near the bromine sigma hole. It begins with a much more extensive set of &rule namelists. One of the weaknesses of this process is that mdgx will often report "atom not found" if there is a misstatement in the rule namelists. It is therefore helpful to construct the file one namelist at a time to ensure that such a mistake is immediately corrected. Aside from that, mdgx prints many helpful messages to guide the user in stating the input in the best possible syntax.

The associated &parmedit input follows from the previous example. The new topology is not that large, as there is only one residue to apply the extra point onto. Here, the extra points get their own atom type because there is no polar hydrogen, no atom with no Lennard-Jones properties. The first extra point is placed using just two frame atoms, at a fixed distance of -0.8 Angstroms from the bromine along the bromine-carbon bond. Subsequent extra points are placed in the benzene ring's Π-system. As in the previous example, the new coordinate file provides a means for visualizing the result:

Bromobenzene with many extra points

The mdout file again provides some useful information for keeping track of what was added. Unique extra points are identified by unique frame geometries, but extra points with similar geometries and different frame atoms are grouped together for the purposes of this summary. Each extra point is shown beside its parent atom.

 Extra points were added to this topology.

     Frame Type    Type Alias    Inherited     Added
     ----------   ------------   ---------   ---------
       Out-3       OutOfPlane           0          12
        FD-2        FixedDis2           0           1

 Details of added extra points follow.
   Example format:
     [Extra Point]   --  
     [Parent Atom]   --  

 Virtual Site ID    0 (   1 instances):
   Frame type:   FD-2 / FixedDis2
   Charge:      -0.25000
   LJ-Sigma:     0.00000
   LJ-Epsilon:   0.00000
   Example:    [ LIG      0 -- SH1      12 ] -> [ LIG      0 -- Br1      11 ]
 Virtual Site ID    0 (   6 instances):
   Frame type:   Out-3 / OutOfPlane
   Charge:      -0.10000
   LJ-Sigma:     0.00000
   LJ-Epsilon:   0.00000
   Examples:   [ LIG      0 -- PL1A     13 ] -> [ LIG      0 -- C1        4 ]
 Virtual Site ID    0 (   6 instances):
   Frame type:   Out-3 / OutOfPlane
   Charge:      -0.10000
   LJ-Sigma:     0.00000
   LJ-Epsilon:   0.00000
   Examples:   [ LIG      0 -- PL1B     14 ] -> [ LIG      0 -- C1        4 ]

There are many other frame types that mdgx can impart to a topology and that pmemd or pmemd.cuda can then run at the best speed of any Amber simulation. While it is not ideal to have this topology modification uniquely supported by mdgx, the combination of topology and coordinate manipulations needed to make this work was a natural extension of features already in the program. Fitting parameters for models with these extra points is its own challenge, and will be elaborated in a subsequent tutorial. However, mdgx has all of the tools to make this work.

"How's that for maxed out?"

Last modified: Jul 25, 2023