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

[Lee-Ping]

In my test case, the charge groups are the same - the molecules are small
enough that each molecule is its own charge group.  My guess is that the
energy discrepancy isn't coming from counting up the charge groups (in my
test case, I made sure that every interaction was counted.)

Regarding whether CHARMM and GROMACS handle charge groups in the same way,
I'm not 100% sure, but my impression is that the charge groups are handled
very similarly.  

-----Original Message-----
From: gmx-developers-bounces at gromacs.org
[mailto:gmx-developers-bounces at gromacs.org] On Behalf Of
chris.neale at utoronto.ca
Sent: Friday, November 27, 2009 9:29 PM
To: gmx-developers at gromacs.org
Subject: Re: [gmx-developers] Question about gmx/g03 interface; also,getting
gmx and charmm to give the same energy.

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.
>



-- 
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.