Proper use of 'cart_hess_read' option in OPTKING

kalju · May 21, 2022, 2:35am

Dear Psi4 Experts,

I have been teaching a Computational Chemistry class using Psi4 (and hope to move couple of examples to Psi4Education soon) this quarter and for the most part it has worked out well. One topic that we are focusing more deeply is the efficiency of different geometry optimization methods, and in this context I was trying to probe the benefits of the ‘cart_hess_read true’ setting for efficiency (e.g. using a DF-MP2 Hessian for CCSD optimizations).

Let’s consider the case of a simple diatomic, CO (req near 1.1). I am intentionally (for teaching purposes) starting the geometry optimization at the region (r = 1.35) where second derivative is small but still positive (the inflection point is near 1.45) and SCF converges well. When I set ‘cart_hess_read false’ and ‘full_hess_every 0’, the minimization succeeds. With Newton-Raphson, the first step overshoots as expected, with RFO not so much and then BFGS corrects. All good!

But if I read in the previously computed Hessian from the disk with ‘cart_hess_read true’ (at the same HF level and at the same geometry), minimization fails as it oscillates back and forth around the minimum. The first step that will be taken is the Newton-Raphson step, this makes sense. But in the next step, the same wrong Hessian will be used again from the disk (unless I force Hessian recalculation with ‘full_hess_every 1’).

I am not sure if the Hessian never gets updated by BFGS here, or if it gets updated and then ignored because I have specified ‘cart_hess_read true’. The ‘full_hess_every -1’ or ‘full_hess_every 0’ does not seem to make difference. I feel that I am missing something and I admit I am not quite following what is happening in driver.py regarding this logic. Below is a pertinent input:

# understanding Hessian utilization in optimization
# Psi4 1.6a1 with oneMKL and libinit2.7 (5-4-3-6-5-4)

memory 2 GB

molecule {
 0 1
C      0.000000    0.000000    0.000000
O      1.350000    0.000000    0.000000
}

set globals {
basis aug-cc-pvtz
scf_type df
}

set scf {
df_scf_guess true
e_convergence 1e-10
}


set hessian_write true

psi4.print_out("\n ~~~  HESSIAN at the HF/aug-cc-pVTZ level ~~~\n")
e_avtz, wfn_avtz = psi4.frequencies('scf', return_wfn=True)

set scf {
guess = read
}

set optking {
step_type nr
cart_hess_read true
print_opt_params true
print_trajectory_xyz_file true
full_hess_every -1
}

psi4.print_out("\n ~~~  OPTIMIZATION at the HF/aug-cc-pVTZ level ~~~\n")
optimize('scf')

So, my question is: is there a good way to read in a not-so-great Hessian computed beforehand at not-so-great geometry with the hope that it will accelerate the convergence if the search is started from that same not-so-great geometry (assuming that the Hessian is positive and large enough not to squish the atoms totally together in the first step.

P.S. This is with Psi4 1.6a1.dev155 as well as with older versions.

Thank you,

Kalju

jmisiewicz · May 21, 2022, 11:06am

@AlexanderH

Mind taking a look at this?

Rollin_King · May 24, 2022, 5:32pm

@kalju , you have described the problem and behavior very clearly.

Please try something like

cart_hess_read true
full_hess_every 20

(or some other large number). This is a bit obscure behavior, but I think the two keywords are not playing nicely together when you (inadvertently) ask to read and to compute the cartesian hessian (the latter of which is already done by the frequencies() line.)

kalju · May 24, 2022, 10:25pm

@Rollin_King

Sure, this gets the geometry to converge in about 25 steps in this specific case because we happen to be in a region close to the minimum at step 20, and the Hessian is now much better.

co_hessian_test_short_20.out:       1 R(1,2)          =      1.085411     0.472854     0.168431     1.253843
co_hessian_test_short_20.out:       1 R(1,2)          =      1.253843    -2.176527    -0.264589     0.989254
co_hessian_test_short_20.out:       1 R(1,2)          =      0.989254     4.060514     0.264589     1.253843
co_hessian_test_short_20.out:       1 R(1,2)          =      1.253843    -2.176527    -0.066147     1.187695
co_hessian_test_short_20.out:       1 R(1,2)          =      1.187695    -1.499733    -0.198441     0.989254
co_hessian_test_short_20.out:       1 R(1,2)          =      0.989254     4.060514     0.198441     1.187695
co_hessian_test_short_20.out:       1 R(1,2)          =      1.187695    -1.499733    -0.049610     1.138085
co_hessian_test_short_20.out:       1 R(1,2)          =      1.138085    -0.716661    -0.148831     0.989254
co_hessian_test_short_20.out:       1 R(1,2)          =      0.989254     4.060514     0.148831     1.138085
co_hessian_test_short_20.out:       1 R(1,2)          =      1.138085    -0.716661    -0.148831     0.989254
co_hessian_test_short_20.out:       1 R(1,2)          =      0.989254     4.060514     0.148831     1.138085
co_hessian_test_short_20.out:       1 R(1,2)          =      1.138085    -0.716661    -0.148831     0.989254
co_hessian_test_short_20.out:       1 R(1,2)          =      0.989254     4.060514     0.148831     1.138085
co_hessian_test_short_20.out:       1 R(1,2)          =      1.138085    -0.716661    -0.148831     0.989254
co_hessian_test_short_20.out:       1 R(1,2)          =      0.989254     4.060514     0.148831     1.138085
co_hessian_test_short_20.out:       1 R(1,2)          =      1.138085    -0.716661    -0.148831     0.989254
co_hessian_test_short_20.out:       1 R(1,2)          =      0.989254     4.060514     0.148831     1.138085
co_hessian_test_short_20.out:       1 R(1,2)          =      1.138085    -0.716661    -0.148831     0.989254
co_hessian_test_short_20.out:       1 R(1,2)          =      0.989254     4.060514     0.148831     1.138085
co_hessian_test_short_20.out:       1 R(1,2)          =      1.138085    -0.716661    -0.038233     1.099852
co_hessian_test_short_20.out:       1 R(1,2)          =      1.099852     0.103425     0.005518     1.105369
co_hessian_test_short_20.out:       1 R(1,2)          =      1.105369    -0.028538    -0.001522     1.103847
co_hessian_test_short_20.out:       1 R(1,2)          =      1.103847     0.007385     0.000394     1.104241
co_hessian_test_short_20.out:       1 R(1,2)          =      1.104241    -0.001946    -0.000104     1.104137

But there are two issues with this approach:

my ultimate goal is to read in the Hessian that I have calculated at the lower level only once and then have BFGS update it during optimization.
with default (RFO) step, trying ‘full_hess_every 20’ just leads to different oscillatory failure as step 20 of RFO happens to be not so great. And the new step 40 is not much better … So, this is not a general solution:

co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.105411    -0.029526    -0.014111     1.091300
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.091300     0.317874     0.141202     1.232503
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.232503    -1.996116    -0.402462     0.830040
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.830040    17.681510     0.512879     1.342919
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.342919    -2.639407    -0.429759     0.913160
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.913160     8.933931     0.497422     1.410582
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.410582    -2.768209    -0.132295     1.278288
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.278288    -2.346032    -0.396883     0.881405
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.881405    11.757617     0.396883     1.278288
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.278288    -2.346032    -0.396883     0.881405
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.881405    11.757617     0.396883     1.278288
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.278288    -2.346032    -0.396883     0.881405
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.881405    11.757616     0.396883     1.278288
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.278288    -2.346032    -0.396883     0.881405
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.881405    11.757616     0.396883     1.278288
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.278288    -2.346032    -0.396883     0.881405
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.881405    11.757616     0.396883     1.278288
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.278288    -2.346032    -0.396883     0.881405
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.881405    11.757616     0.396883     1.278288
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.278288    -2.346032    -0.275881     1.002407
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.002407     3.429886     0.333667     1.336074
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.336074    -2.617369    -0.293047     1.043027
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.043027     1.785026     0.232575     1.275602
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.275602    -2.329190    -0.274742     1.000860
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.000860     3.501376     0.336622     1.337481
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.337481    -2.622056    -0.293325     1.044156
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.044156     1.745141     0.229027     1.273183
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.273183    -2.313661    -0.273684     0.999499
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.999499     3.564836     0.339175     1.338674
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.338674    -2.625964    -0.293557     1.045117
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.045117     1.711433     0.225974     1.271091
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.271091    -2.299951    -0.272744     0.998347
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.998347     3.618952     0.341303     1.339650
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.339650    -2.629119    -0.293743     1.045907
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.045907     1.683899     0.223443     1.269350
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.269350    -2.288335    -0.271943     0.997407
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      0.997407     3.663427     0.343018     1.340425
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.340425    -2.631597    -0.293890     1.046536
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.046536     1.662076     0.221413     1.267948
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.267948    -2.278853    -0.256226     1.011722
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.011722     3.014067     0.300206     1.311928
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.311928    -2.523682    -0.272421     1.039507
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.039507     1.911202     0.228404     1.267911
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.267911    -2.278597    -0.256208     1.011703
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.011703     3.014922     0.300249     1.311952
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.311952    -2.523787    -0.272428     1.039524
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.039524     1.910575     0.228353     1.267877
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.267877    -2.278365    -0.256192     1.011685
co_hessian_test_short_rfo_20.out:           1 R(1,2)          =      1.011685     3.015697     0.300289     1.311974

The practical solution, inspired by the IRC test example, is to move ‘cart_hess_read true’ from the optking scope into globals:

# understanding Hessian utilization in optimization
# Psi4 1.6a1 with oneMKL and libinit2.7 (5-4-3-6-5-4)

memory 2 GB

molecule {
 0 1
C      0.000000    0.000000    0.000000
O      1.350000    0.000000    0.000000
}

set globals {
 basis aug-cc-pvtz
 scf_type df
}

set scf {
 df_scf_guess true
 e_convergence 1e-10
}

set hessian_write true

psi4.print_out("\n ~~~  HESSIAN at the HF/aug-cc-pVTZ level ~~~\n")
e_avtz, wfn_avtz = psi4.frequencies('scf', return_wfn=True)

set scf {
 guess = read
}

set globals {
 cart_hess_read true
}

set optking {
 step_type nr
 print_opt_params true
 print_trajectory_xyz_file true
 full_hess_every -1
}

psi4.print_out("\n ~~~  OPTIMIZATION at the HF/aug-cc-pVTZ level ~~~\n")
optimize('scf')

Now, the calculation converges as expected

            1 R(1,2)          =      1.350000    -2.660250    -0.264589     1.085411
            1 R(1,2)          =      1.085411     0.472854     0.039932     1.125344
            1 R(1,2)          =      1.125344    -0.466916    -0.047825     1.077519
            1 R(1,2)          =      1.077519     0.690318     0.025054     1.102573
            1 R(1,2)          =      1.102573     0.037728     0.001702     1.104275
            1 R(1,2)          =      1.104275    -0.002762    -0.000107     1.104169
            1 R(1,2)          =      1.104169    -0.000240    -0.000009     1.104159

as the Hessian is read in once and gets updated by BFGS in subsequent optimization steps:

read_cartesian_H       =               true
read_cartesian_H       =              false
read_cartesian_H       =              false
read_cartesian_H       =              false
read_cartesian_H       =              false
read_cartesian_H       =              false
read_cartesian_H       =              false

I am not sure if the scope-dependence of ‘cart_hess_read’ option is the desired behavior of Psi4.

Thanks for looking into this!

Rollin_King · May 25, 2022, 12:05am

I am not sure if that is the desired behavior or not either Congratulations on figuring out the scope issue. Apparently, the psi4 driver uses the read_cartesian_H keyword and this was getting hidden from it by the optking scope. I haven’t reviewed that code in awhile. Both the driver and the optimizer are in the middle of a major overhaul, so this kind of thing should be simpler in the future - but glad you figured out how to make it work!

kalju · May 25, 2022, 12:54am

I see. As always, looking forward to what the new features (optimizer in this case) will have to offer!