Previous Thread
Next Thread
Print Thread
Page 1 of 3 1 2 3
#373 10/23/03 10:48 AM
Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
*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
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?

chuan #375 09/06/04 09:42 AM
Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
That would work. You may also want to read the instructions at the beginning of the script.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
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?

tnc #377 11/09/04 05:03 PM
Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
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
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.

tnc #379 11/10/04 02:52 AM
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
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

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.

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?

tnc #382 11/18/04 07:36 PM
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
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

rmv #383 11/18/04 08:24 PM
Joined: Sep 2004
Posts: 33
T
tnc Offline
Forum Member
Offline
Forum Member
T
Joined: Sep 2004
Posts: 33
Thank you for all yours answers and patience.

lennart #384 02/24/05 05:53 AM
Joined: Dec 2003
Posts: 80
L
Forum Member
Offline
Forum Member
L
Joined: Dec 2003
Posts: 80
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
"

label #385 02/24/05 09:01 AM
Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
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)?


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Joined: Oct 2005
Posts: 206
P
Forum Member
Offline
Forum Member
P
Joined: Oct 2005
Posts: 206
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?

Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
"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.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
tnc #388 10/04/07 12:25 PM
Joined: Nov 2006
Posts: 11
Forum Member
Offline
Forum Member
Joined: Nov 2006
Posts: 11
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 Images
15779-run-cpt-md_old.txt (0 Bytes, 744 downloads)
Last edited by vbrahma; 10/04/07 12:33 PM.

Vijaya Brahma Research Scholar Institute of MIcrobial Technology India
vbrahma #389 10/04/07 12:43 PM
Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
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.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
lennart #390 10/04/07 07:43 PM
Joined: Nov 2006
Posts: 11
Forum Member
Offline
Forum Member
Joined: Nov 2006
Posts: 11
Thanks.. I am using RHDO and i am not facing any problem now!


Vijaya Brahma Research Scholar Institute of MIcrobial Technology India
lennart #391 05/22/08 03:27 PM
Joined: Feb 2008
Posts: 3
J
Forum Member
Offline
Forum Member
J
Joined: Feb 2008
Posts: 3
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.

joseph #392 05/22/08 03:33 PM
Joined: Sep 2003
Posts: 4,849
Likes: 6
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,849
Likes: 6
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.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
lennart #393 05/22/08 05:00 PM
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
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.

Page 1 of 3 1 2 3

Moderated by  lennart, rmv 

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

PHP: 7.3.31-1~deb10u1 Page Time: 0.012s Queries: 57 (0.006s) Memory: 0.8546 MB (Peak: 1.0098 MB) Data Comp: Off Server Time: 2022-09-26 12:02:09 UTC
Valid HTML 5 and Valid CSS