Creating a hydrogen mass repartitioned topology

You want to create a hydrogen mass repartitioned topology using CPPTRAJ

In order to speed up simulations, a trick that has been used for some time is to use a redistribution of the masses of the atoms. This allows to increase the time between each step in the simulation as reported by Roitberg and collaborators. The main idea is to reduce the mass of heavy atoms that are connected to hydrogen in order to reduce their vibrational speed, which allows for a faster time step in the molecular dynamics engine. Normally only the solute atoms are hydrogen-mass repartitioned and the solute atoms are kept the same. An illustration of doing h-mass repartition for the adenine nucleobase is presented in the next figure:

Atom name Atom mass before repartition Atom mass after repartition
C5’ 12.01 7.978
H5’ 1.008 3.024
H5’’ 1.008 3.024
C4’ 12.01 9.994
H4’ 1.008 3.024
C3’ 12.01 9.994
H3’ 1.008 3.024
C2’ 12.01 7.978
H2’ 1.008 3.024
H2’’ 1.008 3.024
C1’ 12.01 9.994
H1’ 1.008 3.024
C8 12.01 9.994
H8 1.008 3.024
C2 12.01 7.978
H2 1.008 3.024
N6 14.01 9.978
H61 1.008 3.024
H62 1.008 3.024

From the table, we can see that some mass has been substracted from the heavy atoms and added to the hydrogens. This will allow to run simulations using a higher timestep in the input file:

dt=0.004         ! Time step using H-mass repartition

This increase in time step speeds up the simulation time.

We will use the topology of a 12-mer DNA with the sequence CGCGAATTCGCG as an example (available for download here). In the command line, start CPPTRAJ:

CPPTRAJ: Trajectory Analysis. V6.4.5
   ___  ___  ___  ___
    | \/ | \/ | \/ |
   _|_/\_|_/\_|_/\_|_

| Date/time: 04/05/22 15:49:07
| Available memory: 2.341 GB

>

Load the topology with the command parm:

CPPTRAJ: Trajectory Analysis. V6.4.5
   ___  ___  ___  ___
    | \/ | \/ | \/ |
   _|_/\_|_/\_|_/\_|_

| Date/time: 04/05/22 15:51:09
| Available memory: 2.349 GB

> parm dna.prmtop
[parm dna.prmtop]
     Reading 'dna.prmtop' as Amber Topology
     Radius Set: modified Bondi radii (mbondi)
>

Lets display some information about the topology with the command parminfo:

CPPTRAJ: Trajectory Analysis. V6.4.5
   ___  ___  ___  ___
    | \/ | \/ | \/ |
   _|_/\_|_/\_|_/\_|_

| Date/time: 04/05/22 15:51:09
| Available memory: 2.349 GB

> parm dna.prmtop
[parm dna.prmtop]
     Reading 'dna.prmtop' as Amber Topology
     Radius Set: modified Bondi radii (mbondi)
> parminfo
[parminfo]
     Topology dna.prmtop contains 16074 atoms.
     Title: default_name
     Original filename: dna.prmtop
     5144 residues.
     5122 molecules.
     16110 bonds (15566 to H, 544 other).
     1476 angles (640 with H, 836 other).
     2814 dihedrals (1294 with H, 1520 other).
     Box: Truncated octahedron
     5098 solvent molecules.
     GB radii set: modified Bondi radii (mbondi)
>

This command will print a summary of information contained in the specified topology. We would like to see the difference of a residue when we apply the hydrogen mass repartition command using CPPTRAJ. To get information of a specific residue, we can use the resinfo command:

> resinfo :5
[resinfo :5]
1 residues selected.
#Res Name First Last Natom #Orig #Mol C I
5 DA 125 156 32 5 1

In this case we used the flag :5 to tell the resinfo command that we want information for residue #5 (an adenine). For more information about how to select chains, residues or atoms, refer to the atom mask selection syntax page.

To see the mass value of each individual atom in residue #5, we use atominfo:

> atominfo :5
[atominfo :5]
32 atoms selected.
#Atom Name  #Res Name #Mol Type Charge  Mass    GBradius  El  rVDW   eVDW
125   P     5    DA   1    P     1.1659 30.9700 1.8500    P   2.1000 0.2000
126   OP1   5    DA   1    O2   -0.7761 16.0000 1.5000    O   1.6612 0.2100
127   OP2   5    DA   1    O2   -0.7761 16.0000 1.5000    O   1.6612 0.2100
128   O5'   5    DA   1    OS   -0.4954 16.0000 1.5000    O   1.6837 0.1700
129   C5'   5    DA   1    CJ   -0.0069 12.0100 1.7000    C   1.9080 0.1094
130   H5'   5    DA   1    H1    0.0754 1.0080  1.3000    H   1.3870 0.0157
131   H5''  5    DA   1    H1    0.0754 1.0080  1.3000    H   1.3870 0.0157
132   C4'   5    DA   1    CT    0.1629 12.0100 1.7000    C   1.9080 0.1094
133   H4'   5    DA   1    H1    0.1176 1.0080  1.3000    H   1.3870 0.0157
134   O4'   5    DA   1    OS   -0.3691 16.0000 1.5000    O   1.6837 0.1700
135   C1'   5    DA   1    CT    0.0431 12.0100 1.7000    C   1.9080 0.1094
136   H1'   5    DA   1    H2    0.1838 1.0080  1.3000    H   1.2870 0.0157
137   N9    5    DA   1    N*   -0.0268 14.0100 1.5500    N   1.8240 0.1700
138   C8    5    DA   1    C2    0.1607 12.0100 1.7000    C   1.9080 0.0860
139   H8    5    DA   1    H5    0.1877 1.0080  1.3000    H   1.3590 0.0150
140   N7    5    DA   1    NB   -0.6175 14.0100 1.5500    N   1.8240 0.1700
141   C5    5    DA   1    CB    0.0725 12.0100 1.7000    C   1.9080 0.0860
142   C6    5    DA   1    CA    0.6897 12.0100 1.7000    C   1.9080 0.0860
143   N6    5    DA   1    N2   -0.9123 14.0100 1.5500    N   1.8240 0.1700
144   H61   5    DA   1    H     0.4167 1.0080  1.3000    H   0.6000 0.0157
145   H62   5    DA   1    H     0.4167 1.0080  1.3000    H   0.6000 0.0157
146   N1    5    DA   1    NC   -0.7624 14.0100 1.5500    N   1.8240 0.1700
147   C2    5    DA   1    CQ    0.5716 12.0100 1.7000    C   1.9080 0.0860
148   H2    5    DA   1    H5    0.0598 1.0080  1.3000    H   1.3590 0.0150
149   N3    5    DA   1    NC   -0.7417 14.0100 1.5500    N   1.8240 0.1700
150   C4    5    DA   1    CB    0.3800 12.0100 1.7000    C   1.9080 0.0860
151   C3'   5    DA   1    C7    0.0713 12.0100 1.7000    C   1.9080 0.1094
152   H3'   5    DA   1    H1    0.0985 1.0080  1.3000    H   1.3870 0.0157
153   C2'   5    DA   1    CT   -0.0854 12.0100 1.7000    C   1.9080 0.1094
154   H2'   5    DA   1    HC    0.0718 1.0080  1.3000    H   1.4870 0.0157
155   H2''  5    DA   1    HC    0.0718 1.0080  1.3000    H   1.4870 0.0157
156   O3'   5    DA   1    OS   -0.5232 16.0000 1.5000    O   1.6837 0.1700

This commands shows for residue 5 the atom name, the residue ID (5), the name of the residue (DA), the molecule number of the command (just 1), the type of atom, the charge for each atom, the mass and other properties. Now lets create the hydrogen-mass repartition topology with the hmassrepartition command. Since we only loaded one topology, the hmassrepartition command does not need any additional arguments:

> hmassrepartition
 [hmassrepartition]
     Performing hydrogen mask repartitioning for atoms selected by mask '*'
     Topology is 'dna.prmtop'
     Hydrogen masses will be changed to 3.024000 amu
     Skipping solvent.
     Mask [*] corresponds to 15596 atoms.
     272 hydrogen masses modified.
     196 heavy atom masses modified.

We can see that the command only adjusted the masses of the solute, leaving the solvent intact. If we run the atominfo command again, we can see that the heavy atoms that are directly bonded with hydrogen, their masses are reduced, and the hydrogen atoms are increased:

> atominfo :5
[atominfo :5]
32 atoms selected.
#Atom Name #Res Name #Mol Type Charge Mass GBradius El rVDW eVDW
125   P    5 DA 1 P 1.1659 30.9700 1.8500 P 2.1000 0.2000
126   OP1  5 DA 1 O2 -0.7761 16.0000 1.5000 O 1.6612 0.2100
127   OP2  5 DA 1 O2 -0.7761 16.0000 1.5000 O 1.6612 0.2100
128   O5'  5 DA 1 OS -0.4954 16.0000 1.5000 O 1.6837 0.1700
129   C5'  5 DA 1 CJ -0.0069 7.9780 1.7000 C 1.9080 0.1094
130   H5'  5 DA 1 H1 0.0754 3.0240 1.3000 H 1.3870 0.0157
131   H5'' 5 DA 1 H1 0.0754 3.0240 1.3000 H 1.3870 0.0157
132   C4'  5 DA 1 CT 0.1629 9.9940 1.7000 C 1.9080 0.1094
133   H4'  5 DA 1 H1 0.1176 3.0240 1.3000 H 1.3870 0.0157
134   O4'  5 DA 1 OS -0.3691 16.0000 1.5000 O 1.6837 0.1700
135   C1'  5 DA 1 CT 0.0431 9.9940 1.7000 C 1.9080 0.1094
136   H1'  5 DA 1 H2 0.1838 3.0240 1.3000 H 1.2870 0.0157
137   N9   5 DA 1 N* -0.0268 14.0100 1.5500 N 1.8240 0.1700
138   C8   5 DA 1 C2 0.1607 9.9940 1.7000 C 1.9080 0.0860
139   H8   5 DA 1 H5 0.1877 3.0240 1.3000 H 1.3590 0.0150
140   N7   5 DA 1 NB -0.6175 14.0100 1.5500 N 1.8240 0.1700
141   C5   5 DA 1 CB 0.0725 12.0100 1.7000 C 1.9080 0.0860
142   C6   5 DA 1 CA 0.6897 12.0100 1.7000 C 1.9080 0.0860
143   N6   5 DA 1 N2 -0.9123 9.9780 1.5500 N 1.8240 0.1700
144   H61  5 DA 1 H 0.4167 3.0240 1.3000 H 0.6000 0.0157
145   H62  5 DA 1 H 0.4167 3.0240 1.3000 H 0.6000 0.0157
146   N1   5 DA 1 NC -0.7624 14.0100 1.5500 N 1.8240 0.1700
147   C2   5 DA 1 CQ 0.5716 9.9940 1.7000 C 1.9080 0.0860
148   H2   5 DA 1 H5 0.0598 3.0240 1.3000 H 1.3590 0.0150
149   N3   5 DA 1 NC -0.7417 14.0100 1.5500 N 1.8240 0.1700
150   C4   5 DA 1 CB 0.3800 12.0100 1.7000 C 1.9080 0.0860
151   C3'  5 DA 1 C7 0.0713 9.9940 1.7000 C 1.9080 0.1094
152   H3'  5 DA 1 H1 0.0985 3.0240 1.3000 H 1.3870 0.0157
153   C2'  5 DA 1 CT -0.0854 7.9780 1.7000 C 1.9080 0.1094
154   H2'  5 DA 1 HC 0.0718 3.0240 1.3000 H 1.4870 0.0157
155   H2'' 5 DA 1 HC 0.0718 3.0240 1.3000 H 1.4870 0.0157
156   O3'  5 DA 1 OS -0.5232 16.0000 1.5000 O 1.6837 0.1700

Next, we save the new trajectory with the parmwrite command and the writedata command:

> parmwrite out repartitioned.prmtop
   [parmwrite out repartitioned.prmtop]
      Writing topology 0 (dna.prmtop) to 'repartitioned.prmtop' with format Amber Topology
      Memory used by full exclusion list: 75.872 kB
> writedata
    [writedata]

This new file (repartitioned.prmtop) is now a topology file that has been hydrogen-mass repartitioned and can be used to run simulations using a 4fs timestep.

This process can also be achieved by creating a text file with the name in.cpptraj with the commands:

parm dna.prmtop
hmassrepartition
parmwrite out repartitioned.prmtop
go
quit

and running CPPTRAJ with the command:

cpptraj -i in.cpptraj