Scripts:HS est.py
From UC_Chemistry
| Revision as of 16:10, 1 September 2009 David Rogers (Talk | contribs) (Script:HS est.py moved to Scripts:HS est.py) ← Previous diff |
Revision as of 16:11, 1 September 2009 David Rogers (Talk | contribs) (script->scripts) Next diff → |
||
| Line 49: | Line 49: | ||
| '''3.''' Create a 'weight' file listing the weights to be used with each frame. | '''3.''' Create a 'weight' file listing the weights to be used with each frame. | ||
| - | *# Decide on an appropriate 'lambda' to use with your simulation -- see [[Script:calc_du.py]] | + | *# Decide on an appropriate 'lambda' to use with your simulation -- see [[Scripts:calc_du.py]] |
| - | *# Flag all frames to be excluded from the averaging process (because rmin < lambda) -- see [[Script:rmin2ex.py]] | + | *# Flag all frames to be excluded from the averaging process (because rmin < lambda) -- see [[Scripts:rmin2ex.py]] |
| - | *# Create the weight file by combining the ex.xvg and the model potential re-weighting energies -- see [[Script:mkwt.py]] | + | *# 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 [[Script:mkndhist.py]] | + | '''4.''' Create the Histogram see [[Scripts:mkndhist.py]] |
| <pre>mkndhist.py -n -1 -start $x0 -step $h -w 1 -wf $name.wt \ | <pre>mkndhist.py -n -1 -start $x0 -step $h -w 1 -wf $name.wt \ | ||
| $name.rmin.xvg $N $name.rmin.his</pre> | $name.rmin.xvg $N $name.rmin.his</pre> | ||
Revision as of 16:11, 1 September 2009
The following is our reference implementation of the the hard sphere free energy calculation described in JCP129(13):134505.
Downloading
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.
You can get ucgrad from forcesolve (forcesolve/cg_topol/ucgrad) by downloading or directly using svn:
svn export http://forcesolve.svn.sourceforge.net/svnroot/forcesolve/trunk/cg_topol/ucgrad ucgrad
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
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.
- Create a gromacs index file using make_ndx
- Process the trajectory and index file using g_mindist
- 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.
- Decide on an appropriate 'lambda' to use with your simulation -- see Scripts:calc_du.py
- Flag all frames to be excluded from the averaging process (because rmin < lambda) -- see Scripts:rmin2ex.py
- 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
