-
Notifications
You must be signed in to change notification settings - Fork 7
/
duschinsky_matrix.py
116 lines (99 loc) · 4.35 KB
/
duschinsky_matrix.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
# Example of the calculation of the Duschinsky matrix for methyl peroxy radical
from pyqchem.qchem_core import get_output_from_qchem
from pyqchem.qc_input import QchemInput
from pyqchem.parsers.parser_frequencies import basic_frequencies
from pyqchem.parsers.parser_optimization import basic_optimization
from pyqchem.tools import get_geometry_from_pubchem
from pyqchem.tools.duschinsky import get_duschinsky
from pyqchem import Structure
import numpy as np
molecule = Structure(coordinates=[[ 1.004123, -0.180454, 0.000000],
[-0.246002, 0.596152, 0.000000],
[-1.312366, -0.230256, 0.000000],
[ 1.810765, 0.567203, 0.000000],
[ 1.036648, -0.805445, -0.904798],
[ 1.036648, -0.805445, 0.904798]],
symbols=['C', 'O', 'O', 'H', 'H', 'H'],
multiplicity=2)
basis_set = '6-31+G*'
n_state = 1 # excited state number
# Optimization of ground state geometry
qc_input = QchemInput(molecule,
jobtype='opt',
exchange='b3lyp',
basis=basis_set,
)
output = get_output_from_qchem(qc_input,
processors=4,
force_recalculation=False,
parser=basic_optimization,
)
# Ground state frequencies calculation
qc_input = QchemInput(output['optimized_molecule'],
jobtype='freq',
exchange='b3lyp',
basis=basis_set,
)
gs_output = get_output_from_qchem(qc_input,
processors=4,
force_recalculation=False,
parser=basic_frequencies,
store_full_output=True,
)
# Optimization of the first excited state state geometry
qc_input = QchemInput(molecule,
jobtype='opt',
exchange='b3lyp',
basis=basis_set,
cis_n_roots=10,
cis_singlets=True,
cis_triplets=False,
#cis_convergence=6,
cis_state_deriv=n_state,
extra_rem_keywords={'XC_GRID': '000075000302'}
)
output = get_output_from_qchem(qc_input,
processors=4,
force_recalculation=False,
parser=basic_optimization,
)
# Excited state frequencies calculation
qc_input = QchemInput(output['optimized_molecule'],
jobtype='freq',
exchange='b3lyp',
basis=basis_set,
cis_n_roots=10,
cis_singlets=True,
cis_triplets=False,
#cis_convergence=6,
cis_state_deriv=n_state,
mem_static=200,
extra_rem_keywords={'XC_GRID': '000075000302'}
)
es_output = get_output_from_qchem(qc_input,
processors=4,
force_recalculation=False,
parser=basic_frequencies,
store_full_output=True,
)
# Calculation of Duschinki rotation matrix according to the model where:
#
# q_f = S * q_i + d
#
# where q_i and q_f are the normal mode coordinates of the initial and final electronic states
# S the rotation matrix and d the vector of displacements
print('frequency displacements')
for mode in gs_output['modes']:
print('{:10}'.format(mode['frequency']), mode['displacement'])
# build duschinsky object
duschinsky = get_duschinsky(gs_output, es_output)
# align structures along principal axis of inertia
duschinsky.align_coordinates_pmi()
# compute and print data
with np.printoptions(precision=3, suppress=True, linewidth=100):
s = duschinsky.get_s_matrix()
print('Rotation matrix')
print(s)
d = duschinsky.get_d_vector()
print('vector of displacements')
print(d[None].T)