COSMO cavity construction options in NWChem

[rpseng]

Dear NWChem developers,

First, thank you very much for providing such a powerful package as open source and for being open to discussions and improvements.

For a long time, I have attempted to use NWChem for generating COSMO surface charge densities. However, since the results often did not align with previous findings, I eventually set it aside. Now, I believe it is a good time to revisit the topic and explore potential improvements.

The first point I would like to discuss concerns the total cavity surface area and volume. My findings are described below.

Cavity Area and Volume Comparison

For this comparison, let's use acetone, which runs fairly quickly in NWChem. The input file is given below.

To better highlight potential artifacts, I have selected rsol=1.4 for the NWChem runs.
For this molecule, the surface area and volume should only depend slightly on rsol.

Results for acetone using NWChem options (YK and KS) are shown below along with results from other software:

	YK	KS	GAMESS[^1]	DMol3[^2]
Area A2	104.59	66.575	102.0	102.65
Volume A3	58.7	36.836	85.92	86.42

• 1: lvpp sigma-profile database
• 2: VT sigma database

YK results

The calculated surface segments for YK are shown below:

By making the surface more translucent we can get a better view, since there are many segments actually inside the cavity:

As seen in the results, the YK cavity is essentially a van der Waals cavity.
As a result, it has a larger surface area and a smaller volume compared to the SES surface typically assumed in COSMO calculations.

KS results

The calculated surface segments for KS are shown below:

As seen in the results, surface segments excluded by the solvent are simply deleted. This leads to a smaller surface area and volume. It seems that the closure of concave regions is missing, and I could not find anything related to this in the code. Typically there is a fairly complex section in the code for this surface closure.

Results from other software and further discussion

Results from the lvpp sigma-profile database using GAMESS are shown below:

Results from the VT sigma database using DMol3 are shown below:

As can be seen these packages produce a closed SES surface. Apparently NWChem COSMO cavity construction options are not a SES surface. The YK method currently builds a VDW surface and the KS method apparently leaves an open cavity.

Input files

For the YK method:

start ACETONE
geometry units angstrom noautoz
  C          0.00000        0.00001        0.18289
  O         -0.00000       -0.00003        1.39646
  C          0.00000        1.29064       -0.61124
  C         -0.00000       -1.29062       -0.61121
  H         -0.00001        2.14877        0.07018
  H          0.00000       -2.14871        0.07026
  H          0.88462        1.33518       -1.26604
  H         -0.88462        1.33513       -1.26602
  H         -0.88462       -1.33516       -1.26600
  H          0.88462       -1.33513       -1.26598
end

basis noprint 
  * library 6-31G
end

dft
 XC b3lyp
 direct
 noio
end

cosmo
  rsolv  1.4
  screen ideal
# do_cosmo_ks
  minbem 3
  ificos 1
  
  do_gasphase false
end

task dft energy

For the KS method, just remove the comment for the do_cosmo_ks keyword.

[dmejiar]

@rpseng, given that the surface area generated with the YK scheme is within 2% of the other two packages, I would say that there seems to be a bug in the way that NWChem computes the volume. I took a quick look and it seems that the only volume being considered is the one that can be attributed to the "surviving" surface segments. This means that many contributions to the system's volume would be left out when more than one atom is present in the system. Hence the very small volume as compared to the other two codes.

[rpseng]

There might be a bug in the YK volume calculation, but based on my observations, this seems to be due to the use of different surface definitions. Typically, COSMO codes should compute a Solvent-Excluded Surface (SES), whereas it appears that the YK method is computing a van der Waals (vdW) surface instead.

The SES generally has a lower surface area but a greater volume compared to the vdW surface. Using a different type of surface could pose challenges when integrating it with modern quasi-chemical methods like COSMO-RS or COSMO-SAC. The theory assumes that molecules should make contact at every surface segment, but this is not feasible in regions of the vdW surface that are inaccessible to the solvent.

The issue with the current KS implementation appears to be more pronounced, where the cavities seem to remain open, suggesting a potential inconsistency in the approach.

I could try to work out another surface generation method or fix the KS method, but that can take a while.

[rpseng]

On that topic, if we run the KS method and set rsol 0.0, then it should reduce to a vdW surface. In this case the area and volume calculated are very close to what the YK produced:

molecular surface =    105.409 angstrom**2
molecular volume  =     59.151 angstrom**3

[rpseng]

Another example that makes the issues more clear, triethylamine. Again, the input file follows below.

KS method

The figure below shows the results using the KS method. The small spheres (or dots) represent the actual surface segments being considered:

As seen in the image, the surface segments appear only on part of the molecule, resulting in an open cavity. All segments in portions not accessible to the solvent were deleted.

Consequently, the calculated surface area and volume are significantly smaller than expected. I would expect great impact in the solvation properties calculated as well as potential convergence problems. Perhaps that is why rsolv=0.5 is the default value (usually it is ~1.3), as it reduces the open surface portions.

YK method

Results for the YK method:

Here the cavity is closed, but it corresponds to a van der Waals cavity rather than the expected solvent excluded surface. While the computed surface area seems reasonable, the volume is much smaller than anticipated, indicating a bug.

Input file

start TEA
geometry units angstrom noautoz
  N          0.35095       -0.19083        0.00000
  C          0.80141        0.48442       -1.21210
  C          0.80141        0.48442        1.21210
  C         -1.05008       -0.60654        0.00000
  C          0.65961       -0.34987       -2.48430
  C          0.65961       -0.34987        2.48430
  C         -2.11074        0.50738        0.00000
  H          0.30241        1.47057       -1.36459
  H          1.87002        0.71782       -1.06838
  H          1.87002        0.71782        1.06838
  H          0.30241        1.47057        1.36459
  H         -1.21294       -1.25686       -0.87419
  H         -1.21294       -1.25686        0.87419
  H          1.16710        0.15275       -3.32290
  H          1.11660       -1.34294       -2.34689
  H         -0.39042       -0.49592       -2.78244
  H          1.16710        0.15275        3.32290
  H         -0.39042       -0.49592        2.78244
  H          1.11660       -1.34294        2.34689
  H         -2.03319        1.15374        0.88917
  H         -2.03319        1.15374       -0.88917
  H         -3.12145        0.06765        0.00000
end

basis noprint 
  * library 6-31G
end

dft
 XC b3lyp
 direct
 noio
end

cosmo
  rsolv  1.4
  screen ideal
# do_cosmo_ks
  minbem 3
  ificos 1
  
  do_gasphase false
end

task dft energy

[dmejiar]

@rpseng I have not had time to look at the code in detail but I found a couple of things that might explain part of the behavior that you are seeing.

The first thing to note is that for the KS case, NWChem generates a SAS using Ratm + Rsolv. However, NWChem computes and prints the VdW surface and volume by projecting back the points on the SAS to the VdW surface (it uses only Ratm to compute the volume and surface). I believe this procedure is the one causing the rather small surfaces, since the back projection would lead to open spacew precisely on the surfaces near intersections. Printing of the surface points might also use the same back projection, although I haven't checked if this is really the case.

Another issue is the one I already commented for the volume calculation: NWChem only takes into account the contributions form the surface points, but neglects the inner parts of the atomic spheres.

For YK, the story is a little bit more complicated. The points are indeed located at the VdW surface, but their contribution to the surface (and volume) is scaled by the YK switching function (points around intersecting regions contribute less to the area/volume). As I see it, this procedure should lead to a surface area in between that of VdW and SES, with more aggressive cutoff values perhaps leading to surfaces closer to SES (although this is just an hypothesis). Again, volume calculation neglects the contribution from the inner regions, hence the very small values.

[edoapra]

Thank you very much for this very useful analysis of the cavities produced by NWChem's solvation module.
I agree that the triethylamine case clearly shows the deficiencies of the point charges' surfaces current implemented.
If you can work out another surface generation method, that would be an extremely useful contribution. We can assist you in this process, of course. The only catch is that I would ask you to leave the option of generate cavities using the current methods, too, for keeping the availability of reproducing the current results in the future.

[rpseng]

Thank you for your comment. I’ve managed to generate an SES surface and am now comparing the results. I believe I can submit a PR soon.