about calculation of PMF from SMD trajectories (with pasted smd_config_file)

From: Falgun Shah (fhshah_at_olemiss.edu)
Date: Wed Mar 04 2009 - 12:04:52 CST

Dear NAMD users

I need your help regarding calculation of PMF with NAMD. I have performed
SMD simulation using NAMD (http://www.ks.uiuc.edu/Training/Tutorials/)
little different way then shown in deca alanine stretching tutorial for
computing PMF from deca alanine smd trajectories. i have attached
configuration file for your reference which i used for the calculation of
SMD for protein of my interest. i have all the trajectories with me. In the
tutorial for computing PMF in NAMD for deca alanine, the da_SMD_tcl.out file
is generated which has the information about the time, distance and applied
force. In my case, i have log file which has Time, components of the applied
force as well as the coordinates of the SMD atoms. However, i dont know how
should i implement this information to calculate work using the information
given in the tutorial. Also, the langvein dynamics was set to off during SMD
calculation in my case (please see the config file). does it have a
significant impact on the smd simulation? I have used spring force constant
of 5 kcal/mol/A and 0.02A/ps of constant velocity as shown in config file
for pulling out my ligand. Is this velocity is sufficient enough for
computing PMF. With this information available, how should i go ahead for
calculation of work as well as further average them according to jarzynksi's
averaging scheme (exponential average and 2nd order cumulant). is there any
scripts available for the extracting required information or Do i need to
re-run the SMD simulations as shown in the deca-alanine tutorial. I want to
calculate the PMF using the scripts mentioned in the tutorial. Please let me
know if any important parameters, i am missing in the configuration file
which might be crucial for the calculation of PMF. I would really appreciate
your reply. Thanks.

#############################################################
## JOB DESCRIPTION ##
#############################################################

#constant velocity pulling

#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################

structure
/ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/protein.psf
coordinates
/ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/out_1500.pdb
set outputname protein1_out

set temperature 310

# Continuing a job from the restart files
if {1} {
set inputname
/ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/protein_restart
binCoordinates $inputname.restart.coor
binVelocities $inputname.restart.vel ;#remove temperature if you use
this
extendedsystem $inputname.restart.xsc
firsttimestep 0
}
#############################################################
## SIMULATION PARAMETERS ##
#############################################################

# Input
paraTypeCharmm on
parameters
/ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/par_all27_prot_lipid.inp
#temperature $temperature
#Note: do not use initial velocity temperature if you have specified
.restart.vel file

# Periodic Boundary Conditions
# Note: Do not set the periodic cell basis if you have also specified
# an .xsc restart file!

#if {1} {
#cellBasisVector1 95.5 0. 0.
#cellBasisVector2 0. 75.7 0.
#cellBasisVector3 0. 0. 65.2
#cellOrigin 10.2 1.2 1.3
#}
wrapWater on
wrapAll on

# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 13.
switching on
switchdist 11.
pairlistdist 14.5

# Integrator Parameters
timestep 2.0 ;# 2fs/step
rigidBonds all ;# needed for 2fs steps
nonbondedFreq 1
fullElectFrequency 2
stepspercycle 10

# PME (for full-system periodic electrostatics)
if {1} {
PME yes
PMEGridSizeX 96 ;# 2^5 * 3
PMEGridSizeY 81 ;# 3^3 * 3
PMEGridSizeZ 72 ;# 2^3 * 3^2
}

# Constant Temperature Control

langevin off ;# do langevin dynamics
#langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
#langevinTemp $temperature
#langevinHydrogen no ;# don't couple langevin bath to hydrogens

# ConstantPressure Control (variable volume)
if {1} {
useGroupPressure yes ;# needed for rigidBonds
useFlexibleCell no
useConstantArea no

langevinPiston on
langevinPistonTarget 1.01325 ;# in bar -> 1 atm
langevinPistonPeriod 100.
langevinPistonDecay 50.
langevinPistonTemp 310
}

# Output
outputName $outputname

restartfreq 1000 ;# 500steps = every 1ps
dcdfreq 1000
xstFreq 1000
outputEnergies 1000
outputPressure 1000
binaryoutput no ;# give me the pdb instead of the .coor
binaryrestart yes

#fixed atoms Constraints (set PDB beta-column to 1)
if {1} {
fixedAtoms on
fixedAtomsFile out3_185.ref
fixedAtomsCol B
}

#############################################################
## EXTRA PARAMETERS ##
#############################################################
SMD on
SMDFile /ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/out3_185.ref
SMDk 5
SMDVel 0.00004 # 0.00004/2fs that means 0.02A/ps and 20 A/1ns
SMDDir 0.452 -0.859 -0.240
SMDOutputFreq 10

#############################################################
## EXECUTION SCRIPT ##
#############################################################

# Minimization
#if {1} {
#minimize 100
#reinitvels $temperature
#}
run 1000000 ;# 2000ps 1ps = 5000 steps

-- 
Falgun H shah
PhD candidate (3rd year)
Department of Medicinal Chemistry
2028, Natural Product Center
University of Mississippi
Ph No: 6629151286(O)
          662 801 5667(M)
email: fhshah_at_olemiss.edu
-- 
Falgun H shah
PhD candidate (3rd year)
Department of Medicinal Chemistry
2028, Natural Product Center
University of Mississippi
Ph No: 6629151286(O)
          662 801 5667(M)
email: fhshah_at_olemiss.edu

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:52:26 CST