[gmx-users] CHARMM nonbonded parameters and grompp output

[Mark Abraham]

Justin A. Lemkul wrote:
> Hi all,
> 
> As a number of others have attempted, I am exploring the use of the CHARMM force
> fields in Gromacs.  I have read about a number of difficulties throughout the
> list archive, but I am seeing something that thusfar it seems no one has
> reported, regarding nonbonded parameters.  Let me tell you what I've done so
> far, and perhaps someone can shed some light on what's going on.
> 
> I downloaded the ffcharmm* files from the old User Contribution site, following
> a link I found in the archive.  I created a .hdb file (for ease of use) and
> edited the .tdb files, as they contained some inconsistencies in atom naming
> and formatting.  I converted the CHARMM27 force field parameters from Alex
> MacKerell's site using a script that was also made available through User
> Contributions, giving me ffcharmm.itp, ffcharmmbon.itp, and ffcharmmnb.itp.  So
> far, so good.  For the moment, I am also approximating the Urey-Bradley
> potential by using distance restraints, using scripts that came along with the
> ffcharmm package.

The two packages are not intended to work together in this way. Both 
have to work around the fact that pdb2gmx will not write U-B potentials 
natively. Yuguang Mu worked around this using distance restraints, 
whereas I'm using the implemented U-B bond type that DvdS mentioned. 
Yuguang's scripts predate mine, and I recall when I implemented the 
same, I had to edit the GROMACS 3.3(?) sources to actually make U-B bond 
types available. Presumably that need motivated his choice of 
work-around. That source editing is no longer needed.

Anyway, the point is that using any of Yuguang's awk or shell scripts is 
not supported with the ffcharmm*itp files produced by my conversion 
script. I've no idea what the combination would do, and no interest in 
finding out! :-) The rudimentary documentation in the headers of my 
scripts discuss the necessary issues.

> I was able to produce a topology for hen egg white lysozyme (a decent system to
> test, I thought), with correct disulfides, charges, etc. and life was good.  I
> ran grompp to attempt an in vacuo energy minimization, and saw this among my
> output:
> 
> processing topology...
> Generated 0 of the 8646 non-bonded parameter combinations
> 
> The rest of the grompp output indicated no errors or warnings (aside from a net
> charge, which is OK for now).  I reasoned that there should be _some_ form of
> nonbonded interactions within the protein, correct?

Yes, but that doesn't mean they have to be "generated". grompp allows 
pairwise interactions to be listed explicitly in [nonbond_params], and 
otherwise generates them from the combination rule chosen in [defaults] 
and the atom type parameters in [atomtypes]. Various force fields 
implemented in GROMACS use different combinations of the above to get 
their jobs done. As I note in the header to my ffcharmmnb.itp, the 
former is necessary in CHARMM, because they have an extra constant 
factor of 2 in their 6-12 potential that can't be worked around with any 
of the GROMACS generating options. So I generate all possible 
combinations in [nonbond_params] and grompp finds those.

>  I read about 1-4
> interactions (i.e. OPLSAA generates these by scaling, so they are not
> explicitly included in ffoplsaanb.itp), but such scaling is reportedly not used
> in the CHARMM force field, and thus it was no surprise to find [ nonbond_params
> ] within ffcharmmnb.itp.  This brings me to my question: is grompp not finding
> these parameters?  And if so, why?

"finding" != "generating". You can gmxdump your output .tpr file to see 
what nonbonded parameters are being used. Remember numbering will then 
start from zero.

CHARMM uses some dedicated 1-4 parameters for some atom types. My 
scripts deal with this automatically; I've no idea what Yuguang's 
scripts do here.

> I noticed that ffcharmmnb.itp is formatted much like the GROMOS force field
> files (ffgmxnb.itp, ffG*nb.itp) and saw something that is not mentioned in the
> manual.  The first line under [ atomtypes ] in each of the GROMOS force fields
> has formatting like:
> 
> ;name  at.num      mass        charge   ptype       c6           c12
>     O    8      15.99940       0.000       A   0.22617E-02   0.74158E-06
> 
> The at. num column is not mentioned in the manual (Gromacs version 3.3).  Is it
> needed for proper interpretation of the *nb.itp files?  Such a column is
> missing in my ffcharmmnb.itp file.  The [ nonbond_params ] and [ pairtypes ]
> sections seem to have correct formatting, so I am not sure that this
> (potential) inconsistency among the [ atomtypes ] is causing the issue.

There are a number of optional fields in entries to this directive that 
allow work-arounds for different force field issues. See comment in 
push_at() in src/kernel/toppush.c

> Thanks for reading (yet another) long email from me, and thanks in advance if
> anyone has any ideas as to what's going on.

No problem. I will eventually get around to better documentation for my 
stuff :-)

Mark