separate ligand bond rotations (bug fixed), add ligand parameterizati…

…on scripts

par-scripts/antechamber_script

@@ -0,0 +1,28 @@
+#!/bin/bash -v
+
+#echo "Starting openbabel to add H atoms to the ligand" 
+# Use Babel to add H atoms [ this is already done on converting SDF files to PDB]
+#babel -h -ipdb $ligand.pdb -opdb $ligand.babel-h.pdb >> $molname.err 
+#obabel -h -ipdb $ligand.pdb -opdb $ligand.babel-h.pdb >> $molname.err
+
+#if [ ! -s $ligand.babel-h.pdb ]
+#then
+#        echo "FAILURE: Unable to add H atoms to the ligand. Some likely reasons are invalid pdb format or very non-standard atom name. Check to make sure atom name starts in column 14 of the pdb file. Check the $molname.err file for clues."
+#     exit 1
+#else
+
+# Prepare/read in ligand parameters
+
+# Set initial ligand net charge to zero or whatever is appropriate [filled in by parameterize-it]
+   #ligand_charge=0
+
+   echo "Starting antechamber to create prep file for tleap. Assuming ligand net charge = "  $ligand_charge >> $molname.err
+# Use Amber Antechamber to create tleap prepi file
+   $ACHOME/bin/antechamber -nc $ligand_charge -s 2 -c bcc -i $ligand.pdb -fi pdb -o pdb2pqr.ligand.prepi -fo prepi -at gaff >> $molname.err 
+   if [ ! -s pdb2pqr.ligand.prepi ]
+   then
+      echo "WARNING: Likely inconsistency between the assumed net charge and the number of electrons."
+   fi
+#fi
+
+

par-scripts/check-all

@@ -0,0 +1,19 @@
+#!/bin/bash
+
+rm results
+for drug in `cat successful-drugs`
+#for drug in `echo "oht"`
+do
+	echo $drug
+	cd $drug
+	sander -i ../amber.in -o amber.out -p $drug.prm7 -c $drug.inpcrd -O
+	amber=`awk '/DIHED/ {bond=$3; angle=$6; dih=$9} /VDWAALS/ {vdw=$3; elec=$6} /1-4 VDW/ {vdw14=$4; elec14=$8} END {print bond, angle, dih, vdw+vdw14, elec+elec14}' amber.out`
+	#echo $amber 
+	cp NEWPDB.PDB $drug-test.pdb
+	sed -e "s/:drug;/$drug/g" ../new-exact-energy.txt > exact-energy-$drug.txt
+	~/justin-mrmc/mrmc-7-28-16/mrmc-debug  exact-energy-$drug.txt > exact-energy-$drug.out
+	mrmc=`grep -i "Energy:" exact-energy-$drug.out | awk '($2==0) {print $4,$5,$6,$7,$8,$9,$10}'`
+	cd ..
+	comp=`echo "$drug $amber $mrmc" | awk '{print $1,$7-$2,$8-$3,$9+$10-$4,$11-$5,$12-$6}'`
+	echo $comp >> results
+done

par-scripts/check-collect

@@ -0,0 +1,20 @@
+#!/bin/bash
+
+rm successful-drugs
+rm failed-drugs
+for drug in `cat drug-codes`
+do
+	echo $drug
+	if [ -s $drug/$drug.prm ]; then
+		cp -v $drug/amber99sb-$drug.prm params/
+		cp -v $drug/defs-amber99sb-$drug.txt params/
+		cp -v $drug/$drug.pdb params/
+		cp -v $drug/matches-$drug params/
+		cp -v $drug/$drug.frcmod params/
+		cp -v $drug/$drug.prepi params/
+		#cp -v $drug/$drug.prm7 params/
+		echo $drug >> successful-drugs
+	else
+		echo $drug >> failed-drugs
+	fi
+done

par-scripts/check-parameters

@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rname=$1 #"rnamel" from previous script
+

par-scripts/convert-to-tinker.py

@@ -0,0 +1,163 @@
+#!/usr/bin/env python
+
+import sys
+#to install parmed: python setup.py install --home=~
+#To use parmed: echo $PYTHONPATH=~/lib/python:$PYTHONPATH
+from parmed.amber import *
+#import parmed
+parm = LoadParm(sys.argv[1])
+typestart = int(sys.argv[2])
+classstart = int(sys.argv[3])
+tinkerfname = sys.argv[4]
+tinkerfile = open(tinkerfname,'w')
+
+
+natom = parm.ptr('natom')
+ntypes = parm.ptr('ntypes')
+#need to get all the bonds, angles, and dihedrals
+nbond = parm.ptr('nbonh') + parm.ptr('mbona')
+nangle = parm.ptr('ntheth') + parm.ptr('mtheta')
+ndihedral = parm.ptr('nphih') + parm.ptr('mphia')
+#tinkerfile.write(natom, ntypes, nbond
+
+nextclass = classstart
+atomlinefmt = 'atom {0:d} {1:d} {2} "{3}"  {4:d} {5:.3f} {6:d}\n'
+classname2num = {}
+for iatom in range(0,natom):
+	#tinkerfile.write(parm.atoms[iatom].name, parm.atoms[iatom].type, parm.atoms[iatom].mass, parm.atoms[iatom].charge, parm.atoms[iatom].nb_idx
+	type = iatom + typestart
+	classname = parm.atoms[iatom].type
+	if (classname not in classname2num):
+		iclass = nextclass
+		classname2num[classname] = nextclass
+		#tinkerfile.write('# {0:d} {1}'.format(iclass,classname)
+		nextclass = nextclass + 1
+	else:
+		iclass = classname2num[classname]
+	#todo: the "0" should be the atomic number
+	tinkerfile.write(atomlinefmt.format(type, iclass, parm.atoms[iatom].type, parm.atoms[iatom].name, parm.atoms[iatom].atomic_number, parm.atoms[iatom].mass, 1))
+
+for classname in classname2num:
+	iclass = classname2num[classname]
+	tinkerfile.write('# {0} {1:d}\n'.format(classname,iclass))
+
+bondlinefmt = 'bond {0:d} {1:d} {2:.2f} {3:.4f} #{4} {5}\n'
+bondset = set() #to make sure we don't duplicate
+for ibond in range(0,nbond):
+	iatom = parm.bonds[ibond].atom1
+	jatom = parm.bonds[ibond].atom2
+	iclass = classname2num[iatom.type]
+	jclass = classname2num[jatom.type]
+	k = parm.bonds[ibond].type.k
+	req = parm.bonds[ibond].type.req
+	#todo:really should sort these
+	if ((iclass,jclass) not in bondset):
+		tinkerfile.write(bondlinefmt.format(iclass,jclass,k,req,iatom.name,jatom.name))
+		bondset.add((iclass,jclass))
+
+anglelinefmt = 'angle {0:d} {1:d} {2:d} {3:.2f} {4:.2f} #{5} {6} {7}\n'
+angleset = set()	
+for iangle in range(0,nangle):
+        iatom = parm.angles[iangle].atom1
+        jatom = parm.angles[iangle].atom2
+	katom = parm.angles[iangle].atom3
+        iclass = classname2num[iatom.type]
+        jclass = classname2num[jatom.type]
+	kclass = classname2num[katom.type]
+        k = parm.angles[iangle].type.k
+        theteq = parm.angles[iangle].type.theteq
+        #todo:really should sort these
+        if ((iclass,jclass,kclass) not in angleset):
+                tinkerfile.write(anglelinefmt.format(iclass,jclass,kclass,k,theteq,iatom.name,jatom.name,katom.name))
+        	angleset.add((iclass,jclass,kclass))
+
+#For the dihedrals, we separate identifying the parameters and writing them out, so that we can come up with multiple dihedrals as needed
+ndih = {}
+impdih = {}
+dihset = set()
+for idih in range(0,ndihedral):
+	iatom = parm.dihedrals[idih].atom1
+        jatom = parm.dihedrals[idih].atom2        
+        katom = parm.dihedrals[idih].atom3    
+	latom = parm.dihedrals[idih].atom4
+	iclass = classname2num[iatom.type]      
+        jclass = classname2num[jatom.type]
+        kclass = classname2num[katom.type]
+	lclass = classname2num[latom.type]
+	dihset.add((iclass,jclass,kclass,lclass))
+	per = parm.dihedrals[idih].type.per 
+	classes = (iclass, jclass, kclass, lclass, per)
+	phi_k = parm.dihedrals[idih].type.phi_k
+	per = int(parm.dihedrals[idih].type.per)
+	phase = parm.dihedrals[idih].type.phase
+	dihparams = (phi_k, phase)
+	if (parm.dihedrals[idih].improper):
+		if (classes not in impdih):
+			impdih[classes] = dihparams
+	else: #normal dihedral
+		if (classes not in ndih):
+			ndih[classes] = dihparams
+
+#tinkerfile.write(ndih#
+#tinkerfile.write(impdih
+
+imptorsfmt='imptors {0:d} {1:d} {2:d} {3:d} {4:.3f} {5:.1f} {6:d}\n'
+#assume there is only one improper dihedral for each combination of classes.
+for classes in impdih:
+	iclass,jclass,kclass,lclass,per = classes
+	phi_k,per = impdih[classes]
+	per = int(per)
+	tinkerfile.write(imptorsfmt.format(iclass, jclass, kclass, lclass, phi_k, phase, per))
+
+
+#these do not have \n so we can continue
+torsfmt1='torsion {0:d} {1:d} {2:d} {3:d} '
+torsfmt2=' {0:.3f} {1:.1f} {2:d}' #magnitude, phase, multiplicity
+for classes in dihset:
+	iclass, jclass, kclass, lclass = classes
+	#tinkerfile.write(torsfmt1.format(iclass,jclass,kclass,lclass),
+	firsttorsion = True
+	for per in range(1,7):
+		ndihlookup = (iclass, jclass, kclass, lclass, per)
+		if (ndihlookup in ndih):
+			phi_k, phase = ndih[ndihlookup]
+			if (firsttorsion):
+				tinkerfile.write(torsfmt1.format(iclass,jclass,kclass,lclass))
+				firsttorsion = False
+			tinkerfile.write(torsfmt2.format(phi_k, phase, per))
+	if (not firsttorsion):
+		tinkerfile.write('\n')
+
+
+
+vdwfmt="vdw {0:d} {1:.4f} {2:.4f}\n"
+#tinker's files have charges to only 4 sig digits
+chargefmt="charge {0:d} {1:.4f}\n"
+
+for classname in classname2num:
+	classnum = classname2num[classname]
+	LJtype = parm.LJ_types[classname]-1
+	tinkerfile.write(vdwfmt.format(classnum, parm.LJ_radius[LJtype], parm.LJ_depth[LJtype]))
+
+#Force an integer charge.
+
+totalcharge = 0.0
+for itype in range(0,natom):
+        type = itype + typestart
+        totalcharge += parm.atoms[itype].charge
+
+adjust = -(totalcharge - round(totalcharge,0))/natom
+print "Adjusting charges by {0:.6f} to get total charge of {1:.1f}".format(adjust,round(totalcharge,0))
+
+totalcharge2 = 0.0
+for itype in range(0,natom):
+	type = itype + typestart
+        charge =  parm.atoms[itype].charge+adjust
+        totalcharge2 += charge
+        print parm.atoms[itype].name,charge
+	tinkerfile.write(chargefmt.format(type, charge))
+
+print totalcharge2
+tinkerfile.close()
+
+