[gmx-users] Box dimension size errors in MARTINI soft core simulation
Scott Pendley
scott.pendley at gmail.com
Thu Aug 1 17:35:34 CEST 2013
Justin,
Much appreciation for the suggestions. I found that the large time
constant (tau_p) on the pressure coupling was indeed a significant
contributor. When I dropped that value down to 2.5, the box vectors scaled
down in a more reasonable manner maintaining the cubic simulation box.
With regard to your other questions:
Visualizing is possible but is more technical to do with MARTINI beads. I
am still trying to figure out how to do this on vmd.
The prior npt run also used parrinello-rahman but with a tau-p of 2.0.
This may have contributed to some of the instability seen in this system so
I may switch it to berendsen per your suggestion.
I agree, isotropic coupling would be better for this system.
Some initial runs suggests that turning off the free energy portion does
allow the simulations to run to completion (for short simulations of 10 ns)
but I haven't analyzed them to determine if the sampling time is
insufficient to reveal developing instabilities.
Contrary to logic, at lambda = 0 this simulation instability also occurs.
The LJ changes are composed on contributions from both solvent-solvent
interactions and solvent-solute. I am seeing a ~ 3000 kj/mol different
that grows but plateaus out before the system crashes. I would be hesitant
about ascribing the instability to this term as the decrease in LJ and
total energy would be expected at lamba=0.4 and the energy reached its
maximum prior to the system instability.
Have an amazing day,
Scott
On Tue, Jul 30, 2013 at 10:44 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> On 7/30/13 10:17 AM, Scott Pendley wrote:
>
>> 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
>>
>
> Perfect, thanks.
>
>
> 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.
>>
>>
> I appreciate the sentiment, not sure what I've said to give the bad
> impression. Yesterday was great. Hopefully nothing I have said has come
> across as snarky :)
>
> Back to the issue at hand. A couple of things I notice:
>
> nstxtcout = 1000
>>
>
> Given the frequent .xtc output, you should have a pretty good
> visualization of the trajectory. Does it provide any clues?
>
> Pcoupl = parrinello-rahman
>>
>
> You said you did prior NVT and NPT without free energy options. What
> barostat did you use in the previous NPT? Something forgiving like
> Berendsen?
>
> Pcoupltype = semiisotropic
>>
>
> I don't understand the use of semi-isotropic coupling here. That's
> intended for layered systems like membranes. Sounds to me like you've got
> a solute in a uniform solvent, and thus isotropic pressure coupling is
> appropriate.
>
>
> tau_p = 12.0 12.0 ;parrinello-rahman is more stable
>> with larger tau-p, DdJ, 20130422
>>
>
> Very long tau_p, indeed. Normally a value no larger than 5.0 should be
> necessary. That's not necessarily the source of the problem here, but your
> pressure could be running away from you.
>
>
> ; and set the free energy parameters
>> free-energy = yes
>>
>
> Does turning the free energy code off result in a stable simulation? Your
> first message seemed to indicate this would be the case, but I just want to
> clarify.
>
>
> 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
>>
>
> I don't think there's anything wrong with this approach, but I generally
> do not (de)couple vdW and charges simultaneously. Might be worthwhile to
> try a simple vdw -> none transformation to test.
>
>
> 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
>>
>>
> Nothing necessarily wrong here, though you'll get a lot of output that
> likely won't improve your result and will probably just slow the
> simulations down.
>
> What do results at lambda = 0 give you? In theory, this should be a fully
> interacting system and would help track down any bugs in the free energy
> code itself. If lambda = 0 works, then try a non-zero value.
>
> You mentioned before that the LJ terms are affected over time, but based
> on this .mdp file, they shouldn't since you're not doing any sort of
> slow-growth method. Which terms are affected? By how much? Does the time
> at which you see any major change correspond to any observable weirdness in
> the trajectory?
>
> -Justin
>
>
>> 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>
>>>>>> >
>>>>>> <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/**>
>>>>> <h**ttp://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#**<http://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<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>
>>>>> >
>>>>> <htt**p://lists.gromacs.org/**mailman/**listinfo/gmx-users<http://lists.gromacs.org/mailman/**listinfo/gmx-users>
>>>>> <h**ttp://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/**<http://www.gromacs.org/Support/**>
>>>>>
>>>>> Mailing_Lists/Search<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>
>>>>> <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>
>>>>> <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 <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