Surface-hopping dynamics

In this tutorial we show how to perform nonadiabatic molecular dynamics (NAMD) and analysis of NAMD trajectories with MLatom. The simulations are only possible through the Python API. MLatom currently only supports NAC-free Landau–Zener–Belyev–Lebedev (LZBL) surface hopping.

You can run the TSH with MLatom for the following models:

  • CASSCF through the interface to COLUMBUS

  • ADC(2) through the interface to Turbomole

  • AIQM1/CI

  • MS-ANI

  • Any ML models that provides energies and forces for electronic states of interset.

Running NAMD dynamics

See our paper for more details (please also cite it if you use the corresponding features):

ML-NAMD with AIQM1/MRCI

Note

Please refer to the tutorials of how to use ML for NAMD with MS-ANI and do active learning. The tutorial shown here will need some installations required to perform AIQM1 calculations, while NAMD with ML models will need minimum installations and can be performed online.

Here we show how to use AIQM1/MRCI in propagating LZBL surface-hopping dynamics. Please download the tutorial files (namd_aiqm1.zip).

The script below does all the steps from start to finish:

  • optimizes geometry

  • runs frequency calculations

  • samples initial conditions from the Wigner distribution using Newton-X’s routines. Example can be extended to supports filtering by excitation energy window (refer to the manual for more details)

  • propagates multiple trajectories in parallel (here 16 trajectories for 5 fs with 0.1 fs time step)

  • saves trajectories in the h5md format

  • analyzes results by generating the population plots.

Here is the script:

import mlatom as ml

# Load the initial geometry of a molecule
mol = ml.data.molecule()
mol.charge=1
mol.read_from_xyz_file('cnh4+.xyz')

# Define methods
# .. for NAMD
aiqm1 = ml.models.methods(method='AIQM1',
                        qm_program_kwargs={'save_files_in_current_directory': True,
                                            'read_keywords_from_file':'../materials/mndokw'})
# .. for optimization, frequencies and normal mode calculations
method_optfreq = ml.models.methods(method='B3LYP/Def2SVP', program='pyscf')

# Optimize geometry
geomopt = ml.simulations.optimize_geometry(model=method_optfreq,
                                        initial_molecule=mol)
eqmol = geomopt.optimized_molecule
eqmol.write_file_with_xyz_coordinates('eq.xyz')

# Get frequencies
ml.simulations.freq(model=method_optfreq,
                    molecule=eqmol)
eqmol.dump(filename='eqmol.json', format='json')

# Get initial conditions
init_cond_db = ml.generate_initial_conditions(molecule=eqmol,
                                    generation_method='wigner',
                                    number_of_initial_conditions=16,
                                    initial_temperature=0,
                                    random_seed=1) # To ensure we always get the same initial conditions (should not be used in actual calculations)
init_cond_db.dump('test.json','json')

# Propagate multiple LZBL surface-hopping trajectories in parallel
# .. setup dynamics calculations
namd_kwargs = {
            'model': aiqm1,
            'time_step': 0.25,
            'maximum_propagation_time': 5,
            'hopping_algorithm': 'LZBL',
            'nstates': 3,
            'initial_state': 2, # Numbering of states starts from 0!
            'random_seed': 1 # To ensure we always get the same initial conditions (should not be used in actual calculations)
            }

# .. run trajectories in parallel
dyns = ml.simulations.run_in_parallel(molecular_database=init_cond_db,
                                      task=ml.namd.surface_hopping_md,
                                      task_kwargs=namd_kwargs,
                                      create_and_keep_temp_directories=True)
trajs = [d.molecular_trajectory for d in dyns]

# Dump the trajectories
itraj=0
for traj in trajs:
    itraj+=1
    traj.dump(filename=f"traj{itraj}.h5",format='h5md')

# Analyze the result of trajectories and make the population plot
ml.namd.analyze_trajs(trajectories=trajs, maximum_propagation_time=5)
ml.namd.plot_population(trajectories=trajs, time_step=0.25,
                        max_propagation_time=5, nstates=3, filename=f'pop.png', pop_filename='pop.txt')

And here is the final population plot (your plot will be different because of the random seed in initial conditions and hoppings):

_images/cnh4%2B_aiqm1cis_lznamd_population.png

You will also get the text file with populations pop.txt which should look like:

0.000 0.0 0.0 1.0
0.250 0.0 0.0 1.0
0.500 0.0 0.0 1.0
0.750 0.0 0.0 1.0
...

Download the full file cnh4+_aiqm1cis_lznamd_population.txt.

Multi-state ANI models

Multi-state learning model (MS-ANI) that has unrivaled accuracy for excited state properties (accuracy is often better than for models targeting only ground state!). We demonstrate that this model can be used for trajectory-surface hopping of multiple molecules (not just for a single molecule!) in:

  • Mikołaj Martyka, Lina Zhang, Fuchun Ge, Yi-Fan Hou, Joanna Jankowska, Mario Barbatti, Pavlo O. Dral. Charting electronic-state manifolds across molecules with multi-state learning and gap-driven dynamics via efficient and robust active learning. 2024. Preprint on ChemRxiv: https://doi.org/10.26434/chemrxiv-2024-dtc1w.

Zip with tutorial materials including Jupyter notebook:

msani

ML-NAMD with single-state ML models

In this tutorial, we show an example of running surface-hopping MD with single-state ML models. Please see a separate tutorial on machine learning potentials.

See our paper for more details (please also cite it if you use the corresponding features):

You can download the Jupyter notebook with the required initial conditions and ML models from this paper.

The tutorial calculations are very fast and you should be able to get the final population plot for 5 fs with 0.25 fs time step from 30 trajectories within a minute. Here is the Jupyter notebook code snippet:

import mlatom as ml
import os
import numpy as np

# Read initial conditions
init_cond_db = ml.data.molecular_database.load(filename='materials/init_cond_db_for_pyrazine.json', format='json')

# We need to create a class that accepts the specific arguments shown below and saves the calculated electronic state properties in the molecule object
class mlmodels():
    def __init__(self, nstates = 5):
        folder_with_models = 'materials/lz_models'
        self.models = [None for istate in range(nstates)]
        for istate in range(nstates):
            self.models[istate] = [ml.models.ani(model_file=f'{folder_with_models}/ensemble{ii+1}s{istate}.pt') for ii in range(2)]
            for ii in range(2): self.models[istate][ii].nthreads = 1

    def predict(self,
            molecule=None,
            nstates=5,
            current_state=0,
            calculate_energy=True,
            calculate_energy_gradients=True):

        molecule.electronic_states = [molecule.copy() for ii in range(nstates)]

        for istate in range(nstates):
            moltmp = molecule.electronic_states[istate]
            moltmpens = [moltmp.copy() for ii in range(2)]
            for ii in range(2):
                self.models[istate][ii].predict(molecule=moltmpens[ii], calculate_energy = True, calculate_energy_gradients = True)
            moltmp.energy = np.mean([moltmpens[ii].energy for ii in range(2)])
            moltmp.energy_gradients = np.mean([moltmpens[ii].energy_gradients for ii in range(2)], axis=0)

        molecule.energy = molecule.electronic_states[current_state].energy
        molecule.energy_gradients = molecule.electronic_states[current_state].energy_gradients

models = mlmodels()

# Arguments for running NAMD trajectories
timemax = 5 # fs
namd_kwargs = {
            'model': models,
            'time_step': 0.25, # fs
            'maximum_propagation_time': timemax,
            'dump_trajectory_interval': None,
            'hopping_algorithm': 'LZBL',
            'nstates': 5,
            'random_seed': 1, # making sure that the hopping probabilities are the same (should not be used in actual calculations!)
            'rescale_velocity_direction': 'along velocities',
            'reduce_kinetic_energy': False,
            }

# Run 30 trajectories
dyns = ml.simulations.run_in_parallel(molecular_database=init_cond_db[:30], task=ml.namd.surface_hopping_md, task_kwargs=namd_kwargs)
trajs = [d.molecular_trajectory for d in dyns]
ml.namd.analyze_trajs(trajectories=trajs, maximum_propagation_time=timemax)

# Dump the trajectories
for itraj in range(len(trajs)):
    trajs[itraj].dump(filename=f'traj{itraj+1}.json', format='json')

# Prepare the population plot
ml.namd.plot_population(trajectories=trajs, time_step=0.25,
                        max_propagation_time=timemax, nstates=5, filename=f'pop.png')

Since we used the fixed random seed, you should get the following final population:

_images/pyrazine_lznamd_ml_population.png

Analyzing results

data_analysis.ipynb.

data_analysis