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

[chris.neale at utoronto.ca]

Are the charge groups the same and are they handled the same between  
gromacs and charmm?

Quoting Lee-Ping <leeping at MIT.EDU>:

> 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
>
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.
>