*** Umbrella sampling for the calculation of the potential of mean force
***
The umbrella program is using similar connectivity and
coordinate files (coordinates are in CHARMM format) to run normal Verlet MD.
Output files are standard output that contains info on the simulation
progress, velocity and coordinate files (CHARMM dynamics format).
This program is written to calculate potential of mean force
along reaction coordinate using the umbrella sampling method.
It employs only cartesian
coordinates which makes it convenient for a variety of systems.
This program is based on MD with linear constraints:
The presentation of reaction coordinate with Cartesian coordinates
requires that rigid body motions may effect the value of reaction
coordinate. To avoid this situation we added six linear
constraints on rigid body motions.
At each point along reaction coordinate we imposed umbrella potential
which allows us to run effectvely in vicinity of point of interest .
The application of the method requires considerable human intervention
for more details of this particular implementationof umbrella sampling method
see G. Verkhivker , Ron Elber, Q.H.Gibson
"Microscopic modeling of ligand diffusion through
protein leghemoglobin. Computer simulations and experiments".
J.Amer.Chem.Soc. 1992, in press.
To obtain
the overlapping distributions one needs to histogram the printted value of q.
The center of the distributiuon is extracted from the path structures and
is printted at the beginning of the output as qq=[f]
~ CHAIN TEST FOR CONFORMATIONAL TRANSITION IN VALINE DIPEPTIDE
~
~debug
file conn name=(VAL.WCON) unit=10 read
~ the line below means that coordinates will be read from file f.crd
file rcrd unit=5 read
file wcrd name=(VALUMB.DCD) bina unit=12 wovr
file wvel name=(VALUMB.VCD) bina unit=13 wovr
#ste=1000 #equ=100 #pri=10 #wcr=5 list=0
#sve=20 rand=-3451187 step=0.001 grid=1 #str=2 forc=5.d0
rvrs=-1 temp=300.
rmax=9. epsi=1. cdie v14f=8. e14f=2.
action
PARAMETERS FOR UMBRELLA SAMPLING
file con name=(VAL.WCON) -
connectivity file which contains all parameters
temp temperature
#ste number of integration steps
#pri print each #pri steps
#wcr coordinates are written at step interval #wcr
#eqv velocity scaling at step interval #eqv during equilibration period
#sve velocity scaling at step interval #sve after equilibration period
step step size
#tes period for checking constraints
#equ number of thermalization steps
newv select new velocities each newv steps
rand initial seed for random vel. selection
grid number of segments( usually #str-1)
#str number of structures in file f.crd
forc force constant for umbrella potential
File f.crd should be provided by user and contains the
names of structures along reaction coordinate (in CHARMM format).
In example below the flie f.crd is
1.CRD
2.CRD
So the number of segments is 1 and number of structures is 2.
rvrs=-1/1 the slope of the curvlinear coordinate describing
the reaction coordinate is estimated by finite difference
of cartesian coordinates 2.CRD and 1.CRD (-1 the direction is
from 1 to 2, relevant only if you have more than two structures))
The line with "action" is self explanatory
>From a single simulation a probability density P(q)in the neighborhood of qi
is extracted as a trajectory average of delta function delta(q-q(t))
The equilibrium position is printted at the beginning of the output as qq
umbre> ~ CHAIN TEST FOR CONFORMATIONAL TRANSITION IN VALINE DIPEPTIDE
umbre> ~
umbre> ~debug
umbre> file conn name=(VAL.WCON) unit=10 read
The following unit number was assigned 26
umbre> file rcrd unit=5 read
Input taken from standard file
umbre> file wcrd name=(VALUMB.DCD) bina unit=12 wovr
The following unit number was assigned 27
umbre> file wvel name=(VALUMB.VCD) bina unit=13 wovr
The following unit number was assigned 28
umbre> #ste=1000 #equ=100 #pri=10 #wcr=5 list=0
umbre> #sve=20 rand=-3451187 step=0.001 grid=1 #str=2 forc=5.d0
umbre> rvrs=-1 temp=300.
umbre> rmax=9. epsi=1. cdie v14f=8. e14f=2.
umbre> action
pick> 14 particles picked out of 14
PARAMETERS FOR UMBRELLA SIMULATION:
temperature: 300.0
number of integration steps: 1000
print each 10 steps
number of path segments : 1
data is on unit: 6 coordinates are on unit 27
initial coordinates are read from unit : 5
Coordinate are read in CHAR style
coordinates are written at step interval : 5
velocity scaling at step interval of: 20
force constant for umbrella potential: 5.0
step size: .100E-02
period for checking constraints 500
number of thermalization steps 100
select new velocities each 1000 steps
initial seed for random vel. selection -3451187
debug ? F
number of structures to be compressed: 2
Structure file list is in: f.crd
Path file written in: p.out
******************************
reading structure 1.CRD
number in the PATH file 1
******************************
the number of atoms is 14
******************************
reading structure 2.CRD
number in the PATH file 2
******************************
the number of atoms is 14
initial path coordinates are taken from: p.out
qq = 9.737190508560473E-02
Number of particle pairs 40
eu is .0
current temperature is 312.494
Number of particle pairs 40
q = .0
*** END OF TRAJECTORY INITIALIZATION ***
eu is 2.656308348425582E-03
q = -1.185255918085814E-02
eu is 8.230762257061088E-03
q = -1.555524443719468E-02
eu is 1.239804089336322E-02
q = -6.083927513086563E-03
eu is 1.375920972353991E-02
q = 1.607426627433693E-02
eu is 1.451734530410688E-02
q = 4.516594964775493E-02
eu is 1.694560255778653E-02
q = 7.277373871627585E-02
eu is 2.127083662996143E-02
q = 9.159125648351745E-02
eu is 2.753986128364574E-02
q = 9.862236987288132E-02
eu is 3.762772034718076E-02
q = 9.634040671321991E-02
eu is 5.366596540169265E-02
q = 9.134137206457091E-02
**** PRINTING INFO *****
Kinetic Energy .7441E+01
Potential Energy (1) -35.614
Total Energy -.2817E+02
current temperature is 312.494
eu is 5.365730569839
.
.
. SOME OUTPUT REMOVED
.
.
.
q = 1.15971076319143
eu is 5.38214627509544
q = 1.15336744375706
eu is 5.40213916678386
q = 1.13734259827526
eu is 5.42130024484857
q = 1.11268315961673
eu is 5.43445225711807
q = 1.0841136370278
eu is 5.437423903556
q = 1.05824151113938
eu is 5.42834294862091
q = 1.04069851073807
errors in coordinates constraints
-.61062E-15 -.22204E-14 -.66613E-15 .72164E-15
.45797E-15 -.16879E-14
Calculated constraints -6.106226635438361E-16 -2.220446049250312E-15
-6.661338147750939E-16 7.216449660063517E-16 4.579669976578771E-16
-1.687885942125433E-15
Exact values -.24553489812112 -1.5978205455026 -.826292916040224
-.207097984425299 .1004649990649933 -1.562147175259532E-02
eu is 5.40878337729065
q = 1.03344663162832
errors in coordinates constraints
.55511E-16 .13878E-15 -.13878E-16 .15266E-15
.83267E-16 .10582E-15
Calculated constraints 5.551115123125782E-17 1.387778780781445E-16
-1.387778780781445E-17 1.526556658859590E-16 8.326672684688673E-17
1.058181320345852E-16
Exact values .0 .0 .0 .0 .0 .0
*** checking energy conservation ***
at step 1000 total energy is -19.3474182697219
tot. pot. -31.9164911518526 total kine. 12.56907288213079
*** assigning new velocities at step 1000
irand = 1
errors in coordinates constraints
.16653E-15 -.34694E-17 -.55511E-16 .11102E-15
-.11102E-15 -.45103E-15
Calculated constraints 1.665334536937734E-16 -3.469446951953614E-18
-5.551115123125782E-17 1.110223024625156E-16 -1.110223024625156E-16
-4.510281037539698E-16
Exact values .0 .0 .0 .0 .0 .0
**** PRINTING INFO *****
Kinetic Energy .1014E+02
Potential Energy (1) -31.916
Total Energy -.2178E+02
current temperature is 301.119
**** PRINTING INFO *****