-
Notifications
You must be signed in to change notification settings - Fork 7
/
local_symmetry.py
169 lines (136 loc) · 6.45 KB
/
local_symmetry.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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
# Application of local symmetry analysis of RASCI WF for ethylene dimer
# This example requires posym (https://github.com/abelcarreras/posym)
from pyqchem import Structure, QchemInput, get_output_from_qchem
from pyqchem.parsers.parser_optimization import basic_optimization
from pyqchem.parsers.parser_rasci import parser_rasci
from pyqchem.symmetry import get_state_symmetry
from pyqchem.file_io import write_to_fchk
from pyqchem.utils import get_plane, crop_electronic_structure
from pyqchem.symmetry import get_wf_symmetry
from posym import SymmetryBase
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# define monomer
coor_monomer = [[ 0.6695, 0.0000, 0.0000],
[-0.6695, 0.0000, 0.0000],
[ 1.2321, 0.9289, 0.0000],
[ 1.2321, -0.9289, 0.0000],
[-1.2321, 0.9289, 0.0000],
[-1.2321, -0.9289, 0.0000]]
symbols_monomer = ['C', 'C', 'H', 'H', 'H', 'H']
monomer = Structure(coordinates=coor_monomer,
symbols=symbols_monomer,
charge=0,
multiplicity=1)
# optimization qchem input
qc_input = QchemInput(monomer,
jobtype='opt',
exchange='hf',
basis='sto-3g',
geom_opt_tol_gradient=300,
geom_opt_tol_energy=100,
geom_opt_coords=-1,
geom_opt_tol_displacement=1200)
parsed_data = get_output_from_qchem(qc_input,
processors=4,
parser=basic_optimization)
opt_monomer = parsed_data['optimized_molecule']
print('Optimized monomer structure')
print(opt_monomer)
# Build dimer from monomer
coor_monomer2 = np.array(opt_monomer.get_coordinates())
coor_monomer2[:, 2] += 4.0 # molecules separation
coordinates = opt_monomer.get_coordinates() + coor_monomer2.tolist()
symbols_dimer = symbols_monomer * 2
dimer = Structure(coordinates=coordinates,
symbols=symbols_dimer,
charge=0,
multiplicity=1)
print('Dimer structure')
print(dimer)
# RASCI qchem input
qc_input = QchemInput(dimer,
jobtype='sp',
exchange='hf',
correlation='rasci',
basis='sto-3g',
ras_act=6,
ras_elec=4,
ras_spin_mult=1, # singlets only
ras_roots=8, # calculate 8 states
ras_do_hole=False,
ras_do_part=False,
set_iter=30)
parsed_data, ee = get_output_from_qchem(qc_input,
processors=4,
force_recalculation=False,
parser=parser_rasci,
return_electronic_structure=True,
)
write_to_fchk(ee, filename='dimer.fchk')
symmetry_measures = get_state_symmetry(ee,
parsed_data['excited_states'],
group='D2h',
orientation=(1, 0, 0),
orientation2=(0, 1, 0),
)
# Analysis of diabatic states to use in diabatization
print('\nSymmetry of the electronic states of the dimer')
list_diabatic = []
for i, state in enumerate(parsed_data['excited_states']):
sym_lab = symmetry_measures[i][0]
energy = state['excitation_energy']
print('State {}: {:3}'.format(i+1, sym_lab))
range_frag = list(range(0, 6))
coordinates_frag = ee['structure'].get_coordinates(fragment=range_frag)
center_frag, normal_frag = get_plane(coordinates_frag)
electronic_structure_frag = crop_electronic_structure(ee, range_frag, renormalize=True)
# save test fchk file with new coefficients
write_to_fchk(electronic_structure_frag, filename='fragment.fchk')
molsym = get_wf_symmetry(electronic_structure_frag['structure'],
electronic_structure_frag['basis'],
electronic_structure_frag['coefficients'],
center=center_frag,
orientation=(1, 0, 0),
orientation2=(0, 1, 0),
group='D2h')
molsym.print_alpha_mo_IRD()
molsym.print_beta_mo_IRD()
ir_labels = molsym.IRLab
def get_symmetry_wf(occupation_alpha, occupation_beta):
state_wf = SymmetryBase(group='D2h', rep='Ag')
for orbital_a, occupation in zip(molsym.mo_SOEVs_a, occupation_alpha):
if abs(occupation) > 0.1:
state_orb = SymmetryBase(group='D2h',
rep=pd.Series(np.array(orbital_a),
index=[ "E", "C2", "C2'", "C2''", "sh", "i", "sv", "sd"]))
state_wf = state_wf * state_orb
for orbital_b, occupation in zip(molsym.mo_SOEVs_a, occupation_beta):
if abs(occupation) > 0.1:
state_orb = SymmetryBase(group='D2h',
rep=pd.Series(np.array(orbital_b),
index=[ "E", "C2", "C2'", "C2''", "sh", "i", "sv", "sd"]))
state_wf = state_wf * state_orb
return state_wf
total_contributions = {}
for istate, state in enumerate(parsed_data['excited_states']):
for conf in state['configurations']:
amplitude = conf['amplitude']
state_wf = get_symmetry_wf(conf['occupations']['alpha'], conf['occupations']['beta'])
rep = state_wf.get_ir_representation().values
ir_labels = state_wf.get_ir_representation().keys()
if 'State: {}'.format(istate+1) in total_contributions:
total_contributions['State: {}'.format(istate+1)] += amplitude**2 * rep
else:
total_contributions['State: {}'.format(istate+1)] = amplitude**2 * rep
print('\nLocal symmetry of electronic states on the monomer')
print(' ' + ' '.join(['{:5}'.format(l) for l in ir_labels]))
for istate, contributions in enumerate(total_contributions.items()):
print(contributions[0], ' '.join(['{:5.2f}'.format(abs(c)) for c in contributions[1]]))
plot_data = pd.DataFrame(total_contributions, index=ir_labels)
plot_data.plot(kind="bar", width=1)
plt.title("Local symmetry of electronic states on the monomer")
plt.xlabel("Irreducible representations")
plt.ylabel("Contribution")
plt.show()