Conical intersection optimization

MLatom optimizes minimum-energy conical intersections with the two-state penalty-constrained scheme, which needs only the energies and gradients of the two states and no nonadiabatic couplings. Please cite the paper describing this method:

• Benjamin G. Levine, Joshua D. Coe, Todd J. Martínez. Optimizing Conical Intersections without Derivative Coupling Vectors: Application to Multistate Multireference Second-Order Perturbation Theory (MS-CASPT2). J. Phys. Chem. B 2008, 112, 405–413. DOI: 10.1021/jp0761618.

The optimization can be run on top of any excited-state method or ML model available in MLatom, e.g. OMNI-P2x or an MS-ANI model.

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.

In [ ]:

# 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.

In [2]:

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.

In [3]:

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.

In [4]:

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.

In [5]:

%%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.

In [6]:

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.

In [7]:

# 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

In [8]:

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.

In [9]:

#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
''')

In [10]:

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.

In [11]:

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.

In [12]:

%%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

In [15]:

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.

In [14]:

# 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