CHARMM Development Project
Posted By: lennart run-cpt-md.inp - 10/23/03 10:48 AM
*FILENAME: run-cpt-md.inp
*PURPOSE: setup and run Constant Pressure&Temperature (rectangular box) MD simulation
*AUTHOR: Lennart Nilsson, Karolinska Institutet (October 8, 2003)
*
! Standalone CHARMM input script. PBC MD of protein in box of water
!
! Simple setup with the solute in a box of user specified size
! This file is used both for the initial setup and for the production simulation,
! based on command line argument N (=0, setup; >0, continued simulation)
! Solvate the protein, setup PBC, run simulation.
! Uses unix environment variable CHM_HOME to find root of charmm installation
! and CHM_STREAM to find these scripts
! Uses stream files solvent-box.str
! Reads psf and coordinates from existing files (from gen-prot.inp)
!ARGUMENT: N= run number; N=0 setup and initial MD; N>0, continue with piece N
! 50 ps initial heating followed by 100ps pieces
! charmm < run-cpt-md.inp > run-cpt-md.inp N=0 (or 1,2,....)
! Saves coordinate trajectory, restartfile and snapshot at end of traj. piece
!CONSTANT PRESSURE REFERENCE: Feller, S. E., Zhang, Y., Pastor, R. W., and Brooks, B. R. (1995).
!Constant pressure molecular dynamics simulation: The Langevin piston method.
!J Chem Phys 103, 4613-4621

! get aa definitions and parameters
read rtf card name $CHM_HOME/toppar/top_all22_prot.inp
read para card name $CHM_HOME/toppar/par_all22_prot.inp

! name of the system, used in filenames
set FILE 4pti
if @N .gt. 0 goto get_psf

! Get psf and coordinates, from gen-prot.inp
read psf card name @FILE.psf
read coor card name @FILE.crd

! find out where protein is and put it with center at origin
coor stat sele segi bpti end
! orient molecule with moments of inertia along x,y,z-axes
! everything will be translated/rotated, but the bpti segment determines by how much
coor orie sele segi bpti end
coor stat sele segi bpti end
! Add water; here no ions are added, this should perhaps also be done
! Box size could be autmatically determined here by adding say 8A on each side to the
! dimensions of the protein, we also make the sizes integer numbers (not necessary)
calc XSIZ = INT( ?XMAX - ?XMIN + 16)
calc YSIZ = INT( ?YMAX - ?YMIN + 16)
calc ZSIZ = INT( ?ZMAX - ?ZMIN + 16)
! Let the crystal waters be as they are in segment XWAT
define SOLUTE sele all end
stream $CHM_STREAM/solvent-box.str
! save coordinates and psf
write psf card name @FILE-solv.psf
* bpti with crystal waters solvated with TIP3 waters, setup from top_all22_prot
*
write coor card name @FILE-solv.crd
* bpti with crystal waters solvated with TIP3 waters, setup from top_all22_prot
* in rectangular box: @XSIZ x @YSIZ x @ZSIZ
*
! save size to file so we can reuse it
write title name @FILE-box-size.str
** PBC box definition for @FILE simulation
**
* set XSIZ @XSIZ
* set YSIZ @YSIZ
* set ZSIZ @ZSIZ
* return
*
! setup Periodic Boundary Conditions
! Here we have to use the crystal facility
! Variables XSIZ,YSIZ, and ZSIZ have already been defined
crystal define orthorombic @XSIZ @YSIZ @ZSIZ 90.0 90.0 90.0
crystal build cutoff 14.0

! also specify that water molecules should be recentered as molecules (residues) if they
! get out of primary box, but that the protein molecule is a whole segment
image byresidue sele resn tip3 end xcen 0.0 ycen 0.0 zcen 0.0
image bysegment sele segi bpti end xcen 0.0 ycen 0.0 zcen 0.0

! We need to keep the water molecules rigid, and also want to
! keep covalent X-H bonds fixed, so use SHAKE
shake bonH para ! take reference values for bond lengths from param file
!!!!! MINIMIZATION
! allow water to adjust around protein
cons harm force 50.0 sele segid bpti end
minimize sd nstep 50 -
cdie eps 1.0 fshift vshift cutnb 14.0 ctofnb 12.0 cutim 14.0
cons harm force 20.0 sele segid bpti end
minimze abnr nstep 50
write coor card name @FILE_wmin.crd
* gently mininmized coordinates of 4pti in @XSIZ x @YSIZ x @ZSIZ box of water
*
!!!!! INITIAL DYNAMICS
! turn off harmonic restraints
cons harm force 0.0 sele ALL end

open unit 12 write form name @FILE-0.res
open unit 21 write unform name @FILE-0.cor
open unit 22 write form name @FILE-0.mon

! time to spend on heating may vary.
! 100 to 300 K in a few ps is probably typical
! and then the time for equilibration will have to be judged
! from case to case (5-50ps, or more)
! Use extended system pressure algorithm (see pressure.doc)
! with simple temperature check using veolicty scaling.
! For more precise temperature control Hoover's method may be usd
dynamics start CPT time 0.002 nstep 25000 iseed 31415 -
firstt 48.0 finalt 298.0 teminc 10.0 ihtfrq 500 -
PCONst pmass 400.0 pgamma 20.0 tbath 300.0 PREFerence 1.0 -
ieqfrq 2000 ichecw 1 twindl -5.0 twindh +5.0 iasors 0 -
nprint 100 iprfrq 500 -
atom cdie fshift vshift cutnb 14.0 ctofnb 12.0 -
inbfrq -1 ihbfrq 0 imgfrq -1 - ! -1 for inbfrq/imgfrq means heuristic update
iunwrit 12 iuncrd 21 kunit 22 nsavc 100 !ie update nonbond list when any atom has move
! > 0.5 * (CUTNB-CTOFNB) since last update
! read initial coordinates again and compare
read coor card comp name @FILE_wmin.crd
coor orie rms sele segid bpti end
write coor card name @FILE_0.crd
* 4pti in @XSIZ x @YSIZ x @ZSIZ waterbox
* After initial 50ps CPT dynamics
* RMSD vs minimized: ?RMS
*

! comment out the stop if you want to continue immediately with the next piece
stop
set N 1
goto RUN_DYN
!!!!!!!!!!!!!!!! RESTART AND CONTINUE DYNAMICS
! have to get a few things in place before restarting dynamics
label GET_PSF
read psf card name @FILE-solv.psf
read coor card name @FILE-solv.crd

! setup PBC, this is not saved in the psf
! get variables XSIZ,YSIZ, and ZSIZ
stream @FILE-box-size.str
! at restart the actual boxsize is taken from the restart file
! we just need to setup the transformations with a reasonably similar box
crystal define orthorombic @XSIZ @YSIZ @ZSIZ 90.0 90.0 90.0
crystal build cutoff 14.0
image byresidue sele resn tip3 end
image bysegment sele segi bpti end

shake bonH param

label RUN_DYN

set M @N
decr M by 1
open unit 11 read form name @FILE-@M.res
open unit 12 write form name @FILE-@N.res
open unit 21 write unform name @FILE-@N.cor
open unit 22 write form name @FILE-@N.mon

dynamics restart timestep 0.002 nstep 50000 -
nprint 100 iprfrq 500 -
PCONst pmass 400.0 pgamma 20.0 tbath 300.0 PREFerence 1.0 -
ieqfrq 2000 ichecw 1 twindl -5.0 twindh +5.0 iasors 0 - !keep checking T, may not be needed
atom cdie fshift vshift cutnb 14.0 ctofnb 12.0 -
inbfrq -1 ihbfrq 0 imgfrq -1 - ! heuristic update
iunread 11 iunwrit 12 iuncrd 21 kunit 22 nsavc 100

! read initial coordinates again and compare
read coor card comp name @FILE_wmin.crd
coor orie rms sele segid bpti end
write coor card name @FILE_@N.crd
* 4pti in @XSIZ x @YSIZ x @ZSIZ waterbox
* After initial 50ps + @N * 100ps CPT dynamics
* RMSD vs minimized: ?RMS
*
Posted By: chuan Re: run-cpt-md.inp - 09/06/04 07:26 AM
"
! get aa definitions and parameters
read rtf card name $CHM_HOME/toppar/top_all22_prot.inp
read para card name $CHM_HOME/toppar/par_all22_prot.inp

! name of the system, used in filenames
set FILE 4pti
if @N .gt. 0 goto get_psf

"


In the first few command line above, N doesnt seem to be initiated. Should " N=0 " be added?
Posted By: lennart Re: run-cpt-md.inp - 09/06/04 09:42 AM
That would work. You may also want to read the instructions at the beginning of the script.
Posted By: tnc Re: run-cpt-md.inp - 11/09/04 04:47 PM
I have run a simulation using your script.
The only difference is the fact that I used PME.
The following is my input file.

-----------------------------------

crystal define orthorombic @XSIZ @YSIZ @ZSIZ 90.0 90.0 90.0
crystal build cutoff 14.0
image byresidue sele resn tip3 end
image bysegment sele segi prot end

shake bonH param

label RUN_DYN

set M @N
decr M by 1
open unit 11 read form name @FILE-@M.res
open unit 12 write form name @FILE-@N.res
open unit 21 write unform name @FILE-@N.dcd

dynamics restart timestep 0.002 nstep 250000 -
nprint 500 iprfrq 500 finalt 300.0 -
PCONst pmass 400.0 pgamma 20.0 tbath 300.0 PREFerence 1.0 -
ieqfrq 2000 ichecw 1 iasors 0 - !keep checking T, may not be needed
atom cdie fshift vshift cutnb 14.0 ctofnb 12.0 -
ewald pmewald kappa 0.36 spline order 4 -
inbfrq -1 ihbfrq 0 imgfrq -1 ntrfrq 500 - ! heuristic update
iunread 11 iunwrit 12 iuncrd 21 kunit -1 nsavc 500

-----------------------------------------------

I have got a weird result.
During the simulation, the water box was shrinking along Y axis and being expanded along X and Y axis ( to keep the volume ? ).

When I run a constant volume and PME simulation, this event wasn't observed.

Would you point out the problems of my script?
Posted By: lennart Re: run-cpt-md.inp - 11/09/04 05:03 PM
For a symmetric system you are better off using a symmetric unit cell (CUBIC, RHDO,OCTA), which don't suffer from this particular problem of fairly constant volume but changing shape. In practice you may want to consider the truncated octahedron or rhombic dodecahedron even for an elongated molecule to avoid problems that may arise if (when) the solute rotates. Since these unit cells are quite compact you may even gain in total number of particles.
Posted By: tnc Re: run-cpt-md.inp - 11/09/04 10:06 PM
I'm sorry. Still I can't understand you. Does your answer mean that if anyone use ORTH periodic boundary condition with CPT, the water box shape will always be changed? please let me know the script's problem and how I can accomplish ORTH PB condition simulation. I used your script, solvent-box.str, and I got the box

* PBC BOX DEFINITION FOR INTER SIMULATION
*
SET XSIZ 62.2091
SET YSIZ 56.0373
SET ZSIZ 54.5949
RETURN

and the following is part of the output.

DCNTRL> Constant pressure requested.
Reference pressure tensor (XX,YY,ZZ)= 1.0 1.0 1.0 atm.


I don't know why the water box keeps shrinking along Y axis.
Posted By: rmv Re: run-cpt-md.inp - 11/10/04 02:52 AM
One should not use the ORTHorhombic lattice for CPT dynamics; it can change shape too easily, in undesired ways. An isotropic shape should be used.
Posted By: tnc Re: run-cpt-md.inp - 11/10/04 07:38 PM
Does it depend on the shape of the system(protein)?
Would you explain more?
Thank you.
Posted By: tnc Re: run-cpt-md.inp - 11/18/04 06:23 PM
Nobody answers my question...
I can run CPT with octaheral box and cubic box, BUT..

Is there anybody who can run CPT MD with orthorhombic periodic boundary condition?
Is there anybody who is interested in solving this problem?
Posted By: rmv Re: run-cpt-md.inp - 11/18/04 07:36 PM
As noted above, while it is easily possible to run CPT MD with an ORTHorhobmic lattice, the unit cell doesn't maintain it's shape. Not for proteins, not for for small carbohydates, not for any system. It's a bad choice; one should use an isotropic shape (CUBIC, RHDO, OCTA) for CPT dynamics of a solute in water.

A secondary problem for simulations of a few ns or longer is that the protein can and will rotate during the course of simulation, such that it's long dimension could rotate into one of the shorter dimensions of the lattice, and come into contact with image copies of itself, which is usually highly undesirable.

A final note-- there are rumours (but nothing in the .doc files) of a new RECT lattice type in c31b1 which applies an additional constraint that the ratio of cell edges remains constant during pressure based expansion and contraction with CPT dynamics. This would prevent the collapse along one of the dimensions that commonly occurs with ORTH in CPT dynamics, but is still only suitable for short simulations because of the rotation problem cited above.

Moral: Making a box that "just fits" your protein is a bad idea for an MD simulation, esp. with CPT.
Posted By: tnc Re: run-cpt-md.inp - 11/18/04 08:24 PM
Thank you for all yours answers and patience.
Posted By: label Re: run-cpt-md.inp - 02/24/05 05:53 AM
Dear all ,
I used the script to caculate the density (water + protein) , and find the final density = 0.212099 .
It's not reasonable . The density (water + protein)
should be near 1 (g/cc) .
yeng tseng

http://165.112.184.13/ubbthreads/showflat.php?Cat=0&Number=4274&an=0&page=0
"
scalar mass stat
set m ?STOT
calc rho = 1.66 * @M / ?VOLU
"
Posted By: lennart Re: run-cpt-md.inp - 02/24/05 09:01 AM
Yes you seem to have too little water in you system. How big is your protein? How many water molecules do you have? What does your box look like (edgelengths)?
below given output results

calc XSIZ = INT( ?XMAX - ?XMIN + 16)
RDCMND substituted energy or value "?XMAX" to "20.1479"
RDCMND substituted energy or value "?XMIN" to "-18.9459"
Evaluating: INT(20.1479--18.9459+16)
Parameter: XSIZ <- "55"

CHARMM> calc YSIZ = INT( ?YMAX - ?YMIN + 16)
RDCMND substituted energy or value "?YMAX" to "13.4203"
RDCMND substituted energy or value "?YMIN" to "-12.5572"
Evaluating: INT(13.4203--12.5572+16)
Parameter: YSIZ <- "41"

CHARMM> calc ZSIZ = INT( ?ZMAX - ?ZMIN + 16)
RDCMND substituted energy or value "?ZMAX" to "11.4326"
RDCMND substituted energy or value "?ZMIN" to "-13.2873"
Evaluating: INT(11.4326--13.2873+16)
Parameter: ZSIZ <- "40"

orthogonal lattice of protein +water has a size of 55 41 40angs, how could i find total thickness of the water shell
i wish to extract the system protein +water at different thickness?(commands like show shth ? show shell?)

inorder to chnage the existing script as a cubic system , besides size of the protein i added different integers to find the uniform size of cubic crystal

calc XSIZ = INT( ?XMAX - ?XMIN + 2)
RDCMND substituted energy or value "?XMAX" to "20.1479"
RDCMND substituted energy or value "?XMIN" to "-18.9459"
Evaluating: INT(20.1479--18.9459+2)
Parameter: XSIZ <- "40"

CHARMM> calc YSIZ = INT( ?YMAX - ?YMIN + 15)
RDCMND substituted energy or value "?YMAX" to "13.4203"
RDCMND substituted energy or value "?YMIN" to "-12.5572"
Evaluating: INT(13.4203--12.5572+15)
Parameter: YSIZ <- "40"

CHARMM> calc ZSIZ = INT( ?ZMAX - ?ZMIN + 16)
RDCMND substituted energy or value "?ZMAX" to "11.4326"
RDCMND substituted energy or value "?ZMIN" to "-13.2873"
Evaluating: INT(11.4326--13.2873+16)
Parameter: ZSIZ <- "40"

!now i defined the uniform size
crystal define cubic @xsize @ysize@zsize 90 90 90
crystal facility cutoff 14.0
is this is the way of doing cubic steup?
"coor stat sele protein end" will give you the dimensions of your protein, and
"coor stat" will give the total dimensions. This may help you figure out the thickness of your "water shell" (at least along the x,y,z axes).

If you know that you want a 40x40x40 cubic box all the CALC commands are unnecesary. This suffices:
crystal define cubic 40.0 40.0 40.0 90.0 90.0 90.0
crystal build cutoff 20.0
If your cube is packed with atoms a shorter cutoff can be used; the important fact is that you want 26 images for a cubic system.

Use the POST function of the forums for new posts, not the REPLY.
Posted By: vbrahma Re: run-cpt-md.inp - 10/04/07 12:25 PM
I was using run-cpt-md.inp script, But since ORTH is not recommended I substituted it with another isotropic shape viz OCTA

but I encountered following error:-


CHARMM> ! setup Periodic Boundary Conditions
CHARMM> ! Here we have to use the crystal facility
CHARMM> ! Variables XSIZ,YSIZ, and ZSIZ have already been defined
CHARMM> crystal define OCTAhedral @XSIZ @YSIZ @ZSIZ 109.47 109.47 109.47
Parameter: XSIZ -> "65.0"
Parameter: YSIZ -> "65.0"
Parameter: ZSIZ -> "65.0"
XTLSYM> DIFF,RSMALL 0.3053E-02 0.1000E-09
XTLSYM> XTLTYP=OCTA
XTLABC: 62.546656 -12.508314 62.546656 -12.508314 -12.508314 62.546656

***** LEVEL -4 WARNING FROM <XTLSYM> *****
***** Shape matrix symmetry is not conserved
******************************************
BOMLEV ( -2) IS REACHED - TERMINATING. WRNLEV IS 5

Can someone suggest what went wrong? I am attaching my output script.

Attached File
15779-run-cpt-md_old.txt  (680 downloads)
Posted By: lennart Re: run-cpt-md.inp - 10/04/07 12:43 PM
For the TO you need all the digits in for the angle (see setup-truncatd-octahedron.str). I find the RHDO more robust to use in this respect.
Posted By: vbrahma Re: run-cpt-md.inp - 10/04/07 07:43 PM
Thanks.. I am using RHDO and i am not facing any problem now!
Posted By: joseph Re: run-cpt-md.inp - 05/22/08 03:27 PM
Dear Lennart,

I have a question regarding the temperature control in this script. I can't see neither BERE nor NOSE nor HOOVER keywords. Yet you use "ichecw 1" and "iasors 1".

Does it mean that the temperature is controlled by scaling the velocities when temperatures lies outside TWINDH\L?

Thank you very much.
Posted By: lennart Re: run-cpt-md.inp - 05/22/08 03:33 PM
Yes - from dynamc.doc:
ICHECW 1 The option for checking to see if the average temperature
of the system lies within the allotted temperature window
(between FINALT+TWINDH and FINALT+TWINDL) every
IEQFRQ steps.

For a well behaved system the temperature is usually quite stable and this simple scheme seldom interferes to rescale velocities (if at all); if it has t do many rescaling it works in a similar way as a Berendsen thermostat.
Posted By: rmv Re: run-cpt-md.inp - 05/22/08 05:00 PM
Alternatively, one can set ICHECW 0 IEQFRQ 0 and then add a line such as

HOOVER TMASS 2000. REFT 300. -

and get a true thermostat, coupled to a degree of freedom; with PGAMMA 0, this is rigorous NPT. With a non-zero PGAMMA, it is equivalent to a chain of thermostats.
© CHARMM forums