ROHF/UHF convergence of open-shell doublet VdW diplex

[benj...@gmail.com]

Hi all,

When using Molpro 2015.1 and trying to optimize the geometry of some Van der Waals complexes (N2-NO below) I keep running into ROHF convergence issues. I began with using RKS, B3LYP for the optimization, and also tried UCCSD with aug-cc-pvdz and aug-cc-pvtz basis sets. The calculations never got to the point where the post-HF method mattered b/c the initial single points failed prior to the optimization (these are ignored, but Molpro eventually terminates while calculating the gradients b/c of similar convergence failures).

For other systems I used MCSCF, followed by a single ROHF cycle, to avoid convergence issues, so I tried that again here. Sadly, this did not work (see below input file).

I tried several geometries and have had no problems with any of these in Q-Chem or G16. These use UHF and UB3LYP/UCCSD, but when I tried UHF/UKS in Molpro the calculations still failed.

What thoughts do you have on what is going on and how I can fix these convergence issues?

Thanks,

Benj FitzPatrick

***,geom from a8b
memory,1000,M
gthresh,oneint=1.d-14,twoint=1.d-14,zero=1.d-14

angstrom
symmetry,nosym
geomtyp=xyz
geom={
 N     0.000000     0.000000     0.000000
 N     0.000000     0.000000     1.120000
 N     2.719078     0.000000     2.076060
 O     2.704030     0.167025     3.259679
}

basis={
default,aug-cc-pvdz
}

proc mcscf-dft
 {multi, failsafe, maxit=30,energy=1.d-12,gradient=1.d-7,step=1.d-7; canon,3101.2;
 occ,15;closed,14;}
 {rhf;start,3101.2;save,2101.2;maxit,1;}
 {rks,b3lyp;}
endproc

{optg,proc=mcscf-dft,root=1,displace=symm,numhess=1,maxit=2;
coord,3n,normal;
step,0.1}

proc mcscf-ccsd
 {multi, failsafe, maxit=30,energy=1.d-12,gradient=1.d-7,step=1.d-7; canon,3100.2;
 occ,15;closed,14;}
 {rhf;start,3100.2;save,2100.2;maxit,1;}
 {uccsd,maxit=100}
endproc

{optg,proc=mcscf-ccsd,root=1,maxit=30,displace=symm,numhess=1,hessproc=dft_hess;
coord,3n,normal;
step,0.1}

proc dft_hess
  {rhf}
  {rks,b3lyp}
endproc

PROGRAM * MULTI (Direct Multiconfiguration SCF)       Authors: P.J. Knowles, H.-J. Werner (1984)     S.T. Elbert (1988)


 Number of closed-shell orbitals:  14 (  14 )
 Number of active  orbitals:        1 (   1 )
 Number of external orbitals:      77 (  77 )

 State symmetry 1

 Number of electrons:     1    Spin symmetry=Doublet   Space symmetry=1
 Number of states:        1
 Number of CSFs:          1   (1 determinants, 1 intermediate states)

 Orbital guess generated from atomic densities. Full valence occupancy:   20

 Wavefunction dump at record             3101.2

 Convergence thresholds  0.10E-06 (gradient)  0.10E-11 (energy)  0.10E-06 (step length)

 Number of orbital rotations:     1169   (    14 Core/Active   1078 Core/Virtual   0 Active/Active     77 Active/Virtual)
 Total number of variables:       1170


 ITER. MIC  NCI  NEG     ENERGY(VAR)     ENERGY(PROJ)   ENERGY CHANGE     GRAD(0)  GRAD(ORB)   GRAD(CI)     STEP       TIME

   1   40    7    0    -238.14655608    -238.68106597   -0.53450989    0.47437644 0.00034582 0.00000000  0.28D+00      0.87
   2   27   17    0    -238.02149441    -238.19636021   -0.17486580    0.85278015 0.00000000 0.00000000  0.23D+00      1.45
   3   36   17    0    -238.20536408    -238.20714499   -0.00178090    0.07224975 0.00000003 0.00000000  0.35D-01      2.06
   4   35   17    0    -238.20717140    -238.65449157   -0.44732017    0.00157092 0.00000000 0.00000000  0.25D+00      2.80
   5   30   17    0    -237.97655909    -238.19056610   -0.21400701    0.95773883 0.00000003 0.00000000  0.24D+00      3.41
   6   30   17    0    -238.20329674    -238.20706544   -0.00376870    0.10629393 0.00000003 0.00000000  0.48D-01      3.97
   7   37   17    0    -238.20716552    -238.65540885   -0.44824333    0.00449461 0.00000000 0.00000000  0.25D+00      4.70
   8   28   17    0    -237.97511115    -238.19025748   -0.21514632    0.96034377 0.00000001 0.00000000  0.24D+00      5.28
   9   32   17    0    -238.20316868    -238.20706246   -0.00389378    0.10780526 0.00000006 0.00000000  0.49D-01      5.89
  10   35   17    0    -238.20716849    -238.65538946   -0.44822097    0.00467421 0.00000000 0.00000000  0.25D+00      6.64
  11   30   17    0    -237.97505302    -238.19024640   -0.21519338    0.96045341 0.00000003 0.00000000  0.24D+00      7.26
  12   29   17    0    -238.20316573    -238.20706594   -0.00390021    0.10787771 0.00000003 0.00000000  0.49D-01      7.84
  13   37   17    0    -238.20717224    -238.65532930   -0.44815706    0.00468313 0.00000000 0.00000000  0.25D+00      8.58
  14   27   17    0    -237.97507765    -238.19025122   -0.21517357    0.96041499 0.00000001 0.00000000  0.24D+00      9.15
  15   32   17    0    -238.20316990    -238.20706970   -0.00389980    0.10787026 0.00000005 0.00000000  0.49D-01      9.77
  16   33   17    0    -238.20717592    -238.65526985   -0.44809393    0.00468207 0.00000000 0.00000000  0.25D+00     10.53
  17   30   17    0    -237.97510619    -238.19025837   -0.21515218    0.96036942 0.00000001 0.00000000  0.24D+00     11.13
  18   30   17    0    -238.20317447    -238.20707334   -0.00389887    0.10785563 0.00000003 0.00000000  0.49D-01     11.73
  19   37   17    0    -238.20717947    -238.65521142   -0.44803195    0.00468038 0.00000000 0.00000000  0.25D+00     12.49
  20   27   17    0    -237.97513406    -238.19026384   -0.21512978    0.96032484 0.00000001 0.00000000  0.24D+00     13.11
  21   32   17    0    -238.20317878    -238.20707686   -0.00389808    0.10784410 0.00000005 0.00000000  0.48D-01     13.73
  22   34   17    0    -238.20718290    -238.65515564   -0.44797275    0.00467877 0.00000000 0.00000000  0.25D+00     14.48
  23   29   17    0    -237.97516093    -238.19027062   -0.21510969    0.96028190 0.00000002 0.00000000  0.24D+00     15.10
  24   30   17    0    -238.20318308    -238.20708025   -0.00389717    0.10782992 0.00000003 0.00000000  0.48D-01     15.69
  25   37   17    0    -238.20718620    -238.65510168   -0.44791548    0.00467712 0.00000000 0.00000000  0.25D+00     16.43
  26   27   17    0    -237.97518688    -238.19027571   -0.21508883    0.96024037 0.00000001 0.00000000  0.24D+00     17.05
  27   32   17    0    -238.20318708    -238.20708353   -0.00389644    0.10781930 0.00000005 0.00000000  0.48D-01     17.67
  28   34   17    0    -238.20718939    -238.65505035   -0.44786096    0.00467559 0.00000000 0.00000000  0.25D+00     18.44
  29   29   17    0    -237.97521185    -238.19028108   -0.21506923    0.96020046 0.00000002 0.00000000  0.24D+00     19.05
  30   30   17    0    -238.20319092    -238.20708667   -0.00389575    0.10780730 0.00000003 0.00000000  0.48D-01     19.64

 ** WVFN ****  MAXIMUM NUMBER OF ITERATIONS REACHED

[Hans-Joachim Werner]

Dear Benjamin,
apparently there are two problems: (i) your ROHF does not converge, and (ii) it is very difficult to optimize the geometry of the vdW complex.

(i) Indeed ROHF convergence is difficult to achieve with the standard procedure, and even the second-order mcscf has problems. Molpro 2015 converges with

rhf,maxdiis=20;shift,-0.4,-0.2

to the correct solution, but the convergence depends rather sensitively on the shift parameters and I don’t know if it would work for all other geometries.

In the current version (Molpro 2021.2) we have completely new rohf and mcscf/casscf implementations, and with these robust and much faster convergence is obtained (using rhf,so-sci, which is not available in Molpro2015). I recommend to upgrade your Molpro version, if you can afford it.

(ii) The geometry optimization is another difficulty, since the intermolecular potential may be very flat in some coordinates, while the n2 and no bonds are stiff. Unfortunately, we do not have open-shell ccsd gradients, and so the gradients are computed numerically by finite differences. I would recommend to keep in a first step the N2 and NO distances fixed, and to optimize only the intermolecular coordinates. For these one needs larger displacements than for the N2 and NO bonds in order to get sufficiently accurate gradients. Also, counter-poise corrections and triple excitations should be included. In view of the many displacements needed, it may be easier to compute the PES around the minimum and fit this rather than using the geometry optimizer.

I am sorry not to have a simpler solution.

Best regards
Joachim Werner

> --
> You received this message because you are subscribed to the Google Groups "molpro-user" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to molpro-user...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/8275d32f-fd63-4cb7-89b6-9280da77b4b0n%40googlegroups.com.

[Benj FitzPatrick]

Dr. Werner,

Thank you very much for your thorough response.

Agreed, I would love to upgrade to the newest version, though that is currently out of reach (I'm an independent researcher doing this with my own time and money). This is especially true if UCCSD and UCCSD(T) get analytic gradients in future versions (or the F12 variants).

The level shift is good to know, I have not had much luck with those in the past as you noted their sensitivity to geometry.

I will give your single point strategy + fitting a try (the eventual goal is to look at the A state as well).

Thanks,

Benj