Thank you for your feedback

Best,

Tian

Best Regards,
Kang

I do not easily suggest using scf=dm/xqc, because they are not tested with sobEDA.sh and it is well-known they often converge to unstable wavefunction. You may try to use other ways to solve the unconvergence, such as scf=vshift=xxx, scf=novaracc, other basis sets or theoretical methods, and so on.

Frankly speaking, cleavage of F2 is fairly challenging for sobEDA, the main reason is that when nosymm is used, the orientation of all p orbitals of doublet F (involved in homolytic cleavage) and singlet F+ (involved in heterolytic cleavage) is fully random, however their wavefunctions are needed to construct initial wavefunction of promolecular state, which introduces randomness in physical components given by sobEDA.

The only workaround in the framework of sobEDA is adding a very slightly external perturbation in the calculation of doublet F and singlet F+, the most easiest way I think is adding a very small external field, namely adding "field=z+1" keyword in the F and F+ calculations (I've checked orbital isosurface map using Multiwfn, this treatment indeed removes arbitrariness of orientation of pz orbital for doublet F), this treatment should have negligible impact on total energy, but may make the energy component curves much smoother. I am not sure if this trick really works, but it is worth to try.

Best,

Tian

I performed sobEDA calculations for F2's heterolytic cleavage, considering various bond lengths. Starting from 0.5 AU, I incremented by 0.15 AU until reaching 2.9 AU. Below is a summary of the results. As indicated in the summary.txt, there were errors encountered for several bond lengths during the calculations. Specifically, I have attached the result files for 0.95 and 1.1 AU.

# fragment.txt
2
1 1
1
-1 1
2

# two system.xyz files
--------------------------------------------------------
2
bond length = 0.95 AU
 F    0.000000000000    0.000000000000   -0.475000000000
 F    0.000000000000    0.000000000000    0.475000000000
--------------------------------------------------------
2
bond length = 1.1 AU
 F    0.000000000000    0.000000000000   -0.550000000000
 F    0.000000000000    0.000000000000    0.550000000000
 --------------------------------------------------------

# template.gjf
#p b3lyp/6-31G* em=gd3bj ExtraLinks=L608 nosymm

opted geometry

0 1
[geometry]

-5 5

# summary.txt
# F2 heterolytic cleavage result

--------------------------------------------------------------------
length      E_Tot       E_Els       E_Pauli     E_Orb        E_Ccorr
--------------------------------------------------------------------

  0.50        4955.52      694.12     7276.44    -2972.87      -42.17
  0.65        1547.31     -636.36     4203.55    -1984.73      -35.15
  0.80         262.67     -755.80     2380.40    -1333.16      -28.77
  0.95    -------------------------- Error --------------------------
  1.10        -421.39     -645.94     1236.17     -991.57      -20.05
  1.25    -------------------------- Error --------------------------
  1.40    -------------------------- Error --------------------------
  1.55    -------------------------- Error --------------------------
  1.70    -------------------------- Error --------------------------
  1.85        -462.37     -183.87       48.82     -321.65       -5.68
  2.00    -------------------------- Error --------------------------
  2.15    -------------------------- Error --------------------------
  2.30    -------------------------- Error --------------------------
  2.45        -410.30     -132.81        3.44     -279.21       -1.71
  2.60        -402.32     -132.81        0.81     -269.47       -0.85
  2.75	  -------------------------- Error --------------------------
  2.90        -380.62     -112.67        0.58     -267.95       -0.58

Upon referencing the attached output files, it was observed that the SCF convergence was not achieved using the default DIIS method for the result at 0.95 AU.

# google_drive/0.95AU/fragment1.out
.
.
>>>>>>>>>> Convergence criterion not met.
 SCF Done:  E(RB3LYP) =  -98.9435464416     A.U. after  129 cycles
            NFock=128  Conv=0.25D-07     -V/T= 2.0083
 KE= 9.812919884498D+01 PE=-2.310227613054D+02 EE= 3.395001601885D+01
 Convergence failure -- run terminated.
 .
 .

So, I added the keywords "scf=dm" and "scf=xqc" to the template.gjf file and performed two additional calculations. The EDA results seemed to be obtained normally, but some concerning error codes were also outputted. Can we trust the results despite the error codes, disregarding them?

# google_drive/0.95AU_scf=dm/result.txt
...
E_els = -354.7132107182 Hartree
E_x = -19.387852 Hartree
E_c = -0.859346 Hartree
E_disp = -0.0002765218 Hartree
Frozen state energy: 1) Hartree
(standard_in) 1: syntax error		<-- syntax error?
(standard_in) 1: syntax error		<-- syntax error?

*************************
***** Final results *****
*************************

Total interaction energy:    878.17 kcal/mol

Physical components of interaction energy derived by sobEDA:
...

Below is a plot of the results obtained by adding the "scf=dm" keyword. Using "scf=dm" or "scf=xqc" keywords led to convergence on more IRC points than DIIS, but still fails on some points due to similar errors. Error IRC points were excluded from the plot. The interaction energy plot appears to be not as smooth as expected. Could you advise on the possible reasons for this?
(I also performed homolytic cleavage for F2, but the resulting curve is also not smooth.)

https://drive.google.com/file/d/1NXJtO1 … drive_link

On the other hand, smooth energy curves were successfully obtained for the calculation results of the heterolytic cleavage of H2O2.

https://drive.google.com/file/d/1LSTc66 … drive_link

Detail output files:
https://drive.google.com/drive/folders/ … drive_link

Best Regards,
Kang