[gmx-users] pdb2gmx picks up the wrong .tdb files?

[Justin Lemkul]

On 1/16/20 8:44 AM, András Ferenc WACHA wrote:
> Sorry, now replying to the whole list:
>
> Dear Justin,
>
> I have obtained the base force field from the website of the MacKerell
> Lab
> (http://mackerell.umaryland.edu/download.php?filename=CHARMM_ff_params_files/charmm36-jul2017.ff.tgz).
> I have checked and the merged.n.tdb file is the same in my version and
> in theirs.
>
> Shouldn't pdb2gmx only use the termini databases corresponding to the
> .rtp file from which the given residue has been read? I.e. if it is LYS,
> found in merged.rtp, then only merged.{n,c}.tdb?

The base names should match between .rtp and .tdb but I have never tried 
having multiple types of protein residues in different files.

-Justin

> Kind regards,
>
> Andras
>
> On 1/16/20 2:33 PM, Justin Lemkul wrote:
>>
>> On 1/16/20 8:04 AM, András Ferenc WACHA wrote:
>>> Dear fellow Gromacs users,
>>>
>>> I have developed an extended version of the CHARMM36m force field for
>>> beta-amino acids (https://charmm-betaff.readthedocs.io,
>>> https://gitlab.com/awacha/charmm-beta.ff). For beta-amino acids I found
>>> pdb2gmx works well. However, if I try to use it for natural peptides and
>>> proteins (e.g. 6LYZ from the PDB), the following command:
>>>
>>> gmx pdb2gmx -f 6LYZ.pdb, -ff charmm-beta -ter -ignh
>>>
>>> fails to recognize the correct termini database files. The output is the
>>> following:
>>>
>>> (-------------------------- START of gmx output
>>> -----------------------------------)
>>>
>>> GROMACS:      gmx pdb2gmx, version 2020
>>> Executable:
>>> /opt/spack/opt/spack/linux-archrolling-skylake/gcc-9.2.0/gromacs-2020-ls53yrft6droj2ovhuzdwdufi4an4vir/bin/gmx
>>>
>>> Data prefix:
>>> /opt/spack/opt/spack/linux-archrolling-skylake/gcc-9.2.0/gromacs-2020-ls53yrft6droj2ovhuzdwdufi4an4vir
>>>
>>> Working dir:  /home/wachaandras/gromacs/forcefields
>>> Command line:
>>>     gmx pdb2gmx -f 6lyz.pdb -ff charmm-beta -ter -ignh
>>>
>>> Using the Charmm-beta force field in directory ./charmm-beta.ff
>>>
>>> Opening force field file ./charmm-beta.ff/watermodels.dat
>>>
>>> Select the Water Model:
>>>    1: TIP3P    TIP 3-point, recommended, by default uses CHARMM TIP3 with
>>> LJ on H
>>>    2: TIP4P    TIP 4-point
>>>    3: TIP5P    TIP 5-point
>>>    4: SPC        simple point charge
>>>    5: SPC/E    extended simple point charge
>>>    6: Methanol   Equilibrated methanol
>>>    7: Octanol    Equilibrated octanol
>>>    8: None
>>> 1
>>> going to rename ./charmm-beta.ff/merged.r2b
>>> Opening force field file ./charmm-beta.ff/merged.r2b
>>> Reading 6lyz.pdb...
>>> WARNING: all CONECT records are ignored
>>> Read '', 1001 atoms
>>> Analyzing pdb file
>>> Splitting chemical chains based on TER records or chain id changing.
>>> There are 1 chains and 0 blocks of water and 129 residues with 1001
>>> atoms
>>>
>>>     chain  #res #atoms
>>>     1 'A'   129   1001
>>>
>>> All occupancies are one
>>> Opening force field file ./charmm-beta.ff/atomtypes.atp
>>> Opening force field file ./charmm-beta.ff/betaaminoacids.atp
>>> Opening force field file ./charmm-beta.ff/cyclicbeta.atp
>>> Reading residue database... (Charmm-beta)
>>> Opening force field file ./charmm-beta.ff/betaaminoacids.rtp
>>> Opening force field file ./charmm-beta.ff/cyclicbeta.rtp
>>> Opening force field file ./charmm-beta.ff/extra.rtp
>>> Warning: file does not end with a newline, last line:
>>> ;
>>> Opening force field file ./charmm-beta.ff/merged.rtp
>>> Opening force field file ./charmm-beta.ff/merged.hdb
>>> Opening force field file ./charmm-beta.ff/betaaminoacids.n.tdb
>>> Opening force field file ./charmm-beta.ff/cyclicbeta.n.tdb
>>> Opening force field file ./charmm-beta.ff/merged.n.tdb
>>> Opening force field file ./charmm-beta.ff/betaaminoacids.c.tdb
>>> Opening force field file ./charmm-beta.ff/cyclicbeta.c.tdb
>>> Opening force field file ./charmm-beta.ff/merged.c.tdb
>>> Processing chain 1 'A' (1001 atoms, 129 residues)
>>> Analysing hydrogen-bonding network for automated assignment of histidine
>>>    protonation. 213 donors and 184 acceptors were found.
>>> There are 262 hydrogen bonds
>>> Will use HISD for residue 15
>>> Identified residue LYS1 as a starting terminus.
>>> Identified residue LEU129 as a ending terminus.
>>> 8 out of 8 lines of specbond.dat converted successfully
>>> Special Atom Distance matrix:
>>>                       CYS6   MET12   HIS15   CYS30   CYS64   CYS76
>>> CYS80
>>>                       SG48    SD87  NE2118   SG238   SG513   SG601
>>> SG630
>>>      MET12    SD87   1.205
>>>      HIS15  NE2118   1.804   0.998
>>>      CYS30   SG238   1.450   1.070   2.063
>>>      CYS64   SG513   2.873   1.779   1.755   2.236
>>>      CYS76   SG601   2.740   1.544   1.418   2.132   0.778
>>>      CYS80   SG630   3.006   1.943   1.892   2.387   0.200   0.958
>>>      CYS94   SG724   2.576   1.388   1.329   1.976   0.678   0.211
>>> 0.871
>>>     MET105   SD799   1.866   0.933   1.679   0.896   1.858   1.483
>>> 2.047
>>>     CYS115   SG889   1.625   1.099   2.068   0.204   2.108   2.005
>>> 2.260
>>>     CYS127   SG981   0.203   1.089   1.733   1.332   2.808   2.631
>>> 2.950
>>>                      CYS94  MET105  CYS115
>>>                      SG724   SD799   SG889
>>>     MET105   SD799   1.389
>>>     CYS115   SG889   1.855   0.790
>>>     CYS127   SG981   2.474   1.701   1.505
>>> Linking CYS-6 SG-48 and CYS-127 SG-981...
>>> Linking CYS-30 SG-238 and CYS-115 SG-889...
>>> Linking CYS-64 SG-513 and CYS-80 SG-630...
>>> Linking CYS-76 SG-601 and CYS-94 SG-724...
>>> Select start terminus type for LYS-1
>>>    0: Beta3NH3+
>>>    1: Beta2NH3+
>>>    2: Beta23NH3+
>>>    3: B3_NH2
>>>    4: B2_NH2
>>>    5: B23_NH2
>>>    6: 5TER
>>>    7: None
>>>
>>> (-------------------------- END of gmx output
>>> -----------------------------------)
>>>
>>> Information in the force field relevant to the beta-peptides is stored
>>> in betaaminoacids.{arn,atp,c.tdb,n.tdb,rtp} files. I have kept the
>>> original ones (merged.*) intact. From the output above it seems that
>>> pdb2gmx recognizes and loads all the relevant files (both
>>> betaaminoacids.* and merged.*), but only shows a selection of N-terminal
>>> patches from the betaaminoacids.n.tdb file. More precisely, the 5TER
>>> terminus is from merged.n.tdb, but the other ones before it get
>>> overridden by entries from betaaminoacids.n.tdb. Is there some rule in
>>> pdb2gmx which considers some endgroups equivalent and keeps only one? I
>>> have encountered it neither in the main documentation (section 5.6.5:
>>> pdb2gmx input files), nor in the help of gmx pdb2gmx.
>>>
>>> Could anyone point me to how I should do this correctly?
>> What is the source of these force field files? Based on the .tdb
>> entries, it seems they do not include the "standard" patches that
>> should be included with the official CHARMM force field.
>>
>> -Justin
>>
>>

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

==================================================