*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
*


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden