Single-point calculations

MLatom supports single-point excited-state simulations for different electron structure methods and machine learning models. The most powerful way to use MLatom for such calculations is via Python API, however, some methods also support the calculations via input file.

Python API

Please download the tutorial files (sp_aiqm1.zip).

To start, we provide the initial geometry for \(\text{CNH}_4^+\) in the cnh4+.xyz file, with coordinates in Angstrom, which looks like this:

6
symmetry c1
N        0.051443000     -0.125748000      0.596619000
C        0.067113000     -0.025698000     -0.683445000
H        0.002169000      0.695516000      1.199263000
H        0.087711000     -1.030762000      1.065358000
H        0.027120000      0.954772000     -1.143351000
H        0.120118000     -0.922908000     -1.288953000

Then we can use this geometry to calculate the excited-state properties at the AIQM1/CIS level (in the comments also for AIQM1/MRCI):

#!/usr/bin/env python
# coding: utf-8

import mlatom as ml

# get the molecule object
mol = ml.data.molecule()
mol.read_from_xyz_file('cnh4+.xyz')
mol.charge = 1

# define the method
aiqm1 = ml.models.methods(method='AIQM1')
# To run AIQM1/MRCI calculations, you can uncomment the following line and make sure that mndokw file is available (see below)
#aiqm1 = ml.models.methods(method='AIQM1',qm_program_kwargs={'read_keywords_from_file':'mndokw'})

# calculate electronic state properties
aiqm1.predict(molecule=mol,
                nstates=3, # Number of electronic states to calculate
                current_state=2, # Requests calculating properties for the S2 electronic state (the third one or second excited)
                calculate_energy=True,
                calculate_energy_gradients=[True, True, True], # requests calculating gradients for all three states
                #calculate_energy_gradients=True # requests calculating gradients only for the current state for many methods (but not for AIQM1 yet)
                calculate_nacv=False, # calculate nonadiabatic coupling vectors
                read_density_matrix=False)

# show the energy and gradients of each state
for i in range(len(mol.electronic_states)):
    print(f'Energy of state {i}, {mol.electronic_states[i].energy}')
    print(f'Energy gradients of state {i}')
    print(mol.electronic_states[i].energy_gradients)

# or alternatively:
print('Electronic state energies:', mol.state_energies)
print('Energy gradients in all electronic states')
print(mol.state_gradients)

# save the molecule object in file
mol.dump(filename="mol.json", format="json")

The output produced by the above Python script will look like this:

Energy of state 0, -94.87746605828193
Energy gradients of state 0
[[-2.20979325e-05 -3.09923481e-05  4.19648750e-04]
[ 1.84361461e-05  7.32504232e-05 -9.31292919e-04]
[-2.19042984e-04  4.91462433e-03 -8.93827657e-04]
[ 2.61506161e-04 -4.71725762e-03 -1.64437563e-03]
[ 1.19805843e-04 -2.90049761e-03  1.30834139e-03]
[-1.58607238e-04  2.66087287e-03  1.74150332e-03]]
Energy of state 1, -94.63515481241645
Energy gradients of state 1
[[ 5.04737245e-04  3.19599783e-03 -4.08726944e-02]
[ 4.00318642e-04  2.59615245e-03 -3.32555938e-02]
[ 2.10040873e-04  2.67707897e-05 -1.31297761e-02]
[ 1.06425061e-04  2.01394042e-03 -1.29785531e-02]
[-1.28339570e-04 -1.37049107e-02  4.93587653e-02]
[-1.09318226e-03  5.87204910e-03  5.08778502e-02]]
Energy of state 2, -94.58438542616139
Energy gradients of state 2
[[-3.22384930e-04 -1.95587619e-03  2.50398870e-02]
[ 1.58487817e-03  1.00715972e-02 -1.28916150e-01]
[-9.11050390e-06  3.73598849e-03 -1.03835407e-02]
[ 2.78260000e-04 -2.07710696e-03 -1.08367181e-02]
[ 2.43944556e-03 -6.96100621e-02  5.75290677e-02]
[-3.97108828e-03  5.98354594e-02  6.75674517e-02]]

Electronic state energies: [-94.87746606 -94.63515481 -94.58438543]
Energy gradients in all electronic states
[[[-2.20979325e-05 -3.09923481e-05  4.19648750e-04]
...

All results of the calculations will be dumped into mol.json which you can inspect to get a better idea how the data are stored in MLatom and how to access different properties through its Python API.

If you used AIQM1/MRCI, you would need to select active space and make other adjustments to the keywords passed to the MNDO program (see its documentation). In our case, you would need to use mndokw file which looks like this:

jop=-2 +
iop=-22 immdp=-1 igeom=1 iform=1 nsav15=3 +
icuts=-1 icutg=-1 kitscf=9999 iscf=9 iplscf=9 +
iprint=-1 kprint=-5 lprint=-2 mprint=0 jprint=-1 +
kharge=1 imult=0 nprint=-1 +
kci=5 movo=-1 ici1=2 ici2=1 jci1=1 jci2=1 nciref=0 +
iroot=3 iuvcd=2 +
ncisym=-1 ioutci=3 ipubo=1 +
ncigrd=3 icross=1

The output would be then:

Energy of state 0, -94.88480497779682
Energy gradients of state 0
[[ 1.34242995e-04  9.67083232e-04 -1.23502153e-02]
[-1.78362286e-04 -1.17842971e-03  1.50824778e-02]
[-9.10435874e-05  2.63194636e-03 -2.00336109e-03]
[ 1.56325835e-04 -2.28966888e-03 -2.38851941e-03]
[ 3.31705685e-05 -9.16708517e-04  7.62840787e-04]
[-5.43335121e-05  7.85777537e-04  8.96776412e-04]]
Energy of state 1, -94.63484529385708
Energy gradients of state 1
[[ 0.00050423  0.00321188 -0.04107714]
[ 0.00044946  0.00290478 -0.0372039 ]
[ 0.00017054  0.00049241 -0.01193939]
[ 0.00012002  0.001369   -0.01187431]
[-0.00012679 -0.01405619  0.0502662 ]
[-0.00111746  0.00607812  0.05182854]]
Energy of state 2, -94.56976745482731
Energy gradients of state 2
[[ 6.21429387e-03  3.96821088e-02 -5.07668904e-01]
[-6.24036735e-03 -3.98880933e-02  5.10344195e-01]
[-9.20936201e-05  2.29106816e-03 -1.36149279e-03]
[ 1.28184773e-04 -2.05330005e-03 -1.70182967e-03]
[ 9.81727653e-05 -2.05300286e-03  3.51054010e-05]
[-1.08190451e-04  2.02121919e-03  3.52925482e-04]]

Calculation via input file

Except for pre-trained ML models, the only way to perform single-point excited-state calculations via input file is to do it for methods supported via the MNDO interface, e.g., AIQM1 and OMx series. An example below is given for AIQM1.

For example, if we want to calculate the vertical excitation energies of ethene in 1B1u state, then you can use the following MLatom input vee.inp:

AIQM1
xyzfile=vee.xyz
mndokeywords=mndokw
yestfile=enest.dat

This input requires a vee.xyz file with XYZ geometries of molecules, e.g. (geometries in Å):

6
Ethene 1B1u
C          0.000000      0.000000     0.668188
C          0.000000      0.000000    -0.668188
H          0.000000      0.923274     1.238289
H          0.000000     -0.923274     1.238289
H          0.000000      0.923274    -1.238289
H          0.000000     -0.923274    -1.238289

and mndokw file with MNDO keywords:

iop=-22 immdp=-1 +
igeom=1 iform=1 +
jop=-2 nsav15=3 +
kci=5 ici1=1 ici2=1 iroot=2 ioutci=2 +
movo=-1 nciref=1 mciref=3 levexc=2

After you prepared your input files vee.inp, vee.xyz, and mndokw you can run MLatom as usual:

mlatom vee.inp

The output of vertical excitation energy is saved in the MNDO program outfile mndo.out, where you can see:

State  2,  Mult. 1,  B1u (4),  E-E(1)=  7.910791 eV,  E=   -304.054934 eV