Conversation
…od.at_level()' + prettying up output of 'ExcitedStates.describe_amplitudes()'
There was a name conflict leading to (only) wrong Adc3 equations. Intermediates could be grouped together
moved 'available_kinds' one level lower to 'method' instead of the system/case since it is method dependent adapted AdcMockState to have right attributes and gets the right blocks
Adapted ref data and read-in so that 'available_kinds' is tried to be set from system level (normal ref data) and from adc level for adcc ref data. De facto checks flexible where 'available_kinds' is stored
added 'pole_strength' tot the property blacklist for pp tests added 'is_alpha' attribute to Ex.St. object for the 'describe()' method to work properly
since PP-ADC calculations uses the same intermediates they are now only defined once
PP and IP/EA ADc calcs have different excitation properties Blacklist the non-existing ones
|
Added ip/ea adc2x for consistency with pp basemethods Added tests for IP/EA-ADC:
Modified adcc testdata generation: Modified read-in of reference data and creating MockAdcStates Generated ip/ea reference data for h2o_sto3g and cn_sto3g Removed some intermediates in ip and ea matrix.py to avoid conflicts and reuse code in pp matrix.py Introduced a blacklist for excitation_properties that depend on the type (pp vs ip/ea) All former tests are running |
maxscheurer
left a comment
maxscheurer
left a comment
There was a problem hiding this comment.
Looks nice I think! If you can come up with some minimal test for each method, that'd be great!
| __all__ = ["block"] | ||
|
|
||
| # TODO One thing one could still do to improve timings is implement a "fast einsum" | ||
| # that does not call opt_einsum, but directly dispatches to libadcc. This could | ||
| # lower the call overhead in the applies for the cases where we have only a | ||
| # trivial einsum to do. For the moment I'm not convinced that is worth the | ||
| # effort ... I suppose it only makes a difference for the cheaper ADC variants | ||
| # (ADC(0), ADC(1), CVS-ADC(0-2)-x), but then on the other hand they are not | ||
| # really so much our focus. | ||
|
|
||
|
|
||
| # | ||
| # Dispatch routine | ||
| # | ||
| """ | ||
| `apply` is a function mapping an AmplitudeVector to the contribution of this | ||
| block to the result of applying the ADC matrix. `diagonal` is an | ||
| `AmplitudeVector` containing the expression to the diagonal of the ADC matrix | ||
| from this block. | ||
| """ | ||
| AdcBlock = namedtuple("AdcBlock", ["apply", "diagonal"]) | ||
|
|
||
|
|
||
| def block(ground_state, spaces, order, variant=None, intermediates=None): | ||
| """ | ||
| Gets ground state, potentially intermediates, spaces (ph, pphh and so on) | ||
| and the perturbation theory order for the block, | ||
| variant is "cvs" or sth like that. | ||
|
|
||
| It is assumed largely, that CVS is equivalent to | ||
| mp.has_core_occupied_space, while one would probably want in the long run | ||
| that one can have an "o2" space, but not do CVS. | ||
| """ | ||
| if isinstance(variant, str): | ||
| variant = [variant] | ||
| elif variant is None: | ||
| variant = [] | ||
| reference_state = ground_state.reference_state | ||
| if intermediates is None: | ||
| intermediates = Intermediates(ground_state) | ||
|
|
||
| if ground_state.has_core_occupied_space and "cvs" not in variant: | ||
| raise ValueError("Cannot run a general (non-core-valence approximated)" | ||
| " ADC method on top of a ground state with a " | ||
| "core-valence separation.") | ||
| if not ground_state.has_core_occupied_space and "cvs" in variant: | ||
| raise ValueError("Cannot run a core-valence approximated ADC method on" | ||
| " top of a ground state without a " | ||
| "core-valence separation.") | ||
|
|
||
| fn = "_".join(["block"] + variant + spaces + [str(order)]) | ||
|
|
||
| if fn not in globals(): | ||
| raise ValueError("Could not dispatch: " | ||
| f"spaces={spaces} order={order} variant=variant") | ||
| return globals()[fn](reference_state, ground_state, intermediates) |
There was a problem hiding this comment.
Except for the globals() stuff, this is identical to the other matrix functions, right? I wonder if there's a way to reduce code duplication...
There was a problem hiding this comment.
What could be done is adding a new function in adc_pp/matrix.py or maybe in a new general file:
def get_block_prereqs(ground_state, spaces, order, variant=None,
intermediates=None):
if isinstance(variant, str):
variant = [variant]
elif variant is None:
variant = []
reference_state = ground_state.reference_state
if intermediates is None:
intermediates = Intermediates(ground_state)
if ground_state.has_core_occupied_space and "cvs" not in variant:
raise ValueError("Cannot run a general (non-core-valence approximated) "
"ADC method on top of a ground state with a "
"core-valence separation.")
if not ground_state.has_core_occupied_space and "cvs" in variant:
raise ValueError("Cannot run a core-valence approximated ADC method on "
"top of a ground state without a "
"core-valence separation.")
return ("_".join(["block"] + variant + spaces + [str(order)]),
reference_state)
and then importing it in the ip and ea matrix.py files and using the block routine as:
def block(ground_state, spaces, order, variant=None, intermediates=None):
fn, reference_state = get_block_prereqs(ground_state, spaces, order,
variant=None, intermediates=None)
if fn not in globals():
raise ValueError("Could not dispatch: "
f"spaces={spaces} order={order} variant=variant")
return globals()[fn](reference_state, ground_state, intermediates)
I think more cannot be reduced while keeping this structure. Also the
AdcBlock = namedtuple("AdcBlock", ["apply", "diagonal"])
seems necessary.
I'm not sure about the _all_= line since I think wildcard imports don't occur for the matrices. So if it is needed in general?
| if fn not in globals(): | ||
| raise ValueError("Could not dispatch: " | ||
| f"spaces={spaces} order={order} variant=variant") | ||
| return globals()[fn](reference_state, ground_state, intermediates) |
And thanks for reviewing :) Since I already added some tests, I'm not sure what you have in mind here. The obvious missing tests are testing against Q-Chem. However I struggled with that on the one hand with the interface and on the other hand due to differing different SCF solutions (psi4/pyscf vs qchem energies) that occured for some of the small set of open-shell test molecules leading to error propagation using more or less default inputs. But that should be doable. |
I was talking about testing against Q-Chem. That would be great, because then one can be sure nothing gets broken with the IP/EA feature. The |
added qchem ref data to yml file and added new test file probing ip/ea-adc(2/3) energies and pole strengths against qchem implementation
|
I reuse some of the code that appears in all three matrix files by adding a new function covering that which can be called in the other matrix files Further I added minimal testing against Q-Chem reference data, i.e. h2o-sto3g excitation energies and pole strengths for the first 2-3 states of IP/EA-ADC(2). I dumped the reference values in the qchem yml file |
|
I don't understand the current separation of methods between |
I agree with @jonasleitner. I also don't understand the current separation between |
|
From the top of my head, it can be rationalized as follows: Electronic excitation refers to a single excitation and the corresponding properties, whereas ExcitedStates is kind of a collection of ElectronicExcitations (list of structs vs struct of lists). That way, specific excitations can be addressed individually, for example when computing the gradient or so. Let's put this on the agenda for the next meeting, and I will look thru the code again. Please try to play around with the features of both classes also. I think there is a simple function on ExcitedStates that converts it to a list of ElectronicExcitations. |
I think you are confusing |
|
Ups sorry 😂 I was away from my computer, so didn't check again... |
|
|
||
| self.timer = Timer() | ||
| self.method = method | ||
| if "ip" in method.name: |
There was a problem hiding this comment.
I think this should be handled within the AdcMethod class. Maybe add a adc_type / adc_variant / ... member variable to the class. You can also keep a copy (probably with matching name) on the matrix class for better accessibility.
| split = method.split("-") | ||
| self.__base_method = split[-1] | ||
| split = split[:-1] | ||
| self.__base_method = split.pop() |
There was a problem hiding this comment.
I think base_method should just be adc2. The adc variant (pp/ip/ea) should be stored in a separate variable that needs to be taken into account below (at_level and name) to reassemble the correct method string.
| but not do CVS. | ||
| """ | ||
| fn, reference_state = get_block_prereqs( | ||
| ground_state, spaces, order, variant, intermediates) |
There was a problem hiding this comment.
I think it is a good idea to use a common function to build the function name and to perform some common sanity checks. But it should only return the function name (we have access to the reference_state through the ground state). Probably something like get_block_name is a more fitting name then.
Also a small bug: If intermediates is passed as None the initialization currently happens in get_block_prereqs and is not returned. Probably the lines
if intermediates is None:
intermediates = Intermediates(ground_state)
should be moved back into the block function
| `AmplitudeVector` containing the expression to the diagonal of the ADC matrix | ||
| from this block. | ||
| """ | ||
| AdcBlock = namedtuple("AdcBlock", ["apply", "diagonal"]) |
There was a problem hiding this comment.
It should be possible to import AdcBlock from adc_pp.matrix (same for ip)
| if fn not in globals(): | ||
| raise ValueError("Could not dispatch: " | ||
| f"spaces={spaces} order={order} variant=variant") | ||
| return globals()[fn](reference_state, ground_state, intermediates) |
There was a problem hiding this comment.
Maybe we should raise an explicit exception if cvs is defined for EA/IP, since get_block_prereqs only checks if it is consistently set in the ground_state and the variant.
| # i.e. user-provided, we cannot guarantee for obtaining a particular | ||
| # spin_change in case of a spin_flip calculation. | ||
| spin_change = None | ||
| spin_change = 0 |
There was a problem hiding this comment.
I think it would be better to set spin_change to the default value None if guesses are provided. Since we don't run the guess function in this case we should be fine. (If it can't be removed entirely from here)
| spin_flip=dict(spin_block_symmetrisation="none", spin_change=-1), | ||
| any=dict(spin_block_symmetrisation="none", spin_change=0), | ||
| singlet=dict(spin_block_symmetrisation="symmetric"), | ||
| doublet=dict(spin_block_symmetrisation="none"), |
There was a problem hiding this comment.
Should it not be possible to give guess_kwargs_kind is_alpha to determine the spin change for IP and EA? Then it would not be necessary to forward spin_change all the way from run_adc.
Maybe the spin change can then be fully handled down here? (Also for spin_flip?)
| "ip_adc0": pole_strength_ip_adc0, | ||
| "ip_adc1": pole_strength_ip_adc0, | ||
| "ip_adc2": pole_strength_ip_adc2, | ||
| "ip_adc3": pole_strength_ip_adc3, |
There was a problem hiding this comment.
If I remember correctly adc3 pole strengths are implemented but not used yet, right?
So PP-, IP- and EA-ADC(3) all consistently use ISR2 for properties?
Rename to ip-adcn (see AdcMethod.py)
| valid_bases = ["adc0", "adc1", "adc2", "adc2x", "adc3"] | ||
| valid_bases = ["adc0", "adc1", "adc2", "adc2x", "adc3", | ||
| "ip_adc0", "ip_adc1", "ip_adc2", "ip_adc2x", "ip_adc3", | ||
| "ea_adc0", "ea_adc1", "ea_adc2", "ea_adc2x", "ea_adc3"] |
There was a problem hiding this comment.
Intuitively I would expect ip-adc2 to be a valid method (- vs _). Probably it would be best to accept both (for instance with a replace in AdcMethod.__init__). However, internally we should only use one form. I think ip-adcn to stay consistent with cvs-adcn.
| "ea_adc1": dict(p_p=1, p_pph=None, pph_p=None, pph_pph=None), | ||
| "ea_adc2": dict(p_p=2, p_pph=1, pph_p=1, pph_pph=0), | ||
| "ea_adc2x": dict(p_p=2, p_pph=1, pph_p=1, pph_pph=1), | ||
| "ea_adc3": dict(p_p=3, p_pph=2, pph_p=2, pph_pph=1), |
There was a problem hiding this comment.
See my comment in AdcMethod.py. I think we should consistently use ip-adcn and ea-adcn internally.
4 participants
New methods on user side:
ip_adc{0..3}()
ea_adc{0..3}()
New keywords on user side:
n_doublets, is_alpha
e.g. ip_adc3(scf_result, n_doublets=5, is_alpha=True)
Remarks:
Energies and properties were verified against Q-Chem calculations.
Pole strengths of 3rd order were implemented for completeness. They are currently not used in the calculation since all properties are computed at most at 2nd order even for ADC(3) calculations. Hence they were not verified yet.
No simultaneous calculation of alpha and beta states to simplify the code and error handling and to keep consistency with PP calculations. If alpha and beta would be computed in the same adcc instance, intermediates could be reused. However, these intermediates are not the bottleneck of the calculation and thus recomputed for the other spin.
Quartets not yet implemented since they are pure doubles states which are inaccurate up to 3rd order.
To Do
Testing procedures for IP/EA ADC
Verify 3rd order pole strengths
Observed issue
True for PP/IP/EA
If 'kind="any"' and the user requests a significant amount of states getting somewhere near the number of total accessible states, the Davidson will quickly converge to 0-vector. This is not true if kind is specified.
Example input