[gmx-users] Cannot Generate CHARMM36 Topology for POPC with Virtual Sites

[Chris Ing]

I'm trying to build a charmm36.ff topology file for a POPC with virtual
sites. In GROMACS 4.5.5 and 4.6.3 (with the most recent CHARMM36 March
2014), when I run:

pdb2gmx -f POPC.pdb -o POPC.gro -p POPC.top -i POPC.posre -vsite hydrogen
grompp -f empty.mdp -c POPC.gro -p POPC.top -o test.tpr

I get multiple constraint errors:

ERROR 1 [file POPC.top, line 311]:

  No default Constr. No Conn. types

Which is due to the lack of two constraints not in ffbonded.itp:

MCH3 NTL
MCH3 CTL2

MCH3 is the virtual site of a CH3 group, and NTL refers to the choline
nitrogen in the POPC headgroup and CTL2, in this case, is the carbon
closest to the CH3 group on the POPC tail.

The same question was asked in the gmx-users thread "amber99sb_virtual site
lipids dopc MCH3_N constraint types" by Song Ke, where the reply from Erik
Marklund was to manually add this value "from bond lengths, center of mass
for the group and the moment of inertia.". In lieu of wiping an inch of
dust off my classical mechanics textbook, I wondered if I could copy some
additional constraint values from the OPLS ffbonded.itp file.

The charmm36.ff ffbonded.itp file states that the constrainttypes section
is copied directly from OPLS (the restraint values are slightly adjusted?).
Yet, many entries were not copied from OPLS to CHARMM36:

OPLS:

MNH3   CT_2 2    0.158255
MNH3   MNH3 2    0.080236
MCH3A     C 2    0.166040
MCH3A   C_2 2    0.166040
MCH3A    CW 2    0.164312
MCH3A    CV 2    0.164312
MCH3A    CS 2    0.163448
MCH3A    CA 2    0.164888
MCH3A    CB 2    0.164888
MCH3A     N 2    0.159040
MCH3A     S 2    0.193874
MCH3A MCH3A 2    0.092160
MCH3B    CT 2    0.167031
MCH3B  CT_2 2    0.167031
MCH3B  CT_3 2    0.167031
MCH3B    CO 2    0.167031
MCH3B MCH3B 2    0.091456

CHARMM36:

MNH3     CT3   2    0.159603
MNH3    MNH3   2    0.082619
MCH3     CT1   2    0.168122
MCH3     CT2   2    0.167162
MCH3     CT3   2    0.167354
MCH3       S   2    0.195314
MCH3    MCH3   2    0.093582

Would it be acceptable to use these OPLS constraint values in CHARMM36 to
resolve my problem?:

MCH3A     N 2    0.159040
MCH3B    CT 2    0.167031

Thanks gmx-users/Justin Lemkul!

-Chris

P.S. I had similar missing dihedral/ub errors when trying to build my
charmm36.ff protein topology with virtual sites whenever I used a neutral
terminus. I didn't look into the details because using a charged terminus,
I had no errors.