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

[Ryan Muraglia]

On Fri, Mar 18, 2016 at 7:47 AM, <
gromacs.org_gmx-users-request at maillist.sys.kth.se> wrote:
>
> Message: 4
> Date: Fri, 18 Mar 2016 07:46:48 -0400
> From: Justin Lemkul <jalemkul at vt.edu>
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] Domain decomposition error tied to free
>         energy perturbation
> Message-ID: <56EBEAA8.1050309 at vt.edu>
> Content-Type: text/plain; charset=windows-1252; format=flowed
>
>
>
> 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
>
> ==================================================
>
>
EM does not run when I include the free energy parameters, but does run for
my mutation topology when I exclude them, so I know the issue was tied to
those parameters. Somehow I'd overlooked the couple-intramol flag.

This appears to have fixed the problem. Thank you very much!


-- 
Ryan Muraglia
rmuraglia at gmail.com