Problem with equilibration and pressure

From: Jacob Poehlsgaard (mail_lists_at_hardcoil.net)
Date: Wed Dec 13 2006 - 05:31:49 CST

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 - 05:19:57 CST