Skip to content
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

Nuclear Gradients #124

Open
wants to merge 14 commits into
base: master
Choose a base branch
from
Open

Nuclear Gradients #124

wants to merge 14 commits into from

Conversation

maxscheurer
Copy link
Member

@maxscheurer maxscheurer commented Jan 13, 2021

Implemented together with @Drrehn.

Features

  • MP2
  • ADC(0)
  • ADC(1)
  • ADC(2)
  • ADC(2)-x
  • ADC(3)
  • CVS-ADC(1) (by @iubr)
  • CVS-ADC(2) (by @iubr)
  • CVS-ADC(2)-x (by @iubr)

ToDo

  • consolidate some of the CVS equations for different methods
  • Code documentation (docstrings etc.)
  • Warning about performance
  • Basic documentation (calculations.rst)
  • rebase on Transition dipole moments for asymmetric operators #158
  • finite difference data test generation (preliminary version working)
  • tests:
    • individual components (densities etc.)
    • reference tests (preliminary version working)
  • rebase once Updates for solvent models (new user interface, new schemes) #123 is in (I've added some stuff from this branch here for testing purposes, will be removed)
  • add equation numbers from papers (Rehn 2019, Levchenko 2005)

@lgtm-com
Copy link

lgtm-com bot commented Jan 13, 2021

This pull request introduces 1 alert when merging 6b92a09 into 6e7ac34 - view on LGTM.com

new alerts:

  • 1 for Unused local variable

@lgtm-com
Copy link

lgtm-com bot commented Jan 13, 2021

This pull request introduces 1 alert when merging 7136d03 into 6e7ac34 - view on LGTM.com

new alerts:

  • 1 for Unused local variable

adcc/AmplitudeVector.py Outdated Show resolved Hide resolved
@mfherbst
Copy link
Member

Awesome 👍

@lgtm-com
Copy link

lgtm-com bot commented Jan 16, 2021

This pull request introduces 1 alert when merging 844cc1c into 6e7ac34 - view on LGTM.com

new alerts:

  • 1 for Unused local variable

@lgtm-com
Copy link

lgtm-com bot commented Jan 16, 2021

This pull request introduces 1 alert when merging 10980bf into 6e7ac34 - view on LGTM.com

new alerts:

  • 1 for Unused local variable

@lgtm-com
Copy link

lgtm-com bot commented Jan 17, 2021

This pull request introduces 1 alert when merging 3f82618 into 6e7ac34 - view on LGTM.com

new alerts:

  • 1 for Unused local variable

@lgtm-com
Copy link

lgtm-com bot commented Jan 19, 2021

This pull request introduces 1 alert when merging b92ee1b into 041cdcb - view on LGTM.com

new alerts:

  • 1 for Unused local variable

@lgtm-com
Copy link

lgtm-com bot commented Feb 21, 2021

This pull request introduces 1 alert when merging 1a3c802 into 041cdcb - view on LGTM.com

new alerts:

  • 1 for Unused local variable

@maxscheurer maxscheurer force-pushed the gradients branch 2 times, most recently from 1717c34 to 54a9181 Compare June 3, 2021 12:54
@lgtm-com
Copy link

lgtm-com bot commented Jun 3, 2021

This pull request introduces 1 alert when merging 54a9181 into e4ca2ed - view on LGTM.com

new alerts:

  • 1 for Unused local variable

@maxscheurer maxscheurer force-pushed the gradients branch 2 times, most recently from 6e67a9d to 365cd5a Compare June 3, 2021 13:57
@lgtm-com
Copy link

lgtm-com bot commented Dec 14, 2021

This pull request introduces 4 alerts when merging 92e7535 into e4ca2ed - view on LGTM.com

new alerts:

  • 3 for Unused import
  • 1 for Unused local variable

@lgtm-com
Copy link

lgtm-com bot commented Dec 14, 2021

This pull request introduces 4 alerts when merging df675fa into e4ca2ed - view on LGTM.com

new alerts:

  • 3 for Unused import
  • 1 for Unused local variable

@mfherbst
Copy link
Member

Of course tests need to work and all, but I wonder if a good first aim is to just get the things already implemented merged into master and release a new version?

@maxscheurer
Copy link
Member Author

The tests are horrible to implement... ensuring highly converged eigenvectors is difficult 😞 but in principle all the implemented gradient methods are correct 👍🏻

@iubr
Copy link

iubr commented Dec 15, 2021

Not sure if this is a good suggestion, but maybe we could run a couple of small molecules and save the results in a text file for pytest to compare agains (this is more or less how it's done in veloxchem). I can try out the cvs-adc gradients and find an optimal/good enough setup (I guess the question is of comparing analytical vs. numerical gradients)

@maxscheurer
Copy link
Member Author

The generation of finite-difference reference data works fine already, for both generic and CVS-ADC 😊👍🏻

@mfherbst
Copy link
Member

The tests are horrible to implement... ensuring highly converged eigenvectors is difficult

Can't you lower to tolerance or introduce perturbations in the molecule / orbitals to ensure they are sufficiently distinct (to avoid rotations in invariant subspaces etc.)

@maxscheurer
Copy link
Member Author

You mean the convergence tolerance/atol for asserts?

You mean geometry perturbations to break symmetries? Yes, I can try that for sure 👍🏻

@mfherbst
Copy link
Member

You mean the convergence tolerance/atol for asserts?

Yes. It's not amazing, but it will at least catch forgotten factors of two or switched signs.

@maxscheurer
Copy link
Member Author

@iubr, the CVS-ADC2-x don't seem to be correct right now, at least I get an error for H2O of ~1e-5 (compared to 5-point finite differences). Maybe I made some mistake while cleaning up the code... Which accuracy could you obtain in your tests? I'll try to run your original branch again at some point to get to the bottom of this.

@iubr
Copy link

iubr commented Dec 18, 2021

The ~1e-5 is similar to what I get with my code (in comparison to the 2-point finite difference numerical gradient). I checked it today with a 5-point finite difference and it is a little smaller, but still same order (1e-5). So if there is a bug, it would be in my module as well.
I had discussed this difference with Patrick and we had concluded that it could be due to the tolerance/convergence threshold for the adc step (I couldn't go better than 1e-6 because of the requirement for the HF reference). Of course, this might not be the right explanation for the difference. It's, of course, possible that there is still an undiscovered bug or mistake somewhere in the code/equations. I'll make another check next week. Let's see.

@maxscheurer
Copy link
Member Author

I had discussed this difference with Patrick and we had concluded that it could be due to the tolerance/convergence threshold for the adc step

I converge the SCF to 1e-12, and the ADC procedure is converged to 1e-10 in the energy. Given the fact that all other gradient methods so far can achieve at least 1e-7 accuracy, I'm convinced that one can achieve the same for CVS-ADC2-x . I'll try to add canonical ADC2-x gradients now and see how accurate one can get.

@iubr
Copy link

iubr commented Dec 19, 2021

Thanks a lot!

@iubr
Copy link

iubr commented Feb 11, 2022

Hej Max,

Sorry for the long silence. There have been a couple of things that I had to prioritize, so I couldn't focus on cvs-adc2x, but I have been debugging these last two weeks. I didn't find anything yet, but still have some ideas of what to try (I confirmed the excitation energy PDMs are correct and tested the lambda multipliers against numerical dipole moments -- the error is the same as in cvs-adc2). I am now re-checking all the omega multipliers, but I am not sure how much longer it will take me to get to the bottom of it. I hope not too long. Let me see if I can do it over the weekend, otherwise please go ahead without cvs-adc2x. What was the error in the correlation energy TPDMs?

And thanks a lot for all the help! Super-exciting that you managed to add adc2x and adc3!!

/Iulia

@maxscheurer
Copy link
Member Author

maxscheurer commented Feb 11, 2022

I fixed the TPDMs for CVS-ADC(2)-x in this commit, where the prefactor of g2a.vvvv was wrong.

The "correlation part" of the energy is tested here.

@iubr
Copy link

iubr commented Feb 11, 2022

Thanks! I'll look into it.

@iubr
Copy link

iubr commented Feb 13, 2022

I double checked the equations for the TPDMs and and I do get a factor 2 in the VVVV block for cvs-adc2x. In my module, it's the only way I get the correct excitation energy. I also wrote a stand-alone code to compute the excitation energy from the density matrices and without this factor 2 I do not get the correct energy. I attach the code here as a txt file in case you'd like to have a look cvs_adc2x_dms_from_scratch.txt
. I still have to have a more careful look at the adcc version to understand what's going on, but next week is quite crazy for me up until Thursday. It could also be great to have a short meeting -- maybe you spot something I don't in the equations. What do you think?

Regarding cvs-adc2x updates, over the weekend I still haven't managed to figure out the reason why analytical vs. numerical agree only up to 1e-5 and not 1e-7 (as for the other adc orders). Maybe it's best you go ahead without cvs-adc2x for the moment and we add it once this mystery is solved (which I hope will not be that much longer).

@maxscheurer
Copy link
Member Author

Thanks for the script, I'll have a look.

@iubr
Copy link

iubr commented Feb 13, 2022

Thanks! If you'd like an overview of the density matrices, they are from https://doi.org/10.1063/5.0058221, Appendix C, Table II, 3rd column (ADC matrix terms).

@iubr
Copy link

iubr commented Feb 21, 2022

Hi Max,

I just found the bug! It was a sign in one of the terms for the right-hand-side of the CV block, lambda multipliers. See here.
I tried to merge the latest version of adcc/gradients branch into my cvs_gradients branch, but for some reason I cannot get it running on my computer, so I fixed the bug using the old version of the code. If you could make the fix at your end and test cvs-adc2x again, it would be super great!! Thanks a lot!! I'll make a few more tests on my end just to make sure there's nothing else.

/Iulia

@maxscheurer
Copy link
Member Author

That's great news! I will try it out later. 👍🏻😊

@iubr
Copy link

iubr commented Feb 21, 2022

Thank you!! Let me know how it goes :D
(I still haven't figured out the puzzle of the factor 2 in the VVVV block between my module and adcc, but I'm planning to look into it again tomorrow)

@maxscheurer
Copy link
Member Author

I found the reason for the weird factor 2 behavior in this
block of code. The blocks oooo and vvvv were added again in the case of CVS, such that they existed twice in the TwoParticleDensityMatrix object. 😁

I found some other inconsistencies in the formation of g2a_hf and g2a_oresp, with these fixes altogether, everything works!!! 💪

@maxscheurer maxscheurer force-pushed the gradients branch 2 times, most recently from 06ee1b2 to 51115ed Compare February 21, 2022 15:10
@mfherbst
Copy link
Member

Cool! Great work.

@iubr
Copy link

iubr commented Feb 21, 2022

Wonderful ^-^!! Thanks a lot @maxscheurer!!

Copy link
Member

@mfherbst mfherbst left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a few comments ... obviously I have not at all checked the equations.

Comment on lines +269 to +270
# with pytest.raises(InputError):
# obtain_guesses_by_inspection(matrix1, n_guesses=40, kind="any")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could change this to a test for the warning.

Comment on lines +53 to +58
molecules = [
Molecule("h2o", 0, 1),
]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
molecules = [
Molecule("h2o", 0, 1),
]
molecules = [Molecule("h2o")]

adcc/testdata/dump_fdiff_gradient.py Show resolved Hide resolved
g2_ao_1, g2_ao_2: relaxed two-particle density matrices
"""
natoms = self.mol.natom()
Gradient = {}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why capitalised?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shamelessly copied from some psi4numpy tutorial, will change this 😄

# Build Total OEI Gradient
Gradient["OEI"] = Gradient["T"] + Gradient["V"] + Gradient["S"]

# Build TE contributions
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TEI contributions?

@dataclass(frozen=True)
class GradientComponents:
natoms: int
nuc: np.ndarray
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nuc is a bit short, maybe add a comment here to explain

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good idea, I think it would be good to make all attribute names a bit more verbose.

nuc: np.ndarray
overlap: np.ndarray
hcore: np.ndarray
two_electron: np.ndarray
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

comment: Two electron integrals?

adcc/gradients/__init__.py Show resolved Hide resolved
return self.excitation_or_mp.reference_state

@property
def _energy(self):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why the underscore? If someone uses it it's sort of their fault, no? If you add a warning I think that's fine.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point. Will rename it to correlation_energy or so to make it more clear what exactly is returned.

Comment on lines +337 to +348
lam = conjugate_gradient(A, rhs=rhs, x0=x0, Pinv=Pinv,
explicit_symmetrisation=None,
callback=default_print)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Which tolerance is used and why? Should for sure be related to the convergence used for the ADC excitation vectors ... or maybe aim for relative residual decrease?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it should rather be related to the SCF convergence criterion, because we compute the response of the HF ground state (kind of 😄). In addition, MP states don't have an additional convergence criterion, so I will make this depend on the SCF conv tol. 👍

maxscheurer and others added 14 commits May 6, 2023 10:03
fix dot product issue for bare tensors

add working gradient implementation

add gradient provider to pyscf

flake

fix excitation view tests
lgtm and coveralls (hopefully)

more fdiff

fix
hijacked AmplitudeVector for CVS-ADC orbital response.

CVS-ADC0 working!

cvs-adc1 working
added cvs-adc2 2PDMs,started orbital response (need to check if the RHS is correct)

cvs-adc2 gradient working.

cvs-adc2x working
… vs. cvs-adc rhs.

fixed 1/sqrt(2)

additional comments, RHS

switched to using symmetrise and antisymmetrise to construct density matrices
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants