looser symmetry tolerance for optimization

[loriab]

The List 21 tied with psi4/psi4#232

[loriab]

Alternatively, adjusting this https://github.com/psi4/psi4private/blob/master/src/bin/optking/set_params.cc#L422 Opt_params.symm_tol parameter may be all that's needed. Having a test case where this goes wrong in WebMO would help in finding an acceptable approach.

[polikw]

I will look into this and find an example.

[CDSherrill]

If I recall correctly, methane might be one such?

David

On Mon, Feb 8, 2016 at 3:36 PM, polikw [email protected] wrote:

I will look into this and find an example.

—
Reply to this email directly or view it on GitHub
#4 (comment).

[polikw]

I can reproduce this error in WebMO for
molecule: H2O with idealized comprehensive clean-up
calculation: Optimize and Frequencies
basis: STO-3G, 3-21G, or 6-31G
I also get it for
molecule: H2CO with idealized comprehensive clean-up
calculation: Optimize and Frequencies
basis: 3-21G

Two points:

1. This error occurs for me when optimize is followed by frequencies, through not for optimize alone, suggesting a possible mismatch between tolerances in these two modules.
2. Interestingly, H2CO runs fine with STO-3G or 6-31G basis sets, though not with 3-21G

I am attaching the archive of my jobs, which you can import from Job Manager: New Job: Import Job : Import Type: WebMO archive. From this you can inspect the input and output files in detail. IMPORTANT: You must FIRST unzip this *.tar.zip file to to a *.tar file, and THEN import the *.tar file into WebMO (because GitHub apparently does not support attaching tar files!).

archive21.tar.zip

[polikw]

Seeing David's comment, I just verified that

• molecule: CH4 with idealized comprehensive clean-up
• calculation: Optimize and Frequencies
• basis: STO-3G, 3-21G, or 6-31G
  also all fail.

For the details, convert the attached *.tar.zip to a *.tar file and import into WebMO.
archive27.tar.zip

[loriab]

Ok, I think I have the fix.

original input:

memory 1024 mb

molecule {
0 1
C 0.00000000 0.00000000 0.00000000
H 1.09000000 0.00000000 0.00000000
H -0.36333333 0.83908239 0.59332085
H -0.36333333 0.09428973 -1.02332709
H -0.36333333 -0.93337212 0.43000624
}

set globals {
scf_type pk
basis STO-3G

reference rhf
}

optimize('scf')
frequencies('scf')
oe = OEProp()
oe.add("MULLIKEN_CHARGES")
oe.compute()
print_variables()

working input:

memory 1024 mb

molecule mol {
0 1
C 0.00000000 0.00000000 0.00000000
H 1.09000000 0.00000000 0.00000000
H -0.36333333 0.83908239 0.59332085
H -0.36333333 0.09428973 -1.02332709
H -0.36333333 -0.93337212 0.43000624
}

set globals {
scf_type pk
basis STO-3G

reference rhf
}

mol.update_geometry()
mol.symmetrize(1e-3)

optimize('scf')
e, wfn = frequencies('scf', return_wfn=True)
oe = OEProp(wfn)
oe.add("MULLIKEN_CHARGES")
oe.compute()
print_variables()

The main change are the mol.update_geometry() and mol.symmetrize(1e-3) lines. The latter pre-symmetrizes the input structure (from Cs to C2v) so that when optking symmetrizes the molecule for the second pass and gets C2v, this doesn't look like a change in point group. Note that you need a handle on the molecule in order to do these operations, so have to declare with molecule anyname { rather than just molecule {. So that fixes the opt and freq. In an unrelated issue, the OEProp part fails. The solution to this is just requesting return_wfn=True, collecting e, wfn, and passing OEProp(wfn) the wavefunction object.

So there's no actual changes to psi4. @bwb314, would you want to try incorporating these changes into the WebMO templates?

[CDSherrill]

Can you check how non-symmetric the geometry was that WebMO was providing
in some of the test cases? I thought it was very close to symmetric. If
we fail for things that close, we might still want some change to Psi4
proper, either a changed tolerance or an automatic symmetrization with some
tight value. I am concerned that some user might input s perfectly good
geometry and still get these kinds of WebMO-induced symmetry change
failures because some internal Psi4 tolerance is set to 1E-16 or something.

On Thursday, April 14, 2016, Lori A. Burns [email protected] wrote:

Ok, I think I have the fix.

original input:

memory 1024 mb

molecule {
0 1
C 0.00000000 0.00000000 0.00000000
H 1.09000000 0.00000000 0.00000000
H -0.36333333 0.83908239 0.59332085
H -0.36333333 0.09428973 -1.02332709
H -0.36333333 -0.93337212 0.43000624
}

set globals {
scf_type pk
basis STO-3G

reference rhf
}

optimize('scf')
frequencies('scf')
oe = OEProp()
oe.add("MULLIKEN_CHARGES")
oe.compute()
print_variables()

working input:

memory 1024 mb

molecule mol {
0 1
C 0.00000000 0.00000000 0.00000000
H 1.09000000 0.00000000 0.00000000
H -0.36333333 0.83908239 0.59332085
H -0.36333333 0.09428973 -1.02332709
H -0.36333333 -0.93337212 0.43000624
}

set globals {
scf_type pk
basis STO-3G

reference rhf
}

mol.update_geometry()
mol.symmetrize(1e-3)

optimize('scf')
e, wfn = frequencies('scf', return_wfn=True)
oe = OEProp(wfn)
oe.add("MULLIKEN_CHARGES")
oe.compute()
print_variables()

The main change are the mol.update_geometry() and mol.symmetrize(1e-3)
lines. The latter pre-symmetrizes the input structure (from Cs to C2v) so
that when optking symmetrizes the molecule for the second pass and gets
C2v, this doesn't look like a change in point group. Note that you need a
handle on the molecule in order to do these operations, so have to declare
with molecule anyname { rather than just molecule {. So that fixes the
opt and freq. In an unrelated issue, the OEProp part fails. The solution to
this is just requesting return_wfn=True, collecting e, wfn, and passing
OEProp(wfn) the wavefunction object.

So there's no actual changes to psi4. @bwb314 https://github.com/bwb314,
would you want to try incorporating these changes into the WebMO templates?

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#4 (comment)

[loriab]

No 1e-16. The general Molecule uses #define DEFAULT_SYM_TOL 1.0E-8, while optking resymmetrizes to 0.01.

That said, I don't really understand what's going on with the methane. By Avogadro, it looks perfectly symmetric (1–2 decimal places). But the methane needs symmetrizing with tol<=0.66 and printing out the atom distances I think gives the below as distinct atom distances.

start:   [0.0, 1.941999, 2.059801, 2.378454, 3.363642, 3.883999, 4.119603]
opt end: [0.0, 2.046667, 2.363288, 3.342193, 4.093334]

[CDSherrill]

This is puzzling. Suggest looking back at H2O case -- I think that one
looked really darn close to symmetric.

On Thu, Apr 14, 2016 at 9:32 PM, Lori A. Burns [email protected]
wrote:

No 1e-16. The general Molecule uses #define DEFAULT_SYM_TOL 1.0E-8, while
optking resymmetrizes to 0.01.

That said, I don't really understand what's going on with the methane. By
Avogadro, it looks perfectly symmetric (1–2 decimal places). But the
methane needs symmetrizing with tol<=0.66 and printing out the atom
distances I think gives the below as distinct atom distances.

start: [0.0, 1.941999, 2.059801, 2.378454, 3.363642, 3.883999, 4.119603]
opt end: [0.0, 2.046667, 2.363288, 3.342193, 4.093334]

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#4 (comment)

[loriab]

By calculator, it's symmetric. But with only writing 8 decimal places in the input, it mayn't be possible to satisfy the default symmetry detection of 1e-8. The calibration of tol in symmetrize(tol), I'm not quite following, but in practice, it's behaving exactly as it should. Add a symmetrize(0.01) to any structure that originates with the WebMO molecule builder, and all should be good to go. This is the same thing that we do for structures downloaded from PubChem. Attached is a healed water output.
polik2.out.gz

[polikw]

Please forgive the naivety of my comments, since I do not have a detailed understanding of the inner workings of psi4. Perhaps they are useful, through maybe not...

It seems to me that two different parts of the program have different symmetry detection tolerances, first when the symmetry is initially determined, and second with the optimizer so as to be more efficient. Making these tolerances the same would seem to be a good start.

However, I would not be surprised if these different sections of code even used different algorithms for symmetry detection, since they were presumably written by different people. Thus it might make sense either to make the second tolerance slightly looser than the first, or to symmetrize the molecule to a very high tolerance before starting the process (which I think is Lori's solution).

Thank you for working on this issue!

1. Does the optimizer and the symmetry detection routine

[bwb314]

@loriab I just checked your fix in. All of the test cases provided in this thread can now be optimized successfully using the WebMO interface.

[loriab]

I think it is actually behaving rationally, though there's several things going on. In any psi4 run, the input geometry symmetry is taken as given, subject to a 1e-8 tolerance. After an optimization pass, the projected geometry is loosely symmetrized to 0.01 to clean up any numerical error. At this point if the point group changes to a higher symmetry, this is regarded as a proper error and the message you've seen gets thrown.

• The inputs in these WebMO test cases, we'll call "infinitesimally asymmetric".
• If you run a single point or frequency calculation on an infinitesimally asymmetric geom, it runs fine (and uses the lower point group).
• If you exaggerate the asymmetry slightly and demand very tight geometry convergence, as in the input below, it will still detect as the lower point group and proceed for many optimization cycles. But eventually it reaches a point where the optimizer-projected geometry truly detects as the higher point group (to 1e-8) and the "Point group changed!" error throws.
• The WebMO problem is simply that the infinitesimally asymmetric geom inspires the optimizer's first projected geometry to be at the higher point group, so the situation in the previous bullet happens immediately.
• So either pre-symmetrize, as @bwb314 has already implemented. Or possibly writing out another decimal or two to the input files would allow the initial symmetry detection to occur as we had all expected.
• Oh, and the symmetry detection machinery is the same, thank goodness.

>>> cat polik13.in 
memory 1024 mb

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H 1.05000000 0.00000000 0.00000000
#H -0.35000000 0.00000000 0.98994949
H -0.351000000 0.00000000 0.98994949
}

set globals {
scf_type pk
basis STO-3G

reference rhf
MAX_ENERGY_G_CONVERGENCE 10
RMS_FORCE_G_CONVERGENCE 9
RMS_DISP_G_CONVERGENCE 8

}

mol.update_geometry()
print mol.schoenflies_symbol()
mol.print_out()
#mol.symmetrize(0.01)
print mol.schoenflies_symbol()
mol.print_out()

optimize('scf')
e, wfn = frequencies('scf', return_wfn=True)
oe = OEProp(wfn)
oe.add("MULLIKEN_CHARGES")
oe.compute()
print_variables()

#frequency('scf')

>>>  ../../objdir-pn/bin/psi4 -l ../../share/ polik13.in 
cs
cs

An error has occurred. Traceback:
<class 'p4xcpt.ValidationError'>: Point group changed! You should restart using the last geometry in the output, after carefully making sure all symmetry-dependent input, such as DOCC, is correct.:   File "<string>", line 78, in <module>

  File "../../share//python/driver.py", line 1112, in optimize
    raise ValidationError("""Point group changed! You should restart """


>>> grep '~' polik13.out 
  --------------------------------------------------------------------------------------------- ~
   Step     Total Energy     Delta E     MAX Force     RMS Force      MAX Disp      RMS Disp    ~
  --------------------------------------------------------------------------------------------- ~
    Convergence Criteria    1.00e-10 *    3.00e-04 *    1.00e-09 *    1.20e-03 *    1.00e-08 *  ~
  --------------------------------------------------------------------------------------------- ~
      1     -74.95354732   -7.50e+01      6.48e-02      6.06e-02      2.68e-01      1.97e-01    ~
      2     -74.96321368   -9.67e-03      3.49e-02      3.01e-02      1.17e-01      7.47e-02    ~
      3     -74.96585509   -2.64e-03      4.68e-03      3.60e-03      1.83e-02      1.11e-02    ~
      4     -74.96589941   -4.43e-05      1.02e-03      6.19e-04      3.76e-03      2.18e-03    ~
      5     -74.96590111   -1.70e-06      2.08e-04 *    1.71e-04      4.16e-04 *    4.02e-04    ~
      6     -74.96590119   -8.14e-08      3.05e-05 *    2.50e-05      6.59e-05 *    5.57e-05    ~
      7     -74.96590119   -1.94e-09      7.85e-07 *    6.36e-07      1.99e-06 *    1.60e-06    ~
      8     -74.96590119   -3.84e-13 *    5.68e-07 *    4.56e-07      1.43e-06 *    1.15e-06    ~
      9     -74.96590119   -2.70e-13 *    4.10e-07 *    3.29e-07      1.03e-06 *    8.24e-07    ~
     10     -74.96590119   -2.84e-14 *    2.94e-07 *    2.36e-07      7.42e-07 *    5.93e-07    ~
     11     -74.96590119   -2.84e-14 *    2.17e-07 *    1.74e-07      5.47e-07 *    4.37e-07    ~
     12     -74.96590119    2.84e-14 *    1.64e-07 *    1.31e-07      4.13e-07 *    3.30e-07    ~
     13     -74.96590119   -1.14e-13 *    1.23e-07 *    9.88e-08      3.11e-07 *    2.48e-07    ~
     14     -74.96590119   -9.95e-14 *    9.28e-08 *    7.44e-08      2.34e-07 *    1.87e-07    ~
     15     -74.96590119    9.95e-14 *    6.99e-08 *    5.60e-08      1.76e-07 *    1.40e-07    ~
     16     -74.96590119    2.84e-14 *    5.25e-08 *    4.21e-08      1.32e-07 *    1.06e-07    ~
     17     -74.96590119   -1.42e-14 *    3.95e-08 *    3.17e-08      9.96e-08 *    7.95e-08    ~
     18     -74.96590119   -2.84e-14 *    2.97e-08 *    2.38e-08      7.49e-08 *    5.98e-08    ~
     19     -74.96590119    4.26e-14 *    2.24e-08 *    1.79e-08      5.64e-08 *    4.50e-08    ~
     20     -74.96590119   -7.11e-14 *    1.68e-08 *    1.35e-08      4.24e-08 *    3.39e-08    ~
     21     -74.96590119   -7.11e-14 *    1.27e-08 *    1.02e-08      3.19e-08 *    2.55e-08    ~
     22     -74.96590119    4.26e-14 *    9.54e-09 *    7.64e-09      2.40e-08 *    1.92e-08    ~
     23     -74.96590119    4.26e-14 *    7.18e-09 *    5.75e-09      1.81e-08 *    1.44e-08    ~
     24     -74.96590119   -1.42e-14 *    5.40e-09 *    4.33e-09      1.36e-08 *    1.09e-08    ~
     25     -74.96590119    0.00e+00 *    4.06e-09 *    3.25e-09      1.02e-08 *    8.17e-09 *  ~
     26     -74.96590119   -1.42e-14 *    3.06e-09 *    2.45e-09      7.70e-09 *    6.15e-09 *  ~
     27     -74.96590119   -2.84e-14 *    2.30e-09 *    1.84e-09      5.79e-09 *    4.62e-09 *  ~