EDA fatal error in some path

Questions
#1

I am attempting to run an EDA calculation on a relatively small molecule but have been receiving this error for all but one job:

Q-Chem fatal error occurred in module libalmo\libalmo\eda\eda2.C, line 872:

ortho fragment occupied subspace optimization failed to converge

Please submit a crash report at Q-Chem Crash Reporter
I am not sure where this module is located to investigate the error or how to fix it. The input file for this error is given here:
$molecule
0 1

0 1
C 1.492774000000 0.886198000000 -0.553983000000
C 1.979076000000 1.355305000000 0.660238000000
C 2.200056000000 2.724221000000 0.795290000000
C 1.961897000000 3.584854000000 -0.274929000000
C 1.509068000000 3.085906000000 -1.493403000000
C 1.286007000000 1.718895000000 -1.646765000000
H 2.162524000000 0.677079000000 1.491703000000
H 2.562443000000 3.117497000000 1.743653000000
H 2.143561000000 4.651938000000 -0.160769000000
H 1.335549000000 3.758082000000 -2.332044000000
H 0.958730000000 1.309709000000 -2.602263000000
S 2.365972000000 -1.097400000000 -1.427567000000
O 1.913809000000 -1.195625000000 -2.830977000000
O 2.610447000000 -2.331489000000 -0.655899000000
C 3.948522000000 -0.269675000000 -1.455845000000
H 3.867117000000 0.668077000000 -2.009004000000
H 4.278824000000 -0.107398000000 -0.426869000000
H 4.605774000000 -0.981356000000 -1.969105000000

0 1
Cu 0.304898000000 -0.654795000000 -0.526992000000
N -0.739355000000 -2.164521000000 -1.478981000000
H -0.667103000000 -1.862059000000 -2.450858000000
C -2.123191000000 -1.986766000000 -1.033673000000
H -2.202037000000 -2.406840000000 -0.021592000000
H -2.825366000000 -2.536324000000 -1.677358000000
C -2.477236000000 -0.516965000000 -1.009246000000
H -2.374137000000 -0.089665000000 -2.016298000000
H -3.526058000000 -0.386243000000 -0.703981000000
N -1.563078000000 0.200934000000 -0.122398000000
H -1.693633000000 -0.136235000000 0.831942000000
C -0.276307000000 -3.543555000000 -1.396478000000
H -0.958897000000 -4.229921000000 -1.918565000000
H 0.718461000000 -3.615172000000 -1.842074000000
H -0.209249000000 -3.844971000000 -0.345158000000
C -1.720236000000 1.644705000000 -0.149740000000
H -2.732909000000 1.952230000000 0.150228000000
H -0.994436000000 2.108210000000 0.527520000000
H -1.534139000000 2.016700000000 -1.162907000000
I 0.190167000000 -1.822834000000 2.379833000000
$end

$rem
JOBTYPE EDA
SCF_PRINT_FRGM TRUE
EDA2 1
EDA_PCT_A TRUE
EDA_CLS_DISP TRUE
FRGM_METHOD GIA
FRGM_LPCORR RS_EXACT_SCF
SYMMETRY FALSE
SYM_IGNORE TRUE
METHOD PW6B95
BASIS def2-SVP
DFT_D D3_BJ
EDA_COVP TRUE
EDA_PRINT_COVP 2
MAKE_CUBE_FILES TRUE
PLOTS TRUE
MAX_SCF_CYCLES 200
SCF_ALGORITHM DIIS_GDM
MAX_SCF_CYCLES 150
MEM_TOTAL 30000
MEM_STATIC 2000
SOLVENT_METHOD smd
$end

$plots
grid_points 100 100 100
$end

$smx
solvent acetonitrile
$end
I did have one successful job, similar in structure to the one above, but I am not sure why this one job was able to finish with normal termination. Here is the input file for the successful one:
$molecule
0 1

0 1
C 1.502625720900 0.856278098300 -0.536104639000
C 1.977992003400 1.329122086000 0.681761888800
C 2.186961752000 2.699335149200 0.819572826200
C 1.944858317100 3.561151524600 -0.249687563900
C 1.497843861200 3.060378739800 -1.469231261300
C 1.289412812100 1.691050548900 -1.625665175600
H 2.141007548000 0.645048006500 1.512638957400
H 2.539567754600 3.095289811500 1.770760683900
H 2.115438878300 4.629811920800 -0.132388573500
H 1.317499388800 3.732685840300 -2.306512700000
H 0.970552218900 1.278138882200 -2.582201112600
S 2.366243934200 -1.093964620400 -1.441915502000
O 1.898418228600 -1.142285799400 -2.844615178600
O 2.634713873600 -2.359213971600 -0.732917275900
C 3.952173668100 -0.266105707100 -1.464835725500
H 3.868008538300 0.688124957900 -1.988648510100
H 4.293735094100 -0.133105141200 -0.435189541900
H 4.604664327700 -0.962211436500 -2.004397036400

0 1
Cu 0.331353968100 -0.699520794400 -0.453797380900
N -0.742877647300 -2.153219705300 -1.492523679600
H -0.686529317200 -1.857494903200 -2.466332348300
C -2.111748160800 -1.977279335200 -1.011285792000
H -2.144673595700 -2.368746728600 0.015024088100
H -2.830663833500 -2.550792559000 -1.615862690600
C -2.483226149100 -0.511079101200 -1.008979141000
H -2.407701741300 -0.105071924800 -2.027723222100
H -3.529809694500 -0.389486008400 -0.688277084500
N -1.561927895800 0.234188651600 -0.159917397100
H -1.638613987200 -0.096645236900 0.803187055200
C -0.270500619800 -3.525529780400 -1.386277276400
H -0.972505258500 -4.236867770800 -1.847595803500
H 0.702469773200 -3.609390458500 -1.876092519700
H -0.144274153800 -3.773398486000 -0.326042086900
C -1.730384693000 1.672410440700 -0.209338835900
H -2.739118167300 1.985229614100 0.102440913800
H -0.995006225300 2.153343383000 0.445852620800
H -1.563466634200 2.030481597300 -1.231503990300
Br 0.213539063300 -1.669283393500 1.899745220900
$end

$rem
JOBTYPE EDA
SCF_PRINT_FRGM TRUE
EDA2 1
EDA_PCT_A TRUE
EDA_CLS_DISP TRUE
FRGM_METHOD GIA
FRGM_LPCORR RS_EXACT_SCF
SYMMETRY FALSE
SYM_IGNORE TRUE
METHOD PW6B95
BASIS def2-SVP
DFT_D D3_BJ
EDA_COVP TRUE
EDA_PRINT_COVP 2
MAKE_CUBE_FILES TRUE
PLOTS TRUE
MAX_SCF_CYCLES 200
SCF_ALGORITHM DIIS
MAX_SCF_CYCLES 150
MEM_TOTAL 30000
MEM_STATIC 2000
SOLVENT_METHOD smd
$end

$plots
grid_points 100 100 100
$end

$smx
solvent acetonitrile
$end

#2

When using def2-SVP basis for a system containing the iodine atom, you will need to add “ECP def2-ECP” in the $rem section.

#3

Hello Yuezhi and thanks for the response,

I have attempted to run the same molecules with def2-ecp and even changed my model chemistry. Now my issues are these disparate energy values for both my iodine molecules. I am able to get reasonable values when running the calculations with basis set 6-311g** but the energy in def2-svp are exaggerated.

In addition, I am attempting to run the previously optimized, in Gaussian 16, starting materials in Q-chem but now imaginary frequencies are populated in my aryl iodine geometry. These calculations are using the same model chemistry from my EDA calculations. I am currently adjusting the hessian algorithm from the default one and also varying the method to see if maybe there is an issue with pw6b95 describing the iodine atom.

Has anyone dealt with this problem and if so how did you resolve the issue?

Thanks

#4

6-311G** has all-electron basis functions for iodine, so your first problem sounds to me like there’s an ECP issue with the EDA code (?). If the absolute Hartree energy is wildly off, it could be that the EDA code is not adjusting the nuclear charges to reflect the electrons that are deleted when the ECP is added - I have seen that behavior before, in other parts of the code that hadn’t yet been made compliant with ECPs.

Not sure about 2nd problem. Note that default XC quadrature grids for different functionals are different in Q-Chem vs Gaussian so you do need to re-optimize, even when the model chemistry (functional / basis set) is the same.

#5

I am not sure if the ECP is causing the error. I have a Pople (6-311G**) with lanl2dz and an Ahlrich (def2-SVP) with lanl2dz. Both jobs terminate normally but my Pople value gives -31.9 kcal/mol for EDA total versus -3000 kcal/mol for Ahlrich. For a similar structure with Bromine, my values are -29.6 and -31.8 kcal/mol for Pople and Ahlrich respectively.

The optimization for the Iodine structure is now giving me an error related to the hessian update. The error message can be found below:

ERROR Hessian Appears to have all zero or negative eigenvalues

Q-Chem fatal error occurred in module 0, line 9:

OPTIMIZE fatal error

I also ran the fragments, the reactant geometries, as optimization on qchem. After optimization, the frequency calculation determines that my bromine reactants are minima but my aryl Iodine structure have two imaginary frequencies. I think that the Iodine atom is not being described but I would like a second opinion.

Thanks

#6

That’s consistent with a problem in the EDA code for ECPs, because Ahlrichs basis sets do not use an ECP for Br but do for I. Therefore both def2-SVP and 6-311G** are all-electron for Br (and seem to behave comparably), but for iodine the Pople basis is all-electron but the Ahlrichs basis set uses an ECP. Feels like a bug that you should report to Q-Chem support.

#7

I did a test calculation for CF3I…Cl- with def2-SVP/def2-ECP and printed out the nuclear charge for iodine in the function that EDA2 code uses to compute electrostatic interactions. The result seems correct (28 electrons were properly subtracted).

I’m a little confused by this description: “I have a Pople (6-311G**) with lanl2dz and an Ahlrich (def2-SVP) with lanl2dz”. Since 6-311G** is all-electron for iodine, lanl2dz should not be invoked. While for def2-SVP, the ECP to use should not be lanl2dz either.