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
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
