search

PCMSolver / pcmsolver

An API for the Polarizable Continuum Model
http://pcmsolver.readthedocs.io/
GNU Lesser General Public License v3.0
32 stars 21 forks source link

Psi4 and PCMSolver v1.3 and UFF radii #201

Open loriab opened 1 year ago

loriab commented 1 year ago

Hi @robertodr, don't worry, this isn't a new error. This is an update on that old trouble of PCMSolver v1.3 failing on several of the Psi4 tests; you and I both tried to solve it a couple years ago to no avail. Psi4 since has been using v1.2.1 plus 3 parsing commits, informally dubbed v1.2.1.1.

Anyways, I've been poking at these problems preliminary to putting pcmsolver on conda-forge. Is that ok in general and do you want to be a maintainer on the recipe? I got Psi4 using v1.2.3 with very small tweaks (~1.e-8) to test cases after you introduced the PEDRA pruning. (And I added the citation printing.) So continuing with the v1.2.x branch is one option. Then today, I gave v1.3.0 a try again backing out the 1.0 -> 1.1 UFF radii conversion https://github.com/PCMSolver/pcmsolver/blob/master/src/utils/Atom.cpp#L182 and happily, that heals all the Psi4+PCMSolver tests. Note that the radii change was influencing energies in the third decimal place (see error log below). Does that seem right order-of-magnitude? So updating Psi4 ref values and switching to v1.3.x branch is a viable option.

I quite understand that PCMSolver is under minimal development, and I'm not pressuring otherwise. I wanted to see if you had any advice for which branch to pursue for Psi4 and for c-f packaging. Normally, I'd say v1.3.x, of course, but I don't know if it got as much exercise downstream as v1.2.x. I'd be glad of any advice.

Start testing: Mar 14 11:10 EDT
----------------------------------------------------------
525/533 Testing: pcmsolver-dft
525/533 Test: pcmsolver-dft
Command: "/psi/toolchainconda/envs/defenv10/bin/python3.10" "/psi/gits/hrw-detangle/tests/runtest.py" "/psi/gits/hrw-detangle/tests/pcmsolver/dft/input.dat" "/psi/gits/hrw-detangle/objdir_defenv10/testresults.log" "/psi/gits/hrw-detangle" "/psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-dft/output.dat" "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/share/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/"
Directory: /psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-dft
"pcmsolver-dft" start time: Mar 14 11:10 EDT
Output:
----------------------------------------------------------
RKS-PCM B3LYP, total algorithm
    Nuclear repulsion energy (PCM, total algorithm).......................................PASSED
    Total energy (PCM, total algorithm)...................................................FAILED
Traceback (most recent call last):
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4", line 338, in <module>
    exec(content)
  File "<string>", line 41, in <module>
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/qcdb/testing.py", line 104, in _mergedapis_compare_values
    return qcel.testing.compare_values(expected, computed, **kwargs)
  File "/psi/toolchainconda/envs/defenv10/lib/python3.10/site-packages/qcelemental/testing.py", line 178, in compare_values
    return return_handler(allclose, label, message, return_message, quiet)
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/p4util/testing.py", line 204, in _psi4_true_raise_handler
    raise TestComparisonError(message)

psi4.driver.p4util.exceptions.TestComparisonError:  Total energy (PCM, total algorithm): computed value (-56.55572018797) does not match (-56.55724704841) to atol=1e-09 by difference (0.00152686044).

Printing out the relevant lines from the Psithon --> Python processed input file:
        tmp.write('\n'.join(['SECT toplevel 5 False', 'TAG F KW 2', 'STR UNITS 1 True', 'ANGSTROM', 'INT CODATA 1 False', '2010', 'SECT MEDIUM 1 True', 'TAG F KW 8', 'STR SOLVERTYPE 1 True', 'IEFPCM', 'STR SOLVENT 1 True', 'WATER', 'BOOL NONEQUILIBRIUM 1 False', 'False', 'BOOL MATRIXSYMM 1 False', 'True', 'DBL CORRECTION 1 False', '0.0', 'STR DIAGONALINTEGRATOR 1 False', 'COLLOCATION', 'DBL DIAGONALSCALING 1 False', '1.07', 'DBL PROBERADIUS 1 False', '1.8897261245650618', 'SECT GREEN 1 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT GREEN 0 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT CAVITY 0 True', 'TAG F KW 10', 'STR RADIISET 1 True', 'UFF', 'STR TYPE 1 True', 'GEPOL', 'BOOL SCALING 1 True', 'False', 'DBL AREA 1 True', '1.0713194477591061', 'STR MODE 1 True', 'IMPLICIT', 'STR NPZFILE 1 False', '', 'DBL MINRADIUS 1 False', '100.0', 'INT_ARRAY ATOMS 0 False', 'DBL_ARRAY RADII 0 False', 'DBL_ARRAY SPHERES 0 False', 'SECT MOLECULE 0 False', 'TAG F KW 2', 'DBL_ARRAY GEOMETRY 0 False', 'BOOL MEP 1 False', 'True', 'SECT CHARGEDISTRIBUTION 0 False', 'TAG F KW 2', 'DBL_ARRAY MONOPOLES 0 False', 'DBL_ARRAY DIPOLES 0 False', 'SECT MMFQ 0 False', 'TAG F KW 3', 'BOOL NONPOLARIZABLE 1 False', 'False', 'INT SITESPERFRAGMENT 1 False', '3', 'DBL_ARRAY SITES 0 False']))
    core.set_local_option('PCM', 'PCMSOLVER_PARSED_FNAME', '@pcmsolver.10298.d0840893')
    print('RKS-PCM B3LYP, total algorithm')
    energy_scf1 = energy('b3lyp')
    compare_values(nucenergy, NH3.nuclear_repulsion_energy(), 10, "Nuclear repulsion energy (PCM, total algorithm)")
--> compare_values(totalenergy, energy_scf1, 9, "Total energy (PCM, total algorithm)")
    compare_values(polenergy, psi4.core.variable("PCM POLARIZATION ENERGY"), 6, "Polarization energy (PCM, total algorithm)")
    core.set_global_option("PCM_SCF_TYPE", "separate")
    print('RKS-PCM B3LYP, separate algorithm')
    energy_scf2 = energy('b3lyp')
    compare_values(totalenergy, energy_scf2, 9, "Total energy (PCM, separate algorithm)")

!----------------------------------------------------------------------------------!
!                                                                                  !
!         Total energy (PCM, total algorithm): computed value (-56.55572018797)    !
!     does not match (-56.55724704841) to atol=1e-09 by difference                 !
!     (0.00152686044).                                                             !
!                                                                                  !
!----------------------------------------------------------------------------------!

Exit Status: overall (1)
<end of output>
Test time =   3.54 sec
----------------------------------------------------------
Test Failed.
"pcmsolver-dft" end time: Mar 14 11:10 EDT
"pcmsolver-dft" time elapsed: 00:00:03
----------------------------------------------------------

526/533 Testing: pcmsolver-dipole
526/533 Test: pcmsolver-dipole
Command: "/psi/toolchainconda/envs/defenv10/bin/python3.10" "/psi/gits/hrw-detangle/tests/runtest.py" "/psi/gits/hrw-detangle/tests/pcmsolver/dipole/input.dat" "/psi/gits/hrw-detangle/objdir_defenv10/testresults.log" "/psi/gits/hrw-detangle" "/psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-dipole/output.dat" "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/share/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/"
Directory: /psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-dipole
"pcmsolver-dipole" start time: Mar 14 11:10 EDT
Output:
----------------------------------------------------------
    Energy for displacement 0.............................................................FAILED
Traceback (most recent call last):
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4", line 338, in <module>
    exec(content)
  File "<string>", line 64, in <module>
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/qcdb/testing.py", line 104, in _mergedapis_compare_values
    return qcel.testing.compare_values(expected, computed, **kwargs)
  File "/psi/toolchainconda/envs/defenv10/lib/python3.10/site-packages/qcelemental/testing.py", line 178, in compare_values
    return return_handler(allclose, label, message, return_message, quiet)
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/p4util/testing.py", line 204, in _psi4_true_raise_handler
    raise TestComparisonError(message)

psi4.driver.p4util.exceptions.TestComparisonError:  Energy for displacement 0: computed value (-76.0177183642) does not match (-76.0203919443) to atol=1e-08 by difference (0.0026735801).

Printing out the relevant lines from the Psithon --> Python processed input file:
        dm_z_5point = (8.0*energies[0] - 8.0*energies[1] - energies[2] + energies[3]) / (12.0*pert)

        for val in range(len(lambdas)):
            print_out("Perturbation strength = %7.4f, %s computed energy = %16.10f\n" % (lambdas[val], m, energies[val]))
-->         compare_values(ref_energies[m][val], energies[val], 8, 'Energy for displacement %d' % val)
        print_out("Total dipoles\n")
        print_out("%s 3-point stencil: mu(z) = %10.10f ea0, %10.10f Debye\n" % (m, dm_z_3point, dm_z_3point*psi4.constants.dipmom_au2debye))
        print_out("%s 5-point stencil: mu(z) = %10.10f ea0, %10.10f Debye\n" % (m, dm_z_5point, dm_z_5point*psi4.constants.dipmom_au2debye))

        core.set_global_option("PERTURB_H", "false")

!----------------------------------------------------------------------------------!
!                                                                                  !
!         Energy for displacement 0: computed value (-76.0177183642) does not      !
!     match (-76.0203919443) to atol=1e-08 by difference (0.0026735801).           !
!                                                                                  !
!----------------------------------------------------------------------------------!

Exit Status: overall (1)
<end of output>
Test time =   3.47 sec
----------------------------------------------------------
Test Failed.
"pcmsolver-dipole" end time: Mar 14 11:10 EDT
"pcmsolver-dipole" time elapsed: 00:00:03
----------------------------------------------------------

527/533 Testing: pcmsolver-scf
527/533 Test: pcmsolver-scf
Command: "/psi/toolchainconda/envs/defenv10/bin/python3.10" "/psi/gits/hrw-detangle/tests/runtest.py" "/psi/gits/hrw-detangle/tests/pcmsolver/scf/input.dat" "/psi/gits/hrw-detangle/objdir_defenv10/testresults.log" "/psi/gits/hrw-detangle" "/psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-scf/output.dat" "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/share/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/"
Directory: /psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-scf
"pcmsolver-scf" start time: Mar 14 11:10 EDT
Output:
----------------------------------------------------------
    Nuclear repulsion energy (PCM, total algorithm).......................................PASSED
    Total energy (PCM, total algorithm)...................................................FAILED
Traceback (most recent call last):
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4", line 338, in <module>
    exec(content)
  File "<string>", line 45, in <module>
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/qcdb/testing.py", line 104, in _mergedapis_compare_values
    return qcel.testing.compare_values(expected, computed, **kwargs)
  File "/psi/toolchainconda/envs/defenv10/lib/python3.10/site-packages/qcelemental/testing.py", line 178, in compare_values
    return return_handler(allclose, label, message, return_message, quiet)
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/p4util/testing.py", line 204, in _psi4_true_raise_handler
    raise TestComparisonError(message)

psi4.driver.p4util.exceptions.TestComparisonError:  Total energy (PCM, total algorithm): computed value (-55.454581206064) does not match (-55.455942603983) to atol=1e-10 by difference (0.001361397919).

Printing out the relevant lines from the Psithon --> Python processed input file:
        tmp.write('\n'.join(['SECT toplevel 5 False', 'TAG F KW 2', 'STR UNITS 1 True', 'ANGSTROM', 'INT CODATA 1 False', '2010', 'SECT MEDIUM 1 True', 'TAG F KW 8', 'STR SOLVERTYPE 1 True', 'IEFPCM', 'STR SOLVENT 1 True', 'WATER', 'BOOL NONEQUILIBRIUM 1 False', 'False', 'BOOL MATRIXSYMM 1 False', 'True', 'DBL CORRECTION 1 False', '0.0', 'STR DIAGONALINTEGRATOR 1 False', 'COLLOCATION', 'DBL DIAGONALSCALING 1 False', '1.07', 'DBL PROBERADIUS 1 False', '1.8897261245650618', 'SECT GREEN 1 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT GREEN 0 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT CAVITY 0 True', 'TAG F KW 10', 'STR RADIISET 1 True', 'UFF', 'STR TYPE 1 True', 'GEPOL', 'BOOL SCALING 1 True', 'False', 'DBL AREA 1 True', '1.0713194477591061', 'STR MODE 1 True', 'IMPLICIT', 'STR NPZFILE 1 False', '', 'DBL MINRADIUS 1 False', '100.0', 'INT_ARRAY ATOMS 0 False', 'DBL_ARRAY RADII 0 False', 'DBL_ARRAY SPHERES 0 False', 'SECT MOLECULE 0 False', 'TAG F KW 2', 'DBL_ARRAY GEOMETRY 0 False', 'BOOL MEP 1 False', 'True', 'SECT CHARGEDISTRIBUTION 0 False', 'TAG F KW 2', 'DBL_ARRAY MONOPOLES 0 False', 'DBL_ARRAY DIPOLES 0 False', 'SECT MMFQ 0 False', 'TAG F KW 3', 'BOOL NONPOLARIZABLE 1 False', 'False', 'INT SITESPERFRAGMENT 1 False', '3', 'DBL_ARRAY SITES 0 False']))
    core.set_local_option('PCM', 'PCMSOLVER_PARSED_FNAME', '@pcmsolver.10304.2cddce17')
    print_out('PK-RHF-PCM, total algorithm')
    energy_scf, wfn = energy('scf', return_wfn=True)
    compare_values(nucenergy, NH3.nuclear_repulsion_energy(), 10, "Nuclear repulsion energy (PCM, total algorithm)")
--> compare_values(totalenergy, energy_scf, 10, "Total energy (PCM, total algorithm)")
    compare_values(polenergy, wfn.variable("PCM POLARIZATION ENERGY"), 6, "Polarization energy (PCM, total algorithm)")
    core.set_global_option("PCM_SCF_TYPE", "separate")
    print_out('PK-RHF-PCM, separate algorithm')
    energy_scf, wfn = energy('scf', return_wfn=True)
    compare_values(totalenergy, energy_scf, 10, "Total energy (PCM, separate algorithm)")

!----------------------------------------------------------------------------------!
!                                                                                  !
!         Total energy (PCM, total algorithm): computed value (-55.454581206064)   !
!     does not match (-55.455942603983) to atol=1e-10 by difference                !
!     (0.001361397919).                                                            !
!                                                                                  !
!----------------------------------------------------------------------------------!

Exit Status: overall (1)
<end of output>
Test time =   1.88 sec
----------------------------------------------------------
Test Failed.
"pcmsolver-scf" end time: Mar 14 11:10 EDT
"pcmsolver-scf" time elapsed: 00:00:01
----------------------------------------------------------

528/533 Testing: pcmsolver-opt-fd
528/533 Test: pcmsolver-opt-fd
Command: "/psi/toolchainconda/envs/defenv10/bin/python3.10" "/psi/gits/hrw-detangle/tests/runtest.py" "/psi/gits/hrw-detangle/tests/pcmsolver/opt-fd/input.dat" "/psi/gits/hrw-detangle/objdir_defenv10/testresults.log" "/psi/gits/hrw-detangle" "/psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-opt-fd/output.dat" "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/share/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/"
Directory: /psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-opt-fd
"pcmsolver-opt-fd" start time: Mar 14 11:10 EDT
Output:
----------------------------------------------------------
Optimizer: Optimization complete!
    Nuclear repulsion energy..............................................................FAILED
Traceback (most recent call last):
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4", line 338, in <module>
    exec(content)
  File "<string>", line 38, in <module>
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/qcdb/testing.py", line 104, in _mergedapis_compare_values
    return qcel.testing.compare_values(expected, computed, **kwargs)
  File "/psi/toolchainconda/envs/defenv10/lib/python3.10/site-packages/qcelemental/testing.py", line 178, in compare_values
    return return_handler(allclose, label, message, return_message, quiet)
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/p4util/testing.py", line 204, in _psi4_true_raise_handler
    raise TestComparisonError(message)

psi4.driver.p4util.exceptions.TestComparisonError:  Nuclear repulsion energy: computed value (8.91975) does not match (8.92442) to atol=0.001 by difference (-0.00467).

Printing out the relevant lines from the Psithon --> Python processed input file:
    parsedFile = os.path.join(os.getcwd(), '@pcmsolver.10306.5faf08d7')
    with open(parsedFile, 'w') as tmp:
        tmp.write('\n'.join(['SECT toplevel 5 False', 'TAG F KW 2', 'STR UNITS 1 True', 'ANGSTROM', 'INT CODATA 1 False', '2010', 'SECT MEDIUM 1 True', 'TAG F KW 8', 'STR SOLVERTYPE 1 True', 'IEFPCM', 'STR SOLVENT 1 True', 'WATER', 'BOOL NONEQUILIBRIUM 1 False', 'False', 'BOOL MATRIXSYMM 1 False', 'True', 'DBL CORRECTION 1 False', '0.0', 'STR DIAGONALINTEGRATOR 1 False', 'COLLOCATION', 'DBL DIAGONALSCALING 1 False', '1.07', 'DBL PROBERADIUS 1 False', '1.8897261245650618', 'SECT GREEN 1 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT GREEN 0 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT CAVITY 0 True', 'TAG F KW 10', 'STR RADIISET 1 True', 'UFF', 'STR TYPE 1 True', 'GEPOL', 'BOOL SCALING 1 True', 'False', 'DBL AREA 1 True', '1.0713194477591061', 'STR MODE 1 True', 'IMPLICIT', 'STR NPZFILE 1 False', '', 'DBL MINRADIUS 1 False', '100.0', 'INT_ARRAY ATOMS 0 False', 'DBL_ARRAY RADII 0 False', 'DBL_ARRAY SPHERES 0 False', 'SECT MOLECULE 0 False', 'TAG F KW 2', 'DBL_ARRAY GEOMETRY 0 False', 'BOOL MEP 1 False', 'True', 'SECT CHARGEDISTRIBUTION 0 False', 'TAG F KW 2', 'DBL_ARRAY MONOPOLES 0 False', 'DBL_ARRAY DIPOLES 0 False', 'SECT MMFQ 0 False', 'TAG F KW 3', 'BOOL NONPOLARIZABLE 1 False', 'False', 'INT SITESPERFRAGMENT 1 False', '3', 'DBL_ARRAY SITES 0 False']))
    core.set_local_option('PCM', 'PCMSOLVER_PARSED_FNAME', '@pcmsolver.10306.5faf08d7')
    thisenergy = optimize('scf', molecule=h2o, dertype='energy')
--> compare_values(nucenergy, h2o.nuclear_repulsion_energy(), 3, "Nuclear repulsion energy")   
    compare_values(refenergy, thisenergy, 7, "Reference energy")                               
    clean()
    core.set_global_option("G_CONVERGENCE", "gau_tight")
    thisenergy = optimize('scf', molecule=h2o)
    compare_values(nucenergy, h2o.nuclear_repulsion_energy(), 3, "Nuclear repulsion energy")   

!----------------------------------------------------------------------------------!
!                                                                                  !
!         Nuclear repulsion energy: computed value (8.91975) does not match        !
!     (8.92442) to atol=0.001 by difference (-0.00467).                            !
!                                                                                  !
!----------------------------------------------------------------------------------!

Exit Status: overall (1)
<end of output>
Test time =  15.54 sec
----------------------------------------------------------
Test Failed.
"pcmsolver-opt-fd" end time: Mar 14 11:11 EDT
"pcmsolver-opt-fd" time elapsed: 00:00:15
----------------------------------------------------------

529/533 Testing: pcmsolver-ccsd-pte
529/533 Test: pcmsolver-ccsd-pte
Command: "/psi/toolchainconda/envs/defenv10/bin/python3.10" "/psi/gits/hrw-detangle/tests/runtest.py" "/psi/gits/hrw-detangle/tests/pcmsolver/ccsd-pte/input.dat" "/psi/gits/hrw-detangle/objdir_defenv10/testresults.log" "/psi/gits/hrw-detangle" "/psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-ccsd-pte/output.dat" "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/share/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/"
Directory: /psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-ccsd-pte
"pcmsolver-ccsd-pte" start time: Mar 14 11:11 EDT
Output:
----------------------------------------------------------
    Nuclear repulsion energy (PCM, total algorithm).......................................PASSED
    SCF Total energy (PCM, total algorithm)...............................................FAILED
Traceback (most recent call last):
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4", line 338, in <module>
    exec(content)
  File "<string>", line 41, in <module>
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/qcdb/testing.py", line 104, in _mergedapis_compare_values
    return qcel.testing.compare_values(expected, computed, **kwargs)
  File "/psi/toolchainconda/envs/defenv10/lib/python3.10/site-packages/qcelemental/testing.py", line 178, in compare_values
    return return_handler(allclose, label, message, return_message, quiet)
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/p4util/testing.py", line 204, in _psi4_true_raise_handler
    raise TestComparisonError(message)

psi4.driver.p4util.exceptions.TestComparisonError:  SCF Total energy (PCM, total algorithm): computed value (-199.452249572074) does not match (-199.515074015114) to atol=1e-10 by difference (0.062824443040).

Printing out the relevant lines from the Psithon --> Python processed input file:
    with open(parsedFile, 'w') as tmp:
        tmp.write('\n'.join(['SECT toplevel 5 False', 'TAG F KW 2', 'STR UNITS 1 True', 'ANGSTROM', 'INT CODATA 1 False', '2010', 'SECT MEDIUM 1 True', 'TAG F KW 8', 'STR SOLVERTYPE 1 True', 'IEFPCM', 'STR SOLVENT 1 True', 'WATER', 'BOOL NONEQUILIBRIUM 1 False', 'False', 'BOOL MATRIXSYMM 1 False', 'True', 'DBL CORRECTION 1 False', '0.0', 'STR DIAGONALINTEGRATOR 1 False', 'COLLOCATION', 'DBL DIAGONALSCALING 1 False', '1.07', 'DBL PROBERADIUS 1 False', '1.8897261245650618', 'SECT GREEN 1 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT GREEN 0 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT CAVITY 0 True', 'TAG F KW 10', 'STR RADIISET 1 True', 'UFF', 'STR TYPE 1 True', 'GEPOL', 'BOOL SCALING 1 True', 'False', 'DBL AREA 1 True', '3.5710648258636875', 'STR MODE 1 True', 'IMPLICIT', 'STR NPZFILE 1 False', '', 'DBL MINRADIUS 1 False', '100.0', 'INT_ARRAY ATOMS 0 False', 'DBL_ARRAY RADII 0 False', 'DBL_ARRAY SPHERES 0 False', 'SECT MOLECULE 0 False', 'TAG F KW 2', 'DBL_ARRAY GEOMETRY 0 False', 'BOOL MEP 1 False', 'True', 'SECT CHARGEDISTRIBUTION 0 False', 'TAG F KW 2', 'DBL_ARRAY MONOPOLES 0 False', 'DBL_ARRAY DIPOLES 0 False', 'SECT MMFQ 0 False', 'TAG F KW 3', 'BOOL NONPOLARIZABLE 1 False', 'False', 'INT SITESPERFRAGMENT 1 False', '3', 'DBL_ARRAY SITES 0 False']))
    core.set_local_option('PCM', 'PCMSOLVER_PARSED_FNAME', '@pcmsolver.10335.0ad295f2')
    energy_pte, wfn = energy('ccsd', return_wfn=True)
    compare_values(Mg.nuclear_repulsion_energy(), nucenergy, 10, "Nuclear repulsion energy (PCM, total algorithm)")
--> compare_values(scf_totalenergy, variable("SCF Total energy"), 10, "SCF Total energy (PCM, total algorithm)")
    compare_values(scf_polenergy, wfn.reference_wavefunction().variable("PCM POLARIZATION ENERGY"), 8, "SCF Polarization energy (PCM, total algorithm)")
    compare_values(pte_totalenergy, energy_pte, 10, "CCSD Total energy (PCM PTE algorithm)")
    compare_values(pte_correnergy, variable("CURRENT CORRELATION ENERGY"), 10, "CCSD Correlation energy (PCM PTE algorithm)")

!----------------------------------------------------------------------------------!
!                                                                                  !
!         SCF Total energy (PCM, total algorithm): computed value                  !
!     (-199.452249572074) does not match (-199.515074015114) to atol=1e-10 by      !
!     difference (0.062824443040).                                                 !
!                                                                                  !
!----------------------------------------------------------------------------------!

Exit Status: overall (1)
<end of output>
Test time =   2.15 sec
----------------------------------------------------------
Test Failed.
"pcmsolver-ccsd-pte" end time: Mar 14 11:11 EDT
"pcmsolver-ccsd-pte" time elapsed: 00:00:02
----------------------------------------------------------

530/533 Testing: pcmsolver-ghost
530/533 Test: pcmsolver-ghost
Command: "/psi/toolchainconda/envs/defenv10/bin/python3.10" "/psi/gits/hrw-detangle/tests/runtest.py" "/psi/gits/hrw-detangle/tests/pcmsolver/ghost/input.dat" "/psi/gits/hrw-detangle/objdir_defenv10/testresults.log" "/psi/gits/hrw-detangle" "/psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-ghost/output.dat" "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/share/psi4" "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/"
Directory: /psi/gits/hrw-detangle/objdir_defenv10/tests/pcmsolver/pcmsolver-ghost
"pcmsolver-ghost" start time: Mar 14 11:11 EDT
Output:
----------------------------------------------------------
    Nuclear repulsion energy (PCM, total algorithm).......................................PASSED
    Total energy (PCM, total algorithm)...................................................FAILED
Traceback (most recent call last):
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/bin/psi4", line 338, in <module>
    exec(content)
  File "<string>", line 47, in <module>
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/qcdb/testing.py", line 104, in _mergedapis_compare_values
    return qcel.testing.compare_values(expected, computed, **kwargs)
  File "/psi/toolchainconda/envs/defenv10/lib/python3.10/site-packages/qcelemental/testing.py", line 178, in compare_values
    return return_handler(allclose, label, message, return_message, quiet)
  File "/psi/gits/hrw-detangle/objdir_defenv10/stage/lib/psi4/driver/p4util/testing.py", line 204, in _psi4_true_raise_handler
    raise TestComparisonError(message)

psi4.driver.p4util.exceptions.TestComparisonError:  Total energy (PCM, total algorithm): computed value (-55.801284771319) does not match (-55.802798786575) to atol=1e-10 by difference (0.001514015256).

Printing out the relevant lines from the Psithon --> Python processed input file:
    with open(parsedFile, 'w') as tmp:
        tmp.write('\n'.join(['SECT toplevel 5 False', 'TAG F KW 2', 'STR UNITS 1 True', 'ANGSTROM', 'INT CODATA 1 False', '2010', 'SECT MEDIUM 1 True', 'TAG F KW 8', 'STR SOLVERTYPE 1 True', 'IEFPCM', 'STR SOLVENT 1 True', 'WATER', 'BOOL NONEQUILIBRIUM 1 False', 'False', 'BOOL MATRIXSYMM 1 False', 'True', 'DBL CORRECTION 1 False', '0.0', 'STR DIAGONALINTEGRATOR 1 False', 'COLLOCATION', 'DBL DIAGONALSCALING 1 False', '1.07', 'DBL PROBERADIUS 1 False', '1.8897261245650618', 'SECT GREEN 1 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT GREEN 0 False', 'TAG F KW 13', 'STR TYPE 1 False', 'VACUUM', 'STR DER 1 False', 'DERIVATIVE', 'DBL EPS 1 False', '1.0', 'DBL EPSDYN 1 False', '1.0', 'STR PROFILE 1 False', 'LOG', 'DBL EPS1 1 False', '1.0', 'DBL EPSDYN1 1 False', '1.0', 'DBL EPS2 1 False', '1.0', 'DBL EPSDYN2 1 False', '1.0', 'DBL CENTER 1 False', '188.9726124565062', 'DBL WIDTH 1 False', '9.44863062282531', 'DBL_ARRAY INTERFACEORIGIN 3 False', '0.0', '0.0', '0.0', 'INT MAXL 1 False', '50', 'SECT CAVITY 0 True', 'TAG F KW 10', 'STR RADIISET 1 True', 'UFF', 'STR TYPE 1 True', 'GEPOL', 'BOOL SCALING 1 True', 'False', 'DBL AREA 1 True', '1.0713194477591061', 'STR MODE 1 True', 'IMPLICIT', 'STR NPZFILE 1 False', '', 'DBL MINRADIUS 1 False', '100.0', 'INT_ARRAY ATOMS 0 False', 'DBL_ARRAY RADII 0 False', 'DBL_ARRAY SPHERES 0 False', 'SECT MOLECULE 0 False', 'TAG F KW 2', 'DBL_ARRAY GEOMETRY 0 False', 'BOOL MEP 1 False', 'True', 'SECT CHARGEDISTRIBUTION 0 False', 'TAG F KW 2', 'DBL_ARRAY MONOPOLES 0 False', 'DBL_ARRAY DIPOLES 0 False', 'SECT MMFQ 0 False', 'TAG F KW 3', 'BOOL NONPOLARIZABLE 1 False', 'False', 'INT SITESPERFRAGMENT 1 False', '3', 'DBL_ARRAY SITES 0 False']))
    core.set_local_option('PCM', 'PCMSOLVER_PARSED_FNAME', '@pcmsolver.10337.bbd42454')
    energy_scf, wfn = energy('scf', return_wfn=True)
    compare_values(nucenergy, NH3.nuclear_repulsion_energy(), 10, "Nuclear repulsion energy (PCM, total algorithm)")
--> compare_values(totalenergy, energy_scf, 10, "Total energy (PCM, total algorithm)")
    compare_values(polenergy, wfn.variable("PCM POLARIZATION ENERGY"), 6, "Polarization energy (PCM, total algorithm)")
    dimer = geometry("""
    H            0.000000000000    -0.759061990794     0.521953018286
    O            0.000000000000     0.000000000000    -0.065775570547
    H            0.000000000000     0.759061990794     0.521953018286

!----------------------------------------------------------------------------------!
!                                                                                  !
!         Total energy (PCM, total algorithm): computed value (-55.801284771319)   !
!     does not match (-55.802798786575) to atol=1e-10 by difference                !
!     (0.001514015256).                                                            !
!                                                                                  !
!----------------------------------------------------------------------------------!

Exit Status: overall (1)
<end of output>
Test time =   3.52 sec
----------------------------------------------------------
Test Failed.
"pcmsolver-ghost" end time: Mar 14 11:11 EDT
"pcmsolver-ghost" time elapsed: 00:00:03
----------------------------------------------------------
robertodr commented 1 year ago

Hei Lori! Thanks for pursuing the conda-forge packaging. I'd be happy to be included in the maintainers, though, as you wrote, the library is seeing only minimal maintenance these days. The differences you report do appear to be consistent with a change in radii, so it's OK to update the reference values in Psi4. Can you refresh my memory as to why the UFF radii should not be multiplied by 1.1?

As for which stable line to package: I'm in the (slow) process of finish up a new patch release of 1.3.x in #200 If you start the package to follow line 1.2.x, I think when 1.3.1 is out the bots will update it. So, I think it's fine to get 1.2.3 on conda-forge at first. Once that is in I'll mint 1.3.1 and Psi4 can work with the 1.3.x line.

loriab commented 1 year ago

Great, thanks!

The differences you report do appear to be consistent with a change in radii, so it's OK to update the reference values in Psi4. Can you refresh my memory as to why the UFF radii should not be multiplied by 1.1?

I really hadn't thought beyond "aha, this is the root of the discrepancy". Are you saying that the energies in psi4 shouldn't be changing that much and that instead there should be a compensatory 1.1 scaling somewhere in psi? Also, I noticed that 1.3.x introduced a new interface where one sets an empty object, then sets options by API. Psi is still using the old pcmsolver_new_v1112, so I didn't don't if that could also be a complication.

finish up a new patch release of 1.3.x

Good news on the 1.3.1 release. Note that I've got some cmake updates at https://github.com/PCMSolver/pcmsolver/compare/release/1.2...loriab:pcmsolver:v123_plus_ming (includes #193 ) that it'd be nice to get into the new release once they settle down.