New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Charge-aware frozen core #1350

Merged
merged 12 commits into from Jan 17, 2019

Conversation

Projects
None yet
9 participants
@PeterKraus
Copy link
Contributor

PeterKraus commented Nov 10, 2018

Description

This is my attempt to resolve #1271 . I would like to use this PR to hash out the FC functionality more thoroughly. The questions I have are:

  • what is the ideal expected behaviour of freeze_core True
  • do we want to implement "previous shells" freezing, using eg. freeze_core {1,2,...}
  • how about "ridiculous" cases such as Li2+ or F2-, how much hand-holding do we want the code to do and at which point we give up?

Todos

Notable points (developer or user-interest) that this PR has or will accomplish.

Checklist

Status

  • Ready for review
  • Ready for merge
@PeterKraus

This comment has been minimized.

Copy link
Contributor

PeterKraus commented Nov 10, 2018

To give my 2 cents to the questions above:

  • freeze_core True should be a best-effort solution, always freezing as much as possible, but ensuring there are valence electrons (ie. for Ca2+, [Ne] core is frozen instead of [Ar], while for Ca0, [Ar] is frozen; similarly for Ar- the whole [Ar] shell should be frozen). This one is charge aware. This option should never crash the code.
  • freeze_core N where N = [1,2,3...] will always freeze the N-th previous full shell, ie for N=1 the current behaviour in master. Validation error when no valence electrons remain, to avoid the non-helpful segfault in #1271.
  • @JonathonMisiewicz's idea of per-atom frozen cores (or in my view even better: per-fragment frozen orbitals) is a good one, and I'd be keen to get that sorted, but it'd probably need changes to molecule {} or somewhere else.
@hokru

This comment has been minimized.

Copy link
Contributor

hokru commented Nov 10, 2018

There is also the old-school way via orbital energy thresholds (below -3.5 eV works quite well). Could indirectly solve specific wants.

per atom/fragment declaration could be useful perhaps and would be, I think, a stand-out feature.

@JonathonMisiewicz

This comment has been minimized.

Copy link
Contributor

JonathonMisiewicz commented Nov 10, 2018

I agree with Peter on all points, so far. I have an unexpectedly long to-do list today, so I'll refrain from looking over the code until we're agreed on Peter's questions.

@PeterKraus

This comment has been minimized.

Copy link
Contributor

PeterKraus commented Dec 7, 2018

Just to prevent things falling in between the cracks, @loriab and I have discussed this issue on Slack a while ago, and the consensus was that while adding FC spec to Molecule might be attractive because it's easier to do, it really ought to belong to the Wavefunction object.

To paraphrase @loriab's comments, as the Wavefunction owns a Basisset which owns a Molecule anyway, it is the most logical place to plug this logic in. It should process options like frozen_docc Int or freeze_core True (or any possible freeze_core N spec where N is the N-th previous noble gas, or potential per-fragment frozen_docc). For the latter case (freeze_core True), it should take into account the fragment charge, or any electrons already frozen by ECP's, so that atoms such as Ca2+ with FC have some valence electrons.

@dgasmith dgasmith added this to the Psi4 1.3 milestone Jan 1, 2019

@dgasmith dgasmith requested review from jturney , andysim and robertodr Jan 3, 2019

@andysim

This comment has been minimized.

Copy link
Member

andysim commented Jan 3, 2019

I'm late to look at this one - sorry about that. I fully agree with the design ideas behind this. Yes, True should just err on the side of freezing questionable electrons (and, yes, it should make sure it leaves some!). Exiting gracefully when too many electrons have been frozen by the user is also important. I really like the idea of the fragment based approach, but I think it's too involved to be considered in the current PR.

@andysim

andysim approved these changes Jan 3, 2019

Copy link
Member

andysim left a comment

Great stuff - this will give much more sane behavior for corner cases. Sorry for being so slow to look at it. LGTM

@PeterKraus

This comment has been minimized.

Copy link
Contributor

PeterKraus commented Jan 3, 2019

I've just noticed that this is meant for the 1.3 release. I'll have another look at it and will implement both freeze_core True/False and freeze_core 0/-1/-2/-3 options, where False == 0 witout FC, True is the "smart handling" (previous rare gas of each atom, taking into account charge), and -1/-2/-3 would be "strict" N-th previous rare gas, ignoring charge but throwing a reasonable error when all electrons are frozen. I will have it done by Monday, if that's OK?

@andysim

This comment has been minimized.

Copy link
Member

andysim commented Jan 3, 2019

That would be excellent, if you have time. I think a little bit of documentation and testing would be needed in that case, to alert users of the new syntax.

@loriab

This comment has been minimized.

Copy link
Member

loriab commented Jan 3, 2019

Side note, as I haven’t looked at this yet today: I was going to add some qcelemental calls to enable lookup of levels of frozen core per element. Useful for this or not so much since c-side or other reasons?

@PeterKraus PeterKraus force-pushed the PeterKraus:fc_valence_fix branch from 810b787 to c110f9d Jan 4, 2019

@PeterKraus

This comment has been minimized.

Copy link
Contributor

PeterKraus commented Jan 4, 2019

OK, from my point of view this one is ready for a review. The old syntax (freeze_core true/false) is equivalent to freeze_core 1/0, with 1 being positive-charge-aware; it simply un-freezes another full shell in that case.

The new syntax of freeze_core -1/-2/-3 attempts to freeze Nth previous noble gas on each atom, and throws an exception if no valence electrons remain in the fragment (due to charge).

Show resolved Hide resolved psi4/src/read_options.cc Outdated
Show resolved Hide resolved tests/dfmp2-ecp/input.dat Outdated
// Freeze the number of core electrons strictly corresponding to the
// requested previous n-th shell.
for (int A = 0; A < mymol->natom(); A++) {
double Z = mymol->Z(A);

This comment has been minimized.

@susilehtola

susilehtola Jan 9, 2019

Contributor

Isn't this an int?

This comment has been minimized.

@PeterKraus

PeterKraus Jan 9, 2019

Contributor

Doesn't look like that to me:

/// Nuclear charge of atom
const double& Z(int atom) const;

This comment has been minimized.

@susilehtola

susilehtola Jan 9, 2019

Contributor

Huh... why is that

This comment has been minimized.

@loriab

loriab Jan 16, 2019

Member

Because in the beginning there was the thought that Z would be effective nuclear charge and float in case any MM atoms were around. At least that's my understanding.

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 9, 2019

@JonathonMisiewicz the test in question is dfmp2-ecp which has been unreliable for some time, way before #1456 was merged. See #1433.

@JonathonMisiewicz

This comment has been minimized.

Copy link
Contributor

JonathonMisiewicz commented Jan 9, 2019

@JonathonMisiewicz the test in question is dfmp2-ecp which has been unreliable for some time, way before #1456 was merged. See #1433.

I can't tell if we're disagreeing about facts or best practices.

Peter added new tests of SCF energies to dfmp2-ecp. A few days ago, his new tests passed. Before the force-push, his tests failed. This tipped us off that some other PR changed the energy to which frozen core ECP computations converge. Peter was able to reproduce the change in energies on a non-PR branch, confirming it.

The first problem is that unless I missed something, there shouldn't have been a PR that would change those energies. Apparently, there is one, and we don't know how many computations it afflicts. Probably just frozen core ECP, but there is a question mark on that one.

Second, one of these two answers has to be wrong. I cannot believe this is a case of computations being insufficiently converged or landing on different states. If Psi is giving or was giving wrong answers, then at the absolute minimum, we have an obligation to let people know. And if this test is just meant to detect changes in energies without worrying about why, that needs to be said in the test itself.

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 9, 2019

What I mean is that ECPs appear to have been broken for a while in Psi4. Even a month ago, before any of the SAD guess stuff was included, the dfmp2-ecp job was mysterious. Namely, it failed in the initial version of the new SAD guess, but mysteriously was working again with the reorganized but functionally equivalent version. However, at the same time, the dfmp2-ecp job fails to converge altogether when started from the CORE guess, which really makes no sense since it's dealing with a very weakly interacting noble gas dimer, which should converge straight away.

Are the failures/differences related to the presence of ECPs, or do they also occur in the case of no ECPs?

@PeterKraus

This comment has been minimized.

Copy link
Contributor

PeterKraus commented Jan 9, 2019

We can argue about whether the dfmp2-ecp test was or wasn't broken before PR #1350 was made, but the following point remains:

  • the original part of the test is not touched by this PR, and by looking at git blame for the input.dat file, hasn't been touched in quite some time.
  • as part of this PR, i have added a completely new section of tests, that are Kr--Kr interaction energies with DF-MP2/3-21G (no ECP's).
  • a previous version of this PR, from around 2019-01-04, passed both the new and the old test.
  • a later version of this PR, from until and excluding 2019-01-09, did not pass the new part of the test, but still passed the old test
  • the current version of this PR has modified the new part of the test, after rebasing from master; it passes both the old test and the new test to 7 DP (I need to tighten convergence on the energy to pass it to 8 DP).

A good diff to look at the issue is here:
9e6329a#diff-32d6050d16bb770a64f7f94a67f26101
which compares the output.ref of the test for the version that passed on 2019-01-04 when the test was originally written, and passes the test now (locally) on 2019-01-09 with updated values.

The puzzling thing is: what caused the change of the Kr--Kr interaction energies (and Kr MP2/3-21G energies converged to 1e-9)? What change would cause that, but not affect the first part of the test?

PeterKraus added some commits Jan 9, 2019

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 9, 2019

This is very odd, indeed. With Q-Chem I get the following results

HF MP2
Kr2 dimer -5478.3777659931 -5478.43678797
Kr atom -2739.1975667825 -2739.22713486
interaction energy 0.017368 0.017482

which are close to the original df values, so the new values are wrong.

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 9, 2019

Rerunning the Kr --- Kr 3-21G calculation in Psi4 master gives me 0.017482 with exact integrals, and 0.017535 with density fitting.

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 9, 2019

Btw since the 3-21G case doesn't even have ECPs, it should go in its own test.
Where you could check both the conventional SCF/MP2 as well as the density-fitted ones.

Split tests to dfmp2-ecp and dfmp2-fc. Replaced 3-21g with def2-svp i…
…n dfmp2-fc to get around fitting conditions change.
@PeterKraus

This comment has been minimized.

Copy link
Contributor

PeterKraus commented Jan 9, 2019

To summarise, it seems that running DF-MP2/3-21G is a bad idea for these small systems. The value of the interaction energy one gets, even with E_CONVERGENCE and D_CONVERGENCE at 1.0e-12, differs in the 4th DP based on which DF_FITTING_CONDITION is used. This caused the change in the energy during this PR (df_fitting_condition 1.0e-10 -> 1.0e-12)

The def2-svp basis set is not affected nearly as much, as its fitting basis is better suited for the purpose. Switching the basis set in the test and tightening the energy/density convergence criteria to 1e-10 passes agreement with 1e-12 reference to 8DP's.

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 9, 2019

Well, one wouldn't normally even run DF with basis sets as tiny as 3-21G. IIRC conventional integrals are faster since you get better screening..

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 9, 2019

I think it'd be good to also check freeze_core with conventional integrals. The results should be very close and they're both almost instantaneous calculations.

susilehtola added a commit to susilehtola/psi4 that referenced this pull request Jan 10, 2019

Restore previous DF_FITTING_CONDITION default
The value for ```DF_FITTING_CONDITION``` used by the code was changed in psi4#1446, which resulted in a change of numerical values. This PR restores the previously used threshold. See discussion in psi4#1350 and on Slack.

@susilehtola susilehtola referenced this pull request Jan 10, 2019

Merged

Restore previous DF_FITTING_CONDITION default #1472

3 of 6 tasks complete
@loriab

loriab approved these changes Jan 16, 2019

Copy link
Member

loriab left a comment

LGTM!

@CDSherrill is hunting for HFS's lore of frozen core, but that's a separate implementation.

number of frozen orbitals can be attained by using the keywords
|globals__num_frozen_docc| (gives the total number of orbitals to freeze,
program picks the lowest-energy orbitals) or |globals__frozen_docc| (gives
the number of orbitals to freeze per irreducible representation) -*/

This comment has been minimized.

@loriab

andysim added a commit that referenced this pull request Jan 16, 2019

Restore previous DF_FITTING_CONDITION default (#1472)
* Restore previous DF_FITTING_CONDITION default

The value for ```DF_FITTING_CONDITION``` used by the code was changed in #1446, which resulted in a change of numerical values. This PR restores the previously used threshold. See discussion in #1350 and on Slack.

* Found a few more places where DF_FITTING_CONDITION was not used.

@loriab loriab merged commit 336bee0 into psi4:master Jan 17, 2019

2 checks passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
psi4.psi4 #20190109.6 succeeded
Details

@loriab loriab referenced this pull request Jan 17, 2019

Merged

Finish SAD reorganization and make it the default guess #1458

8 of 10 tasks complete
@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 18, 2019

Doesn't seem to be working too great; I get on my machine

$ psi4 -n 8 
	SAPT0 energy with Ca2+, without ECP: computed value (-0.095615678) does not match (-0.095616042) to 8 digits.
Traceback (most recent call last):
  File "/home/work/psi4/install.susi/bin/psi4", line 287, in <module>
    exec(content)
  File "<string>", line 38, in <module>
  File "/home/work/psi4/install.susi/lib/psi4/driver/p4util/util.py", line 230, in compare_values
    raise TestComparisonError(message)

TestComparisonError: 	SAPT0 energy with Ca2+, without ECP: computed value (-0.095615678) does not match (-0.095616042) to 8 digits.

and the result is the same regardless of the guess...

@loriab

This comment has been minimized.

Copy link
Member

loriab commented Jan 18, 2019

Thanks, your fix does the trick. I also get the above, but I'm guessing that the reference values just weren't constructed defensively. Varying convergence now ...

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 18, 2019

I bet this is since #1472 was merged.

Is there no non-df version of SAPT?

@susilehtola

This comment has been minimized.

Copy link
Contributor

susilehtola commented Jan 18, 2019

So, the reference value should just be regenerated with the correct threshold.

And/or switch to conventional integrals.

@loriab

This comment has been minimized.

Copy link
Member

loriab commented Jan 18, 2019

I should know this. My best recollection is that there is SAPT with conventional integrals but there might be terms that are always constructed through DF. Unless I'm mixing that up with Molpro's CCSD-F12.

Yes, ref energies were probably constructed with default ~6 convergence but checking to 8.

@PeterKraus

This comment has been minimized.

Copy link
Contributor

PeterKraus commented Jan 18, 2019

Apologies. I could swear I ran all references at 1.0e-12, but clearly I must have messed them up for sapt-ecp.

@loriab

This comment has been minimized.

Copy link
Member

loriab commented Jan 18, 2019

No worries. Thanks for checking up on the refs. As I recall, you ran a couple dozen refs to examine the fitting metric issue.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment