CHARMm Principles



7       Setting Stochastic Boundaries

 Stochastic boundaries allow you to simulate infinite systems by using standard molecular dynamics techniques for a local region of the molecule while other parts of the structure are treated stochastically. For example, in studying the solvated active site of a protein, you can use molecular dynamics to study the active site and simulate solvent molecules by setting a stochastic boundary.

This chapter explains

This chapter explains the basic features of stochastic boundary molecular dynamics and the general procedure for running a calculation. A scripted example of how to set up a stochastic boundary is included at the end of the chapter.


Basic features

Stochastic boundary molecular dynamics uses elements of both Langevin dynamics and Newtonian dynamics. The goal of the method is to eliminate atoms distant from, for example, an active site, allowing detailed studies of a spatially localized portion of the reacting molecular system.

The basic features of the model are similar to Langevin dynamics. Langevin dynamics represents an approximation method for eliminating unimportant or uninteresting degrees of freedom from dynamics simulation. The effects of eliminated degrees of freedom on the system are simulated by boundary (mean) and stochastic forces. To use this model, you apply boundary forces and stochastic forces. The kinetic bath effects (dissipation and fluctuation) balance each other to maintain thermal equilibrium.

A solvent boundary force is computed from the analytic deformable boundary model. Structure boundary forces are derived from atomic mean-square displacements. In addition to these mean forces, stochastic forces exist that are analogous to those from Langevin dynamics.

This stochastic boundary method also requires you to partition the structure around a localized region. This is to formally isolate the region in which you are interested. This partitioning is schematically illustrated in the following figure:

Newtonian dynamics is applied for the reaction region and Langevin dynamics for the buffer region.


General procedure

The steps for running a stochastic boundary molecular dynamics calculation are as follows:

1.   Compute the required mean forces for input into the reduced particle equations of motion.

2.   Set up and partition the full system into a reaction region, a buffer region, and a reservoir region.

3.   Generate a stochastic boundary potential and read it into CHARMm.

4.   Set up the mapping CHARMm uses to connect the table entries with the boundary constrained atoms.

5.   Thermalize and equilibrate the reduced particle system.

6.   Carry out production dynamics simulation and analyze the trajectories for desired properties.


Example: Setting up a stochastic boundary

 This script is an example of how to set up a stochastic boundary. Because the script is too long to include in its entirety, critical fragments are presented and discussed. For clarity, these script fragments differ somewhat from those in actual scripts.

This example uses the following files:

SBMD1.STR

SBMD2.STR

SBMD3.STR

SBMD4.STR

SBMD5.STR

SBMD6.STR

SBMD7.STR

SBMD8.STR

SBMD9.STR

SBMD1.STR reads the initial coordinates and sets up the PSF. For convenience, the region of interest is translated to (0,0,0). A 15 Å sphere of water centered at (0,0,0) is appended and waters that overlap the protein are removed. The fragment of the script that codes for these activities follows:


open unit 1 form read name "$CHM_DATA/SPHERE15.CRD"
read sequ coor resi unit 1
rewind unit 1
generate SOLV nodi
read coor card append unit 1
close unit 1
! delete waters which overlap protein
dele atom sele .byres. (segid SOLV .and. type OH2 .and. -
  ((.not. segid SOLV .and. .not. (hydrogen .or. lone)) -
  .around. 2.8 )) end
! delete all solvent atoms more than 12A from the center
! of the reaction zone
dele atom sele .byres. (.not. (point 0.0 0.0 0.0 cut 12) -
  .and. (segid SOLV .and. type OH2)) end
SBMD2.STR performs a brief minimization to remove any bad contacts.

SBMD3.STR performs the following tasks:

Note that the dynamics statements include the keywords LANGEVIN and RBUF required for this process:


*****  SBMD3.STR
*  (A)  set-up stochastic boundary for bulk water
*  (B)  begin thermalization of the solvent within the reaction ! zone
*
! (A) set up solvent boundary - for bulk waters and crystal
! waters only within the 11 A reaction zone.
! The boundary potential file wat11.pot was generated
! using supplementary software (see sbmd.doc for more
! details).
open unit 1 form read name WAT11.POT
sbound read unit 1
close unit 1
sbound set xref 0. yref 0. zref 0. -
  assign 1 sele (( .byres. (point 0.0 0.0 0.0 cut 11)) -
  .and. resn TIP3 .and. type O*) end
! set up solvent friction
scalar FBETA set 62.0 sele resn TIP3 .and. type OH2 end
! fix all protein atoms and all crystal waters outside
! the reaction zone
cons fix sele segid LYSZ .or.(.byres. (segid SOLV .and. resn 
! TIP3 .and. .not. (point 0.0 0.0 0.0 cut @2))) end
! shake bonds to hydrogens.  TIP3 bonds are automatically
! shaken
shake bonh tolerance 1.0e-06 param
! AT LAST! Run dynamics on the reaction zone waters.
! First without temperature scaling to get a higher 
! temperature.
open unit 3 form write name "LYSO_3A.RST"
dynamics langevin nstep 100 timest 0.001 -
  ilbfrq 1 iseed 685963 iasvel 1 firstt 300 -
  tbath 300 rbuf 9.0 nprint 10 iprfrq 100 -
  ieqfrq 0 twindh 10.0 twindl -10.0 -
  iunwrite 3 nsavc 100 nsavv 0
open unit 3 form    read name "LYSO_3A.RST"
open unit 4 form   write name "LYSO_3B.RST"
open unit 1 unform write name "LYSO_3.DCD"
! bring waters to final temperature, 300K
dynamics langevin restart nstep 100 timest 0.001 - 
  ilbfrq 10 iseed 364569 iasvel 1 firstt 300 -
  tbath 300 rbuf 9 nprint 100 iprfrq 100 -
  ieqfrq 0 twindh 10.0 twindl -10.0 -
  isvfrq 100 iunwrite 4 iunread 3 nsavc 100 nsavv 0  iuncrd 1
return
SBMD4.STR does another overlay of waters onto the system to fill any holes that might have occurred as the first waters settle into the protein surface.

SBMD5.STR performs additional equilibration dynamics on the water, still keeping the protein fixed.

SBMD6.STR defines the reaction and reservoir zones of the system. 


! define reaction region
scalar xcomp set 1.0 -
  sele (.byres. (point 0.0 0.0 0.0 cut 9)) -
 .and. .not. ((type N .or. type CA .or. type C .or. type O) -
 .and. .not. point 0.0 0.0 0.0 cut 10) end
! save reaction region flags in storage array 1
scalar xcomp store 1 ! This stores the selection in array "1".
! define initial buffer region
! RECALL 1 refers to the reaction region stored in array "1".
scalar ycomp set 1.0 -
  sele (.byres. (point 0.0 0.0 0.0 cut 11)) -
 .and. .not. recall 1 end
! save buffer region flags in storage array 2
scalar ycomp store 2
! define Langevin atoms
! RECALL 2 in the selection selects the buffer region
! stored in storage array "2".
scalar zcomp set 1.0 -
  sele recall 2 .and. .not. (hydrogen .or. lone .or. resn TIP3) end
! save Langevin region flags in storage array 3
scalar zcomp store 3
! define reservoir region atoms if any
scalar wcomp set 1.0 -
  sele .not. (recall 1 .or. recall 2) end
! write boundary flags, file flags.pro
open unit 1 form write name "FLAGS.PRO"
write coor card comp unit 1
* HEW - TRP 62 SBMD
*  Reduced region partitioning for 11A region about Trp 62.
*  Column 1:  Reaction region atoms, 9A .byres. partitioning
*  Column 2:  Buffer region atoms, .byres. within 11A but not 
* in 9A plus all mainchain atoms within 10A
*  Column 3:  Protein Langevin atoms
*  Column 4:  Reservoir region atoms
*
close unit 1
A large part of this script involves partitioning the system into different parts and applying constraints. The remainder (not included here) involves setting up various constraints on the reaction and buffer regions.



Last updated September 18, 1998 at 04:03PM PDT.
Copyright © 1997, Molecular Simulations, Inc. All rights reserved.