diff --git a/fanpy/wfn/cc/ap1rog_generalized.py b/fanpy/wfn/cc/ap1rog_generalized.py index ccda9ac8..1d255333 100644 --- a/fanpy/wfn/cc/ap1rog_generalized.py +++ b/fanpy/wfn/cc/ap1rog_generalized.py @@ -1,24 +1,40 @@ -"""APG1ro wavefunction with single and double excitations.""" +"""AP1roG wavefunction with generalized singles and paired double excitations.""" from fanpy.tools import slater from fanpy.wfn.cc.pccd_ap1rog import PCCD class AP1roGSDGeneralized(PCCD): - r"""AP1roG wavefunction with single and double excitations, broken spin symmetry. + r"""AP1roG wavefunction with generalized singles and paired doubles. + + The wavefunction is parameterized as .. math:: - \left| {{\Psi }_{APG1roSD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a,b\in virt}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\prod\limits_{i=1}^{N/2\;}{\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{\bar{i};a}} - \hat{\tau }_{i\bar{i}}^{ia}} \right)\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{i;a}}\hat{\tau }_{i\bar{i}}^{a\bar{i}}} - \right)}\left| {{\Phi }_{0}} \right\rangle } - - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. - + \left| \Psi_{\mathrm{AP1roGSDGeneralized}} \right\rangle + = \prod_{\mu \in \mathcal{E}} + \left( 1 + t_\mu\, \tau_\mu \right) \left| \Phi_0 \right\rangle, + \qquad + \mathcal{E} = \left\{ \tau_{i}^{a},\, \tau_{\bar{i}}^{a},\, + \tau_{i}^{\bar{a}},\, \tau_{\bar{i}}^{\bar{a}},\, + \tau_{i\bar{i}}^{a\bar{a}} \right\}, + + where :math:`i, j, k, \ldots` index occupied spin orbitals and + :math:`a, b, c, \ldots` index virtual spin orbitals. The excitation + pool :math:`\mathcal{E}` combines paired double excitations + :math:`\tau_{i\bar{i}}^{a\bar{a}}` inherited from PCCD with + generalized single excitations that permit spin-flip and + spin-symmetry-breaking transitions, in contrast to + :class:`AP1roGSDSpin` where singles are restricted to + spin-conserving form. The single excitation manifold can be further + constrained through the ``s_type`` keyword inherited from the PCCD + base class. + + The reference wavefunction (``refwfn``) must be a seniority-0 Slater + determinant. + + See Ref. [1]_ for the full theoretical background. + Attributes ---------- nelec : int @@ -42,6 +58,11 @@ class AP1roGSDGeneralized(PCCD): dictionary, the keys are tuples with the indices of annihilation and creation operators, and the values are the excitation operators that allow to excite from the annihilation to the creation operators. + s_type : str + Option controlling the screening of single excitations during + overlap evaluation. Inherited from :class:`PCCD`. + See PCCD.assign_s_type for available options and their definitions. + Properties ---------- @@ -82,6 +103,8 @@ class AP1roGSDGeneralized(PCCD): Assign the reference wavefunction. assign_params(self, params=None, add_noise=False) Assign the parameters of the CC wavefunction. + assign_s_type(self, s_type): + Assign the option of seniority-breaking condition used for single excitations. get_ind(self, exop) : int Return the parameter index that corresponds to a given excitation operator. get_exop(self, ind) : list of int @@ -99,6 +122,12 @@ class AP1roGSDGeneralized(PCCD): Assign the excitation operators that can excite from the given indices to be annihilated to the given indices to be created. + References + ---------- + .. [1] P. B. Gaikwad, T. D. Kim, M. Richer, R. A. Lokhande, G. Sánchez-Díaz; + P. A. Limacher, P. W. Ayers and R. A. Miranda-Quintana, "Coupled-cluster-inspired + geminal wavefunctions," *J. Chem. Phys.* **160**, 144108 (2024). + https://doi.org/10.1063/5.0196561 """ def assign_ranks(self, ranks=None): diff --git a/fanpy/wfn/cc/ap1rog_generalized_NEW.py b/fanpy/wfn/cc/ap1rog_generalized_NEW.py deleted file mode 100644 index 547dcae5..00000000 --- a/fanpy/wfn/cc/ap1rog_generalized_NEW.py +++ /dev/null @@ -1,171 +0,0 @@ -"""APG1ro wavefunction with single and double excitations.""" - -from fanpy.tools import slater -from fanpy.wfn.cc.pccd_ap1rog_NEW import PCCD - - -class AP1roGSDGeneralized(PCCD): - r"""AP1roG wavefunction with single and double excitations, broken spin symmetry. - - .. math:: - - \left| {{\Psi }_{APG1roSD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a,b\in virt}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\prod\limits_{i=1}^{N/2\;}{\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{\bar{i};a}} - \hat{\tau }_{i\bar{i}}^{ia}} \right)\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{i;a}}\hat{\tau }_{i\bar{i}}^{a\bar{i}}} - \right)}\left| {{\Phi }_{0}} \right\rangle } - - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. - - Attributes - ---------- - nelec : int - Number of electrons. - nspin : int - Number of spin orbitals (alpha and beta). - dtype : {np.float64, np.complex128} - Data type of the wavefunction. - params : np.ndarray - Parameters of the wavefunction. - memory : float - Memory available for the wavefunction. - ranks : list of ints - Ranks of the excitation operators. - exops : list of list of int - Excitation operators given as lists of ints. The first half of indices correspond to - indices to be annihilated, the second half correspond to indices to be created. - refwfn : int - Reference wavefunction upon which the CC operator will act. - exop_combinations : dict - dictionary, the keys are tuples with the indices of annihilation and creation - operators, and the values are the excitation operators that allow to excite from the - annihilation to the creation operators. - - Properties - ---------- - nparams : int - Number of parameters. - nspatial : int - Number of spatial orbitals - param_shape : tuple of int - Shape of the parameters. - spin : int - Spin of the wavefunction. - seniority : int - Seniority of the wavefunction. - template_params : np.ndarray - Default parameters of the wavefunction. - nexops : int - Number of excitation operators. - nranks : int - Number of allowed ranks. - - Methods - ------- - __init__(self, nelec, nspin, memory=None, ngem=None, orbpairs=None, params=None) - Initialize the wavefunction. - assign_nelec(self, nelec) - Assign the number of electrons. - assign_nspin(self, nspin) - Assign the number of spin orbitals. - assign_dtype(self, dtype) - Assign the data type of the parameters. - assign_memory(self, memory=None) - Assign memory available for the wavefunction. - assign_ranks(self, ranks=None) - Assign the allowed excitation ranks. - assign_exops(self, exops=None) - Assign the allowed excitation operators. - assign_refwfn(self, refwfn=None) - Assign the reference wavefunction. - assign_params(self, params=None, add_noise=False) - Assign the parameters of the CC wavefunction. - get_ind(self, exop) : int - Return the parameter index that corresponds to a given excitation operator. - get_exop(self, ind) : list of int - Return the excitation operator that corresponds to a given parameter index. - product_amplitudes(self, inds, deriv=None) : float - Return the product of the CC amplitudes of the coefficients corresponding to - the given indices. - load_cache(self) - Load the functions whose values will be cached. - clear_cache(self) - Clear the cache. - get_overlap(self, sd, deriv=None) : float - Return the overlap of the wavefunction with a Slater determinant. - generate_possible_exops(self, a_inds, c_inds): - Assign the excitation operators that can excite from the given indices to be annihilated - to the given indices to be created. - - """ - - def assign_ranks(self, ranks=None): - """Assign the ranks of the excitation operators. - - Parameters - ---------- - ranks : {int, list, None} - Ranks of the allowed excitation operators. - - Raises - ------ - TypeError - If 'ranks' is not None. - - """ - if ranks is not None: - raise TypeError("Only the default: ranks=[1, 2] is allowed") - else: - self.ranks = [1, 2] - - def assign_exops(self, indices=None): - """Assign the excitation operators that will be used to construct the CC operator. - - Parameters - ---------- - indices : None - The allowed excitation operators are solely defined by the occupied and virtual - orbitals of the given reference Slater determinant. - - Raises - ------ - TypeError - If `indices` is not None. - - Notes - ----- - The excitation operators are given as a list of lists of ints. - Each sub-list corresponds to an excitation operator. - In each sub-list, the first half of indices corresponds to the indices of the - spin-orbitals to annihilate, and the second half corresponds to the indices of the - spin-orbitals to create. - [a1, a2, ..., aN, c1, c2, ..., cN] - - """ - if indices is not None: - raise TypeError( - "Only the excitation operators constructed by default from " - "the given reference Slater determinant are allowed" - ) - else: - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - ex_to = [i for i in range(self.nspin) if i not in ex_from] - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt in ex_to: - exop = [occ_alpha, virt] - exops[tuple(exop)] = counter - counter += 1 - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt in ex_to: - exop = [occ_alpha + self.nspatial, virt] - exops[tuple(exop)] = counter - counter += 1 - super().assign_exops(None) - shift = len(exops) - for exop, ind in self.exops.items(): - exops[exop] = ind + shift - self.exops = exops diff --git a/fanpy/wfn/cc/ap1rog_spin.py b/fanpy/wfn/cc/ap1rog_spin.py index a2e7ba29..56047723 100644 --- a/fanpy/wfn/cc/ap1rog_spin.py +++ b/fanpy/wfn/cc/ap1rog_spin.py @@ -1,23 +1,37 @@ -"""APG1ro wavefunction with single and double excitations.""" +"""AP1roG wavefunction with single and paired-double excitations.""" from fanpy.tools import slater from fanpy.wfn.cc.pccd_ap1rog import PCCD class AP1roGSDSpin(PCCD): - r"""AP1roG wavefunction with single and double excitations, correct spin symmetry. + r"""AP1roG wavefunction with single and paired-double excitations, both preserving spin symmetry. + + NOTE: The excitation operator pool in this wavefunction consists of paired doubles + inherited from PCCD together with spin-conserving single excitations. .. math:: - \left| {{\Psi }_{APG1roSD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a,b\in virt}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\prod\limits_{i=1}^{N/2\;}{\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{\bar{i};a}} - \hat{\tau }_{i\bar{i}}^{ia}} \right)\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{i;a}}\hat{\tau }_{i\bar{i}}^{a\bar{i}}} - \right)}\left| {{\Phi }_{0}} \right\rangle } + \left| \Psi _{\mathrm{AP1roGSDSpin}} \right\rangle = \prod _{\mu \in \mathcal{E}} + \left( 1+t_{\mu }\tilde{\tau }_{\mu }\right)| \Phi _{0} \rangle - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. + where, + + .. math:: + + \mathcal{E} = \left\{\tau_i^a,\tau_{\bar{i}}^{\bar{a}},\tau_{i\bar{i}}^{a\bar{a}}\right\} + + The excitation operator pool :math:\mathcal{E} consists of + - spin-conserving single excitations + :math:(\tau_i^a,\tau_{\bar{i}}^{\bar{a}}) + - paired double excitations + :math:(\tau_{i\bar{i}}^{a\bar{a}}) + generated from the reference determinant. + + The effective single excitation operators :math:\tilde{\tau} may additionally be + constrained during overlap evaluation through the `s_type` option implemented in the PCCD wavefunction class. + + The reference wavefunction can only be a single Slater determinant with seniority-0. Attributes ---------- @@ -42,6 +56,11 @@ class AP1roGSDSpin(PCCD): dictionary, the keys are tuples with the indices of annihilation and creation operators, and the values are the excitation operators that allow to excite from the annihilation to the creation operators. + s_type : str + Option controlling the screening of single excitations during + overlap evaluation. Inherited from :class:`PCCD`. + See PCCD.assign_s_type for available options and their definitions. + Properties ---------- @@ -82,6 +101,8 @@ class AP1roGSDSpin(PCCD): Assign the reference wavefunction. assign_params(self, params=None, add_noise=False) Assign the parameters of the CC wavefunction. + assign_s_type(self, s_type): + Assign the option of seniority-breaking condition used for single excitations. get_ind(self, exop) : int Return the parameter index that corresponds to a given excitation operator. get_exop(self, ind) : list of int @@ -99,6 +120,12 @@ class AP1roGSDSpin(PCCD): Assign the excitation operators that can excite from the given indices to be annihilated to the given indices to be created. + References + ---------- + .. [1] P. B. Gaikwad, T. D. Kim, M. Richer, R. A. Lokhande, G. Sánchez-Díaz; + P. A. Limacher, P. W. Ayers and R. A. Miranda-Quintana, "Coupled-cluster-inspired + geminal wavefunctions," *J. Chem. Phys.* **160**, 144108 (2024). + https://doi.org/10.1063/5.0196561 """ def assign_ranks(self, ranks=None): @@ -143,6 +170,10 @@ def assign_exops(self, indices=None): spin-orbitals to create. [a1, a2, ..., aN, c1, c2, ..., cN] + + This method constructs only a static pool of excitation operators derived from + the reference determinant. The overlap routines later generate compatible combinations + of these excitation operators during determinant connections. """ if indices is not None: raise TypeError( diff --git a/fanpy/wfn/cc/ap1rog_spin_NEW.py b/fanpy/wfn/cc/ap1rog_spin_NEW.py deleted file mode 100644 index fd6fd0a9..00000000 --- a/fanpy/wfn/cc/ap1rog_spin_NEW.py +++ /dev/null @@ -1,172 +0,0 @@ -"""APG1ro wavefunction with single and double excitations.""" - -from fanpy.tools import slater -from fanpy.wfn.cc.pccd_ap1rog_NEW import PCCD - - -class AP1roGSDSpin(PCCD): - r"""AP1roG wavefunction with single and double excitations, correct spin symmetry. - - .. math:: - - \left| {{\Psi }_{APG1roSD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a,b\in virt}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\prod\limits_{i=1}^{N/2\;}{\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{\bar{i};a}} - \hat{\tau }_{i\bar{i}}^{ia}} \right)\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{i;a}}\hat{\tau }_{i\bar{i}}^{a\bar{i}}} - \right)}\left| {{\Phi }_{0}} \right\rangle } - - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. - - Attributes - ---------- - nelec : int - Number of electrons. - nspin : int - Number of spin orbitals (alpha and beta). - dtype : {np.float64, np.complex128} - Data type of the wavefunction. - params : np.ndarray - Parameters of the wavefunction. - memory : float - Memory available for the wavefunction. - ranks : list of ints - Ranks of the excitation operators. - exops : list of list of int - Excitation operators given as lists of ints. The first half of indices correspond to - indices to be annihilated, the second half correspond to indices to be created. - refwfn : int - Reference wavefunction upon which the CC operator will act. - exop_combinations : dict - dictionary, the keys are tuples with the indices of annihilation and creation - operators, and the values are the excitation operators that allow to excite from the - annihilation to the creation operators. - - Properties - ---------- - nparams : int - Number of parameters. - nspatial : int - Number of spatial orbitals - param_shape : tuple of int - Shape of the parameters. - spin : int - Spin of the wavefunction. - seniority : int - Seniority of the wavefunction. - template_params : np.ndarray - Default parameters of the wavefunction. - nexops : int - Number of excitation operators. - nranks : int - Number of allowed ranks. - - Methods - ------- - __init__(self, nelec, nspin, memory=None, ngem=None, orbpairs=None, params=None) - Initialize the wavefunction. - assign_nelec(self, nelec) - Assign the number of electrons. - assign_nspin(self, nspin) - Assign the number of spin orbitals. - assign_dtype(self, dtype) - Assign the data type of the parameters. - assign_memory(self, memory=None) - Assign memory available for the wavefunction. - assign_ranks(self, ranks=None) - Assign the allowed excitation ranks. - assign_exops(self, exops=None) - Assign the allowed excitation operators. - assign_refwfn(self, refwfn=None) - Assign the reference wavefunction. - assign_params(self, params=None, add_noise=False) - Assign the parameters of the CC wavefunction. - get_ind(self, exop) : int - Return the parameter index that corresponds to a given excitation operator. - get_exop(self, ind) : list of int - Return the excitation operator that corresponds to a given parameter index. - product_amplitudes(self, inds, deriv=None) : float - Return the product of the CC amplitudes of the coefficients corresponding to - the given indices. - load_cache(self) - Load the functions whose values will be cached. - clear_cache(self) - Clear the cache. - get_overlap(self, sd, deriv=None) : float - Return the overlap of the wavefunction with a Slater determinant. - generate_possible_exops(self, a_inds, c_inds): - Assign the excitation operators that can excite from the given indices to be annihilated - to the given indices to be created. - - """ - - def assign_ranks(self, ranks=None): - """Assign the ranks of the excitation operators. - - Parameters - ---------- - ranks : {int, list, None} - Ranks of the allowed excitation operators. - - Raises - ------ - TypeError - If 'ranks' is not None. - - """ - if ranks is not None: - raise TypeError("Only the default: ranks=[1, 2] is allowed") - else: - self.ranks = [1, 2] - - def assign_exops(self, indices=None): - """Assign the excitation operators that will be used to construct the CC operator. - - Parameters - ---------- - indices : None - The allowed excitation operators are solely defined by the occupied and virtual - orbitals of the given reference Slater determinant. - - Raises - ------ - TypeError - If `indices` is not None. - - Notes - ----- - The excitation operators are given as a list of lists of ints. - Each sub-list corresponds to an excitation operator. - In each sub-list, the first half of indices corresponds to the indices of the - spin-orbitals to annihilate, and the second half corresponds to the indices of the - spin-orbitals to create. - [a1, a2, ..., aN, c1, c2, ..., cN] - - """ - if indices is not None: - raise TypeError( - "Only the excitation operators constructed by default from " - "the given reference Slater determinant are allowed" - ) - else: - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - virt_alphas = [i for i in range(self.nspin) if (i not in ex_from) and slater.is_alpha(i, self.nspatial)] - virt_betas = [i for i in range(self.nspin) if (i not in ex_from) and not slater.is_alpha(i, self.nspatial)] - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt_alpha in virt_alphas: - exop = [occ_alpha, virt_alpha] - exops[tuple(exop)] = counter - counter += 1 - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt_beta in virt_betas: - exop = [occ_alpha + self.nspatial, virt_beta] - exops[tuple(exop)] = counter - counter += 1 - super().assign_exops(None) - shift = len(exops) - for exop, ind in self.exops.items(): - exops[exop] = ind + shift - self.exops = exops diff --git a/fanpy/wfn/cc/apg1ro_d.py b/fanpy/wfn/cc/apg1ro_d.py index 8a515666..a9a8af2f 100644 --- a/fanpy/wfn/cc/apg1ro_d.py +++ b/fanpy/wfn/cc/apg1ro_d.py @@ -5,16 +5,32 @@ class APG1roD(PCCD): - r"""APG1ro wavefunction with only double excitations. + r"""APG1ro wavefunction with generalized double excitations. + + The wavefunction is parameterized as .. math:: - \[\left| {{\Psi }_{APG1roD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a,b\in virt}^{{}} - {{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} \right)}\left| {{\Phi }_{0}} \right\rangle \] + \left| \Psi_{\mathrm{APG1roD}} \right\rangle + = \prod_{\mu \in \mathcal{E}} + \left( 1 + t_\mu\, \tau_\mu \right) \left| \Phi_0 \right\rangle, + \qquad + \mathcal{E} = \left\{ \tau_{i\bar{i}}^{ab} \right\}, + + where :math:`i, j, k, \ldots` index occupied spin orbitals, and + :math:`a, b, c, \ldots` index virtual spin orbitals. The excitation pool + :math:`\mathcal{E}` consists of generalized double excitations that + annihilate a paired occupied pair :math:`(i, \bar{i})` and create + electrons in arbitrary virtual spin orbitals :math:`(a, b)`. This + covers both pair-preserving excitations (:math:`b = \bar{a}`) + and pair-breaking excitations (:math:`b \neq = \bar{a}`); + because the virtual creation indices are unrestricted, + the ansatz may break spin symmetry. - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. + The reference wavefunction (``refwfn``) must be a seniority-0 Slater + determinant. + + See Ref. [1]_ for the full theoretical background. Attributes ---------- @@ -96,6 +112,12 @@ class APG1roD(PCCD): Assign the excitation operators that can excite from the given indices to be annihilated to the given indices to be created. + References + ---------- + .. [1] P. B. Gaikwad, T. D. Kim, M. Richer, R. A. Lokhande, G. Sánchez-Díaz; + P. A. Limacher, P. W. Ayers and R. A. Miranda-Quintana, "Coupled-cluster-inspired + geminal wavefunctions," *J. Chem. Phys.* **160**, 144108 (2024). + https://doi.org/10.1063/5.0196561 """ def assign_exops(self, indices=None): diff --git a/fanpy/wfn/cc/apg1ro_d_NEW.py b/fanpy/wfn/cc/apg1ro_d_NEW.py deleted file mode 100644 index 802c95f6..00000000 --- a/fanpy/wfn/cc/apg1ro_d_NEW.py +++ /dev/null @@ -1,141 +0,0 @@ -"""APG1ro wavefunction with only double excitations.""" - -from fanpy.tools import slater -from fanpy.wfn.cc.pccd_ap1rog_NEW import PCCD - - -class APG1roD(PCCD): - r"""APG1ro wavefunction with only double excitations. - - .. math:: - - \[\left| {{\Psi }_{APG1roD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a,b\in virt}^{{}} - {{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} \right)}\left| {{\Phi }_{0}} \right\rangle \] - - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. - - Attributes - ---------- - nelec : int - Number of electrons. - nspin : int - Number of spin orbitals (alpha and beta). - dtype : {np.float64, np.complex128} - Data type of the wavefunction. - params : np.ndarray - Parameters of the wavefunction. - memory : float - Memory available for the wavefunction. - ranks : list of ints - Ranks of the excitation operators. - exops : list of list of int - Excitation operators given as lists of ints. The first half of indices correspond to - indices to be annihilated, the second half correspond to indices to be created. - refwfn : int - Reference wavefunction upon which the CC operator will act. - exop_combinations : dict - dictionary, the keys are tuples with the indices of annihilation and creation - operators, and the values are the excitation operators that allow to excite from the - annihilation to the creation operators. - - Properties - ---------- - nparams : int - Number of parameters. - nspatial : int - Number of spatial orbitals - param_shape : tuple of int - Shape of the parameters. - spin : int - Spin of the wavefunction. - seniority : int - Seniority of the wavefunction. - template_params : np.ndarray - Default parameters of the wavefunction. - nexops : int - Number of excitation operators. - nranks : int - Number of allowed ranks. - - Methods - ------- - __init__(self, nelec, nspin, memory=None, ngem=None, orbpairs=None, params=None) - Initialize the wavefunction. - assign_nelec(self, nelec) - Assign the number of electrons. - assign_nspin(self, nspin) - Assign the number of spin orbitals. - assign_dtype(self, dtype) - Assign the data type of the parameters. - assign_memory(self, memory=None) - Assign memory available for the wavefunction. - assign_ranks(self, ranks=None) - Assign the allowed excitation ranks. - assign_exops(self, exops=None) - Assign the allowed excitation operators. - assign_refwfn(self, refwfn=None) - Assign the reference wavefunction. - assign_params(self, params=None, add_noise=False) - Assign the parameters of the CC wavefunction. - get_ind(self, exop) : int - Return the parameter index that corresponds to a given excitation operator. - get_exop(self, ind) : list of int - Return the excitation operator that corresponds to a given parameter index. - product_amplitudes(self, inds, deriv=None) : float - Return the product of the CC amplitudes of the coefficients corresponding to - the given indices. - load_cache(self) - Load the functions whose values will be cached. - clear_cache(self) - Clear the cache. - get_overlap(self, sd, deriv=None) : float - Return the overlap of the wavefunction with a Slater determinant. - generate_possible_exops(self, a_inds, c_inds): - Assign the excitation operators that can excite from the given indices to be annihilated - to the given indices to be created. - - """ - - def assign_exops(self, indices=None): - """Assign the excitation operators that will be used to construct the CC operator. - - Parameters - ---------- - indices : None - The allowed excitation operators are solely defined by the occupied and virtual - orbitals of the given reference Slater determinant. - - Raises - ------ - TypeError - If `indices` is not None. - - Notes - ----- - The excitation operators are given as a list of lists of ints. - Each sub-list corresponds to an excitation operator. - In each sub-list, the first half of indices corresponds to the indices of the - spin-orbitals to annihilate, and the second half corresponds to the indices of the - spin-orbitals to create. - [a1, a2, ..., aN, c1, c2, ..., cN] - - """ - if indices is not None: - raise TypeError( - "Only the excitation operators constructed by default from " - "the given reference Slater determinant are allowed" - ) - else: - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - ex_to = [i for i in range(self.nspin) if i not in ex_from] - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt1 in ex_to: - for virt2 in ex_to[ex_to.index(virt1) + 1 :]: - exop = [occ_alpha, occ_alpha + self.nspatial, virt1, virt2] - exops[tuple(exop)] = counter - counter += 1 - self.exops = exops diff --git a/fanpy/wfn/cc/apg1ro_sd.py b/fanpy/wfn/cc/apg1ro_sd.py index 42887475..2d8eda26 100644 --- a/fanpy/wfn/cc/apg1ro_sd.py +++ b/fanpy/wfn/cc/apg1ro_sd.py @@ -5,19 +5,35 @@ class APG1roSD(APG1roD): - r"""APG1ro wavefunction with single and double excitations. + r"""APG1ro wavefunction with generalized single and double excitations. + + The wavefunction is parameterized as .. math:: - \left| {{\Psi }_{APG1roSD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a,b\in virt}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\prod\limits_{i=1}^{N/2\;}{\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{\bar{i};a}} - \hat{\tau }_{i\bar{i}}^{ia}} \right)\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{i;a}}\hat{\tau }_{i\bar{i}}^{a\bar{i}}} - \right)}\left| {{\Phi }_{0}} \right\rangle } + \left| \Psi_{\mathrm{APG1roSD}} \right\rangle + = \prod_{\mu \in \mathcal{E}} + \left( 1 + t_\mu\, \tau_\mu \right) \left| \Phi_0 \right\rangle, + \qquad + \mathcal{E} = \left\{ \tau_{i}^{a}, \tau_{i\bar{i}}^{ab} \right\}, + + where :math:`i, j, k, \ldots` index occupied spin orbitals and + :math:`a, b, c, \ldots` index virtual spin orbitals. The excitation + pool :math:`\mathcal{E}` extends APG1roD with generalized single + excitations :math:`\tau_{i}^{a}` alongside the generalized double + excitations :math:`\tau_{i\bar{i}}^{ab}` inherited from APG1roD. + The doubles cover both pair-preserving (:math:`b = \bar{a}`) and + pair-breaking (:math:`b \neq \bar{a}`) excitations. Unlike AP1roGSDSpin, + the single excitations are not restricted to spin-conserving form, + allowing spin-flip excitations and broken spin symmetry. The single + excitation manifold can be further constrained through the ``s_type`` + keyword inherited from the PCCD base class. + + The reference wavefunction (``refwfn``) must be a seniority-0 Slater + determinant. + + See Ref. [1]_ for the full theoretical background. - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. Attributes ---------- @@ -42,6 +58,11 @@ class APG1roSD(APG1roD): dictionary, the keys are tuples with the indices of annihilation and creation operators, and the values are the excitation operators that allow to excite from the annihilation to the creation operators. + s_type : str + Option controlling the screening of single excitations during + overlap evaluation. Inherited from :class:`PCCD`. + See PCCD.assign_s_type for available options and their definitions. + Properties ---------- @@ -82,6 +103,8 @@ class APG1roSD(APG1roD): Assign the reference wavefunction. assign_params(self, params=None, add_noise=False) Assign the parameters of the CC wavefunction. + assign_s_type(self, s_type): + Assign the option of seniority-breaking condition used for single excitations. get_ind(self, exop) : int Return the parameter index that corresponds to a given excitation operator. get_exop(self, ind) : list of int @@ -99,6 +122,12 @@ class APG1roSD(APG1roD): Assign the excitation operators that can excite from the given indices to be annihilated to the given indices to be created. + References + ---------- + .. [1] P. B. Gaikwad, T. D. Kim, M. Richer, R. A. Lokhande, G. Sánchez-Díaz; + P. A. Limacher, P. W. Ayers and R. A. Miranda-Quintana, "Coupled-cluster-inspired + geminal wavefunctions," *J. Chem. Phys.* **160**, 144108 (2024). + https://doi.org/10.1063/5.0196561 """ def assign_ranks(self, ranks=None): diff --git a/fanpy/wfn/cc/apg1ro_sd_NEW.py b/fanpy/wfn/cc/apg1ro_sd_NEW.py deleted file mode 100644 index 8b9b624f..00000000 --- a/fanpy/wfn/cc/apg1ro_sd_NEW.py +++ /dev/null @@ -1,171 +0,0 @@ -"""APG1ro wavefunction with single and double excitations.""" - -from fanpy.tools import slater -from fanpy.wfn.cc.apg1ro_d_NEW import APG1roD - - -class APG1roSD(APG1roD): - r"""APG1ro wavefunction with single and double excitations. - - .. math:: - - \left| {{\Psi }_{APG1roSD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a,b\in virt}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\prod\limits_{i=1}^{N/2\;}{\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{\bar{i};a}} - \hat{\tau }_{i\bar{i}}^{ia}} \right)\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a\in virt}^{{}}{{{t}_{i;a}}\hat{\tau }_{i\bar{i}}^{a\bar{i}}} - \right)}\left| {{\Phi }_{0}} \right\rangle } - - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. - - Attributes - ---------- - nelec : int - Number of electrons. - nspin : int - Number of spin orbitals (alpha and beta). - dtype : {np.float64, np.complex128} - Data type of the wavefunction. - params : np.ndarray - Parameters of the wavefunction. - memory : float - Memory available for the wavefunction. - ranks : list of ints - Ranks of the excitation operators. - exops : list of list of int - Excitation operators given as lists of ints. The first half of indices correspond to - indices to be annihilated, the second half correspond to indices to be created. - refwfn : int - Reference wavefunction upon which the CC operator will act. - exop_combinations : dict - dictionary, the keys are tuples with the indices of annihilation and creation - operators, and the values are the excitation operators that allow to excite from the - annihilation to the creation operators. - - Properties - ---------- - nparams : int - Number of parameters. - nspatial : int - Number of spatial orbitals - param_shape : tuple of int - Shape of the parameters. - spin : int - Spin of the wavefunction. - seniority : int - Seniority of the wavefunction. - template_params : np.ndarray - Default parameters of the wavefunction. - nexops : int - Number of excitation operators. - nranks : int - Number of allowed ranks. - - Methods - ------- - __init__(self, nelec, nspin, memory=None, ngem=None, orbpairs=None, params=None) - Initialize the wavefunction. - assign_nelec(self, nelec) - Assign the number of electrons. - assign_nspin(self, nspin) - Assign the number of spin orbitals. - assign_dtype(self, dtype) - Assign the data type of the parameters. - assign_memory(self, memory=None) - Assign memory available for the wavefunction. - assign_ranks(self, ranks=None) - Assign the allowed excitation ranks. - assign_exops(self, exops=None) - Assign the allowed excitation operators. - assign_refwfn(self, refwfn=None) - Assign the reference wavefunction. - assign_params(self, params=None, add_noise=False) - Assign the parameters of the CC wavefunction. - get_ind(self, exop) : int - Return the parameter index that corresponds to a given excitation operator. - get_exop(self, ind) : list of int - Return the excitation operator that corresponds to a given parameter index. - product_amplitudes(self, inds, deriv=None) : float - Return the product of the CC amplitudes of the coefficients corresponding to - the given indices. - load_cache(self) - Load the functions whose values will be cached. - clear_cache(self) - Clear the cache. - get_overlap(self, sd, deriv=None) : float - Return the overlap of the wavefunction with a Slater determinant. - generate_possible_exops(self, a_inds, c_inds): - Assign the excitation operators that can excite from the given indices to be annihilated - to the given indices to be created. - - """ - - def assign_ranks(self, ranks=None): - """Assign the ranks of the excitation operators. - - Parameters - ---------- - ranks : {int, list, None} - Ranks of the allowed excitation operators. - - Raises - ------ - TypeError - If 'ranks' is not None. - - """ - if ranks is not None: - raise TypeError("Only the default: ranks=[1, 2] is allowed") - else: - self.ranks = [1, 2] - - def assign_exops(self, indices=None): - """Assign the excitation operators that will be used to construct the CC operator. - - Parameters - ---------- - indices : None - The allowed excitation operators are solely defined by the occupied and virtual - orbitals of the given reference Slater determinant. - - Raises - ------ - TypeError - If `indices` is not None. - - Notes - ----- - The excitation operators are given as a list of lists of ints. - Each sub-list corresponds to an excitation operator. - In each sub-list, the first half of indices corresponds to the indices of the - spin-orbitals to annihilate, and the second half corresponds to the indices of the - spin-orbitals to create. - [a1, a2, ..., aN, c1, c2, ..., cN] - - """ - if indices is not None: - raise TypeError( - "Only the excitation operators constructed by default from " - "the given reference Slater determinant are allowed" - ) - else: - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - ex_to = [i for i in range(self.nspin) if i not in ex_from] - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt in ex_to: - exop = [occ_alpha, virt] - exops[tuple(exop)] = counter - counter += 1 - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt in ex_to: - exop = [occ_alpha + self.nspatial, virt] - exops[tuple(exop)] = counter - counter += 1 - super().assign_exops(None) - shift = len(exops) - for exop, ind in self.exops.items(): - exops[exop] = ind + shift - self.exops = exops diff --git a/fanpy/wfn/cc/apset1rog_d.py b/fanpy/wfn/cc/apset1rog_d.py index 99d47656..64f98ff2 100644 --- a/fanpy/wfn/cc/apset1rog_d.py +++ b/fanpy/wfn/cc/apset1rog_d.py @@ -5,18 +5,35 @@ class APset1roGD(PCCD): - r"""APset1roG wavefunction with only double excitations. + r"""APset1roG wavefunction with set-restricted double excitations. + + The reference wavefunction must be a seniority-0 Slater determinant. .. math:: - \left| {{\Psi }_{APset1roGD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{\begin{smallmatrix}a\in A \\ b\in B \\ A\bigcap B=\varnothing - \end{smallmatrix}}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\left| {{\Phi }_{0}} \right\rangle + \left| \Psi_{\mathrm{APset1roGD}} \right\rangle + = \prod_{\mu \in \mathcal{E}} + \left( 1 + t_\mu\, \tau_\mu \right) \left| \Phi_0 \right\rangle, + \qquad + \mathcal{E} = \left\{ \tau_{i\bar{i}}^{ab} + \mid a \in A,\, b \in B,\, A \cap B = \varnothing \right\}, + + where :math:`i, j, k, \ldots` index occupied spin orbitals and + :math:`a, b, c, \ldots` index virtual spin orbitals. The excitation + pool :math:`\mathcal{E}` consists of double excitations that annihilate + an occupied pair :math:`(i, \bar{i})` and create electrons into two + disjoint sets of virtual spin orbitals :math:`A` and :math:`B`. By + default, :math:`A` and :math:`B` correspond to the virtual alpha and + beta spin orbitals respectively, preserving spin projection symmetry; + custom disjoint sets may be supplied through the ``indices`` argument + of :meth:`assign_exops`. The ``s_type`` keyword has no effect on this + ansatz as the excitation pool contains only double excitations. + The reference wavefunction (``refwfn``) must be a seniority-0 Slater + determinant. + + See Ref. [1]_ for the full theoretical background. - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. Attributes ---------- @@ -98,6 +115,12 @@ class APset1roGD(PCCD): Assign the excitation operators that can excite from the given indices to be annihilated to the given indices to be created. + References + ---------- + .. [1] P. B. Gaikwad, T. D. Kim, M. Richer, R. A. Lokhande, G. Sánchez-Díaz; + P. A. Limacher, P. W. Ayers and R. A. Miranda-Quintana, "Coupled-cluster-inspired + geminal wavefunctions," *J. Chem. Phys.* **160**, 144108 (2024). + https://doi.org/10.1063/5.0196561 """ def assign_exops(self, indices=None): diff --git a/fanpy/wfn/cc/apset1rog_d_NEW.py b/fanpy/wfn/cc/apset1rog_d_NEW.py deleted file mode 100644 index f6466611..00000000 --- a/fanpy/wfn/cc/apset1rog_d_NEW.py +++ /dev/null @@ -1,180 +0,0 @@ -"""APset1roG wavefunction with only double excitations.""" - -from fanpy.tools import slater -from fanpy.wfn.cc.pccd_ap1rog_NEW import PCCD - - -class APset1roGD(PCCD): - r"""APset1roG wavefunction with only double excitations. - - .. math:: - - \left| {{\Psi }_{APset1roGD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{\begin{smallmatrix}a\in A \\ b\in B \\ A\bigcap B=\varnothing - \end{smallmatrix}}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\left| {{\Phi }_{0}} \right\rangle - - - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. - - Attributes - ---------- - nelec : int - Number of electrons. - nspin : int - Number of spin orbitals (alpha and beta). - dtype : {np.float64, np.complex128} - Data type of the wavefunction. - params : np.ndarray - Parameters of the wavefunction. - memory : float - Memory available for the wavefunction. - ranks : list of ints - Ranks of the excitation operators. - exops : list of list of int - Excitation operators given as lists of ints. The first half of indices correspond to - indices to be annihilated, the second half correspond to indices to be created. - refwfn : int - Reference wavefunction upon which the CC operator will act. - exop_combinations : dict - dictionary, the keys are tuples with the indices of annihilation and creation - operators, and the values are the excitation operators that allow to excite from the - annihilation to the creation operators. - - Properties - ---------- - nparams : int - Number of parameters. - nspatial : int - Number of spatial orbitals - param_shape : tuple of int - Shape of the parameters. - spin : int - Spin of the wavefunction. - seniority : int - Seniority of the wavefunction. - template_params : np.ndarray - Default parameters of the wavefunction. - nexops : int - Number of excitation operators. - nranks : int - Number of allowed ranks. - - Methods - ------- - __init__(self, nelec, nspin, memory=None, ngem=None, orbpairs=None, params=None) - Initialize the wavefunction. - assign_nelec(self, nelec) - Assign the number of electrons. - assign_nspin(self, nspin) - Assign the number of spin orbitals. - assign_dtype(self, dtype) - Assign the data type of the parameters. - assign_memory(self, memory=None) - Assign memory available for the wavefunction. - assign_ranks(self, ranks=None) - Assign the allowed excitation ranks. - assign_exops(self, exops=None) - Assign the allowed excitation operators. - assign_refwfn(self, refwfn=None) - Assign the reference wavefunction. - assign_params(self, params=None, add_noise=False) - Assign the parameters of the CC wavefunction. - get_ind(self, exop) : int - Return the parameter index that corresponds to a given excitation operator. - get_exop(self, ind) : list of int - Return the excitation operator that corresponds to a given parameter index. - product_amplitudes(self, inds, deriv=None) : float - Return the product of the CC amplitudes of the coefficients corresponding to - the given indices. - load_cache(self) - Load the functions whose values will be cached. - clear_cache(self) - Clear the cache. - get_overlap(self, sd, deriv=None) : float - Return the overlap of the wavefunction with a Slater determinant. - generate_possible_exops(self, a_inds, c_inds): - Assign the excitation operators that can excite from the given indices to be annihilated - to the given indices to be created. - - """ - - def assign_exops(self, indices=None): - """Assign the excitation operators that will be used to construct the CC operator. - - Parameters - ---------- - indices : {list of list of ints, None} - List of indices of the two disjoint sets of virtual orbitals to which one will excite - the occupied orbitals. - The default uses alpha/beta separation. - - Raises - ------ - TypeError - If `indices` is not a list of list of ints. - If `indices` is not None. - ValueError - If one tries to excite to occupied spin-orbitals. - If the lists of creation operators are not disjoint. - - Notes - ----- - The excitation operators are given as a list of lists of ints. - Each sub-list corresponds to an excitation operator. - In each sub-list, the first half of indices corresponds to the indices of the - spin-orbitals to annihilate, and the second half corresponds to the indices of the - spin-orbitals to create. - In previous assign_exops methods if one provides a non-default option for indices the first - sublist corresponds to indices of annihilation operators and the second sublist to indices - of creation operators, and we do not check if there are common elements between these - sublists. In this case both sublists correspond to indices of creation operators, - and they must be disjoint. It is assumed that one will excite from occupied alpha - spin-orbitals to indices given in the first sublist, and from occupied beta - spin-orbitals to indices given in the second sublist. - Takes care of any repetition in the sublists, and sorts them before generating the - excitation operators. - - """ - if indices is None: - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - virt_alphas = [i for i in range(self.nspin) if (i not in ex_from) and slater.is_alpha(i, self.nspatial)] - virt_betas = [i for i in range(self.nspin) if (i not in ex_from) and not slater.is_alpha(i, self.nspatial)] - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt_alpha in virt_alphas: - for virt_beta in virt_betas: - exop = [occ_alpha, occ_alpha + self.nspatial, virt_alpha, virt_beta] - exops[tuple(exop)] = counter - counter += 1 - self.exops = exops - - elif isinstance(indices, list): - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - if len(indices) != 2: - raise TypeError("`indices` must have exactly 2 elements") - for inds in indices: - if not isinstance(inds, list): - raise TypeError("The elements of `indices` must be lists of non-negative ints") - elif not all(isinstance(ind, int) for ind in inds): - raise TypeError("The elements of `indices` must be lists of non-negative ints") - elif not all(ind >= 0 for ind in inds): - raise ValueError("All `indices` must be lists of non-negative ints") - if not set(ex_from).isdisjoint(inds): - raise ValueError("`indices` cannot correspond to occupied spin-orbitals") - if not set(indices[0]).isdisjoint(indices[1]): - raise ValueError("The sets of creation operators must be disjoint") - indices = [list(set(indices[0])), list(set(indices[1]))] - indices[0].sort() - indices[1].sort() - for occ_alpha in ex_from[: len(ex_from) // 2]: - for i in indices[0]: - for j in indices[1]: - exop = [occ_alpha, occ_alpha + self.nspatial, i, j] - exops[tuple(exop)] = counter - counter += 1 - self.exops = exops diff --git a/fanpy/wfn/cc/apset1rog_sd.py b/fanpy/wfn/cc/apset1rog_sd.py index 05a06d89..72134b1c 100644 --- a/fanpy/wfn/cc/apset1rog_sd.py +++ b/fanpy/wfn/cc/apset1rog_sd.py @@ -5,21 +5,38 @@ class APset1roGSD(APset1roGD): - r"""APset1roG wavefunction with single and double excitations. + r"""APset1roG wavefunction with set-restricted single and double excitations. - .. math:: + The wavefunction is parameterized as - \left| {{\Psi }_{APset1roGSD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{\begin{smallmatrix}a\in A \\b\in B\\A\bigcap B=\varnothing - \end{smallmatrix}}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\prod\limits_{i=1}^{N/2\;}{\left( 1+\sum\limits_{b\in B}^{{}}{{{t}_{\bar{i};b}} - \hat{\tau }_{i\bar{i}}^{ib}} \right)\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a\in A}^{{}}{{{t}_{i;a}}\hat{\tau }_{i\bar{i}}^{a\bar{i}}} - \right)}\left| {{\Phi }_{0}} \right\rangle } + .. math:: + \left| \Psi_{\mathrm{APset1roGSD}} \right\rangle + = \prod_{\mu \in \mathcal{E}} + \left( 1 + t_\mu\, \tau_\mu \right) \left| \Phi_0 \right\rangle, + \qquad + \mathcal{E} = \left\{ \tau_{i}^{a},\, \tau_{\bar{i}}^{b},\, + \tau_{i\bar{i}}^{ab} + \mid a \in A,\, b \in B,\, A \cap B = \varnothing \right\}, + + where :math:`i, j, k, \ldots` index occupied spin orbitals and + :math:`a, b, c, \ldots` index virtual spin orbitals. The excitation + pool :math:`\mathcal{E}` extends APset1roGD with set-restricted + spin-conserving single excitations :math:`\tau_{i}^{a}` and + :math:`\tau_{\bar{i}}^{b}` alongside the double excitations + :math:`\tau_{i\bar{i}}^{ab}` inherited from APset1roGD. The disjoint + sets :math:`A` and :math:`B` default to virtual alpha and beta spin + orbitals respectively, preserving spin projection symmetry; custom + disjoint sets may be supplied through the ``indices`` argument of + :meth:`assign_exops`. The single excitation manifold can be further + constrained through the ``s_type`` keyword inherited from the PCCD + base class. + + The reference wavefunction (``refwfn``) must be a seniority-0 Slater + determinant. + + See Ref. [1]_ for the full theoretical background. - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. Attributes ---------- @@ -44,6 +61,11 @@ class APset1roGSD(APset1roGD): dictionary, the keys are tuples with the indices of annihilation and creation operators, and the values are the excitation operators that allow to excite from the annihilation to the creation operators. + s_type : str + Option controlling the screening of single excitations during + overlap evaluation. Inherited from :class:`PCCD`. + See PCCD.assign_s_type for available options and their definitions. + Properties ---------- @@ -84,6 +106,8 @@ class APset1roGSD(APset1roGD): Assign the reference wavefunction. assign_params(self, params=None, add_noise=False) Assign the parameters of the CC wavefunction. + assign_s_type(self, s_type): + Assign the option of seniority-breaking condition used for single excitations. get_ind(self, exop) : int Return the parameter index that corresponds to a given excitation operator. get_exop(self, ind) : list of int @@ -101,6 +125,12 @@ class APset1roGSD(APset1roGD): Assign the excitation operators that can excite from the given indices to be annihilated to the given indices to be created. + References + ---------- + .. [1] P. B. Gaikwad, T. D. Kim, M. Richer, R. A. Lokhande, G. Sánchez-Díaz; + P. A. Limacher, P. W. Ayers and R. A. Miranda-Quintana, "Coupled-cluster-inspired + geminal wavefunctions," *J. Chem. Phys.* **160**, 144108 (2024). + https://doi.org/10.1063/5.0196561 """ def assign_ranks(self, ranks=None): diff --git a/fanpy/wfn/cc/apset1rog_sd_NEW.py b/fanpy/wfn/cc/apset1rog_sd_NEW.py deleted file mode 100644 index 7c675391..00000000 --- a/fanpy/wfn/cc/apset1rog_sd_NEW.py +++ /dev/null @@ -1,218 +0,0 @@ -"""APset1roG wavefunction with single and double excitations.""" - -from fanpy.tools import slater -from fanpy.wfn.cc.apset1rog_d_NEW import APset1roGD - - -class APset1roGSD(APset1roGD): - r"""APset1roG wavefunction with single and double excitations. - - .. math:: - - \left| {{\Psi }_{APset1roGSD}} \right\rangle =\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{\begin{smallmatrix}a\in A \\b\in B\\A\bigcap B=\varnothing - \end{smallmatrix}}^{{}}{{{t}_{i;ab}}\hat{\tau }_{i\bar{i}}^{ab}} - \right)}\prod\limits_{i=1}^{N/2\;}{\left( 1+\sum\limits_{b\in B}^{{}}{{{t}_{\bar{i};b}} - \hat{\tau }_{i\bar{i}}^{ib}} \right)\prod\limits_{i=1}^{N/2\;} - {\left( 1+\sum\limits_{a\in A}^{{}}{{{t}_{i;a}}\hat{\tau }_{i\bar{i}}^{a\bar{i}}} - \right)}\left| {{\Phi }_{0}} \right\rangle } - - - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. - - Attributes - ---------- - nelec : int - Number of electrons. - nspin : int - Number of spin orbitals (alpha and beta). - dtype : {np.float64, np.complex128} - Data type of the wavefunction. - params : np.ndarray - Parameters of the wavefunction. - memory : float - Memory available for the wavefunction. - ranks : list of ints - Ranks of the excitation operators. - exops : list of list of int - Excitation operators given as lists of ints. The first half of indices correspond to - indices to be annihilated, the second half correspond to indices to be created. - refwfn : int - Reference wavefunction upon which the CC operator will act. - exop_combinations : dict - dictionary, the keys are tuples with the indices of annihilation and creation - operators, and the values are the excitation operators that allow to excite from the - annihilation to the creation operators. - - Properties - ---------- - nparams : int - Number of parameters. - nspatial : int - Number of spatial orbitals - param_shape : tuple of int - Shape of the parameters. - spin : int - Spin of the wavefunction. - seniority : int - Seniority of the wavefunction. - template_params : np.ndarray - Default parameters of the wavefunction. - nexops : int - Number of excitation operators. - nranks : int - Number of allowed ranks. - - Methods - ------- - __init__(self, nelec, nspin, memory=None, ngem=None, orbpairs=None, params=None) - Initialize the wavefunction. - assign_nelec(self, nelec) - Assign the number of electrons. - assign_nspin(self, nspin) - Assign the number of spin orbitals. - assign_dtype(self, dtype) - Assign the data type of the parameters. - assign_memory(self, memory=None) - Assign memory available for the wavefunction. - assign_ranks(self, ranks=None) - Assign the allowed excitation ranks. - assign_exops(self, exops=None) - Assign the allowed excitation operators. - assign_refwfn(self, refwfn=None) - Assign the reference wavefunction. - assign_params(self, params=None, add_noise=False) - Assign the parameters of the CC wavefunction. - get_ind(self, exop) : int - Return the parameter index that corresponds to a given excitation operator. - get_exop(self, ind) : list of int - Return the excitation operator that corresponds to a given parameter index. - product_amplitudes(self, inds, deriv=None) : float - Return the product of the CC amplitudes of the coefficients corresponding to - the given indices. - load_cache(self) - Load the functions whose values will be cached. - clear_cache(self) - Clear the cache. - get_overlap(self, sd, deriv=None) : float - Return the overlap of the wavefunction with a Slater determinant. - generate_possible_exops(self, a_inds, c_inds): - Assign the excitation operators that can excite from the given indices to be annihilated - to the given indices to be created. - - """ - - def assign_ranks(self, ranks=None): - """Assign the ranks of the excitation operators. - - Parameters - ---------- - ranks : {int, list, None} - Ranks of the allowed excitation operators. - - Raises - ------ - TypeError - If 'ranks' is not None. - - """ - if ranks is not None: - raise TypeError("Only the default: ranks=[1, 2] is allowed") - else: - self.ranks = [1, 2] - - def assign_exops(self, indices=None): - """Assign the excitation operators that will be used to construct the CC operator. - - Parameters - ---------- - indices : {list of list of ints, None} - List of indices of the two disjoint sets of virtual orbitals to which one will excite - the occupied orbitals. - The default uses alpha/beta separation. - - Raises - ------ - TypeError - If `indices` is not a list of list of ints. - If `indices` is not None. - ValueError - If one tries to excite to occupied spin-orbitals. - If the lists of creation operators are not disjoint. - - Notes - ----- - The excitation operators are given as a list of lists of ints. - Each sub-list corresponds to an excitation operator. - In each sub-list, the first half of indices corresponds to the indices of the - spin-orbitals to annihilate, and the second half corresponds to the indices of the - spin-orbitals to create. - In previous assign_exops methods if one provides a non-default option for indices the first - sublist corresponds to indices of annihilation operators and the second sublist to indices - of creation operators, and we do not check if there are common elements between these - sublists. In this case both sublists correspond to indices of creation operators, - and they must be disjoint. It is assumed that one will excite from occupied alpha - spin-orbitals to indices given in the first sublist, and from occupied beta - spin-orbitals to indices given in the second sublist. - Takes care of any repetition in the sublists, and sorts them before generating the - excitation operators. - - """ - if indices is None: - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - virt_alphas = [i for i in range(self.nspin) if (i not in ex_from) and slater.is_alpha(i, self.nspatial)] - virt_betas = [i for i in range(self.nspin) if (i not in ex_from) and not slater.is_alpha(i, self.nspatial)] - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt_alpha in virt_alphas: - exop = [occ_alpha, virt_alpha] - exops[tuple(exop)] = counter - counter += 1 - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt_beta in virt_betas: - exop = [occ_alpha + self.nspatial, virt_beta] - exops[tuple(exop)] = counter - counter += 1 - super().assign_exops(None) - shift = len(exops) - for exop, ind in self.exops.items(): - exops[exop] = ind + shift - self.exops = exops - - elif isinstance(indices, list): - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - if len(indices) != 2: - raise TypeError("`indices` must have exactly 2 elements") - for inds in indices: - if not isinstance(inds, list): - raise TypeError("The elements of `indices` must be lists of non-negative ints") - elif not all(isinstance(ind, int) for ind in inds): - raise TypeError("The elements of `indices` must be lists of non-negative ints") - elif not all(ind >= 0 for ind in inds): - raise ValueError("All `indices` must be lists of non-negative ints") - if not set(ex_from).isdisjoint(inds): - raise ValueError("`indices` cannot correspond to occupied spin-orbitals") - if not set(indices[0]).isdisjoint(indices[1]): - raise ValueError("The sets of annihilation operators must be disjoint") - indices = [list(set(indices[0])), list(set(indices[1]))] - indices[0].sort() - indices[1].sort() - for occ_alpha in ex_from[: len(ex_from) // 2]: - for i in indices[0]: - exop = [occ_alpha, i] - exops[tuple(exop)] = counter - counter += 1 - for occ_alpha in ex_from[: len(ex_from) // 2]: - for j in indices[1]: - exop = [occ_alpha + self.nspatial, j] - exops[tuple(exop)] = counter - counter += 1 - super().assign_exops(None) - shift = len(exops) - for exop, ind in self.exops.items(): - exops[exop] = ind + shift - self.exops = exops diff --git a/fanpy/wfn/cc/pccd_ap1rog.py b/fanpy/wfn/cc/pccd_ap1rog.py index 23f4d25b..0f37949a 100644 --- a/fanpy/wfn/cc/pccd_ap1rog.py +++ b/fanpy/wfn/cc/pccd_ap1rog.py @@ -1,4 +1,4 @@ -"""Pair Coupled Cluster Doubles-AP1roG wavefunction.""" +"""Pair Coupled Cluster Doubles (pCCD)/AP1roG wavefunction.""" import numpy as np from fanpy.tools import slater @@ -7,20 +7,30 @@ class PCCD(BaseCC): - r"""Pair CC doubles-AP1roG wavefunction. + r"""Pair coupled-cluster doubles (pCCD/AP1roG) wavefunction. + + The wavefunction is parameterized as .. math:: - \begin{align} - & \left| {{\Psi }_{AP1roG}} \right\rangle =\prod\limits_{i} - {\left( a_{i}^{\dagger }a_{{\bar{i}}}^{\dagger }+ - \sum\nolimits_{a}{t_{i}^{a}a_{a}^{\dagger }a_{{\bar{a}}}^{\dagger }} - \right)}\left| \theta \right\rangle \\ & =\prod\limits_{i} - {\left( 1+\sum\nolimits_{a}{t_{i}^{a}a_{a}^{\dagger }a_{{\bar{a}}}^{\dagger }{{a}_ - {{\bar{i}}}}{{a}_{i}}} \right)}\left| {{\Phi }_{0}} \right\rangle \\ \end{align} + \left| \Psi_{\mathrm{pCCD}} \right\rangle + = \prod_{\mu \in \mathcal{E}} + \left( 1 + t_\mu\, \tau_\mu \right) \left| \Phi_0 \right\rangle, + \qquad + \mathcal{E} = \left\{ \tau_{i\bar{i}}^{a\bar{a}} \right\}, + + where :math:`i, j, k, \ldots` index occupied spin orbitals and + :math:`a, b, c, \ldots` index virtual spin orbitals. The excitation + pool :math:`\mathcal{E}` consists exclusively of pair-preserving double + excitations that annihilate an occupied pair :math:`(i, \bar{i})` and + create electrons in a virtual pair :math:`(a, \bar{a})`, equivalent to + the AP1roG geminal formalism via the one-reference orbital construction + of Limacher et al.[1_] + + The reference wavefunction (``refwfn``) must be a seniority-0 Slater + determinant. - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. + See Ref. [2]_ for the full theoretical background. Attributes ---------- @@ -45,6 +55,17 @@ class PCCD(BaseCC): dictionary, the keys are tuples with the indices of annihilation and creation operators, and the values are the excitation operators that allow to excite from the annihilation to the creation operators. + s_type : str + Option indicating how the single excitations are allowed to affect the seniority + of the occupied and/or virtual spaces. + free: no restrictions imposed. + sen-o: single excitations must break an occupied pair. + (only excite from spin-orbital i if its complement is occupied) + sen-v: single excitations cannot form a pair. + (only excite to a virtual spin-orbital a if its complement is empty) + sen-ov: single excitations must break an occupied pair and cannot form extra pairs + (only excite from spin-orbital i if its complement is occupied and only excite to a virtual spin-orbital a if its complement is empty) + Properties ---------- @@ -85,6 +106,8 @@ class PCCD(BaseCC): Assign the reference wavefunction. assign_params(self, params=None, add_noise=False) Assign the parameters of the CC wavefunction. + assign_s_type(self, s_type): + Assign the option of seniority-breaking condition used for single excitations. get_ind(self, exop) : int Return the parameter index that corresponds to a given excitation operator. get_exop(self, ind) : list of int @@ -101,7 +124,21 @@ class PCCD(BaseCC): generate_possible_exops(self, a_inds, c_inds): Assign the excitation operators that can excite from the given indices to be annihilated to the given indices to be created. + _olp_double_derivative(self, sd): int + Returns double derivative of overlap for a given sd + References + ---------- + .. [1] P. A. Limacher, P. W. Ayers, P. A. Johnson, S. De Baerdemacker, D. Van Neck, and + P. Bultinck, "A new mean-field method suitable for strongly correlated electrons: + Computationally facile antisymmetric products of nonorthogonal geminals," + *J. Chem. Theory Comput.*, **9**, 1394-1401 (2013). + https://doi.org/10.1021/ct300902c + + .. [2] P. B. Gaikwad, T. D. Kim, M. Richer, R. A. Lokhande, G. Sánchez-Díaz, + P. A. Limacher, P. W. Ayers and R. A. Miranda-Quintana, "Coupled-cluster-inspired + geminal wavefunctions," *J. Chem. Phys.* **160**, 144108 (2024). + https://doi.org/10.1063/5.0196561 """ def __init__( @@ -115,6 +152,7 @@ def __init__( params=None, exop_combinations=None, refresh_exops=None, + s_type="sen-o", ): """Initialize the wavefunction. @@ -152,6 +190,7 @@ def __init__( super().__init__( nelec, nspin, memory=memory, params=params, exop_combinations=exop_combinations, refresh_exops=refresh_exops ) + self.assign_s_type(s_type=s_type) self.assign_ranks(ranks=ranks) self.assign_exops(indices=indices) self.assign_refwfn(refwfn=refwfn) @@ -224,6 +263,10 @@ def assign_exops(self, indices=None): spin-orbitals to create. [a1, a2, ..., aN, c1, c2, ..., cN] + # For pCCD: Pair excitations only; no single excitations are generated. + # Consequently, s_type-based seniority filtering is inactive for + # the standard PCCD excitation space. + """ if indices is not None: raise TypeError( @@ -279,6 +322,42 @@ def assign_refwfn(self, refwfn=None): # TODO: check that refwfn has the right number of spin-orbs self.refwfn = refwfn + def assign_s_type(self, s_type): + """Assign the seniority option for single excitations. + Notes + ----- + For doubles-only wavefunctions, these seniority restrictions + do not affect the overlap evaluation because the excitation operator + space contains only pair-double excitations and no singles. + The restrictions become relevant only when single excitations are + included in the excitation operator space. + + Parameters + ---------- + s_type : str + Option indicating how the single excitations are allowed to affect the seniority + of the occupied and/or virtual spaces. + free: no restrictions imposed. + sen-o: single excitations must break an occupied pair. + (only excite from spin-orbital i if its complement is occupied) + sen-v: single excitations cannot form a pair. + (only excite to a virtual spin-orbital a if its complement is empty) + sen-ov: single excitations must break an occupied pair and cannot form extra pairs. + (only excite from spin-orbital i if its complement is occupied AND + only excite to a virtual spin-orbital a if its complement is empty) + Raises + ------ + ValueError + If s_type is none of 'free', 'sen-o', 'sen-v', 'sen-ov'. + + """ + allowed = ["free", "sen-o", "sen-v", "sen-ov"] + + if s_type not in allowed: + raise ValueError("Invalid s_type") + + self.s_type = s_type + def _olp(self, sd): r"""Calculate the matrix element of the CC operator between the Slater determinants. @@ -325,8 +404,262 @@ def temp_olp(sd1, sd2): indices_multi = self.exop_combinations[tuple(a_inds + c_inds)] # FIXME: filter out rows whose excitation operators does not have annihilator that is # doubly occupied + occ_indices = set(slater.occ_indices(sd2)) + + # ------------------------------------------------------------------------------ + # NOTE: + # The seniority filtering options (s_type = "sen-o", "sen-v", "sen-ov") + # only affect excitation operators that contain single excitations + # (i.e. len(exop) == 2). + # + # Pure PCCD/AP1roG defines only pair-double excitation operators: + # [i, i_bar, a, a_bar] + # which have len(exop) == 4 and automatically preserve seniority. + # + # Therefore, for standard PCCD wavefunctions, the s_type logic is + # effectively dormant unless single excitations are introduced into + # the excitation operator space. + # ------------------------------------------------------------------------------ + + # "Break occupied pairs" + if self.s_type == "sen-o": + for exc_order in indices_multi: + indices_sign = indices_multi[exc_order] + selected_rows = [] + for row_ind, row in enumerate(indices_sign): + # occupied orbitals but have its spin composite also occupied + # AND its composite CANNOT participate in other excitation operators + trash = set([]) + skip_row = False + for exop in (self.ind_exops[i] for i in row[:-1]): + if len(exop) == 2: + if exop[0] in trash: + # skip + skip_row = True + break + if exop[0] < self.nspatial: + if exop[0] + self.nspatial not in occ_indices: + # skip + skip_row = True + break + else: + # add to trash + trash.add(exop[0]) + trash.add(exop[0] + self.nspatial) + if exop[0] >= self.nspatial: + if exop[0] - self.nspatial not in occ_indices: + # skip + skip_row = True + break + else: + # add to trash + trash.add(exop[0]) + trash.add(exop[0] - self.nspatial) + # FIXME: not sure + else: + for j in exop[: len(exop) // 2]: + if j in trash: + # skip + skip_row = True + break + else: + trash.add(j) + if skip_row: + break + + if not skip_row: + selected_rows.append(row_ind) + indices_multi[exc_order] = indices_sign[selected_rows] + # FIXME: + # print(selected_rows) + # print(indices_multi[exc_order]) + # "Don't form virtual pairs" + elif self.s_type == "sen-v": + for exc_order in indices_multi: + indices_sign = indices_multi[exc_order] + selected_rows = [] + for row_ind, row in enumerate(indices_sign): + # occupied orbitals but have its spin composite also occupied + # AND its composite CANNOT participate in other excitation operators + trash = set([]) + skip_row = False + for exop in (self.ind_exops[i] for i in row[:-1]): + if len(exop) == 2: + if exop[1] in trash: + # skip + skip_row = True + break + if exop[1] < self.nspatial: + if exop[1] + self.nspatial in occ_indices: + # skip + skip_row = True + break + else: + # add to trash + trash.add(exop[1]) + trash.add(exop[1] + self.nspatial) + if exop[1] >= self.nspatial: + if exop[1] - self.nspatial in occ_indices: + # skip + skip_row = True + break + else: + # add to trash + trash.add(exop[1]) + trash.add(exop[1] - self.nspatial) + # FIXME: not sure + else: + for j in exop[len(exop) // 2 :]: + if j in trash: + # skip + skip_row = True + break + else: + trash.add(j) + if skip_row: + break + + if not skip_row: + selected_rows.append(row_ind) + + indices_multi[exc_order] = indices_sign[selected_rows] + # FIXME: + # print(selected_rows) + # print(indices_multi[exc_order]) + # "Break occupied pairs AND don't form virtual pairs" + elif self.s_type == "sen-ov": + for exc_order in indices_multi: + indices_sign = indices_multi[exc_order] + selected_rows = [] + for row_ind, row in enumerate(indices_sign): + # occupied orbitals but have its spin composite also occupied + # AND its composite CANNOT participate in other excitation operators + trash = set([]) + skip_row = False + for exop in (self.ind_exops[i] for i in row[:-1]): + if len(exop) == 2: + if exop[0] in trash or exop[1] in trash: + # skip + skip_row = True + break + if exop[0] < self.nspatial: + if exop[0] + self.nspatial not in occ_indices: + # skip + skip_row = True + break + else: + # add to trash + trash.add(exop[0]) + trash.add(exop[0] + self.nspatial) + if exop[0] >= self.nspatial: + if exop[0] - self.nspatial not in occ_indices: + # skip + skip_row = True + break + else: + # add to trash + trash.add(exop[0]) + trash.add(exop[0] - self.nspatial) + if exop[1] < self.nspatial: + if exop[1] + self.nspatial in occ_indices: + # skip + skip_row = True + break + else: + # add to trash + trash.add(exop[1]) + trash.add(exop[1] + self.nspatial) + if exop[1] >= self.nspatial: + if exop[1] - self.nspatial in occ_indices: + # skip + skip_row = True + break + else: + # add to trash + trash.add(exop[1]) + trash.add(exop[1] - self.nspatial) + # FIXME: not sure + else: + for j in exop: + if j in trash: + # skip + skip_row = True + break + else: + trash.add(j) + if skip_row: + break + + if not skip_row: + selected_rows.append(row_ind) + + indices_multi[exc_order] = indices_sign[selected_rows] + # FIXME: + # print(selected_rows) + # print(indices_multi[exc_order]) + + amplitudes = self.product_amplitudes_multi(indices_multi) + # print(amplitudes) + val = sign * amplitudes + return val + + if isinstance(self.refwfn, CIWavefunction): + val = 0 + for refsd in self.refwfn.sd_vec: + val += temp_olp(sd, refsd) * self.refwfn.get_overlap(refsd) + return val + else: + return temp_olp(sd, self.refwfn) + + def _olp_deriv(self, sd): + """Calculate the derivative of the overlap with the Slater determinant. + + Parameters + ---------- + sd1 : int + Occupation vector of the left Slater determinant given as a bitstring. + sd2 : int + Occupation vector of the right Slater determinant given as a bitstring. + + Returns + ------- + olp : {float, complex} + Derivative of the overlap with respect to the given parameter. + + """ + + def temp_olp(sd1, sd2): + if sd1 == sd2: + return np.zeros(self.nparams) + # FIXME: this should definitely be vectorized + c_inds, a_inds = slater.diff_orbs(sd1, sd2) + if isinstance(a_inds, np.ndarray): + a_inds = a_inds.tolist() + if isinstance(c_inds, np.ndarray): + c_inds = c_inds.tolist() + # NOTE: Indices of the annihilation (a_inds) and creation (c_inds) operators + # that need to be applied to sd2 to turn it into sd1 + + # get sign + sign = slater.sign_excite(sd2, a_inds, c_inds) + + val = np.zeros(self.nparams) + if tuple(a_inds + c_inds) not in self.exop_combinations: + self.generate_possible_exops(a_inds, c_inds) + + # FIXME: sometimes exop contains virtual orbitals in annihilators may need to explicitly + # excite + indices_multi = self.exop_combinations[tuple(a_inds + c_inds)] + # FIXME: filter out rows whose excitation operators do not have annihilator that is + # doubly occupied + + # NOTE: + # Retaining the original pccd_ap1rog_NEW implementation using + # `occ_indices(sd)` here (rather than `sd2`). + # The correct determinant to use for occupancy filtering during + # excitation validation requires further investigation occ_indices = set(slater.occ_indices(sd)) - # print(indices_multi) + for exc_order in indices_multi: indices_sign = indices_multi[exc_order] selected_rows = [] @@ -337,24 +670,30 @@ def temp_olp(sd1, sd2): skip_row = False for exop in (self.ind_exops[i] for i in row[:-1]): if len(exop) == 2: + # skip because annihilator was used in a single excitation as the + # opposite spin if exop[0] in trash: # skip skip_row = True break if exop[0] < self.nspatial: + # skip because corresponding beta orbital is not occupied if exop[0] + self.nspatial not in occ_indices: # skip skip_row = True break + # this annihilator and its beta component cannot be used again else: # add to trash trash.add(exop[0]) trash.add(exop[0] + self.nspatial) if exop[0] >= self.nspatial: + # skip because corresponding alpha orbital is not occupied if exop[0] - self.nspatial not in occ_indices: # skip skip_row = True break + # this annihilator and its alpha component cannot be used again else: # add to trash trash.add(exop[0]) @@ -362,10 +701,14 @@ def temp_olp(sd1, sd2): # FIXME: not sure else: for j in exop[: len(exop) // 2]: + # skip because annihilator was used before as part of a single + # excitation or its opposite spin component if j in trash: # skip skip_row = True break + # not necessary because same orbital is not annihilated multiple + # times by construction else: trash.add(j) if skip_row: @@ -379,39 +722,42 @@ def temp_olp(sd1, sd2): # print(selected_rows) # print(indices_multi[exc_order]) - amplitudes = self.product_amplitudes_multi(indices_multi) + amplitudes = self.product_amplitudes_multi(indices_multi, deriv=True) # print(amplitudes) val = sign * amplitudes return val if isinstance(self.refwfn, CIWavefunction): - val = 0 + val = np.zeros(self.nparams) for refsd in self.refwfn.sd_vec: val += temp_olp(sd, refsd) * self.refwfn.get_overlap(refsd) return val else: return temp_olp(sd, self.refwfn) - def _olp_deriv(self, sd): - """Calculate the derivative of the overlap with the Slater determinant. + + def _olp_double_derivative(self, sd): + r"""Calculate the double derivative of the overlap with the Slater determinant. + + .. math:: + + H_{ab} = \frac{\partial^2}{\partial t_a \partial t_b} + \left[ \langle \mathrm{SD} \mid \Psi_{\mathrm{CC}} \rangle \right] Parameters ---------- - sd1 : int + sd : int Occupation vector of the left Slater determinant given as a bitstring. - sd2 : int - Occupation vector of the right Slater determinant given as a bitstring. - + Returns ------- - olp : {float, complex} - Derivative of the overlap with respect to the given parameter. + olp double derivative: 2D numpy array with Hessian elements :math:`H_{ab}`. """ - def temp_olp(sd1, sd2): + def temp_olp_double_deriv(sd1, sd2): if sd1 == sd2: - return np.zeros(self.nparams) + return np.zeros((self.nparams, self.nparams)) # FIXME: this should definitely be vectorized c_inds, a_inds = slater.diff_orbs(sd1, sd2) if isinstance(a_inds, np.ndarray): @@ -424,7 +770,7 @@ def temp_olp(sd1, sd2): # get sign sign = slater.sign_excite(sd2, a_inds, c_inds) - val = np.zeros(self.nparams) + val = np.zeros((self.nparams, self.nparams)) if tuple(a_inds + c_inds) not in self.exop_combinations: self.generate_possible_exops(a_inds, c_inds) @@ -433,7 +779,15 @@ def temp_olp(sd1, sd2): indices_multi = self.exop_combinations[tuple(a_inds + c_inds)] # FIXME: filter out rows whose excitation operators do not have annihilator that is # doubly occupied - occ_indices = set(slater.occ_indices(sd2)) + + # NOTE: + # Retaining the original pccd_ap1rog_NEW implementation + # of the whole _olp_double_deriv function alogn with + # `occ_indices(sd)` line. + # The correct determinant to use (`sd` or `sd2`) for occupancy filtering during + # excitation validation requires further investigation + occ_indices = set(slater.occ_indices(sd)) + for exc_order in indices_multi: indices_sign = indices_multi[exc_order] selected_rows = [] @@ -492,19 +846,16 @@ def temp_olp(sd1, sd2): selected_rows.append(row_ind) indices_multi[exc_order] = indices_sign[selected_rows] - # FIXME: - # print(selected_rows) - # print(indices_multi[exc_order]) - amplitudes = self.product_amplitudes_multi(indices_multi, deriv=True) - # print(amplitudes) - val = sign * amplitudes + hessian = self.product_amplitudes_multi_double_derivative(indices_multi) + val = sign * hessian return val if isinstance(self.refwfn, CIWavefunction): - val = np.zeros(self.nparams) + val = np.zeros(self.nparams, self.nparams) for refsd in self.refwfn.sd_vec: - val += temp_olp(sd, refsd) * self.refwfn.get_overlap(refsd) + val += temp_olp_double_deriv(sd, refsd) * self.refwfn.get_overlap(refsd) return val else: - return temp_olp(sd, self.refwfn) + return temp_olp_double_deriv(sd, self.refwfn) + diff --git a/fanpy/wfn/cc/pccd_ap1rog_NEW.py b/fanpy/wfn/cc/pccd_ap1rog_NEW.py deleted file mode 100644 index 5c46438f..00000000 --- a/fanpy/wfn/cc/pccd_ap1rog_NEW.py +++ /dev/null @@ -1,789 +0,0 @@ -"""Pair Coupled Cluster Doubles-AP1roG wavefunction.""" - -import numpy as np -from fanpy.tools import slater -from fanpy.wfn.cc.base import BaseCC -from fanpy.wfn.ci.base import CIWavefunction - - -class PCCD(BaseCC): - r"""Pair CC doubles-AP1roG wavefunction. - - .. math:: - - \begin{align} - & \left| {{\Psi }_{AP1roG}} \right\rangle =\prod\limits_{i} - {\left( a_{i}^{\dagger }a_{{\bar{i}}}^{\dagger }+ - \sum\nolimits_{a}{t_{i}^{a}a_{a}^{\dagger }a_{{\bar{a}}}^{\dagger }} - \right)}\left| \theta \right\rangle \\ & =\prod\limits_{i} - {\left( 1+\sum\nolimits_{a}{t_{i}^{a}a_{a}^{\dagger }a_{{\bar{a}}}^{\dagger }{{a}_ - {{\bar{i}}}}{{a}_{i}}} \right)}\left| {{\Phi }_{0}} \right\rangle \\ \end{align} - - In this case the reference wavefunction can only be a single Slater determinant with - seniority 0. - - Attributes - ---------- - nelec : int - Number of electrons. - nspin : int - Number of spin orbitals (alpha and beta). - dtype : {np.float64, np.complex128} - Data type of the wavefunction. - params : np.ndarray - Parameters of the wavefunction. - memory : float - Memory available for the wavefunction. - ranks : list of ints - Ranks of the excitation operators. - exops : list of list of int - Excitation operators given as lists of ints. The first half of indices correspond to - indices to be annihilated, the second half correspond to indices to be created. - refwfn : int - Reference wavefunction upon which the CC operator will act. - exop_combinations : dict - dictionary, the keys are tuples with the indices of annihilation and creation - operators, and the values are the excitation operators that allow to excite from the - annihilation to the creation operators. - s_type : str - Option indicating how the single excitations are allowed to affect the seniority - of the occupied and/or virtual spaces. - free: no restrictions imposed. - sen-o: single excitations must break an occupied pair. - (only excite from spin-orbital i if its complement is occupied) - sen-v: single excitations cannot form a pair. - (only excite to a virtual spin-orbital a if its complement is empty) - sen-ov: single excitations must break an occupied pair and cannot form extra pairs - - Properties - ---------- - nparams : int - Number of parameters. - nspatial : int - Number of spatial orbitals - param_shape : tuple of int - Shape of the parameters. - spin : int - Spin of the wavefunction. - seniority : int - Seniority of the wavefunction. - template_params : np.ndarray - Default parameters of the wavefunction. - nexops : int - Number of excitation operators. - nranks : int - Number of allowed ranks. - - Methods - ------- - __init__(self, nelec, nspin, memory=None, ngem=None, orbpairs=None, params=None) - Initialize the wavefunction. - assign_nelec(self, nelec) - Assign the number of electrons. - assign_nspin(self, nspin) - Assign the number of spin orbitals. - assign_dtype(self, dtype) - Assign the data type of the parameters. - assign_memory(self, memory=None) - Assign memory available for the wavefunction. - assign_ranks(self, ranks=None) - Assign the allowed excitation ranks. - assign_exops(self, exops=None) - Assign the allowed excitation operators. - assign_refwfn(self, refwfn=None) - Assign the reference wavefunction. - assign_params(self, params=None, add_noise=False) - Assign the parameters of the CC wavefunction. - get_ind(self, exop) : int - Return the parameter index that corresponds to a given excitation operator. - get_exop(self, ind) : list of int - Return the excitation operator that corresponds to a given parameter index. - product_amplitudes(self, inds, deriv=None) : float - Return the product of the CC amplitudes of the coefficients corresponding to - the given indices. - load_cache(self) - Load the functions whose values will be cached. - clear_cache(self) - Clear the cache. - get_overlap(self, sd, deriv=None) : float - Return the overlap of the wavefunction with a Slater determinant. - generate_possible_exops(self, a_inds, c_inds): - Assign the excitation operators that can excite from the given indices to be annihilated - to the given indices to be created. - assign_s_type(self, s_type): - Assign the s_type option. - _olp_double_derivative(self, sd): int - Returns double derivative of overlap for a given sd - """ - - def __init__( - self, - nelec, - nspin, - memory=None, - ranks=None, - indices=None, - refwfn=None, - params=None, - exop_combinations=None, - refresh_exops=None, - s_type="sen-o", - ): - """Initialize the wavefunction. - - Parameters - ---------- - nelec : int - Number of electrons. - nspin : int - Number of spin orbitals. - dtype : {float, complex, np.float64, np.complex128, None} - Numpy data type. - Default is `np.float64`. - memory : {float, int, str, None} - Memory available for the wavefunction. - Default does not limit memory usage (i.e. infinite). - ranks : {int, list of int, None} - Ranks of the excitation operators (in increasing order). - If int is provided, it takes all the ranks lower than that. - Default is None, which is equivalent to taking ranks=self.nelec. - indices : {list of list of ints, None} - List of lists containing the indices ot the spin-orbitals to annihilate and create. - The first sub-list contains indices of orbitals to annihilate. - The second sub-list contains indices of orbitals to create. - Default generates all possible indices according to the given ranks. - refwfn: {int, None} - Reference wavefunction upon which the CC operator will act. - params : {np.ndarray, BaseCC, None} - 1-vector of CC amplitudes. - exop_combinations : dict - dictionary, the keys are tuples with the indices of annihilation and creation - operators, and the values are the excitation operators that allow to excite from the - annihilation to the creation operators. - - """ - super().__init__( - nelec, nspin, memory=memory, params=params, exop_combinations=exop_combinations, refresh_exops=refresh_exops - ) - self.assign_s_type(s_type=s_type) - self.assign_ranks(ranks=ranks) - self.assign_exops(indices=indices) - self.assign_refwfn(refwfn=refwfn) - - # mapping from indices to exops - self.ind_exops = sorted(self.exops.keys(), key=lambda x: self.exops[x]) - - def assign_nelec(self, nelec): - """Assign the number of electrons. - - Parameters - ---------- - nelec : int - Number of electrons. - - Raises - ------ - TypeError - If number of electrons is not an integer. - ValueError - If number of electrons is not a positive number. - If number of electrons is odd. - - """ - super().assign_nelec(nelec) - if self.nelec % 2 != 0: - raise ValueError("Odd number of electrons is not supported") - - def assign_ranks(self, ranks=None): - """Assign the ranks of the excitation operators. - - Parameters - ---------- - ranks : None - Ranks of the allowed excitation operators. Set by default to [2]. - - Raises - ------ - ValueError - If a value that is not the default is provided. - If the maximum rank is greater than the number of electrons. - - """ - if ranks is not None: - raise ValueError("Only the default, rank = 2, is allowed") - if self.nelec <= 1: - raise ValueError("Only wavefunctions with more than 1 electron can be considered") - self.ranks = [2] - - def assign_exops(self, indices=None): - """Assign the excitation operators that will be used to construct the CC operator. - - Parameters - ---------- - indices : None - The allowed excitation operators are solely defined by the occupied and virtual - orbitals of the given reference Slater determinant. - - Raises - ------ - TypeError - If `indices` is not None. - - Notes - ----- - The excitation operators are given as a list of lists of ints. - Each sub-list corresponds to an excitation operator. - In each sub-list, the first half of indices corresponds to the indices of the - spin-orbitals to annihilate, and the second half corresponds to the indices of the - spin-orbitals to create. - [a1, a2, ..., aN, c1, c2, ..., cN] - - """ - if indices is not None: - raise TypeError( - "Only the excitation operators constructed by default from " - "the given reference Slater determinant are allowed" - ) - else: - exops = {} - counter = 0 - ex_from = slater.occ_indices(self.refwfn) - ex_to = [i for i in range(self.nspin) if i not in ex_from] - for occ_alpha in ex_from[: len(ex_from) // 2]: - for virt_alpha in ex_to[: len(ex_to) // 2]: - exop = [occ_alpha, occ_alpha + self.nspatial, virt_alpha, virt_alpha + self.nspatial] - exops[tuple(exop)] = counter - counter += 1 - self.exops = exops - - def assign_refwfn(self, refwfn=None): - """Assign the reference wavefunction upon which the CC operator will act. - - Parameters - ---------- - refwfn: {int, None} - Seniority 0 wavefunction that will be modified by the CC operator. - Default is the ground-state Slater determinant. - - Raises - ------ - TypeError - If refwfn is not a int instance. - ValueError - If refwfn does not have the right number of electrons. - If refwfn does not have the right number of spin orbitals. - If refwfn is not a seniority-0 wavefunction. - - """ - if refwfn is None: - self.refwfn = slater.ground(nocc=self.nelec, norbs=self.nspin) - else: - if not isinstance(refwfn, int): - raise TypeError("refwfn must be a int object") - if slater.total_occ(refwfn) != self.nelec: - raise ValueError("refwfn must have {} electrons".format(self.nelec)) - if not all( - [i + self.nspatial in slater.occ_indices(refwfn) for i in slater.occ_indices(refwfn)[: self.nspatial]] - + [i - self.nspatial in slater.occ_indices(refwfn) for i in slater.occ_indices(refwfn)[self.nspatial :]] - ): - raise ValueError("refwfn must be a seniority-0 wavefuntion") - # TODO: check that refwfn has the right number of spin-orbs - self.refwfn = refwfn - - def assign_s_type(self, s_type): - """Assign the seniority option for single excitations. - - Parameters - ---------- - s_type : str - Option indicating how the single excitations are allowed to affect the seniority - of the occupied and/or virtual spaces. - free: no restrictions imposed. - sen-o: single excitations must break an occupied pair. - (only excite from spin-orbital i if its complement is occupied) - sen-v: single excitations cannot form a pair. - (only excite to a virtual spin-orbital a if its complement is empty) - sen-ov: single excitations must break an occupied pair and cannot form extra pairs - - Raises - ------ - ValueError - If s_type is none of 'free', 'sen-o', 'sen-v', 'sen-ov'. - - """ - self.s_type = s_type - - def _olp(self, sd): - r"""Calculate the matrix element of the CC operator between the Slater determinants. - - .. math:: - - \[\left\langle {{m}_{1}} - \right|\prod\limits_{u,v}{\left( 1+t_{u}^{v}\hat{\tau }_{u}^{v} \right)} - \left| {{m}_{2}} \right\rangle \] - - Parameters - ---------- - sd1 : int - Occupation vector of the left Slater determinant given as a bitstring. - sd2 : int - Occupation vector of the right Slater determinant given as a bitstring. - - Returns - ------- - olp : {float, complex} - Matrix element of the CC operator between the given Slater determinant. - - """ - - def temp_olp(sd1, sd2): - if sd1 == sd2: - return 1.0 - c_inds, a_inds = slater.diff_orbs(sd1, sd2) - if isinstance(a_inds, np.ndarray): - a_inds = a_inds.tolist() - if isinstance(c_inds, np.ndarray): - c_inds = c_inds.tolist() - # NOTE: Indices of the annihilation (a_inds) and creation (c_inds) operators - # that need to be applied to sd2 to turn it into sd1 - - # get sign - sign = slater.sign_excite(sd2, a_inds, c_inds) - - val = 0.0 - if tuple(a_inds + c_inds) not in self.exop_combinations: - self.generate_possible_exops(a_inds, c_inds) - - # FIXME: sometimes exop contains virtual orbitals in annihilators may need to explicitly - # excite - indices_multi = self.exop_combinations[tuple(a_inds + c_inds)] - # print("indices_multi", indices_multi) - # FIXME: filter out rows whose excitation operators does not have annihilator that is - # doubly occupied - occ_indices = set(slater.occ_indices(sd2)) - # virt_indices = set(slater.virt_indices(sd, self.nspin)) - # print(indices_multi) - # "Break occupied pairs" - if self.s_type == "sen-o": - for exc_order in indices_multi: - indices_sign = indices_multi[exc_order] - selected_rows = [] - for row_ind, row in enumerate(indices_sign): - # occupied orbitals but have its spin composite also occupied - # AND its composite CANNOT participate in other excitation operators - trash = set([]) - skip_row = False - for exop in (self.ind_exops[i] for i in row[:-1]): - if len(exop) == 2: - if exop[0] in trash: - # skip - skip_row = True - break - if exop[0] < self.nspatial: - if exop[0] + self.nspatial not in occ_indices: - # skip - skip_row = True - break - else: - # add to trash - trash.add(exop[0]) - trash.add(exop[0] + self.nspatial) - if exop[0] >= self.nspatial: - if exop[0] - self.nspatial not in occ_indices: - # skip - skip_row = True - break - else: - # add to trash - trash.add(exop[0]) - trash.add(exop[0] - self.nspatial) - # FIXME: not sure - else: - for j in exop[: len(exop) // 2]: - if j in trash: - # skip - skip_row = True - break - else: - trash.add(j) - if skip_row: - break - - if not skip_row: - selected_rows.append(row_ind) - - indices_multi[exc_order] = indices_sign[selected_rows] - # FIXME: - # print(selected_rows) - # print(indices_multi[exc_order]) - # "Don't form virtual pairs" - elif self.s_type == "sen-v": - for exc_order in indices_multi: - indices_sign = indices_multi[exc_order] - selected_rows = [] - for row_ind, row in enumerate(indices_sign): - # occupied orbitals but have its spin composite also occupied - # AND its composite CANNOT participate in other excitation operators - trash = set([]) - skip_row = False - for exop in (self.ind_exops[i] for i in row[:-1]): - if len(exop) == 2: - if exop[1] in trash: - # skip - skip_row = True - break - if exop[1] < self.nspatial: - if exop[1] + self.nspatial in occ_indices: - # skip - skip_row = True - break - else: - # add to trash - trash.add(exop[1]) - trash.add(exop[1] + self.nspatial) - if exop[1] >= self.nspatial: - if exop[1] - self.nspatial in occ_indices: - # skip - skip_row = True - break - else: - # add to trash - trash.add(exop[1]) - trash.add(exop[1] - self.nspatial) - # FIXME: not sure - else: - for j in exop[len(exop) // 2 :]: - if j in trash: - # skip - skip_row = True - break - else: - trash.add(j) - if skip_row: - break - - if not skip_row: - selected_rows.append(row_ind) - - indices_multi[exc_order] = indices_sign[selected_rows] - # FIXME: - # print(selected_rows) - # print(indices_multi[exc_order]) - # "Break occupied pairs AND don't form virtual pairs" - elif self.s_type == "sen-ov": - for exc_order in indices_multi: - indices_sign = indices_multi[exc_order] - selected_rows = [] - for row_ind, row in enumerate(indices_sign): - # occupied orbitals but have its spin composite also occupied - # AND its composite CANNOT participate in other excitation operators - trash = set([]) - skip_row = False - for exop in (self.ind_exops[i] for i in row[:-1]): - if len(exop) == 2: - if exop[0] in trash or exop[1] in trash: - # skip - skip_row = True - break - if exop[0] < self.nspatial: - if exop[0] + self.nspatial not in occ_indices: - # skip - skip_row = True - break - else: - # add to trash - trash.add(exop[0]) - trash.add(exop[0] + self.nspatial) - if exop[0] >= self.nspatial: - if exop[0] - self.nspatial not in occ_indices: - # skip - skip_row = True - break - else: - # add to trash - trash.add(exop[0]) - trash.add(exop[0] - self.nspatial) - if exop[1] < self.nspatial: - if exop[1] + self.nspatial in occ_indices: - # skip - skip_row = True - break - else: - # add to trash - trash.add(exop[1]) - trash.add(exop[1] + self.nspatial) - if exop[1] >= self.nspatial: - if exop[1] - self.nspatial in occ_indices: - # skip - skip_row = True - break - else: - # add to trash - trash.add(exop[1]) - trash.add(exop[1] - self.nspatial) - # FIXME: not sure - else: - for j in exop: - if j in trash: - # skip - skip_row = True - break - else: - trash.add(j) - if skip_row: - break - - if not skip_row: - selected_rows.append(row_ind) - - indices_multi[exc_order] = indices_sign[selected_rows] - # FIXME: - # print(selected_rows) - # print(indices_multi[exc_order]) - - amplitudes = self.product_amplitudes_multi(indices_multi) - # print(amplitudes) - val = sign * amplitudes - return val - - if isinstance(self.refwfn, CIWavefunction): - val = 0 - for refsd in self.refwfn.sd_vec: - val += temp_olp(sd, refsd) * self.refwfn.get_overlap(refsd) - return val - else: - return temp_olp(sd, self.refwfn) - - def _olp_deriv(self, sd): - """Calculate the derivative of the overlap with the Slater determinant. - - Parameters - ---------- - sd1 : int - Occupation vector of the left Slater determinant given as a bitstring. - sd2 : int - Occupation vector of the right Slater determinant given as a bitstring. - - Returns - ------- - olp : {float, complex} - Derivative of the overlap with respect to the given parameter. - - """ - - def temp_olp(sd1, sd2): - if sd1 == sd2: - return np.zeros(self.nparams) - # FIXME: this should definitely be vectorized - c_inds, a_inds = slater.diff_orbs(sd1, sd2) - if isinstance(a_inds, np.ndarray): - a_inds = a_inds.tolist() - if isinstance(c_inds, np.ndarray): - c_inds = c_inds.tolist() - # NOTE: Indices of the annihilation (a_inds) and creation (c_inds) operators - # that need to be applied to sd2 to turn it into sd1 - - # get sign - sign = slater.sign_excite(sd2, a_inds, c_inds) - - val = np.zeros(self.nparams) - if tuple(a_inds + c_inds) not in self.exop_combinations: - self.generate_possible_exops(a_inds, c_inds) - - # FIXME: sometimes exop contains virtual orbitals in annihilators may need to explicitly - # excite - indices_multi = self.exop_combinations[tuple(a_inds + c_inds)] - # FIXME: filter out rows whose excitation operators do not have annihilator that is - # doubly occupied - occ_indices = set(slater.occ_indices(sd)) - for exc_order in indices_multi: - indices_sign = indices_multi[exc_order] - selected_rows = [] - for row_ind, row in enumerate(indices_sign): - # occupied orbitals but have its spin composite also occupied - # AND its composite CANNOT participate in other excitation operators - trash = set([]) - skip_row = False - for exop in (self.ind_exops[i] for i in row[:-1]): - if len(exop) == 2: - # skip because annihilator was used in a single excitation as the - # opposite spin - if exop[0] in trash: - # skip - skip_row = True - break - if exop[0] < self.nspatial: - # skip because corresponding beta orbital is not occupied - if exop[0] + self.nspatial not in occ_indices: - # skip - skip_row = True - break - # this annihilator and its beta component cannot be used again - else: - # add to trash - trash.add(exop[0]) - trash.add(exop[0] + self.nspatial) - if exop[0] >= self.nspatial: - # skip because corresponding alpha orbital is not occupied - if exop[0] - self.nspatial not in occ_indices: - # skip - skip_row = True - break - # this annihilator and its alpha component cannot be used again - else: - # add to trash - trash.add(exop[0]) - trash.add(exop[0] - self.nspatial) - # FIXME: not sure - else: - for j in exop[: len(exop) // 2]: - # skip because annihilator was used before as part of a single - # excitation or its opposite spin component - if j in trash: - # skip - skip_row = True - break - # not necessary because same orbital is not annihilated multiple - # times by construction - else: - trash.add(j) - if skip_row: - break - - if not skip_row: - selected_rows.append(row_ind) - - indices_multi[exc_order] = indices_sign[selected_rows] - # FIXME: - # print(selected_rows) - # print(indices_multi[exc_order]) - - amplitudes = self.product_amplitudes_multi(indices_multi, deriv=True) - # print(amplitudes) - val = sign * amplitudes - return val - - if isinstance(self.refwfn, CIWavefunction): - val = np.zeros(self.nparams) - for refsd in self.refwfn.sd_vec: - val += temp_olp(sd, refsd) * self.refwfn.get_overlap(refsd) - return val - else: - return temp_olp(sd, self.refwfn) - - - def _olp_double_derivative(self, sd): - r"""Calculate the double derivative of the overlap with the Slater determinant. - - .. math:: - - H_{ab} = \frac{\partial^2}{\partial t_a \partial t_b} - \left[ \langle \mathrm{SD} \mid \Psi_{\mathrm{CC}} \rangle \right] - - Parameters - ---------- - sd : int - Occupation vector of the left Slater determinant given as a bitstring. - - Returns - ------- - olp double derivative: 2D numpy array with Hessian elements :math:`H_{ab}`. - - """ - - def temp_olp_double_deriv(sd1, sd2): - if sd1 == sd2: - return np.zeros((self.nparams, self.nparams)) - # FIXME: this should definitely be vectorized - c_inds, a_inds = slater.diff_orbs(sd1, sd2) - if isinstance(a_inds, np.ndarray): - a_inds = a_inds.tolist() - if isinstance(c_inds, np.ndarray): - c_inds = c_inds.tolist() - # NOTE: Indices of the annihilation (a_inds) and creation (c_inds) operators - # that need to be applied to sd2 to turn it into sd1 - - # get sign - sign = slater.sign_excite(sd2, a_inds, c_inds) - - val = np.zeros((self.nparams, self.nparams)) - if tuple(a_inds + c_inds) not in self.exop_combinations: - self.generate_possible_exops(a_inds, c_inds) - - # FIXME: sometimes exop contains virtual orbitals in annihilators may need to explicitly - # excite - indices_multi = self.exop_combinations[tuple(a_inds + c_inds)] - # FIXME: filter out rows whose excitation operators do not have annihilator that is - # doubly occupied - occ_indices = set(slater.occ_indices(sd)) - for exc_order in indices_multi: - indices_sign = indices_multi[exc_order] - selected_rows = [] - for row_ind, row in enumerate(indices_sign): - # occupied orbitals but have its spin composite also occupied - # AND its composite CANNOT participate in other excitation operators - trash = set([]) - skip_row = False - for exop in (self.ind_exops[i] for i in row[:-1]): - if len(exop) == 2: - # skip because annihilator was used in a single excitation as the - # opposite spin - if exop[0] in trash: - # skip - skip_row = True - break - if exop[0] < self.nspatial: - # skip because corresponding beta orbital is not occupied - if exop[0] + self.nspatial not in occ_indices: - # skip - skip_row = True - break - # this annihilator and its beta component cannot be used again - else: - # add to trash - trash.add(exop[0]) - trash.add(exop[0] + self.nspatial) - if exop[0] >= self.nspatial: - # skip because corresponding alpha orbital is not occupied - if exop[0] - self.nspatial not in occ_indices: - # skip - skip_row = True - break - # this annihilator and its alpha component cannot be used again - else: - # add to trash - trash.add(exop[0]) - trash.add(exop[0] - self.nspatial) - # FIXME: not sure - else: - for j in exop[: len(exop) // 2]: - # skip because annihilator was used before as part of a single - # excitation or its opposite spin component - if j in trash: - # skip - skip_row = True - break - # not necessary because same orbital is not annihilated multiple - # times by construction - else: - trash.add(j) - if skip_row: - break - - if not skip_row: - selected_rows.append(row_ind) - - indices_multi[exc_order] = indices_sign[selected_rows] - - hessian = self.product_amplitudes_multi_double_derivative(indices_multi) - val = sign * hessian - return val - - if isinstance(self.refwfn, CIWavefunction): - val = np.zeros((self.nparams, self.nparams)) - for refsd in self.refwfn.sd_vec: - val += temp_olp_double_deriv(sd, refsd) * self.refwfn.get_overlap(refsd) - return val - else: - return temp_olp_double_deriv(sd, self.refwfn) - diff --git a/tests/test_wfn_cc_pccd_ap1rog.py b/tests/test_wfn_cc_pccd_ap1rog.py index 8ebe8206..66a295b2 100644 --- a/tests/test_wfn_cc_pccd_ap1rog.py +++ b/tests/test_wfn_cc_pccd_ap1rog.py @@ -1,5 +1,6 @@ """Test fanpy.wavefunction.cc.pccd_ap1rog.""" import pytest +import numpy as np from fanpy.tools import slater from fanpy.wfn.cc.pccd_ap1rog import PCCD @@ -64,3 +65,130 @@ def test_assign_refwfn(): test.assign_refwfn(0b11000011) test.assign_refwfn() assert test.refwfn == (0b00110011) + + +def test_assign_refwfn_sen_0_check(): + """Test seniority-0 check in PCCD.assign_refwfn.""" + test = TempPCCD() + test.assign_nelec(4) + test.assign_nspin(16) + with pytest.raises(ValueError): + test.assign_refwfn(0b11000011) + with pytest.raises(ValueError): + test.assign_refwfn(0b0000010100000011) + test.assign_refwfn() + assert test.refwfn == (0b0000001100000011) + + +def test_assign_s_type(): + """Test PCCD.assign_s_type.""" + test = TempPCCD() + + test.assign_s_type("free") + assert test.s_type == "free" + + test.assign_s_type("sen-o") + assert test.s_type == "sen-o" + + test.assign_s_type("sen-v") + assert test.s_type == "sen-v" + + test.assign_s_type("sen-ov") + assert test.s_type == "sen-ov" + + +def test_assign_s_type_invalid(): + """Test invalid s_type.""" + test = TempPCCD() + + with pytest.raises(ValueError): + test.assign_s_type("bad-option") + + +def test_init_s_type(): + """Test initialization of s_type.""" + test = PCCD(4, 8, s_type="sen-v") + assert test.s_type == "sen-v" + + +def tests_type_effct_on_pCCD_overlap(): + """Test s_type effect on pCCD overlap. + pCCD doesn't have singles so it shouldn't enter + into the sen-x logics of singles. + """ + + sd = 0b10100011 + + test_free = PCCD(4, 8, s_type="free") + test_seno = PCCD(4, 8, s_type="sen-o") + + olp_free = test_free.get_overlap(sd) + olp_seno = test_seno.get_overlap(sd) + + assert olp_free == olp_seno + + +def test_olp_double_derivative_shape(): + """Test overlap Hessian shape.""" + + test = PCCD(4, 8) + + hess = test._olp_double_derivative(test.refwfn) + + assert hess.shape == (test.nparams, test.nparams) + + +def test_olp_double_derivative_symmetric(): + """Test overlap Hessian symmetry.""" + + test = PCCD(4, 8) + + sd = 0b11001100 + + hess = test._olp_double_derivative(sd) + + assert np.allclose(hess, hess.T) + + +def test_olp_double_derivative_zero_diagonal(): + """Test Hessian diagonal vanishes.""" + + test = PCCD(4, 8) + + sd = 0b11001100 + + hess = test._olp_double_derivative(sd) + + assert np.allclose(np.diag(hess), 0) + + +def test_olp_double_derivative_finite_difference(): + """Test overlap Hessian with finite differences.""" + + test = PCCD(4, 8) + + sd = 0b11001100 + + h = 1e-7 + + analytic = test._olp_double_derivative(sd) + + numerical = np.zeros_like(analytic) + + orig = test.params.copy() + + for j in range(test.nparams): + + test.params = orig.copy() + test.params[j] += h + plus = test._olp_deriv(sd) + + test.params = orig.copy() + test.params[j] -= h + minus = test._olp_deriv(sd) + + numerical[:, j] = (plus - minus) / (2 * h) + + test.params = orig + + assert np.allclose(analytic, numerical, atol=1e-5)