PCM Solver Geometry Optimization

[theleftcoast]

Hello,

I have been using PSI4 to investigate non-covalent interactions for associating systems (hydrogen bonding and mixed influence). All of the calculations I have done so far have been in vacuum, but most of the applications I am interested in are in solution. To get an understanding of behavior in solution, I have started to explore PCMSolver. The first trial system was the water-water complex and my work flow was as follows.

1. Optimize geometry in vacuum (MP2 with aug-cc-pVTZ).
2. Calculate thermodynamic properties in vacuum.
3. Re-optimize geometry in solution (MP2/PCM with aug-cc-pVTZ).
4. Calculate thermodynamic properties in solution.

The hardest step seems to be re-optimizing the geometry with PCM. I used the finite difference approach (dertype=‘energy’), and it took forever for the optimization to converge (59 iterations, even when starting with an optimized geometry in vacuum). Even though it took a while, the resulting geometry seems reasonable. The calculated hydrogen bond length shortened from 1.929 angstroms to 1.834 angstroms. The relative orientation between the water molecules also changed a bit (see the xyz coordinates shown below). Do these changes in geometry seem like a reasonable result for a hydrogen bonding complex? Should I expect bond lengths to shorten a bit for most complexes when re-optimizing the geometry in water?

In Vacuum (MP2 with aug-cc-pVTZ):

O           -1.478528905013    -0.073068979734     0.004135826010
H           -1.843856649950     0.812352837817    -0.016261037228
H           -0.521931038934     0.057289764657    -0.000318208491
O            1.406567514031     0.069106936961    -0.003826657324
H            1.755671715105    -0.386555474341     0.765611279704
H            1.752195469767    -0.420206634134    -0.753938765385

In Solution (MP2/PCM with aug-cc-pVTZ, solvent was water):

O           -1.442781654230     0.039594940486     0.059699655839
H           -1.685244466118     0.942498597239    -0.160539782306
H           -0.472335854044     0.031272046541     0.016166508974
O            1.358465075784    -0.042253313013    -0.059803175369
H            1.787088761920    -0.206022777191     0.785355001878
H            1.708656845205    -0.725557565457    -0.639338798480

I am also interested to know if the input parameters I used configure PCMSolver are reasonable for investigating hydrogen bonded geometries in solution. Most parameters have been left at the default values suggested in the PCMSolver documentation. However, I did change ‘RadiiSet = UFF’ and left ‘Scaling = False’. What is the impact of setting ‘Scaling = True’? Sorry if these are basic questions, but I am just beginning to learn about the Polarizable Continuum Model. Any references you can point me to are appreciated.

set {
basis aug-cc-pVTZ
opt_coordinates cartesian
scf_type pk
memory 16000mb
pcm true
pcm_scf_type total
}

pcm = {
Units = Angstrom
Medium {
SolverType = IEFPCM
Solvent = Water
}

Cavity {
RadiiSet = UFF
Type = GePol
Scaling = False
Area = 0.3
Mode = Implicit
}
}

optimize(‘MP2’, molecule=water_water, dertype=‘energy’)

I know this is a long post, but I’m really trying to understand the basics before moving ahead to the other associating complexes in my investigation (30 systems in total). Thanks in advance for any assistance the community can provide. It is very much appreciated!

John

[robertodr]

Hi!
Thanks for your interest in Psi4 and PCMSolver (and apologies for the late response!)

1. The standard reference for details about the polarizable continuum model is the 2005 Chemical Reviews paper by Tomasi, Mennucci and Cammi. Another good resource is the Continuum solvation models in chemical physics book
2. As I’m sure you’re aware, PCMSolver uses Bondi radii scaled by 1.2 by default. The Bondi radii of hydrogen and oxygen are 1.2 Ang and 1.52 Ang, respectively. Scaling by 1.2 results in 1.44 Ang and 1.824 Ang, respectively. UFF radii without scaling will give 1.4430 Ang and 1.75 Ang, respectively. So you will obtain a roughly similar cavity. The cavity can be visualized from the cavity.off (using GeomView) or cavity.npz (using the plot_cavity.py script in the PCMSolver repository) files. That famous program defaults to UFF scaled by 1.1 (and also uses slightly different values for the solvent permittivities) If you want to use the same set of radii with PCMSolver, you will have to specify the list of spheres explicitly.[^1] To conclude, I think you’re using perfectly reasonable values for your cavity
3. MP2 with PCM in Psi4 is implemented in the perturbation-to-the-energy (PTE) formalism. The mutual solute-solvent polarization is self-consistent at the SCF level, but remains frozen at the correlated level. The self-consistent scheme at the correlated level would require, as shown in the paper I referenced, to compute the MP2 density and iterate the calculation of the MP2 energy correction. Such a scheme has been in the works and almost ready for CCSD for forever (my bad…) but I doubt it will be necessary for what you are interested in.
4. Geometries optimizations with analytical gradients in Psi4 do not include the PCM contributions. But you’re correctly using the numerical gradient. This should get you a geometry optimization in solution, BUT
   a. I have not tested this explicitly myself;
   b. I am not familiar with the numerical gradient geometry optimization machinery in Psi4 to confidently say that everything is indeed in the right place. @loriab certainly knows more than me in this respect.
5. The geometry optimization in vacuo was run using the analytical gradient, I assume.
   a. How many iterations does that take?
   b. How many iterations does the same optimization run with the numerical gradient take?
   c. How much does a single-point calculation take in vacuo?
   d. How much does a single-point with PCM take?
   e. Could you post the final cavity.off and/or cavity.npz? It could be that at some point during the optimization two disjoint molecular cavities are produced and that would definitely make converge harder with PCM.
   f. The Psi4 manual section on PCM shows only examples using the PK algorithm, but DF and CD also work. This has been tested and the tests are in a pull request in the pipeline on GitHub. It could speed up your calculations.
   g. It might also be helpful in speeding up the geometry optimizations to first optimize with cc-pVTZ and then restart from that geometry using aug-cc-pVTZ.
6. I cannot offer any advice on your very first question. The resulting geometry seems reasonable to me too.

Hope this helps!

[^1] Note, however, that if you’re after perfect agreement between PCM results between different implementations you will have an extremely hard time. Each program/library uses not only different defaults (which one can fiddle with rather easily!), but rather different implementations of the cavity generation, Green’s functions, integral operators matrices and linear equation solver.

[theleftcoast]

Hello,

That was extremely helpful!

1. I’ll take a look at the 2005 Chemical Reviews paper. This looks like a great resource, so thank you for the recommendation.
2. Since ‘RadiiSet = UFF’ and ‘Scaling = False’ will provide a reasonable cavity, I’ll continue to use these settings for future optimizations.
3. Thank you for the reference. Increasing my awareness of how PCMSolver is implemented is helpful while I’m learning. I am sure the current implementation is more than sufficient for my project.
4. Glad that that numerical gradients should work for geometry optimization.
5. Correct, the water_water geometry optimization in vacuo was run with analytical gradients. I’ll re-run the water_water geometry optimization with numeric gradients and compare. My plan is to follow up on points ‘a’-‘e’ in a future reply. I am also running a nitrogen_water geometry optimization with PCM-MP2/aug-cc-pVTZ. Once that finishes, I can post those results too. The point you bring up in ‘e’ seems like a good check to include after any geometry optimization with PCMSolver enabled. I found the plot_cavity.py script and will review the cavities for both water_water and nitrogen_water.
6. Sounds good. I’ll get a better feel for how the geometry changes as I optimize more complexes.

John

[kalju]

The shortening you observed makes sense to me because the dipole-dipole interaction contributes significantly to hydrogen bond. In a very simplified picture, the medium affects the hydrogen bonding through two factors:

1. The medium polarizes the electron density of water molecules. The contribution of this effect is relatively easy to estimate. The dipole moment of H2O in water is about 2.9D while a singe unpolarized water molecule in the gas phase has a dipole moment of 1.85 D. One can calculate the approximate dipole-dipole interaction energy using a classical expression (given, for example in “Principles of Biomolecular Recognition” ) and estimate the contribution due to the polarizing effect of the medium.
2. On the the other hand, the medium also attenuates any Coulombic interactions, including the dipole-dipole interaction. In the case of long range interactions (such as ions in dilute aqueous solution), this is well accounted for by the dielectric constant (about 80 for liquid water). In the case of intimately bound pairs, like you have, there is no easy classical way to say how much the interaction is weakened.

Thus, the shortening that you observed may well reflect the increase of the dipole moment of two water monomers in the presence of the polarizing medium. I do not think it is a general phenomenon, though. You could have molecules where the polarization of monomer charge densities is not so prominent as in H2O, or the dielectric attenuation is bigger because the two monomers are further away. With heterodimeric complexes, the dispersion contribution can play interesting tricks: a repulsive London dispersion is possible if one solute has a larger refractive index than the solvent and the other has smaller refractive index than the solvent. I am not sure if PCMSolver accounts for the latter effect.

[theleftcoast]

Hello,

Quick update…water-water, water-N2, water-CO2, water-CO, water-H2, and water-H2S all successfully converge using the following “recipe.” The final optimized geometry for each complex was contained inside a single molecular cavity.

set {
basis aug-cc-pVTZ
opt_coordinates both
dynamic_level 3
scf_type pk
memory 16000mb
geom_maxiter 200
pcm true
pcm_scf_type total
}

pcm = {
Units = Angstrom
Medium {
SolverType = IEFPCM
Solvent = Water
}
Cavity {
RadiiSet = UFF
Type = GePol
Scaling = False
Area = 0.3
Mode = Implicit
}
}

optimize(‘MP2’, molecule=dimer, dertype=‘energy’)

It turns out that weakly bound complexes (like water-N2) would often never converge without setting dynamic_level > 0. For instance, water-N2 moved towards a reasonable geometry in the first ~30 iterations, but would make no real progress towards a lower energy state for the next 150+ iterations. After trial and error, it looks like dynamic_level 3 works best with convergence typically obtained in 5-15 iterations.

Optimization have been very slow for larger complexes (methanol-benzene, water-dimethyl ether, etc), but that is understandable when using numerical gradient since the number of displacements grows rapidly with the number of atoms.

Bond lengths for water-N2, water-CO2, water-CO, water-H2, and water-H2S did not shorten as much as for water-water.

Thank you again to @robertodr and @kalju for your help!

John

[amuralid]

“” defaults to UFF scaled by 1.1 (and also uses slightly different values for the solvent permittivities) If you want to use the same set of radii with PCMSolver, you will have to specify the list of spheres explicitly""

Hi @robertodr: Thanks for the detailed comment!
Quich question: Is there simple way to set scaling factor to 1.1*UFF?
Atleast for a small list of atoms (like C-S-N-O-H) The explicit option seems to not be user friendly… even sphere centers have to be set…

[robertodr]

The input style is quite deliberately explicit and verbose, because I could not make many assumptions on the format of the data the library would be passed. You can work in so-called “atoms” mode which is in between “implicit” and “explicit” in user-friendliness. Say you want to change the radius on atom 10 in your geometry to 1.8, while keeping UFF without scaling on all the others:

cavity
{
	type = gepol
	area = 0.6
	radiiset = uff
    mode = atoms
    scaling = false
    atoms = [10]
    radii = [1.8]
}

See also here: Input description — PCMSolver documentation