Scripts:HS est.py

From UC_Chemistry

Revision as of 16:11, 1 September 2009 David Rogers (Talk | contribs) (script->scripts) ← Previous diff		Current revision David Rogers (Talk | contribs) (ucgrad info. moved to QCT_Bayesian_Statistics)
Line 3:		Line 3:
	== Downloading ==		== Downloading ==
-	[[Media:HS_est.tgz | HS_est.tgz]]	+	Download from [[Media:HS_est.tgz | HS_est.tgz]].
-	In order to use this script, you must have python with numpy installed on your system and have access to our own utilities directory, ucgrad.	+	In order to use this script, you will have to have our helper library, ucgrad, working -- see [[QCT_Bayesian_Statistics]] for details.
-		+
-	You can get ucgrad from [http://forcesolve.sourceforge.net forcesolve] (forcesolve/cg_topol/ucgrad) by downloading or directly using svn:	+
-	<pre><nowiki>	+
-	svn export http://forcesolve.svn.sourceforge.net/svnroot/forcesolve/trunk/cg_topol/ucgrad ucgrad	+
-	</nowiki></pre>	+
-		+
-	To access it using ''from ucgrad import *'' within python, you have three options:	+
-	* Add it to your python site-packages directory (/usr/lib64/python2.4/site-packages/ on my system).	+
-	* Set the variable PYTHONPATH to wherever you keep ucgrad (e.g. $HOME/src/ucgrad).	+
-	* Make sure it is in the current directory.	+
-		+
-	I use the first and second ones, having added a symlink in site-packages.	+
	== Use ==		== Use ==

---

Current revision

The following is our reference implementation of the the hard sphere free energy calculation described in JCP129(13):134505.

Downloading

Download from HS_est.tgz.

In order to use this script, you will have to have our helper library, ucgrad, working -- see QCT_Bayesian_Statistics for details.

Use

This program works on "histogram" files containing a list of numbers -- representing the count in each bin. A single example case is included and can be run with:

./HS_est.py -p HS.par -o out.expot

Compare out.expot with the included HS.expot and spreadsheet.

Generation of Input Files

The most complicated part of the usage is generating the input histogram files from simulation data. My procedure is as follows:

1. Convert your trajectory to gromacs .trr file format (VMD can usually do this)

name=run
NF=1000

vmd -dispdev text -psf system.psf -dcd $name.dcd <<.
set sel [atomselect top "protein or name SOD or name CLA or name CLM"]
animate write pdb $name.pdb beg 0 end $(($NF-1)) sel \$sel
animate write trr $name.trr beg 0 end $(($NF-1))
.

Of course, substitute your file names and trajectory lengths. You may also have to split up the trajectory using catdcd before processing if VMD chokes because your DCD is larger than your available memory (a violation of the UNIX streams ideal). cat and trjcat -settime can put the resulting pdb-s and trr-s together, respectively.

2. Create a file listing 'rmin' for every frame.

• 1. Create a gromacs index file using make_ndx
  2. Process the trajectory and index file using g_mindist
  3. If rmins between multiple groups were done, write your own script to find the min(rmin) for each frame

3. Create a 'weight' file listing the weights to be used with each frame.

• 1. Decide on an appropriate 'lambda' to use with your simulation -- see Scripts:calc_du.py
  2. Flag all frames to be excluded from the averaging process (because rmin < lambda) -- see Scripts:rmin2ex.py
  3. Create the weight file by combining the ex.xvg and the model potential re-weighting energies -- see Scripts:mkwt.py

4. Create the Histogram see Scripts:mkndhist.py

mkndhist.py -n -1 -start $x0 -step $h -w 1 -wf $name.wt \
                $name.rmin.xvg $N $name.rmin.his