Crash during QMMM calculations

Hello,
Recently I was doing QMMM single point calculations (version of Q-Chem: 6.2). My calculations consist of 3 consecutive jobs (3 inputs in one file): first job was ground state calculation with PCM, second one was ROKS with PCM and the third one was also a ground state but with frozen solvent field (I used rf_ptss_save in ROKS job). Calculations crashed near the beginning of SCF at 3rd job:

 There are 1240 shells and 3015 basis functions
 Adding 8412 external point charges to one-electron Hamiltonian
 Nucleus-charge energy    =    16.8387980807 hartrees
 Charge-charge energy     =     0.0000000000 hartrees
 A cutoff of  1.0D-12 yielded 215867 shell pairs
 There are   1122352 function pairs
 Smallest overlap matrix eigenvalue = 1.87E-06
 Guess MOs from SCF MO coefficient file
 Reading MOs from coefficient file
 Reading MOs from coefficient file
 Long-range K will be added via erf
 Coulomb attenuation parameter = 0.33 bohr**(-1)
 A restricted hybrid HF-DFT SCF calculation will be
 performed using Pulay DIIS extrapolation
 Polarizable Continuum solvation model will be applied
 Exchange:     0.1900 Hartree-Fock + 0.3500 B88 + 0.4600 muB88 + LR-HF
 Correlation:  0.1900 VWN5 + 0.8100 LYP
 Using SG-3 standard quadrature grid
 Dispersion:   Grimme D3
 SCF converges when DIIS error is below 1.0E-08
 Using Q-Chem read-in guess as SCF_GUESS READ specified.

 C-PCM solvent model [f = (eps-1)/eps], solve by conjugate gradient
 using 48 threads for integral computing
 -------------------------------------------------------
 OpenMP Integral computing Module                
 -------------------------------------------------------
 ---------------------------------------
  Cycle       Energy         DIIS Error
 ---------------------------------------
    1   -5360.6077490752      2.29E-04
### Frozen reaction-field calculation using stored ESP to compute ASC ###
Reading ESP for state 0 (0 is MP-GS) from file...
... done. Computing ASC and self-interaction...
... done. Updating ref-RF self-interaction (-11.1539 eV) ...
... done. Separating ASC into fast and slow part...

 Q-Chem fatal error occurred in module include/RefCount.h, line 48:

 RefCount::Constructor called with null pointer argument

I want to ask what do you think, what could be causing this crash?

My input: https://sharetext.io/2a81b0c4

Sincerely,
Tiopirfur

I have tried to look into this but it would be helpful if you could find an example that’s not so large.

Hello,
I’m searching for smaller example and I’ll let you know if I find one.
My system consists of acetonitrile molecules and solute molecule. I tried to repeat calculations with only 2 acetonitrile molecules (one QM and one MM) and solute. These calculations completed without crash, but during SCF in 3rd job I got this warning:

 ---------------------------------------
    1    -995.0625872374      5.53E-04
### Frozen reaction-field calculation using stored ESP to compute ASC ###
Reading ESP for state 0 (0 is MP-GS) from file...
... done. Computing ASC and self-interaction...
... done. Updating ref-RF self-interaction (-2.2007 eV) ...
... done. Separating ASC into fast and slow part...
... done. Computing ESP of SCF density and SCF-fRF interaction ...
... done. Initial SCF-fRF interaction is -1.6214 eV.
Now carrying out SCF calculation in frozen ASC:
Adding current polarization work from ref-RF self-interaction (1.1004 eV)
    2    -995.7223315939      6.48E-04
    3    -995.7760577175      7.97E-05
    4    -995.7767451445      4.01E-05
    5    -995.7769512443      9.88E-06
    6    -995.7769656611      5.85E-06
    7    -995.7769707209      2.23E-06
    8    -995.7769714761      8.61E-07
    9    -995.7769716147      3.85E-07
   10    -995.7769716427      1.53E-07
   11    -995.7769716488      5.30E-08
   12    -995.7769716496      2.61E-08
   13    -995.7769716497      1.17E-08

 Q-Chem warning in module libscrf/libscrf/pcm/PCM_CLASS/pcm_setman_ss.C, line 706:

 SS-PCM is not implemented for CG solver, sorry ~JMM


### RF-ptSS: Non-eq. corr. for transition from EqS-Ref state ####
    => ptSS (<i0|Ri^f - R0^f|i0>) term:    0.000000 eV           
   14    -995.7769716501      4.74E-09 Convergence criterion met
 ---------------------------------------

I don’t know if this warning is related with my crash, but I thought it’s worth sharing.

Input for these calculations:

$molecule
0 1
O   2.114000082015991   -1.9329999685287476   -0.7620000243186951   125   3   5   0   0
N   -3.7760000228881836   -0.17499999701976776   0.5339999794960022   79   25   31   35   0
C   1.621000051498413   -0.6340000033378601   -1.1119999885559082   37   1   15   16   30
O   3.575000047683716   -2.8980000019073486   0.6850000023841858   124   5   0   0   0
C   3.0829999446868896   -1.843999981880188   0.26100000739097595   123   1   4   6   0
C   3.2780001163482666   -0.492000013589859   0.5659999847412109   13   5   7   15   0
C   4.1519999504089355   0.15000000596046448   1.4900000095367432   11   6   8   9   0
H   4.804999828338623   -0.4309999942779541   2.127000093460083   12   7   0   0   0
C   4.145999908447266   1.5260000228881836   1.5479999780654907   11   7   10   11   0
H   4.807000160217285   2.0230000019073486   2.246999979019165   12   9   0   0   0
C   3.315999984741211   2.306999921798706   0.7319999933242798   11   9   12   13   0
H   3.3389999866485596   3.384999990463257   0.8019999861717224   12   11   0   0   0
C   2.444999933242798   1.6710000038146973   -0.18199999630451202   11   11   14   15   0
H   1.7949999570846558   2.2639999389648438   -0.8130000233650208   12   13   0   0   0
C   2.440999984741211   0.3109999895095825   -0.24899999797344208   13   3   6   13   0
C   1.8329999446868896   -0.4320000112056732   -2.6050000190734863   1   3   17   18   36
H   2.888000011444092   -0.5759999752044678   -2.8259999752044678   6   16   0   0   0
H   1.5549999475479126   0.5740000009536743   -2.9100000858306885   6   16   0   0   0
H   -4.265999794006348   -2.1510000228881836   1.0429999828338623   18   35   0   0   0
H   -5.421000003814697   -0.9160000085830688   1.562999963760376   18   35   0   0   0
H   -3.880000114440918   -1.0959999561309814   2.4159998893737793   18   35   0   0   0
C   -0.6460000276565552   0.49799999594688416   -1.1950000524520874   11   23   30   37   0
C   -1.9429999589920044   0.628000020980835   -0.7940000295639038   11   22   24   25   0
H   -2.5409998893737793   1.4290000200271606   -1.1950000524520874   12   23   0   0   0
C   -2.502000093460083   -0.2919999957084656   0.13199999928474426   78   2   23   26   0
C   -1.6759999990463257   -1.3370000123977661   0.621999979019165   11   25   27   28   0
H   -2.062999963760376   -2.049999952316284   1.3309999704360962   12   26   0   0   0
C   -0.3790000081062317   -1.4429999589920044   0.2070000022649765   11   26   29   30   0
H   0.2370000034570694   -2.24399995803833   0.5839999914169312   12   28   0   0   0
C   0.16599999368190765   -0.5350000262260437   -0.7070000171661377   11   3   22   28   0
C   -4.611000061035156   0.9269999861717224   0.07400000095367432   80   2   32   33   34
H   -4.130000114440918   1.878000020980835   0.29100000858306885   18   31   0   0   0
H   -5.561999797821045   0.8830000162124634   0.5889999866485596   18   31   0   0   0
H   -4.78000020980835   0.847000002861023   -0.9990000128746033   18   31   0   0   0
C   -4.369999885559082   -1.1449999809265137   1.444000005722046   80   2   19   20   21
H   1.246999979019165   -1.149999976158142   -3.177000045776367   6   16   0   0   0
H   -0.24899999797344208   1.2120000123977661   -1.899999976158142   12   22   0   0   0
C   4.420928001403809   2.3830971717834473   4.794969081878662   -1   39   40   41   42
H   5.3337273597717285   2.0750107765197754   4.2416157722473145   -4   38   0   0   0
H   4.0875630378723145   3.2775940895080566   4.226583480834961   -4   38   0   0   0
H   3.666041374206543   1.5751404762268066   4.686923980712891   -4   38   0   0   0
C   4.7958550453186035   2.593752384185791   6.205134868621826   -2   38   43   0   0
N   5.115570545196533   2.7171120643615723   7.319411277770996   -3   42   0   0   0
C   -5.660177707672119   0.0014052585465833545   -3.7105185985565186   -1   45   46   47  48
H   -6.3960862159729   -0.8015046119689941   -3.9298226833343506   -4   44   0   0   0
H   -5.950205326080322   0.23080207407474518   -2.662862777709961   -4   44   0   0   0
H   -4.610422611236572   -0.3616585433483124   -3.733290672302246   -4   44   0   0   0
C   -5.9191460609436035   1.085862636566162   -4.668871879577637   -2   44   49   0   0
N   -6.045229434967041   2.0538878440856934   -5.320666790008545   -3   48   0   0   0
$end

$rem
jobtype             sp
method              cam-b3lyp
dft_d               d3_bj
basis               mixed
QM_MM_INTERFACE     JANUS
FORCE_FIELD         OPLSAA
!FORCEMAN_PRINT      4
USER_CONNECT        TRUE
PRINT_GENERAL_BASIS true
SCF_CONVERGENCE     8
MAX_SCF_CYCLES      300
scf_final_print     true
MEM_STATIC          8500
MEM_TOTAL           85000
XC_GRID             3
SOLVENT_METHOD      PCM
$end

$qm_atoms
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
$end

$force_field_params
NumAtomTypes 4
AtomType -1   -0.080    3.300     0.066
AtomType -2    0.460    3.300     0.066
AtomType -3   -0.560    3.200     0.170
AtomType -4    0.060    2.500     0.015
Bond  -1  -2    390.0       1.470
Bond  -2  -3    650.0       1.157
Bond  -1  -4    340.0       1.090
Angle  -4  -1  -2    35.00     108.50
Angle  -1  -2  -3   150.00     180.00
Angle  -4  -1  -4    33.00     107.80
$end

$pcm
Theory              CPCM
Solver              CG
Method              SWIG
!CGthresh            5
!SwitchThresh        6
Radii               FF
!SASradius           1.4
!vdwScale            1.0
NoMatrix
Precond
UseMultipole
$end

$solvent
Dielectric          36.6
OpticalDielectric   1.7956
$end

$basis
O 1
def2-TZVPP
***
N 2
def2-TZVPP
***
C 3
def2-TZVPP
***
O 4
def2-TZVPP
***
C 5
def2-TZVPP
***
C 6
def2-TZVPP
***
C 7
def2-TZVPP
***
H 8
def2-TZVPP
***
C 9
def2-TZVPP
***
H 10
def2-TZVPP
***
C 11
def2-TZVPP
***
H 12
def2-TZVPP
***
C 13
def2-TZVPP
***
H 14
def2-TZVPP
***
C 15
def2-TZVPP
***
C 16
def2-TZVPP
***
H 17
def2-TZVPP
***
H 18
def2-TZVPP
***
H 19
def2-TZVPP
***
H 20
def2-TZVPP
***
H 21
def2-TZVPP
***
C 22
def2-TZVPP
***
C 23
def2-TZVPP
***
H 24
def2-TZVPP
***
C 25
def2-TZVPP
***
C 26
def2-TZVPP
***
H 27
def2-TZVPP
***
C 28
def2-TZVPP
***
H 29
def2-TZVPP
***
C 30
def2-TZVPP
***
C 31
def2-TZVPP
***
H 32
def2-TZVPP
***
H 33
def2-TZVPP
***
H 34
def2-TZVPP
***
C 35
def2-TZVPP
***
H 36
def2-TZVPP
***
H 37
def2-TZVPP
***
C 38
def2-SVP
***
H 39
def2-SVP
***
H 40
def2-SVP
***
H 41
def2-SVP
***
C 42
def2-SVP
***
N 43
def2-SVP
***
$end

@@@

$molecule
    read
$end

$rem
jobtype             sp
method              cam-b3lyp
dft_d               d3_bj
basis               mixed
QM_MM_INTERFACE     JANUS
FORCE_FIELD         OPLSAA
USER_CONNECT        TRUE
PRINT_GENERAL_BASIS true
SCF_GUESS           read
SCF_CONVERGENCE     8
MAX_SCF_CYCLES      300
scf_final_print     true
MEM_STATIC          8500
MEM_TOTAL           85000
XC_GRID             3
IQMOL_FCHK          True
ROKS                True
UNRESTRICTED        False
SOLVENT_METHOD      PCM
$end

$qm_atoms
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
$end

$force_field_params
NumAtomTypes 4
AtomType -1   -0.080    3.300     0.066
AtomType -2    0.460    3.300     0.066
AtomType -3   -0.560    3.200     0.170
AtomType -4    0.060    2.500     0.015
Bond  -1  -2    390.0       1.470
Bond  -2  -3    650.0       1.157
Bond  -1  -4    340.0       1.090
Angle  -4  -1  -2    35.00     108.50
Angle  -1  -2  -3   150.00     180.00
Angle  -4  -1  -4    33.00     107.80
$end

$pcm
Theory              CPCM
Solver              CG
Method              SWIG
!CGthresh            5
!SwitchThresh        6
Radii               FF
!SASradius           1.4
!vdwScale            1.0
NoMatrix
Precond
UseMultipole
rf_ptss_save        true
$end

$solvent
Dielectric          36.6
OpticalDielectric   1.7956
$end

$basis
O 1
def2-TZVPP
***
N 2
def2-TZVPP
***
C 3
def2-TZVPP
***
O 4
def2-TZVPP
***
C 5
def2-TZVPP
***
C 6
def2-TZVPP
***
C 7
def2-TZVPP
***
H 8
def2-TZVPP
***
C 9
def2-TZVPP
***
H 10
def2-TZVPP
***
C 11
def2-TZVPP
***
H 12
def2-TZVPP
***
C 13
def2-TZVPP
***
H 14
def2-TZVPP
***
C 15
def2-TZVPP
***
C 16
def2-TZVPP
***
H 17
def2-TZVPP
***
H 18
def2-TZVPP
***
H 19
def2-TZVPP
***
H 20
def2-TZVPP
***
H 21
def2-TZVPP
***
C 22
def2-TZVPP
***
C 23
def2-TZVPP
***
H 24
def2-TZVPP
***
C 25
def2-TZVPP
***
C 26
def2-TZVPP
***
H 27
def2-TZVPP
***
C 28
def2-TZVPP
***
H 29
def2-TZVPP
***
C 30
def2-TZVPP
***
C 31
def2-TZVPP
***
H 32
def2-TZVPP
***
H 33
def2-TZVPP
***
H 34
def2-TZVPP
***
C 35
def2-TZVPP
***
H 36
def2-TZVPP
***
H 37
def2-TZVPP
***
C 38
def2-SVP
***
H 39
def2-SVP
***
H 40
def2-SVP
***
H 41
def2-SVP
***
C 42
def2-SVP
***
N 43
def2-SVP
***
$end

@@@

$molecule
    read
$end

$rem
jobtype             sp
method              cam-b3lyp
dft_d               d3_bj
basis               mixed
QM_MM_INTERFACE     JANUS
FORCE_FIELD         OPLSAA
USER_CONNECT        TRUE
PRINT_GENERAL_BASIS true
SCF_GUESS           read
SCF_CONVERGENCE     8
MAX_SCF_CYCLES      300
scf_final_print     true
MEM_STATIC          8500
MEM_TOTAL           85000
XC_GRID             3
SOLVENT_METHOD      PCM
$end

$qm_atoms
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
$end

$force_field_params
NumAtomTypes 4
AtomType -1   -0.080    3.300     0.066
AtomType -2    0.460    3.300     0.066
AtomType -3   -0.560    3.200     0.170
AtomType -4    0.060    2.500     0.015
Bond  -1  -2    390.0       1.470
Bond  -2  -3    650.0       1.157
Bond  -1  -4    340.0       1.090
Angle  -4  -1  -2    35.00     108.50
Angle  -1  -2  -3   150.00     180.00
Angle  -4  -1  -4    33.00     107.80
$end

$pcm
Theory              CPCM
Solver              CG
Method              SWIG
!CGthresh            5
!SwitchThresh        6
Radii               FF
!SASradius           1.4
!vdwScale            1.0
NoMatrix
Precond
UseMultipole
rf_ptss_read        true
$end

$solvent
Dielectric          36.6
OpticalDielectric   1.7956
$end

$basis
O 1
def2-TZVPP
***
N 2
def2-TZVPP
***
C 3
def2-TZVPP
***
O 4
def2-TZVPP
***
C 5
def2-TZVPP
***
C 6
def2-TZVPP
***
C 7
def2-TZVPP
***
H 8
def2-TZVPP
***
C 9
def2-TZVPP
***
H 10
def2-TZVPP
***
C 11
def2-TZVPP
***
H 12
def2-TZVPP
***
C 13
def2-TZVPP
***
H 14
def2-TZVPP
***
C 15
def2-TZVPP
***
C 16
def2-TZVPP
***
H 17
def2-TZVPP
***
H 18
def2-TZVPP
***
H 19
def2-TZVPP
***
H 20
def2-TZVPP
***
H 21
def2-TZVPP
***
C 22
def2-TZVPP
***
C 23
def2-TZVPP
***
H 24
def2-TZVPP
***
C 25
def2-TZVPP
***
C 26
def2-TZVPP
***
H 27
def2-TZVPP
***
C 28
def2-TZVPP
***
H 29
def2-TZVPP
***
C 30
def2-TZVPP
***
C 31
def2-TZVPP
***
H 32
def2-TZVPP
***
H 33
def2-TZVPP
***
H 34
def2-TZVPP
***
C 35
def2-TZVPP
***
H 36
def2-TZVPP
***
H 37
def2-TZVPP
***
C 38
def2-SVP
***
H 39
def2-SVP
***
H 40
def2-SVP
***
H 41
def2-SVP
***
C 42
def2-SVP
***
N 43
def2-SVP
***
$end

Sincerely,
Tiopirfur

It means you’re having to diagonalize the solvent response matrix rather than using an iterative solver. That matrix is probably enormous (e.g., 423,445 surface tesserae in my calculation), so that will be incredibly slow and I’m surprised that it even works.

This job has been running on my system for 90+ hours and is in the SCF iterations for the 3rd job, without crashing. (This is Q-Chem v. 6.3.)

I can reproduce your crash, finally after 4 days on 48 threads. The fuller context of the error message that I get is the following:

### Frozen reaction-field calculation using stored ESP to compute ASC ###
Reading ESP for state 0 (0 is MP-GS) from file...
... done. Computing ASC and self-interaction...
... done. Updating ref-RF self-interaction (-11.1539 eV) ...
... done. Separating ASC into fast and slow part...

 Q-Chem warning in module libmdc/NewQAlloc.C, line 469:

 QAllocGeneral() returns a null pointer when allocating 1434445344200 Bytes memory


 Q-Chem fatal error occurred in module include/RefCount.h, line 48:

 RefCount::Constructor called with null pointer argument

Note how the calculation is trying to allocate 1 Tb of memory; that is why it crashes. I think this job is just too large. You could try some of the accelerators for PCM (to not store matrices, in addition to CG solver, and maybe multipole approximations), but this may simply be too large of a cavity surface.

Hello,
Big thank you for your investigation. I have to be more carefull with jobs sizes in the future.
Sincerely,
Tiopirfur

You can also reduce HPoints and HeavyPoints in $pcm section (reducing the discretization grid), but this is so enormous I’m not sure that will help. The options that don’t store matrices are your best bet although those can be slower because they are iterative.