[gmx-developers] Question about gmx/g03 interface; also, getting gmx and charmm to give the same energy.

[Lee-Ping]

Hi there,

I have a GROMACS development project that currently requires the use of both
CHARMM and GROMACS on the same molecular system.  I have not been able to
get the CHARMM and GROMACS energies to agree exactly (for a periodic system
of ~600 atoms, the energy difference is ~10kJ/mol), and for me this is a
cause for concern.

Background:

I've been developing a force-matching method that matches MM energies/forces
to QM (or QMMM) energies/forces by adjusting selected force field
parameters.  The method has the same goal as ffscan, but it is very
different from ffscan in implementation.  There are many helpful features in
the method, and I hope to incorporate it into GROMACS someday.

Conceptually, the algorithm uses a Newton-Raphson algorithm to minimize a
least-squares residual function derived from the difference between the MM
and QM(QMMM) energies/forces.  The MM-QM mean energy gap is subtracted out.

Currently, I use Q-Chem for computing the quantum forces.  When QMMM
energies and forces are desired, I use CHARMM for its Q-Chem interface.
However, I've found that CHARMM and GROMACS give different MM energies when
PBC is turned on.  Consequently, the MM(chm)-MM(gmx) energy difference
contributes to the residual, which is not the desired behavior.

I am considering two things:
 
Option 1:  Eliminate CHARMM from the method by writing a GROMACS interface
to Q-Chem.
 
Out of all the QM methods available to GROMACS, I'm the most familiar with
G03.  I would be comfortable adding a GROMACS/Q-Chem interface by modifying
the G03 interface, but I'm still not 100% sure how the algorithm works after
having read the source code.  

My main question is:  Which program is responsible for calculating the VdW
interactions between the QM and MM atoms?  It seems from reading the
qm_gaussian.c source code that G03 is responsible for this task.  However, I
know for a fact that Q-Chem _cannot_ do this.  

The possible answers are that GROMACS calculates the VdW interaction (in
which case I'd be home free), or G03 calculates it (then I'd have to
implement the QM-MM VdW interaction into GROMACS; I'm not a Q-Chem
developer).  

If somebody could offer advice on how to proceed, it would be awesome.

Option 2:  Investigate why the MM energies in CHARMM and GROMACS do not
agree.
 
I've traced the energy difference to the evaluation of shifted/switched
Coulomb/VdW interactions.  Using a test system of one acetaldehyde molecule
in 210 SPC/E water molecules (in a 2.0nm box), I called GMX with the
following options; PBC is turned off to isolate the source of error.

nstlist     = 1
ns_type     = grid
rlist       = 99.9
coulombtype = shift
vdwtype     = shift
rcoulomb    = 0.80
rvdw        = 0.80
pbc         = no

I used the corresponding CHARMM options (I already carefully converted the
GMX coordinates/topologies to CHARMM format):

NBONDS ATOM FSHIFT CDIE VDW VFSHIFT CUTNB 999.0 CTOFNB 8.0

I then evaluated the potential energies for 3,000 configurations.  The
non-constant energy difference between CHARMM and GROMACS is about 13kJ/mol,
and most of this comes from the Coulomb energy (a smaller part comes from
the VdW energy).  When rcoulomb and rvdw are varied, the error appears to
grow.

As a control, the non-constant energy difference is only 0.01kJ/mol if the
full nonbonded potential is used; here, I use the following options for
CHARMM and GROMACS.

nstlist     = 0
ns_type     = simple
rlist       = 0.0
coulombtype = cut-off
vdwtype     = cut-off
rcoulomb    = 0.0
rvdw        = 0.0
pbc         = no

NBONDS ATOM FSHIFT CDIE VDW VFSHIFT CUTNB 999.0 CTOFNB 998.0

To those well versed in both CHARMM and GROMACS, I'm wondering: Why do the
shifted/switched nonbonded energies differ by so much between the two
programs?  Am I using incorrect run parameters, or do the software packages
calculate the energy in fundamentally different ways?  (I know about the
force shifting functional form in GROMACS, but I can't find the
corresponding functional form in the CHARMM source code.)

Your input is very much appreciated.  Thanks,

- Lee-Ping Wang