Non-Covalent Bond Energy In Solution + Potential PCMSolver Bug

1

Hello!

Many systems with strong non-covalent bonds (water-water, methanol-methanol, water-CO2, etc) can be easily analyzed as a dimer (n=2) ‘in vacuum’. In solution, the density is much higher so associating molecules usually exist in larger clusters (n>2). This means that cooperative effects can significantly change the binding energy. For example, in Ab Initio and DFT Studies on Methanol-Water Clusters the authors show significant increases of binding energy per hydrogen bond for small clusters of methanol and water. I started with the optimized geometries presented in their paper for n=2 and n=4 and calculated binding energies per bond. The result is similar to what is shown in the paper, but there are some differences because I used different model chemistries.

wB97M-V/def2-QZVPPD (no CP)

  • H2O (n=2), -4.91 kcal/mol per bond
  • H2O (n=4), -7.04 kcal/mol per bond
  • Cooperative effects increased binding energy by 43%

MP2/aug-cc-pVTZ (CP)

  • H2O (n=2), -4.70kcal/mol per bond
  • H2O (n=4), -6.73 kcal/mol per bond
  • Cooperative effects increased binding energy by 43%

Clustering and the resulting cooperative effects can clearly change bind energies. The reference above shows that for methanol-methanol, water-water, and water-methanol, the average bond energy increases significantly.

Another important element of chemistry in solution is that the media surrounding associating clusters is polarizeable. I used the same optimized geometries as above and computed binding energies with PCMSolver enabled.

PCM+wB97M-V/def2-QZVPPD (no CP)

  • H2O (n=2), -3.43 kcal/mol per bond
  • H2O (n=4), -3.975 kcal/mol per bond
  • Cooperative effects increased binding energy by 16%

It seems like cooperative effects increase average binding energy for clusters while polarizable media effects reduce average binding energy for clusters ‘in solution’ relative to binding energies ‘in vacuum’.

I have a couple of questions:

  • Is this competing increase/decrease in binding energy from cooperative/polarization effects generally seen in other non-covalent interactions?
  • Has the impact of opposite acting effects of cooperative/polarization effects been discussed in the literature anywhere?

These two effects seem to cancel each other out somewhat. I can’t find any references so far which treat the interplay between these two effects. I’m sure cluster size is important as well. I am mostly looking for studies on small non-covalent complexes (thing S66 and A24 database type complexes) rather than on larger fragments in biochemical applications.

I also ran across a potential PCMSolver bug. I can optimize geometries and perform energy calculations just fine in most cases. Unfortunately, when I try to do energy calculations using counterpoise correction, Psi4 runs into errors. It looks very similar to the issue described in this post.

Thank you in advance for any suggestions or insight you might have for me!

John

2

yes, implicit solvation models will typically reduce interaction energies. There is a charge screening effect from the solvation model and will e.g. most notably reduce strong Coulomb electrostatic interactions, but generally all energy components are affected.
How much polarisation/induction is affected is not well known, I think, but this will have a direct influence on the degree of cooperativity.

One difficulty in studying the interplay you mentioned is to correctly account for SCRF artefacts like cavitation and outlying charges errors.

Could you upload an input example for the crashing PCM calculations?

1 Like
3

Hello,

Thank you for the reply. Here is an input file that works. DFT calculations with and without PCM on to compare results.

molecule water_water_vac {
0 1
H	1.822414	0.767503	-0.267574
O	1.375604	-0.00251	0.093561
H	1.83885	-0.746499	-0.299913
--
0 1
H	-0.570257	-0.005819	-0.014274
O	-1.530212	0.002849	-0.118763
H	-1.854144	-0.017895	0.783375
}

set memory 100000mb
set dft_spherical_points 590
set dft_radial_points 99

set basis def2-QZVPPD
wb97x_v = energy('wb97x-v', molecule=water_water_vac, bsse_type='NoCP')
print("wb97x-v noCP, (kcal/mol) = {}".format(wb97x_v*627.509))

set pcm true
set pcm_scf_type total

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

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

set basis def2-QZVPPD
wb97x_v = energy('wb97x-v', molecule=water_water_vac, bsse_type='NoCP')
print("wb97x-v noCP PCM, (kcal/mol) = {}".format(wb97x_v*627.509))

Here is an input file that causes a crash.

molecule water_water_vac {
0 1
H	1.822414	0.767503	-0.267574
O	1.375604	-0.00251	0.093561
H	1.83885	-0.746499	-0.299913
--
0 1
H	-0.570257	-0.005819	-0.014274
O	-1.530212	0.002849	-0.118763
H	-1.854144	-0.017895	0.783375
}

set memory 100000mb
set dft_spherical_points 590
set dft_radial_points 99
set pcm true
set pcm_scf_type total

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

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

set basis def2-QZVPPD
wb97x_v = energy('wb97x-v', molecule=water_water_vac, bsse_type='CP')
print("wb97x-v noCP PCM, (kcal/mol) = {}".format(wb97x_v*627.509))
4

I found an open issue about PCM and ghost-atoms (used in the CP procedure), where the latter do not get assigned a correct radius.

Please follow the link below for updates:

I’ll try to further the issue along.

5

The pcmsolver bug with ghost atoms was fixed in the master branch.

6

Hello,

Sorry for not replying sooner…it is good to hear that the bug has been addressed. Thank you very much for the help.

John