Re: Problem with equilibration and pressure

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Wed Dec 13 2006 - 07:56:37 CST

Hi Jacob,
I can't guarantee that this is the problem, but you might want to try
the same tests using constraints instead of fixed atoms, with a spring
constant around 5ish. You also might want to try an identically sized
water box on the methyl simulation. Because fixed atoms don't get moved
with the rest of the system during pressure control, they can
occasionally cause pressure equilibration problems (I've only seen this
happen in very large systems, but it is possible that's what is going
on here too).
Peter

Jacob Poehlsgaard wrote:
> Hi all
>
> I have two small RNA systems that I'm trying to equilibrate. The
> difference between them is a single methyl group and a slightly larger
> water box on the "control" without the methyl.
>
> As pr. the the discussion on this list I'm running 250,000 steps with
> pressure control and then 250,000 without to see if things stay stable.
> I'm using a thermostat at 310K all the time.
>
> The problem is that the two systems behave differently. The system
> without the methyl group starts having pressure fluctuations, but they
> seem to be centered around 1 bar, so I'm thinking this is what I want in
> a system without pressure control. The system with the methyl group
> however has fluctuations of the same magnitude, but they're centered
> around 50 bars of pressure, and rapidly find that level after I switch
> off the pressure control.
>
> Control graph: http://www.hardcoil.net/NAMD/control.gif
> Methylated graph: http://www.hardcoil.net/NAMD/methyl.gif
>
> I'm having trouble understanding the difference, surely my parameters
> for the methylated residues can't be THAT bad as to mess with the
> global pressure?
>
> Another thing that puzzles me: If I plot the pressure field, I get a
> HUGE spike 10,000s of bars in the beginning of both runs, even though
> I use restart files from a previous run with pressure control. The
> same happens for the temperature. Is this normal? I thought NAMD was
> supposed to get the velocities and such from the .vel file to avoid
> such spikes in the beginning
>
> Any suggestions or comments are appreciated.
>
> Here's one of the configuration scripts, they are identical in the two
> runs, apart from filenames:
>
> #############################################################
> ## JOB DESCRIPTION ##
> #############################################################
>
> # Equilibration run
> # Uses restart files from _rel run
> # Restraints on terminals
> # Constant temp
> # constant pressure for first half
>
>
> #############################################################
> ## ADJUSTABLE PARAMETERS ##
> #############################################################
>
> set basename capsite_c
>
> structure ${basename}_i.psf
> coordinates ${basename}_iw_restraints.pdb
> #coords are ignored as we use restart files. They just have to be there.
>
> binCoordinates output/3_relaxation/${basename}_rel.coor
> binVelocities output/3_relaxation/${basename}_rel.vel
> extendedSystem output/3_relaxation/${basename}_rel.xsc
>
> firsttimestep 0
> set temperature 310
>
> #############################################################
> ## SIMULATION PARAMETERS ##
> #############################################################
>
> # Input
> paraTypeCharmm on
> parameters /people/disk2/jacobp/MD/common/par_all27_prot_na.prm
> parameters /people/disk2/jacobp/MD/common/par_all27_modified_RNA.prm
>
> wrapWater on
> wrapAll on
>
> ##############################################
> #PME (for full-system periodic electrostatics)
> PME yes
> PMEGridSizeX 72
> PMEGridSizeY 60
> PMEGridSizeZ 72
>
> # Force-Field Parameters
> exclude scaled1-4
> 1-4scaling 1.0
> cutoff 12.
> switching on
> switchdist 10.
> pairlistdist 13.5
>
> # Integrator Parameters
> timestep 1.0 ;# 1fs/step
> rigidBonds none
> nonbondedFreq 1
> fullElectFrequency 1
> stepspercycle 20
>
> # Output
> outputName /scratch/${basename}_equil
> restartfreq 2000
> dcdfreq 5000
> xstFreq 5000
> outputEnergies 5000
> outputPressure 5000
>
> #############################################################
> ## Enemble Controls ##
> #############################################################
>
> ##############################
> # Constant Temperature Control
> langevin on ;# do langevin dynamics
> langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
> langevinTemp $temperature
> langevinHydrogen no ;# don't couple langevin bath to hydrogens
>
> ##############################
> # Constant pressure
> useGroupPressure no ;
> useFlexibleCell no ;# no for water box, yes for membrane
> useConstantArea no ;# no for water box, yes for membrane
> langevinPiston on
> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
> langevinPistonPeriod 100.
> langevinPistonDecay 50.
> langevinPistonTemp $temperature
>
> #############################################################
> ## Fixing and restraining ##
> #############################################################
> ## We're fixing the terminals using same coordinats as before
>
>
> # Fixed Atoms Constraint (set PDB occupancy-column to 1)
> # Fix terminals during entire run
> if {1} {
> fixedAtoms on
> fixedAtomsFile ${basename}_rel_restraints.pdb
> fixedAtomsCol O
> }
>
> #############################################################
> ## EXECUTION SCRIPT ##
> #############################################################
>
> #Run
> run 250000
> langevinPiston off
> run 250000
>
> Jacob Poehlsgaard
> Ph.D. Student, BMB
> SDU, Denmark

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:42:55 CST