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

[Gerrit Groenhof]

Hallo,

The LJ are computed by gromacs. The trick is very simple: in the  
electrostatic embedding, the coulomb charges are set to zero. Coulomb  
interaction between MM and QM atoms are computed by the QM program.

Is the difference between MM energies occuring in the QM/MM, or also  
in MM calculations.

Gerrit



On 28 Nov 2009, at 03:08, Lee-Ping wrote:

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

--
Gerrit Groenhof
MPI biophysical chemistry
Goettingen
Germany
http://wwwuser.gwdg.de/~ggroenh/