Previous Thread
Next Thread
Print Thread
Trying to Perform NVT Dynamics
#37830 04/07/20 09:40 PM
Joined: Mar 2005
Posts: 63
R
rossi Offline OP
Forum Member
OP Offline
Forum Member
R
Joined: Mar 2005
Posts: 63
Hello,

I have a problem with NVT dynamics. Something is going wrong, and I am not sure what I am doing incorrectly. At the end of the NVT run, the resulting pdb file has water molecules that are splay away from the RHDO crystl parameters that were provided. This is a simple system of solvated deca-alanine being prepared for a tutorial

! SETUP CRYSTAL (DEFINE, BUILD), IMAGE CENTERING W. MODIFIED PSF
! WE USE THE SAME PARAMS AS IN SOLVATION
set greaterval = 39.4379626 ! The actual value is read from the restart file
crystal define rhdo @greaterval @greaterval @greaterval 60. 90. 60.
crystal build noper 0

IMAGE BYSEG XCEN 0.0 YCEN 0.0 ZCEN 0.0 SELE SEGID HLX END
IMAGE BYRES XCEN 0.0 YCEN 0.0 ZCEN 0.0 SELE RESNAME TIP3 END
IMAGE BYRES XCEN 0.0 YCEN 0.0 ZCEN 0.0 SELE RESNAME SOD END
IMAGE BYRES XCEN 0.0 YCEN 0.0 ZCEN 0.0 SELE RESNAME CLA END

open unit 31 read card name /scratch/anr11010/charmm/alanine/ala10_alpha_helix_production-3.rst !restart file that will be read This file is derived from a fairly good equilibrated NPT dynamics run

! set up nonbond parameters -- same as for heating
nbond inbfrq -1 imgfrq -1 atom vatom cdie eps 1.0 -
elec ewald pmew fftx 48 ffty 48 fftz 48 kappa .34 spline order 6 -
vdw vswitch cutnb 16.0 cutim 16.0 ctofnb 12. ctonnb 10. wmin 1.0

calc pmass = ?stot * 0.02 ! from the tutorial in the text
calc tmass = stot * 0.20 ! from the tutorial in the text

calc pmass = int ( ?stot / 50.0 ) ! from the turorial
calc tmass = @pmass * 10 ! from the tutorial

! Constrain the hlx to center of simulation cell. ! It's a simple deca-alanine solvated and ions added according to the CHARMM tutorial
MMFP
GEO sphere RCM -
xref 0.0 yref 0.0 zref 0.0 -
force 10.0 droff 0.0 -
select segid hlx end
END

! set up shake
shake bonh param sele all end

! Dynamics Command
DYNA LEAP VERLET RESTRT CPT -
NSTEP 100000 TIMESTEP 0.001 -
IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 5000 -
IUNREA 31 IUNWRI -1 IUNCRD -1 IUNVEL -1 KUNIT -1 -
NPRINT 100 NSAVC 1000 NSAVV 0 ISVFRQ 1000 IHBFRQ 0 INBFRQ 25 -
HOOVER PCONS PMASS 0.0 PREF 1.0 REFT 300.0 TMASS @tmass PGAMMA 0 -
IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0 TWINDL 0.0 -
echeck 100.0

WHAT AM I DOING WRONG. THE ABOVE SCRIPT HAS BEEN MOSTLY GLEANED FROM THE CHARMM TUTORIAL PAGES.


As I stated above, the resulting structure is distorted after 100,000 steps. H2O molecules are flying away, although the output says the volume is constant. The pressure doesn't seem good. It's erratic.

Thanks so much for looking at this.

Kind regards,

Angelo

Re: Trying to Perform NVT Dynamics
rossi #37831 04/08/20 01:50 AM
Joined: Sep 2003
Posts: 8,522
Likes: 2
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,522
Likes: 2
I see nothing obvious in the input; perhaps you've missed something in your careful reading of the output log file.

I generally start a new simulation when I switch ensembles.


Rick Venable
computational chemist

Re: Trying to Perform NVT Dynamics
rmv #37948 06/23/20 11:52 AM
Joined: Mar 2005
Posts: 63
R
rossi Offline OP
Forum Member
OP Offline
Forum Member
R
Joined: Mar 2005
Posts: 63
Good Morning.

There was a problem in my input file that led to spurious results. It wasn't any problem with the CHARMM program.

But I wanted to point out the following.

This is what the CHARMM tutorial says about NVT dynamics:

"PCONS: Turns on pressure reporting. The pressure of the system will be held constant unless PMASS is set to 0
(in this case, pressure statistics are still displayed, but constant pressure is not enforced and an NVT ensemble is obtained)."

This is a correct, but albeit somewhat facile statement.

Below is a snippet of a CHARMM script for a solvated decaalanine solvated in a rhombic dodecahedron simulation cell. This was part of a tutorial for a QM and MD workshop presented in May. The decaalanine was constructed as an idealized alpha-helix and obtained from a CHARMM script in one of the earlier versions of CHARMM. The solvation was performed as given in the Full Example.


! NVT Dynamics
dyna leap cpt restart - ! restart the run from our heating restart, use CPT dynamics
timestep 0.001 nstep 1000000 nprint 1000 - ! do 1.0 ns in 1 fs time step
iunrea 31 iunwri 41 iuncrd 51 isvfrq 1000 nsavc 1000 - ! units for reading/writing restarts and coordinate trajectory
hoover reft 300.0 tmass @tmass - ! Hoover thermostat, constant temp - 300 K, tmass as above
pcons pref 1.0 pmass 0.5 pgamma 0 - ! NVT simulation, pmass 0.5 (need pmass~0), collision frequency 0
ihtfrq 0 ieqfrq 0 ntrfrq 1000 - ! temperature under Hoover's control, cancel rotation/translation every 1000 steps
iasors 1 iasvel 1 iscvel 0 ichecw 0 - ! since we're restarting and ihtfrq = ieqfrq = 0, these don't really matter
echeck 100. ! check energy diff as before


The above partial script is similar to that provided in the tutorial, but modified for use in this simulation. Notice that pmass is NOT equal to zero for a perfect NVT simulation. When pmass=0, the system becomes extremely unstable. To approach the value of pmass=0.50 required four separate 1 ns simulations:

1. Step 1
! calculate masses for the pressure and temperature pistons (pmass & tmass)
scalar wmain = mass
scalar mass stat select segid HLX end

calc pmass = ?stot * 0.02
calc tmass = ?stot * 0.20

or

calc pmass = int ( ?stot / 50.0 )
calc tmass = @pmass * 10

2. Step 2, Restart from Step 1
set pmass = 10.0

3. Step 3, Restart from Step 2
set pmass = 1.0

4. Step 4, Restart from Step 3
set pmass = 0.5

With pmass = 0.5, the cell parameters/volume are quite stable for 1 ns, but for a larger system and longer simulation time (~ 50 ns) instabilities may arise.

I apologize for initially thinking that CHARMM has a problem, when it was my input error in the image command.

Regards,

Angelo

Re: Trying to Perform NVT Dynamics
rossi #37949 06/23/20 07:06 PM
Joined: Sep 2003
Posts: 8,522
Likes: 2
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,522
Likes: 2
I've never done anything like that, and never had to.

Using PMASS 0.5 is not NVT, it is NPT with an extremely low piston mass, and could be unstable. That would appear to be an error or misinformation from that tutorial.

For every NVT simulation I have run (and I've run many) I do the following:

  • Equilibrate the system with NPT using a more typical PMASS, such as step 1 above.
  • Calculate the average volume from a stable portion (no major volume changes) of the simulation.
  • Identify and save a coordinate set from the analyzed interval with the volume closest to the average.
  • Start a new simulation, using the saved coordinate set, and using the unit cell sizes corresponding to that coordinate set.
  • If needed, employ a heating and 3-step Verlet equilibration phase instead of a thermostat for the first run, and then enable a thermostat for subsequent runs.


Rick Venable
computational chemist


Moderated by  BRBrooks, lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.4
(Release build 20200307)
Responsive Width:

PHP: 5.6.33-0+deb8u1 Page Time: 0.014s Queries: 22 (0.006s) Memory: 0.9174 MB (Peak: 0.9984 MB) Data Comp: Off Server Time: 2021-01-17 13:21:13 UTC
Valid HTML 5 and Valid CSS