SS(V)PE input question

Hello, I have a question about how to set up calculations to use use SS(V)PE electrostatics that I suspect @jherbert can answer. Specifically, I’d like to do geometry optimizations with this solvent model.

According to the manual I can use SS(V)PE by setting

$rem
   solvent_method = isosvp
$end

and this appears to work as intended. I’m wondering whether there is any difference between the above and

$rem
   solvent_method = pcm
$end

$pcm
   theory=ssvpe
$end

?

If there is a subtle difference between the two, is one more appropriate for a geometry optimization than another?

Finally, I am interested in seeing the dimensions of the isodensity cavity that the model constructs. Is there a way to output this or to read the data somewhere in the output file?

Thank you!

This is admittedly confusing although it is discussed in the (very long) chapter on solvation models: 11.2 Chemical Solvent Models‣ Chapter 11 Molecules in Complex Environments: Solvent Models, QM/MM and QM/EFP Features, Density Embedding ‣ Q-Chem 5.4 User’s Manual.

The differences between these two approaches comes down to the definition of the solute cavity. SOLVENT_METHOD = ISOSVP implements the SS(V)PE model developed by Dan Chipman, using a cavity that is defined by an isocontour of the solute’s electron density (as done in Chipman’s papers). That’s a nice, clean, one-parameter cavity definition that I like a lot, but with the important limitation that it lacks analytic gradients. You can consult my review (https://doi.org/10.1002/wcms.1519) for a discussion of why those gradients are challenging to formulate. However, it sounds like this is a non-starter for what you want to do.

Methods that do have gradients include SOLVENT_METHOD = PCM, = SM8, and = SMD. The two SMx models are packaged up with nonelectrostatics terms in addition to electrostatics, and in the SMD case the electrostatics actually come from PCM. (There are also nonelectrostatic interactions for the isodensity SS(V)PE; that model is called CMIRS: https://doi.org/10.1021/acs.jctc.6b00644, but it inherits the lack of gradients.)

Therefore, for geometry optimizations let’s talk about PCM. These methods use union-of-atomic-spheres cavity constructions, with atomic radii that are parameters of the model. Gradients are available. There are several different PCMs, the choice of which is stipulated using the “Theory” keyword in the $pcm input section:

$pcm
Theory ssvpe
$end

(Other options are iefpcm and cpcm, and note that there is no “=” as in your example.) The differences between these approaches are somewhat subtle; see the review cited above or else here: Redirecting). SS(V)PE and IEF-PCM differ only by symmetrization of certain matrices that makes SS(V)PE about 2x faster (probably not noticeable except for very large molecules as the PCM time is minor overhead), but this speedup sacrifices some formal properties when van der Waals radii are used to construct the cavity. IEF-PCM uses asymmetric matrices but preserves those properties and is what I recommend. However, C-PCM is equivalent to IEF-PCM in the high-dielectric limit and is simpler and faster; “high-dielectric limit” operationally means epsilon >~ 30-35, i.e., acetonitrile or anything more polar than that.

Summary: for geometry optimizations, I recommend SOLVENT_METHOD = PCM with Theory set to either IEF-PCM for nonpolar solvents or C-PCM for polar solvents.

For “dimensions of isodensity cavity”, I assume (?) that you mean just printing out the discretization points that define it. There’s probably a print option to do that; I will poke around the manual and see if I can find it, but if you’re in a hurry check out this section of the manual: 11.2.5 Isodensity Implementation of SS(V)PE‣ 11.2 Chemical Solvent Models ‣ Chapter 11 Molecules in Complex Environments: Solvent Models, QM/MM and QM/EFP Features, Density Embedding ‣ Q-Chem 5.4 User’s Manual. This is where isodensity implementation is described.

Thank you, John! I guess I was conflating SS(V)PE in general with the isodensity implementation of it. I agree with you that a find the physical consistency of using an isodensity cavity very appealing. If I’m understanding correctly, the only way to get that is with solvent_method=isosvp, whereas setting theory ssvpe will still use empirical radii, it just formulates the underlying electrostatics a little differently.

As an aside, I know from your review and webinar that the new CMIRS model uses SS(V)PE electrostatics. Is that the isodensity variety, or are the radii parameterized?

(1) Yes, only SOLVENT_METHOD = ISOSVP uses an isodensity cavity, but lacks gradients.
(2) SOLVENT_METHOD = PCM uses (scaled) empirical radii for all values of Theory in the $pcm section.
(3) Setting Theory to ssvpe in the $pcm section gives you the same underlying equations as the ISOSVP method but different cavity surface (empirical radii vs. isodensity)
(4) CMIRS uses isodensity cavity, i.e., SOLVENT_METHOD=ISOSVP. It is possibly worth considering whether one could “get away with” empirical radii, which has not been examined in literature to the best of my knowledge. There is some reason to believe that the separability of electrostatics and nonelectrostatics that CMIRS is built upon, leading to a much smaller number of parameters as compared to SMx models, might not survive such a switch. However, this is speculation on my part.

Regarding your question about cavity surface dimension: cavity surface area and volume are printed out, e.g., in the sample job SVP_H2O.in:

Properties
--------------------
Cavity Surface Area =                 149.7698546355  (  41.93983 ANG**2)
Cavity Volume =                       171.0434878910  (  25.34603 ANG**3)
Charge Penetration From Exp Fit =      -0.0592567381
Total Surface Charge =                 -0.0046785901
Charge Rule = -Q*EPS/(EPS-1) =          0.0047390448
Solute Charge - Charge Rule =          -0.0000604547

If you are looking to print the discretization charges on the grid, subroutines to do that are present in the code but I was unable to find a user-facing set of input options that would activate them. If this is what you want, please contact Q-Chem support. It’s a trivial thing to add but unless I’m missing something, it would need to be pushed out in a future release.

Thanks so much @jherbert ; your response covers everything I need right now. Thanks very much for clarifying.