Topic Options
#394 - 10/23/03 06:49 AM run-pbc-md.inp
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4741
Loc: ~ 59N, 15E
*FILENAME: run-pbc-md.inp
*PURPOSE: setup and run Periodic Boundary Condition (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-pbc-md.inp > run-pbc-md.inp N=0 (or 1,2,....)
! Saves coordinate trajectory, restartfile and snapshot at end of traj. piece

! 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 use the image facility, but crystal works as well and
! for a simple minimum image setup there is also the pbound method.
! Variables XSIZ,YSIZ, and ZSIZ have already been defined
open unit 10 read form name $CHM_STREAM/ortho.img
read image card unit 10
! 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
image bysegment sele segi bpti end

! 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)
! even if you want to run NVE, it can be advantageous to run a few ps using CPT
! (see run-cpt-md.inp) to adjust the box volume
dynamics start time 0.002 nstep 50 - !25000 iseed 31415 - ! a new iseed gives a different trajectory
firstt 48.0 finalt 298.0 teminc 10.0 ihtfrq 500 -
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 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

! setup PBC, this is not saved in the psf
! get variables XSIZ,YSIZ, and ZSIZ
stream @FILE-box-size.str
open unit 10 read form name $CHM_STREAM/ortho.img
read image card unit 10
image byresidue sele resn tip3 end
image bysegment sele segi bpti end

shake bonH param

label RUN_DYN

! initial 50 ps heating MD run, no langevin dynamics with the heating
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 50 - !50000 -
nprint 100 iprfrq 500 -
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 dynamics + @N * 100ps
* RMSD vs minimized: ?RMS
*
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#30569 - 08/27/12 02:41 PM Re: run-pbc-md.inp [Re: lennart]
BeWier Offline
Forum Member

Registered: 01/10/11
Posts: 18
Dear Lennard,

Do you have a paper just using methods listed in run-pbc-md.inp? I believe this will help me understand PBC MD simulation much more easily and I also can cite your paper to support my method, because I did use the script (run-pbc-md.inp) to run my simulations.

Also, my enzyme has its highest catalytic activity at PH=3 condition, so the whole system is positive charged. What kind of counterions shall I add?

Thanks a lot,
Bewier

Top
#30576 - 08/29/12 07:08 AM Re: run-pbc-md.inp [Re: lennart]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4741
Loc: ~ 59N, 15E
No, I don't have such a paper.
Your question is quite general, and you need to read some general papers about MD simulations; the two CHARMM papers in J Comp Chem 1983 and 2009 are a good starting point. There is nothing special in the script you used. Welcome back with more specific questions if there is something that you don't quite understand. Note that it is important to have this understanding before you begin your simulations...
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top

Moderator:  chmgr, John Legato, petrella