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
[email protected] | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
--
Gromacs Users mailing list
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
mail to [email protected].
