Optimization and normal modes#
In this section we will see how to visualize geometry optimization trajectories, as well as molecular vibration modes using py3Dmol.
Geometry optimization#
import veloxchem as vlx
import py3Dmol as p3d
from matplotlib import pyplot as plt
import numpy as np
import ipywidgets
import h5py
Let’s assume we have run the geometry optimization and saved all the steps in xyz format in a file. To animate the optimization trajectory we first need a routine which reads the configurations from file.
def read_xyz_file(file_name):
"""Reads all the configurations from an xyz geometry optimization file."""
xyz_file = open(file_name, "r")
data = xyz_file.read()
xyz_file.close()
xyz_file = open(file_name, "r")
i = 0
steps = []
energies = []
for line in xyz_file:
if "Energy" in line:
steps.append(i)
i += 1
parts = line.split()
energy = float(parts[-1])
energies.append(energy)
xyz_file.close()
return data, steps, energies
Using this routine, we can plot the energies during optimization, as well as animate the trajectory.
xyz_file_name = '../../data/md/kahweol_optim.xyz'
xyz_data, steps, energies = read_xyz_file(xyz_file_name)
# Plot the energies
plt.figure(figsize=(7,4))
plt.plot(steps, energies,'o--')
plt.xlabel('Iteration')
plt.ylabel('Energy (H)')
plt.title("Kahweol XTB optimization")
plt.tight_layout(); plt.show()
# and animate the optimization
viewer = p3d.view(width=600, height=300)
viewer.addModelsAsFrames(xyz_data)
viewer.animate({"loop": "forward"})
viewer.setViewStyle({"style": "outline", "width": 0.05})
viewer.setStyle({"stick": {}, "sphere": {"scale": 0.25}})
viewer.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
If instead we would like to see each configuration individually, we can create an interactive widget with a slider. As in the previous section, we need a routine which reads an individual configuration at a time and returns the corresponding py3Dmol viewer.
def read_xyz_index(file_name, index=0):
"""Reads one configuration defined by its index from the xyz optimization file."""
xyz_file = open(file_name, "r")
current_index = 0
data = ""
# Read the number of atoms from the first line
line = xyz_file.readline()
natm = line.split()[0]
data += line
for line in xyz_file:
if index == current_index:
data += line
if 'Energy' in line:
parts = line.split()
energy = float(parts[-1])
if natm in line:
parts = line.split()
if len(parts) == 1:
if index == current_index:
break
current_index += 1
xyz_file.close()
return data, energy
def return_viewer(file_name, step=0):
xyz_data_i, energy_i = read_xyz_index(file_name, step)
# Uncomment if you would also like to see the energy plot
# or comment out if you would like to skip this.
plt.figure(figsize=(7,4))
plt.plot(steps, energies, 'o--')
plt.plot(step, energy_i, 'o', markersize=15)
plt.xlabel('Iteration')
plt.ylabel('Energy (H)')
plt.title("Kahweol XTB optimization")
plt.tight_layout(); plt.show()
viewer = p3d.view(width=600, height=300)
viewer.addModel(xyz_data_i)
viewer.setViewStyle({"style": "outline", "width": 0.05})
viewer.setStyle({"stick": {}, "sphere": {"scale": 0.25}})
viewer.show()
total_steps = steps[-1] # number of configurations in the trajectory file
# Note that the slider only works in a Jupyter notebook.
ipywidgets.interact(return_viewer, file_name=xyz_file_name,
step=ipywidgets.IntSlider(min=0, max=total_steps, step=1, value=3))
<function __main__.return_viewer(file_name, step=0)>
Vibrational normal modes#
Now let’s see how to animate and visualize vibrational normal modes. For this, we will need to determine the molecular Hessian and perform a vibrational analysis.
xtb_drv = vlx.XtbDriver()
xtb_vibanalysis_drv = vlx.VibrationalAnalysis(xtb_drv)
In order to perform a vibrational analysis we need a Molecule object which stores the optimized geometry. To create the Molecule object we will read the last geometry from the xyz file we used in the previous subsection.
# read the last configuration from file
kahweol_xyz, opt_energy = read_xyz_index(file_name=xyz_file_name, index=total_steps)
# create the Mlecul object
xtb_opt_kahweol = vlx.Molecule.read_xyz_string(kahweol_xyz)
# perform vibrational analysis -- diagonalize Hessian, extract frequencies and normal modes
# and calculate IR intensities.
xtb_vibanalysis_drv.compute(xtb_opt_kahweol)
Show code cell output
Vibrational Analysis Driver
=============================
The following will be computed:
- Vibrational frequencies and normal modes
- Force constants
- IR intensities
XTB Hessian Driver
====================
Hessian Type : Numerical
Numerical Method : Symmetric Difference Quotient
Finite Difference Step Size : 0.001 a.u.
*** Time spent in Hessian calculation: 10.80 sec ***
Free Energy Analysis
======================
Note: Rotational symmetry is set to 1 regardless of true symmetry
No Imaginary Frequencies
Free energy contributions calculated at @ 298.15 K:
Zero-point vibrational energy: 263.1952 kcal/mol
H (Trans + Rot + Vib = Tot): 1.4812 + 0.8887 + 10.5627 = 12.9326 kcal/mol
S (Trans + Rot + Vib = Tot): 43.1588 + 34.2859 + 62.5822 = 140.0269 cal/mol/K
TS (Trans + Rot + Vib = Tot): 12.8678 + 10.2223 + 18.6589 = 41.7490 kcal/mol
Ground State Electronic Energy : E0 = -68.32150638 au ( -42872.3925 kcal/mol)
Free Energy Correction (Harmonic) : ZPVE + [H-TS]_T,R,V = 0.37350635 au ( 234.3788 kcal/mol)
Gibbs Free Energy (Harmonic) : E0 + ZPVE + [H-TS]_T,R,V = -67.94800003 au ( -42638.0138 kcal/mol)
Vibrational Analysis
======================
* Info * The 5 dominant normal modes are printed below.
Vibrational Mode 11
----------------------------------------------------
Harmonic frequency: 209.80 cm**-1
Reduced mass: 1.2470 amu
Force constant: 0.0323 mdyne/A
IR intensity: 211.0478 km/mol
Normal mode:
X Y Z
1 O 0.0092 -0.0315 -0.0211
2 O -0.0272 0.0106 0.0491
3 O 0.0320 -0.0118 -0.0044
4 C -0.0009 -0.0003 0.0124
5 C 0.0003 -0.0028 -0.0095
6 C -0.0177 -0.0096 0.0062
7 C -0.0022 0.0100 0.0127
8 C -0.0094 -0.0226 0.0175
9 C 0.0083 0.0095 -0.0136
10 C -0.0183 -0.0189 0.0036
11 C 0.0073 0.0026 0.0175
12 C -0.0209 -0.0126 -0.0245
13 C -0.0281 -0.0055 -0.0235
14 C 0.0177 0.0137 -0.0120
15 C 0.0223 0.0116 -0.0120
16 C -0.0082 -0.0035 -0.0105
17 C -0.0146 0.0184 0.0439
18 C 0.0252 0.0105 0.0019
19 C 0.0178 0.0085 -0.0007
20 C 0.0278 0.0045 0.0037
21 C 0.0309 0.0008 -0.0067
22 C 0.0072 0.0013 0.0129
23 C 0.0160 -0.0118 0.0089
24 H -0.0030 -0.0160 -0.0104
25 H -0.0152 -0.0190 0.0039
26 H -0.0055 0.0286 0.0099
27 H 0.0023 0.0110 0.0217
28 H -0.0159 -0.0469 0.0001
29 H -0.0002 -0.0345 0.0501
30 H -0.0121 -0.0111 0.0153
31 H 0.0209 0.0028 0.0398
32 H -0.0174 -0.0343 -0.0329
33 H -0.0348 0.0000 -0.0347
34 H -0.0312 0.0169 -0.0255
35 H -0.0300 -0.0082 -0.0347
36 H 0.0241 0.0145 -0.0121
37 H 0.0436 0.0028 -0.0102
38 H 0.0182 0.0142 -0.0289
39 H -0.0281 0.0940 -0.0073
40 H -0.0938 -0.0690 0.0065
41 H 0.0772 -0.0474 -0.0284
42 H -0.0060 0.0932 0.0140
43 H 0.0057 0.0187 0.1101
44 H 0.0330 0.0227 0.0158
45 H -0.5900 0.2696 -0.4528
46 H 0.0270 0.0107 0.0157
47 H -0.0057 0.0047 0.0236
48 H -0.3196 0.1764 -0.3725
49 H 0.0126 -0.0210 0.0134
Vibrational Mode 119
----------------------------------------------------
Harmonic frequency: 2925.13 cm**-1
Reduced mass: 1.0756 amu
Force constant: 5.4222 mdyne/A
IR intensity: 75.6487 km/mol
Normal mode:
X Y Z
1 O 0.0002 0.0002 -0.0003
2 O 0.0001 0.0001 -0.0002
3 O -0.0000 -0.0000 -0.0000
4 C 0.0001 0.0003 -0.0003
5 C -0.0002 0.0000 0.0016
6 C -0.0475 -0.0027 0.0619
7 C 0.0004 0.0036 0.0003
8 C 0.0002 0.0005 -0.0008
9 C 0.0001 0.0000 0.0001
10 C -0.0018 -0.0003 0.0019
11 C 0.0001 -0.0001 -0.0010
12 C -0.0001 0.0000 -0.0012
13 C -0.0012 -0.0035 -0.0021
14 C -0.0000 0.0000 0.0002
15 C -0.0001 0.0001 0.0001
16 C 0.0000 -0.0000 0.0000
17 C 0.0007 0.0012 0.0020
18 C -0.0000 0.0000 -0.0000
19 C 0.0000 -0.0000 0.0000
20 C -0.0000 -0.0000 0.0000
21 C 0.0000 0.0000 0.0000
22 C -0.0000 0.0000 -0.0000
23 C -0.0000 0.0000 -0.0000
24 H 0.0028 0.0005 -0.0213
25 H 0.6034 0.0385 -0.7855
26 H 0.0078 0.0064 0.0083
27 H -0.0224 -0.0533 0.0080
28 H -0.0012 -0.0066 0.0093
29 H -0.0003 0.0011 0.0002
30 H -0.0020 0.0006 0.0135
31 H 0.0002 -0.0003 -0.0006
32 H 0.0062 -0.0015 0.0081
33 H -0.0022 -0.0019 0.0001
34 H 0.0328 0.0055 0.0516
35 H -0.0362 0.0306 -0.0059
36 H 0.0001 -0.0001 -0.0026
37 H 0.0000 0.0000 -0.0007
38 H 0.0013 -0.0019 -0.0007
39 H -0.0006 -0.0001 -0.0002
40 H -0.0001 -0.0002 0.0004
41 H -0.0003 0.0005 0.0001
42 H -0.0079 -0.0107 -0.0241
43 H -0.0025 -0.0037 0.0019
44 H 0.0004 -0.0002 0.0001
45 H 0.0047 0.0029 0.0002
46 H 0.0003 0.0006 -0.0002
47 H -0.0000 -0.0002 0.0000
48 H 0.0004 -0.0003 0.0010
49 H 0.0002 -0.0001 0.0001
Vibrational Mode 125
----------------------------------------------------
Harmonic frequency: 2979.68 cm**-1
Reduced mass: 1.0740 amu
Force constant: 5.6181 mdyne/A
IR intensity: 86.4704 km/mol
Normal mode:
X Y Z
1 O -0.0000 0.0002 -0.0002
2 O -0.0000 0.0000 -0.0000
3 O 0.0000 0.0000 0.0000
4 C -0.0012 -0.0026 0.0004
5 C 0.0001 0.0002 -0.0013
6 C -0.0030 -0.0016 0.0024
7 C -0.0241 -0.0565 0.0074
8 C 0.0018 0.0090 -0.0001
9 C -0.0001 0.0000 0.0000
10 C -0.0001 -0.0000 0.0002
11 C -0.0135 -0.0283 -0.0026
12 C -0.0122 -0.0039 -0.0123
13 C -0.0108 0.0057 -0.0070
14 C 0.0002 -0.0004 -0.0007
15 C 0.0109 -0.0162 -0.0132
16 C -0.0013 -0.0003 0.0012
17 C -0.0005 -0.0008 0.0001
18 C -0.0002 0.0003 -0.0001
19 C -0.0000 -0.0001 0.0000
20 C -0.0000 -0.0001 0.0001
21 C 0.0000 0.0001 -0.0000
22 C 0.0001 0.0002 -0.0000
23 C -0.0000 0.0000 -0.0000
24 H -0.0022 -0.0012 0.0160
25 H 0.0266 0.0019 -0.0331
26 H 0.0124 -0.0209 0.0304
27 H 0.2875 0.7255 -0.1240
28 H -0.0019 -0.0153 0.0285
29 H -0.0189 -0.0941 -0.0261
30 H -0.0289 -0.0084 0.1649
31 H 0.2033 0.3600 -0.1353
32 H 0.0930 -0.0342 0.1351
33 H 0.0573 0.0852 0.0140
34 H 0.0445 0.0173 0.0781
35 H 0.0911 -0.0885 0.0054
36 H 0.0003 0.0004 0.0040
37 H -0.0053 -0.0071 0.0935
38 H -0.1392 0.2018 0.0680
39 H -0.0019 -0.0002 0.0006
40 H 0.0101 -0.0136 -0.0032
41 H 0.0063 0.0170 -0.0100
42 H -0.0003 -0.0004 0.0001
43 H 0.0049 0.0092 -0.0011
44 H 0.0015 -0.0023 0.0011
45 H 0.0008 0.0000 0.0008
46 H 0.0007 0.0014 -0.0006
47 H -0.0004 -0.0026 0.0006
48 H 0.0003 -0.0000 0.0001
49 H 0.0000 0.0001 -0.0000
Vibrational Mode 128
----------------------------------------------------
Harmonic frequency: 2999.40 cm**-1
Reduced mass: 1.0753 amu
Force constant: 5.6994 mdyne/A
IR intensity: 72.5337 km/mol
Normal mode:
X Y Z
1 O -0.0001 0.0003 0.0002
2 O 0.0001 0.0001 -0.0001
3 O 0.0000 0.0000 0.0000
4 C -0.0006 -0.0018 0.0006
5 C -0.0000 -0.0011 0.0028
6 C 0.0001 -0.0004 0.0001
7 C -0.0027 -0.0065 0.0010
8 C -0.0122 -0.0609 0.0157
9 C -0.0001 -0.0002 0.0000
10 C -0.0000 -0.0018 0.0002
11 C -0.0019 -0.0039 -0.0002
12 C -0.0040 -0.0370 0.0181
13 C 0.0072 -0.0107 -0.0023
14 C 0.0001 0.0000 -0.0000
15 C 0.0002 0.0000 -0.0035
16 C -0.0010 0.0004 -0.0001
17 C 0.0022 0.0041 -0.0009
18 C -0.0009 0.0011 -0.0004
19 C -0.0000 -0.0000 -0.0000
20 C -0.0000 -0.0001 0.0000
21 C -0.0001 -0.0000 -0.0000
22 C -0.0000 0.0000 -0.0000
23 C -0.0000 0.0000 -0.0000
24 H 0.0039 0.0009 -0.0365
25 H 0.0011 -0.0003 -0.0018
26 H 0.0014 -0.0018 0.0040
27 H 0.0330 0.0836 -0.0156
28 H 0.0310 0.2286 -0.3558
29 H 0.1182 0.5297 0.1645
30 H -0.0041 -0.0013 0.0219
31 H 0.0287 0.0534 -0.0195
32 H -0.2448 0.0602 -0.3385
33 H 0.2999 0.3955 0.1081
34 H 0.0357 0.0043 0.0534
35 H -0.1315 0.1198 -0.0130
36 H 0.0002 -0.0000 -0.0008
37 H -0.0048 0.0006 0.0457
38 H 0.0012 -0.0014 -0.0025
39 H 0.0096 0.0024 0.0024
40 H 0.0037 -0.0049 -0.0011
41 H -0.0013 -0.0010 0.0004
42 H 0.0018 0.0036 0.0020
43 H -0.0291 -0.0518 0.0065
44 H 0.0087 -0.0106 0.0061
45 H 0.0018 -0.0010 -0.0023
46 H 0.0008 0.0015 -0.0006
47 H -0.0000 -0.0001 0.0000
48 H 0.0003 -0.0003 -0.0001
49 H -0.0001 0.0001 -0.0000
Vibrational Mode 130
----------------------------------------------------
Harmonic frequency: 3009.03 cm**-1
Reduced mass: 1.0598 amu
Force constant: 5.6535 mdyne/A
IR intensity: 69.9235 km/mol
Normal mode:
X Y Z
1 O -0.0001 0.0000 0.0002
2 O -0.0000 -0.0000 0.0001
3 O -0.0000 -0.0001 0.0000
4 C 0.0001 0.0000 -0.0013
5 C 0.0001 0.0002 -0.0014
6 C 0.0004 0.0001 -0.0004
7 C 0.0022 0.0016 0.0020
8 C -0.0021 -0.0021 -0.0350
9 C 0.0000 -0.0003 -0.0002
10 C 0.0000 -0.0001 -0.0009
11 C 0.0024 0.0030 -0.0085
12 C -0.0018 0.0002 -0.0025
13 C -0.0005 -0.0008 -0.0018
14 C 0.0001 0.0006 -0.0027
15 C -0.0024 0.0119 -0.0563
16 C -0.0021 -0.0043 0.0022
17 C -0.0014 -0.0024 -0.0004
18 C 0.0001 -0.0001 0.0001
19 C -0.0001 0.0001 -0.0000
20 C -0.0000 -0.0001 0.0001
21 C 0.0000 0.0001 -0.0000
22 C -0.0001 -0.0004 0.0001
23 C -0.0000 0.0001 -0.0000
24 H -0.0024 -0.0020 0.0190
25 H -0.0050 -0.0004 0.0068
26 H -0.0236 -0.0092 -0.0261
27 H -0.0047 -0.0126 0.0035
28 H -0.0397 -0.2731 0.3656
29 H 0.0649 0.2972 0.0659
30 H -0.0119 0.0043 0.0774
31 H -0.0227 -0.0393 0.0113
32 H 0.0214 -0.0074 0.0314
33 H 0.0022 0.0044 -0.0003
34 H 0.0135 0.0033 0.0218
35 H -0.0072 0.0062 -0.0014
36 H 0.0009 0.0012 0.0168
37 H -0.0807 0.0135 0.7877
38 H 0.1099 -0.1555 -0.0869
39 H 0.0089 -0.0000 0.0035
40 H -0.0017 -0.0006 0.0016
41 H 0.0188 0.0482 -0.0304
42 H 0.0036 0.0024 0.0099
43 H 0.0146 0.0265 -0.0039
44 H -0.0016 0.0019 -0.0009
45 H 0.0015 -0.0001 -0.0028
46 H 0.0005 0.0006 -0.0004
47 H 0.0011 0.0043 -0.0012
48 H 0.0000 -0.0000 -0.0002
49 H 0.0007 -0.0004 0.0002
* Info * Full vibrational analysis results written to: vib-results.out
{'molecule_xyz_string': '49\n\nO -4.444356148300 -1.519231405700 0.133818803000\nO -5.717728473400 0.942316237100 -0.066968962400\nO 5.074944380100 0.621839630700 0.095361864000\nC -1.218663248700 -0.830595956300 -0.247223992800\nC -0.477629871500 0.497184226800 -0.529673918500\nC -2.930011813600 0.109407860600 1.091647841000\nC -1.787420846300 -0.894913983900 1.174595923900\nC -2.515711138500 -0.855946194500 -1.095040548800\nC 0.968464285400 0.594394878100 0.023051359100\nC -3.654907965800 -0.372192923700 -0.180394141400\nC -0.338906858500 -2.036979464400 -0.583399668300\nC -1.380813476500 1.698304225800 -0.180085905800\nC -2.322111192600 1.507002611300 1.012160927000\nC 1.723598180200 -0.686836785300 -0.399658653300\nC 1.033616297100 -1.956732648500 0.075560115900\nC 1.068852139200 0.833440819900 1.537692826600\nC -4.585050097800 0.655858178700 -0.853861880100\nC 1.691896058500 1.770006746900 -0.620577515400\nC 3.162347118000 -0.539762374900 -0.024338881700\nC 3.026560322800 1.832733770700 -0.661355332100\nC 3.761395103300 0.686565534800 -0.201084928000\nC 4.186018037300 -1.411360657400 0.402952793600\nC 5.318365378800 -0.649258410900 0.464736569000\nH -0.335036180700 0.523105135000 -1.619153336000\nH -3.599176954100 0.065966348100 1.957690684700\nH -1.080856355100 -0.645807578400 1.956928694700\nH -2.177294049000 -1.897906384700 1.358785220600\nH -2.433801681300 -0.246547664800 -1.993148048300\nH -2.750703334400 -1.876555715600 -1.397466956800\nH -0.203126601800 -2.084320493000 -1.667436130400\nH -0.851890375400 -2.949981497800 -0.273589578900\nH -1.983411213200 1.924016483600 -1.061330676900\nH -0.767945077200 2.580213930300 0.010339706900\nH -1.779556090200 1.670505612800 1.945822615500\nH -3.112392348700 2.258845075600 0.965603116300\nH 1.702548680000 -0.706998814100 -1.502646038800\nH 0.940906440200 -1.969012908100 1.160909209200\nH 1.629806822600 -2.823430346800 -0.216840792700\nH 2.103739193800 1.037127297700 1.803288498700\nH 0.473184324800 1.693615965100 1.827260296200\nH 0.749302119400 -0.028513839400 2.112157924200\nH -4.890341941000 0.251895622800 -1.829136457400\nH -4.081362417900 1.607147612800 -1.012173100000\nH 1.101176933700 2.597977027300 -0.978171475100\nH -4.939036199800 -1.335514044100 0.940761097800\nH 3.561011194900 2.688548100800 -1.042815340900\nH 4.105646417200 -2.452204647200 0.640894394400\nH -6.244800187700 0.136371594000 -0.009353562100\nH 6.323262711400 -0.890385787700 0.749405341100\n',
'gibbs_free_energy': -67.94800002675646,
'free_energy_summary': 'Note: Rotational symmetry is set to 1 regardless of true symmetry\nNo Imaginary Frequencies\n\nFree energy contributions calculated at @ 298.15 K:\nZero-point vibrational energy: 263.1952 kcal/mol\nH (Trans + Rot + Vib = Tot): 1.4812 + 0.8887 + 10.5627 = 12.9326 kcal/mol\nS (Trans + Rot + Vib = Tot): 43.1588 + 34.2859 + 62.5822 = 140.0269 cal/mol/K\nTS (Trans + Rot + Vib = Tot): 12.8678 + 10.2223 + 18.6589 = 41.7490 kcal/mol\n\nGround State Electronic Energy : E0 = -68.32150638 au ( -42872.3925 kcal/mol)\nFree Energy Correction (Harmonic) : ZPVE + [H-TS]_T,R,V = 0.37350635 au ( 234.3788 kcal/mol)\nGibbs Free Energy (Harmonic) : E0 + ZPVE + [H-TS]_T,R,V = -67.94800003 au ( -42638.0138 kcal/mol)\n',
'vib_frequencies': array([ 44.72005936, 50.43202835, 80.29612201, 92.89651977,
134.984309 , 147.48381763, 162.65837666, 190.61455171,
194.86464907, 207.19570858, 209.8044969 , 267.87742003,
271.09026223, 271.46471476, 294.51188597, 301.76664605,
319.20782579, 333.41623758, 357.09521052, 367.61007014,
377.12228617, 381.43564697, 416.84275546, 434.14389834,
447.29954585, 454.30077945, 475.53017159, 506.42525822,
523.94203141, 527.58038219, 561.01028189, 592.24413266,
604.82349516, 628.26916045, 640.70079584, 654.45980158,
685.0265179 , 744.39863576, 749.04124238, 778.29835625,
779.33707256, 789.81218231, 819.16445809, 832.40053367,
834.17279371, 857.90367398, 883.7022011 , 891.85996186,
893.82602348, 900.72653606, 910.02251838, 931.62773841,
942.7653207 , 949.4105173 , 968.35812188, 996.11537957,
1000.20049451, 1023.94672842, 1030.89566603, 1043.35500065,
1048.27578544, 1063.15808523, 1066.85065601, 1078.74911169,
1085.87566787, 1089.80512104, 1090.06042283, 1113.32218997,
1116.47820536, 1140.79505188, 1144.19822605, 1152.23175987,
1155.30529029, 1167.54092271, 1174.26419979, 1178.95014652,
1188.70382454, 1196.77869308, 1198.37992748, 1216.86168022,
1221.92760712, 1240.31460513, 1242.42215869, 1267.06313994,
1271.31008962, 1278.74976758, 1283.39045403, 1294.95859071,
1301.17372806, 1302.69697265, 1310.24123595, 1322.32233513,
1332.62156273, 1339.90403463, 1342.42685345, 1349.37976487,
1352.30313949, 1358.76048189, 1359.87983732, 1374.09796183,
1376.6821786 , 1403.84703714, 1446.48750935, 1459.27496971,
1467.69049341, 1475.43724488, 1477.89578385, 1484.7827952 ,
1493.83012469, 1496.95054015, 1502.40874397, 1512.89843416,
1515.31350985, 1569.09596513, 1637.47804163, 2833.71660634,
2884.03523884, 2891.2902753 , 2925.13216031, 2954.30151087,
2962.49946225, 2971.25806859, 2976.74734794, 2977.66578056,
2979.67889968, 2986.27283685, 2996.2852152 , 2999.4037443 ,
3007.88102936, 3009.02623374, 3030.1820168 , 3036.13720774,
3056.33459653, 3063.95088458, 3076.54023087, 3086.53386038,
3101.10118265, 3166.31368842, 3184.03828961, 3525.83374657,
3530.45886643]),
'normal_modes': array([[-6.36953767e-03, -4.51460058e-02, 7.26002800e-04, ...,
-5.55247418e-02, 8.71853593e-02, 2.19678739e-01],
[ 6.89364212e-02, 3.49147830e-03, 7.91868171e-02, ...,
-8.47277016e-02, -3.21027534e-02, 2.86827525e-01],
[ 9.16692932e-02, -9.02908843e-02, 4.74877224e-02, ...,
2.05998176e-02, -1.47411969e-01, -1.57985015e-01],
...,
[ 5.40542987e-06, 9.35994682e-06, -4.95876433e-07, ...,
-3.68279575e-01, 8.87679018e-02, -1.04349014e-01],
[-3.63625808e-03, 1.38208280e-03, 5.74177359e-03, ...,
1.09668393e-05, -1.23141149e-07, 2.91248057e-06],
[ 3.09676007e-02, -1.20119885e-02, -5.16017281e-02, ...,
1.05282412e-05, -1.87391374e-07, 4.04924113e-06]]),
'reduced_masses': array([4.84831735, 5.08371945, 3.62737455, 4.47435104, 2.76420264,
3.69929211, 4.12642354, 1.98132933, 2.78855431, 3.64304255,
1.24699694, 2.59018842, 2.88532347, 2.85879802, 3.13021897,
2.9888019 , 2.56943163, 2.41912084, 2.81918047, 3.38265503,
2.34207592, 2.71434118, 3.91178041, 2.44527809, 2.59159129,
2.55934104, 1.4136785 , 3.49807953, 3.44611558, 2.83549559,
3.32208912, 4.32686357, 5.2628915 , 3.59253348, 4.49404951,
3.90965461, 3.36850457, 2.64260648, 1.3389677 , 2.79384215,
1.90389555, 3.01611341, 2.9630427 , 1.48555217, 2.22414857,
2.13687412, 1.90690435, 1.38664093, 2.4473202 , 2.29138065,
2.13057412, 2.20706847, 1.96280709, 1.71547878, 1.99126239,
2.02184943, 1.83255064, 1.70491891, 2.12677113, 2.02600072,
1.91015432, 1.80319905, 1.75782426, 1.92011631, 1.501021 ,
1.99825599, 3.06697029, 1.94746491, 1.62414165, 1.37992291,
1.74847076, 1.96496408, 1.59258227, 1.88171632, 1.46247758,
2.4200758 , 1.94061649, 1.46786471, 1.56932745, 1.97665948,
2.13255183, 2.37845309, 2.06709969, 2.07179518, 1.74112541,
1.35505129, 1.72815058, 1.57442149, 1.41251539, 1.78975838,
1.64652632, 1.30086909, 1.49118764, 1.45382665, 1.76128006,
1.42208981, 1.95813151, 1.9336556 , 1.38917325, 1.37327749,
1.35346505, 1.23005572, 7.17261487, 1.08327152, 1.09563733,
1.09370416, 5.57401328, 1.10786379, 1.08750666, 1.08687993,
1.0931636 , 1.07541699, 1.07426352, 9.80494491, 8.95849704,
1.07191195, 1.07450382, 1.07208372, 1.07556775, 1.08754167,
1.09113618, 1.06990815, 1.09751732, 1.07513784, 1.07398771,
1.05189515, 1.07329215, 1.07525687, 1.05458276, 1.05977493,
1.07037099, 1.08597648, 1.0795633 , 1.04842436, 1.06862175,
1.07634202, 1.08882258, 1.08116125, 1.09394433, 1.06567768,
1.06520664]),
'force_constants': array([5.71276057e-03, 7.61806523e-03, 1.37794469e-02, 2.27498838e-02,
2.96747243e-02, 4.74086521e-02, 6.43245727e-02, 4.24149552e-02,
6.23871961e-02, 9.21458957e-02, 3.23403871e-02, 1.09510008e-01,
1.24931681e-01, 1.24125351e-01, 1.59967078e-01, 1.60357709e-01,
1.54253265e-01, 1.58446006e-01, 2.11807419e-01, 2.69328771e-01,
1.96252684e-01, 2.32679022e-01, 4.00469295e-01, 2.71547557e-01,
3.05501724e-01, 3.11218471e-01, 1.88346280e-01, 5.28579890e-01,
5.57373810e-01, 4.65003789e-01, 6.16031841e-01, 8.94180048e-01,
1.13431053e+00, 8.35492486e-01, 1.08692260e+00, 9.86630610e-01,
9.31326903e-01, 8.62767691e-01, 4.42620785e-01, 9.97112893e-01,
6.81308851e-01, 1.10852518e+00, 1.17146779e+00, 6.06460957e-01,
9.11855614e-01, 9.26629742e-01, 8.77386653e-01, 6.49841722e-01,
1.15198553e+00, 1.09530078e+00, 1.03956384e+00, 1.12862812e+00,
1.02786258e+00, 9.11053053e-01, 1.10014703e+00, 1.18200235e+00,
1.08014076e+00, 1.05319477e+00, 1.33168167e+00, 1.29943344e+00,
1.23671550e+00, 1.20085232e+00, 1.17878052e+00, 1.31649338e+00,
1.04279085e+00, 1.39829580e+00, 2.14714292e+00, 1.42220307e+00,
1.19281926e+00, 1.05808442e+00, 1.34868695e+00, 1.53703777e+00,
1.25240745e+00, 1.51129281e+00, 1.18814952e+00, 1.98184685e+00,
1.61561277e+00, 1.23869381e+00, 1.32786176e+00, 1.72450528e+00,
1.87603437e+00, 2.15580039e+00, 1.87996619e+00, 1.95971796e+00,
1.65799525e+00, 1.30550070e+00, 1.67706320e+00, 1.55554655e+00,
1.40900976e+00, 1.78949897e+00, 1.66541106e+00, 1.34016571e+00,
1.56025725e+00, 1.53783682e+00, 1.87007886e+00, 1.52561746e+00,
2.10979476e+00, 2.10336760e+00, 1.51358816e+00, 1.52772059e+00,
1.51134867e+00, 1.42828411e+00, 8.84213472e+00, 1.35913261e+00,
1.39054812e+00, 1.40278653e+00, 7.17308305e+00, 1.43900552e+00,
1.42983061e+00, 1.43498285e+00, 1.45382324e+00, 1.45026271e+00,
1.45333610e+00, 1.42231048e+01, 1.41526043e+01, 5.07133366e+00,
5.26573900e+00, 5.28034541e+00, 5.42224326e+00, 5.59249723e+00,
5.64216469e+00, 5.56515782e+00, 5.72988063e+00, 5.61650679e+00,
5.61808730e+00, 5.52690068e+00, 5.67720392e+00, 5.69944180e+00,
5.62150004e+00, 5.65347962e+00, 5.79057915e+00, 5.89811790e+00,
5.94155533e+00, 5.79897106e+00, 5.95935761e+00, 6.04147002e+00,
6.16934747e+00, 6.38628931e+00, 6.53434446e+00, 7.80548278e+00,
7.82251520e+00]),
'ir_intensities': array([1.47168373e+00, 1.21276911e+00, 3.44051135e-01, 5.22795905e-01,
1.42065811e+00, 1.28175381e+00, 1.52309925e-01, 6.57910704e-01,
1.85988562e+00, 6.18982408e+01, 2.11047847e+02, 7.14047524e+00,
7.85562705e+00, 4.33704700e+00, 8.98609566e+00, 3.00354608e+00,
1.49681686e+00, 2.81983763e+00, 5.14912610e-01, 4.35122216e+00,
3.45035726e+00, 3.47837357e+00, 1.46578858e+00, 5.87040164e+00,
4.56635204e+00, 9.28865038e+00, 1.60138436e+01, 4.33099259e-01,
2.62401040e+00, 3.97402802e+00, 3.48504774e+00, 2.46783948e+00,
1.08744162e+00, 2.82554408e+00, 9.99716576e-01, 5.78708509e+00,
1.15035596e+01, 7.32443880e+00, 7.91866234e+00, 3.31641467e+00,
6.17371001e+00, 2.64369000e+00, 6.28805656e+00, 8.94590437e+00,
4.69883262e+00, 1.87074672e+00, 2.02855918e+00, 2.18361602e+00,
4.86102605e+00, 3.48667200e+00, 1.40051292e+01, 2.92311285e+00,
3.38260409e+00, 9.00192773e+00, 1.27951473e+01, 7.46259019e+00,
1.10212582e+01, 6.32025040e+00, 2.72999022e+01, 1.15664965e+01,
6.45147625e+00, 2.23544913e+01, 3.62386706e+01, 6.86185204e+01,
3.90981649e+00, 2.82245348e+01, 4.18326304e+01, 1.14235203e+01,
3.35891528e+00, 8.96243787e-01, 8.25707583e+00, 2.38913722e+00,
3.33210420e+00, 1.27495530e+01, 1.75322585e+01, 1.88863132e+00,
2.46405209e+01, 1.01398918e+01, 9.06679064e+00, 1.72780927e+01,
5.29932552e+00, 1.97090177e+01, 1.11124390e+01, 3.91942966e+01,
2.22562217e+01, 9.93752721e+00, 1.86046778e+01, 2.86266928e+00,
6.66854252e-01, 7.78169799e-01, 3.99775165e-01, 6.89666028e+01,
1.00081243e+01, 1.96968930e+01, 4.98251478e+00, 4.10671379e-01,
1.67135718e+01, 1.60231645e+01, 2.56694488e+00, 2.71799582e-01,
9.98741633e-01, 5.91697325e+00, 7.51504511e+00, 9.73321477e-01,
1.22468411e+00, 2.75804880e+00, 1.96912247e+01, 8.13362206e-01,
1.27036049e+00, 6.51090311e-01, 2.37293484e-01, 2.74688509e+00,
1.98545273e+00, 1.69299787e+00, 8.87917107e+00, 4.19955744e+01,
4.09377866e+01, 5.98580419e+01, 7.56487074e+01, 3.46790445e+01,
6.84130693e+00, 1.71842834e+01, 5.99381847e+01, 4.03264054e+01,
8.64704481e+01, 4.42256357e+01, 4.49844500e+00, 7.25337388e+01,
1.01394705e+01, 6.99234645e+01, 2.44918203e+01, 3.37072547e+01,
2.39799233e+01, 8.18353596e+00, 3.78711767e+01, 1.16436287e+01,
5.49630803e+01, 5.58580246e+00, 2.37515103e+01, 1.20886048e+01,
9.87890558e+00])}
Now we have the frequencies, normal modes and IR intensities. In order to animate a specific normal mode, we need to write a routine which selects the mode of interest and returns it in the format required by py3Dmol. We would also like to plot the IR spectrum and, for this, we also need a routine which adds a Gaussian or Lorentzian broadening to the calculated IR spectrum.
# To animate the normal mode we will need both the geometry and the displacements
def get_normal_mode(molecule, normal_mode):
elements = molecule.get_labels()
coords = molecule.get_coordinates_in_angstrom()
natm = molecule.number_of_atoms()
vib_xyz = "%d\n\n" % natm
nm = normal_mode.reshape(natm, 3)
for i in range(natm):
# add coordinates:
vib_xyz += elements[i] + " %15.7f %15.7f %15.7f " % (coords[i,0], coords[i,1], coords[i,2])
# add displacements:
vib_xyz += "%15.7f %15.7f %15.7f\n" % (nm[i,0], nm[i,1], nm[i,2])
return vib_xyz
# Broadening function
def add_broadening(list_ex_energy, list_osci_strength, line_profile='Lorentzian', line_param=10, step=10):
x_min = np.amin(list_ex_energy) - 50
x_max = np.amax(list_ex_energy) + 50
x = np.arange(x_min, x_max, step)
y = np.zeros((len(x)))
#print(x)
#print(y)
# go through the frames and calculate the spectrum for each frame
for xp in range(len(x)):
for e, f in zip(list_ex_energy, list_osci_strength):
if line_profile == 'Gaussian':
y[xp] += f * np.exp(-(
(e - x[xp]) / line_param)**2)
elif line_profile == 'Lorentzian':
y[xp] += 0.5 * line_param * f / (np.pi * (
(x[xp] - e)**2 + 0.25 * line_param**2))
return x, y
After adding the broadening, we can plot the spectrum and then animate a specific normal mode, selected based on its index.
wvn, ir = xtb_vibanalysis_drv.vib_frequencies, xtb_vibanalysis_drv.ir_intensities
wvng, irg = add_broadening(wvn, ir, line_profile='Gaussian', line_param=10, step=2)
# Plot the IR spectra
plt.figure(figsize=(7,4))
# Plot the IR spectrum
plt.plot(wvng, irg, label='IR spectrum')
plt.xlabel('Wavenumber (cm**-1)')
#plt.axis(xmin=3200, xmax=3500)
#plt.axis(ymin=-0.2, ymax=1.5)
plt.ylabel('IR intensity (km/mol)')
plt.title("Calculated IR sepctrum, Kahweol")
plt.legend()
plt.tight_layout(); plt.show()
# Visualize the vibrational mode
index = 64
print("\n\n\n Normal mode %d: %.2f cm-1, %.2f km/mol." % (index,
xtb_vibanalysis_drv.vib_frequencies[index-1],
xtb_vibanalysis_drv.ir_intensities[index-1]))
normal_mode = get_normal_mode(xtb_opt_kahweol, xtb_vibanalysis_drv.normal_modes[index-1])
view = p3d.view(viewergrid=(1,1), width=600, height=300)
view.addModel(normal_mode, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}})
view.setViewStyle({"style": "outline", "width": 0.05})
view.setStyle({"stick": {}, "sphere": {"scale": 0.25}})
view.animate({'loop': 'backAndForth'})
view.zoomTo()
view.show()
Normal mode 64: 1078.75 cm-1, 68.62 km/mol.
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
To make a more interactive normal mode selection, we can create a slider using the ipywidgets Python module for Jupyter notebooks. To make use of the slider widget, we first need a routine which can plot the IR spectrum and animate a specific normal mode selected based on its index.
def vibration_viewer(index):
freq = xtb_vibanalysis_drv.vib_frequencies[index-1]
ir_intens = xtb_vibanalysis_drv.ir_intensities[index-1]
# Plot the IR spectrum
plt.figure(figsize=(7,4))
plt.plot(wvng, irg)
plt.plot(freq, ir_intens, 'o', markersize=10)
plt.xlabel('Wavenumber (cm**-1)')
plt.ylabel('IR intensity (km/mol)')
plt.title("Calculated IR sepctrum, Kahweol")
plt.tight_layout(); plt.show()
normal_mode = get_normal_mode(xtb_opt_kahweol, xtb_vibanalysis_drv.normal_modes[index-1])
view = p3d.view(viewergrid=(1,1), width=600, height=300)
view.addModel(normal_mode, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}})
view.setViewStyle({"style": "outline", "width": 0.05})
view.setStyle({"stick": {}, "sphere": {"scale": 0.25}})
view.animate({'loop': 'backAndForth'})
view.zoomTo()
view.show()
no_norm_modes = xtb_vibanalysis_drv.vib_frequencies.shape[0] # number of vibrational modes
# Note that the slider only works in a Jupyter notebook.
ipywidgets.interact(vibration_viewer,
index=ipywidgets.IntSlider(min=1, max=no_norm_modes, step=1, value=92))
<function __main__.vibration_viewer(index)>