Scripts:HS est.py

From UC_Chemistry

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