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],
  32.                            requires_grad=False, device=device, dtype=torch.float32)
  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],
  61.                            requires_grad=False, device=device, dtype=torch.float32)
  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. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top