Coverage for cclib/parser/turbomoleparser.py : 35%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2#
3# Copyright (c) 2020, the cclib development team
4#
5# This file is part of cclib (http://cclib.github.io) and is distributed under
6# the terms of the BSD 3-Clause License.
8"""Parser for Turbomole output files."""
10import re
12import numpy
14from cclib.parser import logfileparser
15from cclib.parser import utils
17class AtomBasis:
18 def __init__(self, atname, basis_name, inputfile):
19 self.symmetries=[]
20 self.coefficients=[]
21 self.atname=atname
22 self.basis_name=basis_name
24 self.parse_basis(inputfile)
26 def parse_basis(self, inputfile):
27 i=0
28 line=inputfile.next()
30 while(line[0]!="*"):
31 (nbasis_text, symm)=line.split()
32 self.symmetries.append(symm)
34 nbasis=int(nbasis_text)
35 coeff_arr=numpy.zeros((nbasis, 2), float)
37 for j in range(0, nbasis, 1):
38 line=inputfile.next()
39 (e1_text, e2_text)=line.split()
40 coeff_arr[j][0]=float(e1_text)
41 coeff_arr[j][1]=float(e2_text)
43 self.coefficients.append(coeff_arr)
44 line=inputfile.next()
46class Turbomole(logfileparser.Logfile):
47 """A Turbomole log file."""
49 def __init__(self, *args, **kwargs):
50 super(Turbomole, self).__init__(logname="Turbomole", *args, **kwargs)
52 def __str__(self):
53 """Return a string representation of the object."""
54 return "Turbomole output file %s" % (self.filename)
56 def __repr__(self):
57 """Return a representation of the object."""
58 return 'Turbomole("%s")' % (self.filename)
60 def normalisesym(self, label):
61 """Normalise the symmetries used by Turbomole.
63 The labels are standardized except for the first character being lowercase.
64 """
65 # TODO more work could be required, but we don't have any logfiles
66 # with non-C1 symmetry.
67 return label[0].upper() + label[1:]
69 def before_parsing(self):
70 self.geoopt = False # Is this a GeoOpt? Needed for SCF targets/values.
71 self.periodic_table = utils.PeriodicTable()
73 @staticmethod
74 def split_molines(inline):
75 """Splits the lines containing mocoeffs (each of length 20)
76 and converts them to float correctly.
77 """
78 line = inline.replace("D", "E")
79 f1 = line[0:20]
80 f2 = line[20:40]
81 f3 = line[40:60]
82 f4 = line[60:80]
84 if(len(f4) > 1):
85 return [float(f1), float(f2), float(f3), float(f4)]
86 if(len(f3) > 1):
87 return [float(f1), float(f2), float(f3)]
88 if(len(f2) > 1):
89 return [float(f1), float(f2)]
90 if(len(f1) > 1):
91 return [float(f1)]
93 def extract(self, inputfile, line):
94 """Extract information from the file object inputfile."""
96 ## This information is in the control file.
97 # $rundimensions
98 # dim(fock,dens)=1860
99 # natoms=20
100 # nshell=40
101 # nbf(CAO)=60
102 # nbf(AO)=60
103 # dim(trafo[SAO<-->AO/CAO])=60
104 # rhfshells=1
105 if line[3:10]=="natoms=":
106 self.natom=int(line[10:])
108 if line[3:11] == "nbf(AO)=":
109 nmo = int(line.split('=')[1])
110 self.set_attribute('nbasis', nmo)
111 self.set_attribute('nmo', nmo)
113 # Extract the version number and optionally the build number.
114 searchstr = ": TURBOMOLE"
115 index = line.find(searchstr)
116 if index > -1:
117 line = line[index + len(searchstr):]
118 tokens = line.split()
119 package_version = tokens[0][1:].replace("-", ".")
120 self.metadata["package_version"] = package_version
121 self.metadata["legacy_package_version"] = package_version
122 if tokens[1] == "(":
123 revision = tokens[2]
124 self.metadata["package_version"] = "{}.r{}".format(package_version, revision)
126 ## Atomic coordinates in job.last:
127 # +--------------------------------------------------+
128 # | Atomic coordinate, charge and isotop information |
129 # +--------------------------------------------------+
130 #
131 #
132 # atomic coordinates atom shells charge pseudo isotop
133 # -2.69176330 -0.00007129 -0.44712612 c 3 6.000 0 0
134 # -1.69851645 -0.00007332 2.06488947 c 3 6.000 0 0
135 # 0.92683848 -0.00007460 2.49592179 c 3 6.000 0 0
136 # 2.69176331 -0.00007127 0.44712612 c 3 6.000 0 0
137 # 1.69851645 -0.00007331 -2.06488947 c 3 6.000 0 0
138 #...
139 # -7.04373606 0.00092244 2.74543891 h 1 1.000 0 0
140 # -9.36352819 0.00017229 0.07445322 h 1 1.000 0 0
141 # -0.92683849 -0.00007461 -2.49592179 c 3 6.000 0 0
142 # -1.65164853 -0.00009927 -4.45456858 h 1 1.000 0 0
143 if 'Atomic coordinate, charge and isotop information' in line:
144 while 'atomic coordinates' not in line:
145 line = next(inputfile)
147 atomcoords = []
148 atomnos = []
149 line = next(inputfile)
150 while len(line) > 2:
151 atomnos.append(self.periodic_table.number[line.split()[3].capitalize()])
152 atomcoords.append([utils.convertor(float(x), "bohr", "Angstrom")
153 for x in line.split()[:3]])
154 line = next(inputfile)
156 self.append_attribute('atomcoords', atomcoords)
157 self.set_attribute('atomnos', atomnos)
158 self.set_attribute('natom', len(atomcoords))
160 # Frequency values in aoforce.out
161 # mode 7 8 9 10 11 12
162 #
163 # frequency 53.33 88.32 146.85 171.70 251.75 289.44
164 #
165 # symmetry a a a a a a
166 #
167 # IR YES YES YES YES YES YES
168 # |dDIP/dQ| (a.u.) 0.0002 0.0000 0.0005 0.0004 0.0000 0.0000
169 # intensity (km/mol) 0.05 0.00 0.39 0.28 0.00 0.00
170 # intensity ( % ) 0.05 0.00 0.40 0.28 0.00 0.00
171 #
172 # RAMAN YES YES YES YES YES YES
173 #
174 # 1 c x 0.00000 0.00001 0.00000 -0.01968 -0.04257 0.00001
175 # y -0.08246 -0.08792 0.02675 -0.00010 0.00000 0.17930
176 # z 0.00001 0.00003 0.00004 -0.10350 0.11992 -0.00003
177 # ....
178 #
179 # reduced mass(g/mol) 3.315 2.518 2.061 3.358 3.191 2.323
181 if 'NORMAL MODES and VIBRATIONAL FREQUENCIES (cm**(-1))' in line:
182 has_raman = False
183 while line[7:11] != 'mode':
184 line = next(inputfile)
185 if line.startswith(" differential RAMAN cross sections"):
186 has_raman = True
187 vibfreqs, vibsyms, vibirs, vibdisps, vibrmasses = [], [], [], [], []
188 while 'all done ****' not in line:
190 if line.strip().startswith('mode'):
191 self.skip_line(inputfile, 'b')
192 line = next(inputfile)
193 assert line.strip().startswith('frequency')
194 freqs = [float(i.replace('i', '-')) for i in line.split()[1:]]
195 vibfreqs.extend(freqs)
196 self.skip_lines(inputfile, ['b'])
197 line = next(inputfile)
198 assert line.strip().startswith('symmetry')
199 syms = line.split()[1:]
200 vibsyms.extend(syms)
202 self.skip_lines(inputfile, ['b', 'IR', 'dDIP/dQ'])
203 line = next(inputfile)
204 assert line.strip().startswith('intensity (km/mol)')
205 irs = [utils.float(f) for f in line.split()[2:]]
206 vibirs.extend(irs)
208 self.skip_lines(inputfile, ['intensity %', 'b', 'RAMAN'])
209 if has_raman:
210 self.skip_lines(
211 inputfile,
212 ['(par,par)', '(ort,ort)', '(ort,unpol)', 'depol. ratio']
213 )
214 line = next(inputfile)
215 assert not line.strip()
216 line = next(inputfile)
217 x, y, z = [], [], []
218 atomcounter = 0
219 while atomcounter < self.natom:
220 x.append([float(i) for i in line.split()[3:]])
221 line = next(inputfile)
222 y.append([float(i) for i in line.split()[1:]])
223 line = next(inputfile)
224 z.append([float(i) for i in line.split()[1:]])
225 line = next(inputfile)
226 atomcounter += 1
228 for j in range(len(x[0])):
229 disps = []
230 for i in range(len(x)):
231 disps.append([x[i][j], y[i][j], z[i][j]])
232 vibdisps.append(disps)
234 line = next(inputfile)
235 assert line.startswith('reduced mass(g/mol)')
236 rmasses = [utils.float(f) for f in line.split()[2:]]
237 vibrmasses.extend(rmasses)
239 line = next(inputfile)
241 self.set_attribute('vibfreqs', vibfreqs)
242 self.set_attribute('vibsyms', vibsyms)
243 self.set_attribute('vibirs', vibirs)
244 self.set_attribute('vibdisps', vibdisps)
245 self.set_attribute('vibrmasses', vibrmasses)
247 # In this section we are parsing mocoeffs and moenergies from
248 # the files like: mos, alpha and beta.
249 # $scfmo scfconv=6 format(4d20.14)
250 # # SCF total energy is -382.3457535740 a.u.
251 # #
252 # 1 a eigenvalue=-.97461484059799D+01 nsaos=60
253 # 0.69876828353937D+000.32405121159405D-010.87670894913921D-03-.85232349313288D-07
254 # 0.19361534257922D-04-.23841194890166D-01-.81711001390807D-020.13626356942047D-02
255 # ...
256 # ...
257 # $end
258 if (line.startswith('$scfmo') or line.startswith('$uhfmo')) and line.find('scfconv') > 0:
259 if line.strip().startswith('$uhfmo_alpha'):
260 self.unrestricted = True
262 # Need to skip the first line to start with lines starting with '#'.
263 line = next(inputfile)
264 while line.strip().startswith('#') and not line.find('eigenvalue') > 0:
265 line = next(inputfile)
267 moenergies = []
268 mocoeffs = []
270 while not line.strip().startswith('$'):
271 info = re.match(r".*eigenvalue=(?P<moenergy>[0-9D\.+-]{20})\s+nsaos=(?P<count>\d+).*", line)
272 eigenvalue = utils.float(info.group('moenergy'))
273 orbital_energy = utils.convertor(eigenvalue, 'hartree', 'eV')
274 moenergies.append(orbital_energy)
275 single_coeffs = []
276 nsaos = int(info.group('count'))
278 while(len(single_coeffs) < nsaos):
279 line = next(inputfile)
280 single_coeffs.extend(Turbomole.split_molines(line))
282 mocoeffs.append(single_coeffs)
283 line = next(inputfile)
285 max_nsaos = max([len(i) for i in mocoeffs])
286 for i in mocoeffs:
287 while len(i) < max_nsaos:
288 i.append(numpy.nan)
290 if not hasattr(self, 'mocoeffs'):
291 self.mocoeffs = []
293 if not hasattr(self, 'moenergies'):
294 self.moenergies = []
296 self.mocoeffs.append(mocoeffs)
297 self.moenergies.append(moenergies)
299 # Parsing the scfenergies, scfvalues and scftargets from job.last file.
300 # scf convergence criterion : increment of total energy < .1000000D-05
301 # and increment of one-electron energy < .1000000D-02
302 #
303 # ...
304 # ...
305 # current damping : 0.700
306 # ITERATION ENERGY 1e-ENERGY 2e-ENERGY NORM[dD(SAO)] TOL
307 # 1 -382.34543727790 -1396.8009423 570.56292464 0.000D+00 0.556D-09
308 # Exc = -57.835278090846 N = 69.997494722
309 # max. resid. norm for Fia-block= 2.782D-05 for orbital 33a
310 # ...
311 # ...
312 # current damping : 0.750
313 # ITERATION ENERGY 1e-ENERGY 2e-ENERGY NORM[dD(SAO)] TOL
314 # 3 -382.34575357399 -1396.8009739 570.56263988 0.117D-03 0.319D-09
315 # Exc = -57.835593208072 N = 69.999813370
316 # max. resid. norm for Fia-block= 7.932D-06 for orbital 33a
317 # max. resid. fock norm = 8.105D-06 for orbital 33a
318 #
319 # convergence criteria satisfied after 3 iterations
320 #
321 #
322 # ------------------------------------------
323 # | total energy = -382.34575357399 |
324 # ------------------------------------------
325 # : kinetic energy = 375.67398458525 :
326 # : potential energy = -758.01973815924 :
327 # : virial theorem = 1.98255043001 :
328 # : wavefunction norm = 1.00000000000 :
329 # ..........................................
330 if 'scf convergence criterion' in line:
331 total_energy_threshold = utils.float(line.split()[-1])
332 one_electron_energy_threshold = utils.float(next(inputfile).split()[-1])
333 scftargets = [total_energy_threshold, one_electron_energy_threshold]
334 self.append_attribute('scftargets', scftargets)
335 iter_energy = []
336 iter_one_elec_energy = []
337 while 'convergence criteria satisfied' not in line:
338 if 'ITERATION ENERGY' in line:
339 line = next(inputfile)
340 info = line.split()
341 iter_energy.append(utils.float(info[1]))
342 iter_one_elec_energy.append(utils.float(info[2]))
343 line = next(inputfile)
345 assert len(iter_energy) == len(iter_one_elec_energy), \
346 'Different number of values found for total energy and one electron energy.'
347 scfvalues = [[x - y, a - b] for x, y, a, b in
348 zip(iter_energy[1:], iter_energy[:-1], iter_one_elec_energy[1:], iter_one_elec_energy[:-1])]
349 self.append_attribute('scfvalues', scfvalues)
350 while 'total energy' not in line:
351 line = next(inputfile)
353 scfenergy = utils.convertor(utils.float(line.split()[4]), 'hartree', 'eV')
354 self.append_attribute('scfenergies', scfenergy)
356 # **********************************************************************
357 # * *
358 # * RHF energy : -74.9644564256 *
359 # * MP2 correlation energy (doubles) : -0.0365225363 *
360 # * *
361 # * Final MP2 energy : -75.0009789619 *
362 # ...
363 # * Norm of MP1 T2 amplitudes : 0.0673494687 *
364 # * *
365 # **********************************************************************
366 # OR
367 # **********************************************************************
368 # * *
369 # * RHF energy : -74.9644564256 *
370 # * correlation energy : -0.0507799360 *
371 # * *
372 # * Final CCSD energy : -75.0152363616 *
373 # * *
374 # * D1 diagnostic : 0.0132 *
375 # * *
376 # **********************************************************************
377 if 'C C S D F 1 2 P R O G R A M' in line:
378 while 'ccsdf12 : all done' not in line:
379 if 'Final MP2 energy' in line:
380 mp2energy = [utils.convertor(utils.float(line.split()[5]), 'hartree', 'eV')]
381 self.append_attribute('mpenergies', mp2energy)
383 if 'Final CCSD energy' in line:
384 self.append_attribute(
385 'ccenergies',
386 utils.convertor(utils.float(line.split()[5]), 'hartree', 'eV')
387 )
389 line = next(inputfile)
391 # *****************************************************
392 # * *
393 # * SCF-energy : -74.49827196840999 *
394 # * MP2-energy : -0.19254365976227 *
395 # * total : -74.69081562817226 *
396 # * *
397 # * (MP2-energy evaluated from T2 amplitudes) *
398 # * *
399 # *****************************************************
400 if 'm p g r a d - program' in line:
401 while 'ccsdf12 : all done' not in line:
402 if 'MP2-energy' in line:
403 line = next(inputfile)
404 if 'total' in line:
405 mp2energy = [utils.convertor(utils.float(line.split()[3]), 'hartree', 'eV')]
406 self.append_attribute('mpenergies', mp2energy)
407 line = next(inputfile)
409 def deleting_modes(self, vibfreqs, vibdisps, vibirs, vibrmasses):
410 """Deleting frequencies relating to translations or rotations"""
411 i = 0
412 while i < len(vibfreqs):
413 if vibfreqs[i] == 0.0:
414 # Deleting frequencies that have value 0 since they
415 # do not correspond to vibrations.
416 del vibfreqs[i], vibdisps[i], vibirs[i], vibrmasses[i]
417 i -= 1
418 i += 1
420 def after_parsing(self):
421 if hasattr(self, 'vibfreqs'):
422 self.deleting_modes(self.vibfreqs, self.vibdisps, self.vibirs, self.vibrmasses)
425class OldTurbomole(logfileparser.Logfile):
426 """A Turbomole output file. Code is outdated and is not being used."""
428 def __init__(self, *args):
429 # Call the __init__ method of the superclass
430 super(Turbomole, self).__init__(logname="Turbomole", *args)
432 def __str__(self):
433 """Return a string representation of the object."""
434 return "Turbomole output file %s" % (self.filename)
436 def __repr__(self):
437 """Return a representation of the object."""
438 return 'Turbomole("%s")' % (self.filename)
440 def atlist(self, atstr):
441 # turn atstr from atoms section into array
443 fields=atstr.split(',')
444 list=[]
445 for f in fields:
446 if(f.find('-')!=-1):
447 rangefields=f.split('-')
448 start=int(rangefields[0])
449 end=int(rangefields[1])
451 for j in range(start, end+1, 1):
452 list.append(j-1)
453 else:
454 list.append(int(f)-1)
455 return(list)
457 def normalisesym(self, label):
458 """Normalise the symmetries used by Turbomole."""
459 return ans
461 def before_parsing(self):
462 self.geoopt = False # Is this a GeoOpt? Needed for SCF targets/values.
464 def split_molines(self, inline):
465 line=inline.replace("D", "E")
466 f1=line[0:20]
467 f2=line[20:40]
468 f3=line[40:60]
469 f4=line[60:80]
471 if(len(f4)>1):
472 return( (float(f1), float(f2), float(f3), float(f4)) )
473 if(len(f3)>1):
474 return( (float(f1), float(f2), float(f3)) )
475 if(len(f2)>1):
476 return( (float(f1), float(f2)) )
477 if(len(f1)>1):
478 return([float(f1)])
479 return
481 def extract(self, inputfile, line):
482 """Extract information from the file object inputfile."""
484 if line[3:11]=="nbf(AO)=":
485 nmo=int(line[11:])
486 self.nbasis=nmo
487 self.nmo=nmo
488 if line[3:9]=="nshell":
489 temp=line.split('=')
490 homos=int(temp[1])
492 if line[0:6] == "$basis":
493 print("Found basis")
494 self.basis_lib=[]
495 line = inputfile.next()
496 line = inputfile.next()
498 while line[0] != '*' and line[0] != '$':
499 temp=line.split()
500 line = inputfile.next()
501 while line[0]=="#":
502 line = inputfile.next()
503 self.basis_lib.append(AtomBasis(temp[0], temp[1], inputfile))
504 line = inputfile.next()
505 if line == "$ecp\n":
506 self.ecp_lib=[]
508 line = inputfile.next()
509 line = inputfile.next()
511 while line[0] != '*' and line[0] != '$':
512 fields=line.split()
513 atname=fields[0]
514 ecpname=fields[1]
515 line = inputfile.next()
516 line = inputfile.next()
517 fields=line.split()
518 ncore = int(fields[2])
520 while line[0] != '*':
521 line = inputfile.next()
522 self.ecp_lib.append([atname, ecpname, ncore])
524 if line[0:6] == "$coord":
525 if line[0:11] == "$coordinate":
526# print "Breaking"
527 return
529# print "Found coords"
530 self.atomcoords = []
531 self.atomnos = []
532 atomcoords = []
533 atomnos = []
535 line = inputfile.next()
536 if line[0:5] == "$user":
537# print "Breaking"
538 return
540 while line[0] != "$":
541 temp = line.split()
542 atsym=temp[3].capitalize()
543 atomnos.append(self.table.number[atsym])
544 atomcoords.append([utils.convertor(float(x), "bohr", "Angstrom")
545 for x in temp[0:3]])
546 line = inputfile.next()
547 self.atomcoords.append(atomcoords)
548 self.atomnos = numpy.array(atomnos, "i")
550 if line[14:32] == "atomic coordinates":
551 atomcoords = []
552 atomnos = []
554 line = inputfile.next()
556 while len(line) > 2:
557 temp = line.split()
558 atsym = temp[3].capitalize()
559 atomnos.append(self.table.number[atsym])
560 atomcoords.append([utils.convertor(float(x), "bohr", "Angstrom")
561 for x in temp[0:3]])
562 line = inputfile.next()
564 if not hasattr(self,"atomcoords"):
565 self.atomcoords = []
567 self.atomcoords.append(atomcoords)
568 self.atomnos = numpy.array(atomnos, "i")
570 if line[0:6] == "$atoms":
571 print("parsing atoms")
572 line = inputfile.next()
573 self.atomlist=[]
574 while line[0]!="$":
575 temp=line.split()
576 at=temp[0]
577 atnosstr=temp[1]
578 while atnosstr[-1] == ",":
579 line = inputfile.next()
580 temp=line.split()
581 atnosstr=atnosstr+temp[0]
582# print "Debug:", atnosstr
583 atlist=self.atlist(atnosstr)
585 line = inputfile.next()
587 temp=line.split()
588# print "Debug basisname (temp):",temp
589 basisname=temp[2]
590 ecpname=''
591 line = inputfile.next()
592 while(line.find('jbas')!=-1 or line.find('ecp')!=-1 or
593 line.find('jkbas')!=-1):
594 if line.find('ecp')!=-1:
595 temp=line.split()
596 ecpname=temp[2]
597 line = inputfile.next()
599 self.atomlist.append( (at, basisname, ecpname, atlist))
601# I have no idea what this does, so "comment" out
602 if line[3:10]=="natoms=":
603# if 0:
605 self.natom=int(line[10:])
607 basistable=[]
609 for i in range(0, self.natom, 1):
610 for j in range(0, len(self.atomlist), 1):
611 for k in range(0, len(self.atomlist[j][3]), 1):
612 if self.atomlist[j][3][k]==i:
613 basistable.append((self.atomlist[j][0],
614 self.atomlist[j][1],
615 self.atomlist[j][2]))
616 self.aonames=[]
617 counter=1
618 for a, b, c in basistable:
619 ncore=0
620 if len(c) > 0:
621 for i in range(0, len(self.ecp_lib), 1):
622 if self.ecp_lib[i][0]==a and \
623 self.ecp_lib[i][1]==c:
624 ncore=self.ecp_lib[i][2]
626 for i in range(0, len(self.basis_lib), 1):
627 if self.basis_lib[i].atname==a and self.basis_lib[i].basis_name==b:
628 pa=a.capitalize()
629 basis=self.basis_lib[i]
631 s_counter=1
632 p_counter=2
633 d_counter=3
634 f_counter=4
635 g_counter=5
636# this is a really ugly piece of code to assign the right labels to
637# basis functions on atoms with an ecp
638 if ncore == 2:
639 s_counter=2
640 elif ncore == 10:
641 s_counter=3
642 p_counter=3
643 elif ncore == 18:
644 s_counter=4
645 p_counter=4
646 elif ncore == 28:
647 s_counter=4
648 p_counter=4
649 d_counter=4
650 elif ncore == 36:
651 s_counter=5
652 p_counter=5
653 d_counter=5
654 elif ncore == 46:
655 s_counter=5
656 p_counter=5
657 d_counter=6
659 for j in range(0, len(basis.symmetries), 1):
660 if basis.symmetries[j]=='s':
661 self.aonames.append("%s%d_%d%s" % \
662 (pa, counter, s_counter, "S"))
663 s_counter=s_counter+1
664 elif basis.symmetries[j]=='p':
665 self.aonames.append("%s%d_%d%s" % \
666 (pa, counter, p_counter, "PX"))
667 self.aonames.append("%s%d_%d%s" % \
668 (pa, counter, p_counter, "PY"))
669 self.aonames.append("%s%d_%d%s" % \
670 (pa, counter, p_counter, "PZ"))
671 p_counter=p_counter+1
672 elif basis.symmetries[j]=='d':
673 self.aonames.append("%s%d_%d%s" % \
674 (pa, counter, d_counter, "D 0"))
675 self.aonames.append("%s%d_%d%s" % \
676 (pa, counter, d_counter, "D+1"))
677 self.aonames.append("%s%d_%d%s" % \
678 (pa, counter, d_counter, "D-1"))
679 self.aonames.append("%s%d_%d%s" % \
680 (pa, counter, d_counter, "D+2"))
681 self.aonames.append("%s%d_%d%s" % \
682 (pa, counter, d_counter, "D-2"))
683 d_counter=d_counter+1
684 elif basis.symmetries[j]=='f':
685 self.aonames.append("%s%d_%d%s" % \
686 (pa, counter, f_counter, "F 0"))
687 self.aonames.append("%s%d_%d%s" % \
688 (pa, counter, f_counter, "F+1"))
689 self.aonames.append("%s%d_%d%s" % \
690 (pa, counter, f_counter, "F-1"))
691 self.aonames.append("%s%d_%d%s" % \
692 (pa, counter, f_counter, "F+2"))
693 self.aonames.append("%s%d_%d%s" % \
694 (pa, counter, f_counter, "F-2"))
695 self.aonames.append("%s%d_%d%s" % \
696 (pa, counter, f_counter, "F+3"))
697 self.aonames.append("%s%d_%d%s" % \
698 (pa, counter, f_counter, "F-3"))
699 elif basis.symmetries[j]=='g':
700 self.aonames.append("%s%d_%d%s" % \
701 (pa, counter, f_counter, "G 0"))
702 self.aonames.append("%s%d_%d%s" % \
703 (pa, counter, f_counter, "G+1"))
704 self.aonames.append("%s%d_%d%s" % \
705 (pa, counter, f_counter, "G-1"))
706 self.aonames.append("%s%d_%d%s" % \
707 (pa, counter, g_counter, "G+2"))
708 self.aonames.append("%s%d_%d%s" % \
709 (pa, counter, g_counter, "G-2"))
710 self.aonames.append("%s%d_%d%s" % \
711 (pa, counter, g_counter, "G+3"))
712 self.aonames.append("%s%d_%d%s" % \
713 (pa, counter, g_counter, "G-3"))
714 self.aonames.append("%s%d_%d%s" % \
715 (pa, counter, g_counter, "G+4"))
716 self.aonames.append("%s%d_%d%s" % \
717 (pa, counter, g_counter, "G-4"))
718 break
719 counter=counter+1
721 if line=="$closed shells\n":
722 line = inputfile.next()
723 temp = line.split()
724 occs = int(temp[1][2:])
725 self.homos = numpy.array([occs-1], "i")
727 if line == "$alpha shells\n":
728 line = inputfile.next()
729 temp = line.split()
730 occ_a = int(temp[1][2:])
731 line = inputfile.next() # should be $beta shells
732 line = inputfile.next() # the beta occs
733 temp = line.split()
734 occ_b = int(temp[1][2:])
735 self.homos = numpy.array([occ_a-1,occ_b-1], "i")
737 if line[12:24]=="OVERLAP(CAO)":
738 line = inputfile.next()
739 line = inputfile.next()
740 overlaparray=[]
741 self.aooverlaps=numpy.zeros( (self.nbasis, self.nbasis), "d")
742 while line != " ----------------------\n":
743 temp=line.split()
744 overlaparray.extend(map(float, temp))
745 line = inputfile.next()
746 counter=0
748 for i in range(0, self.nbasis, 1):
749 for j in range(0, i+1, 1):
750 self.aooverlaps[i][j]=overlaparray[counter]
751 self.aooverlaps[j][i]=overlaparray[counter]
752 counter=counter+1
754 if ( line[0:6] == "$scfmo" or line[0:12] == "$uhfmo_alpha" ) and line.find("scf") > 0:
755 temp = line.split()
757 if temp[1][0:7] == "scfdump":
758# self.logger.warning("SCF not converged?")
759 print("SCF not converged?!")
761 if line[0:12] == "$uhfmo_alpha": # if unrestricted, create flag saying so
762 unrestricted = 1
763 else:
764 unrestricted = 0
766 self.moenergies=[]
767 self.mocoeffs=[]
769 for spin in range(unrestricted + 1): # make sure we cover all instances
770 title = inputfile.next()
771 while(title[0] == "#"):
772 title = inputfile.next()
774# mocoeffs = numpy.zeros((self.nbasis, self.nbasis), "d")
775 moenergies = []
776 moarray=[]
778 if spin == 1 and title[0:11] == "$uhfmo_beta":
779 title = inputfile.next()
780 while title[0] == "#":
781 title = inputfile.next()
783 while(title[0] != '$'):
784 temp=title.split()
786 orb_symm=temp[1]
788 try:
789 energy = float(temp[2][11:].replace("D", "E"))
790 except ValueError:
791 print(spin, ": ", title)
793 orb_en = utils.convertor(energy,"hartree","eV")
795 moenergies.append(orb_en)
796 single_mo = []
798 while(len(single_mo)<self.nbasis):
799 self.updateprogress(inputfile, "Coefficients", self.cupdate)
800 title = inputfile.next()
801 lines_coeffs=self.split_molines(title)
802 single_mo.extend(lines_coeffs)
804 moarray.append(single_mo)
805 title = inputfile.next()
807# for i in range(0, len(moarray), 1):
808# for j in range(0, self.nbasis, 1):
809# try:
810# mocoeffs[i][j]=moarray[i][j]
811# except IndexError:
812# print "Index Error in mocoeffs.", spin, i, j
813# break
815 mocoeffs = numpy.array(moarray,"d")
816 self.mocoeffs.append(mocoeffs)
817 self.moenergies.append(moenergies)
819 if line[26:49] == "a o f o r c e - program":
820 self.vibirs = []
821 self.vibfreqs = []
822 self.vibsyms = []
823 self.vibdisps = []
825# while line[3:31] != "**** force : all done ****":
827 if line[12:26] == "ATOMIC WEIGHTS":
828#begin parsing atomic weights
829 self.vibmasses=[]
830 line=inputfile.next() # lines =======
831 line=inputfile.next() # notes
832 line=inputfile.next() # start reading
833 temp=line.split()
834 while(len(temp) > 0):
835 self.vibmasses.append(float(temp[2]))
836 line=inputfile.next()
837 temp=line.split()
839 if line[5:14] == "frequency":
840 if not hasattr(self,"vibfreqs"):
841 self.vibfreqs = []
842 self.vibfreqs = []
843 self.vibsyms = []
844 self.vibdisps = []
845 self.vibirs = []
847 temp=line.replace("i","-").split()
849 freqs = [utils.float(f) for f in temp[1:]]
850 self.vibfreqs.extend(freqs)
852 line=inputfile.next()
853 line=inputfile.next()
855 syms=line.split()
856 self.vibsyms.extend(syms[1:])
858 line=inputfile.next()
859 line=inputfile.next()
860 line=inputfile.next()
861 line=inputfile.next()
863 temp=line.split()
864 irs = [utils.float(f) for f in temp[2:]]
865 self.vibirs.extend(irs)
867 line=inputfile.next()
868 line=inputfile.next()
869 line=inputfile.next()
870 line=inputfile.next()
872 x=[]
873 y=[]
874 z=[]
876 line=inputfile.next()
877 while len(line) > 1:
878 temp=line.split()
879 x.append(map(float, temp[3:]))
881 line=inputfile.next()
882 temp=line.split()
883 y.append(map(float, temp[1:]))
885 line=inputfile.next()
886 temp=line.split()
887 z.append(map(float, temp[1:]))
888 line=inputfile.next()
890# build xyz vectors for each mode
892 for i in range(0, len(x[0]), 1):
893 disp=[]
894 for j in range(0, len(x), 1):
895 disp.append( [x[j][i], y[j][i], z[j][i]])
896 self.vibdisps.append(disp)
898# line=inputfile.next()
900 def after_parsing(self):
902 # delete all frequencies that correspond to translations or rotations
904 if hasattr(self,"vibfreqs"):
905 i = 0
906 while i < len(self.vibfreqs):
907 if self.vibfreqs[i]==0.0:
908 del self.vibfreqs[i]
909 del self.vibdisps[i]
910 del self.vibirs[i]
911 del self.vibsyms[i]
912 i -= 1
913 i += 1