Errors with Hirshfeld analysis

joannaz · March 20, 2024, 8:32pm

Hello!

I’ve been trying to calculate iterative Hirshfeld population analysis, but I can’t seem to get it right, I always encounter some error.
Earlier in my calculations I’ve been using my basis set, but when I tried using that here, the initial SCF (in the small basis set) didn’t converge.
I thought that either the basis set, or the scf algorithm may have been an issue, so I tried changing both these options.

Change of scf algorithm

input:
$molecule
1 1
C 0.0258960 0.2641517 4.4241795
C 0.5731417 0.3397343 3.1524444
C 1.1789075 1.5074027 2.6842218
C 1.2474938 2.6082257 3.5365011
C 0.7083110 2.5465084 4.8154006
C 0.0955875 1.3782768 5.2524793
O 0.5816658 -0.7266037 2.2771752
C -0.4668564 -1.6269225 2.2902853
C -1.6668918 -1.2868826 1.6603040
C -2.6795458 -2.2460435 1.6198420
C -2.4935740 -3.5013451 2.1854606
C -1.2879083 -3.8176386 2.7997627
C -0.2670310 -2.8762879 2.8538145
P -1.8112513 0.3327653 0.8206085
Cu -0.0381248 1.0203589 -0.4448364
N -0.5566158 2.3401630 -1.9672813
N -0.9750547 3.5529921 -2.0845767
N -1.1037814 3.7772164 -3.3860643
C -0.7670998 2.7009901 -4.1293286
C -0.4087852 1.7604311 -3.1949169
H -1.4252914 4.6877947 -3.6988487
C 0.0604180 0.3755233 -3.2798888
N 0.3363839 -0.1936124 -2.1192049
C 0.7742588 -1.4835477 -2.0897407
C 0.9307928 -2.2344525 -3.2840557
C 0.6294628 -1.5950645 -4.5083746
C 0.1993744 -0.2981007 -4.5122631
C 1.3748332 -3.5761512 -3.1993926
C 1.6487082 -4.1376504 -1.9816210
C 1.4975625 -3.3815391 -0.7967823
C 1.0737945 -2.0813011 -0.8448385
P 1.6788263 1.5586298 0.9316819
H 0.9626357 -1.4860150 0.0590112
H 1.7210182 -3.8402860 0.1659449
H 1.9880474 -5.1710107 -1.9204992
H 1.4930340 -4.1514582 -4.1178739
H 0.7427272 -2.1457400 -5.4427106
H -0.0387183 0.2117850 -5.4439759
H -0.8081272 2.6987488 -5.2123181
H -3.1036289 0.1780650 0.2655092
H -2.1723036 1.1795114 1.8989563
H 0.6901122 -3.0948398 3.3261046
H -1.1397816 -4.8024247 3.2417760
H -3.2957322 -4.2370265 2.1422693
H -3.6263364 -2.0117999 1.1329914
H 1.7178515 3.5295717 3.1922962
H 0.7656228 3.4142728 5.4709673
H -0.3297674 1.3263406 6.2542664
H -0.4450057 -0.6571740 4.7632395
H 2.8523852 0.7643310 0.9194434
H 2.3009238 2.8292296 0.9005737
$end

$rem
jobtype sp
method camb3lyp
dft_d d3_bj
basis gen
basis2 sto-3g
purecart 11
PRINT_GENERAL_BASIS false
gen_scfman false
MAX_SCF_CYCLES 500
SCF_CONVERGENCE 9
scf_final_print true
scf_guess autosad
MEM_STATIC 2000
MEM_TOTAL 80000
hirshiter true
$end

(I’m omitting the $basis section)
Similarly to earlier case, the SCF didn’t converge (withing 500 iterations).

Change of the basis set

$molecule part is the same as above

$rem
jobtype sp
method camb3lyp
dft_d d3_bj
basis def2-svp
purecart 11
PRINT_GENERAL_BASIS false
gen_scfman false
MAX_SCF_CYCLES 500
SCF_CONVERGENCE 9
scf_final_print true
MEM_STATIC 2000
MEM_TOTAL 80000
hirshiter true
$end

This input yielded the following error:
"Q-Chem fatal error occurred in module gesman/frgmchild.C, line 676:

Error executing Q-Chem for fragment 19"

Change of both scf algorithm and basis set

One more version I tried:
$rem
jobtype sp
method camb3lyp
dft_d d3_bj
basis def2-svp
purecart 11
PRINT_GENERAL_BASIS false
gen_scfman true
MAX_SCF_CYCLES 500
SCF_CONVERGENCE 9
scf_final_print true
scf_guess sap
MEM_STATIC 2000
MEM_TOTAL 80000
hirshiter true
$end

This resulted in a different error:
"Performing Hirshfeld-I population analysis.

Q-Chem fatal error occurred in module libmdc/newfileman.C, line 376:

FileMan error: End of file reached prematurely reading (204) bytes in file UNKNOWN FILE
Path: /net/ascratch/people/plgjoannaz/slurm_jobdir/8128024/tmp/qchem2397561/1621.0"

Could someone please explain what’s going wrong in these cases here? And what can I do to make it work?

I’ll be grateful for any help.

jherbert · March 20, 2024, 9:06pm

Does this SCF calculation converge in the absence of the Hirshfeld keywords? Sounds like you are using an exotic basis set so it’s not clear that it does; also for molecules of this size, I would suggest setting THRESH = 12.

joannaz · March 20, 2024, 9:19pm

Yes, it does, in less then 30 SCF iterations. I’m using DZP basis set rewritten in Q-Chem format, as it was unavailable.

jherbert · March 21, 2024, 12:02am

In that case I can’t really diagnose without a complete input file. Please use the preformatted text buttom, </>, otherwise this GUI will eat the special characters in your $basis section.

joannaz · March 21, 2024, 11:05am

Here is my original input:

$molecule
1 1
C     0.0258960    0.2641517    4.4241795 
C     0.5731417    0.3397343    3.1524444 
C     1.1789075    1.5074027    2.6842218 
C     1.2474938    2.6082257    3.5365011 
C     0.7083110    2.5465084    4.8154006 
C     0.0955875    1.3782768    5.2524793 
O     0.5816658   -0.7266037    2.2771752 
C    -0.4668564   -1.6269225    2.2902853 
C    -1.6668918   -1.2868826    1.6603040 
C    -2.6795458   -2.2460435    1.6198420 
C    -2.4935740   -3.5013451    2.1854606 
C    -1.2879083   -3.8176386    2.7997627 
C    -0.2670310   -2.8762879    2.8538145 
P    -1.8112513    0.3327653    0.8206085 
Cu   -0.0381248    1.0203589   -0.4448364 
N    -0.5566158    2.3401630   -1.9672813 
N    -0.9750547    3.5529921   -2.0845767 
N    -1.1037814    3.7772164   -3.3860643 
C    -0.7670998    2.7009901   -4.1293286 
C    -0.4087852    1.7604311   -3.1949169 
H    -1.4252914    4.6877947   -3.6988487 
C     0.0604180    0.3755233   -3.2798888 
N     0.3363839   -0.1936124   -2.1192049 
C     0.7742588   -1.4835477   -2.0897407 
C     0.9307928   -2.2344525   -3.2840557 
C     0.6294628   -1.5950645   -4.5083746 
C     0.1993744   -0.2981007   -4.5122631 
C     1.3748332   -3.5761512   -3.1993926 
C     1.6487082   -4.1376504   -1.9816210 
C     1.4975625   -3.3815391   -0.7967823 
C     1.0737945   -2.0813011   -0.8448385 
P     1.6788263    1.5586298    0.9316819 
H     0.9626357   -1.4860150    0.0590112 
H     1.7210182   -3.8402860    0.1659449 
H     1.9880474   -5.1710107   -1.9204992 
H     1.4930340   -4.1514582   -4.1178739 
H     0.7427272   -2.1457400   -5.4427106 
H    -0.0387183    0.2117850   -5.4439759 
H    -0.8081272    2.6987488   -5.2123181 
H    -3.1036289    0.1780650    0.2655092 
H    -2.1723036    1.1795114    1.8989563 
H     0.6901122   -3.0948398    3.3261046 
H    -1.1397816   -4.8024247    3.2417760 
H    -3.2957322   -4.2370265    2.1422693 
H    -3.6263364   -2.0117999    1.1329914 
H     1.7178515    3.5295717    3.1922962 
H     0.7656228    3.4142728    5.4709673 
H    -0.3297674    1.3263406    6.2542664 
H    -0.4450057   -0.6571740    4.7632395 
H     2.8523852    0.7643310    0.9194434 
H     2.3009238    2.8292296    0.9005737 
$end

$rem
jobtype            sp 
method             camb3lyp 
dft_d              d3_bj
basis              gen
basis2             sto-3g
purecart           11 
PRINT_GENERAL_BASIS false
gen_scfman         false 
MAX_SCF_CYCLES     500
SCF_CONVERGENCE    9
scf_final_print    true
MEM_STATIC         2000
MEM_TOTAL          80000
hirshiter          true
$end

$basis
cu 0
   s  6  1.0 
  441087.25070     -0.21814173104E-03
  66112.021187     -0.16921935882E-02
  15047.011425     -0.88137131631E-02
  4263.4273084     -0.35954659517E-01
  1396.3815797     -0.11742970473
  511.96055788     -0.28844267427
   s  2  1.0
  203.45426948     -0.42678898925
  82.792337027     -0.33044128545
   s  1  1.0
  20.854285634       1.0000000000
   s  1  1.0
  9.0410679584       1.0000000000
   s  1  1.0
  2.7518135173       1.0000000000
   s  1  1.0
  1.0434856515       1.0000000000
   s  1  1.0
 0.11172292442       1.0000000000
   s  1  1.0
 0.41041020330E-01   1.0000000000
   p  3  1.0
  2530.0965671      0.19137941259E-02
  600.09792954      0.15797670964E-01
  194.08204479      0.76271259615E-01
   p  3  1.0
  73.671821377      0.23881452309
  30.447369690      0.44980015897
  13.122714875      0.39337682391
   p  1  1.0
  5.5214839972       1.0000000000
   p  1  1.0
  2.1457922130       1.0000000000
   p  1  1.0
 0.76797488702       1.0000000000
   d  3  1.0
  47.313743703      0.32399760431E-01
  13.154688449      0.16822706530
  4.3662885749      0.38494429603
   d  1  1.0
  1.4122065936       1.0000000000
   d  1  1.0
 0.38840713001       1.0000000000
   p  1  1.0
 0.15506500000       1.0000000000
****
p 0
   s  5  1.0
  43629.146384      0.60290694324E-03
  6544.9101761      0.46551468038E-02
  1489.8784631      0.23786703379E-01
  422.38300193      0.90866640549E-01
  138.69747966      0.24743794714
   s  2  1.0
  50.865386699      0.45854747443
  19.826936819      0.30750563013
   s  1  1.0
  4.5933029207       1.0000000000
   s  1  1.0
  1.7570187795       1.0000000000
   s  1  1.0
 0.34358804778       1.0000000000
   s  1  1.0
 0.12451452677       1.0000000000
   p  4  1.0
  215.97067949      0.99703587825E-02
  50.428986942      0.70614266037E-01
  15.653604009      0.25773443022
  5.4085976775      0.49153067625
   p  1  1.0
  1.9189363936       1.0000000000
   p  1  1.0
 0.45420822084       1.0000000000
   p  1  1.0
 0.13190937083       1.0000000000
   d  1  1.0
 0.45000000000       1.0000000000
****
h 0
   s  3  1.0
  13.010701000      0.19682158000E-01
  1.9622572000      0.13796524000
 0.44453796000      0.47831935000
   s  1  1.0
 0.12194962000       1.0000000000
   p  1  1.0
 0.80000000000       1.0000000000
****
o 0 
   s  5  1.0
  6773.3747000      0.17266130000E-02
  1016.7970000      0.13246484000E-01
  231.26738000      0.66026157000E-01
  65.008454000      0.23528361000
  20.499808000      0.54976236000
   s  1  1.0
  6.7160664000       1.0000000000
   s  1  1.0
  1.1471572000       1.0000000000
   s  1  1.0
 0.33423251000       1.0000000000
   p  3  1.0
  17.694264000      0.43590224000E-01
  3.8536027000      0.23178087000
  1.0467465000      0.51455969000
   p  1  1.0
 0.27586043000       1.0000000000
   d  1  1.0
  1.2000000000       1.0000000000

****
n 0
   s  5  1.0
  5071.9892000      0.17067238000E-02
  761.41791000      0.13087696000E-01
  173.18418000      0.65098252000E-01
  48.670390000      0.23051603000
  15.331448000      0.53300474000
   s  1  1.0
  5.0186710000       1.0000000000
   s  1  1.0
 0.83554772000       1.0000000000
   s  1  1.0
 0.24626814000       1.0000000000
   p  3  1.0
  13.550722000      0.40630840000E-01
  2.9178682000      0.22085799000
 0.79831252000      0.51860142000
   p  1  1.0
 0.21900125000       1.0000000000
   d  1  1.0
  1.0000000000       1.0000000000

****
f 0
   s  5  1.0
  8709.5462000      0.13935739000E-02
  1307.4134000      0.10695073000E-01
  297.36441000      0.53389308000E-01
  83.602154000      0.19112479000
  26.385121000      0.44987537000
   s  1  1.0
  8.6499174000       1.0000000000
   s  1  1.0
  1.5043291000       1.0000000000
   s  1  1.0
 0.43408405000       1.0000000000
   p  3  1.0
  22.665994000      0.34028872000E-01
  4.9760341000      0.17865315000
  1.3476543000      0.38497323000
   p  1  1.0
 0.34775168000       1.0000000000
   d  1  1.0
  1.4000000000       1.0000000000

****
b 0
   s  5  1.0
  2410.2061000      0.17549915000E-02
  361.86992000      0.13436405000E-01
  82.308593000      0.66369025000E-01
  23.109942000      0.23025297000
  7.2517259000      0.51498212000
   s  1  1.0
  2.3666469000       1.0000000000
   s  1  1.0
 0.35987078000       1.0000000000
   s  1  1.0
 0.11162889000       1.0000000000
   p  3  1.0
  6.0017073000     -0.35545320000E-01
  1.2401097000     -0.19834505000
 0.33668024000     -0.50478729000
   p  1  1.0
 0.95590862000E-01   1.0000000000
   d  1  1.0
 0.50000000000       1.0000000000
****
c 0
   s  5  1.0
  3623.8613000      0.16339191000E-02
  544.04621000      0.12521701000E-01
  123.74338000      0.62113914000E-01
  34.763209000      0.21817729000
  10.933333000      0.49800431000
   s  1  1.0
  3.5744765000       1.0000000000
   s  1  1.0
 0.57483245000       1.0000000000
   s  1  1.0
 0.17303640000       1.0000000000
   p  3  1.0
  9.4432819000      0.37895451000E-01
  2.0017986000      0.20818177000
 0.54629718000      0.50474166000
   p  1  1.0
 0.15202684000       1.0000000000
   d  1  1.0
 0.80000000000       1.0000000000
****
$end

kcarter-fenk · March 21, 2024, 4:53pm

Regarding the job where the SCF converged, but iterative Hirshfeld failed: I can reproduce this by setting SCF_GUESS = SAP (as in your input file). File 1621 is written only in the SCF_GUESS = AUTOSAD routines. Any time Hirshfeld charges or their analogues such as CM5 or Hirshfeld-I are desired, AUTOSAD guess is required as it generates the neutral atomic densities for stockholder partitioning.

It at least used to be the case that Hirshfeld and related atomic charges would force the calculation path through the AUTOSAD routine, and then replace the AUTOSAD density with the user-requested guess if it differs from AUTOSAD. I will take a look at this and either restore this behavior or otherwise figure out why it isn’t doing this for Hirshfeld-I charges.

My suggestion for now would be to run your alternate basis set calculation (Def2-SVP) with the AUTOSAD guess since that will likely converge and write the required files to disk for Hirshfeld-I.

If AUTOSAD guess isn’t converging your SCF, a hacky fix would be to run your exact input with an AUTOSAD guess and SCF_MAX_CYCLES = 1, then in a second job (in the same input file, delimited by @@@) just run the same exact calculation with your preferred guess (just not READ because the MOs written to disk will be a bad guess). This way, the AUTOSAD files will be in scratch and your second calculation will be able to read them. Just be careful not to change details of the method or basis set between the two jobs.

joannaz · March 22, 2024, 10:38am

Okay, so I ran the calculations as you suggested, and once again I got the error:
“Q-Chem fatal error occurred in module gesman/frgmchild.C, line 676:
Error executing Q-Chem for fragment 19”

After rerunning it with “SCF_PRINT_FRGM = true” to obtain more information, I obtained the following:

"Q-Chem fatal error occurred in module qparser/MoleculeInput.C, line 805:
Invalid charge/multiplicity combination in MoleculeInput::getNElectrons!
Please submit a crash report at Q-Chem Crash Reporter

Q-Chem fatal error occurred in module gesman/frgmchild.C, line 676:
Error executing Q-Chem for fragment 19"

So could it maybe be a bug in the program?

The second way that you suggested unfortunately doesn’t work either - since SCF doesn’t converge within that one step, Q-Chem crashes.

kcarter-fenk · March 22, 2024, 11:39am

Which version of Q-Chem are you using? This was a bug with iterative Hirshfeld on transition metal atoms until Q-Chem v6.0.2.

joannaz · March 25, 2024, 2:05pm

Thank you! That was the reason. I had been using Q-Chem 6.0.1 instead of 6.1. After switching to the newer version it works fine, as long as I’m using built-in basis sets. Hance my further questions:

Is there anything I can do to make it work with my own basis set?
More importantly, can I calculate Hirshfeld population analysis in the excited state? I tried adding CIS_STATE_DERIV=1 along with other TD-DFT controlling keywords, but I also got pop. analysis for the ground state.

jherbert · March 25, 2024, 4:49pm

If it doesn’t work with BASIS=GEN, you might be able to get it to work via a user-defined basis set (BASIS=USER1 or USER2, then you need to define user1.bas in $QCAUX/basis). This is described in the manual, in the chapter on basis sets.

For excited-state Hirshfeld, it probably doesn’t exist. You could contact Q-Chem support with a feature request but no guarantees on the timescale.

joannaz · March 26, 2024, 10:06pm

Okay, thank you very much!

On an unrelated to Hirshfeld note, some time ago I had asked a question about a weird behaviour of difference densities in ROKS (Unexpected ROKS difference density - #4 by joannaz) and I didn’t get any answer. Do either of you perhaps have any idea what could be going on there?

Thank you in advance.