Skip to content

Davidson works when doubles folding in ADC(2)#176

Open
TianLin-228 wants to merge 7 commits intoadc-connect:masterfrom
TianLin-228:davidson_foldingMatrix
Open

Davidson works when doubles folding in ADC(2)#176
TianLin-228 wants to merge 7 commits intoadc-connect:masterfrom
TianLin-228:davidson_foldingMatrix

Conversation

@TianLin-228
Copy link

@TianLin-228 TianLin-228 commented Aug 16, 2024

Implement a modified Davidson procedure for doubles-folding. It solves pseudo-eigenvalue problem ωi ui = A(ωi) ui, where A(ωi) is an effective ADC matrix which depends on the eigenvalue ωi. The modified procedure is based on the implementation stated in the paper ( J. Chem. Phys. 158, 124121 (2023)).

The algorithm consists of state-specific macro and micro iterations and further DIIS acceleration. Micro-iterations use Davidson iterations to diagonalize the effective matrix Aeffi) with the fixed eigenvalue ωi, giving a new eigenvalue ω'i. After convergence, a new macro-iteration begins with the effective matrix Aeff'i) . This process continues until the eigenvectors meet a threshold (default=1e-3), after which DIIS further converges the individual roots.

NB: This procedure only works for single excitations.

Changes:

  1. In file “davidson.py”:
  • New functions named “davidson_folded_DIIS“ and “eigsh_folded” are added, which implements the main procedure of the modified Davidson method. It consists of loops for macro iterations, in which the davidson_iterations are performed as micro iterations, and loops for DIIS iterations.

  • Changes in function “davidson_iterations” :

    • For the folding case, when it breaks at convergence or otherwise, the eigenvectors with the number of n_block are calculated by the rvecs and guesses, and the subspace_vectors is updated for the next macro iteration;
    • A different convergence criteria for the folding case as following:
      ω0i - ω1i > ωn-1i - ωni , where ω0i is the initial guess of ωi and ωni is the Ritz value of the current micro iteration.
  • Add new Class DIIS to implement diis acceleration.

  • There are two new Classes: Class FoldedDavidsonState inherits from Class DavidsonState and contains more properties about the DavidsonState in doubles-folding case; Class DIIS is used when performing extrapolation on the corrected vector.

  • Add new properties for Class DavidsonState: “folded” is the flag to distinguish normal or modified davidson procedures when running “davidson_iterations”, and “DIIS_iter” sums the total number of performed DIIS iterations for all required states.

  • In the function “default_print”, there are new identifiers for the folding case.

  1. In the file “AdcMatrix.py”, the matrix-vector product of the folded ADC matrix with singles component is defined in the new Class AdcMatrixFolded.

  2. Workflow - two options of guesses vector in doubles-folding case:
    guesses_fold = “adc1”, which takes adc1 results (excitation_vector and excitation_energy_uncorrected) as initial guess vectors and guess eigenvalues.
    Or obtaining guess vectors by inspecting the diagonal matrix elements as default.

Testing:

test_davidson_folded.py: testing davidson_folded_DIIS in two cases of given guesses and taking guesses as adc1 results.
test_workflow.py: I added two tests for new parameters “fold”, “guesses_fold” and “omegas”.

@TianLin-228 TianLin-228 marked this pull request as ready for review August 19, 2024 13:06
@apapapostolou
Copy link
Contributor

@TianLin-228 could you update this against the current master?

@TianLin-228
Copy link
Author

Updating brings the branch up to date with the current master. The updates include:

  1. davidson.py: involving the use of n_block

  2. davidson_test.py: Expanded testing for davidson_folded_DIIS by adding two cases:

     - Testing with  given guess inputs.
    
     - Testing by using ADC1 results as the guess.
    

"""
diag = super().diagonal().pphh
e = diag.ones_like()
u2 = self.block_apply("pphh_ph", u1.ph) / (e * self.omega - diag)
Copy link
Contributor

Choose a reason for hiding this comment

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

You should explicitly symmetrise the doubles vector. This can be done, for example, using the IndexSymmetrisation class from solver/explicit_symmetrisation.py. You can find usage examples in respondo/MatrixWrapper.py.

print(" Total solver time: ", strtime(soltime), file=file)
elif identifier == "restart":
print("=== Restart ===", file=file)
# For folded matrix
Copy link
Contributor

Choose a reason for hiding this comment

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

The printout for the folded davidson solver looks a bit messy, the alignment is off, for example. I also noticed that the style is different compared to the standard davidson solver in adcc. Was there a specific reason for changing the formatting? If not, it might be better to stick with the earlier style for better readability.

assert res.converged
assert res.eigenvalues[:n_states] == pytest.approx(ref_triplets[:n_states])

def test_adc2_singlets_folded(self):
Copy link
Contributor

Choose a reason for hiding this comment

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

I would put it in an other class. Like TestSolverDavidsonFolded.
You should also think about testing the correct symmetry of the unfolded eigenvector.

@apapapostolou
Copy link
Contributor

@TianLin-228 what is the current status here? Can you push your local changes?

@TianLin-228
Copy link
Author

@TianLin-228 what is the current status here? Can you push your local changes?

Hi, sorry for delayed update!

I've been suffering from the symmetry test, and there is another one triggered by max_subspace. I’d also like to discuss the tests for unfold in AdcMatrix_test.py. This is also to keep a record.

block_orders=matrix.block_orders,
intermediates=matrix.intermediates)
if matrix.method.name != "adc2":
raise TypeError("The matrix should be ADC2 Matrix for doubles-folding.")
Copy link
Contributor

Choose a reason for hiding this comment

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

I would call it:
"Doubles folding can currently only be applied to an ADC(2) matrix."

"""
return AmplitudeVector(ph=super().diagonal().ph)

def u2_fold(self, other):
Copy link
Contributor

Choose a reason for hiding this comment

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

Please rename it so that it is a bit more meaningful:
get_doubles_component(self, in_ampl)

u2 = self.block_apply("pphh_ph", other.ph) / (e * self.omega - diag)
return self.isymm.symmetrise(AmplitudeVector(pphh=u2))

def matvec(self, other):
Copy link
Contributor

Choose a reason for hiding this comment

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

same here, you can rename other to in_ampl to keep it more in line with the AdcMatrix implementation.

return AmplitudeVector(ph=self.block_apply("ph_ph", other.ph)
+ self.block_apply("ph_pphh", u2.pphh))

def unfold(self, other):
Copy link
Contributor

Choose a reason for hiding this comment

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

unfold_solution(self, in_ampl)

from .preconditioner import JacobiPreconditioner
from .SolverStateBase import EigenSolverStateBase
from .explicit_symmetrisation import IndexSymmetrisation
from .DIIS import DIIS
Copy link
Contributor

Choose a reason for hiding this comment

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

You forgot to add the DIIS.py file :)

# a block_view and then wrap that in an AdcMatrixProjected.


class AdcMatrixFolded(AdcMatrix):
Copy link
Contributor

Choose a reason for hiding this comment

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

I would move it into a separate file within the solver directory, since it’s currently quite solver-specific (for example, the unfold_solution function clearly refers to solving the eigenvalue problem handled in the solver) and specific to ADC(2). I’d also rename it to Adc2MatrixFolded.

However, we should consider how to generalize it in the future. :)

super().__init__(matrix.method, matrix.ground_state,
block_orders=matrix.block_orders,
intermediates=matrix.intermediates)
if matrix.method.name != "adc2":
Copy link
Contributor

Choose a reason for hiding this comment

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

As an idea for a more general approach: matrix.method.level % 2 !=0, because in general the ADC(4) matrix can also be folded.

- odiag["ph"].to_ndarray())) < 1e-12
assert not hasattr(fdiag, "pphh")

def test_u2_fold(self, system: str):
Copy link
Contributor

Choose a reason for hiding this comment

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

please rename it to test_get_doubles_component

diff = diff_manual - diff_matvec
assert np.max(np.abs(diff.to_ndarray())) < 1e-12

def test_unfold_random_vec(self,system:str):
Copy link
Contributor

Choose a reason for hiding this comment

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

I think that either this or the test_u2_fold are necessary as they essentially do the same.

"""
diag = super().diagonal().pphh
e = diag.ones_like()
u2 = self.block_apply("pphh_ph", other.ph) / (e * self.omega - diag)
Copy link
Contributor

@frieschneider frieschneider Oct 30, 2025

Choose a reason for hiding this comment

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

to get the correct symmetry: try to explicitly symmetrise it:

u2 = (
self.block_apply("pphh_ph", other.ph) / (e * self.omega - diag)
).antisymmetrise(0, 1).antisymmetrise(2, 3).symmetrise([(0, 1), (2, 3)])

@frieschneider
Copy link
Contributor

and please update it against the master :)

from adcc.AmplitudeVector import AmplitudeVector
print(ref_singlets_vec_s,type(ref_singlets_vec_s),type(res.eigenvectors[0]),res.eigenvectors[0])
ref_vec = AmplitudeVector(ph=ref_singlets_vec_s,pphh=ref_singlets_vec_d)
assert res.eigenvectors[0].ph.describe_symmetry() == ref_vec.ph.describe_symmetry()
Copy link
Contributor

Choose a reason for hiding this comment

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

Instead of using describe_symmetry to compare the symmetries, you can use the assert_equal_symmetry function from projection_test.py. It provides a more robust comparison.

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.

3 participants