[gmx-users] Box dimension size errors in MARTINI soft core simulation
Scott Pendley
scott.pendley at gmail.com
Tue Jul 30 16:17:43 CEST 2013
Justin,
Thank you for the continued help. I have copied the text of the mdp file
below. It is a little long, so fair warning. The system is a MARTINI
coarse grain representation of 1,2-dipalmitoyl-phosphatidyl-glycerole (LHG)
surrounded by 700 octanol molecules. Being a neutral MARTINI
representation, all represented beads or atoms have a charge of 0.0.
Octanol parameters were taken from MARTINI v 2.0 lipid parameters taken
from the official website. 1,2-dipalmitoyl-phosphatidyl-glycerole
parameters were created in house (and seem to work well in aqueous
simulations). From your note yesterday, it sounded like it may have been a
rough day yesterday. I hope today treats you better.
Cheers,
Scott
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 02 feb 2013 by DdJ
;
; for use with GROMACS 4.5/4.6
;
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable
; with dt=40 fs, some (especially rings and polarizable water) require
20-30 fs.
; Note that time steps of 40 fs and larger may create local heating or
; cooling in your system. Although the use of a heat bath will globally
; remove this effect, it is advised to check consistency of
; your results for somewhat smaller time steps in the range 20-30 fs.
; Time steps exceeding 40 fs should not be used; time steps smaller
; than 20 fs are also not required unless specifically stated in the itp
file.
integrator = sd
dt = 0.002
nsteps = 10000000
nstcomm = 10
comm-grps =
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 100
nstxtcout = 1000
xtc_precision = 100
xtc-grps =
energygrps = LHG OCO
; NEIGHBOURLIST and MARTINI
; Due to the use of shifted potentials, the noise generated
; from particles leaving/entering the neighbour list is not so large,
; even when large time steps are being used. In practice, once every
; ten steps works fine with a neighborlist cutoff that is equal to the
; non-bonded cutoff (1.2 nm). However, to improve energy conservation
; or to avoid local heating/cooling, you may increase the update frequency
; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option
; is computationally less expensive and leads to improved energy
conservation
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.4
; MARTINI and NONBONDED
; Standard cut-off schemes are used for the non-bonded interactions
; in the Martini model: LJ interactions are shifted to zero in the
; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2 nm.
; The treatment of the non-bonded cut-offs is considered to be part of
; the force field parameterization, so we recommend not to touch these
; values as they will alter the overall balance of the force field.
; In principle you can include long range electrostatics through the use
; of PME, which could be more realistic in certain applications
; Please realize that electrostatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit
screening
; in aqueous environment. In this case, the use of PME is more realistic.
;
; For use in combination with the Verlet-pairlist algorithm implemented
; in Gromacs 4.6 a straight cutoff in combination with the potential
; modifiers can be used. Although this will change the potential shape,
; preliminary results indicate that forcefield properties do not change a
lot
; when the LJ cutoff is reduced to 1.1 nm. Be sure to test the effects for
; your particular system. The advantage is a gain of speed of 50-100%.
coulombtype = Shift ;Reaction_field (for use with
Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)
;cutoff-scheme = verlet
;coulomb-modifier = Potential-shift
;vdw-modifier = Potential-shift
;epsilon_rf = 0 ; epsilon_rf = 0 really means epsilon_rf =
infinity
;verlet-buffer-drift = 0.005
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale
(V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10-4 - 10-5 bar-1. Note that, for equilibration purposes,
; the Berendsen thermostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = LHG OCO
tau_t = 1.0 1.0
ref_t = 300 300
Pcoupl = parrinello-rahman
Pcoupltype = semiisotropic
tau_p = 12.0 12.0 ;parrinello-rahman is more stable
with larger tau-p, DdJ, 20130422
compressibility = 3e-4 3e-4
ref_p = 1.0 1.0
gen_vel = no
gen_temp = 320
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
unconstrained_start = no
lincs_order = 4
lincs_warnangle = 30
; and set the free energy parameters
free-energy = yes
couple-moltype = LHG
init-lambda = 0.4
; these 'soft-core' parameters make sure we never get overlapping
; charges as lambda goes to 0
sc-power = 1
sc-sigma = 0.3
sc-alpha = 1.0
; we still want the molecule to interact with itself at lambda=0
couple-intramol = no
couple-lambda0 = vdw-q
couple-lambda1 = none
foreign-lambda = 0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8
0.85 0.9 0.95 1
On Mon, Jul 29, 2013 at 4:39 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> On 7/29/13 2:12 PM, Scott Pendley wrote:
>
>> Thank you, Justin. I am using gromacs version 4.5.5 and have attached
>> .mdp
>> file. I followed your advise and pointer to trouble shooting the system
>>
>
> The mailing list does not accept attachments. Please either copy and
> paste its contents or provide a link from which it can be downloaded.
>
>
> and decomposed the energy to find potential sources for the problem. The
>> simulation ran for 978 ps, the total energy changed from approximately
>> -13000 kj/mol to -15000 kj/mol and plateau'ed out around 900 ps without a
>> big energy blow-out. The majority of the change in energy was from the LJ
>> contribution which is the expected source with the disappearing atoms of
>> my
>> ligand. Of more concern is the change in box coordinates/shape. The
>> simulation cell starts with dimensions 9.45nm x 9.45nm x 9.45nm and ends
>> at
>> 2.8 x 2.8 x 45.8 nm which does seem to give rise to the error that I
>> originally noted. While some change in volume and cell dimensions is
>> expected, these changes seems to be a focused solely on expanding along
>> the
>> Z-coordinate. Is there any way to constrain the ratio of x:y:z
>> coordinates
>> during npt simulations or do you know of something that I may be missing
>> in
>> my simulations?
>>
>>
> The magnitude of change you're seeing suggests extreme instability and
> thus there is no real "hack" to simply "make it work." Isotropic pressure
> coupling is normally sufficient for this purpose, but something seems to be
> going very wrong in your system. The complete .mdp file, and a description
> of what your ligand is (or at least how big) would be helpful. You've got
> a very large box for what one would normally consider a "ligand."
>
> -Justin
>
> Thank you,
>>
>> Scott
>>
>>
>> On Wed, Jul 24, 2013 at 4:59 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>
>>
>>>
>>> On 7/24/13 4:33 PM, Scott Pendley wrote:
>>>
>>> I am fairly new to gromacs and I am trying to run a thermodynamic
>>>> integration simulation of a ligand disappearing in a box of octanol at a
>>>> single set lambda point. I have previous successful nvt and npt runs of
>>>> this system. When I have added the free energy portions to the input
>>>> file,
>>>> I get the following error:
>>>>
>>>> Fatal error:
>>>> One of the box vectors has become shorter than twice the cut-off length
>>>> or
>>>> box_yy-|box_zy| or box_zz has become smaller than the cut-off.
>>>> For more information and tips for troubleshooting, please check the
>>>> GROMACS
>>>> website at http://www.gromacs.org/****Documentation/Errors<http://www.gromacs.org/**Documentation/Errors>
>>>> <http://**www.gromacs.org/Documentation/**Errors<http://www.gromacs.org/Documentation/Errors>
>>>> >
>>>>
>>>>
>>>> This seems unusual. The box dimensions are 9.45 nm x 9.45 x 9.45 nm so
>>>> that is fairly large even accounting for some shrinkage with a
>>>> disappearing
>>>> ligand.
>>>>
>>>>
>>>> The available information suggests the system has become unstable and
>>> is
>>> imploding. See general troubleshooting information at
>>> http://www.gromacs.org/****Documentation/Terminology/**<http://www.gromacs.org/**Documentation/Terminology/**>
>>> Blowing_Up#Diagnosing_an_****Unstable_System<http://www.**
>>> gromacs.org/Documentation/**Terminology/Blowing_Up#**
>>> Diagnosing_an_Unstable_System<http://www.gromacs.org/Documentation/Terminology/Blowing_Up#Diagnosing_an_Unstable_System>
>>> >
>>>
>>> .
>>>
>>>
>>> Cutoffs in the input file are set as follows: rlist =1.4, rvdw=1.2,
>>>
>>>> rcoulomb=1.2. Doubling any of them would still be less than 3 nm which
>>>> is
>>>> significantly smaller than the box size. Is there anything I am missing
>>>> or
>>>> any suggestions that others can give me?
>>>>
>>>>
>>>> I wouldn't mess with the cutoffs; they're an essential part of the
>>> force
>>> field. For further diagnostics, please consider the points above and
>>> provide your .mdp file and Gromacs version.
>>>
>>> -Justin
>>>
>>> --
>>> ==============================****====================
>>>
>>>
>>> Justin A. Lemkul, Ph.D.
>>> Postdoctoral Fellow
>>>
>>> Department of Pharmaceutical Sciences
>>> School of Pharmacy
>>> Health Sciences Facility II, Room 601
>>> University of Maryland, Baltimore
>>> 20 Penn St.
>>> Baltimore, MD 21201
>>>
>>> jalemkul at outerbanks.umaryland.****edu <jalemkul at outerbanks.**
>>> umaryland.edu <jalemkul at outerbanks.umaryland.edu>> | (410)
>>> 706-7441
>>>
>>> ==============================****====================
>>>
>>> --
>>> gmx-users mailing list gmx-users at gromacs.org
>>> http://lists.gromacs.org/****mailman/listinfo/gmx-users<http://lists.gromacs.org/**mailman/listinfo/gmx-users>
>>> <htt**p://lists.gromacs.org/mailman/**listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
>>> >
>>> * Please search the archive at http://www.gromacs.org/**
>>> Support/Mailing_Lists/Search<h**ttp://www.gromacs.org/Support/**
>>> Mailing_Lists/Search<http://www.gromacs.org/Support/Mailing_Lists/Search>>before
>>> posting!
>>>
>>> * Please don't post (un)subscribe requests to the list. Use the www
>>> interface or send it to gmx-users-request at gromacs.org.
>>> * Can't post? Read http://www.gromacs.org/****Support/Mailing_Lists<http://www.gromacs.org/**Support/Mailing_Lists>
>>> <http://**www.gromacs.org/Support/**Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>
>>> >
>>>
>>>
>>>
>>>
> --
> ==============================**====================
>
> Justin A. Lemkul, Ph.D.
> Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 601
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalemkul at outerbanks.umaryland.**edu <jalemkul at outerbanks.umaryland.edu> | (410)
> 706-7441
>
> ==============================**====================
> --
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
> * Please search the archive at http://www.gromacs.org/**
> Support/Mailing_Lists/Search<http://www.gromacs.org/Support/Mailing_Lists/Search>before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>
>
More information about the gromacs.org_gmx-users
mailing list