Tutorial 5 – SHARC MD with NN Model

In Tutorial 3, we showed how to train a NN model for predicting energies, forces and nonadiabatic couplings. Here we demonstrate, how to use the trained model, when running surface hopping molecular dynamics simulations with SHARC 3.0.

Some Text about SHARC …

%source $SHARC/sharcvars.sh

UsageError: Line magic function `%source` not found.

import os, sys
import spainn
from spainn.interface.sharcinterface import SHARC_NN

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[2], line 2
      1 import os, sys
----> 2 import spainn
      3 from spainn.interface.sharcinterface import SHARC_NN

ModuleNotFoundError: No module named 'spainn'

Molecule out of Dataset

In the following we create an Atoms object from the symbols and positions of our target molecule.

# molecular symbols of target structure (same as database)
symbols = 'CNHHHH'
# positions of atoms in target structure (not included in database)
positions = np.array(
    [[ 0.0000,  0.0000,  0.0000 ],
     [ 2.4321,  0.0000,  0.0000 ],
     [-1.0111,  1.7951,  0.0000 ],
     [ 3.4373,  1.6202, -0.2566 ],
     [ 3.4373, -1.6202,  0.2566 ],
     [-1.0111, -1.7951,  0.0000 ]]
)
# create ase Atoms object
target_mol = ase.Atoms(symbols=symbols, positions=positions)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 4
      2 symbols = 'CNHHHH'
      3 # positions of atoms in target structure (not included in database)
----> 4 positions = np.array(
      5     [[ 0.0000,  0.0000,  0.0000 ],
      6      [ 2.4321,  0.0000,  0.0000 ],
      7      [-1.0111,  1.7951,  0.0000 ],
      8      [ 3.4373,  1.6202, -0.2566 ],
      9      [ 3.4373, -1.6202,  0.2566 ],
     10      [-1.0111, -1.7951,  0.0000 ]]
     11 )
     12 # create ase Atoms object
     13 target_mol = ase.Atoms(symbols=symbols, positions=positions)

NameError: name 'np' is not defined

Define the calculator

In this step, we define a calculator for predicting for the target geometry the properties. Therefore, we have to specify the path to the trained model (see Tutorial 3) as well as the neighbor list for building the representation.

from schnetpack.transform import MatScipyNeighborList
from spainn.interface import NacCalculator

# set up calculator for predictions
calc = NacCalculator(
    model_file=os.path.join(os.getcwd(), 'train', 'best_model_E_F_C'),
    neighbor_list=MatScipyNeighborList(cutoff=10.0)
)
# apply calculator to target molecule
target_mol.calc = calc

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[4], line 1
----> 1 from schnetpack.transform import MatScipyNeighborList
      2 from spainn.interface import NacCalculator
      4 # set up calculator for predictions

ModuleNotFoundError: No module named 'schnetpack'

In the next step, we get the predictions from our model, by calling the get_properties function for the set of properties, we are interested in (here energies and nonadiabatic couplings). Note, the NacCalculator always returns nonadiabatic couplings \(\mathbf{C}_{ij}\) rather than the smoothed nonadiabatic couplings \(\mathbf{C}_{ij}^s\) employed in the training. Nevertheless, the predictions, need to be requested with the same property key as used in the training.

# predict energies and nacs from NN model
pred = target_mol.get_properties(['energy', 'smooth_nacs'])
pred

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 2
      1 # predict energies and nacs from NN model
----> 2 pred = target_mol.get_properties(['energy', 'smooth_nacs'])
      3 pred

NameError: name 'target_mol' is not defined

Mean absolute error for energies

# reference energies from (6,4)-CASSCF/cc-pVDZ calculations
ref = {'energy': np.array([-94.69603647, -94.43388056, -94.40410243])}

mae = np.mean(
    np.abs(np.asarray(pred['energy']).flatten()
           - np.asarray(ref['energy']).flatten()
          ))

print('MAE: '+str(mae))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 2
      1 # reference energies from (6,4)-CASSCF/cc-pVDZ calculations
----> 2 ref = {'energy': np.array([-94.69603647, -94.43388056, -94.40410243])}
      4 mae = np.mean(
      5     np.abs(np.asarray(pred['energy']).flatten()
      6            - np.asarray(ref['energy']).flatten()
      7           ))
      9 print('MAE: '+str(mae))

NameError: name 'np' is not defined

Visualization of prediction quality for energies

# define plot
plt.figure(figsize=(5, 5))
plt.scatter(ref['energy'], pred['energy'], s=100.0)
plt.rcParams["figure.autolayout"] = True
plt.rcParams["font.size"] = 14

# get limits of axis and aplly same for x and y-axis (square plot)
limits = [*plt.xlim(), *plt.ylim()]
lower_lim = min(limits)
upper_lim = max(limits)
plt.xlim(lower_lim, upper_lim)
plt.ylim(lower_lim, upper_lim)

# plot y=x line to guide eye
plt.plot([lower_lim-2.5, upper_lim+2.5], [lower_lim-2.5, upper_lim+2.5], c='k', alpha=0.2, linewidth=8, linestyle='-')

plt.xlabel('E$_i$(ref) / Ha')
plt.ylabel('E$_i$(NN) / Ha')
plt.title('MAE: '+str(round(mae, 4)))

plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 2
      1 # define plot
----> 2 plt.figure(figsize=(5, 5))
      3 plt.scatter(ref['energy'], pred['energy'], s=100.0)
      4 plt.rcParams["figure.autolayout"] = True

NameError: name 'plt' is not defined