i was trying to examine the PES of azobenzene while performing a 180° rotation along the dihedral spanning the azo-group in the middle.
i ran into similar problems as the person in another thread where some steps would not converge. i wanted to investigate what was happening, so i modified the generated python code
to ignore any convergence errors and continue doing whatever it was doing and outputting the results. in order to see what it attempts
to do. i subsequently visualized all geometries that were written to the console.
i was using the following code to do that
import sys
sys.path.insert(1,'/usr/lib')
import psi4
from psi4 import *
from psi4.core import *
from psi4.driver.diatomic import anharmonicity
from psi4.driver.gaussian_n import *
from psi4.driver.aliases import *
from psi4.driver.driver_cbs import *
from psi4.driver.wrapper_database import database, db, DB_RGT, DB_RXN
from psi4.driver.wrapper_autofrag import auto_fragments
from psi4.driver.constants.physconst import *
psi4_io = core.IOManager.shared_object()
psi4_io.set_default_path("/tmp/")
geometry("""
0 1
H
H 1 0.74
""","blank_molecule_psi4_yo")
core.set_memory_bytes(4000000000)
core.efp_init()
Azobenzene = geometry("""
units angstrom
C -2.39739 -0.78827 0.00735
C -2.34437 0.60598 -0.01640
C -1.10938 1.25611 -0.02197
C 0.07499 0.51774 -0.00380
C 0.02012 -0.87819 0.01973
C 3.58030 1.10357 0.00379
C 4.76468 0.36518 0.02130
C 6.05268 2.40957 -0.00736
C 5.99966 1.01532 0.01569
C 3.63517 2.49951 -0.01902
C 4.87006 3.15035 -0.02457
H -3.26195 1.18396 -0.03042
H -1.06845 2.33955 -0.04037
C -1.21477 -1.52904 0.02529
H -3.35614 -1.29539 0.01183
H 2.72434 3.08743 -0.03255
H 4.91048 4.23401 -0.04238
H 7.01143 2.91670 -0.01186
H 6.91725 0.43733 0.02915
H 4.72375 -0.71826 0.03913
H -1.25519 -2.61268 0.04365
H 0.93094 -1.46610 0.03384
N 1.25619 1.16647 -0.00911
N 2.39910 0.45483 0.00905
""","Azobenzene")
core.IO.set_default_namespace("Azobenzene")
steps = 18
step_size = 10
ecp = {}
for counter in range (0, steps):
core.set_global_option("BASIS", "sto-3g")
core.set_global_option("GEOM_MAXITER", 100)
core.set_global_option("SCF_TYPE", "df")
core.set_global_option("GUESS", "sad")
core.set_global_option("INTS_TOLERANCE", 1.0E-8)
dihedral = 180.1 - (1.0 * counter * step_size)
dihedral_string = "4 23 24 6 %.1f" % dihedral
core.set_local_option("OPTKING", "FROZEN_DIHEDRAL", "")
core.set_local_option("OPTKING", "FIXED_DIHEDRAL", dihedral_string)
try:
optimize('scf')
except RuntimeError:
psi4.core.print_out("WARNING: Did not converge")
ecp[counter] = -1337.0
continue
except ConvergenceError:
psi4.core.print_out("WARNING: Did not converge")
ecp[counter] = -1337.0
continue
core.set_global_option("BASIS", "cc-pVDZ")
core.set_global_option("SCF_TYPE", "df")
core.set_global_option("GUESS", "sad")
try:
ecp[counter] = energy('scf')
except RuntimeError:
psi4.core.print_out("WARNING: SPE Did not converge")
ecp[counter] = -1338.0
continue
except ConvergenceError:
psi4.core.print_out("WARNING: Did not converge")
ecp[counter] = -1337.0
continue
psi4.core.print_out("\n")
psi4.core.print_out("Energies for different dihedral angles with sto-3g \n\n")
psi4.core.print_out(" Dihedral [degree] E_int [kcal/mol] \n")
psi4.core.print_out("-----------------------------------------------------\n")
for counter in range (0, steps):
dihedral = 180.1 - (1.0 * counter * step_size)
e = ecp[counter] * psi4.constants.hartree2kcalmol
psi4.core.print_out(" %3.1f %10.6f\n" % (dihedral, e))
afterwards extracted all geometry frames generated in the process and rendered a video from it that can be seen here:
https://rainbow-lollipop.de/azobenzene.avi
what i expected it to do at every loop step according to docs and tutorials:
* create a geometry where the dihedral is adjusted at a specific angle
* optimize that geometry to the nearest local minimum without changing the dihedral but everything else
* calculate the energy of that minimum
what i expected the results to look like:
* i expected that the benzene rings would rotate along the C3-N23 and C4-N24 axes to a local minimum at each step
* i expected to get a pointy energy curve as i know that i cut the hypersurface in the proximity of a conical intersection.
what actually happened:
* optimize seems to perform a step-by-step approximation of the whole rotation multiple times
* therefore the energy constantly rises, causing the optimization to not converge
* it seems as if the whole geometry was fixed
apparently i misinterpret how the optimize command works with constraints. can someone shine light on this issue?
how would i do something like “rotate dihedral for 10° at each step, optimize, give me single point energy”
bunch of other questions:
* is it possible to perform a search for conical intersections between two PES of the same molecule at different excitations? If so, how would one approach that in psi4?
* which algorithm would you suggest to investigate conical intersections? casscf seems appropriate but manual
active space selection is cumbersome. are there any more novel possibilites that do not need as much interaction?