Azobenzene PES scan along dihedral does not converge and does not behave as intended

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
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()
0 1
H 1 0.74

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
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)
    except RuntimeError:
        psi4.core.print_out("WARNING: Did not converge")
        ecp[counter] = -1337.0
    except ConvergenceError:
        psi4.core.print_out("WARNING: Did not converge")
        ecp[counter] = -1337.0
    core.set_global_option("BASIS", "cc-pVDZ")
    core.set_global_option("SCF_TYPE", "df")
    core.set_global_option("GUESS", "sad")
        ecp[counter] = energy('scf')
    except RuntimeError:
        psi4.core.print_out("WARNING: SPE Did not converge")
        ecp[counter] = -1338.0
    except ConvergenceError:
        psi4.core.print_out("WARNING: Did not converge")
        ecp[counter] = -1337.0
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")
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:

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?

If this is still a concern for you:

  1. As far as I’m aware, Psi has nothing to deal with conical intersections
  2. Ignoring failures of SCF and geometry convergence is not going to be helpful in investigating whatever-you-have. I can take a look at this if you give me a plain input file, not a modified auto-generated one. I tried working backwards to what the original input file was, but I was getting bizarre results, so I’d like to double-check what I did.