SCF energy gradient gives NAN values for organometallic compound

Dear Q-Chem community,

I am attempting to optimize the structure of an organometallic compounds containing Zr, O, C, and H using the BP86+D3/def2-SVP method with effective core potentials (def2-ecp). Without the D3 corrections my structure converges well. With the D3 method, however, the first SCF cycle converges but the analytical gradient of the SCF energy reports “-nan” values and the optimization fails. Can anyone help me to understand why the D3 method might cause this error, and what I might do to resolve this issue? I would greatly appreciate any suggestions!

A complete input file will be required in order to diagnose this. Please use the formatted text box </> to avoid losing special characters. It is best if you can find a “minimal example” that illustrates the problem, stripping away any unnecessary keywords.

$molecule
1 1
  H    0.2120610    3.2282449    3.0585973
  H   -2.9049710   -1.4198724    3.0723946
  C    0.1770789    2.6608870    2.0963942
  C   -2.3935563   -1.1770121    2.1091306
  O   -0.9693456    2.4622594    1.5862972
  O    1.2883650    2.2946394    1.6031293
  O   -1.6380731   -2.0712838    1.6168961
  O   -2.6400812   -0.0406933    1.5977379
 Zr   -1.6456562    1.0973507    0.0080725
 Zr    1.7813494    0.8685712    0.0104477
 Zr   -0.1289242   -1.9876044    0.0227133
  O    0.0031180   -0.0058661    0.7124107
  O    0.1196316    1.7603047   -1.0663431
  H    0.1754953    2.4825208   -1.7146750
  O   -1.5865918   -0.7956702   -1.0521337
  H   -2.2227191   -1.1138609   -1.7147017
  C   -2.9194207    1.9936849   -2.0150017
  C   -2.6055122    3.1362900   -1.2163871
  C   -3.3442929    3.0427473    0.0004953
  C   -4.1175944    1.8477699   -0.0433548
  C   -3.8586010    1.2003343   -1.2875894
  O    1.4890124   -1.0022862   -1.0450477
  H    2.0774711   -1.3964367   -1.7111739
  O    2.6138249   -0.3879645    1.6081907
  C    3.1489666    1.6161845   -2.0124468
  C    3.9876469    0.7069928   -1.2976805
  C    4.3293935    1.3064488   -0.0489312
  C    3.7089346    2.5867812    0.0094411
  C    2.9809620    2.7808761   -1.2011479
  O    1.3538384   -2.2697734    1.6131813
  C   -0.2121088   -3.5436645   -1.9938685
  C   -1.4042990   -3.8105268   -1.2530790
  C   -1.0312414   -4.4092554   -0.0130414
  C    0.3884775   -4.5183109    0.0145398
  C    0.8960168   -3.9865998   -1.2074801
  H   -2.5336375    1.7801114   -3.0213500
  H   -1.9384869    3.9631153   -1.4964567
  H   -3.3103545    3.7578311    0.8319065
  H   -4.7853186    1.4840494    0.7475047
  H   -4.3327001    0.2714939   -1.6332598
  C    2.2202794   -1.4855464    2.1110846
  H    2.7366948    1.4605369   -3.0189906
  H    4.3396693   -0.2708035   -1.6541183
  H    4.9527209    0.8569050    0.7340689
  H    3.7676840    3.2927388    0.8473018
  H    2.4184058    3.6857715   -1.4682778
  H   -0.1602128   -3.1093455   -3.0018594
  H   -2.4341631   -3.6214516   -1.5857619
  H   -1.7157869   -4.7172482    0.7870857
  H    0.9867424   -4.9283429    0.8378938
  H    1.9547065   -3.9574768   -1.4990952
  H    2.6948858   -1.8008840    3.0721471
$end

$rem
   BASIS = def2-SVP
   ECP = def2-ecp
   JOB_TYPE = Optimization
   METHOD = BP86
   DFT_D = D3
$end

Thank you for your help! I have confirmed that the above minimal example causes the same issue–the SCF energy gradient gives all “-nan” values.

Hi Ari,
(1) I can reproduce this error with the latest version of Q-Chem (pre-release v. 6.3).
(2) It appears to be problem with the dispersion correction, because commenting out DFT_D removes the problem.
(3) It’s specifically a problem for Zr, because replacing Zr with Ti removes the problem.
(3) Those changes isolate the problem but it’s not what you want. As a workaround, I found that DFT_D = D4 seems to work (at least through one optimization step, no NaN).

I will submit a bug ticket about D3 (which affects both the D3(0) and D3(BJ) versions, I checked both).

I just checked and the D4 method worked well. Thank you for helping with this!