Previous Thread
Next Thread
Print Thread
Page 1 of 3 1 2 3
run-cpt-md.inp
#373 10/23/03 10:48 AM
Joined: Sep 2003
Posts: 4,829
Likes: 3
lennart Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 4,829
Likes: 3
*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
Re: run-cpt-md.inp
lennart #374 09/06/04 07:26 AM
Joined: May 2004
Posts: 222
C
Forum Member
Offline
Forum Member
C
Joined: May 2004
Posts: 222
"
! 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?

Re: run-cpt-md.inp
chuan #375 09/06/04 09:42 AM
Joined: Sep 2003
Posts: 4,829
Likes: 3
lennart Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 4,829
Likes: 3
That would work. You may also want to read the instructions at the beginning of the script.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: run-cpt-md.inp
lennart #376 11/09/04 04:47 PM
Joined: Sep 2004
Posts: 33
T
tnc Offline
Forum Member
Offline
Forum Member
T
Joined: Sep 2004
Posts: 33
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?

Re: run-cpt-md.inp
tnc #377 11/09/04 05:03 PM
Joined: Sep 2003
Posts: 4,829
Likes: 3
lennart Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 4,829
Likes: 3
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.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: run-cpt-md.inp
lennart #378 11/09/04 10:06 PM
Joined: Sep 2004
Posts: 33
T
tnc Offline
Forum Member
Offline
Forum Member
T
Joined: Sep 2004
Posts: 33
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.

Re: run-cpt-md.inp
tnc #379 11/10/04 02:52 AM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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.


Rick Venable
computational chemist

Re: run-cpt-md.inp
rmv #380 11/10/04 07:38 PM
Joined: Sep 2004
Posts: 33
T
tnc Offline
Forum Member
Offline
Forum Member
T
Joined: Sep 2004
Posts: 33
Does it depend on the shape of the system(protein)?
Would you explain more?
Thank you.

Re: run-cpt-md.inp
lennart #381 11/18/04 06:23 PM
Joined: Sep 2004
Posts: 33
T
tnc Offline
Forum Member
Offline
Forum Member
T
Joined: Sep 2004
Posts: 33
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?

Re: run-cpt-md.inp
tnc #382 11/18/04 07:36 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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.


Rick Venable
computational chemist

Page 1 of 3 1 2 3

Moderated by  lennart, rmv 

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

PHP: 7.3.31-1~deb10u1 Page Time: 0.007s Queries: 35 (0.005s) Memory: 0.7788 MB (Peak: 0.8753 MB) Data Comp: Off Server Time: 2021-11-30 23:44:54 UTC
Valid HTML 5 and Valid CSS