Conical intersection optimization
Download tutorial Jupyter notebook with required files.
Conical Intersection Optimization with MLatom¶
By Mikolaj Martyka, 6 November 2025.
This short tutorial demonstrates how to perform conical-intersection (CI) optimizations using the ci_opt class in MLatom.
We’ll walk through two examples:
A small π-system (fulvene) optimized with the MS-ANI neural network model.
Cyclohexadiene (CHD), optimized with a semi-empirical ODM2 model.
The CI optimization method implemented here follows the penalty-constraint approach of Levine et al. (J. Phys. Chem. B 112, 405–413 (2008)), where the geometry is optimized to minimize the average energy of two electronic states while penalizing the energy gap between them.
This method is useful for locating conical intersections with models that cannot easily predict nonadiabatic couplings (NACs), although more sophisticated algorithms exist for higher accuracy.
# First, lets import MLatom. Make sure you have version >= 3.19.1
import mlatom as ml
Optimization of CI of fulvene with MS-ANI model¶
Next, we create a simple molecule, fulvene, directly from a SMILES string. This is convenient for quick demonstrations without needing an input file.
fulvene = ml.molecule.from_smiles_string("C=C1\C=C/C=C1")
We can visualize the generated molecular structure to confirm that it was read correctly.
fulvene.view()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
Now we initialize an MS-ANI model. This is a multistate version of the ANI neural network potential, capable of describing excited states. The model file is loaded from disk, and its DOI (10.1038/s41524-025-01636-z) provides details on the architecture and training.
mlmodel = ml.models.msani(model_file='fulvene_msani.pt')
model loaded from fulvene_msani.pt
With the model ready, we can perform a conical-intersection optimization using ml.ci_opt. This will adjust the geometry to minimize the average energy between two electronic states while applying a penalty to the energy gap — following the method of Levine et al. (J. Phys. Chem. B 112, 405–413 (2008)). The result should be a structure near the S₀/S₁ intersection seam.
%%time
_ = ml.ci_opt(molecule=fulvene, model=mlmodel)
CPU times: user 12.1 s, sys: 343 ms, total: 12.4 s Wall time: 37.4 s
After optimization, let’s visualize the optimized geometry. Comparing it to the starting structure can reveal how the molecule distorts toward the intersection, by twisting the C=CH2 bond to be more perpendicular.
print(f'Dihedral angle: {fulvene.dihedral_angle(6, 0, 1, 2):.2f} degree')
fulvene.view()
Dihedral angle: 59.17 degree
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
We can inspect the energy gaps between the lowest two electronic states. A small S₀–S₁ gap indicates that the optimizer successfully approached the conical intersection.
# observe the small energy gap between S0 and S1
print(f'Energy gap between S0 and S1: {fulvene.energy_gaps[0][1]*ml.constants.hartree2eV:.2f} eV')
Energy gap between S0 and S1: 0.04 eV
print(fulvene.get_xyz_string())
12 C 2.2261782498407 0.0288249015972 -0.0224619628967 C 0.7534574980195 0.0023262310261 -0.0178462961598 C -0.1048341377048 1.1136448374900 -0.0985504717530 C -1.4814337972701 0.6433604187912 -0.0699358298845 C -1.4564188747425 -0.7199309909833 0.0461200770069 C -0.0634519070593 -1.1395907859531 0.0668748382233 H 2.7698346864724 0.5658982042672 0.7270175692633 H 2.7841724590479 -0.4876574080361 -0.7758438945271 H 0.2192085967770 2.1270107716469 -0.1991350439951 H -2.3521232761106 1.2677188871543 -0.1125228801202 H -2.3035671275105 -1.3755162489931 0.0936247706131 H 0.2979994192141 -2.1403643284423 0.1656517896772
Optimization of cyclohexadiene¶
Let’s move to a more complex system — cyclohexadiene (CHD). Here, we load the geometry from an .xyz file representing an S₂→S₁ hopping point from a nonadiabatic dynamics trajectory. For this algorithm, providing a good starting guess is very important.
#load cyclohexadiene. It is important to make a reasonable starting guess in this case, so we take an S2->S1 hopping geometry from NAMD.
chd = ml.molecule.from_xyz_string('''14
C 1.3352823257446 -0.7696992754936 -0.0487115159631
C 1.2886266708374 0.7333662509918 0.0120503930375
C 0.0400483943522 1.4208787679672 0.0668141320348
C -1.1751360893250 0.8934385776520 -0.3302345275879
C -1.1844748258591 -0.8462789058685 0.3520835340023
C 0.0853536650538 -1.4187078475952 -0.1039694100618
H 2.2638611793518 1.2157090902328 0.1544666290283
H -0.2900925278664 2.2322068214417 0.7879427075386
H -1.2981988191605 0.6020563244820 -1.3704229593277
H -1.9767922163010 -1.4304752349854 -0.1777630448341
H -0.1033942699432 -2.3882086277008 -0.6410480141640
H 2.3468582630157 -1.2727781534195 0.1109914630651
H -2.0999982357025 1.2376796007156 0.2713582217693
H -1.4057227373123 -0.4323001801968 1.4544125795364
''')
print(f'Distance between atoms 4 and 5: {chd.bond_length(3, 4):.3f} Angstrom')
chd.view()
Distance between atoms 4 and 5: 1.869 Angstrom
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
MECI's can also be optimized with quantum chemical methods, here we will use the semi-empirical ODM2 method, which provides a good description of excited states, for a relatively low cost. We will need to use the MRCI variant, which is defined in the keywords file mndokw.
odm2 = ml.methods(method='odm2', read_keywords_from_file='mndokw') # load the odm2 method
As mentioned, we will optimize the S2/S1 MECI, which is of physical importance in this system.
%%time
_ = ml.ci_opt( molecule=chd, model=odm2, state_0=1, state_1=2, maximum_number_of_steps=400) # conical intersection between higher excited states can also be optimized
CPU times: user 20.4 s, sys: 272 ms, total: 20.6 s Wall time: 1min
print(f'Distance between atoms 4 and 5: {chd.bond_length(3, 4):.3f} Angstrom')
chd.view()
Distance between atoms 4 and 5: 1.907 Angstrom
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
# observe the small energy gap between S1 and S2
print(f'Energy gap between S1 and S2: {chd.energy_gaps[1][2]*ml.constants.hartree2eV:.2f} eV')
Energy gap between S1 and S2: 0.04 eV