Add pseudopotentials from Gaussian output

This commit is contained in:
Anthony Scemama 2021-12-20 18:16:13 +01:00
parent 427514e36e
commit 2d33329c1c
1 changed files with 69 additions and 26 deletions

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
# resultsFile is a library which allows to read output files of quantum
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@ -18,11 +18,11 @@
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
@ -35,7 +35,7 @@ import re
gaussianFile_defined_vars = [ "date", "version", "machine", "memory", "disk",\
"cpu_time", "author", "title", "units", "methods", "options", \
"spin_restrict", "conv_threshs", "energies", \
"ee_pot_energies", \
"ee_pot_energies", "pseudo", \
"Ne_pot_energies", "pot_energies", \
"kin_energies", "point_group", "num_elec", \
"charge", "multiplicity","nuclear_energy","dipole","geometry",\
@ -103,7 +103,7 @@ class gaussianFile(resultsFile.resultsFileX):
for x in dip:
self._dipole.append(float(x))
return self._dipole
def get_multiplicity(self):
if self._multiplicity is None:
self._multiplicity = int(self.options[3][0].split(',')[1])
@ -229,7 +229,7 @@ class gaussianFile(resultsFile.resultsFileX):
def get_spin_restrict(self):
if self._spin_restrict is None:
method = self.methods[0]
self._spin_restrict = True
self._spin_restrict = True
if method.startswith("U"): self._spin_restrict = False
return self._spin_restrict
@ -266,7 +266,7 @@ class gaussianFile(resultsFile.resultsFileX):
def get_ee_pot_energies(self):
if self._ee_pot_energies is None:
self._ee_pot_energies = []
self._ee_pot_energies = []
for i,e in enumerate(self.kin_energies):
self._ee_pot_energies.append(self.energies[i]\
-self.nuclear_energy\
@ -299,7 +299,7 @@ class gaussianFile(resultsFile.resultsFileX):
def get_point_group(self):
if self._point_group is None:
self._point_group = self.options[4][-1].split()[0].split('=')[1]
return self._point_group
return self._point_group
def get_geometry(self):
if self._geometry is None:
@ -323,7 +323,7 @@ class gaussianFile(resultsFile.resultsFileX):
temp.charge = float(buffer[1])
temp.coord = (float(buffer[3]), float(buffer[4]), float(buffer[5]))
self._geometry.append(temp)
pos += 1
pos += 1
for i,line in enumerate(self.options[3][1:]):
buffer = line.split(',')
self._geometry[i].name = buffer[0]
@ -454,21 +454,21 @@ class gaussianFile(resultsFile.resultsFileX):
elif basis[k][0] == "G":
if not Polar:
mylist = [ "ZZZZ", "ZZZY", "YYZZ", "YYYZ", "YYYY", "ZZZX",
"ZZXY", "YYXZ","YYYX", "XXZZ", "XXYZ",
"ZZXY", "YYXZ","YYYX", "XXZZ", "XXYZ",
"XXYY", "XXXZ","XXXY", "XXXX" ]
else:
mylist = [ "G0", "G1+", "G1-", "G2+", "G2-", "G3+", "G3-",
"G4+", "G4-" ]
elif basis[k][0] == "H":
if not Polar:
mylist = [ "ZZZZZ", "YZZZZ", "YYZZZ", "YYYZZ", "YYYYZ", "YYYYY",
mylist = [ "ZZZZZ", "YZZZZ", "YYZZZ", "YYYZZ", "YYYYZ", "YYYYY",
"XZZZZ", "XYZZZ", "XYYZZ", "XYYYZ", "XYYYY", "XXZZZ", "XXYZZ", "XXYYZ",
"XXYYY", "XXXZZ", "XXXYZ", "XXXYY", "XXXXZ", "XXXXY", "XXXXX" ]
"XXYYY", "XXXZZ", "XXXYZ", "XXXYY", "XXXXZ", "XXXXY", "XXXXX" ]
else:
mylist = [ "H0", "H1+", "H1-", "H2+", "H2-", "H3+", "H3-",
"H4+", "H4-", "H5+", "H5-" ]
elif basis[k][0] == "I":
if not Polar:
if not Polar:
mylist = [ "ZZZZZZ", "YZZZZZ", "YYZZZZ", "YYYZZZ",
"YYYYZZ", "YYYYYZ", "YYYYYY", "XZZZZZ", "XYZZZZ", "XYYZZZ",
"XYYYZZ", "XYYYYZ", "XYYYYY", "XXZZZZ", "XXYZZZ", "XXYYZZ",
@ -510,7 +510,7 @@ class gaussianFile(resultsFile.resultsFileX):
def get_mo_types(self):
if self._mo_types is None:
self.get_mo_sets()
return self._mo_types
return self._mo_types
def get_mo_sets(self):
if self._mo_sets is None:
@ -585,7 +585,7 @@ class gaussianFile(resultsFile.resultsFileX):
line = line.replace('**********','-99999999.')
line = line.split()
if line[0].lower() == "unable":
pos +=1
pos +=1
line = line.split()
eigval = [ float(k) for k in line ]
pos += 1
@ -616,7 +616,7 @@ class gaussianFile(resultsFile.resultsFileX):
if self._mo_sets == {}:
self._mo_sets = None
return self._mo_sets
def get_mulliken_mo(self):
if self._mulliken_mo is None:
pass
@ -632,6 +632,49 @@ class gaussianFile(resultsFile.resultsFileX):
pass
return self._lowdin_ao
def get_pseudo(self):
if self._pseudo is None:
try:
self.find_string("Pseudopotential Parameters")
except IndexError:
return None
self._pos += 5
pos = self._pos
try:
self.find_next_string("=======")
except IndexError:
return None
end = self._pos
pseudo_read = []
text = self.text
bdict = {"S": 0, "P": 1, "D": 2, "F": 3, "G": 4, "H": 5}
while pos < end:
ecp = {}
center, _, valence = [int(x) for x in text[pos].split()]
ecp["atom"] = center
ecp["zcore"] = valence
ecp["lmax"] = 0
pos += 1
l = 0
while text[pos][:35].strip() == "":
# Read angular momentum
line = text[pos].lstrip()
buffer = line.split()
if buffer[0] in "SPDFGH":
contraction = []
l = bdict[buffer[0]]
ecp["lmax"] = max(ecp["lmax"], l)
else:
# Read gaussian parameters
power, expo, coef = int(buffer[0]), float(buffer[1]), float(buffer[2])
contraction.append( (coef, power, expo) )
ecp[str(l)] = contraction
pos += 1
pseudo_read.append(ecp)
self._pseudo = pseudo_read
return self._pseudo
def get_mulliken_atom(self):
if self._mulliken_atom is None:
@ -671,14 +714,14 @@ class gaussianFile(resultsFile.resultsFileX):
self.find_string("alpha electrons")
pos = self._pos
self._num_alpha = int(self.text[pos].split()[0])
return self._num_alpha
return self._num_alpha
def get_num_beta(self):
if self._num_beta is None:
self.find_string("beta electrons")
pos = self._pos
self._num_beta = int(self.text[pos].split()[3])
return self._num_beta
return self._num_beta
def get_determinants_mo_type(self):
if self._determinants_mo_type is None:
@ -890,7 +933,7 @@ class gaussianFile(resultsFile.resultsFileX):
for motype in self.mo_types:
occ[motype] = [ mo.eigenvalue for mo in self.mulliken_mo ]
self._occ_num = occ
return self._occ_num
return self._occ_num
def get_num_states(self):
if self._num_states is None:
@ -923,7 +966,7 @@ class gaussianFile(resultsFile.resultsFileX):
self._file = fortranBinary(filename,"rb")
def get_form(self):
return self._file.form
return self._file.form
def set_form(self,form):
self._file.form = form
@ -988,7 +1031,7 @@ def write_vec(res, file):
monum = monum[-2:]
fields.append ( monum )
linum = "%4d"%(line)
linum = linum[-3:]
linum = linum[-3:]
fields.append ( linum )
for n in range(5):
fields.append ( "%15.8E"%v.pop(0) )
@ -1041,8 +1084,8 @@ def write_data(res, file):
print(" $END", file=file)
resultsFile.fileTypes.append(gaussianFile)
if __name__ == '__main__':
resultsFile.main(gaussianFile)