[gmx-users] dialanine using charmm27 ff

Sandeep Somani ssomani at gmail.com
Sat Sep 10 02:23:40 CEST 2011


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.

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.

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.

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.

Thanks for your help
Sandeep

================final rtp entries============================
[ ACE ]
 [ atoms ]
       CT3     CT3     -0.270  0
       H1     HA      0.090   1
       H2     HA      0.090   2
       H3     HA      0.090   3
       C       C       0.510   4
       O       O       -0.510  5
 [ bonds ]
       C       CT3
       C       +N
       CT3     H1
       CT3     H2
       CT3     H3
       O       C
 [ impropers ]
       C       CT3     +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
       CT      CT3     -0.110  2
       H1     HA      0.090   3
       H2     HA      0.090   4
       H3     HA      0.090   5
 [ bonds ]
       -C      N
       N       HN
       N       CT
       CT      H1
       CT      H2
       CT      H3
 [ impropers ]
       N       -C      CT     HN
;      -C      CT     N       -O ; commented to match with charmm psf
==============final rtp entries==================

==============final pdb======================
REMARK  NONE *

REMARK   DATE:     8/10/11      5:30:23      CREATED BY USER: ss2029

ATOM      1  CT3 ACE     1      -2.160   0.537   0.930  1.00  0.00      DIAL
ATOM      2  H1  ACE     1      -2.466   0.003   0.005  1.00  0.00      DIAL
ATOM      3  H2  ACE     1      -2.562   1.572   0.910  1.00  0.00      DIAL
ATOM      4  H3  ACE     1      -2.562   0.001   1.816  1.00  0.00      DIAL
ATOM      5  C   ACE     1      -0.672   0.582   1.009  1.00  0.00      DIAL
ATOM      6  O   ACE     1      -0.105   1.128   1.954  1.00  0.00      DIAL
ATOM      7  N   ALA     2       0.000   0.000   0.000  1.00  0.00      DIAL
ATOM      8  HN  ALA     2      -0.453  -0.444  -0.769  1.00  0.00      DIAL
ATOM      9  CA  ALA     2       1.459   0.000   0.000  1.00  0.00      DIAL
ATOM     10  HA  ALA     2       1.812  -0.497   0.897  1.00  0.00      DIAL
ATOM     11  CB  ALA     2       1.949  -0.834  -1.207  1.00  0.00      DIAL
ATOM     12  HB1 ALA     2       1.514  -1.854  -1.155  1.00  0.00      DIAL
ATOM     13  HB2 ALA     2       1.628  -0.373  -2.167  1.00  0.00      DIAL
ATOM     14  HB3 ALA     2       3.056  -0.932  -1.214  1.00  0.00      DIAL
ATOM     15  C   ALA     2       2.096   1.401   0.000  1.00  0.00      DIAL
ATOM     16  O   ALA     2       1.425   2.432   0.000  1.00  0.00      DIAL
ATOM     17  N   ALA     3       3.451   1.458   0.000  1.00  0.00      DIAL
ATOM     18  HN  ALA     3       3.954   0.594   0.000  1.00  0.00      DIAL
ATOM     19  CA  ALA     3       4.276   2.664   0.000  1.00  0.00      DIAL
ATOM     20  HA  ALA     3       4.065   3.236  -0.897  1.00  0.00      DIAL
ATOM     21  CB  ALA     3       3.864   3.539   1.207  1.00  0.00      DIAL
ATOM     22  HB1 ALA     3       2.776   3.756   1.155  1.00  0.00      DIAL
ATOM     23  HB2 ALA     3       4.063   3.014   2.167  1.00  0.00      DIAL
ATOM     24  HB3 ALA     3       4.407   4.508   1.214  1.00  0.00      DIAL
ATOM     25  C   ALA     3       5.791   2.399   0.000  1.00  0.00      DIAL
ATOM     26  O   ALA     3       6.597   3.328   0.000  1.00  0.00      DIAL
ATOM     27  N   CT3     4       6.175   1.110   0.000  1.00  0.00      DIAL
ATOM     28  HN  CT3     4       5.528   0.351   0.000  1.00  0.00      DIAL
ATOM     29  CT  CT3     4       7.566   0.777   0.000  1.00  0.00      DIAL
ATOM     30  H1  CT3     4       8.048   1.201   0.907  1.00  0.00      DIAL
ATOM     31  H2  CT3     4       7.683  -0.327   0.000  1.00  0.00      DIAL
ATOM     32  H3  CT3     4       8.048   1.201  -0.907  1.00  0.00      DIAL
TER      33      CT3      4
END
==============final pdb======================












On Fri, Sep 9, 2011 at 11:35 AM, Mark Abraham <mark.abraham at anu.edu.au>wrote:

>
>
> On 09/09/11, *Sandeep Somani *<ssomani at gmail.com> wrote:
>
> Hi Mark
>
> Thanks for the rtp definitions. I added the rtp defs to aminoacids.rtp and
> renamed residue and atom names in pdb file.
>
> But pdb2gmx is still complaining about the terminal groups. It seems to
> have picked up ACE, but not CT3.
> pdb2gmx error message and modified pdb file are listed below. Any idea ? Do
> I need to modify .c.tdb ? If so, do you have that entry as well :) ?
>
>
> Ah, no. To get rid of the warning, you need to teach GROMACS that CT3 is a
> protein residue. (It knew about ACE from other force fields.)
> Copy residuetypes.dat from the share/top directory to your working
> directory, and add CT3 as Protein at the end.
>
> To get rid of the error, you need to tell pdb2gmx not to try to change your
> termini. Use pdb2gmx -ter, and choose "None" for both N and C termini.
>
> Mark
>
>
>
> Thanks
> Sandeep
>
> =============pdb2gmx error=======================================
> Back Off! I just backed up topol.top to ./#topol.top.3#
> Processing chain 1 (32 atoms, 4 residues)
> There are 3 donors and 3 acceptors
> There are 4 hydrogen bonds
> Identified residue ACE1 as a starting terminus.
> Warning: Residue CT32 in chain has different type (Other) from starting
> residue ACE1 (Protein).
> Identified residue ALA2 as a ending terminus.
> 8 out of 8 lines of specbond.dat converted successfully
> Start terminus ACE-1: NH3+
> End terminus ALA-2: COO-
>
> -------------------------------------------------------
> Program pdb2gmx_d, VERSION 4.5.4
> Source code file: pdb2top.c, line: 1070
>
> Fatal error:
> atom N not found in buiding block 1ACE while combining tdb and rtp
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
> -------------------------------------------------------
> ================pdb2gmx error====================================
>
> ================modified pdb=====================================
>  REMARK  NONE *
>
> REMARK   DATE:     8/10/11      5:30:23      CREATED BY USER: ss2029
>
> *ATOM      1  CT3 ACE     1      -2.160   0.537   0.930  1.00  0.00
>  DIAL*
> *ATOM      2  HT1 ACE     1      -2.466   0.003   0.005  1.00  0.00
>  DIAL*
> *ATOM      3  HT2 ACE     1      -2.562   1.572   0.910  1.00  0.00
>  DIAL*
> *ATOM      4  HT3 ACE     1      -2.562   0.001   1.816  1.00  0.00
>  DIAL*
> *ATOM      5  C   ACE     1      -0.672   0.582   1.009  1.00  0.00
>  DIAL*
> *ATOM      6  O   ACE     1      -0.105   1.128   1.954  1.00  0.00
>  DIAL*
> ATOM      7  N   ALA     1       0.000   0.000   0.000  1.00  0.00
>  DIAL
> ATOM      8  HN  ALA     1      -0.453  -0.444  -0.769  1.00  0.00
>  DIAL
> ATOM      9  CA  ALA     1       1.459   0.000   0.000  1.00  0.00
>  DIAL
> ATOM     10  HA  ALA     1       1.812  -0.497   0.897  1.00  0.00
>  DIAL
> ATOM     11  CB  ALA     1       1.949  -0.834  -1.207  1.00  0.00
>  DIAL
> ATOM     12  HB1 ALA     1       1.514  -1.854  -1.155  1.00  0.00
>  DIAL
> ATOM     13  HB2 ALA     1       1.628  -0.373  -2.167  1.00  0.00
>  DIAL
> ATOM     14  HB3 ALA     1       3.056  -0.932  -1.214  1.00  0.00
>  DIAL
> ATOM     15  C   ALA     1       2.096   1.401   0.000  1.00  0.00
>  DIAL
> ATOM     16  O   ALA     1       1.425   2.432   0.000  1.00  0.00
>  DIAL
> ATOM     17  N   ALA     2       3.451   1.458   0.000  1.00  0.00
>  DIAL
> ATOM     18  HN  ALA     2       3.954   0.594   0.000  1.00  0.00
>  DIAL
> ATOM     19  CA  ALA     2       4.276   2.664   0.000  1.00  0.00
>  DIAL
> ATOM     20  HA  ALA     2       4.065   3.236  -0.897  1.00  0.00
>  DIAL
> ATOM     21  CB  ALA     2       3.864   3.539   1.207  1.00  0.00
>  DIAL
> ATOM     22  HB1 ALA     2       2.776   3.756   1.155  1.00  0.00
>  DIAL
> ATOM     23  HB2 ALA     2       4.063   3.014   2.167  1.00  0.00
>  DIAL
> ATOM     24  HB3 ALA     2       4.407   4.508   1.214  1.00  0.00
>  DIAL
> ATOM     25  C   ALA     2       5.791   2.399   0.000  1.00  0.00
>  DIAL
> ATOM     26  O   ALA     2       6.597   3.328   0.000  1.00  0.00
>  DIAL
> *ATOM     27  N   CT3     2       6.175   1.110   0.000  1.00  0.00
>  DIAL*
> *ATOM     28  HN  CT3     2       5.528   0.351   0.000  1.00  0.00
>  DIAL*
> *ATOM     29  CT  CT3     2       7.566   0.777   0.000  1.00  0.00
>  DIAL*
> *ATOM     30  HT1 CT3     2       8.048   1.201   0.907  1.00  0.00
>  DIAL*
> *ATOM     31  HT2 CT3     2       7.683  -0.327   0.000  1.00  0.00
>  DIAL*
> *ATOM     32  HT3 CT3     2       8.048   1.201  -0.907  1.00  0.00
>  DIAL*
> TER      33      CT3      2
> END
> ================modified pdb=====================================
>
>
>
>
> On Fri, Sep 9, 2011 at 2:52 AM, Mark Abraham <Mark.Abraham at anu.edu.au>wrote:
>
>> On 9/09/2011 8:05 AM, Sandeep Somani wrote:
>>
>>> Hi
>>>
>>> I am trying to set up a simulation for dialanine using charmm27 ff but am
>>> getting errors in pdb2gmx due to missing definitions in rtp file.
>>>
>>> I created the pdb file (below) using Charmm for the sequence
>>> ACE-ALA-ALA-CT3.
>>>
>>> pdb2gmx recognizes ALA atoms but not those of the terminal groups ACE and
>>> CT3.
>>>
>>
>> Yes, CHARMM in GROMACS has lacked these for some time. I asked over a year
>> ago for them to be added, but that hasn't happened. You can add
>>
>> [ ACE ]
>>  [ atoms ]
>>        CT3     CT3     -0.270  0
>>        HT1     HA      0.090   1
>>        HT2     HA      0.090   2
>>        HT3     HA      0.090   3
>>        C       C       0.510   4
>>        O       O       -0.510  5
>>  [ bonds ]
>>        C       CT3
>>        C       +N
>>        CT3     HT31
>>        CT3     HT32
>>        CT3     HT33
>>        O       C
>>  [ impropers ]
>>        C       CT3     +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
>>        CT      CT3     -0.110  2
>>        HT1     HA      0.090   3
>>        HT2     HA      0.090   4
>>        HT3     HA      0.090   5
>>  [ bonds ]
>>        -C      N
>>        N       HN
>>        N       CAT
>>        CT      HT1
>>        CT      HT2
>>        CT      HT3
>>
>>  [ impropers ]
>>        N       -C      CAT     HN
>>        -C      CAT     N       -O
>>
>> to the end of aminoacids.rtp in a local copy of the charmm27.ff folder to
>> make this work.
>>
>>
>>
>>> Does someone have the appropriate definitions for these residues ?
>>> ( I have tried http://swissparam.ch/ to generate an itp file for the
>>> molecule, but upon using it in GMX it does not give the same energy as that
>>> from charmm. )
>>>
>>
>> That will likely be a separate issue.
>>
>> Mark
>> --
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
>> Please search the archive at http://www.gromacs.org/**
>> Support/Mailing_Lists/Search<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-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>
>>
>
>
>
> --
> gmx-users mailing list    gmx-users at 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-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110910/8241c038/attachment.html>


More information about the gromacs.org_gmx-users mailing list