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

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Wed Mar 04 2009 - 12:52:47 CST

On Wed, 4 Mar 2009, Falgun Shah wrote:

FS> Dear NAMD users
FS>
FS> I need your help regarding calculation of PMF with NAMD. I have performed
FS> SMD simulation using NAMD (http://www.ks.uiuc.edu/Training/Tutorials/)
FS> little different way then shown in deca alanine stretching tutorial for
FS> computing PMF from deca alanine smd trajectories. i have attached

FS> configuration file for your reference which i used for the calculation of
FS> SMD for protein of my interest. i have all the trajectories with me. In the
FS> tutorial for computing PMF in NAMD for deca alanine, the da_SMD_tcl.out file
FS> is generated which has the information about the time, distance and applied
FS> force. In my case, i have log file which has Time, components of the applied
FS> force as well as the coordinates of the SMD atoms. However, i dont know how
FS> should i implement this information to calculate work using the information
FS> given in the tutorial. Also, the langvein dynamics was set to off during SMD

how about checking out the literature that describe the
method in detail?

- Izrailev, Stepaniants, Isralewitz, Kosztin, Lu, Molnar,
  Wriggers, Schulten. Computational Molecular Dynamics: Challenges,
  Methods, Ideas, volume 4 of Lecture Notes in Computational Science and
  Engineering, pp. 39-65. Springer-Verlag, Berlin, 1998.

- Park, Schulten, J. Chem. Phys. 120 (13), 5946 (2004)

- Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)

axel.

FS> calculation in my case (please see the config file). does it have a
FS> significant impact on the smd simulation? I have used spring force constant
FS> of 5 kcal/mol/A and 0.02A/ps of constant velocity as shown in config file
FS> for pulling out my ligand. Is this velocity is sufficient enough for
FS> computing PMF. With this information available, how should i go ahead for
FS> calculation of work as well as further average them according to jarzynksi's
FS> averaging scheme (exponential average and 2nd order cumulant). is there any
FS> scripts available for the extracting required information or Do i need to
FS> re-run the SMD simulations as shown in the deca-alanine tutorial. I want to
FS> calculate the PMF using the scripts mentioned in the tutorial. Please let me
FS> know if any important parameters, i am missing in the configuration file
FS> which might be crucial for the calculation of PMF. I would really appreciate
FS> your reply. Thanks.
FS>
FS> #############################################################
FS> ## JOB DESCRIPTION ##
FS> #############################################################
FS>
FS> #constant velocity pulling
FS>
FS> #############################################################
FS> ## ADJUSTABLE PARAMETERS ##
FS> #############################################################
FS>
FS> structure
FS> /ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/protein.psf
FS> coordinates
FS> /ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/out_1500.pdb
FS> set outputname protein1_out
FS>
FS> set temperature 310
FS>
FS> # Continuing a job from the restart files
FS> if {1} {
FS> set inputname
FS> /ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/protein_restart
FS> binCoordinates $inputname.restart.coor
FS> binVelocities $inputname.restart.vel ;#remove temperature if you use
FS> this
FS> extendedsystem $inputname.restart.xsc
FS> firsttimestep 0
FS> }
FS> #############################################################
FS> ## SIMULATION PARAMETERS ##
FS> #############################################################
FS>
FS> # Input
FS> paraTypeCharmm on
FS> parameters
FS> /ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/par_all27_prot_lipid.inp
FS> #temperature $temperature
FS> #Note: do not use initial velocity temperature if you have specified
FS> .restart.vel file
FS>
FS> # Periodic Boundary Conditions
FS> # Note: Do not set the periodic cell basis if you have also specified
FS> # an .xsc restart file!
FS>
FS> #if {1} {
FS> #cellBasisVector1 95.5 0. 0.
FS> #cellBasisVector2 0. 75.7 0.
FS> #cellBasisVector3 0. 0. 65.2
FS> #cellOrigin 10.2 1.2 1.3
FS> #}
FS> wrapWater on
FS> wrapAll on
FS>
FS>
FS> # Force-Field Parameters
FS> exclude scaled1-4
FS> 1-4scaling 1.0
FS> cutoff 13.
FS> switching on
FS> switchdist 11.
FS> pairlistdist 14.5
FS>
FS>
FS> # Integrator Parameters
FS> timestep 2.0 ;# 2fs/step
FS> rigidBonds all ;# needed for 2fs steps
FS> nonbondedFreq 1
FS> fullElectFrequency 2
FS> stepspercycle 10
FS>
FS>
FS> # PME (for full-system periodic electrostatics)
FS> if {1} {
FS> PME yes
FS> PMEGridSizeX 96 ;# 2^5 * 3
FS> PMEGridSizeY 81 ;# 3^3 * 3
FS> PMEGridSizeZ 72 ;# 2^3 * 3^2
FS> }
FS>
FS>
FS> # Constant Temperature Control
FS>
FS> langevin off ;# do langevin dynamics
FS> #langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
FS> #langevinTemp $temperature
FS> #langevinHydrogen no ;# don't couple langevin bath to hydrogens
FS>
FS>
FS> # ConstantPressure Control (variable volume)
FS> if {1} {
FS> useGroupPressure yes ;# needed for rigidBonds
FS> useFlexibleCell no
FS> useConstantArea no
FS>
FS> langevinPiston on
FS> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
FS> langevinPistonPeriod 100.
FS> langevinPistonDecay 50.
FS> langevinPistonTemp 310
FS> }
FS>
FS> # Output
FS> outputName $outputname
FS>
FS> restartfreq 1000 ;# 500steps = every 1ps
FS> dcdfreq 1000
FS> xstFreq 1000
FS> outputEnergies 1000
FS> outputPressure 1000
FS> binaryoutput no ;# give me the pdb instead of the .coor
FS> binaryrestart yes
FS>
FS> #fixed atoms Constraints (set PDB beta-column to 1)
FS> if {1} {
FS> fixedAtoms on
FS> fixedAtomsFile out3_185.ref
FS> fixedAtomsCol B
FS> }
FS>
FS>
FS>
FS> #############################################################
FS> ## EXTRA PARAMETERS ##
FS> #############################################################
FS> SMD on
FS> SMDFile /ptmp/r1146/smd/smd8/3216.sequoia/3225.sequoia/out3_185.ref
FS> SMDk 5
FS> SMDVel 0.00004 # 0.00004/2fs that means 0.02A/ps and 20 A/1ns
FS> SMDDir 0.452 -0.859 -0.240
FS> SMDOutputFreq 10
FS>
FS> #############################################################
FS> ## EXECUTION SCRIPT ##
FS> #############################################################
FS>
FS> # Minimization
FS> #if {1} {
FS> #minimize 100
FS> #reinitvels $temperature
FS> #}
FS> run 1000000 ;# 2000ps 1ps = 5000 steps
FS>
FS>
FS>

-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.

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