Coverage for cclib/parser/adfparser.py : 98%
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) 2018, 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 ADF output files"""
10import itertools
11import re
13import numpy
15from cclib.parser import logfileparser
16from cclib.parser import utils
19class ADF(logfileparser.Logfile):
20 """An ADF log file"""
22 def __init__(self, *args, **kwargs):
24 # Call the __init__ method of the superclass
25 super(ADF, self).__init__(logname="ADF", *args, **kwargs)
27 def __str__(self):
28 """Return a string representation of the object."""
29 return "ADF log file %s" % (self.filename)
31 def __repr__(self):
32 """Return a representation of the object."""
33 return 'ADF("%s")' % (self.filename)
35 def normalisesym(self, label):
36 """Use standard symmetry labels instead of ADF labels.
38 To normalise:
39 (1) any periods are removed (except in the case of greek letters)
40 (2) XXX is replaced by X, and a " added.
41 (3) XX is replaced by X, and a ' added.
42 (4) The greek letters Sigma, Pi, Delta and Phi are replaced by
43 their lowercase equivalent.
44 """
45 greeks = ['Sigma', 'Pi', 'Delta', 'Phi']
46 for greek in greeks:
47 if label.startswith(greek):
48 return label.lower()
50 ans = label.replace(".", "")
51 if ans[1:3] == "''":
52 temp = ans[0] + '"'
53 ans = temp
55 l = len(ans)
56 if l > 1 and ans[0] == ans[1]: # Python only tests the second condition if the first is true
57 if l > 2 and ans[1] == ans[2]:
58 ans = ans.replace(ans[0]*3, ans[0]) + '"'
59 else:
60 ans = ans.replace(ans[0]*2, ans[0]) + "'"
61 return ans
63 def normalisedegenerates(self, label, num, ndict=None):
64 """Generate a string used for matching degenerate orbital labels
66 To normalise:
67 (1) if label is E or T, return label:num
68 (2) if label is P or D, look up in dict, and return answer
69 """
71 if not ndict:
72 ndict = {
73 'P': {0: "P:x", 1: "P:y", 2: "P:z"},
74 'D': {0: "D:z2", 1: "D:x2-y2", 2: "D:xy", 3: "D:xz", 4: "D:yz"}
75 }
77 if label in ndict:
78 if num in ndict[label]:
79 return ndict[label][num]
80 else:
81 return "%s:%i" % (label, num+1)
82 else:
83 return "%s:%i" % (label, num+1)
85 def before_parsing(self):
87 # Used to avoid extracting the final geometry twice in a GeoOpt
88 self.NOTFOUND, self.GETLAST, self.NOMORE = list(range(3))
89 self.finalgeometry = self.NOTFOUND
91 # Used for calculating the scftarget (variables names taken from the ADF manual)
92 self.accint = self.SCFconv = self.sconv2 = None
94 # keep track of nosym and unrestricted case to parse Energies since it doens't have an all Irreps section
95 self.nosymflag = False
96 self.unrestrictedflag = False
98 SCFCNV, SCFCNV2 = list(range(2)) # used to index self.scftargets[]
99 maxelem, norm = list(range(2)) # used to index scf.values
101 def extract(self, inputfile, line):
102 """Extract information from the file object inputfile."""
104 # If a file contains multiple calculations, currently we want to print a warning
105 # and skip to the end of the file, since cclib parses only the main system, which
106 # is usually the largest. Here we test this by checking if scftargets has already
107 # been parsed when another INPUT FILE segment is found, although this might
108 # not always be the best indicator.
109 if line.strip() == "(INPUT FILE)" and hasattr(self, "scftargets"):
110 self.logger.warning("Skipping remaining calculations")
111 inputfile.seek(0, 2)
112 return
114 # We also want to check to make sure we aren't parsing "Create" jobs,
115 # which normally come before the calculation we actually want to parse.
116 if line.strip() == "(INPUT FILE)":
117 while True:
118 self.updateprogress(inputfile, "Unsupported Information", self.fupdate)
119 line = next(inputfile) if line.strip() == "(INPUT FILE)" else None
120 if line and not line[:6] in ("Create", "create"):
121 break
122 line = next(inputfile)
124 version_searchstr = "Amsterdam Density Functional (ADF)"
125 if version_searchstr in line:
126 startidx = line.index(version_searchstr) + len(version_searchstr)
127 trimmed_line = line[startidx:].strip()[:-1]
128 # The package version is normally a year with revision number
129 # (such as 2013.01), but it may also be a random string (such as a
130 # version control branch name).
131 match = re.search(r"([\d\.]{4,7})", trimmed_line)
132 if match:
133 package_version = match.groups()[0]
134 # Use YYYY.MM as a short version.
135 self.metadata["legacy_package_version"] = package_version
136 else:
137 # This isn't as well-defined, but the field shouldn't be left
138 # empty. Grab whatever is there and parse it out in the
139 # following lines.
140 package_version = trimmed_line.strip()
141 # More detailed information can be found before "A D F", even if
142 # the above package version isn't numeric.
143 self.skip_line(inputfile, 's')
144 line = next(inputfile)
145 # Get the contents between the star border.
146 tokens = line.split()[1:-1]
147 assert len(tokens) >= 1
148 if tokens[0] == "Build":
149 package_version += "+{}".format(tokens[1])
150 else:
151 assert tokens[0][0] == "r"
152 # If a year-type version has already been parsed (YYYY(.nn)),
153 # it should take precedence, otherwise use the more detailed
154 # version first.
155 if match:
156 package_version = '{}dev{}'.format(package_version, tokens[0][1:])
157 else:
158 year = tokens[1].split("-")[0]
159 self.metadata["package_version_description"] = package_version
160 package_version = '{}dev{}'.format(year, tokens[0][1:])
161 self.metadata["legacy_package_version"] = year
162 self.metadata["package_version_date"] = tokens[1]
163 self.metadata["package_version"] = package_version
165 # In ADF 2014.01, there are (INPUT FILE) messages, so we need to use just
166 # the lines that start with 'Create' and run until the title or something
167 # else we are sure is is the calculation proper. It would be good to combine
168 # this with the previous block, if possible.
169 if line[:6] == "Create":
170 while line[:5] != "title" and "NO TITLE" not in line:
171 line = inputfile.next()
173 if line[1:10] == "Symmetry:":
174 info = line.split()
175 if info[1] == "NOSYM":
176 self.nosymflag = True
178 # Use this to read the subspecies of irreducible representations.
179 # It will be a list, with each element representing one irrep.
180 if line.strip() == "Irreducible Representations, including subspecies":
182 self.skip_line(inputfile, 'dashes')
184 self.irreps = []
185 line = next(inputfile)
186 while line.strip() != "":
187 self.irreps.append(line.split())
188 line = next(inputfile)
190 if line[4:13] == 'Molecule:':
191 info = line.split()
192 if info[1] == 'UNrestricted':
193 self.unrestrictedflag = True
195 if line[1:6] == "ATOMS":
196 # Find the number of atoms and their atomic numbers
197 # Also extract the starting coordinates (for a GeoOpt anyway)
198 # and the atommasses (previously called vibmasses)
199 self.updateprogress(inputfile, "Attributes", self.cupdate)
201 self.atomcoords = []
203 self.skip_lines(inputfile, ['header1', 'header2', 'header3'])
205 atomnos = []
206 atommasses = []
207 atomcoords = []
208 coreelectrons = []
209 line = next(inputfile)
210 while len(line) > 2: # ensure that we are reading no blank lines
211 info = line.split()
212 element = info[1].split('.')[0]
213 atomnos.append(self.table.number[element])
214 atomcoords.append(list(map(float, info[2:5])))
215 coreelectrons.append(int(float(info[5]) - float(info[6])))
216 atommasses.append(float(info[7]))
217 line = next(inputfile)
218 self.atomcoords.append(atomcoords)
220 self.set_attribute('natom', len(atomnos))
221 self.set_attribute('atomnos', atomnos)
222 self.set_attribute('atommasses', atommasses)
223 self.set_attribute('coreelectrons', coreelectrons)
225 if line[1:10] == "FRAGMENTS":
226 header = next(inputfile)
228 self.frags = []
229 self.fragnames = []
231 line = next(inputfile)
232 while len(line) > 2: # ensure that we are reading no blank lines
233 info = line.split()
235 if len(info) == 7: # fragment name is listed here
236 self.fragnames.append("%s_%s" % (info[1], info[0]))
237 self.frags.append([])
238 self.frags[-1].append(int(info[2]) - 1)
240 elif len(info) == 5: # add atoms into last fragment
241 self.frags[-1].append(int(info[0]) - 1)
243 line = next(inputfile)
245 # Extract charge
246 if line[1:11] == "Net Charge":
248 charge = int(line.split()[2])
249 self.set_attribute('charge', charge)
251 line = next(inputfile)
252 if len(line.strip()):
253 # Spin polar: 1 (Spin_A minus Spin_B electrons)
254 # (Not sure about this for higher multiplicities)
255 mult = int(line.split()[2]) + 1
256 else:
257 mult = 1
258 self.set_attribute('mult', mult)
260 if line[1:22] == "S C F U P D A T E S":
261 # find targets for SCF convergence
263 if not hasattr(self, "scftargets"):
264 self.scftargets = []
266 self.skip_lines(inputfile, ['e', 'b', 'numbers'])
268 line = next(inputfile)
269 self.SCFconv = float(line.split()[-1])
270 line = next(inputfile)
271 self.sconv2 = float(line.split()[-1])
273 # In ADF 2013, the default numerical integration method is fuzzy cells,
274 # although it used to be Voronoi polyhedra. Both methods apparently set
275 # the accint parameter, although the latter does so indirectly, based on
276 # a 'grid quality' setting. This is translated into accint using a
277 # dictionary with values taken from the documentation.
278 if "Numerical Integration : Voronoi Polyhedra (Te Velde)" in line:
279 self.integration_method = "voronoi_polyhedra"
280 if line[1:27] == 'General Accuracy Parameter':
281 # Need to know the accuracy of the integration grid to
282 # calculate the scftarget...note that it changes with time
283 self.accint = float(line.split()[-1])
284 if "Numerical Integration : Fuzzy Cells (Becke)" in line:
285 self.integration_method = 'fuzzy_cells'
286 if line[1:19] == "Becke grid quality":
287 self.grid_quality = line.split()[-1]
288 quality2accint = {
289 'BASIC': 2.0,
290 'NORMAL': 4.0,
291 'GOOD': 6.0,
292 'VERYGOOD': 8.0,
293 'EXCELLENT': 10.0,
294 }
295 self.accint = quality2accint[self.grid_quality]
297 # Half of the atomic orbital overlap matrix is printed since it is symmetric,
298 # but this requires "PRINT Smat" to be in the input. There are extra blank lines
299 # at the end of the block, which are used to terminate the parsing.
300 #
301 # ====== smat
302 #
303 # column 1 2 3 4
304 # row
305 # 1 1.00000000000000E+00
306 # 2 2.43370854175315E-01 1.00000000000000E+00
307 # 3 0.00000000000000E+00 0.00000000000000E+00 1.00000000000000E+00
308 # ...
309 #
310 if "====== smat" in line:
312 # Initialize the matrix with Nones so we can easily check all has been parsed.
313 overlaps = [[None] * self.nbasis for i in range(self.nbasis)]
315 self.skip_line(inputfile, 'blank')
317 line = inputfile.next()
318 while line.strip():
320 colline = line
321 assert colline.split()[0] == "column"
322 columns = [int(i) for i in colline.split()[1:]]
324 rowline = inputfile.next()
325 assert rowline.strip() == "row"
327 line = inputfile.next()
328 while line.strip():
330 i = int(line.split()[0])
331 vals = [float(col) for col in line.split()[1:]]
332 for j, o in enumerate(vals):
333 k = columns[j]
334 overlaps[k-1][i-1] = o
335 overlaps[i-1][k-1] = o
337 line = inputfile.next()
339 line = inputfile.next()
341 # Now all values should be parsed, and so no Nones remaining.
342 assert all([all([x is not None for x in ao]) for ao in overlaps])
344 self.set_attribute('aooverlaps', overlaps)
346 if line[1:11] == "CYCLE 1":
348 self.updateprogress(inputfile, "QM convergence", self.fupdate)
350 newlist = []
351 line = next(inputfile)
353 if not hasattr(self, "geovalues"):
354 # This is the first SCF cycle
355 self.scftargets.append([self.sconv2*10, self.sconv2])
356 elif self.finalgeometry in [self.GETLAST, self.NOMORE]:
357 # This is the final SCF cycle
358 self.scftargets.append([self.SCFconv*10, self.SCFconv])
359 else:
360 # This is an intermediate SCF cycle in a geometry optimization,
361 # in which case the SCF convergence target needs to be derived
362 # from the accint parameter. For Voronoi polyhedra integration,
363 # accint is printed and parsed. For fuzzy cells, it can be inferred
364 # from the grid quality setting, as is done somewhere above.
365 if self.accint:
366 oldscftst = self.scftargets[-1][1]
367 grdmax = self.geovalues[-1][1]
368 scftst = max(self.SCFconv, min(oldscftst, grdmax/30, 10**(-self.accint)))
369 self.scftargets.append([scftst*10, scftst])
371 while line.find("SCF CONVERGED") == -1 and line.find("SCF not fully converged, result acceptable") == -1 and line.find("SCF NOT CONVERGED") == -1:
372 if line[4:12] == "SCF test":
373 if not hasattr(self, "scfvalues"):
374 self.scfvalues = []
376 info = line.split()
377 newlist.append([float(info[4]), abs(float(info[6]))])
378 try:
379 line = next(inputfile)
380 except StopIteration: # EOF reached?
381 self.logger.warning("SCF did not converge, so attributes may be missing")
382 break
384 if line.find("SCF not fully converged, result acceptable") > 0:
385 self.logger.warning("SCF not fully converged, results acceptable")
387 if line.find("SCF NOT CONVERGED") > 0:
388 self.logger.warning("SCF did not converge! moenergies and mocoeffs are unreliable")
390 if hasattr(self, "scfvalues"):
391 self.scfvalues.append(newlist)
393 # Parse SCF energy for SP calcs from bonding energy decomposition section.
394 # It seems ADF does not print it earlier for SP calculations.
395 # Geometry optimization runs also print this, and we want to parse it
396 # for them, too, even if it repeats the last "Geometry Convergence Tests"
397 # section (but it's usually a bit different).
398 if line[:21] == "Total Bonding Energy:":
400 if not hasattr(self, "scfenergies"):
401 self.scfenergies = []
403 energy = utils.convertor(float(line.split()[3]), "hartree", "eV")
404 self.scfenergies.append(energy)
406 if line[51:65] == "Final Geometry":
407 self.finalgeometry = self.GETLAST
409 # Get the coordinates from each step of the GeoOpt.
410 if line[1:24] == "Coordinates (Cartesian)" and self.finalgeometry in [self.NOTFOUND, self.GETLAST]:
412 self.skip_lines(inputfile, ['e', 'b', 'title', 'title', 'd'])
414 atomcoords = []
415 line = next(inputfile)
416 while list(set(line.strip())) != ['-']:
417 atomcoords.append(list(map(float, line.split()[5:8])))
418 line = next(inputfile)
420 if not hasattr(self, "atomcoords"):
421 self.atomcoords = []
422 self.atomcoords.append(atomcoords)
424 # Don't get any more coordinates in this case.
425 # KML: I think we could combine this with optdone (see below).
426 if self.finalgeometry == self.GETLAST:
427 self.finalgeometry = self.NOMORE
429 # There have been some changes in the format of the geometry convergence information,
430 # and this is how it is printed in older versions (2007.01 unit tests).
431 #
432 # ==========================
433 # Geometry Convergence Tests
434 # ==========================
435 #
436 # Energy old : -5.14170647
437 # new : -5.15951374
438 #
439 # Convergence tests:
440 # (Energies in hartree, Gradients in hartree/angstr or radian, Lengths in angstrom, Angles in degrees)
441 #
442 # Item Value Criterion Conv. Ratio
443 # -------------------------------------------------------------------------
444 # change in energy -0.01780727 0.00100000 NO 0.00346330
445 # gradient max 0.03219530 0.01000000 NO 0.30402650
446 # gradient rms 0.00858685 0.00666667 NO 0.27221261
447 # cart. step max 0.07674971 0.01000000 NO 0.75559435
448 # cart. step rms 0.02132310 0.00666667 NO 0.55335378
449 #
450 if line[1:27] == 'Geometry Convergence Tests':
452 if not hasattr(self, "geotargets"):
453 self.geovalues = []
454 self.geotargets = numpy.array([0.0, 0.0, 0.0, 0.0, 0.0], "d")
456 if not hasattr(self, "scfenergies"):
457 self.scfenergies = []
459 self.skip_lines(inputfile, ['e', 'b'])
461 energies_old = next(inputfile)
462 energies_new = next(inputfile)
463 self.scfenergies.append(utils.convertor(float(energies_new.split()[-1]), "hartree", "eV"))
465 self.skip_lines(inputfile, ['b', 'convergence', 'units', 'b', 'header', 'd'])
467 values = []
468 for i in range(5):
469 temp = next(inputfile).split()
470 self.geotargets[i] = float(temp[-3])
471 values.append(float(temp[-4]))
473 self.geovalues.append(values)
475 # This is to make geometry optimization always have the optdone attribute,
476 # even if it is to be empty for unconverged runs.
477 if not hasattr(self, 'optdone'):
478 self.optdone = []
480 # After the test, there is a message if the search is converged:
481 #
482 # ***************************************************************************************************
483 # Geometry CONVERGED
484 # ***************************************************************************************************
485 #
486 if line.strip() == "Geometry CONVERGED":
487 self.skip_line(inputfile, 'stars')
488 self.optdone.append(len(self.geovalues) - 1)
490 # Here is the corresponding geometry convergence info from the 2013.01 unit test.
491 # Note that the step number is given, which it will be prudent to use in an assertion.
492 #
493 #----------------------------------------------------------------------
494 #Geometry Convergence after Step 3 (Hartree/Angstrom,Angstrom)
495 #----------------------------------------------------------------------
496 #current energy -5.16274478 Hartree
497 #energy change -0.00237544 0.00100000 F
498 #constrained gradient max 0.00884999 0.00100000 F
499 #constrained gradient rms 0.00249569 0.00066667 F
500 #gradient max 0.00884999
501 #gradient rms 0.00249569
502 #cart. step max 0.03331296 0.01000000 F
503 #cart. step rms 0.00844037 0.00666667 F
504 if line[:31] == "Geometry Convergence after Step":
506 stepno = int(line.split()[4])
508 # This is to make geometry optimization always have the optdone attribute,
509 # even if it is to be empty for unconverged runs.
510 if not hasattr(self, 'optdone'):
511 self.optdone = []
513 # The convergence message is inline in this block, not later as it was before.
514 if "** CONVERGED **" in line:
515 if not hasattr(self, 'optdone'):
516 self.optdone = []
517 self.optdone.append(len(self.geovalues) - 1)
519 self.skip_line(inputfile, 'dashes')
521 current_energy = next(inputfile)
522 energy_change = next(inputfile)
523 constrained_gradient_max = next(inputfile)
524 constrained_gradient_rms = next(inputfile)
525 gradient_max = next(inputfile)
526 gradient_rms = next(inputfile)
527 cart_step_max = next(inputfile)
528 cart_step_rms = next(inputfile)
530 if not hasattr(self, "scfenergies"):
531 self.scfenergies = []
533 energy = utils.convertor(float(current_energy.split()[-2]), "hartree", "eV")
534 self.scfenergies.append(energy)
536 if not hasattr(self, "geotargets"):
537 self.geotargets = numpy.array([0.0, 0.0, 0.0, 0.0, 0.0], "d")
539 self.geotargets[0] = float(energy_change.split()[-2])
540 self.geotargets[1] = float(constrained_gradient_max.split()[-2])
541 self.geotargets[2] = float(constrained_gradient_rms.split()[-2])
542 self.geotargets[3] = float(cart_step_max.split()[-2])
543 self.geotargets[4] = float(cart_step_rms.split()[-2])
545 if not hasattr(self, "geovalues"):
546 self.geovalues = []
548 self.geovalues.append([])
549 self.geovalues[-1].append(float(energy_change.split()[-3]))
550 self.geovalues[-1].append(float(constrained_gradient_max.split()[-3]))
551 self.geovalues[-1].append(float(constrained_gradient_rms.split()[-3]))
552 self.geovalues[-1].append(float(cart_step_max.split()[-3]))
553 self.geovalues[-1].append(float(cart_step_rms.split()[-3]))
555 if line.find('Orbital Energies, per Irrep and Spin') > 0 and not hasattr(self, "mosyms") and self.nosymflag and not self.unrestrictedflag:
556 #Extracting orbital symmetries and energies, homos for nosym case
557 #Should only be for restricted case because there is a better text block for unrestricted and nosym
559 self.mosyms = [[]]
561 self.moenergies = [[]]
563 self.skip_lines(inputfile, ['e', 'header', 'd', 'label'])
565 line = next(inputfile)
566 info = line.split()
568 if not info[0] == '1':
569 self.logger.warning("MO info up to #%s is missing" % info[0])
571 #handle case where MO information up to a certain orbital are missing
572 while int(info[0]) - 1 != len(self.moenergies[0]):
573 self.moenergies[0].append(99999)
574 self.mosyms[0].append('A')
576 homoA = None
578 while len(line) > 10:
579 info = line.split()
580 self.mosyms[0].append('A')
581 self.moenergies[0].append(utils.convertor(float(info[2]), 'hartree', 'eV'))
582 if info[1] == '0.000' and not hasattr(self, 'homos'):
583 self.set_attribute('homos', [len(self.moenergies[0]) - 2])
584 line = next(inputfile)
586 self.moenergies = [numpy.array(self.moenergies[0], "d")]
588 if line[1:29] == 'Orbital Energies, both Spins' and not hasattr(self, "mosyms") and self.nosymflag and self.unrestrictedflag:
589 #Extracting orbital symmetries and energies, homos for nosym case
590 #should only be here if unrestricted and nosym
592 self.mosyms = [[], []]
593 moenergies = [[], []]
595 self.skip_lines(inputfile, ['d', 'b', 'header', 'd'])
597 homoa = 0
598 homob = None
600 line = next(inputfile)
601 while len(line) > 5:
602 info = line.split()
603 if info[2] == 'A':
604 self.mosyms[0].append('A')
605 moenergies[0].append(utils.convertor(float(info[4]), 'hartree', 'eV'))
606 if info[3] != '0.00':
607 homoa = len(moenergies[0]) - 1
608 elif info[2] == 'B':
609 self.mosyms[1].append('A')
610 moenergies[1].append(utils.convertor(float(info[4]), 'hartree', 'eV'))
611 if info[3] != '0.00':
612 homob = len(moenergies[1]) - 1
613 else:
614 print(("Error reading line: %s" % line))
616 line = next(inputfile)
618 self.moenergies = [numpy.array(x, "d") for x in moenergies]
620 self.set_attribute('homos', [homoa, homob])
622 # Extracting orbital symmetries and energies, homos.
623 if line[1:29] == 'Orbital Energies, all Irreps' and not hasattr(self, "mosyms"):
625 self.symlist = {}
626 self.mosyms = [[]]
627 self.moenergies = [[]]
629 self.skip_lines(inputfile, ['e', 'b', 'header', 'd'])
631 homoa = None
632 homob = None
634 #multiple = {'E':2, 'T':3, 'P':3, 'D':5}
635 # The above is set if there are no special irreps
636 names = [irrep[0].split(':')[0] for irrep in self.irreps]
637 counts = [len(irrep) for irrep in self.irreps]
638 multiple = dict(list(zip(names, counts)))
639 irrepspecies = {}
640 for n in range(len(names)):
641 indices = list(range(counts[n]))
642 subspecies = self.irreps[n]
643 irrepspecies[names[n]] = dict(list(zip(indices, subspecies)))
645 line = next(inputfile)
646 while line.strip():
647 info = line.split()
648 if len(info) == 5: # this is restricted
649 #count = multiple.get(info[0][0],1)
650 count = multiple.get(info[0], 1)
651 for repeat in range(count): # i.e. add E's twice, T's thrice
652 self.mosyms[0].append(self.normalisesym(info[0]))
653 self.moenergies[0].append(utils.convertor(float(info[3]), 'hartree', 'eV'))
655 sym = info[0]
656 if count > 1: # add additional sym label
657 sym = self.normalisedegenerates(info[0], repeat, ndict=irrepspecies)
659 try:
660 self.symlist[sym][0].append(len(self.moenergies[0])-1)
661 except KeyError:
662 self.symlist[sym] = [[]]
663 self.symlist[sym][0].append(len(self.moenergies[0])-1)
665 if info[2] == '0.00' and not hasattr(self, 'homos'):
666 self.homos = [len(self.moenergies[0]) - (count + 1)] # count, because need to handle degenerate cases
667 line = next(inputfile)
668 elif len(info) == 6: # this is unrestricted
669 if len(self.moenergies) < 2: # if we don't have space, create it
670 self.moenergies.append([])
671 self.mosyms.append([])
672# count = multiple.get(info[0][0], 1)
673 count = multiple.get(info[0], 1)
674 if info[2] == 'A':
675 for repeat in range(count): # i.e. add E's twice, T's thrice
676 self.mosyms[0].append(self.normalisesym(info[0]))
677 self.moenergies[0].append(utils.convertor(float(info[4]), 'hartree', 'eV'))
679 sym = info[0]
680 if count > 1: # add additional sym label
681 sym = self.normalisedegenerates(info[0], repeat)
683 try:
684 self.symlist[sym][0].append(len(self.moenergies[0])-1)
685 except KeyError:
686 self.symlist[sym] = [[], []]
687 self.symlist[sym][0].append(len(self.moenergies[0])-1)
689 if info[3] == '0.00' and homoa is None:
690 homoa = len(self.moenergies[0]) - (count + 1) # count because degenerate cases need to be handled
692 if info[2] == 'B':
693 for repeat in range(count): # i.e. add E's twice, T's thrice
694 self.mosyms[1].append(self.normalisesym(info[0]))
695 self.moenergies[1].append(utils.convertor(float(info[4]), 'hartree', 'eV'))
697 sym = info[0]
698 if count > 1: # add additional sym label
699 sym = self.normalisedegenerates(info[0], repeat)
701 try:
702 self.symlist[sym][1].append(len(self.moenergies[1])-1)
703 except KeyError:
704 self.symlist[sym] = [[], []]
705 self.symlist[sym][1].append(len(self.moenergies[1])-1)
707 if info[3] == '0.00' and homob is None:
708 homob = len(self.moenergies[1]) - (count + 1)
710 line = next(inputfile)
712 else: # different number of lines
713 print(("Error", info))
715 if len(info) == 6: # still unrestricted, despite being out of loop
716 self.set_attribute('homos', [homoa, homob])
718 self.moenergies = [numpy.array(x, "d") for x in self.moenergies]
720 # Section on extracting vibdisps
721 # Also contains vibfreqs, but these are extracted in the
722 # following section (see below)
723 if line[1:28] == "Vibrations and Normal Modes":
725 self.vibdisps = []
727 self.skip_lines(inputfile, ['e', 'b', 'header', 'header', 'b', 'b'])
729 freqs = next(inputfile)
730 while freqs.strip() != "":
731 minus = next(inputfile)
732 p = [[], [], []]
733 for i in range(len(self.atomnos)):
734 broken = list(map(float, next(inputfile).split()[1:]))
735 for j in range(0, len(broken), 3):
736 p[j//3].append(broken[j:j+3])
737 self.vibdisps.extend(p[:(len(broken)//3)])
738 self.skip_lines(inputfile, ['b', 'b'])
739 freqs = next(inputfile)
740 self.vibdisps = numpy.array(self.vibdisps, "d")
742 if line[1:24] == "List of All Frequencies":
743 # Start of the IR/Raman frequency section
744 self.updateprogress(inputfile, "Frequency information", self.fupdate)
746 # self.vibsyms = [] # Need to look into this a bit more
747 self.vibirs = []
748 self.vibfreqs = []
749 for i in range(8):
750 line = next(inputfile)
751 line = next(inputfile).strip()
752 while line:
753 temp = line.split()
754 self.vibfreqs.append(float(temp[0]))
755 self.vibirs.append(float(temp[2])) # or is it temp[1]?
756 line = next(inputfile).strip()
757 self.vibfreqs = numpy.array(self.vibfreqs, "d")
758 self.vibirs = numpy.array(self.vibirs, "d")
759 if hasattr(self, "vibramans"):
760 self.vibramans = numpy.array(self.vibramans, "d")
762 #******************************************************************************************************************8
763 #delete this after new implementation using smat, eigvec print,eprint?
764 # Extract the number of basis sets
765 if line[1:49] == "Total nr. of (C)SFOs (summation over all irreps)":
766 nbasis = int(line.split(":")[1].split()[0])
767 self.set_attribute('nbasis', nbasis)
769 # now that we're here, let's extract aonames
771 self.fonames = []
772 self.start_indeces = {}
773 self.atombasis = [[] for frag in self.frags] # parse atombasis in the case of trivial SFOs
775 self.skip_line(inputfile, 'blank')
777 note = next(inputfile)
778 symoffset = 0
780 self.skip_line(inputfile, 'blank')
781 line = next(inputfile)
782 if len(line) > 2: # fix for ADF2006.01 as it has another note
783 self.skip_line(inputfile, 'blank')
784 line = next(inputfile)
785 self.skip_line(inputfile, 'blank')
787 self.nosymreps = []
788 while len(self.fonames) < self.nbasis:
790 symline = next(inputfile)
791 sym = symline.split()[1]
792 line = next(inputfile)
793 num = int(line.split(':')[1].split()[0])
794 self.nosymreps.append(num)
796 #read until line "--------..." is found
797 while line.find('-----') < 0:
798 line = next(inputfile)
800 line = next(inputfile) # the start of the first SFO
802 while len(self.fonames) < symoffset + num:
803 info = line.split()
805 #index0 index1 occ2 energy3/4 fragname5 coeff6 orbnum7 orbname8 fragname9
806 if not sym in list(self.start_indeces.keys()):
807 #have we already set the start index for this symmetry?
808 self.start_indeces[sym] = int(info[1])
810 orbname = info[8]
811 orbital = info[7] + orbname.replace(":", "")
813 fragname = info[5]
814 frag = fragname + info[9]
816 coeff = float(info[6])
818 # parse atombasis only in the case that all coefficients are 1
819 # and delete it otherwise
820 if hasattr(self, 'atombasis'):
821 if coeff == 1.:
822 ibas = int(info[0]) - 1
823 ifrag = int(info[9]) - 1
824 iat = self.frags[ifrag][0]
825 self.atombasis[iat].append(ibas)
826 else:
827 del self.atombasis
829 line = next(inputfile)
830 while line.strip() and not line[:7].strip(): # while it's the same SFO
831 # i.e. while not completely blank, but blank at the start
832 info = line[43:].split()
833 if len(info) > 0: # len(info)==0 for the second line of dvb_ir.adfout
834 frag += "+" + fragname + info[-1]
835 coeff = float(info[-4])
836 if coeff < 0:
837 orbital += '-' + info[-3] + info[-2].replace(":", "")
838 else:
839 orbital += '+' + info[-3] + info[-2].replace(":", "")
840 line = next(inputfile)
841 # At this point, we are either at the start of the next SFO or at
842 # a blank line...the end
844 self.fonames.append("%s_%s" % (frag, orbital))
845 symoffset += num
847 # blankline blankline
848 next(inputfile)
849 next(inputfile)
851 if line[1:32] == "S F O P O P U L A T I O N S ,":
852 #Extract overlap matrix
854# self.fooverlaps = numpy.zeros((self.nbasis, self.nbasis), "d")
856 symoffset = 0
858 for nosymrep in self.nosymreps:
860 line = next(inputfile)
861 while line.find('===') < 10: # look for the symmetry labels
862 line = next(inputfile)
864 self.skip_lines(inputfile, ['b', 'b'])
866 text = next(inputfile)
867 if text[13:20] != "Overlap": # verify this has overlap info
868 break
870 self.skip_lines(inputfile, ['b', 'col', 'row'])
872 if not hasattr(self, "fooverlaps"): # make sure there is a matrix to store this
873 self.fooverlaps = numpy.zeros((self.nbasis, self.nbasis), "d")
875 base = 0
876 while base < nosymrep: # have we read all the columns?
878 for i in range(nosymrep - base):
880 self.updateprogress(inputfile, "Overlap", self.fupdate)
881 line = next(inputfile)
882 parts = line.split()[1:]
883 for j in range(len(parts)):
884 k = float(parts[j])
885 self.fooverlaps[base + symoffset + j, base + symoffset + i] = k
886 self.fooverlaps[base + symoffset + i, base + symoffset + j] = k
888 #blank, blank, column
889 for i in range(3):
890 next(inputfile)
892 base += 4
894 symoffset += nosymrep
895 base = 0
897# The commented code below makes the atombasis attribute based on the BAS function in ADF,
898# but this is probably not so useful, since SFOs are used to build MOs in ADF.
899# if line[1:54] == "BAS: List of all Elementary Cartesian Basis Functions":
900#
901# self.atombasis = []
902#
903# # There will be some text, followed by a line:
904# # (power of) X Y Z R Alpha on Atom
905# while not line[1:11] == "(power of)":
906# line = inputfile.next()
907# dashes = inputfile.next()
908# blank = inputfile.next()
909# line = inputfile.next()
910# # There will be two blank lines when there are no more atom types.
911# while line.strip() != "":
912# atoms = [int(i)-1 for i in line.split()[1:]]
913# for n in range(len(atoms)):
914# self.atombasis.append([])
915# dashes = inputfile.next()
916# line = inputfile.next()
917# while line.strip() != "":
918# indices = [int(i)-1 for i in line.split()[5:]]
919# for i in range(len(indices)):
920# self.atombasis[atoms[i]].append(indices[i])
921# line = inputfile.next()
922# line = inputfile.next()
924 if line[48:67] == "SFO MO coefficients":
926 self.mocoeffs = [numpy.zeros((self.nbasis, self.nbasis), "d")]
927 spin = 0
928 symoffset = 0
929 lastrow = 0
931 # Section ends with "1" at beggining of a line.
932 while line[0] != "1":
933 line = next(inputfile)
935 # If spin is specified, then there will be two coefficient matrices.
936 if line.strip() == "***** SPIN 1 *****":
937 self.mocoeffs = [numpy.zeros((self.nbasis, self.nbasis), "d"),
938 numpy.zeros((self.nbasis, self.nbasis), "d")]
940 # Bump up the spin.
941 if line.strip() == "***** SPIN 2 *****":
942 spin = 1
943 symoffset = 0
944 lastrow = 0
946 # Next symmetry.
947 if line.strip()[:4] == "=== ":
948 sym = line.split()[1]
949 if self.nosymflag:
950 aolist = list(range(self.nbasis))
951 else:
952 aolist = self.symlist[sym][spin]
953 # Add to the symmetry offset of AO ordering.
954 symoffset += lastrow
956 # Blocks with coefficient always start with "MOs :".
957 if line[1:6] == "MOs :":
958 # Next line has the MO index contributed to.
959 monumbers = [int(n) for n in line[6:].split()]
961 self.skip_lines(inputfile, ['occup', 'label'])
963 # The table can end with a blank line or "1".
964 row = 0
965 line = next(inputfile)
966 while not line.strip() in ["", "1"]:
967 info = line.split()
969 if int(info[0]) < self.start_indeces[sym]:
970 #check to make sure we aren't parsing CFs
971 line = next(inputfile)
972 continue
974 self.updateprogress(inputfile, "Coefficients", self.fupdate)
975 row += 1
976 coeffs = [float(x) for x in info[1:]]
977 moindices = [aolist[n-1] for n in monumbers]
978 # The AO index is 1 less than the row.
979 aoindex = symoffset + row - 1
980 for i in range(len(monumbers)):
981 self.mocoeffs[spin][moindices[i], aoindex] = coeffs[i]
982 line = next(inputfile)
983 lastrow = row
985 # **************************************************************************
986 # * *
987 # * Final excitation energies from Davidson algorithm *
988 # * *
989 # **************************************************************************
990 #
991 # Number of loops in Davidson routine = 20
992 # Number of matrix-vector multiplications = 24
993 # Type of excitations = SINGLET-SINGLET
994 #
995 # Symmetry B.u
996 #
997 # ... several blocks ...
998 #
999 # Normal termination of EXCITATION program part
1000 if line[4:53] == "Final excitation energies from Davidson algorithm":
1002 while line[1:9] != "Symmetry" and "Normal termination" not in line:
1003 line = next(inputfile)
1004 symm = self.normalisesym(line.split()[1])
1006 # Excitation energies E in a.u. and eV, dE wrt prev. cycle,
1007 # oscillator strengths f in a.u.
1008 #
1009 # no. E/a.u. E/eV f dE/a.u.
1010 # -----------------------------------------------------
1011 # 1 0.17084 4.6488 0.16526E-01 0.28E-08
1012 # ...
1013 while line.split() != ['no.', 'E/a.u.', 'E/eV', 'f', 'dE/a.u.'] and "Normal termination" not in line:
1014 line = next(inputfile)
1016 self.skip_line(inputfile, 'dashes')
1018 etenergies = []
1019 etoscs = []
1020 etsyms = []
1021 line = next(inputfile)
1022 while len(line) > 2:
1023 info = line.split()
1024 etenergies.append(utils.convertor(float(info[2]), "eV", "wavenumber"))
1025 etoscs.append(float(info[3]))
1026 etsyms.append(symm)
1027 line = next(inputfile)
1029 # There is another section before this, with transition dipole moments,
1030 # but this should just skip past it.
1031 while line[1:53] != "Major MO -> MO transitions for the above excitations":
1032 line = next(inputfile)
1034 # Note that here, and later, the number of blank lines can vary between
1035 # version of ADF (extra lines are seen in 2013.01 unit tests, for example).
1036 self.skip_line(inputfile, 'blank')
1037 excitation_occupied = next(inputfile)
1038 header = next(inputfile)
1039 while not header.strip():
1040 header = next(inputfile)
1041 header2 = next(inputfile)
1042 x_y_z = next(inputfile)
1043 line = next(inputfile)
1044 while not line.strip():
1045 line = next(inputfile)
1047 # Before we start handeling transitions, we need to create mosyms
1048 # with indices; only restricted calcs are possible in ADF.
1049 counts = {}
1050 syms = []
1051 for mosym in self.mosyms[0]:
1052 if list(counts.keys()).count(mosym) == 0:
1053 counts[mosym] = 1
1054 else:
1055 counts[mosym] += 1
1056 syms.append(str(counts[mosym]) + mosym)
1058 etsecs = []
1059 printed_warning = False
1060 for i in range(len(etenergies)):
1062 etsec = []
1063 info = line.split()
1064 while len(info) > 0:
1066 match = re.search('[^0-9]', info[1])
1067 index1 = int(info[1][:match.start(0)])
1068 text = info[1][match.start(0):]
1069 symtext = text[0].upper() + text[1:]
1070 sym1 = str(index1) + self.normalisesym(symtext)
1072 match = re.search('[^0-9]', info[3])
1073 index2 = int(info[3][:match.start(0)])
1074 text = info[3][match.start(0):]
1075 symtext = text[0].upper() + text[1:]
1076 sym2 = str(index2) + self.normalisesym(symtext)
1078 try:
1079 index1 = syms.index(sym1)
1080 except ValueError:
1081 if not printed_warning:
1082 self.logger.warning("Etsecs are not accurate!")
1083 printed_warning = True
1085 try:
1086 index2 = syms.index(sym2)
1087 except ValueError:
1088 if not printed_warning:
1089 self.logger.warning("Etsecs are not accurate!")
1090 printed_warning = True
1092 etsec.append([(index1, 0), (index2, 0), float(info[4])])
1094 line = next(inputfile)
1095 info = line.split()
1097 etsecs.append(etsec)
1099 # Again, the number of blank lines between transition can vary.
1100 line = next(inputfile)
1101 while not line.strip():
1102 line = next(inputfile)
1104 if not hasattr(self, "etenergies"):
1105 self.etenergies = etenergies
1106 else:
1107 self.etenergies += etenergies
1109 if not hasattr(self, "etoscs"):
1110 self.etoscs = etoscs
1111 else:
1112 self.etoscs += etoscs
1114 if not hasattr(self, "etsyms"):
1115 self.etsyms = etsyms
1116 else:
1117 self.etsyms += etsyms
1119 if not hasattr(self, "etsecs"):
1120 self.etsecs = etsecs
1121 else:
1122 self.etsecs += etsecs
1124 if "M U L L I K E N P O P U L A T I O N S" in line:
1125 if not hasattr(self, "atomcharges"):
1126 self.atomcharges = {}
1127 while line[1:5] != "Atom":
1128 line = next(inputfile)
1129 self.skip_line(inputfile, 'dashes')
1130 mulliken = []
1131 line = next(inputfile)
1132 while line.strip():
1133 mulliken.append(float(line.split()[2]))
1134 line = next(inputfile)
1135 self.atomcharges["mulliken"] = mulliken
1137 # Dipole moment is always printed after a point calculation,
1138 # and the reference point for this is always the origin (0,0,0)
1139 # and not necessarily the center of mass, as explained on the
1140 # ADF user mailing list (see cclib/cclib#113 for details).
1141 #
1142 # =============
1143 # Dipole Moment *** (Debye) ***
1144 # =============
1145 #
1146 # Vector : 0.00000000 0.00000000 0.00000000
1147 # Magnitude: 0.00000000
1148 #
1149 if line.strip()[:13] == "Dipole Moment":
1151 self.skip_line(inputfile, 'equals')
1153 # There is not always a blank line here, for example when the dipole and quadrupole
1154 # moments are printed after the multipole derived atomic charges. Still, to the best
1155 # of my knowledge (KML) the values are still in Debye.
1156 line = next(inputfile)
1157 if not line.strip():
1158 line = next(inputfile)
1160 assert line.split()[0] == "Vector"
1161 dipole = [float(d) for d in line.split()[-3:]]
1163 reference = [0.0, 0.0, 0.0]
1164 if not hasattr(self, 'moments'):
1165 self.moments = [reference, dipole]
1166 else:
1167 try:
1168 assert self.moments[1] == dipole
1169 except AssertionError:
1170 self.logger.warning('Overwriting previous multipole moments with new values')
1171 self.moments = [reference, dipole]
1173 # Molecular response properties.
1174 if line.strip()[1:-1].strip() == "RESPONSE program part":
1176 while line.strip() != "Normal termination of RESPONSE program part":
1178 if "THE DIPOLE-DIPOLE POLARIZABILITY TENSOR:" in line:
1179 if not hasattr(self, 'polarizabilities'):
1180 self.polarizabilities = []
1181 polarizability = numpy.empty(shape=(3, 3))
1182 self.skip_lines(inputfile, ['b', 'FREQUENCY', 'coordinates'])
1183 # Ordering of rows/columns is Y, Z, X.
1184 ordering = [1, 2, 0]
1185 indices = list(itertools.product(ordering, ordering))
1186 for i in range(3):
1187 tokens = next(inputfile).split()
1188 for j in range(3):
1189 polarizability[indices[(i*3)+j]] = tokens[j]
1190 self.polarizabilities.append(polarizability)
1192 line = next(inputfile)
1194 if line[:24] == ' Buffered I/O statistics':
1195 self.metadata['success'] = True