Bjoern Olausson

gauss2mol2 PDF Print E-mail
  
Tuesday, 03 March 2009 15:21

  1. #!/usr/bin/python
  2. # -*- coding: utf-8 -*-
  3. #
  4. #--------------------------------------------------------------------------------
  5. #gauss2mol2.py v1.1, Copyright Bjoern Olausson
  6. #--------------------------------------------------------------------------------
  7. #This program is free software; you can redistribute it and/or modify
  8. #it under the terms of the GNU General Public License as published by
  9. #the Free Software Foundation; either version 2 of the License, or
  10. #(at your option) any later version.
  11. #
  12. #This program is distributed in the hope that it will be useful,
  13. #but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. #MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15. #GNU General Public License for more details.
  16. #
  17. #To view the license visit
  18. #http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
  19. #or write to
  20. #Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  21. #--------------------------------------------------------------------------------
  22. #--------------------------------------------------------------------------------
  23. #
  24. #This script uses antechamber to read gaussion log files, and generates mol2 files with resp charges
  25. #Antechamber does not preserve the coordinates an correct atomtypes.
  26. #So in this script you can choose to replace the charges in your original mol2 file with the gaussian charges
  27. #(taken from the antechamber generated mol2 file) or to replace the coordinates and atom types in the antechamber
  28. #generated mol2 file with the coordinates and atom types from your original mol2 file.
  29. #
  30. #run with --help to see how to use this script.
  31. #
  32.  
  33. import sys, os, shutil, string, glob, subprocess, itertools
  34. from optparse import OptionParser
  35. parser = OptionParser()
  36.  
  37. parser.add_option(
  38. "-f",
  39. action="store",
  40. type="string",
  41. dest="F",
  42. help="""File(s) to process. Original mol2 file(s) must be in the same place as the gaussian log file(s) and must have the same prefix (AAA.mol2, AAA.log). If you want to use wildcards, please wrap your expression in quotes (gauss2mol2.py --crd -f "*.log*" )""",
  43. metavar="FILES"
  44. )
  45.  
  46. parser.add_option(
  47. "-c",
  48. "--chg",
  49. action="store_true",
  50. dest="CHG",
  51. help="Replace the charges in your mol2 with the charges from the antechamber generated mol2 file",
  52. metavar="REPCHG"
  53. )
  54.  
  55. parser.add_option(
  56. "-o",
  57. "--crd",
  58. action="store_true",
  59. dest="CRD",
  60. help="Replace the coordinates and atom types in the antechamber generated mol2 file with the coordinates and atom types of your original mol2 file",
  61. metavar="REPCRD"
  62. )
  63.  
  64. # instruct optparse to parse the program's command line
  65. (options, args) = parser.parse_args()
  66.  
  67. # Checking for required otions
  68. if not options.F:
  69. parser.error("-f must be set")
  70. elif not options.CHG and not options.CRD:
  71. parser.error("--chg or --crd must be set")
  72.  
  73. #Convert Options to a list
  74. F = options.F.split(" ")
  75. CHG = options.CHG
  76. CRD = options.CRD
  77.  
  78. #Change to the current working dir
  79. path = os.getcwd()
  80. os.chdir(path)
  81.  
  82. def linefilter(lines, start, stop) :
  83. lines = iter(lines) #state-preserving iterator
  84. #if 'start' matches BOF, we will start to read at the beginning of the file,
  85. #else we start at the following line that matches 'start'
  86. if start != "BOF" :
  87. for line in lines :
  88. if line == start :
  89. break
  90. for line in lines :
  91. #if 'stop' matches EOF, we will read to the end of file.
  92. #else we stop reading on the line that matches 'stop'
  93. if stop == "EOF" :
  94. yield line
  95. else :
  96. if line ==  stop :
  97. break
  98. else :
  99. yield line
  100.  
  101. #If "*" is found in a parameter, glob it.
  102. for element in F :
  103. if '*' in element :
  104. F.extend(glob.glob(element))
  105. F.remove(element)
  106.  
  107.  
  108. def altercrd(gaufile, orgfile) :
  109. #open the new file
  110. fout = open(IN+"-CRD.mol2", "w")
  111.  
  112. #write all from  beginning to '@<TRIPOS>ATOM' from the gaussian mol2 to the new file
  113. top = linefilter(gaufile.splitlines(), 'BOF', '@<TRIPOS>ATOM')
  114. for line in linefilter(gaufile.splitlines(), 'BOF', '@<TRIPOS>ATOM') :
  115. fout.writelines(line+"
    "
    )
  116.  
  117. #Write the '@<TRIPOS>ATOM' match string cause i is stripped in the linefilter
  118. fout.writelines('@<TRIPOS>ATOM'+"
    "
    )
  119.  
  120. #Now lets write the collumns 0,1 from the gaussian mol2file, the collumn 2,3,4 and 5
  121. #from the original mol2 files and again the collumns 6,7 and 8 from the gaussian mol2 file
  122. #So it should look like this
  123. #
  124. #ORIGINAL FILE EXSAMPLE:
  125. #@<TRIPOS>ATOM
  126. #      1 N1        -18.5878   -6.4372    8.0759 N.pl3     1 <1>        -0.9900
  127. #      2 H1        -18.4216   -6.0346    8.9973 H         1 <1>         0.3600
  128. #@<TRIPOS>BOND
  129.  
  130. #GAUSSIAN FILE EXAMPLE:
  131. #@<TRIPOS>ATOM
  132. #      1 N1          5.1000    0.1990   -0.0330 N.3       1 MOL      -1.0597
  133. #      2 H1          5.9140   -0.3800   -0.1080 H         1 MOL       0.3916
  134. #@<TRIPOS>BOND
  135.  
  136. #Merge table
  137. #G == gaussion mol2 fil
  138. #O == original mol2 files
  139. #@<TRIPOS>ATOM
  140. #      G G          O    O   O O       G G   G
  141. #      G G          O    O   O O       G G   G
  142. #@<TRIPOS>BOND
  143.  
  144. #MERGED FILE EXAMPLE:
  145. #@<TRIPOS>ATOM
  146. #      1 N1        -18.5878   -6.4372    8.0759 N.pl3     1 MOL      -1.0597
  147. #      2 H1        -18.4216   -6.0346    8.9973 H         1 MOL       0.3916
  148. #@<TRIPOS>BOND
  149.  
  150.  
  151. gau = linefilter(gaufile.splitlines(), '@<TRIPOS>ATOM', '@<TRIPOS>BOND')
  152. org = linefilter(orgfile.splitlines(), '@<TRIPOS>ATOM', '@<TRIPOS>BOND')
  153.  
  154.  
  155. for gauline, orgline in itertools.izip(gau, org) :
  156. glist = gauline.split()
  157. olist = orgline.split()
  158.  
  159. anbr = glist[0]
  160. anbrF = anbr.rjust(7," ")+" "
  161.  
  162. atom = glist[1]
  163. atomF = atom.ljust(10," ")+" "
  164.  
  165. x = olist[2]
  166. xF = x.rjust(10," ")+" "
  167.  
  168. y = olist[3]
  169. yF = y.rjust(10," ")+" "
  170.  
  171. z = olist[4]
  172. zF = z.rjust(10," ")+" "
  173.  
  174. atype = olist[5]
  175. atypeF = atype.ljust(10," ")+" "
  176.  
  177. molnbr = glist[6]
  178. molnbrF = molnbr.rjust(10," ")+" "
  179.  
  180. molnam = glist[7]
  181. molnamF = " "+molnam.ljust(10," ")+" "
  182.  
  183. charge = glist[8]
  184. chargeF = charge.rjust(10," ")+" "
  185.  
  186. fout.write(anbrF+atomF+xF+yF+zF+atypeF+molnbrF+molnamF+chargeF+"
    "
    )
  187.  
  188. #Here we write the rest of the gaussian mol2 file to the new file
  189. #Again we have to write the '@<TRIPOS>BOND' manually cause it is stripped by linefilter
  190. fout.writelines('@<TRIPOS>BOND'+"
    "
    )
  191. for line in linefilter(gaufile.splitlines(), '@<TRIPOS>BOND', 'EOF') :
  192. fout.writelines(line+"
    "
    )
  193. #finally save the new mol2 file.
  194. fout.close()
  195.  
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202.  
  203. def alterchg(gaufile, orgfile) :
  204. sys.stdout.flush()
  205. #open the new file
  206. fout = open(IN+"-CHG.mol2", "w")
  207.  
  208. #write all from  beginning to '@<TRIPOS>ATOM' from the gaussian mol2 to the new file
  209. for line in linefilter(orgfile.splitlines(), 'BOF', '@<TRIPOS>ATOM') :
  210. fout.writelines(line+"
    "
    )
  211.  
  212. #Write the '@<TRIPOS>ATOM' match string cause i is stripped in the linefilter
  213. fout.writelines('@<TRIPOS>ATOM'+"
    "
    )
  214.  
  215. #Now lets write the collumns 0 to 7 from the original mol2file and take collumn 8
  216. #from the gassian files
  217. #So it should look like this
  218. #
  219. #ORIGINAL FILE EXSAMPLE:
  220. #@<TRIPOS>ATOM
  221. #      1 N1        -18.5878   -6.4372    8.0759 N.pl3     1 <1>        -0.9900
  222. #      2 H1        -18.4216   -6.0346    8.9973 H         1 <1>         0.3600
  223. #@<TRIPOS>BOND
  224.  
  225. #GAUSSIAN FILE EXAMPLE:
  226. #@<TRIPOS>ATOM
  227. #      1 N1          5.1000    0.1990   -0.0330 N.3       1 MOL      -1.0597
  228. #      2 H1          5.9140   -0.3800   -0.1080 H         1 MOL       0.3916
  229. #@<TRIPOS>BOND
  230.  
  231. #Merge table
  232. #G == gaussion mol2 fil
  233. #O == original mol2 files
  234. #@<TRIPOS>ATOM
  235. #      O O          O    O   O O       O O   G
  236. #      O O          O    O   O O       O O   G
  237. #@<TRIPOS>BOND
  238.  
  239. #MERGED FILE EXAMPLE:
  240. #@<TRIPOS>ATOM
  241. #      1 N1        -18.5878   -6.4372    8.0759 N.pl3     1 <1>      -1.0597
  242. #      2 H1        -18.4216   -6.0346    8.9973 H         1 <1>       0.3916
  243. #@<TRIPOS>BOND
  244.  
  245.  
  246. gau = linefilter(gaufile.splitlines(), '@<TRIPOS>ATOM', '@<TRIPOS>BOND')
  247. org = linefilter(orgfile.splitlines(), '@<TRIPOS>ATOM', '@<TRIPOS>BOND')
  248.  
  249.  
  250. for gauline, orgline in itertools.izip(gau, org) :
  251. glist = gauline.split()
  252. olist = orgline.split()
  253.  
  254. anbr = olist[0]
  255. anbrF = anbr.rjust(7," ")+" "
  256.  
  257. atom = olist[1]
  258. atomF = atom.ljust(10," ")+" "
  259.  
  260. x = olist[2]
  261. xF = x.rjust(10," ")+" "
  262.  
  263. y = olist[3]
  264. yF = y.rjust(10," ")+" "
  265.  
  266. z = olist[4]
  267. zF = z.rjust(10," ")+" "
  268.  
  269. atype = olist[5]
  270. atypeF = atype.ljust(10," ")+" "
  271.  
  272. molnbr = olist[6]
  273. molnbrF = molnbr.rjust(10," ")+" "
  274.  
  275. molnam = olist[7]
  276. molnamF = " "+molnam.ljust(10," ")+" "
  277.  
  278. charge = glist[8]
  279. chargeF = charge.rjust(10," ")+" "
  280.  
  281. fout.write(anbrF+atomF+xF+yF+zF+atypeF+molnbrF+molnamF+chargeF+"
    "
    )
  282.  
  283. #Here we write the rest of the gaussian mol2 file to the new file
  284. #Again we have to write the '@<TRIPOS>BOND' manually cause it is stripped by linefilter
  285. fout.writelines('@<TRIPOS>BOND'+"
    "
    )
  286. for line in linefilter(orgfile.splitlines(), '@<TRIPOS>BOND', 'EOF') :
  287. fout.writelines(line+"
    "
    )
  288. #finally save the new mol2 file.
  289. fout.close()
  290.  
  291.  
  292. #remove extension from elements, append .log and .mol2, and generate with antechamber
  293. #sybyl conform mol2 files with resp charges calculated from the gaussian log files
  294. for element in F :
  295. IN = os.path.splitext(element)[0]
  296. sys.stdout.write("\rantechamber -pf y -at sybyl -i %s.log -fi gout -o %s-g.mol2 -fo mol2 -c resp
    "
    % (IN, IN))
  297. sys.stdout.flush()   
  298. amberout = subprocess.Popen(["antechamber", "-i", IN+".log", "-fi", "gout", "-o", IN+"-g.mol2", "-fo", "mol2", "-c", "resp", "-pf", "y", "-at", "sybyl"], stderr=subprocess.PIPE, stdout=subprocess.PIPE)
  299.  
  300. stdout, stderr = amberout.communicate()
  301. if stderr :
  302. print stderr
  303. print stdout
  304.  
  305. tmpfiles = ["esout", "punch", "qout", "QOUT"]
  306. for remfile in tmpfiles :
  307. os.remove(remfile)
  308. #open the gaussian and original mol2 file
  309. gaumol2 = open(IN+"-g.mol2", "r").read()
  310. orgmol2 = open(IN+".mol2", "r").read()
  311. if CRD :
  312. altercrd(gaumol2,orgmol2)
  313.  
  314. if CHG :
  315. alterchg(gaumol2,orgmol2)
  316.  
  317. print "If you want to make a amber prep file from your mol2 files, you can do this with"
  318. print " # for i in *.mol2 ; do antechamber -pf y -fi mol2 -i ${i} -fo prepi -o ${i/.mol2/.prep} -at amber ; done"
  319. print "Be shure to check the Atomnames. Numbering should start with 1"
  320. print "If you have missing parameters, run parmcheck to generate a frcmod file with the missing parameters"
  321. print " # for i in *.prep ; do parmchk -i ${i} -f prepi -o ${i/.prep/.frcmod} ; done"
  322.  

Attachments:
FileFile sizeLast Modified
Download this file (gauss2mol2.py)gauss2mol2.py10 Kb07/09/09 10:48
Last Updated ( Thursday, 09 July 2009 10:48 )
 

Add comment


Security code
Refresh

6M Pixel

Banner

www. is deprecated

Banner

Play OGG

Banner