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

Moh-Pou commented 1 year ago

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 commented 11 months ago

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 commented 11 months ago

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 commented 11 months ago

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 commented 10 months ago

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 commented 10 months ago

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 commented 10 months ago

Thanks for the report, and I'm sorry for the trouble. Looking a little closer, https://github.com/loriab/pcmsolver/compare/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, https://github.com/loriab/pcmsolver/compare/v1211...loriab:pcmsolver:v123_plus_ming#diff-8d11e793b6dd5e60de9e3e8d113979ae53e634b6651aed4d1342a81985481048R547
Pruning the tesserae of little contributions that might break pos-def-ness, https://github.com/loriab/pcmsolver/compare/v1211...loriab:pcmsolver:v123_plus_ming#diff-9bb493e48722de80574a99af86291d9a574403d31c862161e5e29f4c0a41b89cL239
Acutally checking pos-def-ness https://github.com/loriab/pcmsolver/compare/v1211...loriab:pcmsolver:v123_plus_ming#diff-a7ccae0968d4e17551af3b0886f63f2b171c3d32bae13d848fe3418bfbb7d9f4R43
changelog: https://github.com/loriab/pcmsolver/compare/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 (https://github.com/loriab/pcmsolver/compare/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 commented 10 months ago

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 commented 9 months ago

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 commented 5 months ago

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 commented 2 days ago

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 commented 1 day ago

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.

PCMSolver / pcmsolver