forked from Masterluke87/psixas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
spec.py
77 lines (52 loc) · 2.62 KB
/
spec.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
import numpy as np
import pdb
import psi4
import pickle
def CalcSpec(mol,func):
psi4.core.print_out("\n\nX-Ray Absorption Spectrum Calculation:\n"+38*"="+"\n\n")
prefix = psi4.core.get_local_option("PSIXAS","PREFIX")
psi4.core.print_out("Using orbitals, occupations from file: {} \n".format(prefix+"_exorbs.npz"))
Ca = np.load(prefix+"_exorbs.npz")["Ca"]
Cb = np.load(prefix+"_exorbs.npz")["Cb"]
occa = np.load(prefix+"_exorbs.npz")["occa"]
occb = np.load(prefix+"_exorbs.npz")["occb"]
epsa = np.load(prefix+"_exorbs.npz")["epsa"]
epsb = np.load(prefix+"_exorbs.npz")["epsb"]
orbitals = np.load(prefix+"_exorbs.npz", allow_pickle=True)["orbitals"]
psi4.core.print_out("Occupation pattern: \n")
printOccupation("Alpha",occa,15)
printOccupation("Beta ",occb,15)
#psi4.core.be_quiet()
wfn = psi4.core.Wavefunction.build(mol,psi4.core.get_global_option('BASIS'))
mints = psi4.core.MintsHelper(wfn.basisset())
psi4.core.reopen_outfile()
D = mints.ao_dipole()
Dx,Dy,Dz = np.asarray(D[0]),np.asarray(D[1]),np.asarray(D[2])
spec = {}
if ("b" in [x['spin'] for x in orbitals]):
psi4.core.print_out("\nBETA orbitals")
orbI = ([c for c,x in enumerate(occb) if x != 1.0][0])
psi4.core.print_out("\nInitial Orbital: {:d}".format(orbI))
Ci = Cb[orbI]
orbF = ([c for c,x in enumerate(occb) if (x != 1.0) and (c!=orbI)])
psi4.core.print_out("\nFinal Orbitals: {} \n\n".format(str(orbF)))
Cf = Cb[orbF]
Mx = Cb.T @ Dx @ Cb
My = Cb.T @ Dy @ Cb
Mz = Cb.T @ Dz @ Cb
spec["En"] = epsb[orbF] - epsb[orbI]
spec["Dx"] = Mx[orbI,orbF]
spec["Dy"] = My[orbI,orbF]
spec["Dz"] = Mz[orbI,orbF]
with open(prefix+'_b.spectrum', 'wb') as handle:
pickle.dump(spec, handle, protocol=pickle.HIGHEST_PROTOCOL)
psi4.core.print_out(("\n{}"+"_b.spectrum written.. \n\n").format(prefix))
if ("a" in [x['spin'] for x in orbitals]):
print ("ALFA orbitals")
def printOccupation(title,occs,width):
psi4.core.print_out("\n{}: \n".format(title))
for i in range(int(len(occs)/width)):
psi4.core.print_out(("|"+"{:^4}|"*width).format(*[x for x in range(i*width,i*width+ width)]))
psi4.core.print_out(("\n|"+"{:3.2f}|"*width+"\n\n").format(*[x for x in occs[i*width:i*width+width]]))
psi4.core.print_out(("|"+"{:^4}|"*(len(occs) % width)).format(*[x for x in range(int(len(occs)/width)*width,len(occs))]))
psi4.core.print_out(("\n|"+"{:3.2f}|"*(len(occs) % width)+"\n\n").format(*[x for x in occs[int(len(occs)/width)*width:]]))