Skip to content

Instantly share code, notes, and snippets.

@platinhom

platinhom/0-PQR-related-process.md

Last active March 15, 2016 20:33
Show Gist options
  • Save platinhom/bcef539bd97bc192e231 to your computer and use it in GitHub Desktop.
Save platinhom/bcef539bd97bc192e231 to your computer and use it in GitHub Desktop.
Download ZIP
To deal with pqr/pqra/pqrta
Raw
combinePQRTA.py
#! /usr/bin/env python
# -*- coding: utf8 -*-
# Author: Platinhom; Last Updated: 2015-11-26
#
# Argv1: A pqr/pqra/pqrt file name, used to parse file prefix.
# Need pqrt and pqra file with same prefix in same place.
import os,sys
fname=sys.argv[1]
fnamelist=os.path.splitext(fname)
frpqra=open(fnamelist[0]+'.pqra','r')
frpqrt=open(fnamelist[0]+'.pqrt','r')
fw=open(fnamelist[0]+'.pqrta','w')
atomarea=[]
for line in frpqra:
if (line[:6]=='ATOM ' or line[:6]=='HETATM'):
aarea=line[78:90]
atomarea.append(aarea)
elif (line[:6]=='REMARK'): fw.write(line)
nowi=0
for line in frpqrt:
if (line[:6]=='ATOM ' or line[:6]=='HETATM'):
fw.write(line.strip()+atomarea[nowi]+'\n')
nowi+=1;
fw.close();frpqra.close();frpqrt.close();
Raw
ExtractXYZR.py
#! /usr/bin/env python
# -*- coding: utf8 -*-
# Author: Platinhom; Last Updated: 2015-09-30
# Usage: Extract XYZR from PQR file and save it in xyzr format.
# Only for pqr format now.
# python ExtractXYZR.py file.pqr
import os,sys
supportFormat=[".pqr",]
if (__name__ == '__main__'):
if (len(sys.argv)!=2):
print "Give the input file!"
exit()
fname=sys.argv[1]
fnamelist=os.path.splitext(fname)
if (fnamelist[1]==".xyzr"):
print "Input file is xyzr format! Exit!"
exit()
if (not fnamelist[1].lower() in supportFormat):
print fnamelist[1][1:]+" format can't be supported now.."
exit()
fxyzr=open(fnamelist[0]+".xyzr",'w')
fr=open(fname)
inlines=fr.readlines();
fr.close();
# Write out the corresponding xyzr file.
for line in inlines:
# Each atom
if (line[:4]=="ATOM" or line[:6]=="HETATM"):
# Extract x, y, z, r from pqr to xyzr file
radius="%10.5f" % float(line[62:70].strip());
xcoor="%10.5f" % float(line[30:38].strip());
ycoor="%10.5f" % float(line[38:46].strip());
zcoor="%10.5f" % float(line[46:54].strip());
xyzrstr=xcoor+ycoor+zcoor+radius+"\n";
fxyzr.write(xyzrstr);
fxyzr.close()
#end main
Raw
HIS-HIE-HID-HIP-rename.py
#! /usr/bin/env python
# -*- coding: utf8 -*-
""" To rename the residues based on protonization"""
import os,sys
def getatomname(line):
return line[12:16]
def getresname(line):
return line[17:20]
def getresnum(line):
return line[22:26]
if (__name__ == "__main__"):
if (len(sys.argv)!=2):
print "Please give input file(pdb/pqr)."
input()
exit()
fname=sys.argv[1]
fnamelist=os.path.splitext(fname)
fwname=fnamelist[0]+"_HM"+fnamelist[1]
fr=open(fname)
fw=open(fwname,"w")
lines=fr.readlines()
startH=False
endH=0
for i in range(len(lines)):
line=lines[i]
if (line[:4]=="ATOM" or line[:6]=="HETATM"):
# His line have been process
if (startH and i<endH):
continue;
elif (startH and i==endH):
startH=False;
restype=getresname(line);
# first line of HIS residue
if (not startH and (restype=="HIS" or restype=="HID"
or restype=="HIE" or restype=="HIP")):
startH=True
nowrid=getresnum(line)
atomnames=[getatomname(line).strip()]
newlines=[line]
for j in range(i+1,i+30):
line2=lines[j]
rid2=getresnum(line2)
if (rid2!=nowrid):
# set up end pos
endH=j
break
if (not (line2[:4]=="ATOM" or line2[:6]=="HETATM")):
# meet the TER end
endH=j
break
newlines.append(line2)
atomnames.append(getatomname(line2).strip())
# Find the right residue
shouldname="HIS"
if ("HD1" in atomnames and "HE2" in atomnames):
shouldname="HIP"
elif ("HD1" in atomnames):
shouldname="HID"
elif ("HE2" in atomnames):
shouldname="HIE"
for line3 in newlines:
fw.write(line3[:17]+shouldname+line3[20:])
continue
else:
startH=False
fw.write(line)
fr.close()
fw.close()
Raw
mol2pqr.sh
#! /bin/bash
# File: mol2pqr.sh
# Author: Zhixiong Zhao, Last Updated: 2015.10.1
#
# Change the input molecule to pqr file with vdw radius.
# Input 1: molecule file in pdb/mol2/mpdb/pqr/ac format
# Input 2: output file format in ac/pdb/mol2/mpdb/xyzr/pqr/pqrt/pqra/pqrta format
# Input 3: charge type: no(default), bcc, mul, gas
# Input 4: atom type: gaff(default), amber, sybyl, bcc, gas
# Input 5: vdw radius set: mbondi(default), mbondi2, mbondi3, bondi, amber6, zap9
# Input 6: set any value to keep pqr file when for conversion to pqrt/pqra/pqrta
#
# Usage: ./mol2pqr.sh input.mol2 pqr bcc gaff mbondi
# Need: python, ambertools, MS_Intersection(for pqra/pqrta)
# Custom: $AMBERHOME enviroment, MS_Intersection parameter,
# Setup the amber environment. You should modify by your own environment
if [ -z $AMBERHOME ];then
source /mnt/home/zhaozx/AmberTools/amber.sh
##if you are in HPCC
#module swap GNU Intel/13.0.1.117; module load OpenMPI/1.4.4; module load Amber/14v19;
if [ -z $AMBERHOME ];then
echo "No AmberTools is installed or its path can't be known!"
exit
fi
fi
probe=1.4
grid=0.2
solbuffer=4.0
if [ -z $1 ];then
echo "Please assign the input file!"
exit
fi
outFormat=`echo $2 | tr 'A-Z' 'a-z'`
if [ -z $outFormat ];then
echo "No output format is given! Can be: "
echo " pdb/mol2/mpdb/pqr/ac/pqrt/pqra/pqrta"
echo " Now using pqr format..."
outFormat=pqr
fi
chargetype=`echo $3 | tr 'A-Z' 'a-z'`
if [[ -z $chargetype || ( $chargetype != "bcc" && $chargetype != "mul" \
&& $chargetype != "gas" && $chargetype != "no" ) ]];then
echo "Please assign the charge method! Can be:"
echo " no: Not to change the charges"
echo " bcc: for AM1-BCC charge"
echo " mul: for Mulliken charge"
echo " gas: for Gasteiger charge"
echo " cm1, cm2, esp, resp charge need more programs.. so don't support here."
echo " Now using no charge or origin charge...."
chargetype=no
fi
# atom type as: gaff(default), amber, sybyl, bcc, gas.(In amber14)
atomtype=`echo $4 | tr 'A-Z' 'a-z'`
if [[ -z $atomtype || ( $atomtype != "gaff" && $atomtype != "amber" && \
$atomtype != "sybyl" && $atomtype != "bcc" && $atomtype != "gas" ) ]];then
echo "Please assign the atom type in molecule! Can be:"
echo "gaff(default), amber, sybyl, bcc, gas.(In amber14)"
echo "Using default gaff...."
atomtype=gaff
fi
# vdw Type as: bondi, mbondi(default), mbondi2, mbondi3, amber6 (In amber14) and zap9.
vdwtype=`echo $5 | tr 'A-Z' 'a-z'`
if [[ -z $vdwtype || ( $vdwtype != "bondi" && $vdwtype != "mbondi" && $vdwtype != "mbondi2" \
&& $vdwtype != "mbondi3" && $vdwtype != "amber6" && $vdwtype != "zap9") ]];then
echo "Please assign the vdw radius type of atom! Can be:"
echo " bondi, mbondi(default), mbondi2, mbondi3, amber6 (In amber14) and zap9"
echo " Using default mbondi...."
vdwtype=mbondi
fi
MoreChangevdw="False"
if [[ $vdwtype = "zap9" ]];then
MoreChangevdw=$vdwtype
# use mbondi instead first
vdwtype=mbondi
fi
#If need changevdw, change value to "True". And change the vdwtype value
#You can change it to any other string to deactivate it.
#changevdw="True"
outPQR="False"
if [[ ${outFormat:0:3} = "pqr" || $outFormat = "xyzr" ]];then
outPQR="True"
fi
#mpdb will use mbondi, can't be changed!
if [ $outPQR = "True" ];then
echo "Partial charge: $chargetype , atom type: $atomtype , vdw radiis: $vdwtype ."
else
echo "Partial charge: $chargetype , atom type: $atomtype ."
fi
ESES=""
if [[ ${outFormat} = "pqra" || ${outFormat} = "pqrta" ]];then
if [ -f ./MS_Intersection ];then
ESES="./MS_Intersection"
elif which MS_Intersection 2>/dev/null; then
ESES="MS_Intersection"
else
echo "No ESES program (MS_Intersection) exists!"
exit
fi
fi
basename=${1%.*}
exdname=${1##*.}
if [ $chargetype = "no" ];then # No need to change charges
if [ $outPQR = "True" ];then #pqr/pqrt/pqra/pqrta
if [[ $exdname = "mol2" || $exdname = "pdb" || $exdname = "ac" ]];then
antechamber -fi $exdname -fo ac -pf y -i $1 -at $atomtype -o ${basename}_gaff.ac
elif [[ $exdname = "mpdb" || $exdname = "pqr" ]];then
antechamber -fi mpdb -fo ac -pf y -i $1 -at $atomtype -o ${basename}_gaff.ac
else
echo "Only support for pdb/mol2/mpdb/pqr/ac files!"
exit
fi
else #pdb/mol2/mpdb/ac
if [[ $exdname = "mol2" || $exdname = "pdb" || $exdname = "ac" ]];then
antechamber -fi $exdname -fo $outFormat -pf y -i $1 -at $atomtype -o ${basename}.$outFormat
elif [[ $exdname = "mpdb" || $exdname = "pqr" ]];then
antechamber -fi mpdb -fo $outFormat -pf y -i $1 -at $atomtype -o ${basename}.$outFormat
else
echo "Only support for pdb/mol2/mpdb/pqr/ac files!"
exit
fi
fi
else # Need to calculate charges.
if [ $outPQR = "True" ];then #pqr/pqrt/pqra/pqrta
if [[ $exdname = "mol2" || $exdname = "pdb" || $exdname = "ac" ]];then
antechamber -fi $exdname -fo ac -pf y -i $1 -c $chargetype -at $atomtype -o ${basename}_gaff.ac
elif [[ $exdname = "mpdb" || $exdname = "pqr" ]];then
antechamber -fi mpdb -fo ac -pf y -i $1 -c $chargetype -at $atomtype -o ${basename}_gaff.ac
else
echo "Only support for pdb/mol2/mpdb/pqr/ac files!"
exit
fi
else #pdb/mol2/mpdb/ac
if [[ $exdname = "mol2" || $exdname = "pdb" || $exdname = "ac" ]];then
antechamber -fi $exdname -fo $outFormat -pf y -i $1 -c $chargetype -at $atomtype -o ${basename}.$outFormat
elif [[ $exdname = "mpdb" || $exdname = "pqr" ]];then
antechamber -fi mpdb -fo $outFormat -pf y -i $1 -c $chargetype -at $atomtype -o ${basename}.$outFormat
else
echo "Only support for pdb/mol2/mpdb/pqr/ac files!"
exit
fi
fi
fi
# Convert to pqr file
if [ $outPQR = "True" ];then
# Using tleap to generate top file
# Leap only know gaff name and loadmol2, so convert back to mol2 file
antechamber -fi ac -fo mol2 -pf y -i ${basename}_gaff.ac -o ${basename}_gaff.mol2
parmchk -i ${basename}_gaff.mol2 -f mol2 -o ${basename}_gaff.frcmod
echo "source leaprc.gaff">> ${basename}_gaff.leapin
echo "loadamberparams ${basename}_gaff.frcmod" >> ${basename}_gaff.leapin
echo "hom=loadmol2 ${basename}_gaff.mol2">> ${basename}_gaff.leapin
echo "set default PBRadii $vdwtype" >> ${basename}_gaff.leapin
echo "saveamberparm hom ${basename}_gaff.top ${basename}_gaff.crd">> ${basename}_gaff.leapin
echo "quit">> ${basename}_gaff.leapin
tleap -f ${basename}_gaff.leapin >/dev/null
# Change the vdw radius by amber sets by parmed.py
# Old version, now it use leap 'set default PBRadii type' instead.
####
#if [ $changevdw = "True" ];then
#echo "changeRadii $vdwtype">>${basename}_gaff.parmedin
#echo "outparm ${basename}_gaff.top">>${basename}_gaff.parmedin
#parmed.py -p ${basename}_gaff.top -i ${basename}_gaff.parmedin -O > /dev/null
#fi
####
# Use ambpdb to convert amber input to pqr.
ambpdb -p ${basename}_gaff.top -pqr < ${basename}_gaff.crd > ${basename}.pqr
# Change to Other radius...
if [ ! $MoreChangevdw = "False" ];then
if [ $MoreChangevdw = "zap9" ];then
echo "#! /usr/bin/env python" > zap9radiusFromPQR.py
echo "import os,sys" >> zap9radiusFromPQR.py
echo "fname=sys.argv[1]" >> zap9radiusFromPQR.py
echo "fnamelist=os.path.splitext(fname)" >> zap9radiusFromPQR.py
echo "fwname=fnamelist[0]+'_zap9'+fnamelist[1]" >> zap9radiusFromPQR.py
echo "fr=open(fname)" >> zap9radiusFromPQR.py
echo "fw=open(fwname,'w')" >> zap9radiusFromPQR.py
echo "for line in fr:" >> zap9radiusFromPQR.py
echo " if (line[:4]=='ATOM' or line[:6]=='HETATM'):" >> zap9radiusFromPQR.py
echo " element=line.split()[-1].upper();" >> zap9radiusFromPQR.py
echo " radius=' 0.000';" >> zap9radiusFromPQR.py
echo " if (element=='H'): radius=' 1.100';" >> zap9radiusFromPQR.py
echo " elif (element=='C'): radius=' 1.870';" >> zap9radiusFromPQR.py
echo " elif (element=='N'): radius=' 1.400';" >> zap9radiusFromPQR.py
echo " elif (element=='O'): radius=' 1.760';" >> zap9radiusFromPQR.py
echo " elif (element=='F'): radius=' 2.400';" >> zap9radiusFromPQR.py
echo " elif (element=='S'): radius=' 2.150';" >> zap9radiusFromPQR.py
echo " elif (element=='CL'): radius=' 1.820';" >> zap9radiusFromPQR.py
echo " newline=line[:62]+radius+line[70:];" >> zap9radiusFromPQR.py
echo " fw.write(newline);" >> zap9radiusFromPQR.py
echo " else:" >> zap9radiusFromPQR.py
echo " fw.write(line);" >> zap9radiusFromPQR.py
echo "fr.close();" >> zap9radiusFromPQR.py
echo "fw.close();" >> zap9radiusFromPQR.py
python zap9radiusFromPQR.py ${basename}.pqr
rm zap9radiusFromPQR.py ${basename}.pqr
mv ${basename}_zap9.pqr ${basename}.pqr
fi
fi
# Convert to xyzr
if [ $outFormat = "xyzr" ];then
echo "#! /usr/bin/env python" > pqr2xyzr.py
echo "import os,sys" >> pqr2xyzr.py
echo "fname=sys.argv[1]" >> pqr2xyzr.py
echo "fnamelist=os.path.splitext(fname)" >> pqr2xyzr.py
echo "fxyzr=open(fnamelist[0]+'.xyzr','w')" >> pqr2xyzr.py
echo "fr=open(fname)" >> pqr2xyzr.py
echo "for line in fr:" >> pqr2xyzr.py
echo " if (line[:4]=='ATOM' or line[:6]=='HETATM'):" >> pqr2xyzr.py
echo " radius='%10.5f' % float(line[62:70].strip());" >> pqr2xyzr.py
echo " xcoor='%10.5f' % float(line[30:38].strip());" >> pqr2xyzr.py
echo " ycoor='%10.5f' % float(line[38:46].strip());" >> pqr2xyzr.py
echo " zcoor='%10.5f' % float(line[46:54].strip());" >> pqr2xyzr.py
echo " xyzrstr=xcoor+ycoor+zcoor+radius+'\n';" >> pqr2xyzr.py
echo " fxyzr.write(xyzrstr);" >> pqr2xyzr.py
echo "fr.close(); fxyzr.close()" >> pqr2xyzr.py
python pqr2xyzr.py ${basename}.pqr
rm pqr2xyzr.py
fi
# Convert to pqrt
if [ ${outFormat:0:4} = "pqrt" ];then
echo "#! /usr/bin/env python" > pqr_ac2pqrt.py
echo "import os,sys" >> pqr_ac2pqrt.py
echo "fname=sys.argv[1]" >> pqr_ac2pqrt.py
echo "fnamelist=os.path.splitext(fname)" >> pqr_ac2pqrt.py
echo "frac=open(sys.argv[2],'r')" >> pqr_ac2pqrt.py
echo "frpqr=open(fnamelist[0]+'.pqr','r')" >> pqr_ac2pqrt.py
echo "fw=open(fnamelist[0]+'.pqrt','w')" >> pqr_ac2pqrt.py
echo "actype=[]" >> pqr_ac2pqrt.py
echo "for line in frac:" >> pqr_ac2pqrt.py
echo " if (line[:6]=='ATOM ' or line[:6]=='HETATM'):" >> pqr_ac2pqrt.py
echo " at=line[64:74]" >> pqr_ac2pqrt.py
echo " actype.append(at)" >> pqr_ac2pqrt.py
echo "nowi=0" >> pqr_ac2pqrt.py
echo "for line in frpqr:" >> pqr_ac2pqrt.py
echo " if (line[:6]=='ATOM ' or line[:6]=='HETATM'):" >> pqr_ac2pqrt.py
echo " fw.write(line.strip()+actype[nowi]+'\n')" >> pqr_ac2pqrt.py
echo " nowi+=1;" >> pqr_ac2pqrt.py
python pqr_ac2pqrt.py ${basename}.pqr ${basename}_gaff.ac
rm pqr_ac2pqrt.py
if [[ ${outFormat} = "pqrt" && $6 = "" ]];then
rm ${basename}.pqr
fi
fi
outPQRA=""
# Convert to pqra
if [[ $outFormat = "pqra" || $outFormat = "pqrta" ]];then
# Come from ESES_ElementArea.py
echo "#! /usr/bin/env python" > pqr2pqra.py
echo "import os,sys" >> pqr2pqra.py
echo "probe=${probe}" >> pqr2pqra.py
echo "grid=${grid}" >> pqr2pqra.py
echo "solbuffer=${solbuffer}" >> pqr2pqra.py
echo "fname=sys.argv[1]" >> pqr2pqra.py
echo "fnamelist=os.path.splitext(fname)" >> pqr2pqra.py
echo "fxyzr=open(fnamelist[0]+'.xyzr','w')" >> pqr2pqra.py
echo "fr=open(fname)" >> pqr2pqra.py
echo "inlines=fr.readlines();" >> pqr2pqra.py
echo "fr.close();" >> pqr2pqra.py
echo "atomtypes=[];" >> pqr2pqra.py
echo "for line in inlines:" >> pqr2pqra.py
echo " if (line[:4]=='ATOM' or line[:6]=='HETATM'):" >> pqr2pqra.py
echo " tmp=line.split();" >> pqr2pqra.py
echo " element=tmp[-1].upper();" >> pqr2pqra.py
echo " atomtypes.append(element);" >> pqr2pqra.py
echo " radius='%10.5f' % float(line[62:70].strip());" >> pqr2pqra.py
echo " xcoor='%10.5f' % float(line[30:38].strip());" >> pqr2pqra.py
echo " ycoor='%10.5f' % float(line[38:46].strip());" >> pqr2pqra.py
echo " zcoor='%10.5f' % float(line[46:54].strip());" >> pqr2pqra.py
echo " xyzrstr=xcoor+ycoor+zcoor+radius+'\n';" >> pqr2pqra.py
echo " fxyzr.write(xyzrstr);" >> pqr2pqra.py
echo "fxyzr.close()" >> pqr2pqra.py
echo "p=os.popen('./MS_Intersection '+fnamelist[0]+'.xyzr '+str(probe)+' '+str(grid)+' '+str(solbuffer),'r')" >> pqr2pqra.py
echo "totalArea='0'" >> pqr2pqra.py
echo "totalVolume='0'" >> pqr2pqra.py
echo "while 1:" >> pqr2pqra.py
echo " line=p.readline();" >> pqr2pqra.py
echo " if 'area:' in line: totalArea=line.split(':')[1].split()[0]" >> pqr2pqra.py
echo " if 'volume:' in line: totalVolume=line.split(':')[1].split()[0]" >> pqr2pqra.py
echo " if not line:break" >> pqr2pqra.py
echo "fa=open('partition_area.txt')" >> pqr2pqra.py
echo "atomareas=[];# tmp save atom area by atom number" >> pqr2pqra.py
echo "typedefault=['H','C','N','O','F','S','P','CL','BR','I'];" >> pqr2pqra.py
echo "typeareas={'H':0.0,'C':0.0,'N':0.0,'O':0.0,'F':0.0,'S':0.0,'P':0.0,'CL':0.0,'BR':0.0,'I':0.0};" >> pqr2pqra.py
echo "atomnum=0;" >> pqr2pqra.py
echo "for line in fa:" >> pqr2pqra.py
echo " tmp=line.split();" >> pqr2pqra.py
echo " atomarea='%12.6f' % float(tmp[1]);" >> pqr2pqra.py
echo " atomareas.append(atomarea);" >> pqr2pqra.py
echo " atype=atomtypes[atomnum];" >> pqr2pqra.py
echo " typeareas[atype]=typeareas.setdefault(atype,0.0)+float(tmp[1]);" >> pqr2pqra.py
echo " atomnum=atomnum+1;" >> pqr2pqra.py
echo "fa.close()" >> pqr2pqra.py
echo "fwname=fnamelist[0]+'.pqra'" >> pqr2pqra.py
echo "fw=open(fwname,'w')" >> pqr2pqra.py
echo "typeused=['H','C','N','O','F','S','P','CL','BR','I'];" >> pqr2pqra.py
echo "for i in typeareas.iterkeys():" >> pqr2pqra.py
echo " if i not in typeused:typeused.append(i);" >> pqr2pqra.py
echo "outputelearea=fnamelist[0]+' Areas: '+totalArea+' Volumes: '+totalVolume+' ';" >> pqr2pqra.py
echo "fw.write('REMARK AREAS '+totalArea+'\n');" >> pqr2pqra.py
echo "fw.write('REMARK VOLUMES '+totalVolume+'\n'); " >> pqr2pqra.py
echo "for element in typedefault:" >> pqr2pqra.py
echo " fw.write('REMARK AREA '+'%2s'%element+' '+'%20.6f'%typeareas.get(element,0.0)+'\n');" >> pqr2pqra.py
echo " outputelearea=outputelearea+element+': '+str(typeareas[element])+' ';" >> pqr2pqra.py
echo "print outputelearea" >> pqr2pqra.py
echo "fr=open(fname)" >> pqr2pqra.py
echo "atomnum=0;" >> pqr2pqra.py
echo "for line in fr:" >> pqr2pqra.py
echo " if (line[:4]=='ATOM' or line[:6]=='HETATM'):" >> pqr2pqra.py
echo " tmp=line.split();" >> pqr2pqra.py
echo " element=tmp[-1].upper();" >> pqr2pqra.py
echo " newline=line.strip('\n')+atomareas[atomnum]+'\n';" >> pqr2pqra.py
echo " fw.write(newline);" >> pqr2pqra.py
echo " atomnum=atomnum+1;" >> pqr2pqra.py
echo " else:" >> pqr2pqra.py
echo " fw.write(line);" >> pqr2pqra.py
echo "fr.close();" >> pqr2pqra.py
echo "fw.close();" >> pqr2pqra.py
outPQRA=`python pqr2pqra.py ${basename}.pqr`
rm pqr2pqra.py ${basename}.xyzr bounding_box.txt grid_info.txt intersection_info.txt partition_area.txt;
if [[ $6 = "" ]];then
rm ${basename}.pqr
fi
fi
# Convert to pqrta
if [ $outFormat = "pqrta" ];then
echo "#! /usr/bin/env python" > combine2pqrta.py
echo "import os,sys" >> combine2pqrta.py
echo "fname=sys.argv[1]" >> combine2pqrta.py
echo "fnamelist=os.path.splitext(fname)" >> combine2pqrta.py
echo "frpqra=open(fnamelist[0]+'.pqra','r')" >> combine2pqrta.py
echo "frpqrt=open(fnamelist[0]+'.pqrt','r')" >> combine2pqrta.py
echo "fw=open(fnamelist[0]+'.pqrta','w')" >> combine2pqrta.py
echo "atomarea=[]" >> combine2pqrta.py
echo "for line in frpqra:" >> combine2pqrta.py
echo " if (line[:6]=='ATOM ' or line[:6]=='HETATM'):" >> combine2pqrta.py
echo " aarea=line[78:90]" >> combine2pqrta.py
echo " atomarea.append(aarea)" >> combine2pqrta.py
echo " elif (line[:6]=='REMARK'): fw.write(line)">> combine2pqrta.py
echo "nowi=0" >> combine2pqrta.py
echo "for line in frpqrt:" >> combine2pqrta.py
echo " if (line[:6]=='ATOM ' or line[:6]=='HETATM'):" >> combine2pqrta.py
echo " fw.write(line.strip()+atomarea[nowi]+'\n')" >> combine2pqrta.py
echo " nowi+=1;" >> combine2pqrta.py
echo "fw.close();frpqra.close();frpqrt.close();" >> combine2pqrta.py
python combine2pqrta.py ${basename}.pqrt ${basename}.pqra
rm combine2pqrta.py ${basename}.pqrt ${basename}.pqra
fi
rm ${basename}_gaff.* leap.log
fi
Raw
read_atomicvalues.py
#! /usr/bin/env python
# -*- coding: utf8 -*-
#
__author__="Zhixiong Zhao"
__date__="2015-11-30"
__doc__=""" A script to extract the atomic based data to pqr file.
You can define new function as read_atomicSingleValue(filename) to define how to read the data file.
And then define a new jobtype in writejob(fname, jobtype=0, jobname=None) function.
write_pqrx(filename,fout=None,atomvalues=None,keys=None) help to combine the pqr and data to a pqrx file."""
import os,sys
from optparse import OptionParser
def read_atomicSingleValue(filename):
# read the data for atoms
# return a list of data based on atoms sequence
# Datafile: One data in a line for corresponding atom.
try:
SingleValues=[];
fr=open(filename)
for line in fr:
tmp=line.strip().split()
if (len(tmp)>0):
SingleValues.append(float(tmp[0]));
fr.close()
return SingleValues
except IOError as e:
print "Can't find the file "+filename+" since: "+str(e)
return None
except ValueError as e:
print "Value error for file: "+filename+" since: "+str(e)
fr.close()
return None
def read_multipole(filename):
# read the multipole data (Bao Wang *.mda) for atoms
# Return (total_dipole, total_quadrupole, dipoles_list, quadrupoles_list) based on atoms sequence
# Datafile: First line " Atomic Multipole Dipole, Quadrapole moment:"
# Then: atomic_dipole atomic_quadrupole for each atom
# Then: " Total Dipole, Quadrapole moment Refer to 0: "
# Then: total_dopole total_quadrupole
try:
fr=open(filename)
finddata=False
findtotal=False
dipoles=[]
qupoles=[]
moldipole=0.0
molqupole=0.0
for line in fr:
tmp=line.strip().split()
if (findtotal):
moldipole=float(tmp[0])
molqupole=float(tmp[1])
findtotal=False;
break;
if (tmp[0]=="Total"):
finddata=False;
findtotal=True;
continue;
if (finddata):
dipoles.append(float(tmp[0]))
qupoles.append(float(tmp[1]))
continue;
if (tmp[0]=="Atomic"):
finddata=True;
continue;
fr.close()
return moldipole,molqupole,dipoles,qupoles
except IOError as e:
print "Can't find the file "+filename+" since: "+str(e)
return None
def write_pqrx(filename,fout=None,atomvalues=None,keys=None):
# Combine a pqr(filename) and given data in a list(atomvalues) based on atomic sequence
# fout is the output file name
# keys is Some information to be written in REMARK part. Such as total_dipole: value pair.
# It's the *Key* function to write out pqrx file
try:
if (not atomvalues):
print "Empty values for atoms!"
return False
fr=open(filename)
fnamelist=os.path.splitext(filename)
if (not fout):
fout=fnamelist[0]+".pqrx"
fw=open(fout,'w')
if (isinstance(keys,dict)):
for key,value in keys.items():
fw.write('REMARK '+key+' '+str(value)+'\n');
count=0
for line in fr:
if (line[:4]=='ATOM' or line[:6]=='HETATM'):
newline=line.strip('\n')+'%12.6f' % atomvalues[count]+'\n';
fw.write(newline);
count+=1;
else:
fw.write(line);
fr.close();
fw.close();
return True;
except Exception as e:
print "Can't successfully deal with file: "+fout+": "+str(e);
fr.close();fw.close();
return False;
def writejob(fname, jobtype=0, jobname=None):
# Define how to deal with your data file and output it
# You can set new function to process your data file! and be loaded here!
# fname is datafile to be read. such as 1ajj.data (for 1ajj.pqr, same prefix)
# jobtype:
# 0: a file saving value at each line
# 1: a file saving multipoles informations
# Jobname: Middle name in output file
fnamelist=os.path.splitext(fname)
if (jobtype==0):
results=read_atomicSingleValue(fname)
if (jobname):
jobname="_"+jobname;
else:
jobname=""
if (not results is None):
return write_pqrx(fnamelist[0]+".pqr",fout=fnamelist[0]+jobname+".pqrx",atomvalues=results)
elif (jobtype==1):
results=read_multipole(fname);
if (not results is None):
write_pqrx(fnamelist[0]+".pqr",fout=fnamelist[0]+"_dipole.pqrx",atomvalues=results[2], keys={"DIPOLE":results[0]})
return write_pqrx(fnamelist[0]+".pqr",fout=fnamelist[0]+"_qupole.pqrx",atomvalues=results[3], keys={"QUPOLE":results[1]})
else:
print "Unrecognized job type: "+str(jobtype);
return False
if (__name__ == '__main__'):
helpdes="""To read value for atom in pqr and finally generate pqrx file"""
parser = OptionParser(description=helpdes)
parser.add_option("-i", "--input", action="store",
dest="input", default="",
help="Read input data from input file")
parser.add_option("-m", "--multi", action="store",
dest="multi", default="",
help="File containing file name without extension, format must be assigned!")
parser.add_option("-f", "--format", action="store",
dest="format", default="",
help="Input file format (file extension)")
parser.add_option("-j", "--job", action="store",
dest="jobtype", default="0",
help="The job type, default 0 for one value per line.")
parser.add_option("-o", "--outfix", action="store",
dest="outfix", default=None,
help="The middle name for the output file.")
(options, args) = parser.parse_args()
if (len(sys.argv)<2):
print "Please assign a calculated result file!"
parser.print_help()
#parser.print_description()
#parser.print_usage()
exit()
# a file saving id (for id.pqr and id.data) and process all the id-based files!
# test.py -m id.txt -format data -j jobid [-o dataname]
if (options.multi !=""):
if (options.format==""):
print "Using multifile mode. But input format (extension) is not giving!"
exit(1)
fm=open(options.multi)
for line in fm:
tmp=line.strip().split()
if (len(tmp)>0):
fname=tmp[0]+"."+options.format;
writejob(fname,int(options.jobtype),options.outfix);
fm.close()
# Single file to be process.
# test.py -i id.txt -j jobid [-o dataname]
elif (options.input != ""):
fname=options.input;
writejob(fname,int(options.jobtype),options.outfix);
Raw
zap9radius4pqr.py
#! /usr/bin/env python
# -*- coding: utf8 -*-
# Author: Platinhom; Last Updated: 2015-09-09
# Replace the radius of atom by ZAP9 strategy
# Only for pqr format.
# python zap9radius4pqr.py file.pqr
import os,sys
if (__name__ == '__main__'):
fname=sys.argv[1]
fnamelist=os.path.splitext(fname)
fwname=fnamelist[0]+"_zap"+fnamelist[1]
fr=open(fname)
fw=open(fwname,'w')
for line in fr:
if (line[:4]=="ATOM" or line[:6]=="HETATM"):
tmp=line.split();
element=tmp[-1].upper();
radius=" 0.000";
if (element=='H'):
radius=" 1.100";
elif (element=='C'):
radius=" 1.870";
elif (element=='N'):
radius=" 1.400";
elif (element=='O'):
radius=" 1.760";
elif (element=='F'):
radius=" 2.400";
elif (element=='S'):
radius=" 2.150";
elif (element=='CL'):
radius=" 1.820";
newline=line[:62]+radius+line[70:];
fw.write(newline);
else:
fw.write(line);
fr.close()
fw.close()
#end main
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment