From: Marc Q. Ma (qma_at_oak.njit.edu)
Date: Wed Apr 13 2005 - 08:59:31 CDT
Dear Lina,
Did you only minimize 2000 steps? That maybe insufficient. To find out,
you can plot the total energy to see if it is still dropping. If it is,
then you need extra steps. I would use at least 60000 steps for
minimization. Make sure you have a fully minimized system.
Also, as you noticed, turning on rigidbonds will also help with
stability. It makes sense to turn on rigidbonds in biomolecular
simulations.
Other parameters that you can fine tune include the cutoff distance and
switchon distance. Since you are using PME and MTS integrator, which
take care of long range interactions periodically, you do not need a
cutoff distance of 12 AA. A longer cutoff value makes it more time
consuming to compute the real part of the electrostatics forces in PME
calculations. I would make it 6.5 to 8. The swichdist should be smaller
than the cutoff value. If you use 6.5 as cutoff, you may use 2-4.5 as
switchdist -- which is used in the smoothing function for nonbonded
medium range force calculations.
If it still does not work, make your timestep .5 or even .1. This put
you on a super conservative safe side regarding stability of
simulations. Of course this would result in more computer time to reach
desired simulation length, the numerical artificat is much much
smaller. When you get a fully thermalized and equilibrated system, you
may use a timestep of 1 for long production runs.
The last thing I would try is to double your PMEGridSize[X,Y,Z], which
can reduce the numerical error related to PME calculations.
Let us know how it goes. Cheers!
Marc
On Apr 13, 2005, at 4:16 AM, Lina Nilsson wrote:
>
>
> Dear all,
> I need to use AMBER force fields in NAMD, but from .crd and .top files
> created with xleap, I can't get a system that is stable in my NAMD
> simulation (tried with a protein I have previously simulated fine
> using Charmm, and with a two amino acid 'test protein.') Do I need to
> do extra minimzation steps when using AMBER? Other ideas????
>
> Details:
> My periodic boxes (water and protein) blow up as soon as I begin to
> increase the temperature above 0 K in NAMD (the energies, volumes etc
> seem to be fine/plateau during mimimization in NAMD, as far as I can
> tell). With the temp turned on, the system volume shrinks rapidly,
> only to then start increasing and then 'blowing up' as atoms start
> moving too fast for the simulation to handle (the temperaure is
> increasing rapidly too, of course).
>
> The NAMD base scripts I use work fine for charmm ff based simulations.
> The densities I get from the xleap prep are low (<0.8 g/cc)compared to
> what I get solvating my charmm systems. Reducing the closeness factor
> in xleap to 0.5 doesn't help. May the low density be causing
> instabilities? Turning on RigidBonds reduces the problem, but does not
> get rid of it, for what that is worth. I am completely stuck on this
> and ANY input is much appreciated!
>
> Lina
>
>
>
> TWO AMINO ACID TEST CASE
>
> XLEAP:
>
> xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff # tried ff03 and then
> ff94
>
> filename = loadpdb PDB.pdb
> solvatebox filename TIP3PBOX 16
> charge filename # --> no charge
> check filename
> # --> no errors. "Unit is ok" Some hydrogens are too close to each
> other
> saveamberparm filename filename.top filename.crd
> savepdb filename filename.pdb
>
> # Other than the water density, the system looks fine in VMD (no
> missing atoms etc).
>
> THEN MINIMIZING IN NAMD in 3 steps (first water only, then holding
> protein backbone still, then all atom minimization):
>
> cwd ./
> amber on
> readexclusions no
> parmfile ff94_xxx.top
> coordinates ff94_xxx.pdb
>
> restartname xxx.rst
> restartfreq 100
> binaryrestart yes
> outputname xxx
> binaryoutput no
> DCDfile xxx.dcd
> DCDfreq 100
> outputEnergies 1
> outputMomenta 1
> outputPressure 1
> XSTfile xxx.xst
> XSTfreq 100
>
> numsteps 2000
> timestep 1.0
>
> fixedAtoms on
> fixedAtomsForces on
> fixedAtomsFile fixedAtoms.pdb
> fixedAtomsCol B
>
> stepspercycle 4
> nonbondedFreq 1
> cutoff 12.0
> switching on
> switchdist 10.0
> pairlistdist 13.0
> exclude scaled1-4
>
> 1-4scaling 0.833333
> fullElectFrequency 4
> PME yes
> PMETolerance 0.000001
> PMEInterpOrder 4
> PMEGridSizeX 36
> PMEGridSizeY 36
> PMEGridSizeZ 36
>
> scnb 2
> temperature 0
> dielectric 1.0
>
> cellBasisVector1 34.572 0 0
> cellBasisVector2 0 33.2880 0
> cellBasisVector3 0 0 33.2840
> cellOrigin -0.1518 0.1325 -0.0739
> wrapWater on
>
> constraints off
> minimization on
>
>
> --> output *eems* ok...
>
> FINALLY, THERMALIZING:
>
>
> cwd ./
> amber on
> readexclusions no
> parmfile ff94_xxx.top
> coordinates minoutput.coor
> extendedSystem minoutput.xsc
>
> restartname xxx.rst
> restartfreq 100
> binaryrestart yes
> outputname xxx
> binaryoutput no
> DCDfile xxx.dcd
> DCDfreq 100
> outputEnergies 100
> outputMomenta 100
> outputPressure 100
> XSTfile xxx.xst
> XSTfreq 100
>
> numsteps 40000
> timestep 1.0
>
> fixedAtoms on
> fixedAtomsForces on
> fixedAtomsFile fixedAtoms.pdb
> fixedAtomsCol B
>
> stepspercycle 4
> nonbondedFreq 1
> cutoff 12.0
> switching on
> switchdist 10.0
> pairlistdist 13.0
> exclude scaled1-4
>
> 1-4scaling 0.833333
> fullElectFrequency 4
> PME yes
> PMETolerance 0.000001
> PMEInterpOrder 4
> PMEGridSizeX 36
> PMEGridSizeY 36
> PMEGridSizeZ 36
>
> scnb 2
> temperature 0
> dielectric 1.0
>
> reassignFreq 1000
> reassignTemp 10.0
> reassignIncr 10.0
> reassignHold 310.0
>
> useGroupPressure no
> useFlexibleCell no
> BerendsenPressure on
> BerendsenPressureTarget 1.01325
> BerendsenPressureCompressibility 0.000045
> BerendsenPressureRelaxationTime 1000
> BerendsenPressureFreq 4
> wrapWater on
> constraints off
>
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:39:20 CST