Webinar 51 by John Herbert: Implicit Solvation Models in Q-Chem, for Ground and Excited States

Webinars
#1

Recording: Webinar 51: Implicit Solvation Models in Q-Chem, for Ground and Excited States - YouTube

Abstract:

This talk will provide an overview of dielectric continuum implicit solvent models available in Q-Chem. For ground-state properties, these include several SMx models along with PCMs. Nonequilibrium versions of PCMs are available that afford corrections for vertical excitation and vertical ionization energies. Finally, a Poisson Equation Solver (PEqS) is available that can provide arbitrary, anisotropic dielectric boundary conditions for electronic structure calculations.

Please post your questions under this topic.

#2

FYI, this webinar was structured around my recent review article that you can find here: https://doi.org/10.1002/wcms.1519. Feel free to email me for a copy if you don’t have access to WIRES-CMS.

#3

@jherbert , thank you for the excellent webinar and review on implicit solvation models. I am very interested in trying out the CMIRS model you described, but am having trouble getting calculations to converge. I am specifically trying to use it to calculate solvation energies of single monatomic ions, so the very small system size might be to blame.

The gas-phase part of the SVP calculation succeeds, but then I get DiagJacobi failed to converge errors when starting the SVP portion. This occurs for small and large basis (def2-SVPD, TZPVD, or QZVPD), for all five solvents available in CMIRS, and for multiple ion types. Can you provide any guidance on how to troubleshoot?

I have already tried lowering the various convergence parameters like CVGLIN, SVP_CHARGE_CONV, and SVP_CAVITY_CONV, and have also tried reducing the no. of Lebdev points via NPTLEB, all with the same result.

Here is the tail of the output.

The Gas Phase SCF has converged
 ---------------------------------------


 Now restarting with SVP Calculation
          --------------
          -SVPE- SOLVENT
          --------------

 ----- NAMELIST $SVP -----
 ISHAPE=   0 ITRCAV=  99 NDRCAV=   0
 INTCAV=   0 NPTLEB=  50 NPTTHE=   8 NPTPHI=  16
 ITRNGR=   2 IROTGR=   2
 LINEQ =   1 IRDRF =   0 IPNRF =   1
 IOPPRD=   0
 IDEFESR=   1

 RHOISO=    0.001000 TOLCAV=  0.1000D-09 RADSPH=    2.645617 TOLCHG=  0.1000D-06
 DIELST=   78.540000 CVGLIN=  0.1000D-06 CSDIAG=    1.104000
 ROTTHE=    0.000000 ROTPHI=    0.000000 ROTCHI=    0.000000
 TRANX =    0.000000 TRANY =    0.000000 TRANZ =    0.000000


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

 DiagJacobi failed to converge in 100 sweeps

And the input:

$molecule
 1 1
 Tl      0.0000000000      0.0000000000      0.0000000000
$end

$rem
   job_type = sp
   basis = def2-qzvpd
   max_scf_cycles = 200
   xc_grid = 3
   resp_charges = false
   symmetry = false
   sym_ignore = true
   method = wb97xv
   plots = true
   make_cube_files = true
   ecp = def2-ecp
   solvent_method = isosvp
   thresh = 14
$end

$plots
   grid_spacing 0.05
   total_density 0
$end

$svp
rhoiso=0.001, dielst=78.54, ipnrf=1, idefesr=1
$end

$pcm_nonels
   a -0.006736
   b 0.032698
   c -1249.6
   d -21.405
   gamma 3.7
   solvrho 0.05
   delta 7.0
   gaulag_n 40
$end

Thank you in advance for any guidance you can offer.

#4

Hi Ryan,
I am able to reproduce this crash with the latest version of Q-Chem (5.4.1). It is coming from the isodensity implementation of the electrostatics [SS(V)PE] model, independent of the nonelectrostatic parts of CMIRS. There are a couple of things going on, which are NOT obvious (took me a while).

(1) The diagonalization error that you see can be made to go away through a combination of IRotGr=0 in the $svp section and increasing the number of Lebedev points used to discretization the isodensity surface (e.g., NPtLeb=974). I need to do more testing to figure out just what it is about this job that creates problems, but when I come up with some general advice I’ll post it here and add it to the manual.

(2) Based on looking at the energetics, I strongly suspect that the implementation of SS(V)PE does not work in the presence of pseudopotentials.

That 2nd part may be a problem for your application, but if it’s not, then I was able to get your job to run by commenting out the ECP=DEF2-ECP line in $rem, and using the following for $svp:

$svp
rhoiso=0.001, dielst=78.54, ipnrf=1, idefesr=1, irotgr=0, itrngr=0, nptleb=1454
$end

I need to check the code to be sure about the pseudopotential part (which may be a relatively easy fix), but I’m pretty sure because if you make the changes to $svp but leave the pseudopotential turned on, you will notice that the SCF energy of -3998 Ha suddenly gets MUCH smaller when solvation model turns out.

Let me know if you have additional problems. For now, the solution is no ECPs and change the $svp. Try to fix the former and get better guidance about the latter into the manual, for a future release.

#5

@jherbert , thank you so much for the thorough reply. Another thing I noticed in my testing is that you must set gen_scfman to false in order for the CMIRS (or possibly just SS(V)PE) calc to run. If you use the example CMIRS input calc from the manual but add gen_scfman = true, I think you’ll find that it fails.

I will do more testing based on your suggestions and report back.

#6

You are correct about GEN_SCFMAN although that is now set to false automatically for SOLVENT_METHOD = ISOSVP. That may not be automatic some in older versions of Q-Chem

#7

For what it’s worth, I think I can confirm your hypothesis about the ECPs. I tried a few different elements and the calc only fails when the ECP is active. The def2-ECP is only defined for elements Rb and higher. So, I reran the above calculation using the SVP settings below, but kept the ECP in the input file. I tried Li, K, Rb, and Tl.

$svp
rhoiso=0.001, irotgr=0, dielst=78.54, nptleb=1454, ipnrf=1, idefesr=1
$end

With Li and K the calculation worked fine. With Rb and Tl (where the ECP is defined) it failed with the following message

 -------------------------------------------------------
  Cycle       Energy         DIIS Error    SVP Error
 -------------------------------------------------------
    1    -635.7961540595      1.61E-02      4.58E-02
    2   -1492.0846809472      5.24E-02      1.06E-02
 SINGLE CENTER EXPANSION VALIDITY CHECK INDICATES FAILURE AT  1454 POINTS. STOP.
 LARGEST FAILURE WAS FOUND FOR SHELL    1 ON RAY    7:
 INNER ISODENSITY POINT IS MORE THAN    0.5917795 AU FROM OUTER ISODENSITY POINT
 APPROX INNER ISODENSITY POINT AT    0.2182068    0.2182068    0.2182068 AU HAS RHO=  0.69459D-03
 ACTUAL OUTER ISODENSITY POINT AT    0.5598708    0.5598708    0.5598708 AU HAS RHO=  0.10000D-02

Interestingly, I was able to complete a calculation with I-, for which the ECP should be active. However as you noted above, the energies during the gas phase calculation and the SVP portion were vastly differently, leading to an absurdly large solvation energy.

Some excerpts from the output are below.

 ---------------------------------------
  Cycle       Energy         DIIS Error
 ---------------------------------------
    1    -297.5436466021      3.46E-03
    2    -297.7515104287      4.11E-03
    3    -297.7775438125      1.77E-03
    4    -297.7861477468      6.07E-04
    5    -297.7873579672      1.95E-04
    6    -297.7874147021      1.95E-05
    7    -297.7874173095      2.73E-06

The Gas Phase SCF has converged
 ---------------------------------------

 -------------------------------------------------------
  Cycle       Energy         DIIS Error    SVP Error
 -------------------------------------------------------
    1    -370.8343519107      1.13E-02      3.00E-02
    2    -497.9388207394      1.68E-02      6.28E-03
    3    -587.1410385074      8.96E-02      8.02E-03
    4    -586.4723711603      9.07E-02      1.02E-05
    5    -593.2094999315      8.51E-02      7.74E-05
    6    -611.7404577073      7.17E-02      2.02E-04
    7    -648.0011424077      4.50E-02      3.58E-04
    8    -778.8755624394      3.17E-02      5.68E-03
    9    -778.3690577913      3.21E-02      3.59E-06
   10    -781.3331891375      3.03E-02      2.15E-05
   11    -801.4124459684      1.72E-02      1.46E-04
   12    -827.7534450373      2.76E-03      1.72E-04
   13    -828.0942178472      2.19E-04      7.73E-06
   14    -828.1878216727      3.18E-05      3.90E-06
   15    -828.2115181717      2.94E-06      1.64E-06
   16    -828.2136833467      4.36E-07      1.28E-06
   17    -828.2136647450      3.98E-07      1.23E-06
   18    -828.2136658201      7.14E-08      1.34E-05
   19    -828.2136666446      3.82E-08      1.78E-07
   20    -828.2137586459      7.56E-07      9.99E-06
   21    -828.2136691835      8.12E-07      1.02E-06
   22    -828.2136682191      5.95E-07      6.54E-06
   23    -828.2136489638      1.21E-06      4.15E-06
   24    -828.2136757219      4.58E-07      3.96E-06
   25    -828.2136724040      1.60E-06      5.70E-06
   26    -828.2136651537      1.53E-06      1.34E-06
   27    -828.2136567014      1.49E-06      3.67E-07
   28    -828.2136665876      1.74E-08      8.17E-06
   29    -828.2136666643      1.34E-08      1.45E-07
   30    -828.2141870055      5.17E-08      7.50E-07
   31    -828.2136665429      1.20E-08      7.76E-07
   32    -828.2136665524      1.04E-08      2.48E-08 Convergence criterion met
 -------------------------------------------------------


--------------------------------------------------------------------------------
DEFESR calculation with single-center isodensity surface

 Total number of surface points = 974
 Total number of solute electrons = 26.000000
 Dispersion Time: 21.24 s (CPU) 13.64 s (wall)
 Exchange Time: 0.15 s (CPU) 0.15 s (wall)
 Extremum Field Time: 16.73 s (CPU) 0.75 s (wall)


--------------------------------------------------------------------------------
The Final SS(V)PE energies and Properties
--------------------------------------------------------------------------------

Energies
--------------------
The Final Solution-Phase Energy =    -828.2136665524
The Solute Internal Energy =         -102.2936344128
The Change in Solute Internal Energy = 195.4937828967  (122674.20029 KCAL/MOL)
The Reaction Field Free Energy =     -725.9200321396  (-455521.69536 KCAL/MOL)
The Dispersion Energy =                21.3976937694  (-127.69320 KCAL/MOL)
The Exchange Energy =                  13.8876983741  ( 390.33844 KCAL/MOL)
Min. Negative Field Energy =           25.2967250168  (   0.00000 KCAL/MOL)
Max. Positive Field Energy =           25.2968284992  (-57529928009.17393 KCAL/MOL)
The Total Solvation Free Energy = -91680306.4379944056  (-57530260594.02376 KCAL/MOL)
#8

Yes, thanks. For the next release I’ll fix this if I find the time or at least disable the job and print a sensible message about why.