• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Jun 27th, 2019 70 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import torchani
2. import torch
3. from simtk import unit
4. import numpy as np
5. import matplotlib.pyplot as plt
6.
7. hartree_to_kJ_mol = 2625.499638
8. nm_to_angstroms = (1.0 * unit.nanometer) / (1.0 * unit.angstrom)
9.
10. device = torch.device('cpu')
11. model = torchani.models.ANI1ccx()
12. species = model.species_to_tensor('CHHHH').to(device).unsqueeze(0)
13.
14. def ANI1ccx_energy(x):
15.     """
16.     Parameters
17.     ----------
18.     x : numpy array, unit'd
19.         coordinates
20.
21.     Returns
22.     -------
23.     energy : float, unit'd
24.         energy, with unit of kJ/mol attached
25.     """
26.
27.     assert(type(x) == unit.Quantity)
28.
29.     x_in_nm = x.value_in_unit(unit.nanometer)
30.
31.     coordinates = torch.tensor([x_in_nm],
33.
34.     # convert from nm to angstroms
35.     coordinates_in_angstroms = coordinates * nm_to_angstroms
36.     _, energy_in_hartree = model((species, coordinates_in_angstroms))
37.
38.     # convert energy from hartrees to kJ/mol
39.     energy_in_kJ_mol = energy_in_hartree * hartree_to_kJ_mol
40.     return float(energy_in_kJ_mol.detach()) * unit.kilojoule_per_mole
41.
42.
43. def ANI1ccx_energy_ensemble(x):
44.     """
45.     Parameters
46.     ----------
47.     x : numpy array, unit'd
48.         coordinates
49.
50.     Returns
51.     -------
52.     energies : numpy array, unit'd
53.         8 energies from the ensemble, with units of kJ/mol attached
54.     """
55.
56.     assert(type(x) == unit.Quantity)
57.
58.     x_in_nm = x.value_in_unit(unit.nanometer)
59.
60.     coordinates = torch.tensor([x_in_nm],
62.
63.     # convert from nm to angstroms
64.     coordinates_in_angstroms = coordinates * nm_to_angstroms
65.
66.     energies_in_kJ_mol = [float(hartree_to_kJ_mol * m((species, coordinates_in_angstroms))[1].detach()) for m in model]
67.
68.     return np.array(energies_in_kJ_mol) * unit.kilojoule_per_mole
69.
70. scale_factors = np.linspace(0,2,1000)
71.
72. tetrahedral_geometry = np.array([
73.     (0, 0, 0),
74.     (1, 1, 1),
75.     (-1, -1, 1),
76.     (-1, 1, -1),
77.     (1, -1, -1),
78. ], dtype=float)
79. tetrahedral_geometry /= np.linalg.norm(tetrahedral_geometry[0] - tetrahedral_geometry[1])
80.
81. # lowest energy found
82. bent_optimum = np.array([
83.     [ 0.        ,  0.        ,  0.        ],
84.     [ 0.64702535,  0.64647233,  0.8811238 ],
85.     [-0.64734113, -0.6475463 ,  0.8808765 ],
86.     [-0.7168874 ,  0.7172707 , -0.33484012],
87.     [ 0.7182246 , -0.71852016, -0.33497408]
88. ]) # in angstroms
89. lowest_energy = ANI1ccx_energy(bent_optimum * unit.angstroms)
90.
91. # scaled so that the average bond length is 1 angstrom
92. normalized_bent_optimum = bent_optimum / np.mean(np.linalg.norm(bent_optimum[1:], axis=1))
93.
94.
95. # energy scan from global optimum (bent geometry), with ensemble uncertainty
96. bent_scan = np.array([normalized_bent_optimum * scale_factor for scale_factor in scale_factors])
97. bent_energies = [ANI1ccx_energy_ensemble(g * unit.angstrom) for g in bent_scan]
98. bent_energy_array = np.array([(e - lowest_energy) / unit.kilojoule_per_mole for e in bent_energies])
99. bent_means = np.mean(bent_energy_array, 1)
100. bent_stdevs = np.std(bent_energy_array, 1)
101.
102.
103. # energy scan from tetrahedral geometry, with ensemble uncertainty
104. tetrahedral_scan = np.array([tetrahedral_geometry * scale_factor for scale_factor in scale_factors])
105. tetra_energies = [ANI1ccx_energy_ensemble(g * unit.angstrom) for g in tetrahedral_scan]
106. tetra_energy_array = np.array([(e - lowest_energy) / unit.kilojoule_per_mole for e in tetra_energies])
107. tetra_means = np.mean(tetra_energy_array, 1)
108. tetra_stdevs = np.std(tetra_energy_array, 1)
109.
110.
111. def average_bond_length_in_snapshot(xyz):
112.     return np.mean(np.linalg.norm(xyz[1:], axis=1))
113.
114.
115. # plot
116. plt.figure()
117. ax = plt.subplot(1,1,1)
118. ax.spines['top'].set_visible(False)
119. ax.spines['right'].set_visible(False)
120.
121. plt.ylim(-20,400)
122. plt.xlim(0.85,1.35)
123.
124. plt.xlabel('C-H bond length (Å)')
125. plt.ylabel('$\Delta E \pm$ stddev (kJ/mol)\nrelative to global minimum')
126. plt.title('methane energy scans: varying C-H bond length')
127.
128. # tetrahedral
129. # x: average C-H bond length
130. tetrahedral_x = [average_bond_length_in_snapshot(xyz) for xyz in tetrahedral_scan]
131. # y: mean energy +/- standard deviation
132. plt.plot(tetrahedral_x, tetra_means, label='scan from tetrahedron')
133. plt.fill_between(tetrahedral_x, tetra_means - tetra_stdevs, tetra_means + tetra_stdevs, alpha=0.5)
134.
135. # bent
136. # x: average C-H bond length
137. bent_x = [average_bond_length_in_snapshot(xyz) for xyz in bent_scan]
138. # y: mean energy +/- standard deviation
139. plt.plot(bent_x, bent_means, label='scan from global minimum (bent)')
140. plt.fill_between(bent_x, bent_means - bent_stdevs, bent_means + bent_stdevs, alpha=0.5)
141.
142. # horizontal line indicating global minimum
143. plt.hlines(0, 0,2,
144.            linestyles='dotted', linewidth=1, label='global minimum (bent)')
145.
146.
147. # clip x and y limits
148. plt.ylim(-150,400)
149. plt.xlim(0.85,1.35)
150.
151.
152. plt.legend(loc='lower left')
153. plt.tight_layout()
154.
155. plt.savefig('methane_energy_scans_with_uncertainties.png', dpi=300, bbox_inches='tight')
156. plt.close()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.
Not a member of Pastebin yet?