Subject: CHARMM Verlet dynamics

From: [email-masked]
Date: Sun, 28 Jul 2002 15:12:53 -0400 (EDT)

Hi all

I've managed to run 1 nanosecond of Verlet dynamics on my solvated 
peptide.  The temperature output was around 300 K, which is what I 
wanted; unfortunately, the waters flew away!  If I watched a movie of the 
dynamics, it looked like my system blew up.  This is not good.

Here were the settings that I used:

DYNA  VERLET  NSTEP 70000 TIMESTEP 0.001 IHBFRQ 10 -
IPRFRQ 1000 IHTFRQ 0 IEQFRQ 50 NTRFRQ 0 NSAVC 100 ISEED  314159 -
NSAVV 0 FIRSTT 300 FINALT 300 IUNCRD 21 CDIE EPS 1.0 GROUP SWITCH -
IUNWRI 20 ISVFRQ 1000 INBFRQ 20 CUTNB 13.0 CTOFNB 12.0 CTONNB 10.0 -
TWINDL -10.0 TWINDH 10.0 IASORS 0 ISCVEL 1 IASVEL 1 VSWITCH ICHECW 0 
IPRFRQ 1000

I noticed that I had ICHECW set to 0.  I thought this might be the 
problem, so I set it to 1 (in other words, to check to see if the temp. 
falls into the window of 290 to 310 K).  But this didn't help.

So I went to the CHARMM documentation and I found the following example:

"For late equilibration and analysis runs:

DYNAMICS LEAP VERLET RESTART  NSTEP 20000 TIMESTEP 0.001 -
    IPRFRQ 1000 IHTFRQ 2000 IEQFRQ 5000(*) NTRFRQ 5000  -
    IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 -
    NPRINT 100 NSAVC 100 NSAVV 0 IHBFRQ 0 INBFRQ 25  -
    hbond-spec  nonbond-spec   -
    FIRSTT 100.0 FINALT 300.0 TEMINC 100.0   -
    IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 1 TWINDH 10.0 TWINDL -10.0

(*)   Window checking should be disabled for the analysis run (i.e. IEQFRQ=0)
      if you want a real microcanonical ensemble."

Now, this has some heating in it, which I don't want to do.

So I reset my IEQFRQ from 50 to 0, changed ICHECW from 0 to 1, changed 
ISCVEL to 0, and changed NTRFRQ.

This is what my settings are now:

DYNA  VERLET  NSTEP 70000 TIMESTEP 0.001 IHBFRQ 10 -
IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 500 NSAVC 100 ISEED  314159 -
NSAVV 0 TSTRUCT 300 FIRSTT 300 FINALT 300 IUNCRD 21 CDIE EPS 1.0 GROUP SWITCH -
IUNWRI 20 ISVFRQ 1000 INBFRQ 20 CUTNB 13.0 CTOFNB 12.0 CTONNB 10.0 -
TWINDL -10.0 TWINDH 10.0 IASORS 0 ISCVEL 0 IASVEL 1 VSWITCH ICHECW 1 IPRFRQ 1000

The temperature after a few ps is around 145 K or so.

Will it go up?  Are my settings "correct"? Does the fact that I have 
harmonically constrained the ends of my peptide cause the waters to fly away?

I appreciate any help you can offer.

Richard