[gmx-developers] MD run blows up!

[Justin Lemkul]

Please do not post usage problems to the development list.  I am CC'ing this 
message over to gmx-users.  Please continue the discussion there.

On 9/14/16 5:58 PM, Carlos ES Bernardes wrote:
> Dear Sir or Madam:
>
>
>
> I’m a new GROMACS user, and I’m interested in the study of solutions of small
> molecules, like ionic liquids.
>
> At this moment I’m trying to ma
>
>
>
> ke a simulation of a molecule called bicalutamide. As a first step, I tested the
> force field (OPLS-AA) with DL-POLY, and everything run without problems.
> However, when I changed to GROMACS, the simulations start to blow up.
>
> To understand what was happening, I made the input files only for one
> bicalutamide molecule and I run a small simulation. The files I used are
> attached to this e-mail. The analysis of the simulation results revealed that
> the problem is related with two hydrogen atoms (one connected to a nitrogen
> atoms and the other belonging to a hydroxyl group), that come very close to each
> other and, when this happens, the system goes crazy! What eludes me, is the fact
> that the parametrization works really good with DL-POLY (being able to reproduce
> the crystal structures parameters and the energetics of this molecule very
> accurately), but fails with GROMACS. Furthermore, when I compare the single
> point energies computed by the two programs for the same starting configuration,
> the following differences are obtained (data in kJ/mol):
>
>
>
> Bond
>
> 	
>
> Angle
>
> 	
>
> Dih.
>
> 	
>
> LJ
>
> 	
>
> Coulomb
>
> 	
>
> Potential
>
> 0.02
>
> 	
>
> -0.04
>
> 	
>
> -0.05
>
> 	
>
> -98.24
>
> 	
>
> 59.97
>
> 	
>
> -158.28
>
>
>
> i.e., the bonds, angles and dihedrals energies computed by the two programs are
> approximately the same, but the Coulomb and LJ are significantly different.
> Initially, I thought this could be related with the 1-4 interactions, but I did
> not find anything suspicious error in the input files (maybe because I’m new
> with GROMACS). However, if this was the problem, it does not explain the strange
> attraction between two highly positively charged hydrogen atoms (~+0.3e each),
> since they are 5 bonds apart from each other.
>
> At this point, I do not know what I can check more, and I’m not sure if the
> problem is related with my input files or if is something intrinsic of GROMACS
> and the routines LINCS or SHAKE (the behavior is similar with both algorithms).
> Can you give me a help on this?
>

This is because of 1-4 interactions.  Your topology defines no [pairs] so 
everything is being treated like normal nonbonded interactions (note that you 
only have Coul(SR) and LJ(SR) as energy terms, when you should also have 1-4 
terms).

In OPLS-AA, the 1-4 interactions (pairs) have charges and LJ scaled by 0.5 
(FudgeLJ and FudgeQQ).  Since you define no pairs, the topology is incorrect and 
leads to overly strong 1-4 interactions.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================