[gmx-users] dialanine using charmm27 ff

[Mark Abraham]

On 10/09/2011 10:23 AM, Sandeep Somani wrote:
> Hi Mark
>
> It worked ! pdb2gmx is now properly processing the pdb, and moreover 
> the potential energy from gmx is within 1.27 kJ/mol of that from 
> Charmm. More on energy comparison shortly.
>
> However, further tinkering of the rtp and pdb file was required to 
> generate the topology:
> 1.
> The original rtp entry which had HT1/2/3 for the methyl hydrogens was 
> giving the error
> "
> Fatal error:
> Atom H1 in residue ACE 1 was not found in rtp entry ACE with 6 atoms
> while sorting atoms.
>  "
> even though neither rtp nor pdb had H1 atom. I do not understand why 
> it was insisting on H1/2/3. It worked upon changing HT1/2/3 to H1/2/3.

Yeah. There's a bizarre atom-name-translation feature that was breaking 
here. That'll teach me not to do actual testing.

>
> 2.
> I changed bond entry  'N CAT'  -> 'N CT' in rtp entry of CT3. Guess 
> CAT was a typo.
>
> 3.
> The improper dihedral in CT3
> '      -C      CT     N       -O ; commented to match with charmm psf '
> was absent from the charmm psf. So I commented it.
>
> fyi - final pdb and rtp entries are listed below.

Here's an updated patch that really does work :)

diff --git a/share/top/charmm27.ff/aminoacids.hdb 
b/share/top/charmm27.ff/aminoacids.hdb
index 84e0859..7a9a4ea 100644
--- a/share/top/charmm27.ff/aminoacids.hdb
+++ b/share/top/charmm27.ff/aminoacids.hdb
@@ -40,6 +40,9 @@ ASPP    4
  1      5       HA      CA      N       C       CB
  2      6       HB      CB      CG      CA
  1      2       HD2     OD2     CG      CB
+CT3     2
+1      1       HN      N       -C      CH3
+3       4       HH3    CH3     N       HN
  CYS2    3
  1      1       HN      N       -C      CA
  1      5       HA      CA      N       C       CB
diff --git a/share/top/charmm27.ff/aminoacids.rtp 
b/share/top/charmm27.ff/aminoacids.rtp
index f9f6823..0b9978f 100644
--- a/share/top/charmm27.ff/aminoacids.rtp
+++ b/share/top/charmm27.ff/aminoacids.rtp
@@ -1672,3 +1672,43 @@
  [ ZN2 ]
   [ atoms ]
         ZN      ZN      2.00    0
+
+[ ACE ]
+ [ atoms ]
+        CH3     CT3     -0.270  0
+        HH31    HA      0.090   1
+        HH32    HA      0.090   2
+        HH33    HA      0.090   3
+        C       C       0.510   4
+        O       O       -0.510  5
+ [ bonds ]
+        C       CH3
+        C       +N
+        CH3     HH31
+        CH3     HH32
+        CH3     HH33
+        O       C
+ [ impropers ]
+        C       CH3     +N      O
+
+[ CT3 ]
+; this can also be done with the .c.tdb, but the atom naming is different
+; and this can matter
+ [ atoms ]
+        N       NH1     -0.470  0
+        HN      H       0.310   1
+        CH3     CT3     -0.110  2
+        HH31    HA      0.090   3
+        HH32    HA      0.090   4
+        HH33    HA      0.090   5
+ [ bonds ]
+        -C      N
+        N       HN
+        N       CH3
+        CH3     HH31
+        CH3     HH32
+        CH3     HH33
+
+ [ impropers ]
+        N       -C      CH3     HN
+        -C      CH3     N       -O
diff --git a/share/top/residuetypes.dat b/share/top/residuetypes.dat
index 05bcbcf..1e8d527 100644
--- a/share/top/residuetypes.dat
+++ b/share/top/residuetypes.dat
@@ -10,6 +10,7 @@ ASP   Protein
  ASP1   Protein
  ASPH   Protein
  ASH    Protein
+CT3    Protein
  CYS    Protein


>
> Now, regarding the energies, I compared the different bonded and 
> non-bonded energy components from charmm with those from gmx.
> Bonds, dihedral and coulomb energies matched to within 0.01 kJ/mol 
> which was great!
> The discrepancy is in angles (U-B) where gmx value is higher by 1.6652 
> kJ/mol and LJ where gmx is lower by 0.459 kJ/mol.

The label CHARMM and GROMACS use for UB varies. The former includes only 
the harmonic distance component of the U-B, the latter includes both.

>
> Any idea on how to fix the angles (not too bothered by LJ) ?
> I checked that the total number of angles in charmm psf and gmx top 
> are the same (=54), but havn't yet compared individual entries.
>
> btw, would you (or someone else here) have a sample charmm (inp, crd, 
> psf) and gmx (mdp, gro, top) files where the two energies show a 
> better match? It would help me with this troubleshooting.

I don't. Ask the authors of the port of this forcefield to GROMACS 
(Bjelkemar et al JCTC)

Mark