Re: Rattle Error During Equilibration of Lipid Bilayer Membrane

From: Dan Strahs (dstrahs_at_pace.edu)
Date: Mon Feb 13 2006 - 11:30:35 CST

Hi:

Two things: I understand that DMPC has approximately the same density
as water; thus, 11,500 waters (or 34,500) requires a cube with a ~70.1
side or 3.44x10^5 A^3. Your system has nearly 38,900 atoms (based on the
200 DMPC/5100 waters); thus, I would expect you need ~3.88x10^5 A^3 or a
cube of side ~72.9. I think your cell may be too small, and NAMD can't
equilibrate the size to escape the inevitable RATTLE error. But this is
all based on comparing your cell to a water cell - I have no experience
with DMPC.

If this isn't helping, I would suggest dumping coordinates at a high
frequency and using VMD to watch the system, particularly the water that
gives the RATTLE error. In the past, I successfully tracked errors in
CHARMM that caused SHAKE errors using this protocol.

Dan Strahs

On Sun, 12 Feb 2006, Arneh Babakhani wrote:

> Hi Morad,
>
> Yep, all of the below is correct. Still no good, still aborting. I
> even went back and did an additional 10,000 steps of minimization. Then
> when I tried to heat it up again, it gradually heated up to 330K, stayed
> constant at the temperature for a while (indicating that heating was
> complete), then aborted again.
>
> Thanks again,
>
> Arneh
>
> Morad Alawneh wrote:
>
> > *Hi Arneh,
> >
> > In this case try do what I did which helped me to get rid of this error:
> >
> > 1) make sure your box is big enough; e.g. let x length be larger by 2A
> > from your atoms for each side, and do it for the other axises
> >
> > 2) make sure these parameters are the same in all your configuration
> > files: switching, switchDist, and cutoff
> >
> > 3) check the following
> > timestep = 2.0
> > fullElectFrequency = 2
> > nonbondedFreq = 1
> >
> > Note: **fullElectFrequency * **nonbondedFreq <= 4.0
> >
> > That is my way, I hope it will work for you
> > *
> >
> > *
> > *
> >
> >
> > /*Morad Alawneh*/
> >
> > *Department of Chemistry and Biochemistry*
> >
> > *C100 BNSN, BYU*
> >
> > *Provo, UT 84602*
> >
> >
> >
> > Arneh Babakhani wrote:
> >
> >> Hi Morad, thank you for your suggestion.
> >>
> >> I didn't explicitly specify the basis vectors; I used the previous
> >> dimensions, from the xsc file. Those dimensions are about
> >> ~71x71x71. I took your suggestion of increasing the PME gridsize to
> >> be >=1.5 times that length(~110x110x110). Indeed, that did help; the
> >> simulation ran for much longer than it did before, but it still
> >> aborted prematurely, with the same error. So although that might have
> >> been part of the problem, it appears there's still something else
> >> going on here.
> >>
> >> In any case, thank you for your help. Would greatly appreciate any
> >> further ideas.
> >>
> >> Regards,
> >>
> >> Arneh
> >>
> >>
> >>
> >> Morad Alawneh wrote:
> >>
> >>> *Hi,
> >>>
> >>> It seems to me that you don't have periodic boundary for your
> >>> system, and once you build that you have to make PMEGridSize to be
> >>> >= 1.5 of the length of that axis.
> >>>
> >>> I hope that will fix your problem.
> >>> *
> >>>
> >>>
> >>>
> >>> /*Morad Alawneh*/
> >>>
> >>> *Department of Chemistry and Biochemistry*
> >>>
> >>> *C100 BNSN, BYU*
> >>>
> >>> *Provo, UT 84602*
> >>>
> >>>
> >>>
> >>> Arneh Babakhani wrote:
> >>>
> >>>> Dear NAMD community,
> >>>>
> >>>> I constructed and solvated a lipid bilayer membrane. This consists
> >>>> of 200 DMPC lipids (100 per monolayer), about 5100 water
> >>>> molecules, for a total of ~40,000 atoms.
> >>>>
> >>>> I conducted 10,000 steps of minimization, ultimately achieving a
> >>>> gradient tolerance of ~2.91. At the end of this minimization, I
> >>>> zeroed the velocities of all atoms. No problem here . . .
> >>>>
> >>>> I then gradually heated up my system, over 2500 steps (2fs
> >>>> timestep), to 330K. Looking at the plot of the temperature, the
> >>>> system does heat up gradually, eventually plateauing out at 330K,
> >>>> towards the end of the run. No problem here . . .
> >>>>
> >>>> I then turn on constant pressure and attempt to run for another
> >>>> 2500 steps. The temperature spikes from 330K to ~375K, by the
> >>>> second step! (I am implementing langevin temp control). The temp.
> >>>> gradually decreases through the course of the simulation, but it
> >>>> never quite gets back down to 330K. NAMD prematurely aborts
> >>>> halfway through the run, at about 350K, giving the following error:
> >>>>
> >>>> ERROR: Constraint failure in RATTLE algorithm for atom 13561!
> >>>>
> >>>> ERROR: Constraint failure; simulation has become unstable.
> >>>>
> >>>> ERROR: Exiting prematurely.
> >>>>
> >>>>
> >>>> I found through some trial & error that by increasing the langevin
> >>>> Damping parameter, I can decrease the magnitude of this initial
> >>>> temperature spike (although I can't get rid of it completely, very
> >>>> strange). By decreasing this spike, the system comes back down to
> >>>> 330K much faster, and the simulation proceeds to completion (with
> >>>> langevin Damping = 7).
> >>>>
> >>>> I then try to run a 1 ns equilibration of my membrane, again using
> >>>> lang. Damping = 7 (since this seemed to work for the constant
> >>>> pressure run). This runs for about 850 ps and then aborts, with
> >>>> the same RATTLE algorithm error. However, during this
> >>>> equilibration, the temp. stays constant, as does the pressure (no
> >>>> spikes), and nothing else is observably wrong. I repeated the
> >>>> attempt for lang. Damping =5, and 9, no joy, still aborts.
> >>>>
> >>>> Would greatly appreciate any ideas? (I've pasted my conf file
> >>>> below for reference)
> >>>>
> >>>> Regards,
> >>>>
> >>>> Arneh
> >>>>
> >>>> #############################################################
> >>>> # 1NS equilibration dmpc lipid bilayer in a Water Box
> >>>> #############################################################
> >>>> structure ../common/solvated.psf
> >>>> coordinates ../common/solvated.pdb
> >>>> bincoordinates lipwaterconstpress2.coor
> >>>> binvelocities lipwaterconstpress2.vel
> >>>> extendedsystem lipwaterconstpress2.xsc
> >>>>
> >>>> set temperature 330
> >>>> set outputname lipwaterequil1
> >>>>
> >>>> firsttimestep 0
> >>>>
> >>>> ## SIMULATION PARAMETERS ##
> >>>>
> >>>> paraTypeCharmm on
> >>>> parameters ../common/par_all27_prot_lipid.inp
> >>>>
> >>>> # Force-Field Parameters
> >>>> exclude scaled1-4
> >>>> 1-4scaling 1.0
> >>>> cutoff 10.
> >>>> switching on
> >>>> switchdist 9.
> >>>> pairlistdist 11.5
> >>>>
> >>>>
> >>>> # Integrator Parameters
> >>>> timestep 2.0 ;# 2fs/step
> >>>> rigidBonds all ;# needed for 2fs steps
> >>>> nonbondedFreq 1
> >>>> fullElectFrequency 2
> >>>> stepspercycle 10
> >>>>
> >>>>
> >>>> # Constant Temperature Control
> >>>> langevin on ;# do langevin dynamics
> >>>> langevinDamping 7 ;# damping coefficient (gamma) of 9/ps
> >>>> langevinTemp $temperature
> >>>> langevinHydrogen off ;# don't couple langevin bath to hydrogens
> >>>>
> >>>> wrapAll on
> >>>>
> >>>> # PME (for full-system periodic electrostatics)
> >>>> PME yes
> >>>> PMEGridSizeX 64
> >>>> PMEGridSizeY 64
> >>>> PMEGridSizeZ 64
> >>>>
> >>>> # Constant Pressure Control (variable volume)
> >>>> useGroupPressure yes ;# needed for rigidBonds
> >>>> useFlexibleCell no
> >>>> useConstantArea no
> >>>>
> >>>> langevinPiston on
> >>>> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
> >>>> langevinPistonPeriod 100.
> >>>> langevinPistonDecay 200.
> >>>> langevinPistonTemp $temperature
> >>>>
> >>>> # Output
> >>>> outputName $outputname
> >>>> restartfreq 500 ;# 500steps = every 1ps
> >>>> dcdfreq 500
> >>>> outputEnergies 100
> >>>> outputPressure 100
> >>>>
> >>>> run 500000 ; # 1ns
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>
> >>
>
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:41:37 CST