/
formchk_interface.py
87 lines (75 loc) · 3.18 KB
/
formchk_interface.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
import numpy as np
class FormchkInterface:
def __init__(self, file_path):
self.file_path = file_path
self.natm = NotImplemented
self.nao = NotImplemented
self.nmo = NotImplemented
self.initialization()
def initialization(self):
self.natm = int(self.key_to_value("Number of atoms"))
self.nao = int(self.key_to_value("Number of basis functions"))
self.nmo = int(self.key_to_value("Number of independent functions"))
def key_to_value(self, key, file_path=None):
if file_path is None:
file_path = self.file_path
flag_read = False
expect_size = -1
vec = []
with open(file_path, "r") as file:
for l in file:
if l[:len(key)] == key:
try:
expect_size = int(l[len(key):].split()[2])
flag_read = True
continue
except IndexError:
try:
return float(l[len(key):].split()[1])
except IndexError:
continue
if flag_read:
try:
vec += [float(i) for i in l.split()]
except ValueError:
break
if len(vec) != expect_size:
raise ValueError("Number of expected size is not consistent with read-in size!")
return np.array(vec)
def total_energy(self, file_path=None):
if file_path is None:
file_path = self.file_path
return self.key_to_value("Total Energy", file_path)
def grad(self, file_path=None):
if file_path is None:
file_path = self.file_path
return self.key_to_value("Cartesian Gradient", file_path).reshape((self.natm, 3))
def dipole(self, file_path=None):
if file_path is None:
file_path = self.file_path
return self.key_to_value("Dipole Moment", file_path)
@staticmethod
def tril_to_symm(tril: np.ndarray):
dim = int(np.floor(np.sqrt(tril.size * 2)))
if dim * (dim + 1) / 2 != tril.size:
raise ValueError("Size " + str(tril.size) + " is probably not a valid lower-triangle matrix.")
indices_tuple = np.tril_indices(dim)
iterator = zip(*indices_tuple)
symm = np.empty((dim, dim))
for it, (row, col) in enumerate(iterator):
symm[row, col] = tril[it]
symm[col, row] = tril[it]
return symm
def hessian(self, file_path=None):
if file_path is None:
file_path = self.file_path
return self.tril_to_symm(self.key_to_value("Cartesian Force Constants", file_path))
def polarizability(self, file_path=None):
if file_path is None:
file_path = self.file_path
# two space after `Polarizability' is to avoid `Polarizability Derivative'
return self.tril_to_symm(self.key_to_value("Polarizability ", file_path))
def dipolederiv(self, file_path=None):
if file_path is None:
file_path = self.file_path
return self.key_to_value("Dipole Derivatives", file_path).reshape(-1, 3)