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 <jalem...@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 <jalem...@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@outerbanks.umaryland.****edu <jalemkul@outerbanks.** >>> umaryland.edu <jalem...@outerbanks.umaryland.edu>> | (410) >>> 706-7441 >>> >>> ==============================****==================== >>> >>> -- >>> gmx-users mailing list gmx-users@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-requ...@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@outerbanks.umaryland.**edu <jalem...@outerbanks.umaryland.edu> | > (410) > 706-7441 > > ==============================**==================== > -- > gmx-users mailing list gmx-users@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-requ...@gromacs.org. > * Can't post? Read > http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists> > -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at 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-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
