Atomic Forces and Geometry Optimizaiton

What this tutorial is about

In this tutorial we will learn how to compute the excited state atomic forces of a molecular structure and how to use these to perform geometry optimization. Let’s start by importing the relevant modules

[1]:

from pyxtp import xtp
from ase.io import write
from ase.build import molecule
from ase.optimize import GoodOldQuasiNewton

Define the molecular structure

We define here a CO molecule in its ground state geometry.

[2]:

atoms = molecule('CO')

Configure the xtp calculator

We now define the xtp calculator and configure it to target particular forces. We can select which forces we want to compute with the method .select_force() of the calculator. We can choose to target different forces :

  • energy='energy' : DFT total energy

  • energy='singlets': BSE singlet energy

  • energy='triplets': BSE triplet energy

  • energy='qp': Quasi particile energy

For singlets, triplets and qp one can also specify which particular level to target to compute the forces. For example energy='singlets', level=0 will target the lowest lying singlet state.

[3]:

# instantiate the calculator
calc = xtp(nthreads=2)

# select the force we want to use
calc.select_force(energy='singlets', level=0, dynamic=False)

# this allows to change all options
# calc.options.dftpackage.functional = 'PBE'
calc.options.dftpackage.basisset = 'def2-svp'
calc.options.dftpackage.auxbasisset = 'aux-def2-svp'
calc.options.gwbse.gw.qp_grid_search_mode = 'dense'
calc.options.gwbse.gw.qp_grid_spacing = 0.001
calc.options.gwbse.gw.qp_grid_steps = 1001

# set up the logger
calc.options.logging_file = 'CO_forces.log'

# set the calculator
atoms.calc = calc

Compute the forces

If you are simply interested in computing the forces, they can easily be accessed through the .get_forces() method

[4]:

atoms.get_forces()

[4]:

array([[-7.53175300e-10, -7.53175300e-10,  3.19229776e-01],
       [ 6.82121026e-10,  6.96331881e-10, -3.19229776e-01]])

Geometry optimization

Geometry optimization can be run by leveraging the intrinsic ASE capabilities. We can for example use the QuasiNewton method implemented in ASE to relax the molecular structre in the excited states we have just specified. We here fix the number of steps to 10 to limit the computational cost.

[5]:

dyn = GoodOldQuasiNewton(atoms, trajectory='test.traj')
dyn.run(fmax=0.05, steps=10)
write('final.xyz', atoms)

                    Step     Time          Energy          fmax
GoodOldQuasiNewton:    0 08:39:42     -113.101204        0.319230
new radius 0.012247
eigenvalues 20.00 20.00 20.00
Corrected Newton step: abs(D) = 0.02
Abs Gbar estimate 0.042646422477687254
GoodOldQuasiNewton:    1 08:40:28     -113.098913        0.319230
energies -113.0989134292243 -113.10120446793485
reject step
new radius 0.006124
eigenvalues 20.00 20.00 20.00
Corrected Newton step: abs(D) = 0.02
Abs Gbar estimate 0.10823086221853342
GoodOldQuasiNewton:    2 08:41:16     -113.100204        0.248237
energies -113.10020359454646 -113.10120446793485
reject step
new radius 0.003062
eigenvalues 20.00 20.00 20.00
Corrected Newton step: abs(D) = 0.02
Abs Gbar estimate 0.15227308033174852
GoodOldQuasiNewton:    3 08:42:04     -113.100741        0.282892
energies -113.1007413294438 -113.10120446793485
reject step
new radius 0.001531
eigenvalues 20.00 20.00 20.00
Corrected Newton step: abs(D) = 0.02
Abs Gbar estimate 0.17710669010008095
GoodOldQuasiNewton:    4 08:42:51     -113.100982        0.300816
energies -113.10098236796148 -113.10120446793485
reject step
new radius 0.000765
eigenvalues 20.00 20.00 20.00
Corrected Newton step: abs(D) = 0.02
Abs Gbar estimate 0.19022661986696826
GoodOldQuasiNewton:    5 08:43:38     -113.101096        0.310054
energies -113.1010958025806 -113.10120446793485
reject step
new radius 0.000383
eigenvalues 20.00 20.00 20.00
Corrected Newton step: abs(D) = 0.02
Abs Gbar estimate 0.19696236601507155
GoodOldQuasiNewton:    6 08:44:24     -113.101151        0.314829
energies -113.10115073498848 -113.10120446793485
energy change; actual: 0.000054 estimated: -0.000171
Energy prediction factor 1.3136347101201804
Force prediction factor -1.6250813122652765
Scale factors 0.364431 0.364431
new radius 0.000139
eigenvalues 7.62 20.00 20.00
Corrected Newton step: abs(D) = 0.06
Abs Gbar estimate 0.20023874313710724
GoodOldQuasiNewton:    7 08:45:11     -113.101131        0.317168
energies -113.10113085428141 -113.10115073498848
energy change; actual: 0.000020 estimated: -0.000062
Energy prediction factor 1.318149213725695
Force prediction factor 0.1283867435690288
Scale factors 0.364431 0.510204
new radius 0.000100
eigenvalues 8.74 20.00 20.00
Corrected Newton step: abs(D) = 0.05
Abs Gbar estimate 0.19931742848174566
GoodOldQuasiNewton:    8 08:45:58     -113.101117        0.316306
energies -113.10111650360759 -113.10113085428141
energy change; actual: 0.000014 estimated: -0.000045
Energy prediction factor 1.3211254389053801
Force prediction factor 0.2542237295882319
Scale factors 0.364431 0.510204
new radius 0.000100
eigenvalues 11.72 20.00 20.00
Corrected Newton step: abs(D) = 0.04
Abs Gbar estimate 0.19800698922037333
GoodOldQuasiNewton:    9 08:46:44     -113.101102        0.315477
energies -113.10110207138129 -113.10111650360759
energy change; actual: 0.000014 estimated: -0.000045
Energy prediction factor 1.3239078563198683
Force prediction factor -0.9962336589739281
Scale factors 0.364431 0.364431
new radius 0.000100
eigenvalues 5.87 20.00 20.00
Corrected Newton step: abs(D) = 0.08
Abs Gbar estimate 0.1980050211877179
GoodOldQuasiNewton:   10 08:47:30     -113.101088        0.315062