Bonded ALMO-EDA for F2

Questions
#1

I’m trying to get the bonded ALMO-EDA for F2, which was successfully run in the PNAS paper [1]. However, the calculations have been failing in QChem 6.0.1. Below is one input and its corresponding error message:

--------------------------------------------------------------
User input:
--------------------------------------------------------------
$molecule
0 1
--  An alpha spin H atom
0 2
F1
--  Another alpha spin H atom. Bonded ALMO-EDA does not need the multiplicities to sum to that of the molecule
0 2
F2 F1 1.3
$end

$rem
JOBTYPE            eda                       # use almo-eda
METHOD             wb97m-v                   # note: bonded eda examples only "EXCHANGE" instead of METHOD. Why?
BASIS              aug-cc-pVTZ              # chosen because PySCF uses 'full' by default
EDA2               10                        # EDA2 = 10 for bonded EDA. Use EDA2 = 2 otherwise
BONDED_EDA         1                         # (1) Perform ALMO-EDA with non-orthogonal CI. (2) Perform ALMO-EDA with spin-projected formalism, which is more efficient
SCF_CONVERGENCE    7                         # convergence in H = 10^SCF_CONVERGENCE
MAX_SCF_CYCLES     500
ROSCF              true                      # required for bonded EDA
SCF_GUESS          fragmo
DISP_FREE_X      revPBE                          # HF recommended for hybrid functionals (B97X-V), revPBE recommended for semi-local functionals (B97M-V)
DISP_FREE_C      PBE                        # put PBE if DISP_FREE_X is revPBE, NONE if DISP_FREE_X is set to an exchange-correlation functional
EDA_BSSE         false
SCF_ALGORITHM      diis
SCFMI_MODE         1                         # required for bonded EDA
SCF_PRINT_FRGM     true
FRZ_RELAX          true
FRZ_RELAX_METHOD   2
FRZ_ORTHO_DECOMP   1
$end

$rem_frgm
scf_convergence   8
scf_algorithm     diis
scf_guess         sad
$end

=====polarizing for fragment clarification=====
String list
        1        0
        0        1
Constructing subrot_ross from passed subrot_ross in u_ross PickSubrot
 Exchange:     0.1500 Hartree-Fock + 1.0000 wB97M-V + LR-HF
 Correlation:  1.0000 wB97M-V
 Using SG-2 standard quadrature grid
 Nonlocal Correlation:  VV10 with C = 0.0100 and b = 6.00 and scale = 1.00000
 Grid used for NLC:  SG-1 standard quadrature
constructed subrot_u_ross object from passed subrot_ross
Constructing subrot_ross from passed subrot_ross in u_ross PickSubrot
 Exchange:     0.1500 Hartree-Fock + 1.0000 wB97M-V + LR-HF
 Correlation:  1.0000 wB97M-V
 Using SG-2 standard quadrature grid
 Nonlocal Correlation:  VV10 with C = 0.0100 and b = 6.00 and scale = 1.00000
 Grid used for NLC:  SG-1 standard quadrature
constructed subrot_u_ross object from passed subrot_ross
ci vectors
  -0.7072   0.7070
   0.7072   0.7070
H(NOCI)
-213.8037102 -0.0344917
-0.0344917 -213.8037102
Ov(NOCI)
   1.0000   0.0002
   0.0002   1.0000
S2(NOCI)
   0.9998  -0.9998
  -0.9998   0.9998
Nuc e: 32.9718108
 ===================================================================
 Theory           State:          Energy (au) (kcal/mol)       <S2>
 -------------------------------------------------------------------
 SF-ALMO-NOCI         1:      -213.8095280827     0.0000     2.0000
 SF-ALMO-NOCI         2:      -213.7978945589     7.3001     0.0000
 ===================================================================

              Active Space Coefficients > 0.10
 -------------------------------------------------------------------
 State  1      Configuration (1 = SA, _ = SB)         Coeff.   <S^2>
 -------------------------------------------------------------------
                                        1 | _      -0.707107   1.000
                                        _ | 1       0.707107   1.000
 -------------------------------------------------------------------
 State  2      Configuration (1 = SA, _ = SB)         Coeff.   <S^2>
 -------------------------------------------------------------------
                                        1 | _       0.707107   1.000
                                        _ | 1       0.707107   1.000

<S^2> of spin-flipped wavefunction before polarization iterations 0.00000000000
Beginning iterations for the determination of the low-spin optimized wavefunction
Energy prior to optimization (guess energy) = -213.797894558943
Begin Timing: Total SCF Calculation
begin iterations for algorithm: DIIS         
the SCF tolerance is set to 1.00e-07
gen_scfman_exception: objective_function: get_vector is not implemented in inheriting class.

 Q-Chem fatal error occurred in module libgscf/gen_scfman/gen_scfman_main.C, line 246:

 Error in gen_scfman


 Please submit a crash report at q-chem.com/reporter

What are some ways I can improve the success of the run? I’ve changed the tolerances, SCF algorithms for both molecule and fragments, turned BSSE on/off, and tried both orthogonal CI and spin-projected formalisms. I’ve also tried it with cc-pVTZ and cc-pVQZ to remove diffuse functions but these also fail.

[1] Levine, D. S. and M. Head-Gordon (2017). “Energy decomposition analysis of single bonds within Kohn–Sham density functional theory.” Proceedings of the National Academy of Sciences of the United States of America 114(48): 12649-12656.

#2

This calculation should run fine with SCF_ALGORITHM set to GDM_LS instead of DIIS with Q-Chem 6.2

#3

Thanks for your response. Even after switching to version 6.2 and using GDM-LS, I’ve got errors. Here’s the input:

$molecule
   0 1
--  An alpha spin H atom
   0 2
   F1
--  Another alpha spin H atom. Bonded ALMO-EDA does not need the multiplicities to sum to that of the molecule
   0 2
   F2 F1 1.3
$end

$rem
   JOBTYPE            eda                       # use almo-eda
   METHOD             wb97m-v                   # note: bonded eda examples only "EXCHANGE" instead of METHOD. Why?
   BASIS              aug-cc-pVTZ              # chosen because PySCF uses 'full' by default
   EDA2               10                        # EDA2 = 10 for bonded EDA. Use EDA2 = 2 otherwise
   BONDED_EDA         1                         # (1) Perform ALMO-EDA with non-orthogonal CI. (2) Perform ALMO-EDA with spin-projected formalism, which is more efficient
   SCF_CONVERGENCE    7                         # convergence in H = 10^SCF_CONVERGENCE
   MAX_SCF_CYCLES     500
   ROSCF              true                      # required for bonded EDA
   SCF_GUESS          fragmo
   DISP_FREE_X      revPBE                          # HF recommended for hybrid functionals (B97X-V), revPBE recommended for semi-local functionals (B97M-V)
   DISP_FREE_C      PBE                        # put PBE if DISP_FREE_X is revPBE, NONE if DISP_FREE_X is set to an exchange-correlation functional
   EDA_BSSE         false
   SCF_ALGORITHM      gdm_ls
   SCFMI_MODE         1                         # required for bonded EDA
   SCF_PRINT_FRGM     true
   FRZ_RELAX          true
   FRZ_RELAX_METHOD   2
   FRZ_ORTHO_DECOMP   1
$end

$rem_frgm
   scf_convergence   8
   scf_algorithm     gdm_ls
   scf_guess         sad
$end

And the error message:

 Q-Chem warning in module libalmo/libalmo/scfmi/subrot_ci_u.C, line 3720:

 haven't worked out this preconditioner, using I instead

  497    -558.2624540621      1.67e-07     00000 BFGS Step
  498    -558.2624540618      1.05e-05     00000 LineSearch Step

 Q-Chem warning in module libalmo/libalmo/scfmi/subrot_ci_u.C, line 3720:

 haven't worked out this preconditioner, using I instead

  499    -558.2624540621      1.43e-07     00000 BFGS Step
  500    -558.2624540620      5.46e-06     00000 LineSearch Step
max iterations reached for algorithm: GDM_LS       
we have reached max iterations on our final algorithm

 Q-Chem fatal error occurred in module libgscf/gen_scfman/gen_scfman_hybrid_algorithm.C, line 402:

 SCF failed to converge


 Please submit a crash report at q-chem.com/reporter

Please advise, thanks!

#4

Hello lleung,

Thank you for the input/output files. Is there a reason that you would like to change the convergence criteria from (7, 6) → (8, 7) (12.7.2 ALMO-EDA for Bonded Interactions‣ 12.7 Additional ALMO-EDA Capabilities ‣ Chapter 12 Fragment-Based Methods ‣ Q-Chem 6.2 User’s Manual)? The solutions seems to oscillate for (8, 7). According to the bondded EDA section (see previous link), you may want to use SCF_ALGORITHM = L_BFGS for BONDED-EDA = 1. The following input should converge without any issues.

$molecule
   0 1
--  An alpha spin H atom
   0 2
   F1
--  Another alpha spin H atom. Bonded ALMO-EDA does not need the multiplicities to sum to that of the molecule
   0 2
   F2 F1 1.3
$end

$rem
   JOBTYPE            eda                       
   METHOD             wb97m-v                  
   BASIS              aug-cc-pVTZ              
   EDA2               10                        
   BONDED_EDA         1                        
   SCF_CONVERGENCE    6                         
   MAX_SCF_CYCLES     500
   ROSCF              true                    
   SCF_GUESS          fragmo
   DISP_FREE_X      revPBE                         
   DISP_FREE_C      PBE                       
   EDA_BSSE         false
   SCF_ALGORITHM      l_bfgs
   SCFMI_MODE         1                        
   SCF_PRINT_FRGM     true
   CHILD_MP           true
   CHILD_MP_ORDERS    1233
   FRZ_RELAX          true
   FRZ_RELAX_METHOD   2
   FRZ_ORTHO_DECOMP   1
   INTEGRAL_SYMMETRY    false
   POINT_GROUP_SYMMETRY false
$end

$rem_frgm
   scf_convergence   7
   scf_algorithm     l_bfgs
   scf_guess         autosad
$end

Note that I change fragment guess to autosad for faster convergence.

1 Like
#5

Actually, bonded_eda = 1 corresponds to the method based on non-orthogonal CI (https://doi.org/10.1021/acs.jctc.6b00571) according to the Q-Chem manual, which is not compatible with DFT methods. In my personal opinion, this method is obsolete and one should use bonded_eda = 2 (based on spin-projected DFT) instead.

I would use the following input instead:

$molecule
   0 1 
--  An alpha spin H atom
   0 2 
   F1  
--  Another alpha spin H atom. Bonded ALMO-EDA does not need the multiplicities to sum to that of the molecule
   0 2 
   F2 F1 1.3 
$end

$rem
   JOBTYPE            eda                       # use almo-eda
   METHOD             wb97m-v                   # note: bonded eda examples only "EXCHANGE" instead of METHOD. Why?
   BASIS              aug-cc-pVTZ              # chosen because PySCF uses 'full' by default
   EDA2               10                        # EDA2 = 10 for bonded EDA. Use EDA2 = 2 otherwise
   BONDED_EDA         2                         # (1) Perform ALMO-EDA with non-orthogonal CI. (2) Perform ALMO-EDA with spin-projected formalism, which is more efficient
   SCF_CONVERGENCE    7                         # convergence in H = 10^SCF_CONVERGENCE
   MAX_SCF_CYCLES     500 
   ROSCF              true                      # required for bonded EDA 
   SCF_GUESS          fragmo
   DISP_FREE_X      revPBE                          # HF recommended for hybrid functionals (B97X-V), revPBE recommended for semi-local functionals (B97M-V)
   DISP_FREE_C      PBE                        # put PBE if DISP_FREE_X is revPBE, NONE if DISP_FREE_X is set to an exchange-correlation functional
   EDA_BSSE         false
   SCF_ALGORITHM      gdm_ls
   SCFMI_MODE         1                         # required for bonded EDA 
   SCF_PRINT_FRGM     true
   FRZ_RELAX          true
   FRZ_RELAX_METHOD   2   
$end

$rem_frgm
   scf_convergence   8   
   scf_algorithm     gdm_ls
   scf_guess         sad 
$end
~
1 Like
#6

Thanks for the helpful suggestions regarding both the SCF tolerance and method. Will incorporate both.