From tinkergpu
Jump to: navigation, search

Modeling books in dropbox


Running Gaussian QM job in tacc

Stampede each node in the Normal queue has 16 core and 32GB memory. Each job can run for 2 days.

There are additional 16 nodes in the “largemem” queue, each with 1TB memory and 32 core.

Create a submission script for each job. Only the job name, name of standard out, number of cores and the last line need to be changed for each job:

For example: qm.job:



# Example SLURM job script to run MPI applications on 

# TACC's Stampede system.


# $Id: job.mpi 1580 2013-01-08 04:10:50Z karl $



#SBATCH -J dmpwat           # Job name

#SBATCH -o dmpwat.%j.out   # Name of stdout output file (%j expands to jobId)

#SBATCH -p largemem             # Queue name

#SBATCH -N 1                  # Total number of nodes requested (16 cores/node)

#SBATCH -n 16                 # Total number of mpi tasks/cores requested. If you want to run 2 QM jobs on one node, this needs to be 8.

#SBATCH -t 48:00:00           # Run time (hh:mm:ss) - 48 hours

#SBATCH -A TG-MCB100057       # <-- Allocation name to charge job against

# Launch the MPI executable named "a.out"

module load gaussian

module list

g09 < DMPwat-mp2.com > DMPwat-mp2.out


To submit the job on stampede, log into stampede

ssh myusername@stampede.tacc.utexas.edu

It is best to run the jobs in work directory:



mkdir DMP

cd DMP

vi qm.job #edit the content as above example


in the *.com file, specify the Nproc and mem accordingly. For example

Nproc=8  and mem=16GB will allow you to run two QM jobs on each node in the "normal' queue

Nproc=16 and mem=32GB will work for one QM job on each node in the "normal' queue

Nproc=16  and mem=300GB will work for two QM jobs on each node in the "largemem' queue


To submit the job

sbatch qm.job

To check the status

squeue –u myusername

To cancel the job

scancel jobid_from_squeue

You can write a script to submit several jobs via a loop as above.


Download from


Installed on bme-pluto (python 2.6)and can run on node100 etc.

python setup.py install --user --install-scripts ~pren/bin

To run:

~pren/bin/propka31 1hpx.pdb 
propka31 -help

Output has pka for each residue:

Group      pKa  model-pKa   ligand atom-type
ASP  15 X     2.90       3.80                      
ASP  30 X     4.37       3.80                      
ASP  31 X     4.02       3.80       

From the above pka, you can estimate that at pH=7 all ASP is negatively charged. At pH=4, the bottom two may be neutral.

Brownian Dynamics and nanoparticle

Current working directory for this project is: /users/gchattree/Research/brownian/
Email at gaurav.chattree@utexas.edu for questions.
Brownian Dynamics Tutorial Using Browndye:

There are two version of Browndye installed.
Strangely you will need them both.

/home/gchattree/Research/browndye/bin/ and /home/muxiao/software/browndye/bin/
I use the * character extensively to represent parts of the file name differs for each simulation. In general * = molecule name. But it varies depdending on context.
Browndye manual is found: http://browndye.ucsd.edu/browndye/doc/manual.html
Browndye tutorial is found: http://browndye.ucsd.edu/browndye/thrombin-example/tutorial12.html
A sample directory can be found: /users/gchattree/Research/brownian/spring-2013/trastuzumab/gold-ninth-A1-1/


1. Download your two pdb's. Give them both succinct names.
2. Go to: http://kryptonite.nbcr.net/pdb2pqr_1.8/ to convert these pdb files to pqr.
-Note, many issues can occur here. Any residue label that has a letter in it, such as 99B, will be ignored. Change these residue labels to a number. 99B -> 999.
-HETATMs are also not usually recognized and have to be added to the pqr manually.
-Be sure to generate the APBS input file as well.
-Download both the *.in file and the *.pqr file
-Overlook the *.pqr file to make sure everything looks okay.

-Delete all the lines containing just TER from the file
3. Now enter the *.pqr file in vim. Run the following command: %s/\([0-9]\)-\([0-9]\)/\1 -\2/
-This adds spaces between numbers in case they got conjugated during pdb2pqr
4. Certain nodes have APBS installed. One of these is bme-earth. One one of these nodes run: apbs *.in
-If this takes much too long or there is a memory error, edit the the dime line to something smaller, such as: dime 61 61 61.
-Check the apbs output in pymol using the APBS visualizer tool. There is a chance that some of the potential map got cut off. In that case increase the size of the cglen line in the *.in file and re run apbs.
-If you want more info about the keywords in the APBS input file, the APBS website has a good faq.
5. Now rename the APBS output file. I usually do something like: mv *.dx *-PE0.dx
6. Now run pqr2xml to create an atoms xml file as follows: /home/gchattree/Research/browndye/bin/pqr2xml < *.pqr > *-atoms.xml
7. Now the annoying part. You must create a pairs.xml file. You must decide which pairs will count towards the reaction occuring. I use pymol's view hydrogen bonding tool to decide on these atoms.
-Example pair file can be found: /users/gchattree/Research/brownian/spring-2013/trastuzumab/gold-ninth-A1-1/pairs.xml
-You can also use the default command to create the file: /home/gchattree/Research/browndye/bin/make_rxn_pairs -mol1 molecule1.pqrxml -mol2 molecule2.pqrxml -ctypes contacts.xml -dist 5.0 > pairs.xml
8. Now you make the reactions file (*-rxns.xml). Basically here you just define what distance between each pair will constitute as binding, and how many bound pairs are needed to count as a reaction (keywords 'distance' and 'nneeded' respectively). Example command: /home/gchattree/Research/browndye/bin/make_rxn_file -pairs pairs.xml -distance 4.0 -nneeded 2 > mol0-mol1-rxns.xml
9. Now copy over the following input file and change the file names as necessary and make other edits as you see fit (such as changing 'n-threads'): /users/gchattree/Research/brownian/spring-2013/trastuzumab/gold-ninth-A1-1/input.xml
10. Now run /home/muxiao/software/browndye/bin/bd_top input.xml to create all the necessary input files.
-must use the latest bd_top binary which is found in muxiao's directory
-If this step takes too long you may have to kill it and alter the size of your system.
-IMPORTANT: You may want to check the *-surface.xml file after this runs. If there seems to be way too few atoms on the surface, surface_spheres may have to be run manually with a different probe size, and then bd_top should be rerun.

--visualize the surface as follows:

a. cat gold_A1-atoms.xml | surface_spheres -probe_radius 4.0 > gold_A1-surface.xml (in this example 4.0 is set as probe size. default is 1.6)

b. make_surface_sphere_list -surface gold_A1-surface.xml -rxn1 R-gold_A1-rxns.xml -spheres gold_A1-atoms.xml > gold_A1-surface-atoms-hard.xml

c. cat gold_A1-surface-atoms-hard.xml | make_soft_spheres > gold_A1-surface-atoms.xml

d. vim gold_A1-surface-atoms.xml

--in vim, run: :%s/<\// <\// and then :%s/>/> /

e. now run: /home/gchattree/Research/brownian/scripts/find_surface.py

-- open find_surface.py in vim first and edit the file names as needed
11. Now run the simulation: /home/gchattree/Research/browndye/bin/nam_simulation *-simulation.xml &
-Wait a few days.

Progress can be measured two ways:
A. Vim the file 'results.xml' file
B.cat results.xml | /home/gchattree/Research/browndye/bin/compute_rate_constant
Option B provides information on the current state of the rate constant that is being calculated.

Output can be visualized as follows:
1. for i in {1..100}; do process_trajectories -traj traj.xml0.xml -index traj.xml0.index.xml -n $i | cat > ptraj.xml; xyz_trajectory -mol0 mol0-atoms.xml -mol1 mol1-atoms.xml -trajf ptraj.xml >> traj.xyz; done
2. vmd traj.xyz
Adjust the above command as needed. As written it adds first 100 end states calculated by the first thread to traj.xyz. All of these states can be viewed as a movie in vmd.


Creating the gold nanoparticle (AuNP) and conjugating molecules to the surface:

1. Go to /home/gchattree/Research/brownian/scripts/
2. Run creategold.py with the arguments diameter size and output file name

-example: python creategold.py 40 gold_out.xyz

-gold_out.xyz should be an AuNP with radius of about 40 angstrom

-This script can take quite a while to run depending on the size of the AuNP

-optional, run either of the creategoldshell files to turn the AuNP into just a shell.

-open the files and adjust the file names as needed and then run it

-the smooth creategoldshell file can take quite a while to run

3. Use openbabel to convert from xyz to pdb

-example: babel -ixyz gold_out.xyz -opdb gold_out.pdb

4. Trim the molecule if it is too large using the trimgold.py script

-Open the file to adjust file names as necessary

-the trim_general.py file works for trimming not just gold molecules. It needs the *-pairs.xml file to make sure it does not trim an atom in a pair.

5. Now run attach.py to attach molecules to the gold particles surface

-Open the file and set all the parameters as necessary. There are many parameters to be set.

-It is currently set up for attaching 9 Trastuzumabs (A1) to a gold nanoparticle of diameter ~18nm

-A2 = Pertuzumab

-Open the output in pymol to make sure everything is okay.

6. Go to: <a href="http://kryptonite.nbcr.net/pdb2pqr_1.8/">http://kryptonite.nbcr.net/pdb2pqr_1.8/</a> to convert the pdb files to pqr.

-All the gold atoms will be stripped during this step

-must be added manually

--How to add manually:

--a. Open the newly created AuNP *.pqr file and the original AuNP *.pdb file using vim's vsp command.

--b. Copy over all of the gold atom lines from the *.pdb file to the *.pqr file

--c. Select all of the gold atom lines in the *.pqr file using Shift+v and then run the following command: s/ 1.00 0.00/ 0.00 1.35/

---The gold atom lines should look like this now:

---HETATM 1 AU LIG 1 11.803 -22.201 31.558 0.00 1.35 Au

--d. The id's of the atoms are all wrong now however. To find the new id's open the script find_newIds.py, edit the parts marked by #set, such as file names and starting numbers, and then run it.

--e. Now do a copy and replace of these new ids over the old ones


Current Status:

The various runs can be found in /brownian/spring-2013/trastuzumab

As you can see the current results suggest no signifant change in association constant using this AuNP conjugation

Problems here involve how the surface is found. The surface probe size could be too big.

The results also don't make intuitive sense. We would expect greater association constants with so many antibodies in close proximity



Where to go now:
A smaller gold shell (smaller diameter) so there is less surface area.
More trim (The more you trim the larger you have to make the probe size for the surface_spheres command)
Faster computer
Manually adjusting inside_points so that it can use multiple cores

Having many receptors attached to like a cell membrane




pymol tips

to show the virus capsid, the whole biological assembly as oppose to a single asymmetric unit:

set all states,1


qchem 5

qchem 5 can be run on node40, 41, 42, 142 and g2-node38 (or node38).

Needed to install openmpi related rpm and the LD path (see below): yum install -y *openmpi*

You can run multiple jobs on each node using 4-8 cores.

setenv QCSCRATCH /opt/scratch/qchemscratch
setenv QCLOCALSCR /scratch/qchemscratch
setenv QCMPI mpich
source /opt/qchem/5/qcenv.csh
setenv LD_LIBRARY_PATH /usr/lib64/openmpi-1.10/lib/:$LD_LIBRARY_PATH

qchem 4.2

Qchem (4.2) can be run on the following nodes (all cores):

node50-57, 8 cores each, 7GB mem. node53, 55 dead now (9/2016).

Also work on these nodes with large memory:


The max memory can be used on the nodes is set to 50GB for now. If need more, please let me know.

Never set the memory in the qchem input more than the actual mem (free -g) on each node.

(mem & disk in MB):

mem_static  10000
mem_total 50000
ao2mo_disk   130000

/opt/qchem/4.2/config/preference MEM_TOTAl is set to 50GB

v 4.2

July 2015 (last version of free upgrade) change 4.1 below to run multithread version. MPICH version is in /opt/qchem/4.2mpi

setenv QCSCRATCH /opt/scratch/qchemscratch
setenv QCLOCALSCR /scratch/qchemscratch
setenv QCMPI mpich
source /opt/qchem/4.2/qcenv.csh

Then (see /home/pren/QuantumMechanics/qchem/ for example inputs)

qchem -nt 6 n60.in n60.out &

v 4.1

User manual: http://www.q-chem.com/qchem-website/manual_4-1.html

setenv QCSCRATCH /opt/scratch/qchemscratch;
setenv QCLOCALSCR /scratch
#setenv PBS_NODEFILE /home/pren/qchemnode
setenv QCMPI mpich
source /opt/qchem/4.1/qcenv.csh

to run multi-thread:

qchem -nt 6 n60.in n60.out &

before 4.1

source the following if you want to use version 4.0 (bash shell):

  1. QCHEM #qchem
    export QC
    export OMP_NUM_THREADS=1
    export QCAUX
    export QCPLATFORM
    export QCSCRATCH
    export QCLOCALSCR=/scratch
    export QCMPI=mpich
    export QCRSH=/usr/bin/ssh
    export PBS_NODEFILE=/home/pren/qchemnode
    export PATH=/work/opt/qchem4.0/bin/:$PATH
    if [ -e $QC/bin/qchem.setup.sh ] ; then


source the following if you want to use verion 3.2:

    setenv QC /opt/qchem-para if (-e "$QC"/bin/qchem.setup) source $QC/bin/qchem.setup;

    setenv QCAUX $QC/aux; setenv QCPLATFORM LINUX_Ix86_64;

    setenv QCSCRATCH /opt/scratch/qchemscratch;

    setenv QCLOCALSCR /scratch

    setenv QCMPI mpich

    setenv QCRSH /usr/bin/ssh setenv noclobber ""

    setenv PBS_NODEFILE /home/pren/qchemnode


To run (example: /home/pren/QuantumMechanics/qchem) qchem n60.in n60.out & qchem -pbs -np 7 2huw-sp.in 2huw-sp.out &


mpdboot debug

Cleanup first:

/usr/local/sbin/rcom "/opt/mpich2-fc13-p26/bin/mpdcleanup"

Do this form nova and This will cleanup your mpd on each node.
ssh nodexxx "/opt/mpich2-fc13-p26/bin/mpdcleanup"

~/.mpd.conf has the user's secretword, chmod 600 ~~/.mpd.conf

df -h make sure disk is not full

iptables stop

hosts file has all the nodes

<span style="font-size: larger;" />make sure you can ssh into the nodes


On the node you have trouble, you can check with "mpdtrace -l" if it is clean (mpdtrace: cannot connect to local....)


Make sure you can ssh Hosts file has all the nodes. The top line should be " localhost.bme.utexas.edu localhost"

Diable firewall:

/sbin/chkconfig --list iptables

Disable NetworkManager !

/sbin/chkconfig --list NetworkManger
/sbin/chkconfig --level 345 NetworkManager off
/etc/init.d/NetworkManager stop
/etc/init.d/network restart

You may need to add to resolv.conf


Make sure the hostname is correct

 less /etc/sysconfig/network
mpd&  ###here you may see error message about hostname ip etc.

You may have to reboot.




mpdboot [ -n <#nodes> ] [ -f <hostsfile> ] [ -h ] [ -r <rshcmd> ] \ [ -u <user> ] [ -m <mpdcmd> ] [ --loccons ] [ --remcons ] \ [ -s ] [ -d ] [ -v ] [ -1 ] [ --ncpus=<ncpus> ]


mpdboot [ --totalnum=<#nodes> ] [ --file=<hostsfile> ] [ --help ] \ [ --rsh=<rshcmd> ] [ --user=<user> ] [ --mpd=<mpdcmd> ] \ [ --loccons ] [ --remcons ] [ --shell ] [ --debug ] \ [ --verbose ] [ -1 ] [ --ncpus=<ncpus> ]


-h, --help
Display help message
-d, --debug
Print debug information
–v, --verbose
Print extra verbose information. Show the rshcmd attempts
-n <#nodes>
Number of nodes in mpd.hosts on which daemons start
-r <rshcmd>
Specify remote shell to start daemons and jobs
-f <hostsfile>
Path/name of file that contains the list of machine names on which daemons start.
Remove a restriction of starting only one mpd per machine
-m <mpdcmd>
Specify the full path name of mpd on the remote hosts
-s, --shell
Specify shell
-u <user>
Specify user
Do not create local MPD consoles
Do not create remote MPD consoles
Indicate how many processors to use on the local machine (the rest runs on other nodes)

/opt/mpich2 /opt/mpich2-fc7 for parallel mpi jobs.

If you have trouble with mpdboot, it is typically becasue disk is full, which is cuased by egregious amount of error message of previous mpi jobs. Oe see the debug lists under MPICH2 above.

Please always do this when you run paallel job using mpi to get rid of the mpi error messages.

"mpirun ....your prgram... 2> /dev/null"

2 means standard error, which is redirected to /dev/null (black hole) instead of system log.

If you have trouble running mpdboot (which is a script starting mpd), you can mannaully start mpd, based on mpich2 installation guide:

 mpd &                 # starts the local daemon  (i.e. node18)
mpdtrace -l         # makes the local daemon print its host and port in the form <host>_<port> 

Then log into second machines,  and do: 
mpd -h <hostname> -p <port> &
Note you can add --ncpus=4 to mpd on both nodes as well.

Here is an example (put /opt/mpich2-fc7/bin in PATH): log into node18

  mpd --ncpus=4 --listenport=55551 &

log into node19

 mpd -h node18 -p 55551 --ncpus=4 & 

log into node21

 mpd -h node18 -p 55551 --ncpus=4 &

Try mpdtrace to see which nodes are in the mpd ring. You should see node18, node19 and node21. Note node 18 is the head node here, and every other node will connect to it. The 55551 is a number I specified but you can change. cpus is also a variables. mpd -h gives mode info.







compute cd spectrum from QM

Orca offer CIS and RDFT but the spectrum part is broken (author notified. will fix in next release) G03 offer TDFT: #B3LYP/6-31+G* TD(NSTATE=60)

The NSTATE can be determined by a quick test run. The initial state can also be changed (see g03 manual).

The output (log file) will list different states (eV/nm) and rotation strength R. /opt/software/CHEM/cdspectrum/ has the software fit the data into a spectrum.

   vi x (get the length R from g03 log file, copy log to x and delete everything but velocity R)
grep Singlet-A td.log | cut -c35-44 > x2 (get the eV)
paste x2 x > td.rot
vi td.rot (add header according to cdpsectrum, see righttd.rot)
/opt/software/CHEM/cdspectrum/plotspecTM < td.rot

see example in:

/home/pren/cd-qm/dpd (gpt file shows conversion frm ev to nm)

GABEDIT can be used to extract and vis the CD results.



  • The small molecule databases for Docking or Virtual Screening have been put in /opt/database directory
     ChemBridge: /opt/database/ChemBridge
     1. EXPRESS-Pick™ Database: EXPRESS-Pick October 2008.sdf  or EXPRESS-Pick October 2008.db
     2. KINASet Database (ATP analogue): KINASet.sdf or KINASet.db
     This is the same as kinaset in CHembridge/ KINASet@: New_custom_Kinaset_3D_05072007 as a sdf or db or zip file (/opt/database/chembridge/kinaset)
     3. NOVACore/PHARMACophore Combined Database: NCL-PCL October 2008.sdf or NCL-PCL October 2008.db
     4. GPCR and KINACore Library databases: GPCR September2008.zip
     5. FOCUSCore Library database: FocusCore_September2008.zip
     Diversity set (non ATP analogue)
     6. NCI(National Cancer Institute)-diversity Set: NCI-Open_09-03.sdf.gz(/opt/database/NCI-diversity)
     7. Maybridge Chemical Company (England): MayBridge.sdf (/opt/database/maybridge)


  • pdbbind opt/reprints/benchmark/pdbbind/all-core-10-2008.xls
  • binding DB (Gilson)

Small Molecule Preparation

If you need to convert the structure from 2D to 3D, add hydrogen and optimize. 1) You can use babel to add hydrogen’s, generate 3D structure, and optimize. An example is in /opt/database/chembridge/DIVERSet-06/2d-3d-test/ 2) Alternatively, you can use openeye software to do this. Chunli can help you with this. 3) Or use my script in /op/database/sdfbatch.pl The first option is probably the best. Let me know if you have questions.







Updated license info can be found :




GOLD Tutorial

Before running GOLD from the command line, it is recommended to follow a tutorial using the graphical interface to gain familarity with how GOLD works.


A tutorial for docking two ligands into a protein while taking flexibility into consideration is located here:


A tutorial for ligand binding pose prediction is located here:



Remove old env in your bashrc or cshrc if you have one (so far I don't think we need it for 2020 version)

Activate GOLD the first time you use it:

/opt/CCDC/CCDC2020/CSD_2020/bin/ccdc_activator -k 570769-C9BD90-4986BB-904E59-D56C9B-025394

If you want to use gradphical interface for activating or setting up GOLD job initially (running job requires no GUI), you need a X grapics on your local computer.

Mac has Xquartz.  On windows, use mobaxterm that comes with X (or you can install VcXsrv & use it with ubuntu). Make sure you use ssh -Y when logging into a node.

You can check that X is working by typing xeyes and hitting enter. You should get a pair of googly eyes (which may be in another window).

/opt/CCDC/CCDC2020/CSD_2020/bin/ccdc_activator_gui & enter key from above (570769-C9BD90-4986BB-904E59-D56C9B-025394)

Running graphical interface (e.g. to generate gold.conf)

ssh -Y nova



Once you set up the job/create the conf file on nova, you can run the job on other nodes (below)

GOLD user guide: https://www.ccdc.cam.ac.uk/support-and-resources/ccdcresources/GOLD_User_Guide.pdf 

Running GOLD from the command line (you can use gold.conf above or mannually creat/edit conf)

/opt/CCDC/CCDC2020/Discovery_2020/bingold_auto gold.conf & 

2018 update

(installation and activation same as 2016)

setenv GOLD_DIR "/opt/CCDC2018/GoldSuite_2018/GOLD"
setenv PVM_ROOT "/opt/CCDC2018/GoldSuite_2018/GOLD/pvm3"
setenv PVM_RSH "/usr/bin/ssh"
setenv CCDC_CSD_LICENCE_FILE "/opt/CCDC2018/CSD_2018/csd/csd_license.dat"
setenv LD_LIBRARY_PATH /opt/CCDC2018/GoldSuite_2018/c_linux-64/lib/:$LD_LIBRARY_PATH

To activate on a node or WS: /opt/CCDC2018/GoldSuite_2018/bin/gold,  Enter site 771,  Confirmation Code: 1f9986 and email.

After installation, you may run the GUI (hermes) on workstation such as water.bme.utexas.edu


Installation: /opt/CCDC/installation/csds-linux-x64.run  Mannual in /op/CCDC

Activate on a node (e.g. node101) by running /"'opt/CCDC/GoldSuite_2016/bin/gold" (using X windows if you do so remotely). Windows will pop up, asking you to enter the site number (771) and purchase number/confirmation code to register online. This confirmation number changes every spring each year (1f9986 for 2018). The code is now available from x-ray lab web site: https://sites.utexas.edu/chem-x-ray/ (right panel)

Once activated, specify the license file: csd_licence.dat /opt/CCDC//CSD_2016/csd/csd_licence.dat 

To run:

setenv GOLD_DIR "/opt/CCDC/GoldSuite_2016/GOLD"
setenv PVM_ROOT "/opt/CCDC/GoldSuite_2016/GOLD/pvm3"
setenv PVM_RSH "/usr/bin/ssh"
setenv CCDC_CSD_LICENCE_FILE "/opt/CCDC/CSD_2016/csd/csd_licence.dat"
setenv LD_LIBRARY_PATH /opt/CCDC/GoldSuite_2016/c_linux-64/lib/:$LD_LIBRARY_PATH
Then start pvm with a host file (first killall pvmd3 on all nodes you will be using and cleam /tmp/pvm*)
/opt/CCDC/GoldSuite_2016/GOLD/pvm3/lib/LINUX64/pvmd3 hostfile
Hostfile looks like this:
Then (assuming gold.conf exists)
/opt/CCDC/GoldSuite_2016/GOLD/bin/parallel_gold_auto  16
You need another gold.hosts file:
node101.bme.utexas.edu 8
no_of_processes 8
node36.bme.utexas.edu 8
no_of_processes 8

 Notes from Dave Broadwell

CSDS licensing, which GOLD now also uses, is node-based - this means that each individual node that GOLD (or any of the CSD software) is to run on must be registered.
You installation and registration on one of your cluster nodes has enabled use on that particular node - in order to run on other nodes they similarly will need to be registered. That can be achieved by running one of the programs (other than GOLD itself - you will need to run one of the programs with a GUI such as ConQuest, Mogul or Hermes), or if you need to register multiple nodes and/or need to do so at the command line we can provide a tool to help with that.
Once registered then hopefully your PVM use should work. I should note, however, that use of PVM to parallelise GOLD is no no longer our recommended method of doing so and PVM will become deprecated for use with GOLD at some point in a future release. There are some entries in our support database for PVM use - see:
In general, we now recommend use of a system such as univa grid and a script to split GOLD jobs and use the grid's scheduling system to run them - for an example script see:
and in particular:


You can split the jobs and combine the results by uinsg the script discussed in the link above (Once all the jobs have been launched, a script is generated that may be used to collate all the results once the dockings have completed. )

Gold Tutorial:

All scripts and example gold.hosts and gold.conf files are available here

First, you need to initialize a file structure for docking. run the makedirs.sh script, making sure that the gold.hosts file( explained above) exists in the current working directory, as well as the gold.conf. The gold.conf file should not contain a ligand line, this is created by the makedirs.sh script.

Secondly, you need to run makerunscript.sh with an argument equal to the total number of cores, This will create a file run.sh that you will run to start job execution

Start up the master daemon on node36 by entering nohup /opt/CCDC/GoldSuite_2016/GOLD/pvm3/lib/pvmd hostfile &, with hostffile being the hostfile descripted above( not the gold.hosts)

To check that the daemon is active, run /opt/CCDC/GoldSuite_2016/GOLD/pvm3/lib/pvm and enter conf in the resulting command prompt. You should see a list of all the hosts provided to docking.

run nohup ./run.sh. This starts docking. Expect docking to take 1-1.5 weeks

to analyze, cat */bestranking.lst to a new file. This collects all ligand scores.

I have provided a script( extracthigh.py) that takes all lines in the collected bestranking file that score above a limit. To use run python extracthigh.py inputfile.txt outputfile.txt minimum score( ex. 90).

revised TI3 usage

The queuing script is ineffecent. Instead, run jobs on each individual node. example code:

for dir in */; do

       cd $dir
   echo $cwd
   sleep 30
   while $check == 0; do
           if  cut -f2| sed -n 1p|awk '{printf "%.0f\n", $1}') -lt $limit ; then
                   ssh compute-0-9 /apps/gold/goldsuite-5.2/bin/gold_auto $cwd/gold.conf
                   echo "9"
           elif  cut -f2| sed -n 2p|awk '{printf "%.0f\n", $1}') -lt $limit ; then
                           ssh compute-0-3 /apps/gold/goldsuite-5.2/bin/gold_auto $cwd/gold.conf
                           echo "3"
               	sleep 360

cd ../ done

notes on gold inputs

Gold will not accept ligand files written in sdf format. They need to be mol2s

For gold.conf file, make sure that path is full, don't use ~ shortcut.


<a href="Research%3ADrugDesign">autodock 4</a>


  • Start SSH
  • Create a public and a private key.
    • Click "Generate New"
    • Change DSA to RSA
    • Private Key Stats
      • Filename: id_rsa (default, but can put whatever in here)
      • Passphase: (Optional, used to access private key)
  • Put the key into the computers that are to be remotely accessed.
  • After putting the keys into those other comptuers, click "Configure"
  • Click on "Profiles" in the console, and enter name of computer.
    • bme-(planetname).bme.utexas.edu (All planets but Mercury)
  • To copy the public key, open console, go to the home directory, and then click the .ssh directory.
    • Open a file called authorized_keys.
    • Copy the public key to that file, make sure it is 1 line. Be sure to add "ssh-rsa " to the front of the key.


TACC portal

Submit a job on lonestar

ssh username@lonestar.tacc.utexas.edu

showq -u #check your job

Put yoyr job in work directory: cdw or cd $WORK

bsub < job_script &

example script  /home1/00433/prentacc/EXAMPLE/run

To use TINKER:

    Before run tinker job, do "module load mkl", or add to your source file.

    Regular tinker binary: /home1/00433/prentacc/tinker_6.2.06_tacc/source

    Osrw tinker bianry: /home1/00433/prentacc/tinker-osrw/source-osrw-2-17s

Example job_script (change BOLD part for each job):

#$ -V       #Inherit the submission environment
#$ -cwd   # Start job in submission directory
#$ -N solv3      # Job Name
#$ -A  TG-MCB100057    # or A-biomod 
#$ -j y     # Combine stderr and stdout
#$ -o $JOB_NAME.o$JOB_ID         # Name of the output file (eg. myMPI.oJobID)
#$ -pe 12way 12  # Requests 12 tasks/node, 12 cores total. Increase the second number for MPI jobs.
#$ -q normal     # Queue name normal
#$ -l h_rt=24:00:00      # Run time (hh:mm:ss) 12 hr max
#$ -M pren@mail.utexas.edu       # Address for email notification
#$ -m be         # Email at Begin and End of job #ibrun ./a.out  # Run the MPI executable named a.out
module load mkl
/home1/00433/prentacc/tinker-osrw/source-osrw-2-17s/dynamic.x solv03-100-100.xyz_2  1000000 2.0 5 2 298

Submit a series of jobs in TACC

you can create a series of jobs that are dependent on previous steps being completed in lonestar or ranger. I found this very useful when you have a few steps but you can submit them at the same time; or when you have a job that would last for a few days but got killed every 24 hours. By defining the dependency between jobs, you don't have to resubmit your jobs manually everyday. In addition, it saves your waiting time a lot because your jobs will be automatically in the waiting list whenever the previous ones are finished. To specify a job dependency, use "-hold_jid" tag in your qsub job script.

login2$ qsub -help | grep hold

  [-hold_jid job_identifier_list]          define jobnet interdependencies
[-hold_jid_ad job_identifier_list]       define jobnet array interdependencies

For example, if you want to run a job first and you have just submitted it, which has a job id of"67121". Assume you have your restart input files in place for the next step, you can specify "-hold_jid" tag in your qsub job script "job2", for example in the following case ("#$ -hold_jid 67121"). Then you can submit the second job right away to the scheduler: qsub < job2. The second job will wait until the first job is finished (or killed) to begin execution or to be placed in the waiting queue. You can also have job3 (day3), job4 (day4), job5 (day5).... which are dependent on the previous ones.

#$ -V
#$ -cwd
#$ -N min3
#$ -A  A-biomod
#$ -j y
#$ -o alphA
#$ -pe 12way 24
#$ -q development
#$ -l h_rt=01:00:00
#$ -m be
#$ -hold_jid 67121



<a href="Research%3ARemd">REMD HOWTO</a>

GROMACS tutorial

<a href="http://www.dddc.ac.cn/embo04/#lecture">http://www.dddc.ac.cn/embo04/#lecture</a> IED for visualize major mode of motion: <a href="http://mccammon.ucsd.edu/ied/">http://mccammon.ucsd.edu/ied/</a>

creating custom residue

Pyrrole as an example: /home/pren/gromacs-top-pyr/readme-pyr

Only in /home/pren/gromacs-top-pyr: backupffoplsaa.rtp
diff top/ffoplsaa.atp /home/pren/gromacs-top-pyr/ffoplsaa.atp
<  opls_542   14.00670  ; N   in pyrrole
>  opls_542   11.00670  ; N   in pyrrole
<  opls_544   12.01100  ; C3  in pyrrole
<  opls_545    1.00800  ; H1  in pyrrole
>  opls_544   9.01100  ; C3  in pyrrole
>  opls_545    4.00800  ; H1  in pyrrole
<  opls_547    1.00800  ; H3  in pyrrole
>  opls_547    4.00800  ; H3  in pyrrole
diff top/ffoplsaabon.itp /home/pren/gromacs-top-pyr/ffoplsaabon.itp
>   CW     NA     H       1   120.000    292.880   ; 
<   NA     CW     CS      1   107.700    585.760   ; wj/nm
>   NA     CW     CS      1   107.700    785.760   ; wj/nm
>   CW     CW     CS      1   107.300    585.760   ; -idem-
>   CS     CW     NA     H       3     00.00200   0.00000 -12.55200   0.00000   0.00000   0.00000 ; TRP chi-4
<   X      CW     CW     X       3     44.97800   0.00000 -44.97800   0.00000   0.00000   0.00000 ; 
> ;  X      CW     CW     X       3     44.97800   0.00000 -44.97800   0.00000   0.00000   0.00000 ; 
>   NA     CW     CW     NA       3     6.3500   -2.450000 -2.01000   0.00000   0.00000   0.00000 ; 
>   CS     CW     CW     NA       3     0.00000   0.000000 0.00000  0.00000   0.00000   0.00000 ; 
>   CS     CW     CW     CS       3     0.00000   0.000000 0.00000  0.00000   0.00000   0.00000 ; 
diff top/ffoplsaanb.itp /home/pren/gromacs-top-pyr/ffoplsaanb.itp
<  opls_542   NA     714.00670    -0.239       A    3.25000e-01  7.11280e-01
>  opls_542   NA     711.00670    -0.239       A    3.25000e-01  7.11280e-01
<  opls_544   CS     612.01100    -0.149       A    3.55000e-01  2.92880e-01
<  opls_545   H      1 1.00800     0.317       A    0.00000e+00  0.00000e+00
>  opls_544   CS     69.01100    -0.149       A    3.55000e-01  2.92880e-01
>  opls_545   H      1 4.00800     0.317       A    0.00000e+00  0.00000e+00
<  opls_547   HA     1 1.00800     0.118       A    2.42000e-01  1.25520e-01
>  opls_547   HA     1 4.00800     0.118       A    2.42000e-01  1.25520e-01
diff top/ffoplsaa.rtp /home/pren/gromacs-top-pyr/ffoplsaa.rtp
<    CG2    CB    CA     C    dih_VAL_chi1_C_C_C_CO
>    GCG2    CB    CA     C    dih_VAL_chi1_C_C_C_CO
> ;up uncharged pyrrole
> [ PYU ]
>  [ atoms ]
>    N1P   opls_542    -0.214    1
>    C1P   opls_543    -0.036    1
>    C2P   opls_543    -0.036    1
>    C3P   opls_544    -0.152    1
>    C4P   opls_544    -0.152    1
>    HNP   opls_545     0.210    1
>    H3P   opls_547     0.190    2
>    H4P   opls_547     0.190    2
>  [ bonds ]
>    N1P   HNP
>    C1P   N1P
>    C2P   N1P
>    C3P   C1P
>    C2P   C4P
>    C3P   C4P
>    H3P   C3P
>    H4P   C4P
>   -C2P   C1P
>  [ impropers ]
>   -C2P   C1P   N1P    C3P    improper_Z_CA_X_Y 
>    C4P   N1P   C2P    C1P    improper_Z_CA_X_Y 
>    C3P   N1P   C1P    C2P    improper_Z_CA_X_Y 
>    H3P   C3P   N1P    C4P    improper_Z_CA_X_Y 
>    H4P   C4P   N1P    C3P    improper_Z_CA_X_Y 
>    HNP   N1P   C4P    C1P    improper_Z_CA_X_Y 
>  ;  C2P   C1P  -C2P    C3P    improper_Z_CA_X_Y
> ; charged pyrrole 
> ; [ atoms ] reduce charge becuase of 1-4 interaction
> ;   N1P   opls_542    -0.380    1
> ;   C1P   opls_543     0.422    1
> ;   C2P   opls_543     0.422    1
> ;   C3P   opls_544    -0.120    1
> ;   C4P   opls_544    -0.120    1
> ;   HNP   opls_545     0.396    1
> ;   H3P   opls_547     0.190    2
> ;   H4P   opls_547     0.190    2
> [ PYC ]
>  [ atoms ]
>    N1P   opls_542    -0.380    1
>    C1P   opls_543     0.422    1
>    C2P   opls_543     0.422    1
>    C3P   opls_544    -0.120    1
>    C4P   opls_544    -0.120    1
>    HNP   opls_545     0.396    1
>    H3P   opls_547     0.190    2
>    H4P   opls_547     0.190    2
>  [ bonds ]
>    N1P   HNP
>    C1P   N1P
>    C2P   N1P
>    C3P   C1P
>    C2P   C4P
>    C3P   C4P
>    H3P   C3P
>    H4P   C4P
>   -C2P   C1P
>  [ impropers ]
>   -C2P   C1P   N1P    C3P    improper_Z_CA_X_Y 
>    C4P   N1P   C2P    C1P    improper_Z_CA_X_Y 
>    C3P   N1P   C1P    C2P    improper_Z_CA_X_Y 
>    H3P   C3P   N1P    C4P    improper_Z_CA_X_Y 
>    H4P   C4P   N1P    C3P    improper_Z_CA_X_Y 
>    HNP   N1P   C4P    C1P    improper_Z_CA_X_Y 
> ;   C2P   C1P  -C2P    C3P    improper_Z_CA_X_Y
Only in /home/pren/gromacs-top-pyr: readme-pyr
torsion of backbone see: /home/pren/pyr/QM/torsion.xls

AMBER Tutorial


Basic MD

set up parameters for new ligands and unnatural resiude

MM-PBSA for binding: http://ambermd.org/tutorials/advanced/tutorial3/


see this page for installation, nodes and additional information  AMBER-GPU

To run AMBER-GPU, first include in your bashrc file:

source /opt/intel/composer_xe_2013.2.146/bin/compilervars.sh intel64

export PATH=$PATH:/usr/local/cuda-6.5/bin
export LD_LIBRARY_PATH=/usr/local/cuda-6.5/lib64:$LD_LIBRARY_PATH
export AMBERHOME=/opt/software/amber14-gpu

alias python=/usr/bin/env\ python2.6
export PYTHONPATH=/usr/bin/python2.6
export PATH=/opt/mpich2-1.4.1p1/bin/mpicc/:/opt/mpich2-1.4.1p1/bin:$PATH
export CUDA_HOME=/usr/local/cuda-6.5

export DO_PARALLEL='mpirun -np 1'

The export DO_PARALLEL statement tells amber how many gpus to use. We have three gpus in sugar, but in order not to double up one gpu, just use one gpu per job.

You can find out if the gpu's are being used by the command: nvidia-smi (assuming you have the above variables in you bashrc/equivalent file) In order to run on a particular gpu, just set the environment variable:

export CUDA_VISIBLE_DEVICES="#"  ;where # is the gpu number (0,1,2 for sugar)

For node105 (possibly 106 and 107 too) use /opt/software/amber14-gpu-earth/ rather than /opt/software/amber14-gpu/

The executables can be found in: /opt/software/amber14-gpu/bin/

An example executable command line is given below(taken from /home/davbell/software/amber/test/cuda/dhfr/Run.dhfr)

$DO_PARALLEL $sander -O -o $output -c md12.x -r restrt -x mdcrd < dummy || goto error

Further tutorials can be found through the AMBER manual: http://ambermd.org/doc12/Amber14.pdf and on the website http://ambermd.org/gpus/



MM-PBSA free energy estimation tutorial

Author: Chengwen Liu
Latest update: Dec.5th, 2017

a. Theory preparation

A very good tutorial can be found in http://ambermd.org/tutorials/advanced/tutorial3/. I will suggest you to first read this tutorial so you can get a basic understanding of MM-PB/GBSA. The following steps will guide you to b) set up the complex system, c) run MD simulations on our GPU cards using pmemd.cuda, and d) estimate the complex binding affinity.

b. System setting up 

System setting up in AMBER can be done in either xleap or tleap program. The former provides a graphical interface thus it is easier for the beginners. Here we will use tleap to set the system up (generate the necessary files for AMBER run).

Copy the leaprc.ff14SB (this is the force field file we will be using) to the current working directory as leaprc.

 cp $AMBERHOME/dat/leap/cmd/leaprc.ff14SB ./leaprc

 Append the following commands to the leaprc file.

Notes: these commands are helping us to generate the complex prmtop and rst7 files in gas phase (COMP.prmtop & COMP.rst7), add Na+ and Cl- ions to the salt concentration we wanted (addions commands), add solvent molecules (TIP3P water here) and generate complex prmtop and rst7 files in explicit water (COMP_wat.prmtop & COMP_wat.rst7), generate the ligand and protein prmtop and rst7 files (LIG.prmtop & LIG.rst7; PRO.prmtop and PRO.rst7)

 COMP = loadpdb COMP.pdb
 saveamberparm COMP COMP.prmtop COMP.rst7
 addions COMP Na+ 0
 addions COMP Cl- 0
 addions COMP Cl-  number (see notes *)
 addions COMP Na+ number (see notes *)
 solvateoct COMP TIP3PBOX 8.0
 saveamberparm COMP COMP_wat.prmtop COMP_wat.rst7
 savepdb COMP COMP_wat.pdb
 LIG = loadpdb LIG.pdb
 saveamberparm LIG LIG.prmtop LIG.rst7
 PRO = loadpdb PRO.pdb
 saveamberparm PRO PRO.prmtop PRO.rst7

Execute tleap command to generate the bunch of prmtop and rst7 files


After we neutralize the system, the number of ions needed to reach the salts concentration should be calculated.*

The above steps help us generate one individual system for the next runs. Here is a python script to automatically produce multiple systems. Modifications are very possibly needed for your own cases.

#!/usr/bin/env python

#set up AMBER input files for use in MM-PBSA free energy calculations

# files need to be prepared:
# 1. leaprc file
# 2. complex pdb file

#Chengwen Liu

import sys,os 

def prepareInputs(complexPdbFile):
  lines = file(complexPdbFile).readlines()
  filename = complexPdbFile.split(".pdb")[0]
  f1 = open("Integrin_%s.pdb"%filename, "w")
  if "I1" in complexPdbFile:
    nAtoms = 1557
  if "I2" in complexPdbFile:
    nAtoms = 1544
  for n in xrange(nAtoms):

  f2 = open("Collagen_%s.pdb"%filename, "w")
  for n in xrange(nAtoms, len(lines)):
  f3 = open("temp", "w")
  f3.write("COMP = loadpdb %s\n"%complexPdbFile)
  f3.write("saveamberparm COMP %s.prmtop %s.rst7\n"%(filename, filename))
  f3.write("addions COMP Na+ 0\n")
  f3.write("addions COMP Cl- 0\n")
  f3.write("solvateoct COMP TIP3PBOX 8.0\n")
  f3.write("saveamberparm COMP %s_wat.prmtop %s_wat.rst7\n"%(filename, filename))
  f3.write("savepdb COMP %s_wat.pdb\n"%filename)
  f3.write("LIG = loadpdb Collagen_%s.pdb\n"%filename)
  f3.write("saveamberparm LIG Collagen_%s.prmtop Collagen_%s.rst7\n"%(filename, filename))
  f3.write("PRO = loadpdb Integrin_%s.pdb\n"%filename)
  f3.write("saveamberparm PRO Integrin_%s.prmtop Integrin_%s.rst7\n"%(filename, filename))
  os.system("cat leap temp >leaprc") 
  #do the HMassRepartition and generate a new prmtop file 
  f4 = open("parmed.in", "w")
  f4.write("parmout %s_wat_heavyH.prmtop\n"%filename)
  os.system("/opt/software/ParmEd-master/parmed.py -p %s_wat.prmtop -i parmed.in"%filename)

currDir = os.getcwd()
files = os.listdir(currDir)
for eachPdb in files:
  if ".pdb" in eachPdb:
    #print eachPdb 
    filename = eachPdb
    os.system("mkdir %s/%s"%(currDir,filename.split(".pdb")[0]))
    os.system("cp %s/leap %s/%s"%(currDir, currDir,filename.split(".pdb")[0]))
    os.chdir("%s/%s"%(currDir, filename.split(".pdb")[0]))
    os.system("cp %s/%s ."%(currDir, filename))



c. Run MD simulations with pmemd.cuda

 Routinely in AMBER, we will run minimization, heating, equilibration, cooling, equilibration and final production run in NVT ensemble. Here are the input and run files I am using in collagen-integrin system.


initial minimization prior to MD 
  imin   = 1,
  maxcyc = 6000,
  ncyc   = 2500,
  cut    = 12


200ps MD 
  imin   = 0,
  irest  = 0,
  ntx    = 1,
  ntb    = 1,
  cut    = 10.0,
  ntc    = 2,
  ntf    = 2,
  tempi  = 0.0,
  temp0  = 360.0,
  ntt    = 3,
  gamma_ln = 1.0,
  nstlim = 100000, dt = 0.002,
  ioutfm = 1,
  ntpr = 250, ntwx = 250, ntwr = 10000


200ps of MD
  imin = 0, irest = 1, ntx = 7,
  ntb = 2, pres0 = 1.0, ntp = 1,
  taup = 2.0,
  cut = 10.0, ntr = 0,
  ntc = 2, ntf = 2,
  tempi = 360.0, temp0 = 360.0,
  ntt = 3, gamma_ln = 1.0,
  nstlim = 100000, dt = 0.002,
  ioutfm = 1,
  ntpr = 250, ntwx = 250, ntwr = 10000


200ps MD 
  imin   = 0,
  irest  = 0,
  ntx    = 1,
  ntb    = 1,
  cut    = 10.0,
  ntc    = 2,
  ntf    = 2,
  tempi  = 360.0,
  temp0  = 300.0,
  ntt    = 3,
  gamma_ln = 1.0,
  nstlim = 100000, dt = 0.002,
  ioutfm = 1,
  ntpr = 250, ntwx = 250, ntwr = 10000


200ps of MD
  imin = 0, irest = 1, ntx = 7,
  ntb = 2, pres0 = 1.0, ntp = 1,
  taup = 2.0,
  cut = 10.0, ntr = 0,
  ntc = 2, ntf = 2,
  tempi = 300.0, temp0 = 300.0,
  ntt = 3, gamma_ln = 1.0,
  nstlim = 100000, dt = 0.002,
  ioutfm = 1,
  ntpr = 250, ntwx = 250, ntwr = 10000


100ns MD NVT 
  imin   = 0,
  irest  = 0,
  ntx    = 1,
  ntb    = 1,
  cut    = 10.0,
  ntc    = 2,
  ntf    = 2,
  tempi  = 300.0,
  temp0  = 300.0,
  ntt    = 3,
  gamma_ln = 1.0,
  nstlim = 25000000, dt = 0.004,
  ioutfm = 1,
  ntpr = 1000, ntwx = 1000, ntwr = 1000

Run script:

export filename=your file name

source ~pren/amber-git/amber17/bashrc.amber16 
export CUDA_VISIBLE_DEVICES="0" #or 1

$AMBERHOME/bin/pmemd.cuda -O -i min.in -o $filename\_min.out -c $filename\.rst7 -p $filename\.prmtop -r $filename\_min.rst7 
$AMBERHOME/bin/pmemd.cuda -O -i heat360.in -o $filename\_heat360.out -c $filename\_min.rst7 -p $filename\.prmtop -r $filename\_heat360.rst7 -x $filename\_heat360.nc 
$AMBERHOME/bin/pmemd.cuda -O -i md200ps360.in -o $filename\_md200ps360.out -c $filename\_heat360.rst7 -p $filename\.prmtop -r $filename\_md200ps360.rst7 -x $filename\_md200ps360.nc 
$AMBERHOME/bin/pmemd.cuda -O -i cool300.in -o $filename\_cool300.out -c $filename\_md200ps360.rst7 -p $filename\.prmtop -r $filename\_cool300.rst7 -x $filename\_cool300.nc 
$AMBERHOME/bin/pmemd.cuda -O -i md200ps300.in -o $filename\_md200ps300.out -c $filename\_cool300.rst7 -p $filename\.prmtop -r $filename\_md200ps300.rst7 -x $filename\_md200ps300.nc 

$AMBERHOME/bin/pmemd.cuda -O -i mdprod.in -o $filename\_md1.out -c $filename\_md200ps300.rst7 -p $filename\_heavyH.prmtop -r $filename\_md1.rst7 -x $filename\_md1.nc 

After the production run, we may want to delete the explicit water molecules in the trajectory file (*.mdcrd or *.nc), which helps us to save the hard disk storage. In this step the ptraj (or cpptraj) program will be used. 

Here is the trajin file

 parm COMP_wat.prmtop
 trajin COMP_wat_md1.nc
 strip :WAT:Na+:Cl-
 trajout COMP_wat_md1_traj.nc
Here is run command
cpptraj -i trajin

d. MM-PB/GBSA calculation

We will use MMPBSA.py to estimate the binding affinity, where the entropy contribution is not computed due to the expensive cost of normal mode analysis calculations for biomolecules. (See theory preparation section tutorials for more details)

pbsa.in file

MMPBSA from JCTC paper -CW
  startframe=1, endframe=25000,
  keep_files=0, interval=50, debug_printlevel=2,
   igb=2, saltcon=0.100
  istrng=0.1, exdi=80, indi=1.0,
  inp=2, cavity_surften=0.0378, cavity_offset=0.5692,
  fillratio=4, scale=2.0,
  linit=1000, prbrad=1.4, radiopt=1,

Run file

export filename=your file name
source  /home/pren/amber-git/amber17/ambercpu/amber.sh
$AMBERHOME/bin/MMPBSA.py -O -i pbsa.in -o results_mmpbsa.dat -sp $filename\.prmtop -cp $filename\.prmtop  -rp Integrin_$filename\.prmtop


e. MM-PB/GBSA results

see a. theory preparation for how to read/understand the MMPBSA results.


CHARMM tutorial

<a href="http://www.ch.embnet.org/MD_tutorial/">http://www.ch.embnet.org/MD_tutorial/</a> <a href="http://www.charmming.org/charmming/">http://www.charmming.org/charmming/</a> <a href="http://www.charmm.org/ubbthreads/ubbthreads.php?ubb=cfrm">http://www.charmm.org/ubbthreads/ubbthreads.php?ubb=cfrm</a>

Movie making

You can use trjconv to convert groamcs trajectory to a xxx.pdb (many frames), and then use the pdb to make movies. One option for trjconv is –pbc, which can be used to keep the molecules in the box. One of these may work for you:

  • -pbc mol. This put the center of mass of molecules in the box.
  • -pbc nojump. Here you use the first frame (-s xxx.tpr) use as reference, and get rid of any jumps out of the box. This only works if the firs frame is whole in the box.

Amber has similar utilities to convert trajectory to pdb

To make movies look slower, you can more frames 9smaller steps between frames), decrease the no of frames per second or increase delay in animated gif if you use unfreez under windows…

  • use vmd movie maker. You can decrease frames per seconds (default is 24 the lowest already; see Format/compression setting).
<a href="http://www.ks.uiuc.edu/Training/Tutorials/vmd/tutorial-html/node3.html#SECTION00034200000000000000">http://www.ks.uiuc.edu/Training/Tutorials/vmd/tutorial-html/node3.html#SECTION00034200000000000000</a> 
<a href="http://ambermd.org/tutorials/basic/tutorial2/section9.htm">http://ambermd.org/tutorials/basic/tutorial2/section9.htm</a>

to turn on ray tracing for each png file:
  viewport 320,240 9decide how big the picture is)
set ray_trace_frames=1
set cache_frames=0
mpng mymov (or go to File/export movie)

Additional pymol commands;<a href="http://pymol.sourceforge.net/newman/user/S0300movies.html">http://pymol.sourceforge.net/newman/user/S0300movies.html</a>

mset 1 x30      # creates a 30 frame movie consisting of state 1 played 30 times.
mset 1 -30      # creates a 30 frame movie: states 1, 2, ... up through 30.
mset 1 -30 -2   # 58 frames: states 1, 2, ... , 29, 30, then 29, 28, ... , 4, 3, down to 2
mset 1 6 5 2 3  # 5 frames: states 1, 6, 5, 2, 3 in that order.



Discovery Studio & Pipeline Pilot

Latest version

DS 3.5
PP 8.5.3



# network access problem from outside lab

The default route ( ip route) shows was the default gateway, in addition to

Use "route del default" twice to remove all the default and then add using "/sbin/route add default gw eth0"

Also removed gateway from ifcfg-eth1.

  • PP installation dir on pharma is /home/qwang/software/Accelrys/PP_install
Web acess: pharma.bme.utexas.edu:9943 (https)
   pharma.bme.utexas.edu:9944 (http)
Admin Accounts: scitegicadmin and/or root (same password)


  • DS installation is done through the PP installation script. A single license file contains both PP and DS is required to install both of them in a integrated installation.

A such license file looks like:

INCREMENT License_Holder msi 2016.120 30-dec-2016 uncounted 62B406C87B7EAC123E26 VENDOR_STRING="200531:218555:5:University of Texas at Austin" HOSTID=ANY TS_OK
INCREMENT scitegic_ppclient msi 2016.120 30-dec-2016 uncounted A2949618D29B97720F7F VENDOR_STRING=599622[25854] HOSTID=ANY TS_OK
INCREMENT DSCollection msi 2016.120 30-dec-2016 uncounted 42B4168820BAED55C170 VENDOR_STRING=614772[25554] HOSTID=ANY TS_OK
INCREMENT DSDeveloperClient msi 2016.120 30-dec-2016 uncounted 527406E85FF8986EF951 VENDOR_STRING=614772[25554] HOSTID=ANY TS_OK
INCREMENT scitegic_core msi 2016.120 30-dec-2016 uncounted 82C4C62864E053BDFFDF VENDOR_STRING=614772[25554] HOSTID=ANY TS_OK

Once the installation is done, there should have two more new directories created (DiscoveryStudio35 and LicensePack) at the same level of your PP install dir (PP_install in our case).


  • Then, another license file is required to activate the DS component. A request form like below sending to Accelrys can get us this file.
COMPANY: <your company's name>
CONTACT: <your name>
EMAIL:   <email address to which license file should be sent>

Enter the Accelrys products for which you are requesting a license.
(include order number if available)
<Product Name> <Order Number>

You should not need to edit any of the remaining portion of the form.

COMPUTER:       pharma.bme.utexas.edu
HOST ID:        MAC address
OS TYPE:        Linux

End Accelrys Software License File Request Form

Such a license file should looks like this:

SERVER pharma.bme.utexas.edu MACaddress 1715 (This line defines the port and server name that we can used to connect to in client machines)
VENDOR msi PORT=1720 (NB: This port number is crucial for client in connecting to the server.)
             (    It is used for the license managing program to communicate between client and server.)
             (    So we need to make sure it is open! The DS installation manual does NOT even mention it.)
             (    We have been stuck at this point for quite a while.)
INCREMENT License_Holder msi 2016.120 30-DEC-2016 999 2D42D7207F6D64AE0EBC "200531:218557:1:University of Texas at Austin"
INCREMENT DS_ADMET_User_MP msi 2016.120 30-DEC-2016 1 ADB287204B698029EC64 "599619[28071]"
INCREMENT DS_MCSS_User_MP msi 2016.120 30-DEC-2016 1 9DC29750CCAE2EF8B007 "599619[28071]"


  • It is important to make sure the system has the older ethX network interface card naming convention rather than the newer emX names. Otherwise, the DS installer will not smart enough to extract the MAC address of the server, thus the license server won't work. Dell machines and Fedora15 or newer usually follow the emX naming scheme, so they need to be changed to ethX. A method that works for our pharma server can be found at: <a href="http://www.sysarchitects.com/em1_to_eth0">http://www.sysarchitects.com/em1_to_eth0</a>
1. Add biosdevname=0 to the kernel boot arguments in /etc/grub.conf.
2. Rename /etc/sysconfig/network-scripts/ifcfg-em1 to /etc/sysconfig/network-scripts/ifcfg-eth0, changing the line
3. Delete /etc/udev/rules.d/70-persistent-net.rules
4. Reboot


  • DS license server

Sometimes, the license server may go down. To restart, go to /home/qwang/software/Accelrys/LicensePack/linux/bin

And run ./lp_server -s to start the server.
            -h for help
            -w for status


Once above things are done, we can download and install a DS client in a linux or windows machine from our server's website by following the following instructions:

  • Go to pharma.bme.utexas.edu:9944, and click on the Discovery Studio link at the right bottom corner. This will bring your to the DS website. Then you can simply click on the Install/Upgrade Discovery Studio Client link to install either windows or linux version. This will install a free version of DS client.


  • To enable the full version, we can click on the 'Enable additional features' button at the right bottom corner of the client window. Then choose connect to a license server option, and enter the port number 1715 and the server name pharma.bme.utexas.edu to appropriate boxes. Then you should be able to tell whether it is connected or not from the window. And the 'Enable additional features' button should have been changed to 'Server: none' now. You would also notice there are much more options appear (but still are grayed menu items) in the left panel of each module.


  • If you click on these grayed options, there should be a pop-up window asking for server name. Just enter pharma.bme.utexas.edu, and then use your username and password to login. Now, you should have the full power of the software!

Schrodinger and GLIDE Docking

We have 20 perpetual Glide licenses on Ti3D cluster; 2 perpetual ligprep license. (see Box Sync/LAB/software/ti3d glide)

We have 50 annual GLIDE licenses on Ren cluster and expires in 2015?

Prepping Protein and Grid

Preparation of the protein and grid can be accomplished using the maestro interface, available at at /opt/Schrodinger_Suites_2013-2/maestro. protein preparation and grid preparation instructions are available at https://www.schrodinger.com//AcrobatFile.php?type=supportdocs&type2=&ident=252. The final output of the grid generation is a .zip file. This file contains protein coordinates as well as other data. No other protein input is necessary.

Ligand Preparation

The ZINC server has a discrepancy in the size of the drug_lime in stock library depending on if one selects .mol2 or .sdf output. The latest molecules appear to be only in the .sdf output. To download the drug-like in stock library, download the file under UNIX SDF at http://zinc.docking.org/subsets/drugs-now. In terminal, run this file using csh. This will download the library as many .sdf files, which can be decompressed using gunzip. Each .sdf file corresponds to part of the library.

Do not attempt to run ligprep using maestro or the bme-nova nodes. We are not licenced for ligprep. To run ligprep, you need to connect to a remote server in Kong Ren's lab( see Pengyu for information). Once connected, one can run ligprep from /opt/schrodinger2013/ligprep. For a complete manual on ligprep, refer to http://isp.ncifcrf.gov/files/isp/uploads/2010/07/lp23_user_manual.pdf. Example of simple usage is:

ligprep -g -epik -n 1:23 -isd inputfile.sdf -omae outputfile.mae.

-g keeps current chiral state( important for ZINC, as not all chiralities are avaialable). -epik uses epik for tautomerization. -n 1:23 tells the server to only run ligprep on molecule 1-23. If one wants to run multiple instances of ligprep, one can run using same input and use -n to divide molecules between processes. -isd is to indicate input file is in .sdf format, and -omae indicates .mae output format. The output file(s) are now ready for docking.

Docking using Glide

Do not attempt to dock using maestro's interface. It cannot submit jobs to the nodes. instead, start glide using the command line interface. Before starting a job, one needs to make a .in file. an example .in file is :


GRIDFILE "/users/mh43854/Desktop/glide_prepped/glide-grid_1.zip"

LIGANDFILE "/users/mh43854/Desktop/glide_prepped/ligprep_4-out.maegz"


I don't know the function of USECOMPMAE. GRIDFILE should point to the prepared protein grid file( a .zip file). LIGANDFILE points to the output of ligprep, and PRECISION selects the degree of precision , either HTVS, SP, or XP, in increasing order of precision.To run glide, connect to an available node, and run the command

  /opt/software/Schrodinger_Suites_2014/glide MY_INFILE.in.

All options are written into the in file, so no more options are neccesary. By default, log and output files are written to the current working directory. output is as one large .mae file readable by maestro, and contains all poses as well as score information.

license server (root required)


 /opt/goldsuite-5.2/bin/gold_licence status
/opt/software/Schrodinger_Suites_2014/licadmin stat

to check what is running.

In bashrc on NOVA:

export SCHRODINGER="/opt/software/Schrodinger_Suites_2014" 
export tmpdir=/scratch

to run of ti3d:

       ssh into one of the compute nodes ( names accessible by ganlia cpu-user)
      export SCHRODINGER="/apps/schrodinger/2010_update_1/"


/opt/software/Schrodinger_Suites_2014/licadmin start -c ./license

You may have to stop GOLD server first:

/opt/goldsuite-5.2/bin/gold_licence stop (this now works on uranus)

The Schrodinger license_key file was generated for bme-nova and use the MAC of em1 ($SHRODINGER/machid)

Client (node50-56)


See  /opt/software/Schrodinger_Suites_2014/schrodinger.hosts

To run Schrodinger on selected nodes requires a HOST file. A sample is shown below:

name:        localhost
schrodinger: /opt/software/Schrodinger_Suites_2014
tmpdir:      /scratch
name:        node50
host:        nova
schrodinger: /opt/software/Schrodinger_Suites_2014
tmpdir:      /scratch
processors:  8
name:        node51
host:        nova
schrodinger: /opt/software/Schrodinger_Suites_2014
tmpdir:      /scratch
processors:  8

I (David Bell) experienced trouble logging into Schrodinger as the root user. It seems like the chosen nodes had trouble communicating when logged in as root. Finally, Schrodinger employs a gui that requires the X display server (not BME-mars).




howto (2019//6/7)

On node 38, follow instructions at https://admiring-tesla-08529a.netlify.com/installs/v132/( installer)


Psi4 is available at /opt/software/CHEM/PSI4/PSI4new.

Most nodes use centos 6, source /opt/software/CHEM/PSI4/.bashrc.psi4.centos6

For centos7(node36) source /opt/software/CHEM/PSI4/.bashrc.psi4.centos7

The only difference is in Glibc usage, centos 6 requires an external newer version of glibc, while centos7 has a newer version pre-installed .

how to calculate energy component of a dimer.

To run PSI4:

a). First, you need to define some parameters and settings in a '.inp' file for PSI4 to recognize.:

memory 10 gb

molecule Pyridine_Ethene {
0 1
N    1.381382191  -0.000233477   0.131463739
C    0.679350788  -1.140239457   0.092079660
H    1.258719604  -2.054962232   0.125883611
C   -0.709722319  -1.193114069   0.006664265
H   -1.214087679  -2.148561633  -0.025308510
C   -1.421613574   0.000133427  -0.040816901
H   -2.500696146   0.000257571  -0.109169733
C   -0.709401197   1.193175380   0.006521983
H   -1.213511626   2.148747839  -0.025528309
C    0.679651666   1.139956229   0.091893028
H    1.259260730   2.054510898   0.125502478
0 1
C    0.016268721   0.666423181   3.116599871
H    0.926742724   1.225908900   2.957485028
H   -0.893278777   1.228827408   3.273560375
C    0.016601405  -0.666264119   3.116732110
H    0.927337096  -1.225346603   2.957721677
H   -0.892686689  -1.229088890   3.273821261
units angstrom
no_reorient  #important for SAPT in psi4, default
symmetry c1  #important for SAPT in psi4, default

set {
basis aug-cc-pVDZ
scf_type DF
freeze_core True

b). Second, you need to call PSI4:


export PSI_SCRATCH=/scratch/YourUserName

if [ ! -d "$PSI_SCRATCH" ]

source /opt/intel/compilerpro- intel64


$my_psi4/bin/psi4 -n 8 -i 2901_26UracilUracilpipi100_TZ.inp -o 2901_26UracilUracilpipi100_TZ.out

By defining 'PSI_SCRATCH', user can tell PSI4 to store the temporary data into this directory. Without it, temporary data willl be stored in the default directory of PSI4, which may not be big enough.

Molpro 2010

Installation Info

a) Making a ~/.molpro dir, and copying the token file into it

b) Configured with the following options:

  ./configure -icc -ifort -mpp -mppbase /opt/mpich2-fc13-p262/include/ 
      -blaspath /opt/intel/compilerpro- 
      -lapackpath /opt/intel/compilerpro- 
      -prefix /home/qwang/local

c) For some unknown reason, the generated CONFIG file contains a invalid path to mpich2, so we need to manually correct it to:

  -I/opt/mpich2-fc13-p262/include -L/opt/mpich2-fc13-p262/lib

d) The final installation step will install the program in a new dir under -prefix, and cp the token file under -prefix/lib automatically

e) The full guide of Molpro can be found at: <a href="http://www.molpro.net/info/current/doc/manual/node4.html">http://www.molpro.net/info/current/doc/manual/node4.html</a>

f) Example of input files, including SAPT calculations can be found at: /home/qwang/software/Molpro/examples

g) To run Molpro, please make sure you source the lib of MKL and Intel compilers (/opt/intel/compilerpro-

An example script to run Molpro would be:


export TMPDIR=/scratch/dir_name_you_prefer
export LD_LIBRARY_PATH=/opt/mpich2-fc13-p262/lib:$LD_LIBRARY_PATH

source /opt/intel/compilerpro- intel64


if [ ! -d "$TMPDIR" ]
mkdir $TMPDIR

$mpich2_home/bin/mpiexec -np 1 $molpro_home/bin/molpro -n 4 -s h2odimer_sapt_hf.com


An example script (if you have trouble, remove your bashrc or cshrc)

export TMPDIR=/scratch/pren
export LD_LIBRARY_PATH=/opt/mpich2-fc13-p262/lib:$LD_LIBRARY_PATH
# this doesn't seem to be necessary
#source /opt/intel/compilerpro- intel64
#$mpich2_home/bin/mpdboot -f mpd.hosts -n 1 --ncpus=8 &
$mpich2_home/bin/mpd  & 
if [ ! -d "$TMPDIR" ]
mkdir $TMPDIR
#$mpich2_home/bin/mpiexec -np 1 $molpro_home/bin/molpro -n 4 -s h2odimer_sapt_hf.com
#$molpro_home/bin/molpro -n 4 -s h2odimer_sapt_hf.com
#$mpich2_home/bin/mpiexec -np 1 $molpro_home/bin/molpro -n 4 -s ch2cl2_pbe0_opt_ion_pot.com >&/dev/null &
$molpro_home/bin/molpro -n 4 -s ion_pot.com > log &

ion_pot.com is the input file has the structures of the molecules and options for calculations.

!about 1.5G mem
1, O,,   -0.525329794  -0.050971084  -0.314516861
2, H,,   -0.942006633   0.747901631   0.011252816
3, H,,    0.403696525   0.059785981  -0.073568368
!ground state
!one electron less


Scripts to compute T-wham free energy:

Matlab script to draw the phi-psi map for dialanine

<a href="http://biomol.bme.utexas.edu/~yues/remd/">http://biomol.bme.utexas.edu/~yues/remd/</a>

<a href="http://biomol.bme.utexas.edu/~yues/remd/energy">http://biomol.bme.utexas.edu/~yues/remd/energy</a> has the energies extracted at different temperatures. <a href="http://biomol.bme.utexas.edu/~yues/remd/dihedral">http://biomol.bme.utexas.edu/~yues/remd/dihedral</a> has the phi psi angles for all snapshots. <a href="http://biomol.bme.utexas.edu/~yues/remd/map">http://biomol.bme.utexas.edu/~yues/remd/map</a> has Jenny's matlab codes (.m files) with T-wham algorithm. All the other files in <a href="http://biomol.bme.utexas.edu/~yues/remd/">http://biomol.bme.utexas.edu/~yues/remd/</a> are trajectories from AMBER.



A compiled version of SAPT2008 can be found at: /home/qwang/software/SAPT2008

Installation Info

a) Compiled with gcc4.1.2 compilers installed at /home/qwang/software/gcc-4.1.2-install

b) Higher version of gcc installed on the nodes does NOT work for ATMOL1024.

c) Examples to run SAPT with ATMOL can be found at /home/qwang/software/SAPT


d) Useful links:

  <a href="http://www.physics.udel.edu/~szalewic/SAPT/manual.html">http://www.physics.udel.edu/~szalewic/SAPT/manual.html</a>
<a href="http://www.theochem.ru.nl/~pwormer/atmol/integw.html">http://www.theochem.ru.nl/~pwormer/atmol/integw.html</a>

Intel Compiler

V13 fortran and c compilers are installed (1/2014). The installation directory for both ifort and icc is "/opt/intel/composer_xe_2013.2.146". The bin in this directory is also linked to composer_xe_2013/

mkl library is included within.

Csh or tcsh:

source /opt/intel/composer_xe_2013.2.146/bin/compilervars.csh intel64



source /opt/intel/composer_xe_2013.2.146/bin/compilervars.sh intel64

Open MM


needs to be done on a node with GPU and Cuda( ex. node 201)

may need to install ccmake(http://www.cmake.org/), depending on version present on machine (or do yum install cmake as root; may also need swig, flex etc)

 download openMM  source from https://simtk.org/project/xml/downloads.xml?group_id=161

make a build folder

source bash script from /opt/software/bash/bashrc.openmm

run ccmake ( ccmake -i (source directory)) from inside build directory.

   press "c" on keyboard

turn on desired features ( such as CUDA, etc.), both cuda directories should point to /usr/local/cuda-6.5

   press g


make install

make PythonInstall

perform tests inside build directory ( make test)

cd into ~pren/openMM/dhfr-test/ and  run simulateMTS.py from within this directory. Expect performance of around 3 ns/day ( observed in node 201).



Softcore implementation


in ccmake , need to set C compiler flag --std=c++0x. This is needed to enable use of tuples

changed /plugins/amoeba/openmmapi/include/openMM/AomebaVdwForce.h to have lambda particle descriptor

changed /plugins/amoeba/openMMapi/src/AmoebaVdwForceimpl.cpp to use soft core Vdw.

changed /plugins/amoeba/plaforms/cuda/src/AmoebaCudaKernals.cpp so that SigmaEpsilon array contains a 3rd variable ( lambda)

Input Set up

  • 1) you need pdb file and xml file that contains parameters. Examples are in ~pren/openMM/dhfr-test/. OpenMM match the amoeba2013.xml with the pdb file and assign parameters by using residue names.
  • 2) or openMM can use a "system" xml file which contains parameters for all particles. Example: /home/pren/OpenMM/Output-2fs/7.jac-dhfr/mutual-1e-5/ see input.xml and input.pdb

A script for doing this (not sure exactly which one of above):

 python processTinkerForceField.py amoebapro13.prm

Need to modify the residueFinal.xml location See /home/pren/OpenMM/Accuracy/7.jac-dhfr/mutual-1e-6

---2014-- Xiao, I am able to run OpenMM/AMOEBA on the GPU (BME-SUGAR) now. You will need my ~/.bashrc file to set up all the env vars and paths.

The openMM source (precompiled binary) is in /opt/software/OpenMM6.1-Linux64/. It is installed into /usr/local/openmm, and /opt/software/python27/appdata/canopy-

The two examples are /home/pren/OpenMM/zzz/simulate.py and simulateMTS.py. The later uses larger time step.

A couple of things we need to work with Lee-Ping

1) the tinker/OpenMM interface that we need to convert tinker xyz/key to OpenMM input (xml). The above examples uses a PDB and amoeba2013.xml and OpenMM will assign the parameters. However for complicated system this would not work.

2)The soft core stuff. TINKER uses “ligand -100 110” to specific a solute (atom index between 100 and 110); VDW-LAMBDA and ELE–LAMBDA to scale the vdw and ele interaction between ligand and the rest. The vdw soft core is implement in ehal.f and ehal3.f (energy) and ehal1.f (force), search vlambda. The electrostatic part just scale the multipoles and polarizability by the “elambda”. This means we turn off ele intramolecular interaction when decouple it so recharging in gas-phase is necessary (but cheap).

Peter suggests: The easiest way for you to do it is probably to start from the existing AmoebaVdwForce class and modify it to use a soft core potential. The relevant code is in amoebaVdwForce2.cu for the CUDA platform, and AmoebaReferenceVdwForce.cpp for the Reference platform. It should be quite straightforward to change.

---2014 from Lee-Ping-- (1) You can use Peter's "pdbfixer" code (https://github.com/SimTk/pdbfixer) to set up any PDB you want (you start with the PDB ID and end up with a simulation-ready PDB file), and simulateAmoeba.py should run it properly.

However, if you keep molecules and residues not supported in the force field XML file it will crash. There is a script to convert TINKER .key/.prm files into OpenMM .xml files (https://www.dropbox.com/s/7urbdgl0bja6hmd/processTinkerForceField.py?dl=0); I believe Mark Friedrichs wrote the code when he implemented AMOEBA into OpenMM.

(2) To use the multiple timestep Langevin integrator, paste the code for MTSVVVRIntegrator() into simulateAmoeba.py.

Then instead of LangevinIntegrator(298.15*U.kelvin, 1.0/U.picosecond, 1.0*U.femtosecond) use MTSVVVRIntegrator(298.15*U.kelvin, 1.0/U.picosecond, 2.0*U.femtosecond, system).

The extra "system" argument is needed because the MTSVVVRIntegrator partitions the forces in the system into different groups (slow and vast), and applies a different time step to the force groups.


Generate structures of orientations

In AMOEBA parameter determining process, we often change the VDW parameters to match the total interaction energy. In this situation, QM data from different orientated dimers are very necessary. Here I will provide a python script to automatically (or semi-automatically) generate dimers with different orientation.

The basic idea of how we get a rotated (x',y',z') coordinate is to first construct a rotation matrix, in which the rotation axis and rotation angle should be specified. Then we simply apply the dot product of the rotation matrix and the original coordinate (x,y,z).

Here is the script used for ChloroBenzene-Water. It could be easily modified to work on your case if you understand how this works.

Python script

 #!/usr/bin/env python

 import numpy as np
 import math
 import sys

 '''Rotate the first Molecule along an Axis for Theta radians'''
 '''Chengwen Liu'''

 def rotMatrix(axis, theta):
   Return the rotation matrix associated with counterclockwise rotation about
   the given axis by theta radians.
   axis = np.asarray(axis)
   axis = axis/math.sqrt(np.dot(axis, axis))
   a = math.cos(theta/2.0)
   b, c, d = -axis*math.sin(theta/2.0)
   aa, bb, cc, dd = a*a, b*b, c*c, d*d
   bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d
   return np.array([[aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],
                    [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],
                    [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]])

 def rotate_water(txyz):
   """rotate water molecule; fix benzenechloride """
   lines = file(txyz).readlines()
   coord = []
   for line in lines[1:]:
     data = line.split()
     coord.append( [float(data[2]), float(data[3]),float(data[4]),])

   #read some points
   Cl_p = np.array(coord[9])
   C6_p = np.array(coord[5])
   C7_p = np.array(coord[6])
   C8_p = np.array(coord[7])

   Wox_p1 = np.array(coord[0])
   Wh1_p1 = np.array(coord[1])
   Wh2_p1 = np.array(coord[2])

   #define some arrays
   C7C6_v = []
   C7C8_v = []
   C6C8_v = []

   for i in xrange(3):
     C7C6_v.append(C6_p[i] - C7_p[i])
     C7C8_v.append(C8_p[i] - C7_p[i])
     C6C8_v.append(C8_p[i] - C6_p[i])


   theta_list = []
   for i in xrange(0,7,1):
   axis1 = C6C8_v

   for theta in theta_list:
     Wox_p2 = np.dot(rotMatrix(axis1,theta), Wox_p1)
     Wh1_p2 = np.dot(rotMatrix(axis1,theta), Wh1_p1)
     Wh2_p2 = np.dot(rotMatrix(axis1,theta), Wh2_p1)

     '''rotate around the axis across Cl atom, so we need a translation'''
     trans_v = (Wox_p2 - Cl_p)/np.linalg.norm(Wox_p2-Cl_p) * (np.linalg.norm(Wox_p2-Cl_p) - np.linalg.norm(Wox_p1 - Cl_p))
     Wox_p2 = Wox_p2 - trans_v
     Wh1_p2 = Wh1_p2 - trans_v
     Wh2_p2 = Wh2_p2 - trans_v

     filename = str("%.1f"%float(theta*180.0/np.pi))
     ifile = open("chlorobenzH20_OP_%sdeg.txyz"%filename,'w') #axis1 out of plane rotation
     ifile.write("%s%12.6f %12.6f %12.6f%s\n"%(lines[1][:6],Wox_p2[0], Wox_p2[1], Wox_p2[2], lines[1][49:-1]))
     ifile.write("%s%12.6f %12.6f %12.6f%s\n"%(lines[2][:6],Wh1_p2[0], Wh1_p2[1], Wh1_p2[2], lines[2][49:-1]))
     ifile.write("%s%12.6f %12.6f %12.6f%s\n"%(lines[3][:6],Wh2_p2[0], Wh2_p2[1], Wh2_p2[2], lines[3][49:-1]))
     for i in xrange(4,len(lines),1):


I rotated along the axis parallel to C6-C8 and across Cl atom (see txyz file)

Initial Tinker file

  1 O    4.553365    -0.050957     0.011596        1     2     3
  2 H    5.042109     0.300726    -0.736075        2     1
  3 H    5.043243     0.270403     0.772051        2     1
  4 C   -3.077858     0.019185     0.016973      402     5    10    11
  5 C   -2.362505     1.222155    -0.004935      403     4     6    12
  6 C   -0.963557     1.214342     0.010611      401     5     7    13
  7 C   -0.282885    -0.007338    -0.024152      404     6     8     9
  8 C   -0.986495    -1.215973     0.010115      405     7
  9 C   -2.385288    -1.197092    -0.005379      401     7    10    14
 10 Cl   1.451131    -0.025174    -0.010370      403     4     9    15
 11 H   -4.164013     0.029438     0.010388      407     4
 12 H   -2.890909     2.171654     0.010734      406     5
 13 H   -0.401633     2.143467    -0.001102      408     6
 14 H   -0.441889    -2.155321    -0.001783      408     9
 15 H   -2.931533    -2.136433     0.009921      406    10

Here are the rotated structures.