Solvation energy depend on cavity tessera area.

Sungwoo-Kang commented 6 years ago

Hello. I greatly appreciate for such a helpful library!

I have figured out that the solvation energy discontinuously bumps as area variable changes. I think it is natural that the solvation energy depends on the number of quadrature points by ~1 kcal/mol, but these bumping behavior seems unnatural, especially the solvation energy bumps about 100 Ha and then recovers as the number of finite elements decreases. I wrote detailed descriptions below. I appreciate any suggestions.

Current Behavior

Solvation energy discontinuously bumps as area variable changes. Following table contains the solvation energy and number of finite elements with varying partitioned surface area input. Test molecule is the bare nuclei of biphenyl molecule. Similar tendency occurs for various radii set, scaling option, compilers (intel 2015, GCC 5.3.0) and optimization options.

Area (input)	Number of finite elements	Solvation energy
0.17-0.21	1840	-478.1969864
0.22-0.26	1512	-379.1572135
0.27-0.32	1328	-378.6571250
0.33-0.35	976	-477.1193199

Running with single hydrogen nucleus shows ~10 mHa order of difference, which does not seems suspicious.

Steps to Reproduce (for bugs)

Run go-pcm.py with following input, with varying {{ area }} from 0.17 to 0.30.

units = angstrom
MEDIUM
{
    solvertype = cpcm
    solvent = water
}

CAVITY
{
    type = gepol
    area = {{ area }}
    radiiset = bondi
    mode = implicit
}

MOLECULE
{
    geometry = [
    0.3731926992,   0.6469064084,   0.0000000000,  6,
   -0.3731920247,  -0.6469053527,   0.0000000000,  6,
    1.7876957010,   0.6875961212,   0.0000000000,  6,
   -1.7876948822,  -0.6875966348,   0.0000000000,  6,
   -0.2987095886,   1.8919872045,   0.0000000000,  6,
    0.2987087796,  -1.8919868417,   0.0000000000,  6,
    2.4890930848,   1.8980892202,   0.0000000000,  6,
   -2.4890928933,  -1.8980892208,   0.0000000000,  6,
    0.3992471290,   3.1045982298,   0.0000000000,  6,
   -0.3992472455,  -3.1045981328,   0.0000000000,  6,
    1.8007582144,   3.1184678710,   0.0000000000,  6,
   -1.8007582843,  -3.1184679570,   0.0000000000,  6,
    2.3666720015,  -0.2395767444,   0.0000000000,  1,
   -2.3666717003,   0.2395758360,   0.0000000000,  1,
   -1.3911689868,   1.9277753687,   0.0000000000,  1,
    1.3911680825,  -1.9277757862,   0.0000000000,  1,
    3.5842514984,   1.8839459229,   0.0000000000,  1,
   -3.5842512870,  -1.8839456274,   0.0000000000,  1,
   -0.1599871452,   4.0461738267,   0.0000000000,  1,
    0.1599870403,  -4.0461737001,   0.0000000000,  1,
    2.3479356426,   4.0664812050,   0.0000000000,  1,
   -2.3479358356,  -4.0664812165,   0.0000000000,  1]
}

Your Environment

Version used: v1.2.1 and master(1c2e836)
Build option: -fc=gfortran --cc=gcc --cxx=g++ --type=relwithdebinfo (GCC) and --fc=ifort --cc=icc --cxx=icxx --omp --type=relwithdebinfo --extra-cc-flags="-O0 -std=c++11" (Intel)

Sungwoo-Kang commented 6 years ago

Uh, are there anything for me to try? I appreciate for any comments and have no idea how to deal with this.

robertodr commented 6 years ago

Ooops, sorry for the slow response! We haven't had time yet to look at this in detail, we will report as soon as we have a lead.

Sungwoo-Kang commented 6 years ago

Thank you for kind reply! I feel sorry that I might have interrupted you.

robertodr commented 5 years ago

This went completely under by radar for a long time, quite embarassing.

I've had a look at I think it all stems from a bad cavity tesselation. These cases are very hard to predict in general and even harder to correct in any meaningful way once they occur. In PR #170 I've put in place some checks to detect problems and stop the calculation before it outputs utter nonsense. For your 4 different tesselations I get:

0.17-0.21: Gauss' theorem check fails with a 0.016459 absolute difference.
0.22-0.26 S matrix is not positive-definite.
0.27-0.32 S matrix is not positive-definite.
0.33-0.35: Gauss' theorem check fails with a 0.0389664 absolute difference.

Whereas I think cases 1 and 4 are now artifactual failure due to the too stringent condition on the Gauss' theorem check, cases 2 and 3 are failures due to artifacts in the cavity formation.

Again, apologies for the extremely slow reply to this issue and for not being of more help!

Sungwoo-Kang commented 5 years ago

Thank you for your attention! I have tested a latest commit with a few molecules, and found the program halts or gives good solvation free energy values. Most difficult case was diphenyl molecule, but it works well enough after the Gauss' theorem check threshold changed to relative one. Thanks!

robertodr commented 5 years ago

Thanks for trying it out! It's not a solution, but at least gives the users a chance to understand what's going wrong. We'll mint a new minor release soon.

PCMSolver / pcmsolver