面跳跃动力学
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):
Lina Zhang, Sebastian Pios, Mikołaj Martyka, Fuchun Ge, Yi-Fan Hou, Yuxinxin Chen, Joanna Jankowska, Lipeng Chen, Mario Barbatti, Pavlo O. Dral. MLatom software ecosystem for surface hopping dynamics in Python with quantum mechanical and machine learning methods. J. Chem. Theory Comput. 2024, 20, 5043–5057. DOI: 10.1021/acs.jctc.4c00468. Preprint on arXiv: https://arxiv.org/abs/2404.06189.
ML-NAMD with AIQM1/MRCI
备注
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
).
以下脚本展示了所有步骤:
优化几何构型
运行频率计算
使用Newton-X的例程从Wigner分布中生成初始条件。该示例可以扩展以支持按激发能量窗口进行筛选(有关更多详情,请参阅 手册 )。
并行传播多条轨迹(此处为16条轨迹,每条轨迹为5 fs,时间步长为0.1 fs)。
将轨迹储存为h5md格式
通过生成群体图来分析结果。
脚本如下:
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')
这是最终的群体图(由于初始条件和跃迁中的随机种子不同,您的图可能会有所不同):

您还将获得带有群体信息的文本文件 pop.txt
,其内容应如下:
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
...
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:
#Import MLatom as a Python package.
import mlatom as ml
#First, we define a MS-ANI model
model = ml.models.msani(model_file='fulvene_tutorial.npz',nstates=2)
the trained MS-ANI model will be saved in fulvene_tutorial.npz
#then, we load the training data as a MLatom molecular database. This is a small dataset for
#demo purposes.
train_data = ml.data.molecular_database.load("tutorial_data.json", format="json")
#Now we can train the model.
model.train(molecular_database=train_data,
property_to_learn='energy',
xyz_derivative_property_to_learn='energy_gradients',
hyperparameters={'max_epochs': 100}) #100 epochs is not enough, only for demo.
#The model can be used to make single-point predictions
#First we need to load a fulvene molecule from .xyz
mol = ml.data.molecule()
mol.read_from_xyz_file('fulvene.xyz')
mol.view()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
#Now we can make a prediction
model.predict(molecule=mol, nstates=2, current_state=0,
calculate_energy=True, calculate_energy_gradients=True)
#And print the results
for state in mol.electronic_states:
print(state.energy)
gap = (mol.electronic_states[1].energy-mol.electronic_states[0].energy)*ml.constants.Hartree2eV
print("The S1-S0 gap is {} eV".format(gap))
-230.56873760870891 -230.49484480298358 The S1-S0 gap is 2.010727281361726 eV
#Finally, we can use the MS-ANI model to run NAMD.
#For that we will load a model trained in the AL loop.
model_namd = ml.models.msani(model_file='fulvene_energy_iteration_19.npz')
model loaded from fulvene_energy_iteration_19.npz
#Load an initial conditions database
init_db = ml.data.molecular_database.load("init_cond_db.json", format="json")
#And define NAMD arguments
timemax = 60 # fs
namd_kwargs = {
'model': model_namd,
'time_step': 0.1, # fs
'maximum_propagation_time': timemax,
'hopping_algorithm': 'LZBL',
'nstates': 2,
'reduce_kinetic_energy': True,
}
#And finally run the NAMD. This should take about a minute.
dyns = ml.simulations.run_in_parallel(molecular_database=init_db[:4], task=ml.namd.surface_hopping_md, task_kwargs=namd_kwargs, create_and_keep_temp_directories=False)
#Now we can plot the results
trajs = [d.molecular_trajectory for d in dyns]
ml.namd.plot_population(trajectories=trajs, time_step=0.1,
max_propagation_time=timemax, nstates=2)
#MLatom trajectories can be analyzed and visualized using simple Python code. An example below.
energies = [[],[]]
time = []
gap = []
for step in trajs[0].steps:
energies[0].append(step.molecule.electronic_states[0].energy)
energies[1].append(step.molecule.electronic_states[1].energy)
time.append(step.time)
gap.append(step.molecule.electronic_states[1].energy-step.molecule.electronic_states[0].energy)
import matplotlib.pyplot as plt
import numpy as np
plt.plot(time, np.transpose(energies))
plt.xlabel("Time (fs)")
plt.ylabel("Energy (Hartree)")
plt.legend(['S0', 'S1'])
<matplotlib.legend.Legend at 0x2adbeae2d940>
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):
Lina Zhang, Sebastian Pios, Mikołaj Martyka, Fuchun Ge, Yi-Fan Hou, Yuxinxin Chen, Joanna Jankowska, Lipeng Chen, Mario Barbatti, Pavlo O. Dral. MLatom software ecosystem for surface hopping dynamics in Python with quantum mechanical and machine learning methods. J. Chem. Theory Comput. 2024, 20, 5043–5057. DOI: 10.1021/acs.jctc.4c00468. Preprint on arXiv: https://arxiv.org/abs/2404.06189.
您可以从本文 下载
具有所需初始条件和ML模型的Jupyter笔记本。
教程中的计算速度非常快,您应该能够在一分钟内从30条轨迹中获得5 fs的最终群体图,时间步长为0.25 fs。以下是Jupyter笔记本的代码片段:
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')
由于我们使用了固定的随机种子,您应该得到以下最终的群体分布:

Analyzing results
Author: Mikołaj Martyka
This tutorial shows how to analyse results of NAMD simulations and extract data from an ensamble of trajectories. The data used in this tutorial is available at this DOI: 10.6084/m9.figshare.28044863 Taken from our paper on active learning for TSH simulations, available at: https://chemrxiv.org/engage/chemrxiv/article-details/66e7e01712ff75c3a177a270
import mlatom as ml
Example 1: Analysis of the hopping geometries in fulvene¶
As a first example we will analyse hopping geometries of fulvene, after excitation to S1. Reference, eg. doi.org/10.12688/openreseurope.13624.2
#Define a function for extracting hopping geometries.
def get_hop_geos(traj):
nstates = len(traj.steps[0].molecule.electronic_states)
local_hop_list = [[ [] for k in range(nstates)] for j in range(nstates)]
for index, geom in enumerate(traj.steps):
if index !=0:
if geom.current_state != traj.steps[index-1].current_state:
hop_to = int(geom.current_state)
hop_from = int(traj.steps[index-1].current_state)
traj.steps[index-1].molecule.hoptime = traj.steps[index-1].time
print
local_hop_list[hop_from][hop_to].append(traj.steps[index-1].molecule)
geoms = []
for k in range(nstates):
for j in range(nstates):
if local_hop_list[k][j] != None:
for idx, hop_geom in enumerate(local_hop_list[k][j]):
hop_geom.hop_from = k
hop_geom.hop_to = j
geoms.append(hop_geom)
return geoms
#Load a fulvene molecule in MLatom format
mol = ml.data.molecule.load("eqmol.json", format='json')
#Check what is the numbering of the atoms in the xyz file
mol.write_file_with_xyz_coordinates("fulvene.xyz")
mol.view()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
#Define lists for the relevant quantities
C_CH2_lens = [] #The bond length between atoms no. 0 and 5
#The four dihedral angles required for eq. (8)
C_CH2_dihedral_Cis1= []
C_CH2_dihedral_Cis2 = []
C_CH2_dihedral_trans1 = []
C_CH2_dihedral_trans2 = []
#List for hop times
hop_times = []
#Iterate over trajectories...
for i in range(50):
#...load them...
traj = ml.data.molecular_trajectory()
traj.load("traj{}.json".format(i+1), format='json')
#...extract hop geoms...
hop_geoms = get_hop_geos(traj)
#...and finally extract the data.
for mol in hop_geoms:
C_CH2_lens.append(mol.bond_length(0, 5))
C_CH2_dihedral_Cis1.append(mol.dihedral_angle(2, 0, 5,10))
C_CH2_dihedral_Cis2.append(mol.dihedral_angle(1, 0, 5,11))
C_CH2_dihedral_trans1.append(mol.dihedral_angle(2, 0, 5,11))
C_CH2_dihedral_trans2.append(mol.dihedral_angle(1, 0, 5,10))
hop_times.append(mol.hoptime)
#Calculate the mean dihedral angle according to Barbatti et al Open. Res. Eur. 2022, 1, 49.
hop_angle = []
for i in range(len(C_CH2_dihedral_Cis1)):
sum1 = abs(C_CH2_dihedral_Cis1[i])+abs(C_CH2_dihedral_Cis2[i])+abs(180-abs(C_CH2_dihedral_trans1[i]))+abs(180-abs(C_CH2_dihedral_trans2[i]))
hop_angle.append(0.25*sum1)
import matplotlib.pyplot as plt
#Plot the results, using the hopping times as a colormap. See fig 5 of 10.26434/chemrxiv-2024-dtc1w-v2 for a plot with more points and explanation.
plt.scatter(C_CH2_lens, hop_angle, c=hop_times)
cbar = plt.colorbar()
cbar.set_label('Time (fs)')
plt.xlabel("C=CH2 distance (A)")
plt.ylabel("Mean C=CH2 dihedral (degrees)")
Text(0, 0.5, 'Mean C=CH2 dihedral (degrees)')
Example 2: Analysis of the quantum yield of cis -> trans photoisomerization of azobenzene¶
Now we will calculate the quantum yield (QY) of azobenzene cis -> trans isomerisation. Reference: 10.26434/chemrxiv-2024-dtc1w-v2, fig. 7 and text.
#list for dihedral angles
dih_fin = []
#Iterate over trajectories
for i in range(50):
#Load the final geometries of each trajectory. Notice that the format and location is different, this is the typical output of MLatom's run_in_parallel.
traj = ml.data.molecular_trajectory()
traj.load("job_surface_hopping_md_{}/traj.h5".format(i+1), format='h5md')
#Take the final geometry in each trajectory...
finmol = traj.steps[-1].molecule
#...and the dihedral angle.
dih_fin.append(finmol.dihedral_angle(11, 23, 22,3))
#View one representative trajectory
traj.view()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
#Plot a histogram of the dihedral angles
plt.hist(dih_fin)
plt.ylabel("Count")
plt.xlabel("Dihedral angle (Deg.)")
Text(0.5, 0, 'Dihedral angle (Deg.)')
#Count the number of trans-azobenzenes at the end, taking the descriptor of dihedral >60 as trans-azobenzene.
counter = 0
for angle in dih_fin:
if angle > 60:
counter+=1
#Print the QY
QY = counter/50*100
print("The quantum yield is {}%!".format(QY))
The quantum yield is 68.0%!