PCMSolver fatal error. S matrix is not positive-definite!

[Moh-Pou]

Getting an error in pcm-enabled energy calculation:

Expected Behavior

There is no error when inputting a very small molecule such as below. The energy is calculated smoothly.

molecule h2o {
N     -0.0000000001    -0.1040380466      0.0000000000
H     -0.9015844116     0.4818470201     -1.5615900098
H     -0.9015844116     0.4818470201      1.5615900098
H      1.8031688251     0.4818470204      0.0000000000
units bohr
}

Current Behavior

When I change the molecule to something slightly bigger, I get the following error. The error reads:

PCMSolver fatal error.
In function computeS at line 51 of file /home/conda/feedstock_root/build_artifacts/pcmsolver-split_1685512869815/work/src/bi_operators/IBoundaryIntegralOperator.cpp
S matrix is not positive-definite!
Consider changing the average area of the cavity finite elements.
Please report this issue: https://github.com/PCMSolver/pcmsolver/issues

Attached please find the outputs.
output.txt
output.log

Steps to Reproduce (for bugs)

Here is the input (named input):

memory 350 GB
molecule comp {
0 1
C    -0.754200    -2.008900     0.801900
H    -2.172800    -5.040100     0.661700
H    -0.579900     3.155700    -1.867900
N    -0.082900    -4.208100     0.405800
H    -1.199800    -0.224000     1.896900
C    -1.522400    -4.148900     0.662600
H     0.761100     1.339300    -2.776700
H    -1.738500     2.434400    -2.916100
S    -2.446300    -2.627900     1.016100
H     1.334400    -2.473100     0.381300
H    -0.703300     0.479700    -3.307700
H    -1.562700     1.885800    -1.145500
C    -1.087600     2.212600    -2.070400
C    -0.229600     0.999400    -2.474900
C     0.289900    -2.796500     0.521000
C     0.803400    -0.044800     1.005500
C    -0.667400    -0.481900     0.981400
C    -1.234600     0.205700    -0.262400
C    -0.100500     0.040700    -1.276600
C     1.173600     0.276500    -0.453800
H     0.912200     0.871200     1.635600
H     1.458200    -0.834900     1.445900
H    -2.196200    -0.244700    -0.607900
H    -1.404300     1.289900    -0.053800
H    -0.108400    -1.005300    -1.669200
H     1.488000     1.345000    -0.540800
H     2.018800    -0.352600    -0.823700
}

set {
  basis aug-cc-pvdz
  e_convergence 10
  d_convergence 10
  scf_type pk
  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
  }
}

# Explicitly ask for numerical
thisenergy = optimize('b3lyp', molecule=comp, dertype='energy')
comp.save_xyz_file('name.xyz',0)

How to run: psi4 input -n 16 output.txt

Context

Your Environment

• Version used: psi4 1.8.1
• Environment name and version (e.g. PHP 5.4 on nginx 1.9.1): I don't think this is applicable.
• Server type and version: I don't think this is applicable.
• Operating System and version: Ubuntu 20.04.6 LTS

[Moh-Pou]

I wanted to gently follow up about this issue. If you could spare some time to review and provide your feedback or guidance, it would be greatly appreciated.

[markperri]

I'm getting the same issue with 1.8.1 on Rocky8. I've tried different values for area (0.1, 0.3, 0.5), but I still get the error.

[welchbr]

I'm getting the same issue albeit with MRCC. Same exact error and I've changed the area from really small (0.1) to large ones like 5.0.

[jgonthier]

Well hello everyone, I am also getting the same issue with Psi4 1.8.2. Tried a number of different cavity areas. I can only speculate but it looks like it's a recent problem, however the code in this repo hasn't changed since 2021. Maybe it's coming from a dependency? Maybe earlier versions of Psi4 would work?

[jgonthier]

In case anybody is interested: using Psi4 1.6 I was able to get the solvation to work. I think it's because Psi4 1.6 uses a previous version of the PCMSolver library.

[loriab]

Thanks for the report, and I'm sorry for the trouble. Looking a little closer, loriab/pcmsolver@v1211...loriab:pcmsolver:v123_plus_ming is the diff for the pcmsolver used with psi4 v1.8 (v1.2.3-ish) vs. earlier (v.1.2.1-ish). Skipping over all the frills, there's only three changes of substance:

• Gauss check on nuclear charges, loriab/pcmsolver@v1211...loriab:pcmsolver:v123_plus_ming#diff-8d11e793b6dd5e60de9e3e8d113979ae53e634b6651aed4d1342a81985481048R547
• Pruning the tesserae of little contributions that might break pos-def-ness, loriab/pcmsolver@v1211...loriab:pcmsolver:v123_plus_ming#diff-9bb493e48722de80574a99af86291d9a574403d31c862161e5e29f4c0a41b89cL239
• Acutally checking pos-def-ness loriab/pcmsolver@v1211...loriab:pcmsolver:v123_plus_ming#diff-a7ccae0968d4e17551af3b0886f63f2b171c3d32bae13d848fe3418bfbb7d9f4R43
• changelog: loriab/pcmsolver@v1211...loriab:pcmsolver:v123_plus_ming#diff-06572a96a58dc510037d5efa622f9bec8519bc1beab13c9f251e97e657a9d4ed

So it looks like the reason you were seeing success with psi4 1.6 is b/c the pos-def check hasn't been implemented yet at pcm ~1.2.1.

I could do a 1.2.1 on c-f and link psi4 to it, but that does seem to be heading backwards. Maybe for the larger mols you're using, the 1.0e-4 filter (loriab/pcmsolver@v1211...loriab:pcmsolver:v123_plus_ming#diff-9bb493e48722de80574a99af86291d9a574403d31c862161e5e29f4c0a41b89cR252) is too small to eliminate all the tesserae contaminating positive-definiteness? I guess I could try to make that an adjustable parameter. Thoughts, @robertodr ? It'd be good to find a strategy by the end of the month so that psi4 v1.9 in early Dec can use the new version.

[jgonthier]

Thank you for the details, Lori! It seems the lack of positive-definiteness in 1.6 was causing some issues in SCF convergence, maybe. I was also looking at a weird system.

[loriab]

Hi Jerome! It seems like the "area" knob isn't enough to solve some of the errors being thrown. I'm not sure if optionally turning off that warning and letting SCF take its chances or if there's another knob that can be exposed. Hoping @robertodr can provide some guidance.

It doesn't look like there's any changes in v1.3 that would address this problem.

[metma99]

Hello:
Did anybody find a solution to this issue? I have installed the latest psi4 version 1.9.1.
Any pointers are highly appreciated.
Best regards,
Markus

[amcisaac]

I am also having this issue (Psi4 1.9.1), have there been any updates? I tried using 0.3 AU, 0.1 AU, and 0.05 AU for the area but to no avail. Is there any guidance on how serious the error is, e.g. can we trust the results?

[loriab]

tl;dr I'll convert the error to warning in the conda package, but it's not best advised. Consider pyddx instead of PCMSolver.

Hi all, I've been in touch with Roberto DiRemegio on this issue. He says the S-not-positive-definite is a real issue potentially producing nonsense if ignored. But he mentions another check that might take a little more physics into account so as to be a better gauge for sensible results. He's fine with me converting these to bold warnings but does implore you all to watch for bad results. (The current state on conda is that the latest pcm package has the main (S matrix) error as warning, but I don't know if you get that package with the latest psi4, as I haven't rebuilt that.) The problem can't really be healed because it's coming from legacy Fortran PCM code that constructs the cavity and wasn't designed for larger molecules. He suggests that Michael Herbst's pyddx+psi4 interface might be the way to go test case.

[loriab]

Update that the latest PCMSolver and Psi4 packages have both the errors converted to warnings. They're based on v1.2.3 .

[proteins247]

Hi all, @loriab, I wanted to raise attention for my issue, #207, where I notice that the cavity generated by PCMSolver using Psi4 1.9.1 is incorrect. Did some of the changes introduced into the most recent updates to PCMSolver introduce some sort of bug that at the very least makes the generated cavity file to be incorrect?

I fear that the update mentioned above that converts errors into warnings doesn't address what seems to be a regression in the software.