Save and load data using HDF5 files
Some of richmol objects can be stored and loaded using HDF5 files.
We start from generating field-free rotational energies and matrix elements of dipole moment and polarizability for camphor molecule. For details, see “Rotational dynamics Quickstart”:
[1]:
from richmol.rot import solve, Solution
from richmol.rot import LabTensor
from richmol.rot import Molecule
from richmol.field import CarTens
WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
[2]:
camphor = Molecule()
camphor.XYZ = ("angstrom",
"O", -2.547204, 0.187936, -0.213755,
"C", -1.382858, -0.147379, -0.229486,
"C", -0.230760, 0.488337, 0.565230,
"C", -0.768352, -1.287324, -1.044279,
"C", -0.563049, 1.864528, 1.124041,
"C", 0.716269, -1.203805, -0.624360,
"C", 0.929548, 0.325749, -0.438982,
"C", 0.080929, -0.594841, 1.638832,
"C", 0.791379, -1.728570, 0.829268,
"C", 2.305990, 0.692768, 0.129924,
"C", 0.730586, 1.139634, -1.733020,
"H", -1.449798, 1.804649, 1.756791,
"H", -0.781306, 2.571791, 0.321167,
"H", 0.263569, 2.255213, 1.719313,
"H", 1.413749, -1.684160, -1.316904,
"H", -0.928638, -1.106018, -2.110152,
"H", -1.245108, -2.239900, -0.799431,
"H", 1.816886, -1.883799, 1.170885,
"H", 0.276292, -2.687598, 0.915376,
"H", -0.817893, -0.939327, 2.156614,
"H", 0.738119, -0.159990, 2.396232,
"H", 3.085409, 0.421803, -0.586828,
"H", 2.371705, 1.769892, 0.297106,
"H", 2.531884, 0.195217, 1.071909,
"H", 0.890539, 2.201894, -1.536852,
"H", 1.455250, 0.830868, -2.487875,
"H", -0.267696, 1.035608, -2.160680)
camphor.dip = [1.21615, -0.30746, 0.01140]
camphor.pol = [[115.80434, -0.58739, 0.03276], \
[-0.58739, 112.28245, 1.36146], \
[0.03276, 1.36146, 108.47809]]
camphor.frame = "diag(pol)"
sol = solve(camphor, Jmin=0, Jmax=10)
dip = LabTensor(camphor.dip, sol)
pol = LabTensor(camphor.pol, sol)
Here is how we can store and read Molecule object
[3]:
# store molecule camphor
camphor.store('camphor.h5', comment="S-Camphor structure from the Supersonic expansion FTMW spectra, Kisiel, et al., PCCP 5, 820 (2003)", replace=True)
# read molecule into camphor2
camphor2 = Molecule()
camphor2.read('camphor.h5')
# obtain rotational solution for camphor2
sol2 = solve(camphor2, Jmin=0, Jmax=10)
# compare solutions for camphor and camphor2
print("\ncompare solutions")
for J in sol.keys():
for sym in sol[J].keys():
for i in range(sol[J][sym].nstates):
enr = sol[J][sym].enr[i]
enr2 = sol2[J][sym].enr[i]
diff = enr - enr2
print(J, "%4s"%sym, i, "%12.6f"%enr, "%12.6f"%enr2, "%12.6f"%diff)
compare solutions
0 A 0 0.000000 0.000000 0.000000
1 A 0 0.076070 0.076070 0.000000
1 A 1 0.084863 0.084863 0.000000
1 A 2 0.087741 0.087741 0.000000
2 A 0 0.227612 0.227612 0.000000
2 A 1 0.234126 0.234126 0.000000
2 A 2 0.242759 0.242759 0.000000
2 A 3 0.269138 0.269138 0.000000
2 A 4 0.269736 0.269736 0.000000
3 A 0 0.453583 0.453583 0.000000
3 A 1 0.457678 0.457678 0.000000
3 A 2 0.474864 0.474864 0.000000
3 A 3 0.497349 0.497349 0.000000
3 A 4 0.500187 0.500187 0.000000
3 A 5 0.548850 0.548850 0.000000
3 A 6 0.548929 0.548929 0.000000
4 A 0 0.752988 0.752988 0.000000
4 A 1 0.755209 0.755209 0.000000
4 A 2 0.783441 0.783441 0.000000
4 A 3 0.801159 0.801159 0.000000
4 A 4 0.808865 0.808865 0.000000
4 A 5 0.854125 0.854125 0.000000
4 A 6 0.854669 0.854669 0.000000
4 A 7 0.924882 0.924882 0.000000
4 A 8 0.924891 0.924891 0.000000
5 A 0 1.125354 1.125354 0.000000
5 A 1 1.126433 1.126433 0.000000
5 A 2 1.167533 1.167533 0.000000
5 A 3 1.180186 1.180186 0.000000
5 A 4 1.195808 1.195808 0.000000
5 A 5 1.235838 1.235838 0.000000
5 A 6 1.237902 1.237902 0.000000
5 A 7 1.306557 1.306557 0.000000
5 A 8 1.306635 1.306635 0.000000
5 A 9 1.397421 1.397421 0.000000
5 A 10 1.397421 1.397421 0.000000
6 A 0 1.570646 1.570646 0.000000
6 A 1 1.571134 1.571134 0.000000
6 A 2 1.625927 1.625927 0.000000
6 A 3 1.633993 1.633993 0.000000
6 A 4 1.660444 1.660444 0.000000
6 A 5 1.693825 1.693825 0.000000
6 A 6 1.699453 1.699453 0.000000
6 A 7 1.764922 1.764922 0.000000
6 A 8 1.765299 1.765299 0.000000
6 A 9 1.855364 1.855364 0.000000
6 A 10 1.855374 1.855374 0.000000
6 A 11 1.966491 1.966491 0.000000
6 A 12 1.966492 1.966492 0.000000
7 A 0 2.088963 2.088963 0.000000
7 A 1 2.089172 2.089172 0.000000
7 A 2 2.157544 2.157544 0.000000
7 A 3 2.162142 2.162142 0.000000
7 A 4 2.201829 2.201829 0.000000
7 A 5 2.227777 2.227777 0.000000
7 A 6 2.239923 2.239923 0.000000
7 A 7 2.300048 2.300048 0.000000
7 A 8 2.301363 2.301363 0.000000
7 A 9 2.390006 2.390006 0.000000
7 A 10 2.390062 2.390062 0.000000
7 A 11 2.500692 2.500692 0.000000
7 A 12 2.500693 2.500693 0.000000
7 A 13 2.632097 2.632097 0.000000
7 A 14 2.632097 2.632097 0.000000
8 A 0 2.680384 2.680384 0.000000
8 A 1 2.680471 2.680471 0.000000
8 A 2 2.761852 2.761852 0.000000
8 A 3 2.764240 2.764240 0.000000
8 A 4 2.818738 2.818738 0.000000
8 A 5 2.837268 2.837268 0.000000
8 A 6 2.859247 2.859247 0.000000
8 A 7 2.911899 2.911899 0.000000
8 A 8 2.915538 2.915538 0.000000
8 A 9 3.001493 3.001493 0.000000
8 A 10 3.001727 3.001727 0.000000
8 A 11 3.111556 3.111556 0.000000
8 A 12 3.111564 3.111564 0.000000
8 A 13 3.242557 3.242557 0.000000
8 A 14 3.242557 3.242557 0.000000
8 A 15 3.394236 3.394236 0.000000
8 A 16 3.394236 3.394236 0.000000
9 A 0 3.344956 3.344956 0.000000
9 A 1 3.344991 3.344991 0.000000
9 A 2 3.438826 3.438826 0.000000
9 A 3 3.439984 3.439984 0.000000
9 A 4 3.509795 3.509795 0.000000
9 A 5 3.521809 3.521809 0.000000
9 A 6 3.556686 3.556686 0.000000
9 A 7 3.600292 3.600292 0.000000
9 A 8 3.608650 3.608650 0.000000
9 A 9 3.689945 3.689945 0.000000
9 A 10 3.690724 3.690724 0.000000
9 A 11 3.799224 3.799224 0.000000
9 A 12 3.799259 3.799259 0.000000
9 A 13 3.929651 3.929651 0.000000
9 A 14 3.929652 3.929652 0.000000
9 A 15 4.080960 4.080960 0.000000
9 A 16 4.080960 4.080960 0.000000
9 A 17 4.252910 4.252910 0.000000
9 A 18 4.252910 4.252910 0.000000
10 A 0 4.082701 4.082701 0.000000
10 A 1 4.082715 4.082715 0.000000
10 A 2 4.188630 4.188630 0.000000
10 A 3 4.189164 4.189164 0.000000
10 A 4 4.273832 4.273832 0.000000
10 A 5 4.280902 4.280902 0.000000
10 A 6 4.331072 4.331072 0.000000
10 A 7 4.364892 4.364892 0.000000
10 A 8 4.381250 4.381250 0.000000
10 A 9 4.455414 4.455414 0.000000
10 A 10 4.457593 4.457593 0.000000
10 A 11 4.563844 4.563844 0.000000
10 A 12 4.563980 4.563980 0.000000
10 A 13 4.693498 4.693498 0.000000
10 A 14 4.693502 4.693502 0.000000
10 A 15 4.844296 4.844296 0.000000
10 A 16 4.844296 4.844296 0.000000
10 A 17 5.015900 5.015900 0.000000
10 A 18 5.015900 5.015900 0.000000
10 A 19 5.208118 5.208118 0.000000
10 A 20 5.208118 5.208118 0.000000
Here is how we can store and read rotational solutions, produced by solve function
[4]:
# store solutions sol
sol.store('camphor.h5', replace=True)
# read solutions into sol2
sol2 = Solution()
sol2.read('camphor.h5')
# compute matrix elements of laboratory-frame tensors using sol2 solutions
dip2 = LabTensor(camphor.dip, sol2)
pol2 = LabTensor(camphor.pol, sol2)
# compare tensor matrix elements obtained with sol and sol2
# full matrix representations of dip and dip2 for x, y, and z components
dip_mat = [dip.tomat(form='full', repres='dense', cart=cart) for cart in 'xyz']
dip2_mat = [dip2.tomat(form='full', repres='dense', cart=cart) for cart in 'xyz']
# full matrix representations of pol and pol2 for xx, xy,, xz, yx, ... etc components
pol_mat = [pol.tomat(form='full', repres='dense', cart=cart+cart2) for cart in 'xyz' for cart2 in 'xyz']
pol2_mat = [pol2.tomat(form='full', repres='dense', cart=cart+cart2) for cart in 'xyz' for cart2 in 'xyz']
# compute and print maximal differences
import numpy
dip_diff = [numpy.max(abs(m-m2)) for m,m2 in zip(dip_mat, dip2_mat)]
pol_diff = [numpy.max(abs(m-m2)) for m,m2 in zip(pol_mat, pol2_mat)]
print("max diffs sol-sol2 in dipole matrix elements:", dip_diff)
print("max diffs sol-sol2 in polarizability matrix elements:", pol_diff)
max diffs sol-sol2 in dipole matrix elements: [0.0, 0.0, 0.0]
max diffs sol-sol2 in polarizability matrix elements: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
… and store/read Cartesian tensor operators
[5]:
# store dip and pol matrix elements
dip.store('camphor.h5', replace=True)
pol.store('camphor.h5', replace=True)
# read matrix elements into dip2 and pol2
dip2 = CarTens(filename='camphor.h5', name='dip')
pol2 = CarTens(filename='camphor.h5', name='pol')
# compare matrix elements, calculated and stored/read form file
# full matrix representations of dip and dip2 for x, y, and z components
dip_mat = [dip.tomat(form='full', repres='dense', cart=cart) for cart in 'xyz']
dip2_mat = [dip2.tomat(form='full', repres='dense', cart=cart) for cart in 'xyz']
# full matrix representations of pol and pol2 for xx, xy,, xz, yx, ... etc components
pol_mat = [pol.tomat(form='full', repres='dense', cart=cart+cart2) for cart in 'xyz' for cart2 in 'xyz']
pol2_mat = [pol2.tomat(form='full', repres='dense', cart=cart+cart2) for cart in 'xyz' for cart2 in 'xyz']
# compute and print maximal differences
import numpy
dip_diff = [numpy.max(abs(m-m2)) for m,m2 in zip(dip_mat, dip2_mat)]
pol_diff = [numpy.max(abs(m-m2)) for m,m2 in zip(pol_mat, pol2_mat)]
print("max diffs dip-dip2:", dip_diff)
print("max diffs pol-pol2:", pol_diff)
max diffs dip-dip2: [0.0, 0.0, 0.0]
max diffs pol-pol2: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
The field-free solutions can also be transformed and stored in a tensor form
[6]:
# convert solutions sol2 into tensor form h and store
h = LabTensor(camphor, sol2)
h.store('camphor.h5', comment="S-Camphor field-free energies (in cm^-1)", replace=True)
# read solutions into h2
h2 = CarTens(filename='camphor.h5', name='h')
# compare matrix elements, calculated and stored/read form file
mat = h.tomat(form='full', repres='dense', cart='0')
mat2 = h2.tomat(form='full', repres='dense', cart='0')
diff = numpy.max(abs(mat-mat2))
print("max diffs:", diff)
max diffs: 0.0
[ ]: