Re: membrane equilibration

From: Jeffrey Potoff (jpotoff_at_chem1.eng.wayne.edu)
Date: Thu Mar 03 2011 - 19:58:11 CST

Hi Jeremy,
What we do to get these things going is use a 0.1-0.5 fs time step in
the canonical ensemble and a slow temperature ramp from 5 K to wherever
we need to get to. The total simulation is about 200000 steps; although
many times we can get away with a faster temperature ramp. After that,
we can run at any temperature using a 2fs timestep and the system is
completely stable. In your case, if you have a structure that was
already stable in Charmm, you might just need a 10000 steps or so at 0.1
fs to bump it out of the rattle error.

Here is an input file I used to get a DMPC bilayer started (initial
coordinates were from CHARMM-GUI):

#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################
structure start.psf
coordinates start.pdb

set TEMP 10.0
reassignFreq 10000
reassignTemp 10
reassignIncr 10
reassignHold 220

set outputname heat
#set inputname heat
#binCoordinates $inputname.restart.coor
#binVelocities $inputname.restart.vel ;# remove the "temperature"
entry if you use this!
#extendedSystem $inputname.restart.xsc
firsttimestep 0

#############################################################
## SIMULATION PARAMETERS ##
#############################################################

# Input
paraTypeCharmm on
parameters par_all27_prot_lipid.inp
temperature $TEMP

# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 14
switching on
switchdist 13.0
pairlistdist 16.0
margin 3.0

# Integrator Parameters
timestep 0.1 ;# 2fs/step
rigidBonds none ;# needed for 2fs steps
nonbondedFreq 1
fullElectFrequency 2
stepspercycle 10

# Constant Temperature Control
langevin on ;# do langevin dynamics
langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
langevinTemp $TEMP
langevinHydrogen on ;# don't couple langevin bath to hydrogens

# Periodic Boundary Conditions
cellBasisVector1 64.2 0. 0.0
cellBasisVector2 0. 64.2 0.0
cellBasisVector3 0.0 0. 64.2
cellOrigin 0.0 0.0 0.0

wrapAll on
# PME (for full-system periodic electrostatics)
PME yes
PMEGridSizeX 64
PMEGridSizeY 64
PMEGridSizeZ 64

# Constant Pressure Control
(vpar_all27_prot_lipid.inppar_all27_prot_lipid.inpariable volume)
useGroupPressure no ;# needed for rigidBonds
useFlexibleCell no
useConstantArea no

langevinPiston off
#langevinPistonTarget 1.0 ;# in bar -> 1 atm
#langevinPistonPeriod 100.
#langevinPistonDecay 50.
#langevinPistonTemp $TEMP

# Output
outputName $outputname

restartfreq 5000 ;# 500steps = every 1ps
dcdfreq 5000
xstFreq 5000
outputEnergies 5000
outputPressure 5000

#############################################################
## EXECUTION SCRIPT ##
#############################################################

# Minimization
   minimize 1000

for { set TEMP 10 } { $TEMP < 220 } { incr TEMP 10 } {
run 10000
reinitvels $TEMP
langevinTemp $TEMP
}

On 3/3/2011 7:30 PM, jeremy adler wrote:
> Hi Dir. Harrison,
>
> these are the last few lines of minimization (1000 steps with
> virtually no change in the energies. If the equilibration in charmm
> was carried out over about 1500ps behaved itslef perfectly well in the
> subsequent production run.. As potentially another solution, do you
> know of a script that will convert charmm restart files into namd
> restart files?
> thanks!
>
> Jeremy
>
> starting system:
> http://www.mediafire.com/file/iq1q1z96ac79nje/step6.6_equilibration.pdb
>
> BRACKET: 2.85511e-09 31.25 -9.12482e+11 -9.07398e+11 -8.99207e+11
> PRESSURE: 998 3.88841e+16 -2.14232e+16 1.12373e+17 3.04533e+16
> -2.00369e+16 2.26124e+16 7.48332e+15 -6.20386e+15 1.91964e+16
> GPRESSURE: 998 3.88841e+16 -2.14232e+16 1.12373e+17 3.04533e+16
> -2.00369e+16 2.26124e+16 7.48332e+15 -6.20386e+15 1.91964e+16
> ENERGY: 998 2346753.6788 468115.1150 4454.7007
> 4848.7842 -31018.7247 99999999.9999 0.0000
> 0.0000 0.0000 99999999.9999 0.0000
> 99999999.9999 99999999.9999 0.0000 99999999.9999
> 99999999.9999 23040.0000 99999999.9999 99999999.9999
>
> BRACKET: 1.76455e-09 16.1875 -9.12482e+11 -9.07398e+11 -9.04265e+11
> PRESSURE: 999 3.88841e+16 -2.14232e+16 1.12373e+17 3.04533e+16
> -2.00369e+16 2.26124e+16 7.48332e+15 -6.20386e+15 1.91964e+16
> GPRESSURE: 999 3.88841e+16 -2.14232e+16 1.12373e+17 3.04533e+16
> -2.00369e+16 2.26124e+16 7.48332e+15 -6.20386e+15 1.91964e+16
> ENERGY: 999 2346695.8646 468093.2401 4454.6479
> 4848.7877 -31019.1667 99999999.9999 0.0000
> 0.0000 0.0000 99999999.9999 0.0000
> 99999999.9999 99999999.9999 0.0000 99999999.9999
> 99999999.9999 23040.0000 99999999.9999 99999999.9999
>
> BRACKET: 1.09055e-09 3.0625 -9.09338e+11 -9.07398e+11 -9.04265e+11
> NEW SEARCH DIRECTION
> INITIAL STEP: 2.5e-07
> GRADIENT TOLERANCE: 1.89282e+08
> TIMING: 1000 CPU: 2751.05, 2.03812/step Wall: 2751.05, 2.03812/step,
> 0 hours remaining, 65648 kB of memory in use.
> PRESSURE: 1000 3.88841e+16 -2.14232e+16 1.12373e+17 3.04533e+16
> -2.00369e+16 2.26124e+16 7.48332e+15 -6.20386e+15 1.91964e+16
> GPRESSURE: 1000 3.88841e+16 -2.14232e+16 1.12373e+17 3.04533e+16
> -2.00369e+16 2.26124e+16 7.48332e+15 -6.20386e+15 1.91964e+16
> ETITLE: TS BOND ANGLE DIHED
> IMPRP ELECT VDW BOUNDARY
> MISC KINETIC TOTAL TEMP
> TOTAL2 TOTAL3 TEMPAVG PRESSURE
> GPRESSURE VOLUME PRESSAVG GPRESSAVG
>
> ENERGY: 1000 2346695.8646 468093.2401 4454.6479
> 4848.7877 -31019.1667 99999999.9999 0.0000
> 0.0000 0.0000 99999999.9999 0.0000
> 99999999.9999 99999999.9999 0.0000 99999999.9999
> 99999999.9999 23040.0000 99999999.9999 99999999.9999
>
> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 1000
> OPENING COORDINATE DCD FILE
> WRITING COORDINATES TO DCD FILE AT STEP 1000
> WRITING COORDINATES TO RESTART FILE AT STEP 1000
> FINISHED WRITING RESTART COORDINATES
> WRITING VELOCITIES TO RESTART FILE AT STEP 1000
> FINISHED WRITING RESTART VELOCITIES
> REINITIALIZING VELOCITIES AT STEP 1000 TO 300 KELVIN.
> TCL: Running for 250000 steps
> ERROR: Constraint failure in RATTLE algorithm for atom 3241!
> ERROR: Constraint failure; simulation has become unstable.
> ERROR: Constraint failure in RATTLE algorithm for atom 5379!
> ERROR: Constraint failure; simulation has become unstable.
> ERROR: Exiting prematurely.
> ==========================================
> WallClock: 2760.421875 CPUTime: 2760.421875 Memory: 65648 kB
> End of program
>
> On Thu, Mar 3, 2011 at 7:10 PM, Chris Harrison <charris5_at_gmail.com
> <mailto:charris5_at_gmail.com>> wrote:
>
> Were you using multiple-time-stepping in CHARMM? I suspect not, and
> that you are in NAMD. If this is the case, you could have some
> mildly pathological interactions that are not updated as often
> and lead to strong velocities that result in RATTLE failures. Some
> solutions are to either determine the specific atom (via it's
> index) and
> resolve whatever may be causing the instability, or you could attempt
> longer minimization/equilibration runs in hopes the unstable
> interaction
> resolves itself. If this is indeed your situation, you'll be more
> efficient if you re-run the crashing simulation, but outputting the
> frames and energies every step, and then analyzing that trajectory for
> the source of your RATTLE failure.
>
> HOWEVER, all the above is just a best guess based on your mention of a
> RATTLE failure. Axel is right that it's difficult to diagnose the
> problem
> without seeing the actual output in your query. It's often useful to
> cut and paste 20 or so lines around the error message, so others
> can see
> exactly the output. With that kind of data the namd-list can be
> much more helpful.
>
> Best,
> Chris
>
>
> --
> Chris Harrison, Ph.D.
> Theoretical and Computational Biophysics Group
> NIH Resource for Macromolecular Modeling and Bioinformatics
> Beckman Institute for Advanced Science and Technology
> University of Illinois, 405 N. Mathews Ave., Urbana, IL 61801
>
> char_at_ks.uiuc.edu <mailto:char_at_ks.uiuc.edu>
> Voice: 217-244-1733
> http://www.ks.uiuc.edu/~char <http://www.ks.uiuc.edu/%7Echar>
> Fax: 217-244-6078
>
>
> jeremy adler <jeremyeadler_at_gmail.com
> <mailto:jeremyeadler_at_gmail.com>> writes:
> > Date: Thu, 3 Mar 2011 18:03:18 -0500
> > From: jeremy adler <jeremyeadler_at_gmail.com
> <mailto:jeremyeadler_at_gmail.com>>
> > To: Axel Kohlmeyer <akohlmey_at_gmail.com <mailto:akohlmey_at_gmail.com>>
> > Cc: namd-l_at_ks.uiuc.edu <mailto:namd-l_at_ks.uiuc.edu>
> > Subject: Re: namd-l: membrane equilibration
> >
> > Hi Dr Kohlmeyer,
> >
> > It seems as though the energies never come down and I get
> "ERROR: Constraint
> > failure in RATTLE algorithm" when the program exits the minimization
> > portion. Ive tried using both fixed atom as well as harmonic
> restraints with
> > the same basic result. I had seen somewhere in the mailing list that
> > periodic boundary conditions might be too small so i made them
> bigger only
> > to have the water collapse into the membrane. Im not exactly
> clear why the
> > energies are so he given they werent nearly as high when i was
> using charmm
> > and it equilibrated fine then (i am using the preequilbrated
> membrane from a
> > charmm run i conducted). I would like to use namd because i have
> access to a
> > bluegene computer and namd functions much better on bluegene
> than charmm.
> >
> >
> > thanks for your help,
> >
> > Jeremy
> >
> > On Thu, Mar 3, 2011 at 2:43 PM, Axel Kohlmeyer
> <akohlmey_at_gmail.com <mailto:akohlmey_at_gmail.com>> wrote:
> >
> > > On Thu, Mar 3, 2011 at 2:02 PM, jeremy adler
> <jeremyeadler_at_gmail.com <mailto:jeremyeadler_at_gmail.com>>
> > > wrote:
> > > > Hi
> > > >
> > > > I was wondering if anyone has a membrane equilbration
> protocol that
> > > works. I
> > > > tried adapting the one in the tutorial with a membrane
> previousely
> > > > equilibrated in charmm to no avail
> > >
> > > jeremy,
> > >
> > > before stating that something doesn't work, you should explain
> > > how it fails for you and why you think that this is not adequate.
> > >
> > > in molecular simulations, there are so many little things that
> > > can matter, it is rare that there is a one-size-fits-all solution.
> > >
> > > cheers,
> > > axel.
> > >
> > > >
> > > > thanks
> > > >
> > > > Jeremy
> > > >
> > >
> > >
> > >
> > > --
> > > Dr. Axel Kohlmeyer
> > > akohlmey_at_gmail.com <mailto:akohlmey_at_gmail.com> http://goo.gl/1wk0
> > >
> > > Institute for Computational Molecular Science
> > > Temple University, Philadelphia PA, USA.
> > >
>
>

-- 
======================================================================
Jeffrey J. Potoff			  jpotoff_at_chem1.eng.wayne.edu
Associate Professor			  Wayne State University		
Department of Chemical Engineering and Materials Science
5050 Anthony Wayne Dr			  Phone:(313)577-9357		
Detroit, MI 48202  	                  Fax:  (313)578-5815
http://potoff1.eng.wayne.edu
======================================================================

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:19:53 CST