[gmx-users] Domain decomposition error tied to free energy perturbation

[Justin Lemkul]

On 3/17/16 8:21 PM, Ryan Muraglia wrote:
> Hello,
>
> I have been attempting to carry out some free energy calculations, but to
> verify the sanity of my parameters, I decided to test them on a structure I
> knew to be stable -- the lysozyme from Lemkul's lysozyme in water tutorial.
>
> I chose the L75A mutation because it is out on the surface to minimize the
> "difficulty of the transformation."
> Using my regular mdp file (even with my mutatation topology generated with
> the pmx package), my minimization runs to completion with no errors.
>
> Once I introduce the following lines to my mdp file:
>
> "
> ; Free energy calculations
> free_energy = yes
> delta_lambda = 0 ; no Jarzynski non-eq
> calc_lambda_neighbors = 1 ; only calculate energy to immediate neighbors
> (suitable for BAR, but MBAR needs all)
> sc-alpha                 = 0.5
> sc-coul                  = no
> sc-power                 = 1.0
> sc-sigma                 = 0.3
> couple-moltype           = Protein_chain_A  ; name of moleculetype to
> decouple
> couple-lambda0           = vdw-q      ; all interactions
> couple-lambda1           = vdw     ; remove electrostatics, only vdW
> couple-intramol          = no
> nstdhdl                  = 100
>
> ; lambda vectors ; decharging only.
> ; init_lambda_state   0   1   2   3   4   5   6   7   8   9   10
> init_lambda_state = 00
> coul_lambdas =        0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
> vdw_lambdas =         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
> bonded_lambdas =      0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ; match
> vdw
> mass_lambdas =        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ; match
> vdw
> temperature_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ; not
> doing simulated tempering
> "
>
> I notice two things:
> 1) Running grompp to generate the tpr file takes much longer
> 2) The minimization fails to run due the following error related to domain
> decomposition:
>
> "
> Fatal error:
> There is no domain decomposition for 4 ranks that is compatible with the
> given box and a minimum cell size of 5.51109 nm
> Change the number of ranks or mdrun option -rdd
> Look in the log file for details on the domain decomposition
> "
>
> I noted that it lists a two-body bonded interaction with a strangely large
> distance:
>
> "
> Initial maximum inter charge-group distances:
>      two-body bonded interactions: 5.010 nm, LJC Pairs NB, atoms 1074 1937
>    multi-body bonded interactions: 0.443 nm, Proper Dih., atoms 1156 1405
> Minimum cell size due to bonded interactions: 5.511 nm
> "
>
> Atom 1074 corresponds to a hydrogen off the beta-carbon of proline 70, and
> atom 1937 refers to a hydrogen on arginine 128. Neither residue is part of
> the protein that is being mutated, and they certainly should not be bonded.
> The [bonds] directive in the topology confirms that there should be no
> interaction between these atoms.

With "couple-intramol = no" (from the manual):

"All intra-molecular non-bonded interactions for moleculetype couple-moltype are 
replaced by exclusions and explicit pair interactions."

So you have a much larger distance for intramolecular interactions, hence DD 
complains and you are more limited in the number of DD cells that can be 
constructed.  Trying to decouple an entire protein chain is (1) not usually 
reasonable and (2) fraught with algorithmic challenges.

> To force the run to begin to get more information on the nature of the
> error, I gave mdrun the -nt 1 option, and got the following warning at the
> beginning of the minimization (which goes on to end prematurely prior to
> reaching the desired Fmax):
>
> "
> WARNING: Listed nonbonded interaction between particles 1 and 195
> at distance 2.271 which is larger than the table limit 2.200 nm.
> "
>
> I'm at a loss in terms of understanding why the addition of my FEP
> parameters is causing this error, and appears to be causing the grompp
> parser to decide that there is a bond where there shouldn't be, forcing the

It's not magically creating bonds; see above.  grompp is taking forever because 
it has to generate a massive list of exclusions and pairs.

> minimimum box size to exceed what makes sense for domain decomposition.
>

If EM fails, that's usually a dead giveaway that either the topology is unsound 
or the initial coordinates are unsuitable in some way.  Without more 
information, it's hard to guess at what's going on.  Does EM proceed without the 
free energy options turned on?

-Justin

> Additional information that may be relevant: I am using the amber99sb
> forcefield with explicit tip3p waters. I am attempting steepest descent
> minimization. rcoulomb and rvdw are both set to 1.2.
>
> Any advice would be greatly appreciated. Thank you!
>
>



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

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

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