From 947c44b43947706c9aa39e24a140d502aed2c2bc Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Tue, 12 Mar 2024 13:47:49 -0500 Subject: [PATCH 01/11] Add fanpt files FanPT files imported from FanCI --- .gitignore | 1 + pyci/fanci/fanpt_wrapper.py | 229 ++++++++++++++++++ pyci/fanpt/__init__.py | 17 ++ pyci/fanpt/base_fanpt_container.py | 298 +++++++++++++++++++++++ pyci/fanpt/fanpt_constant_terms.py | 167 +++++++++++++ pyci/fanpt/fanpt_cont_e_free.py | 259 ++++++++++++++++++++ pyci/fanpt/fanpt_cont_e_param.py | 236 ++++++++++++++++++ pyci/fanpt/fanpt_updater.py | 374 +++++++++++++++++++++++++++++ 8 files changed, 1581 insertions(+) create mode 100644 pyci/fanci/fanpt_wrapper.py create mode 100644 pyci/fanpt/__init__.py create mode 100644 pyci/fanpt/base_fanpt_container.py create mode 100644 pyci/fanpt/fanpt_constant_terms.py create mode 100644 pyci/fanpt/fanpt_cont_e_free.py create mode 100644 pyci/fanpt/fanpt_cont_e_param.py create mode 100644 pyci/fanpt/fanpt_updater.py diff --git a/.gitignore b/.gitignore index befbdc7..f43948b 100644 --- a/.gitignore +++ b/.gitignore @@ -94,3 +94,4 @@ deps/ # Jupyter Notebooks .ipynb_checkpoints +.DS_Store diff --git a/pyci/fanci/fanpt_wrapper.py b/pyci/fanci/fanpt_wrapper.py new file mode 100644 index 0000000..e285a76 --- /dev/null +++ b/pyci/fanci/fanpt_wrapper.py @@ -0,0 +1,229 @@ +""" FANPT wrapper""" + +import numpy as np +import pyci + +from .fanci import FanCI +from ..fanpt import FANPTUpdater, FANPTContainerEParam, FANPTContainerEFree + + +def solve_fanpt( + fanci_wfn, + ham0, + ham1, + guess_params, + fill, + energy_active=True, + resum=False, + ref_sd=0, + final_order=1, + lambda_i=0.0, + lambda_f=1.0, + steps=1, + kwargs=None, + solver_kwargs=None, +): + """[summary] + + Args: + fanci_wfn : FanCI class + FanCI wavefunction. + guess_params : np.ndarray + Initial guess for wave function parameters. + ham0 : pyci.hamiltonian + PyCI Hamiltonian of the ideal system. + ham1 : pyci.hamiltonian + PyCI Hamiltonian of the real system. + energy_active : bool, optional + Whether the energy is an active parameter. It determines which FANPT + method is used. If set to true, FANPTContainerEParam is used. + Defaults to True. + resum : bool, optional + Indicates if we will solve the FANPT equations by re-summing the series. + Defaults to False. + ref_sd : int, optional + Index of the Slater determinant used to impose intermediate normalization. + = 1. Defaults to 0. + final_order : int, optional + Final order of the FANPT calculation. Defaults to 1. + lambda_i : float, optional + Initial lambda value for the solution of the FANPT equations. Defaults to 0.0. + lambda_f : float, optional + Lambda value up to which the FANPT calculation will be performed. Defaults to 1.0. + steps (int, optional): int, optional + Solve FANPT in n stepts between lambda_i and lambda_f. Defaults to 1. + kwargs (dict, optional): + Additional keyword arguments for FanPTContainer class. Defaults to {}. + + Raises: + TypeError: [description] + + Returns: + params: np.ndarray + Solution of the FANPT calculation. + """ + if not isinstance(fanci_wfn, FanCI): + raise TypeError("fanci_wfn must be a FanCI wavefunction") + if kwargs is None: + kwargs = {} + if solver_kwargs is None: + solver_kwargs = {} + + # Check for normalization constraint in FANCI wfn + # Assumes intermediate normalization relative to ref_sd only + if f"<\\psi_{{{ref_sd}}}|\\Psi> - v_{{{ref_sd}}}" in fanci_wfn.constraints: + inorm = True + norm_det = [(ref_sd, 1.0)] + else: + inorm = False + norm_det = None + + # Select FANPT method + if energy_active: + fanptcontainer = FANPTContainerEParam + if not fanci_wfn.mask[-1]: + fanci_wfn.unfreeze_parameter(-1) + else: + fanptcontainer = FANPTContainerEFree + if fanci_wfn.mask[-1]: + fanci_wfn.freeze_parameter(-1) + + if resum: + if energy_active: + raise ValueError( + "The energy parameter must be inactive with the resumation option." + ) + nequation = fanci_wfn.nequation + nactive = fanci_wfn.nactive + steps = 1 + if not inorm and (nequation == nactive): + norm_det = [(ref_sd, 1.0)] + elif inorm and (nequation - 1) == nactive: + fanci_wfn.remove_constraint(f"<\\psi_{{{ref_sd}}}|\\Psi> - v_{{{ref_sd}}}") + inorm = False + else: + raise ValueError( + "The necesary condition of a determined system of equations is not met." + ) + + # Get initial guess for parameters at initial lambda value. + results = fanci_wfn.optimize(guess_params, **solver_kwargs) + guess_params[fanci_wfn.mask] = results.x + params = guess_params + + # Solve FANPT equations + for l in np.linspace(lambda_i, lambda_f, steps, endpoint=False): + fanpt_container = fanptcontainer( + fanci_wfn=fanci_wfn, + params=params, + ham0=ham0, + ham1=ham1, + l=l, + inorm=inorm, + ref_sd=ref_sd, + **kwargs, + ) + + final_l = l + (lambda_f - lambda_i) / steps + print(f'Solving FanPT problem at lambda={final_l}') + + fanpt_updater = FANPTUpdater( + fanpt_container=fanpt_container, + final_order=final_order, + final_l=final_l, + solver=None, + resum=resum, + ) + new_wfn_params = fanpt_updater.new_wfn_params + new_energy = fanpt_updater.new_energy + + # These params serve as initial guess to solve the fanci equations for the given lambda. + fanpt_params = np.append(new_wfn_params, new_energy) + print('Frobenius Norm of parameters: {}'.format(np.linalg.norm(fanpt_params - params))) + print('Energy change: {}'.format(np.linalg.norm(fanpt_params[-1] - params[-1]))) + + # Initialize perturbed Hamiltonian with the current value of lambda using the static method of fanpt_container. + ham = fanpt_container.linear_comb_ham(ham1, ham0, final_l, 1 - final_l) + + # Initialize fanci wfn with the perturbed Hamiltonian. + fanci_wfn = update_fanci_wfn(ham, fanci_wfn, norm_det, fill) + + # Solve the fanci problem with fanpt_params as initial guess. + # Take the params given by fanci and use them as initial params in the + # fanpt calculation for the next lambda. + results = fanci_wfn.optimize(fanpt_params, **solver_kwargs) + + fanpt_params[fanci_wfn.mask] = results.x + params = fanpt_params + + if not energy_active: + fanci_wfn.freeze_parameter([-1]) + + # Output for debugging purposes + results["energy"] = fanpt_params[-1] + results["residuals"] = results.fun + #output = { + # "x": params, + # "variables": { + # "steps": steps, + # "inorm": inorm, + # "norm_det": norm_det, + # "eactive": fanci_wfn.mask[-1], + # "final_l": final_l, + # }, + #} + return results + + +def update_fanci_wfn(ham, fanciwfn, norm_det, fill): + fanci_class = fanciwfn.__class__ + + if isinstance(fanciwfn.wfn, pyci.fullci_wfn): + nocc = (fanciwfn.wfn.nocc_up, fanciwfn.wfn.nocc_dn) + else: + nocc = fanciwfn.wfn.nocc_up + + # Activate energy parameter + fanciwfn.unfreeze_parameter([-1]) + + # FIXME: for FanCI class + #return fanci_class( + # ham, fanciwfn.wfn, fanciwfn.nproj, fanciwfn.nparam, norm_det=norm_det, mask=fanciwfn.mask, fill=fill + #) + + # for fanpy class + # ASSUMING GeneratedFanCI class + return fanci_class( + ham, fanciwfn._fanpy_wfn, fanciwfn._wfn.nocc_up + fanciwfn._wfn.nocc_dn, + nproj=fanciwfn.nproj, wfn=fanciwfn.wfn, fill=fill, seniority=fanciwfn.seniority, + step_print=fanciwfn.step_print, step_save=fanciwfn.step_save, tmpfile=fanciwfn.tmpfile, + mask=fanciwfn._mask, objective_type=fanciwfn.objective_type, norm_det=norm_det, + param_selection=fanciwfn.indices_component_params, + constraints=fanciwfn._constraints + ) + + +def reduce_to_fock(two_int, lambda_val=0): + """Reduce given two electron integrals to that of the correspoding Fock operator. + + Parameters + ---------- + two_int : np.ndarray(K, K, K, K) + Two electron integrals of restricted orbitals. + + """ + fock_two_int = two_int * lambda_val + nspatial = two_int.shape[0] + indices = np.arange(nspatial) + fock_two_int[indices[:, None, None], indices[None, :, None], indices[None, None, :], indices[None, :, None]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, None, :], indices[None, :, None]] + fock_two_int[indices[:, None, None], indices[None, :, None], indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, :, None], indices[None, None, :]] + #fock_two_int[indices[None, :, None], indices[:, None, None], indices[None, None, :], indices[None, :, None]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, None, :], indices[None, :, None]] + #fock_two_int[indices[None, :, None], indices[:, None, None], indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, :, None], indices[None, None, :]] + + #occ_indices = np.arange(nelec // 2) + #fock_two_int[indices[:, None, None], occ_indices[None, :, None], indices[None, None, :], occ_indices[None, :, None]] = two_int[indices[:, None, None], occ_indices[None, :, None], indices[None, None, :], occ_indices[None, :, None]] + #fock_two_int[indices[:, None, None], occ_indices[None, :, None], occ_indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], occ_indices[None, :, None], occ_indices[None, :, None], indices[None, None, :]] + #fock_two_int[occ_indices[None, :, None], indices[:, None, None], indices[None, None, :], occ_indices[None, :, None]] = two_int[indices[:, None, None], occ_indices[None, :, None], indices[None, None, :], occ_indices[None, :, None]] + #fock_two_int[occ_indices[None, :, None], indices[:, None, None], occ_indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], occ_indices[None, :, None], occ_indices[None, :, None], indices[None, None, :]] + + return fock_two_int diff --git a/pyci/fanpt/__init__.py b/pyci/fanpt/__init__.py new file mode 100644 index 0000000..1175807 --- /dev/null +++ b/pyci/fanpt/__init__.py @@ -0,0 +1,17 @@ +r"""FanCI module.""" + + +__all__ = [ + "FANPTContainer", + "FANPTContainerEFree", + "FANPTContainerEParam", + "FANPTConstantTerms", + "FANPTUpdater", +] + + +from .base_fanpt_container import FANPTContainer +from .fanpt_cont_e_free import FANPTContainerEFree +from .fanpt_cont_e_param import FANPTContainerEParam +from .fanpt_updater import FANPTUpdater +from .fanpt_constant_terms import FANPTConstantTerms diff --git a/pyci/fanpt/base_fanpt_container.py b/pyci/fanpt/base_fanpt_container.py new file mode 100644 index 0000000..78bbc4c --- /dev/null +++ b/pyci/fanpt/base_fanpt_container.py @@ -0,0 +1,298 @@ +r"""Base class that contains the elements required to perform a FANPT calculation.""" + +from abc import ABCMeta, abstractmethod +import pyci + +class FANPTContainer(metaclass=ABCMeta): + r"""Container for the matrices and vectors required ot perform a FANPT calculation. + + We assume that the equations to be solved have the following structure: + + A block of nproj equations obtained by projecting the Schrödinger equation in a space + of nproj Slater determinants: + + G_1 = <1|ham(l)|psi(l)> - E * <1|psi(l)> = 0 + G_2 = <2|ham(l)|psi(l)> - E * <2|psi(l)> = 0 + .... + G_nproj = - E * = 0 + + We also have constraint equations (at least one, used to impose intermediate normalization). + It is assumed that the constraint equations only depend on the wavefunction parameters, p_k, + and are independent of the energy, E, and lambda, l. This implies that while all the vector + have nequation elements and all the matrices have nequation rows, except for the coefficient + matrix (dG/dp_k), only the first nproj elements of the vectors and the first nproj rows of + the matrices are non-zero. + + Attributes + ---------- + fanci_wfn : FanCI instance + FanCI wavefunction. + params : np.ndarray + Wavefunction parameters and energy at for the given lambda value. + ham0 : pyci.hamiltonian + PyCI Hamiltonian of the ideal system. + ham1 : pyci.hamiltonian + PyCI Hamiltonian of the real system. + l : float + Lambda value. + ref_sd : int + Index of the Slater determinant used to impose intermediate normalization. + = 1. + inorm : bool + Indicates whether we will work with intermediate normalization or not. + ham : pyci.hamiltonian + PyCI Hamiltonian for the given value of lambda. + ham = l * ham1 + (1 - l) * ham0 + f_pot : pyci.hamiltonian + PyCI Hamiltonian corresponding to the fluctuation potential. + f_pot = ham1 - ham0 + wfn_params : np.ndarray + Wavefunction parameters. + energy : float + Energy for the current value of lambda. + active_energy : bool + Indicates if the energy will be varied in the calculations. + It is False either when the energy is frozen in a E-param calculation + or in any E-free calculation. + ham_ci_op : pyci.sparse_op + PyCI sparse operator corresponding to the perturbed Hamiltonian. + f_pot_ci_op : pyci.sparse_op + PyCI sparse operator corresponding to the fluctuation potential. + ovlp_s : np.ndarray + Overlaps of the wavefunction with the determinants in the "S" space. + d_ovlp_s : np.ndarray + Derivatives of the overlaps of the wavefunction with the determinants in the "S" space + with respect to the active wavefunction parameters. + d_g_lambda : np.ndarray + Derivative of the FANPT equations with respect to the lambda parameter. + numpy array with shape (self.nequations,). + d2_g_lambda_wfnparams : np.ndarray + Derivative of the FANPT equations with respect to lambda and the wavefunction + parameters. + numpy array with shape (self.nequations, len(self.wfn_params_active)). + c_matrix : np.ndarray + Coefficient matrix of the FANPT system of equations. + numpy array with shape (self.nequations, len(self.nactive)). + + Properties + ---------- + nactive : int + Number of active parameters. + nequation : int + Number of equations. + nproj : int + Number of determinants in the projection ("P") space. + + Methods + ------- + __init__(self, fanci_wfn, params, ham0, ham1, l=0, ref_sd=0) + Initialize the FANPT container. + linear_comb_ham(ham1, ham0, a1, a0) + Return a linear combination of two PyCI Hamiltonians. + der_g_lambda(self) + Derivative of the FANPT equations with respect to the lambda parameter. + der2_g_lambda_wfnparams(self) + Derivative of the FANPT equations with respect to lambda and the wavefunction parameters. + gen_coeff_matrix(self) + Generate the coefficient matrix of the linear FANPT system of equations. + """ + + def __init__( + self, + fanci_wfn, + params, + ham0, + ham1, + l=0, + ref_sd=0, + inorm=False, + ham_ci_op=None, + f_pot_ci_op=None, + ovlp_s=None, + d_ovlp_s=None, + ): + r"""Initialize the FANPT container. + + Parameters + ---------- + fanci_wfn : FanCI instance + FanCI wavefunction. + params : np.ndarray + Wavefunction parameters and energy at for the given lambda value. + ham0 : pyci.hamiltonian + PyCI Hamiltonian of the ideal system. + ham1 : pyci.hamiltonian + PyCI Hamiltonian of the real system. + l : float + Lambda value. + ref_sd : int + Index of the Slater determinant used to impose intermediate normalization. + = 1. + ham_ci_op : {pyci.sparse_op, None} + PyCI sparse operator of the perturbed Hamiltonian. + f_pot_ci_op : {pyci.sparse_op, None} + PyCI sparse operator of the fluctuation potential. + ovlp_s : {np.ndarray, None} + Overlaps in the "S" projection space. + d_ovlp_s : {np.ndarray, None} + Derivatives of the overlaps in the "S" projection space. + """ + # Separate parameters for better readability. + self.params = params + self.wfn_params = params[:-1] + self.energy = params[-1] + self.active_energy = fanci_wfn.mask[-1] + self.inorm = inorm + + # Assign ideal and real Hamiltonians. + self.ham1 = ham1 + self.ham0 = ham0 + + # Lambda parameter. + self.l = l + + # Assign FanCI wfn + self.fanci_wfn = fanci_wfn + + # Build the perturbed Hamiltonian. + self.ham = FANPTContainer.linear_comb_ham(self.ham1, self.ham0, self.l, 1 - self.l) + + # Construct the perturbed Hamiltonian and fluctuation potential sparse operators. + if ham_ci_op: + self.ham_ci_op = ham_ci_op + else: + self.ham_ci_op = pyci.sparse_op(self.ham, self.fanci_wfn.wfn, self.fanci_wfn.nproj) + if f_pot_ci_op: + self.f_pot_ci_op = f_pot_ci_op + else: + self.f_pot = FANPTContainer.linear_comb_ham(self.ham1, self.ham0, 1.0, -1.0) + self.f_pot_ci_op = pyci.sparse_op(self.f_pot, self.fanci_wfn.wfn, self.fanci_wfn.nproj) + + if ovlp_s: + self.ovlp_s = ovlp_s + else: + self.ovlp_s = self.fanci_wfn.compute_overlap(self.wfn_params, "S") + if d_ovlp_s: + self.d_ovlp_s = d_ovlp_s + else: + self.d_ovlp_s = self.fanci_wfn.compute_overlap_deriv(self.wfn_params, "S") + + # Update Hamilonian in the fanci_wfn. + self.fanci_wfn._ham = self.ham + + # Assign ref_sd. + if self.inorm: + if f"<\\psi_{{{ref_sd}}}|\\Psi> - v_{{{ref_sd}}}" in self.fanci_wfn.constraints: + self.ref_sd = ref_sd + else: + raise KeyError( + "The normalization of the Slater determinant is not constrained" + "in the FanCI wavefunction." + ) + else: + self.ref_sd = ref_sd + + # Generate the required vectors and matrices. + self.der_g_lambda() + self.der2_g_lambda_wfnparams() + self.gen_coeff_matrix() + + @staticmethod + def linear_comb_ham(ham1, ham0, a1, a0): + r"""Return a linear combination of two PyCI Hamiltonians. + + Parameters + ---------- + ham1 : pyci.hamiltonian + PyCI Hamiltonian of the real system. + ham0 : pyci.hamiltonian + PyCI Hamiltonian of the ideal system. + a1 : float + Coefficient of the real Hamiltonian. + a0 : float + Coefficient of the ideal Hamiltonian. + + Returns + ------- + pyci.hamiltonian + """ + ecore = a1 * ham1.ecore + a0 * ham0.ecore + one_mo = a1 * ham1.one_mo + a0 * ham0.one_mo + two_mo = a1 * ham1.two_mo + a0 * ham0.two_mo + return pyci.hamiltonian(ecore, one_mo, two_mo) + + @property + def nactive(self): + r"""Return the number of active parameters. + + Returns + ------- + self.fanci_wfn.nactive : int + Number of active parameters. + """ + return self.fanci_wfn.nactive + + @property + def nequation(self): + r"""Return the number of equations. + + Returns + ------- + self.fanci_wfn.nequation : int + Number of equations (includes the number of constraints). + """ + return self.fanci_wfn.nequation + + @property + def nproj(self): + r"""Return the number of determinants in the projection "P" space. + + Returns + ------- + self.fanci_wfn.nproj + Number of determinants in the "P" space. + """ + return self.fanci_wfn.nproj + + @abstractmethod + def der_g_lambda(self): + r"""Derivative of the FANPT equations with respect to the lambda parameter. + + dG/dl + + Generates + --------- + d_g_lambda : np.ndarray + Derivative of the FANPT equations with respect to the lambda parameter. + numpy array with shape (self.nequations,). + """ + self.d_g_lambda = None + + @abstractmethod + def der2_g_lambda_wfnparams(self): + r"""Derivative of the FANPT equations with respect to lambda and the wavefunction parameters. + + d^2G/dldp_k + + Generates + --------- + d2_g_lambda_wfnparams : np.ndarray + Derivative of the FANPT equations with respect to lambda and the wavefunction + parameters. + numpy array with shape (self.nequations, len(self.wfn_params_active)). + """ + self.d2_g_lambda_wfnparams = None + + @abstractmethod + def gen_coeff_matrix(self): + r"""Generate the coefficient matrix of the linear FANPT system of equations. + + dG/dp_k + + Generates + --------- + c_matrix : np.ndarray + Coefficient matrix of the FANPT system of equations. + numpy array with shape (self.nequations, len(self.nactive)). + """ + self.c_matrix = None diff --git a/pyci/fanpt/fanpt_constant_terms.py b/pyci/fanpt/fanpt_constant_terms.py new file mode 100644 index 0000000..35dade4 --- /dev/null +++ b/pyci/fanpt/fanpt_constant_terms.py @@ -0,0 +1,167 @@ +r"""Class that generates and contains the constant terms of the FANPT system of equations.""" + +from math import factorial +import numpy as np +from .base_fanpt_container import FANPTContainer + +class FANPTConstantTerms: + r"""Generates and contains the constant terms of the FANPT system of equations. + + If the order is 1: + -dG_n/dl + + If the energy is not an active parameter: + -N * sum_k {d^2G_n/dldp_k * d^(N-1)p_k/dl^(N-1)} + + If the energy is an active parameter: + linear_term = t1 + t2 + t1 = -N * sum_k {d^2G_n/dldp_k * d^(N-1)p_k/dl^(N-1)} + t2 = - sum_k{d^2G_n/dEdp_k * sum_m{C(N,m) * d^mE/dl^m * d^(N-m)p_k/dl^(N-m)}} + + Notation + N: order + sum_k: sum over all active wfn parameters. + sum_m: sum from 1 to N-1. + C(N,m): N-choose-m (binomial coefficient). Calculated using math.factorial because + math.comb is only available starting at Python 3.8. + + Attributes + ---------- + fanpt_container : FANPTContainer + Object containing the FANPT matrices and vectors. + order : int + Order of the current FANPT iteration. + previous_responses : np.ndarray + Previous responses of the FANPT calculations up to order = order -1 + + Methods + ------- + __init__(self, fanpt_container, order=1, previous_responses=None) + Initialize the constant terms. + assign_fanpt_container(self, fanpt_container) + Assign the FANPT container. + assign_order(self, order) + Assign the order. + assign_previous_responses(self, previous_responses) + Assign the previous responses. + gen_constant_terms(self) + Generate the constant terms. + """ + + def __init__(self, fanpt_container, order=1, previous_responses=[]): + r"""Initialize the constant terms. + + Parameters + ---------- + fanpt_container : FANPTContainer + Object containing the FANPT matrices and vectors. + order : int + Order of the current FANPT iteration. + previous_responses : np.ndarray + Previous responses of the FANPT calculations up to order = order -1 + """ + self.assign_fanpt_container(fanpt_container=fanpt_container) + self.assign_order(order=order) + self.assign_previous_responses(previous_responses=previous_responses) + self.gen_constant_terms() + + def assign_fanpt_container(self, fanpt_container): + r"""Assign the FANPT container. + + Parameters + ---------- + fanpt_container : FANPTContainer + FANPTContainer object that contains all the matrices and vectors required to perform + the FANPT calculation. + + Raises + ------ + TypeError + If fanpt_container is not a child of FANPTContainer. + """ + if not isinstance(fanpt_container, FANPTContainer): + raise TypeError("fanpt_container must be a child of FANPTContainer") + self.fanpt_container = fanpt_container + + def assign_order(self, order): + r"""Assign the order. + + Parameters + ---------- + order : int + Order of the current FANPT iteration. + + Raises + ------ + TypeError + If order is not an int. + ValueError + If order is negative. + """ + if not isinstance(order, int): + raise TypeError("order must be an integer.") + if order < 0: + raise ValueError("order must be non-negative") + self.order = order + + def assign_previous_responses(self, previous_responses): + r"""Assign the previous responses. + + Parameters + ---------- + previous_responses : np.ndarray + Previous responses of the FANPT calculations up to order = order -1. + + Raises + ------ + TypeError + If previous_responses is not a numpy array. + If the elements of previous_responses are not numpy arrays. + ValueError + If previous_responses is None and order is not 1. + If the shape of previous_responses is not equal to (order - 1, nactive). + """ + if self.order == 1: + self.previous_responses = previous_responses + else: + if not isinstance(previous_responses, np.ndarray): + raise TypeError("previous_responses must be a numpy array.") + if not all([isinstance(response, np.ndarray) for response in previous_responses]): + raise TypeError("The elements of previous_responses must be numpy arrays.") + if previous_responses.shape != (self.order - 1, self.fanpt_container.nactive): + raise ValueError( + "The shape of previous_responses must be ({}, {}).".format( + self.order - 1, self.fanpt_container.nactive + ) + ) + self.previous_responses = previous_responses + + def gen_constant_terms(self): + r"""Generate the constant terms. + + Returns + ------- + constant_terms : np.ndarray + Constant terms of the FANPT linear system of equations. + numpy array with shape (nequation,). + """ + if self.order == 1: + constant_terms = -self.fanpt_container.d_g_lambda + else: + if self.fanpt_container.active_energy: + r_vector = np.zeros(self.fanpt_container.nactive - 1) + for o in range(1, self.order): + comb = factorial(self.order) / (factorial(o) * factorial(self.order - o)) + r_vector += ( + comb + * self.previous_responses[o - 1][-1] + * self.previous_responses[self.order - o - 1][:-1] + ) + constant_terms = -self.order * np.dot( + self.fanpt_container.d2_g_lambda_wfnparams, self.previous_responses[-1][:-1] + ) - np.dot(self.fanpt_container.d2_g_e_wfnparams, r_vector) + else: + constant_terms = -self.order * np.dot( + self.fanpt_container.d2_g_lambda_wfnparams, self.previous_responses[-1] + ) + self.constant_terms = constant_terms diff --git a/pyci/fanpt/fanpt_cont_e_free.py b/pyci/fanpt/fanpt_cont_e_free.py new file mode 100644 index 0000000..3585a72 --- /dev/null +++ b/pyci/fanpt/fanpt_cont_e_free.py @@ -0,0 +1,259 @@ +r"""Class that contains the elements required to perform a FANPT calculation with implicit E.""" + +import numpy as np +from .fanpt_cont_e_param import FANPTContainerEParam + +class FANPTContainerEFree(FANPTContainerEParam): + r""" + Container for the matrices and vectors required ot perform a FANPT calculation. + + We assume that the equations to be solved have the following structure: + + A block of nproj equations obtained by projecting the Schrödinger equation in a space + of nproj Slater determinants: + + G_1 = <1|ham(l)|psi(l)> - E * <1|psi(l)> = 0 + G_2 = <2|ham(l)|psi(l)> - E * <2|psi(l)> = 0 + .... + G_nproj = - E * = 0 + + We also have constraint equations (at least one, used to impose intermediate normalization). + It is assumed that the constraint equations only depend on the wavefunction parameters, p_k, + and are independent of the energy, E, and lambda, l. This implies that while all the vector + have nequation elements and all the matrices have nequation rows, except for the coefficient + matrix (dG/dp_k), only the first nproj elements of the vectors and the first nproj rows of + the matrices are non-zero. + + Attributes + ---------- + fanci_wfn : FanCI instance + FanCI wavefunction. + params : np.ndarray + Wavefunction parameters and energy at for the given lambda value. + ham0 : pyci.hamiltonian + PyCI Hamiltonian of the ideal system. + ham1 : pyci.hamiltonian + PyCI Hamiltonian of the real system. + l : float + Lambda value. + ref_sd : int + Index of the Slater determinant used to impose intermediate normalization. + = 1. + inorm : bool + Indicates whether we will work with intermediate normalization or not. + ham : pyci.hamiltonian + PyCI Hamiltonian for the given value of lambda. + ham = l * ham1 + (1 - l) * ham0 + f_pot : pyci.hamiltonian + PyCI Hamiltonian corresponding to the fluctuation potential. + f_pot = ham1 - ham0 + wfn_params : np.ndarray + Wavefunction parameters. + energy : float + Energy for the current value of lambda. + active_energy : bool + Indicates if the energy will be varied in the calculations. + It is False either when the energy is frozen in a E-param calculation + or in any E-free calculation. + ham_ci_op : pyci.sparse_op + PyCI sparse operator corresponding to the perturbed Hamiltonian. + f_pot_ci_op : pyci.sparse_op + PyCI sparse operator corresponding to the fluctuation potential. + ovlp_s : np.ndarray + Overlaps of the wavefunction with the determinants in the "S" space. + d_ovlp_s : np.ndarray + Derivatives of the overlaps of the wavefunction with the determinants in the "S" space + with respect to the active wavefunction parameters. + d_g_lambda : np.ndarray + Derivative of the FANPT equations with respect to the lambda parameter. + numpy array with shape (self.nequations,). + super_d_g_lambda : np.array + Derivarive of the FANPT equations with respect to the lambda parameter + as if it was calculated in the E-param way. + numpy array with shape (self.nequations,). + d2_g_lambda_wfnparams : np.ndarray + Derivative of the FANPT equations with respect to lambda and the wavefunction + parameters. + numpy array with shape (self.nequations, len(self.wfn_params_active)). + c_matrix : np.ndarray + Coefficient matrix of the FANPT system of equations. + numpy array with shape (self.nequations, len(self.nactive)). + + Properties + ---------- + nactive : int + Number of active parameters. + nequation : int + Number of equations. + nproj : int + Number of determinants in the projection ("P") space. + + Methods + ------- + __init__(self, fanci_wfn, params, ham0, ham1, l=0, ref_sd=0) + Initialize the FANPT container. + linear_comb_ham(ham1, ham0, a1, a0) + Return a linear combination of two PyCI Hamiltonians. + der_g_lambda(self) + Derivative of the FANPT equations with respect to the lambda parameter. + der2_g_lambda_wfnparams(self) + Derivative of the FANPT equations with respect to lambda and the wavefunction parameters. + gen_coeff_matrix(self) + Generate the coefficient matrix of the linear FANPT system of equations. + """ + + def __init__( + self, + fanci_wfn, + params, + ham0, + ham1, + l=0, + ref_sd=0, + inorm=False, + ham_ci_op=None, + f_pot_ci_op=None, + ovlp_s=None, + d_ovlp_s=None, + ): + r"""Initialize the FANPT container. + + Parameters + ---------- + fanci_wfn : FanCI instance + FanCI wavefunction. + params : np.ndarray + Wavefunction parameters and energy at for the given lambda value. + ham0 : pyci.hamiltonian + PyCI Hamiltonian of the ideal system. + ham1 : pyci.hamiltonian + PyCI Hamiltonian of the real system. + l : float + Lambda value. + ref_sd : int + Index of the Slater determinant used to impose intermediate normalization. + = 1. + ham_ci_op : {pyci.sparse_op, None} + PyCI sparse operator of the perturbed Hamiltonian. + f_pot_ci_op : {pyci.sparse_op, None} + PyCI sparse operator of the fluctuation potential. + ovlp_s : {np.ndarray, None} + Overlaps in the "S" projection space. + d_ovlp_s : {np.ndarray, None} + Derivatives of the overlaps in the "S" projection space. + """ + if fanci_wfn.mask[-1]: + raise TypeError("The energy cannot be an active parameter.") + else: + super().__init__( + fanci_wfn, + params, + ham0, + ham1, + l, + ref_sd, + inorm, + ham_ci_op, + f_pot_ci_op, + ovlp_s, + d_ovlp_s, + ) + + def der_g_lambda(self): + r"""Derivative of the FANPT equations with respect to the lambda parameter. + + dG/dl = - * + + dG/dl = super() - * + + Generates + --------- + d_g_lambda : np.ndarray + Derivative of the FANPT equations with respect to the lambda parameter. + numpy array with shape (self.nequations,). + """ + super().der_g_lambda() + self.super_d_g_lambda = self.d_g_lambda.copy() + if self.inorm: + self.d_g_lambda[: self.nproj] -= ( + self.d_g_lambda[self.ref_sd] * self.ovlp_s[: self.nproj] + ) + else: + self.d_g_lambda[: self.nproj] -= ( + self.d_g_lambda[self.ref_sd] * self.ovlp_s[: self.nproj] / self.ovlp_s[self.ref_sd] + ) + + def der2_g_lambda_wfnparams(self): + r"""Derivative of the FANPT equations with respect to lambda and the wavefunction parameters. + + d^2G/dldp_k = - * + - * + + d^2G/dldp_k = super() - * + - * + + Generates + --------- + d2_g_lambda_wfnparams : np.ndarray + Derivative of the FANPT equations with respect to lambda and the wavefunction + parameters. + numpy array with shape (self.nequations, len(self.wfn_params_active)). + """ + super().der2_g_lambda_wfnparams() + if self.inorm: + self.d2_g_lambda_wfnparams[: self.nproj] -= self.d2_g_lambda_wfnparams[ + self.ref_sd + ] * self.ovlp_s[: self.nproj].reshape(self.nproj, 1) + self.d2_g_lambda_wfnparams[: self.nproj] -= ( + self.super_d_g_lambda[self.ref_sd] * self.d_ovlp_s[: self.nproj] + ) + else: + self.d2_g_lambda_wfnparams[: self.nproj] -= ( + ( + self.d2_g_lambda_wfnparams[self.ref_sd] + - self.super_d_g_lambda[self.ref_sd] + * self.d_ovlp_s[self.ref_sd] + / self.ovlp_s[self.ref_sd] + ) + * self.ovlp_s[: self.nproj].reshape(self.nproj, 1) + / self.ovlp_s[self.ref_sd] + ) + self.d2_g_lambda_wfnparams[: self.nproj] -= ( + self.super_d_g_lambda[self.ref_sd] + * self.d_ovlp_s[: self.nproj] + / self.ovlp_s[self.ref_sd] + ) + + def gen_coeff_matrix(self): + r"""Generate the coefficient matrix of the linear FANPT system of equations. + + dG/dp_k = - E * + - * + + dG/dp_k = super() - * + + Notes + ----- + - There is not a dG/dE column. + - The E that appears in the above equations is just . + + Generates + --------- + c_matrix : np.ndarray + Coefficient matrix of the FANPT system of equations. + numpy array with shape (self.nequations, len(self.nactive)). + """ + super().gen_coeff_matrix() + f_proj = np.empty((self.nproj, self.nactive), order="F") + for f_proj_col, d_ovlp_col in zip(f_proj.transpose(), self.d_ovlp_s.transpose()): + self.ham_ci_op(d_ovlp_col, out=f_proj_col) + if self.inorm: + self.c_matrix[: self.nproj] -= f_proj[self.ref_sd] * self.ovlp_s[: self.nproj].reshape( + self.nproj, 1 + ) + else: + self.c_matrix[: self.nproj] -= ( + (f_proj[self.ref_sd] - self.energy * self.d_ovlp_s[self.ref_sd]) + * self.ovlp_s[: self.nproj].reshape(self.nproj, 1) + / self.ovlp_s[self.ref_sd] + ) diff --git a/pyci/fanpt/fanpt_cont_e_param.py b/pyci/fanpt/fanpt_cont_e_param.py new file mode 100644 index 0000000..0ccb538 --- /dev/null +++ b/pyci/fanpt/fanpt_cont_e_param.py @@ -0,0 +1,236 @@ +r"""Class that contains the elements required to perform a FANPT calculation with explicit E.""" + +import numpy as np +from .base_fanpt_container import FANPTContainer + +class FANPTContainerEParam(FANPTContainer): + r"""Container for the matrices and vectors required ot perform a FANPT calculation. + + We assume that the equations to be solved have the following structure: + + A block of nproj equations obtained by projecting the Schrödinger equation in a space + of nproj Slater determinants: + + G_1 = <1|ham(l)|psi(l)> - E * <1|psi(l)> = 0 + G_2 = <2|ham(l)|psi(l)> - E * <2|psi(l)> = 0 + .... + G_nproj = - E * = 0 + + We also have constraint equations (at least one, used to impose intermediate normalization). + It is assumed that the constraint equations only depend on the wavefunction parameters, p_k, + and are independent of the energy, E, and lambda, l. This implies that while all the vector + have nequation elements and all the matrices have nequation rows, except for the coefficient + matrix (dG/dp_k), only the first nproj elements of the vectors and the first nproj rows of + the matrices are non-zero. + + Attributes + ---------- + fanci_wfn : FanCI instance + FanCI wavefunction. + params : np.ndarray + Wavefunction parameters and energy at for the given lambda value. + ham0 : pyci.hamiltonian + PyCI Hamiltonian of the ideal system. + ham1 : pyci.hamiltonian + PyCI Hamiltonian of the real system. + l : float + Lambda value. + ref_sd : int + Index of the Slater determinant used to impose intermediate normalization. + = 1. + inorm : bool + Indicates whether we will work with intermediate normalization or not. + ham : pyci.hamiltonian + PyCI Hamiltonian for the given value of lambda. + ham = l * ham1 + (1 - l) * ham0 + f_pot : pyci.hamiltonian + PyCI Hamiltonian corresponding to the fluctuation potential. + f_pot = ham1 - ham0 + wfn_params : np.ndarray + Wavefunction parameters. + energy : float + Energy for the current value of lambda. + active_energy : bool + Indicates if the energy will be varied in the calculations. + It is False either when the energy is frozen in a E-param calculation + or in any E-free calculation. + ham_ci_op : pyci.sparse_op + PyCI sparse operator corresponding to the perturbed Hamiltonian. + f_pot_ci_op : pyci.sparse_op + PyCI sparse operator corresponding to the fluctuation potential. + ovlp_s : np.ndarray + Overlaps of the wavefunction with the determinants in the "S" space. + d_ovlp_s : np.ndarray + Derivatives of the overlaps of the wavefunction with the determinants in the "S" space + with respect to the active wavefunction parameters. + d_g_lambda : np.ndarray + Derivative of the FANPT equations with respect to the lambda parameter. + numpy array with shape (self.nequations,). + d2_g_lambda_wfnparams : np.ndarray + Derivative of the FANPT equations with respect to lambda and the wavefunction + parameters. + numpy array with shape (self.nequations, len(self.wfn_params_active)). + d2_g_e_wfnparams : np.ndarray + Derivative of the FANPT equations with respect to the energy and the wavefunction + parameters. + numpy array with shape (self.nequations, len(self.wfn_params_active)). + c_matrix : np.ndarray + Coefficient matrix of the FANPT system of equations. + numpy array with shape (self.nequations, len(self.nactive)). + + Properties + ---------- + nactive : int + Number of active parameters. + nequation : int + Number of equations. + nproj : int + Number of determinants in the projection ("P") space. + + Methods + ------- + __init__(self, fanci_wfn, params, ham0, ham1, l=0, ref_sd=0) + Initialize the FANPT container. + linear_comb_ham(ham1, ham0, a1, a0) + Return a linear combination of two PyCI Hamiltonians. + der_g_lambda(self) + Derivative of the FANPT equations with respect to the lambda parameter. + der2_g_lambda_wfnparams(self) + Derivative of the FANPT equations with respect to lambda and the wavefunction parameters. + der2_g_e_wfnparams(self) + Derivative of the FANPT equations with respect to the energy and the wavefunction + parameters. + gen_coeff_matrix(self) + Generate the coefficient matrix of the linear FANPT system of equations. + """ + + def __init__( + self, + fanci_wfn, + params, + ham0, + ham1, + l=0, + ref_sd=0, + inorm=False, + ham_ci_op=None, + f_pot_ci_op=None, + ovlp_s=None, + d_ovlp_s=None, + ): + r"""Initialize the FANPT container. + + Parameters + ---------- + fanci_wfn : FanCI instance + FanCI wavefunction. + params : np.ndarray + Wavefunction parameters and energy at for the given lambda value. + ham0 : pyci.hamiltonian + PyCI Hamiltonian of the ideal system. + ham1 : pyci.hamiltonian + PyCI Hamiltonian of the real system. + l : float + Lambda value. + ref_sd : int + Index of the Slater determinant used to impose intermediate normalization. + = 1. + ham_ci_op : {pyci.sparse_op, None} + PyCI sparse operator of the perturbed Hamiltonian. + f_pot_ci_op : {pyci.sparse_op, None} + PyCI sparse operator of the fluctuation potential. + ovlp_s : {np.ndarray, None} + Overlaps in the "S" projection space. + d_ovlp_s : {np.ndarray, None} + Derivatives of the overlaps in the "S" projection space. + """ + super().__init__( + fanci_wfn, + params, + ham0, + ham1, + l, + ref_sd, + inorm, + ham_ci_op, + f_pot_ci_op, + ovlp_s, + d_ovlp_s, + ) + self.der2_g_e_wfnparams() + + def der_g_lambda(self): + r"""Derivative of the FANPT equations with respect to the lambda parameter. + + dG/dl = + + Generates + --------- + d_g_lambda : np.ndarray + Derivative of the FANPT equations with respect to the lambda parameter. + numpy array with shape (self.nequations,). + """ + f = np.zeros(self.nequation) + f_proj = f[: self.nproj] + self.f_pot_ci_op(self.ovlp_s, out=f_proj) + self.d_g_lambda = f + + def der2_g_lambda_wfnparams(self): + r"""Derivative of the FANPT equations with respect to lambda and the wavefunction parameters. + + d^2G/dldp_k = + + Generates + --------- + d2_g_lambda_wfnparams : np.ndarray + Derivative of the FANPT equations with respect to lambda and the wavefunction + parameters. + numpy array with shape (self.nequations, len(self.wfn_params_active)). + """ + if self.active_energy: + ncolumns = self.nactive - 1 + else: + ncolumns = self.nactive + f = np.zeros((self.nequation, ncolumns), order="F") + f_proj = f[: self.nproj] + for f_proj_col, d_ovlp_col in zip(f_proj.transpose(), self.d_ovlp_s.transpose()): + self.f_pot_ci_op(d_ovlp_col, out=f_proj_col) + self.d2_g_lambda_wfnparams = f + + def der2_g_e_wfnparams(self): + r"""Derivative of the FANPT equations with respect to the energy and the wavefunction + parameters. + + d^2G/dEdp_k = - + + Generates + --------- + d2_g_e_wfnparams : np.ndarray + Derivative of the FANPT equations with respect to the energy and the wavefunction + parameters. + numpy array with shape (self.nequations, len(self.wfn_params_active)). + """ + if self.active_energy: + ncolumns = self.nactive - 1 + f = np.zeros((self.nequation, ncolumns), order="F") + f[: self.nproj] = -self.d_ovlp_s[: self.nproj] + self.d2_g_e_wfnparams = f + else: + self.d2_g_e_wfnparams = None + + def gen_coeff_matrix(self): + r"""Generate the coefficient matrix of the linear FANPT system of equations. + + dG/dp_k = - E * + + If the energy is active, the last column of the matrix has the form: + + dG/dE = - + + Generates + --------- + c_matrix : np.ndarray + Coefficient matrix of the FANPT system of equations. + numpy array with shape (self.nequations, len(self.nactive)). + """ + self.c_matrix = self.fanci_wfn.compute_jacobian(self.params) diff --git a/pyci/fanpt/fanpt_updater.py b/pyci/fanpt/fanpt_updater.py new file mode 100644 index 0000000..4d52f71 --- /dev/null +++ b/pyci/fanpt/fanpt_updater.py @@ -0,0 +1,374 @@ +r"""FANPT updater.""" + +from functools import partial +from math import factorial +import numpy as np + +import pyci +from .base_fanpt_container import FANPTContainer +from .fanpt_constant_terms import FANPTConstantTerms + +class FANPTUpdater: + r"""Solve the FANPT equations up to a given order, updates the wavefunction parameters and + energy. + + Attributes + ---------- + fanpt_container : FANPTContainer + Object containing the FANPT matrices and vectors. + final_order : int + Final order of the FANPT calculation. + final_l : float + Lambda value up to which the FANPT calculation will be performed. + solver : callable + Solver that will be used to solve the FANPT equations. + resum : bool + Indicates if we will solve the FANPT equations by re-summing the series. + inverse : np.ndarray + Inverse of the coefficient matrix. + A_matrix : np.ndarray + A matrix in the resummation equations. + b_vector : np.ndarray + b vector in the resummation equations. + responses : np.ndarray + Responses up to the specified order. + new_wfn_params : np.ndarray + New wavefunction parameters. + new_energy : float + New energy. + new_ham_op : pyci.sparse_op + PyCI sparse operator of the perturbed Hamiltonian at the new value of lambda. + new_ovlp_s : np.ndarray + Overlaps of the wavefunction in the "S" projection space with the new parameters. + fanpt_e : float + Energy calculated as the sum of the FANPT responses of the variable E. + + Methods + ------- + __init__(self, fanpt_container, final_order=1, final_l=1, solver=None) + Initialize the updater. + assign_fanpt_container(self, fanpt_container) + Assign the FANPT container. + assign_final_l(self, final_l) + Assign the final lambda. + inverse_coeff_matrix(self) + Return the inverse of the matrix of derivatives of the FANPT equations with respect to + the active parameters. + resum_matrix(self) + Return the matrix that appears in the resummation expressions. + resum_vector(self) + Return the vector that appears in the resummation expressions. + fanpt_resum_correction(self): + Return the resummation of all the responses. + assign_final_order(self, final_order) + Assign the final order. + assign_solver(self, solver) + Assign the solver. + get_responses(self) + Find the responses up to the final order. + params_updater(self) + Update the wavefunction parameters with the new responses up to the given value of + final_lambda. + energy_ham_ovlp_updater(self) + Update the energy, Hamiltonian sparse operator, and wavefunction overlaps. + fanpt_e_response(self) + Return the energy calculated as the sum of the FANPT responses of the E variable. + """ + + def __init__(self, fanpt_container, final_order=1, final_l=1.0, solver=None, resum=False): + r"""Initialize the updater. + + Parameters + ---------- + fanpt_container : FANPTContainer + Object containing the FANPT matrices and vectors. + final_order : int + Final order of the FANPT calculation. + final_l : float + Lambda value up to which the FANPT calculation will be performed. + solver : callable + Solver that will be used to solve the FANPT equations. + resum : bool + Indicates if we will solve the FANPT equations by re-summing the series. + """ + self.assign_fanpt_container(fanpt_container=fanpt_container) + self.assign_final_l(final_l=final_l) + self.assign_resum(resum) + if self.resum: + self.inverse_coeff_matrix() + self.resum_matrix() + self.resum_vector() + self.fanpt_resum_correction() + else: + self.assign_final_order(final_order=final_order) + self.assign_solver(solver=solver) + self.get_responses() + self.params_updater() + self.energy_ham_ovlp_updater() + if self.fanpt_container.active_energy: + self.fanpt_e_response() + + def assign_fanpt_container(self, fanpt_container): + r"""Assign the FANPT container. + + Parameters + ---------- + fanpt_container : FANPTContainer + FANPTContainer object that contains all the matrices and vectors required to perform + the FANPT calculation. + + Raises + ------ + TypeError + If fanpt_container is not a child of FANPTContainer. + """ + if not isinstance(fanpt_container, FANPTContainer): + raise TypeError("fanpt_container must be a child of FANPTContainer") + self.fanpt_container = fanpt_container + + def assign_resum(self, resum): + r"""Assign the value of resum. + + Returns + ------ + resum : bool + Indicates if we will solve the FANPT equations by re-summing the series. + + Note + ---- + The resummation can only be performed when the energy is inactive, and for a determined + system of equations, that is, if we have the same number of variables as FANPT equations, + hence the condition: nequation == nactive. + """ + if self.fanpt_container.active_energy: + self.resum = False + elif resum: + if self.fanpt_container.nequation == self.fanpt_container.nactive: + self.resum = resum + else: + self.resum = False + else: + self.resum = False + + def inverse_coeff_matrix(self): + r"""Return the inverse of the matrix of derivatives of the FANPT equations with respect to + the active parameters. + + Returns + ------- + inverse : np.ndarray + Inverse of the matrix of derivatives of the FANPT equations with respect to the + active parameters. + """ + self.inverse = np.linalg.inv(self.fanpt_container.c_matrix) + + def resum_matrix(self): + r"""Return the matrix that appears in the resummation expressions. + + Returns + ------- + A_matrix : np.ndarray + Matrix that appears in the resummation expressions. + + A = [(dG/dp_k)^-1] * [d^2G/dldp_k] + """ + self.A_matrix = np.dot(self.inverse, self.fanpt_container.d2_g_lambda_wfnparams) + + def resum_vector(self): + r"""Return the vector that appears in the resummation expressions. + + Returns + ------- + b_vector : np.ndarray + Vector that appears in the resummation expressions. + + b = [(dG/dp_k)^-1] * dG/dl + """ + self.b_vector = np.dot(self.inverse, self.fanpt_container.d_g_lambda) + + def fanpt_resum_correction(self): + r"""Return the resummation of all the responses. + + Returns + ------- + resum_result_correction : np.ndarray + Sum of all the responses up to infinite order. + + -l * [(l * A + I)^-1] * b + I : identity matrix + """ + self.resum_correction = -( + self.final_l + * np.dot( + np.linalg.inv( + self.final_l * self.A_matrix + np.identity(self.fanpt_container.nequation) + ), + self.b_vector, + ) + ) + + def assign_final_order(self, final_order): + r"""Assign the final order. + + Parameters + ---------- + final_order : int + Order of the current FANPT iteration. + + Raises + ------ + TypeError + If final_order is not an int. + ValueError + If final_order is negative. + """ + if not isinstance(final_order, int): + raise TypeError("final_order must be an integer.") + if final_order < 0: + raise ValueError("final_order must be non-negative.") + self.final_order = final_order + + def assign_final_l(self, final_l): + r"""Assign the final lambda. + + Parameters + ---------- + final_l : float + Lambda value up to which the FANPT calculation will be performed. + + Raises + ------ + TypeError + If final_l is not a float. + ValueError + If final_l is not between self.fanpt_container.l_param and 1. + """ + if not isinstance(final_l, float): + raise TypeError("final_l must be given as a float.") + if not self.fanpt_container.l < final_l <= 1.0: + raise ValueError( + "final_l must be greater than {} and lower or equal than 1.".format( + self.fanpt_container.l + ) + ) + self.final_l = final_l + + def assign_solver(self, solver): + r"""Assign solver.""" + if solver is None: + #self.solver = partial(np.linalg.lstsq, rcond=None) + self.solver = partial(np.linalg.lstsq, rcond=1e-6) + + def get_responses(self): + r"""Find the responses up to the final order. + + Returns + ------- + resp_matrix : np.ndarray + Numpy array with shape (final_order, nactive). + resp_matrix[k] contains the responses of order k+1. + If the energy is active, the final element of each row is the response of the energy. + """ + resp_matrix = np.zeros((self.final_order, self.fanpt_container.nactive)) + for o in range(1, self.final_order + 1): + c_terms = FANPTConstantTerms( + fanpt_container=self.fanpt_container, + order=o, + previous_responses=resp_matrix[: o - 1, :], + ) + constant_terms = c_terms.constant_terms + resp_matrix[o - 1] = self.solver(self.fanpt_container.c_matrix, constant_terms)[0] + self.responses = resp_matrix + + def params_updater(self): + r"""Update the wavefunction parameters with the new responses up to the given value of + final_lambda. + + Returns + ------- + new_wfn : BaseWavefunction + Updated wavefunction with the FANPT responses. + Does not modify the value of self.fanpt_container.wfn (it uses a deepcopy). + + Notes + ----- + It only updates the active wavefunction parameters. + """ + wfn_params = self.fanpt_container.wfn_params + wfn_params = wfn_params.flatten() + if self.resum: + corrections = self.resum_correction + else: + l0 = self.fanpt_container.l + dl = np.array( + [ + (self.final_l - l0) ** order / factorial(order) + for order in range(1, self.final_order + 1) + ] + ) + if self.fanpt_container.active_energy: + wfn_responses = self.responses[:, :-1] + else: + wfn_responses = self.responses + corrections = np.sum(wfn_responses * dl.reshape(self.final_order, 1), axis=0) + active_wfn_indices = np.where(self.fanpt_container.fanci_wfn.mask[:-1])[0] + for c, active_index in zip(corrections, active_wfn_indices): + wfn_params[active_index] += c + self.new_wfn_params = wfn_params + + def energy_ham_ovlp_updater(self): + r""""Update the energy, Hamiltonian sparse operator, and wavefunctoin overlaps. + + Generates + --------- + new_energy : float + E = + new_ham_op : pyci.sparse_op + PyCI sparse operator of the perturbed Hamiltonian at the new value of lambda. + new_ovlp_s : np.ndarray + Overlaps of the wavefunction in the "S" projection space with the new parameters. + + Notes + ----- + This E satisfies the 2n + 1 rule. + """ + new_ham = FANPTContainer.linear_comb_ham( + self.fanpt_container.ham1, self.fanpt_container.ham0, self.final_l, 1 - self.final_l + ) + new_ham_op = pyci.sparse_op( + new_ham, self.fanpt_container.fanci_wfn.wfn, self.fanpt_container.nproj + ) + new_ovlp_s = self.fanpt_container.fanci_wfn.compute_overlap(self.new_wfn_params, "S") + f = np.empty(self.fanpt_container.nproj, dtype=pyci.c_double) + new_ham_op(new_ovlp_s, out=f) + if self.fanpt_container.inorm: + energy = f[self.fanpt_container.ref_sd] + else: + energy = f[self.fanpt_container.ref_sd] / new_ovlp_s[self.fanpt_container.ref_sd] + self.new_energy = energy + self.new_ham_op = new_ham_op + self.new_ovlp_s = new_ovlp_s + + def fanpt_e_response(self): + r"""Return the energy calculated as the sum of the FANPT responses of the E variable. + + Generates + --------- + fanpt_e : float + Energy calculated as the sum of the FANPT responses of the variable E. + E(final_l) = E(l) + sum_k{1/k! * (final_l - l)^k * d^kE/dl^k} + + Notes + ----- + This E does not satisfy the 2n + 1 rule. + """ + e0 = self.fanpt_container.energy + l0 = self.fanpt_container.l + dl = np.array( + [ + (self.final_l - l0) ** order / factorial(order) + for order in range(1, self.final_order + 1) + ] + ) + e_responses = self.responses[:, -1] + self.fanpt_e = e0 + np.sum(e_responses * dl) From e74a4b97ac8c335966e0d89ca5a12893e5bccc46 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Fri, 12 Apr 2024 14:12:34 -0500 Subject: [PATCH 02/11] FanPT Wrapper added --- pyci/fanci/__init__.py | 1 + pyci/fanci/fanpt_wrapper.py | 1 - setup.py | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyci/fanci/__init__.py b/pyci/fanci/__init__.py index 150f7e4..0ee4b68 100644 --- a/pyci/fanci/__init__.py +++ b/pyci/fanci/__init__.py @@ -18,3 +18,4 @@ from .ap1rog import AP1roG from .detratio import DetRatio from .pccds import pCCDS +from .fanpt_wrapper import reduce_to_fock, solve_fanpt diff --git a/pyci/fanci/fanpt_wrapper.py b/pyci/fanci/fanpt_wrapper.py index e285a76..4e1ff9f 100644 --- a/pyci/fanci/fanpt_wrapper.py +++ b/pyci/fanci/fanpt_wrapper.py @@ -6,7 +6,6 @@ from .fanci import FanCI from ..fanpt import FANPTUpdater, FANPTContainerEParam, FANPTContainerEFree - def solve_fanpt( fanci_wfn, ham0, diff --git a/setup.py b/setup.py index 1f0f57d..5d57e2d 100644 --- a/setup.py +++ b/setup.py @@ -61,7 +61,7 @@ } -packages = ["pyci", "pyci.fanci", "pyci.test", "pyci.fanci.test"] +packages = ["pyci", "pyci.fanci", "pyci.fanpt", "pyci.test", "pyci.fanci.test"] package_data = { From a8cc99931b537631faca4f12d71857bc7e670f03 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Fri, 12 Apr 2024 14:31:35 -0500 Subject: [PATCH 03/11] Added masks support in FanCI class --- pyci/fanci/fanci.py | 170 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 144 insertions(+), 26 deletions(-) diff --git a/pyci/fanci/fanci.py b/pyci/fanci/fanci.py index 4d81ac7..a296ffb 100644 --- a/pyci/fanci/fanci.py +++ b/pyci/fanci/fanci.py @@ -4,18 +4,12 @@ """ from abc import ABCMeta, abstractmethod - from collections import OrderedDict - from typing import Any, Callable, Dict, List, Sequence, Tuple, Union - import numpy as np - from scipy.optimize import OptimizeResult, least_squares, root - import pyci - __all__ = [ "FanCI", ] @@ -50,6 +44,15 @@ def nparam(self) -> int: """ return self._nparam + + @property + def nactive(self) -> int: + r""" + Number of active parameters for this FanCI problem. + + """ + return self._nactive + @property def constraints(self) -> Tuple[str]: r""" @@ -58,6 +61,14 @@ def constraints(self) -> Tuple[str]: """ return tuple(self._constraints.keys()) + @property + def mask(self) -> np.ndarray: + r""" + Frozen parameter mask. + + """ + return self._mask_view + @property def ham(self) -> pyci.hamiltonian: r""" @@ -147,6 +158,7 @@ def __init__( norm_param: Sequence[Tuple[int, float]] = None, norm_det: Sequence[Tuple[int, float]] = None, constraints: Dict[str, Tuple[Callable, Callable]] = None, + mask: Sequence[int] = None, fill: str = "excitation", ) -> None: r""" @@ -170,6 +182,10 @@ def __init__( them. constraints : Dict[str, Tuple[Callable, Callable]], optional Pairs of functions (f, dfdx) corresponding to additional constraints. + mask : Sequence[int] or Sequence[bool], optional + List of parameters to freeze. If the list contains ints, then each element corresponds + to a frozen parameter. If the list contains bools, then each element indicates whether + that parameter is active (True) or frozen (False). fill : ('excitation' | 'seniority' | None) Whether to fill the projection ("P") space by excitation level, by seniority, or not at all (in which case ``wfn`` must already be filled). @@ -193,8 +209,26 @@ def __init__( name = f"<\\psi_{{{index}}}|\\Psi> - v_{{{index}}}" constraints[name] = self.make_det_constraint(index, value) + # Generate mask (True for active parameter, False for frozen parameter) + if mask is None: + mask = np.ones(nparam, dtype=np.bool) + else: + mask = np.array(mask) + if mask.dtype == np.bool: + # Check length of boolean mask + if mask.size != nparam: + raise ValueError(f"Mask size is {mask.size}; must have size `nparam={nparam}`") + else: + # Convert integer mask to boolean + ints = mask + mask = np.ones(nparam, dtype=np.bool) + mask[ints] = 0 + # Number of nonlinear equations and active parameters nequation = nproj + len(constraints) + nactive = mask.sum() + if nequation < nactive: + print(f"WARNING: System is underdetermined with dimensions {nequation}, {nactive}. Continuing anyways") # Generate determinant spaces wfn = fill_wavefunction(wfn, nproj, fill) @@ -210,7 +244,10 @@ def __init__( self._nequation = nequation self._nproj = nproj self._nparam = nparam + self._nactive = nactive self._constraints = constraints + self._mask = mask + self._mask_view = mask[...] self._ham = ham self._wfn = wfn self._ci_op = ci_op @@ -220,6 +257,7 @@ def __init__( # Set read-only flag on public array attributes self._pspace.setflags(write=False) self._sspace.setflags(write=False) + self._mask_view.setflags(write=False) def optimize( self, @@ -249,17 +287,30 @@ def optimize( Result of optimization. """ + # Check if system is underdetermined + if self.nequation < self.nactive: + raise ValueError("system is underdetermined") + # Convert x0 to proper dtype array x0 = np.asarray(x0, dtype=pyci.c_double) # Check input x0 length if x0.size != self.nparam: raise ValueError("length of `x0` does not match `param`") - # Use bare functions - f = self.compute_objective - j = self.compute_jacobian + # Prepare objective, Jacobian, x0 + if self.nactive < self.nparam: + # Generate objective, Jacobian, x0 with frozen parameters + x_ref = np.copy(x0) + f = self.mask_function(self.compute_objective, x_ref) + j = self.mask_function(self.compute_jacobian, x_ref) + x0 = np.copy(x0[self.mask]) + else: + # Use bare functions + f = self.compute_objective + j = self.compute_jacobian # Set up initial arguments to optimizer + opt_args = f, x0 opt_kwargs = kwargs.copy() if use_jac: opt_kwargs["jac"] = j @@ -319,6 +370,7 @@ def optimize_stochastic( nocc_up = self._wfn.nocc_up nocc_dn = self._wfn.nocc_dn constraints = self._constraints + mask = self._mask ci_cls = self._wfn.__class__ # Start at sample 1 @@ -361,6 +413,7 @@ def optimize_stochastic( nproj, nparam, constraints=constraints, + mask=mask, fill=fill, ) @@ -401,6 +454,36 @@ def remove_constraint(self, name: str) -> None: # Update nequation self._nequation = self._nproj + len(self._constraints) + def freeze_parameter(self, *params: Sequence[int]) -> None: + r""" + Set a FanCI parameter to be frozen during optimization. + + Parameters + ---------- + params : Sequence[int] + Indices of parameters to freeze. + + """ + for param in params: + self._mask[param] = False + # Update nactive + self._nactive = self._mask.sum() + + def unfreeze_parameter(self, *params: Sequence[int]) -> None: + r""" + Set a FanCI parameter to be active during optimization. + + Parameters + ---------- + params : Sequence[int] + Indices of parameters to unfreeze. + + """ + for param in params: + self._mask[param] = True + # Update nactive + self._nactive = self._mask.sum() + def compute_objective(self, x: np.ndarray) -> np.ndarray: r""" Compute the FanCI objective function. @@ -469,7 +552,7 @@ def compute_jacobian(self, x: np.ndarray) -> np.ndarray: """ # Allocate Jacobian matrix (in transpose memory order) - jac = np.empty((self._nequation, self._nparam), order="F", dtype=pyci.c_double) + jac = np.empty((self._nequation, self._nactive), order="F", dtype=pyci.c_double) jac_proj = jac[: self._nproj] jac_cons = jac[self._nproj :] @@ -486,17 +569,20 @@ def compute_jacobian(self, x: np.ndarray) -> np.ndarray: # d_ovlp = self.compute_overlap_deriv(x[:-1], "S") - # Compute final Jacobian column - # - # dE/d(p_k) = dE/d(p_k) \delta_{nk} c_n - # - ovlp = self.compute_overlap(x[:-1], "P") - ovlp *= -1 - jac_proj[:, -1] = ovlp - # - # Remove final column from jac_proj - # - jac_proj = jac_proj[:, :-1] + # Check is energy parameter is active: + if self._mask[-1]: + # + # Compute final Jacobian column if mask[-1] == True + # + # dE/d(p_k) = dE/d(p_k) \delta_{nk} c_n + # + ovlp = self.compute_overlap(x[:-1], "P") + ovlp *= -1 + jac_proj[:, -1] = ovlp + # + # Remove final column from jac_proj + # + jac_proj = jac_proj[:, :-1] # Iterate over remaining columns of Jacobian and d_ovlp for jac_col, d_ovlp_col in zip(jac_proj.transpose(), d_ovlp.transpose()): @@ -545,15 +631,16 @@ def f(x: np.ndarray) -> float: Constraint function p_{i} - v_{i}. """ - return x[i] - val + return x[i] - val if self._mask[i] else 0 def dfdx(x: np.ndarray) -> np.ndarray: r""" Constraint gradient \delta_{ki}. """ - y = np.zeros(self.nparam, dtype=pyci.c_double) - y[i] = 1 + y = np.zeros(self.nactive, dtype=pyci.c_double) + if self._mask[i]: + y[self._mask[:i].sum()] = 1 return y return f, dfdx @@ -590,13 +677,44 @@ def dfdx(x: np.ndarray) -> np.ndarray: Constraint gradient d(<\psi_{i}|\Psi>)/d(p_{k}). """ - y = np.zeros(self._nparam, dtype=pyci.c_double) + y = np.zeros(self._nactive, dtype=pyci.c_double) d_ovlp = self.compute_overlap_deriv(x[:-1], self._sspace[np.newaxis, i])[0] - y[: -1] = d_ovlp + y[: self._nactive - self._mask[-1]] = d_ovlp return y return f, dfdx + def mask_function(self, f: Callable, x_ref: np.ndarray) -> Callable: + r""" + Generate masked function for optimization with frozen parameters. + + Parameters + ---------- + f : Callable + Initial function. + x_ref : np.ndarray + Full parameter vector including frozen terms. + + Returns + ------- + f : Callable + Masked function. + + """ + if self.nactive == self.nparam: + return f + + def masked_f(x: np.ndarray) -> Any: + r""" + Masked function. + + """ + y = np.copy(x_ref) + y[self._mask] = x + return f(y) + + return masked_f + @abstractmethod def compute_overlap(self, x: np.ndarray, occs_array: Union[np.ndarray, str]) -> np.ndarray: r""" From df40695d5d2fcaa4e5babcdf3854d47ae6aed609 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Mon, 15 Apr 2024 12:51:14 -0500 Subject: [PATCH 04/11] Revert "Added masks support in FanCI class" This reverts commit d6e7c4b84c688f1b9c90e41ed832ba70a01f3b2e. --- pyci/fanci/fanci.py | 170 +++++++------------------------------------- 1 file changed, 26 insertions(+), 144 deletions(-) diff --git a/pyci/fanci/fanci.py b/pyci/fanci/fanci.py index a296ffb..4d81ac7 100644 --- a/pyci/fanci/fanci.py +++ b/pyci/fanci/fanci.py @@ -4,12 +4,18 @@ """ from abc import ABCMeta, abstractmethod + from collections import OrderedDict + from typing import Any, Callable, Dict, List, Sequence, Tuple, Union + import numpy as np + from scipy.optimize import OptimizeResult, least_squares, root + import pyci + __all__ = [ "FanCI", ] @@ -44,15 +50,6 @@ def nparam(self) -> int: """ return self._nparam - - @property - def nactive(self) -> int: - r""" - Number of active parameters for this FanCI problem. - - """ - return self._nactive - @property def constraints(self) -> Tuple[str]: r""" @@ -61,14 +58,6 @@ def constraints(self) -> Tuple[str]: """ return tuple(self._constraints.keys()) - @property - def mask(self) -> np.ndarray: - r""" - Frozen parameter mask. - - """ - return self._mask_view - @property def ham(self) -> pyci.hamiltonian: r""" @@ -158,7 +147,6 @@ def __init__( norm_param: Sequence[Tuple[int, float]] = None, norm_det: Sequence[Tuple[int, float]] = None, constraints: Dict[str, Tuple[Callable, Callable]] = None, - mask: Sequence[int] = None, fill: str = "excitation", ) -> None: r""" @@ -182,10 +170,6 @@ def __init__( them. constraints : Dict[str, Tuple[Callable, Callable]], optional Pairs of functions (f, dfdx) corresponding to additional constraints. - mask : Sequence[int] or Sequence[bool], optional - List of parameters to freeze. If the list contains ints, then each element corresponds - to a frozen parameter. If the list contains bools, then each element indicates whether - that parameter is active (True) or frozen (False). fill : ('excitation' | 'seniority' | None) Whether to fill the projection ("P") space by excitation level, by seniority, or not at all (in which case ``wfn`` must already be filled). @@ -209,26 +193,8 @@ def __init__( name = f"<\\psi_{{{index}}}|\\Psi> - v_{{{index}}}" constraints[name] = self.make_det_constraint(index, value) - # Generate mask (True for active parameter, False for frozen parameter) - if mask is None: - mask = np.ones(nparam, dtype=np.bool) - else: - mask = np.array(mask) - if mask.dtype == np.bool: - # Check length of boolean mask - if mask.size != nparam: - raise ValueError(f"Mask size is {mask.size}; must have size `nparam={nparam}`") - else: - # Convert integer mask to boolean - ints = mask - mask = np.ones(nparam, dtype=np.bool) - mask[ints] = 0 - # Number of nonlinear equations and active parameters nequation = nproj + len(constraints) - nactive = mask.sum() - if nequation < nactive: - print(f"WARNING: System is underdetermined with dimensions {nequation}, {nactive}. Continuing anyways") # Generate determinant spaces wfn = fill_wavefunction(wfn, nproj, fill) @@ -244,10 +210,7 @@ def __init__( self._nequation = nequation self._nproj = nproj self._nparam = nparam - self._nactive = nactive self._constraints = constraints - self._mask = mask - self._mask_view = mask[...] self._ham = ham self._wfn = wfn self._ci_op = ci_op @@ -257,7 +220,6 @@ def __init__( # Set read-only flag on public array attributes self._pspace.setflags(write=False) self._sspace.setflags(write=False) - self._mask_view.setflags(write=False) def optimize( self, @@ -287,30 +249,17 @@ def optimize( Result of optimization. """ - # Check if system is underdetermined - if self.nequation < self.nactive: - raise ValueError("system is underdetermined") - # Convert x0 to proper dtype array x0 = np.asarray(x0, dtype=pyci.c_double) # Check input x0 length if x0.size != self.nparam: raise ValueError("length of `x0` does not match `param`") - # Prepare objective, Jacobian, x0 - if self.nactive < self.nparam: - # Generate objective, Jacobian, x0 with frozen parameters - x_ref = np.copy(x0) - f = self.mask_function(self.compute_objective, x_ref) - j = self.mask_function(self.compute_jacobian, x_ref) - x0 = np.copy(x0[self.mask]) - else: - # Use bare functions - f = self.compute_objective - j = self.compute_jacobian + # Use bare functions + f = self.compute_objective + j = self.compute_jacobian # Set up initial arguments to optimizer - opt_args = f, x0 opt_kwargs = kwargs.copy() if use_jac: opt_kwargs["jac"] = j @@ -370,7 +319,6 @@ def optimize_stochastic( nocc_up = self._wfn.nocc_up nocc_dn = self._wfn.nocc_dn constraints = self._constraints - mask = self._mask ci_cls = self._wfn.__class__ # Start at sample 1 @@ -413,7 +361,6 @@ def optimize_stochastic( nproj, nparam, constraints=constraints, - mask=mask, fill=fill, ) @@ -454,36 +401,6 @@ def remove_constraint(self, name: str) -> None: # Update nequation self._nequation = self._nproj + len(self._constraints) - def freeze_parameter(self, *params: Sequence[int]) -> None: - r""" - Set a FanCI parameter to be frozen during optimization. - - Parameters - ---------- - params : Sequence[int] - Indices of parameters to freeze. - - """ - for param in params: - self._mask[param] = False - # Update nactive - self._nactive = self._mask.sum() - - def unfreeze_parameter(self, *params: Sequence[int]) -> None: - r""" - Set a FanCI parameter to be active during optimization. - - Parameters - ---------- - params : Sequence[int] - Indices of parameters to unfreeze. - - """ - for param in params: - self._mask[param] = True - # Update nactive - self._nactive = self._mask.sum() - def compute_objective(self, x: np.ndarray) -> np.ndarray: r""" Compute the FanCI objective function. @@ -552,7 +469,7 @@ def compute_jacobian(self, x: np.ndarray) -> np.ndarray: """ # Allocate Jacobian matrix (in transpose memory order) - jac = np.empty((self._nequation, self._nactive), order="F", dtype=pyci.c_double) + jac = np.empty((self._nequation, self._nparam), order="F", dtype=pyci.c_double) jac_proj = jac[: self._nproj] jac_cons = jac[self._nproj :] @@ -569,20 +486,17 @@ def compute_jacobian(self, x: np.ndarray) -> np.ndarray: # d_ovlp = self.compute_overlap_deriv(x[:-1], "S") - # Check is energy parameter is active: - if self._mask[-1]: - # - # Compute final Jacobian column if mask[-1] == True - # - # dE/d(p_k) = dE/d(p_k) \delta_{nk} c_n - # - ovlp = self.compute_overlap(x[:-1], "P") - ovlp *= -1 - jac_proj[:, -1] = ovlp - # - # Remove final column from jac_proj - # - jac_proj = jac_proj[:, :-1] + # Compute final Jacobian column + # + # dE/d(p_k) = dE/d(p_k) \delta_{nk} c_n + # + ovlp = self.compute_overlap(x[:-1], "P") + ovlp *= -1 + jac_proj[:, -1] = ovlp + # + # Remove final column from jac_proj + # + jac_proj = jac_proj[:, :-1] # Iterate over remaining columns of Jacobian and d_ovlp for jac_col, d_ovlp_col in zip(jac_proj.transpose(), d_ovlp.transpose()): @@ -631,16 +545,15 @@ def f(x: np.ndarray) -> float: Constraint function p_{i} - v_{i}. """ - return x[i] - val if self._mask[i] else 0 + return x[i] - val def dfdx(x: np.ndarray) -> np.ndarray: r""" Constraint gradient \delta_{ki}. """ - y = np.zeros(self.nactive, dtype=pyci.c_double) - if self._mask[i]: - y[self._mask[:i].sum()] = 1 + y = np.zeros(self.nparam, dtype=pyci.c_double) + y[i] = 1 return y return f, dfdx @@ -677,44 +590,13 @@ def dfdx(x: np.ndarray) -> np.ndarray: Constraint gradient d(<\psi_{i}|\Psi>)/d(p_{k}). """ - y = np.zeros(self._nactive, dtype=pyci.c_double) + y = np.zeros(self._nparam, dtype=pyci.c_double) d_ovlp = self.compute_overlap_deriv(x[:-1], self._sspace[np.newaxis, i])[0] - y[: self._nactive - self._mask[-1]] = d_ovlp + y[: -1] = d_ovlp return y return f, dfdx - def mask_function(self, f: Callable, x_ref: np.ndarray) -> Callable: - r""" - Generate masked function for optimization with frozen parameters. - - Parameters - ---------- - f : Callable - Initial function. - x_ref : np.ndarray - Full parameter vector including frozen terms. - - Returns - ------- - f : Callable - Masked function. - - """ - if self.nactive == self.nparam: - return f - - def masked_f(x: np.ndarray) -> Any: - r""" - Masked function. - - """ - y = np.copy(x_ref) - y[self._mask] = x - return f(y) - - return masked_f - @abstractmethod def compute_overlap(self, x: np.ndarray, occs_array: Union[np.ndarray, str]) -> np.ndarray: r""" From bb285bd763f192b704c2b1ddcb36c9efb9f23169 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Mon, 15 Apr 2024 14:33:35 -0500 Subject: [PATCH 05/11] Remova mask support in FanPT Wrapper --- pyci/fanci/fanci.py | 5 ----- pyci/fanci/fanpt_wrapper.py | 18 ++++-------------- 2 files changed, 4 insertions(+), 19 deletions(-) diff --git a/pyci/fanci/fanci.py b/pyci/fanci/fanci.py index 4d81ac7..c4fe7f5 100644 --- a/pyci/fanci/fanci.py +++ b/pyci/fanci/fanci.py @@ -4,15 +4,10 @@ """ from abc import ABCMeta, abstractmethod - from collections import OrderedDict - from typing import Any, Callable, Dict, List, Sequence, Tuple, Union - import numpy as np - from scipy.optimize import OptimizeResult, least_squares, root - import pyci diff --git a/pyci/fanci/fanpt_wrapper.py b/pyci/fanci/fanpt_wrapper.py index 4e1ff9f..150fdec 100644 --- a/pyci/fanci/fanpt_wrapper.py +++ b/pyci/fanci/fanpt_wrapper.py @@ -80,12 +80,6 @@ def solve_fanpt( # Select FANPT method if energy_active: fanptcontainer = FANPTContainerEParam - if not fanci_wfn.mask[-1]: - fanci_wfn.unfreeze_parameter(-1) - else: - fanptcontainer = FANPTContainerEFree - if fanci_wfn.mask[-1]: - fanci_wfn.freeze_parameter(-1) if resum: if energy_active: @@ -107,8 +101,7 @@ def solve_fanpt( # Get initial guess for parameters at initial lambda value. results = fanci_wfn.optimize(guess_params, **solver_kwargs) - guess_params[fanci_wfn.mask] = results.x - params = guess_params + params = results.x # Solve FANPT equations for l in np.linspace(lambda_i, lambda_f, steps, endpoint=False): @@ -151,9 +144,7 @@ def solve_fanpt( # Take the params given by fanci and use them as initial params in the # fanpt calculation for the next lambda. results = fanci_wfn.optimize(fanpt_params, **solver_kwargs) - - fanpt_params[fanci_wfn.mask] = results.x - params = fanpt_params + params = results.x if not energy_active: fanci_wfn.freeze_parameter([-1]) @@ -167,7 +158,6 @@ def solve_fanpt( # "steps": steps, # "inorm": inorm, # "norm_det": norm_det, - # "eactive": fanci_wfn.mask[-1], # "final_l": final_l, # }, #} @@ -187,7 +177,7 @@ def update_fanci_wfn(ham, fanciwfn, norm_det, fill): # FIXME: for FanCI class #return fanci_class( - # ham, fanciwfn.wfn, fanciwfn.nproj, fanciwfn.nparam, norm_det=norm_det, mask=fanciwfn.mask, fill=fill + # ham, fanciwfn.wfn, fanciwfn.nproj, fanciwfn.nparam, norm_det=norm_det, fill=fill #) # for fanpy class @@ -196,7 +186,7 @@ def update_fanci_wfn(ham, fanciwfn, norm_det, fill): ham, fanciwfn._fanpy_wfn, fanciwfn._wfn.nocc_up + fanciwfn._wfn.nocc_dn, nproj=fanciwfn.nproj, wfn=fanciwfn.wfn, fill=fill, seniority=fanciwfn.seniority, step_print=fanciwfn.step_print, step_save=fanciwfn.step_save, tmpfile=fanciwfn.tmpfile, - mask=fanciwfn._mask, objective_type=fanciwfn.objective_type, norm_det=norm_det, + objective_type=fanciwfn.objective_type, norm_det=norm_det, param_selection=fanciwfn.indices_component_params, constraints=fanciwfn._constraints ) From 111120c4b0eac698a4360581090893093fe7d823 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Tue, 28 May 2024 13:47:28 -0500 Subject: [PATCH 06/11] Mask support removed --- pyci/fanci/ap1rog.py | 4 ++-- pyci/fanci/apig.py | 4 ++-- pyci/fanpt/base_fanpt_container.py | 6 ------ 3 files changed, 4 insertions(+), 10 deletions(-) diff --git a/pyci/fanci/ap1rog.py b/pyci/fanci/ap1rog.py index a6607d6..d593f8f 100644 --- a/pyci/fanci/ap1rog.py +++ b/pyci/fanci/ap1rog.py @@ -109,7 +109,7 @@ def __init__( idx_param_cons=idx_param_cons, param_cons=param_cons, ) - def compute_overlap(self, x: np.ndarray) -> np.ndarray: + def compute_overlap(self, x: np.ndarray, occs_array = None) -> np.ndarray: r""" Compute the FanCI overlap vector. @@ -126,7 +126,7 @@ def compute_overlap(self, x: np.ndarray) -> np.ndarray: """ return self._cext.overlap(x) - def compute_overlap_deriv(self, x: np.ndarray) -> np.ndarray: + def compute_overlap_deriv(self, x: np.ndarray, occs_array = None) -> np.ndarray: r""" Compute the FanCI overlap derivative matrix. diff --git a/pyci/fanci/apig.py b/pyci/fanci/apig.py index 95715e8..589820a 100644 --- a/pyci/fanci/apig.py +++ b/pyci/fanci/apig.py @@ -104,7 +104,7 @@ def __init__( idx_param_cons=idx_param_cons, param_cons=param_cons, ) - def compute_overlap(self, x: np.ndarray) -> np.ndarray: + def compute_overlap(self, x: np.ndarray, occs_array = None) -> np.ndarray: r""" Compute the FanCI overlap vector. @@ -121,7 +121,7 @@ def compute_overlap(self, x: np.ndarray) -> np.ndarray: """ return self._cext.overlap(x) - def compute_overlap_deriv(self, x: np.ndarray) -> np.ndarray: + def compute_overlap_deriv(self, x: np.ndarray, occs_array = None) -> np.ndarray: r""" Compute the FanCI overlap derivative matrix. diff --git a/pyci/fanpt/base_fanpt_container.py b/pyci/fanpt/base_fanpt_container.py index 78bbc4c..c7bd7f3 100644 --- a/pyci/fanpt/base_fanpt_container.py +++ b/pyci/fanpt/base_fanpt_container.py @@ -50,10 +50,6 @@ class FANPTContainer(metaclass=ABCMeta): Wavefunction parameters. energy : float Energy for the current value of lambda. - active_energy : bool - Indicates if the energy will be varied in the calculations. - It is False either when the energy is frozen in a E-param calculation - or in any E-free calculation. ham_ci_op : pyci.sparse_op PyCI sparse operator corresponding to the perturbed Hamiltonian. f_pot_ci_op : pyci.sparse_op @@ -141,7 +137,6 @@ def __init__( self.params = params self.wfn_params = params[:-1] self.energy = params[-1] - self.active_energy = fanci_wfn.mask[-1] self.inorm = inorm # Assign ideal and real Hamiltonians. @@ -167,7 +162,6 @@ def __init__( else: self.f_pot = FANPTContainer.linear_comb_ham(self.ham1, self.ham0, 1.0, -1.0) self.f_pot_ci_op = pyci.sparse_op(self.f_pot, self.fanci_wfn.wfn, self.fanci_wfn.nproj) - if ovlp_s: self.ovlp_s = ovlp_s else: From cb8db5efe499d1c76e728b6c4620392b7fed8ffa Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Mon, 3 Jun 2024 11:42:59 -0500 Subject: [PATCH 07/11] Solved conflicts with nactive parameter in FanPT container Added Control of Energy as Active Parameter Small corrections active_energy argument full removed Revert "active_energy argument full removed" This reverts commit 1bc07c2e18b19bea68e8cd4e77792d7904888213. Revert "Small corrections" This reverts commit 712b527f63e4317795d209204de4cfea1eb49dec. Revert "Added Control of Energy as Active Parameter" This reverts commit 9623fd46095343b4f1ba9004f5faca60bd66517d. Revert "Solved conflicts with nactive parameter in FanPT container" This reverts commit 3ab67bd08ac64f15074608bb6253b5dfc8038182. Removing nactive parameter Corrected number of columns Small bug solved Mask removed in fanpt_updater Updating update_fanci_wfn Updating update_fanci_wfn to PyCI Corrected update_fanci_wfn More updates in update_fanci_wfn Solving class conflicts Solving bug Corrections to active_energy support Bug solved on gen_coeff_matrix --- pyci/fanci/fanpt_wrapper.py | 28 ++++++---------------------- pyci/fanpt/base_fanpt_container.py | 26 +++++++++++++++----------- pyci/fanpt/fanpt_cont_e_free.py | 5 ++++- pyci/fanpt/fanpt_cont_e_param.py | 9 +++++++-- pyci/fanpt/fanpt_updater.py | 5 ++--- 5 files changed, 34 insertions(+), 39 deletions(-) diff --git a/pyci/fanci/fanpt_wrapper.py b/pyci/fanci/fanpt_wrapper.py index 150fdec..9371148 100644 --- a/pyci/fanci/fanpt_wrapper.py +++ b/pyci/fanci/fanpt_wrapper.py @@ -80,6 +80,8 @@ def solve_fanpt( # Select FANPT method if energy_active: fanptcontainer = FANPTContainerEParam + else: + fanptcontainer = FANPTContainerEFree if resum: if energy_active: @@ -87,11 +89,11 @@ def solve_fanpt( "The energy parameter must be inactive with the resumation option." ) nequation = fanci_wfn.nequation - nactive = fanci_wfn.nactive + nparams = len(fanci_wfn.wfn_params) steps = 1 - if not inorm and (nequation == nactive): + if not inorm and (nequation == nparams): norm_det = [(ref_sd, 1.0)] - elif inorm and (nequation - 1) == nactive: + elif inorm and (nequation - 1) == nparams: fanci_wfn.remove_constraint(f"<\\psi_{{{ref_sd}}}|\\Psi> - v_{{{ref_sd}}}") inorm = False else: @@ -146,9 +148,6 @@ def solve_fanpt( results = fanci_wfn.optimize(fanpt_params, **solver_kwargs) params = results.x - if not energy_active: - fanci_wfn.freeze_parameter([-1]) - # Output for debugging purposes results["energy"] = fanpt_params[-1] results["residuals"] = results.fun @@ -172,23 +171,8 @@ def update_fanci_wfn(ham, fanciwfn, norm_det, fill): else: nocc = fanciwfn.wfn.nocc_up - # Activate energy parameter - fanciwfn.unfreeze_parameter([-1]) - - # FIXME: for FanCI class - #return fanci_class( - # ham, fanciwfn.wfn, fanciwfn.nproj, fanciwfn.nparam, norm_det=norm_det, fill=fill - #) - - # for fanpy class - # ASSUMING GeneratedFanCI class return fanci_class( - ham, fanciwfn._fanpy_wfn, fanciwfn._wfn.nocc_up + fanciwfn._wfn.nocc_dn, - nproj=fanciwfn.nproj, wfn=fanciwfn.wfn, fill=fill, seniority=fanciwfn.seniority, - step_print=fanciwfn.step_print, step_save=fanciwfn.step_save, tmpfile=fanciwfn.tmpfile, - objective_type=fanciwfn.objective_type, norm_det=norm_det, - param_selection=fanciwfn.indices_component_params, - constraints=fanciwfn._constraints + ham, nocc, fanciwfn.nproj, fanciwfn.wfn ) diff --git a/pyci/fanpt/base_fanpt_container.py b/pyci/fanpt/base_fanpt_container.py index c7bd7f3..950f919 100644 --- a/pyci/fanpt/base_fanpt_container.py +++ b/pyci/fanpt/base_fanpt_container.py @@ -50,6 +50,10 @@ class FANPTContainer(metaclass=ABCMeta): Wavefunction parameters. energy : float Energy for the current value of lambda. + active_energy : bool + Indicates if the energy will be varied in the calculations. + It is False either when the energy is frozen in a E-param calculation + or in any E-free calculation. ham_ci_op : pyci.sparse_op PyCI sparse operator corresponding to the perturbed Hamiltonian. f_pot_ci_op : pyci.sparse_op @@ -99,6 +103,7 @@ def __init__( params, ham0, ham1, + active_energy, l=0, ref_sd=0, inorm=False, @@ -119,6 +124,10 @@ def __init__( PyCI Hamiltonian of the ideal system. ham1 : pyci.hamiltonian PyCI Hamiltonian of the real system. + active_energy : bool + Indicates if the energy will be varied in the calculations. + It is False either when the energy is frozen in a E-param calculation + or in any E-free calculation. l : float Lambda value. ref_sd : int @@ -139,6 +148,12 @@ def __init__( self.energy = params[-1] self.inorm = inorm + self.active_energy = active_energy + if active_energy: + self.nactive = len(params) + else: + self.nactive = len(params) - 1 + # Assign ideal and real Hamiltonians. self.ham1 = ham1 self.ham0 = ham0 @@ -215,17 +230,6 @@ def linear_comb_ham(ham1, ham0, a1, a0): two_mo = a1 * ham1.two_mo + a0 * ham0.two_mo return pyci.hamiltonian(ecore, one_mo, two_mo) - @property - def nactive(self): - r"""Return the number of active parameters. - - Returns - ------- - self.fanci_wfn.nactive : int - Number of active parameters. - """ - return self.fanci_wfn.nactive - @property def nequation(self): r"""Return the number of equations. diff --git a/pyci/fanpt/fanpt_cont_e_free.py b/pyci/fanpt/fanpt_cont_e_free.py index 3585a72..8250060 100644 --- a/pyci/fanpt/fanpt_cont_e_free.py +++ b/pyci/fanpt/fanpt_cont_e_free.py @@ -108,6 +108,7 @@ def __init__( params, ham0, ham1, + active_energy=False, l=0, ref_sd=0, inorm=False, @@ -142,7 +143,7 @@ def __init__( d_ovlp_s : {np.ndarray, None} Derivatives of the overlaps in the "S" projection space. """ - if fanci_wfn.mask[-1]: + if active_energy: raise TypeError("The energy cannot be an active parameter.") else: super().__init__( @@ -150,6 +151,7 @@ def __init__( params, ham0, ham1, + active_energy, l, ref_sd, inorm, @@ -244,6 +246,7 @@ def gen_coeff_matrix(self): numpy array with shape (self.nequations, len(self.nactive)). """ super().gen_coeff_matrix() + self.c_matrix = self.c_matrix[:, :self.nactive] f_proj = np.empty((self.nproj, self.nactive), order="F") for f_proj_col, d_ovlp_col in zip(f_proj.transpose(), self.d_ovlp_s.transpose()): self.ham_ci_op(d_ovlp_col, out=f_proj_col) diff --git a/pyci/fanpt/fanpt_cont_e_param.py b/pyci/fanpt/fanpt_cont_e_param.py index 0ccb538..44a7ea1 100644 --- a/pyci/fanpt/fanpt_cont_e_param.py +++ b/pyci/fanpt/fanpt_cont_e_param.py @@ -110,6 +110,7 @@ def __init__( params, ham0, ham1, + active_energy=True, l=0, ref_sd=0, inorm=False, @@ -144,11 +145,13 @@ def __init__( d_ovlp_s : {np.ndarray, None} Derivatives of the overlaps in the "S" projection space. """ + super().__init__( fanci_wfn, params, ham0, ham1, + active_energy, l, ref_sd, inorm, @@ -185,14 +188,16 @@ def der2_g_lambda_wfnparams(self): d2_g_lambda_wfnparams : np.ndarray Derivative of the FANPT equations with respect to lambda and the wavefunction parameters. - numpy array with shape (self.nequations, len(self.wfn_params_active)). + numpy array with shape (self.nequations, self.nactive). """ if self.active_energy: ncolumns = self.nactive - 1 else: ncolumns = self.nactive + f = np.zeros((self.nequation, ncolumns), order="F") f_proj = f[: self.nproj] + for f_proj_col, d_ovlp_col in zip(f_proj.transpose(), self.d_ovlp_s.transpose()): self.f_pot_ci_op(d_ovlp_col, out=f_proj_col) self.d2_g_lambda_wfnparams = f @@ -208,7 +213,7 @@ def der2_g_e_wfnparams(self): d2_g_e_wfnparams : np.ndarray Derivative of the FANPT equations with respect to the energy and the wavefunction parameters. - numpy array with shape (self.nequations, len(self.wfn_params_active)). + numpy array with shape (self.nequations, self.nactive). """ if self.active_energy: ncolumns = self.nactive - 1 diff --git a/pyci/fanpt/fanpt_updater.py b/pyci/fanpt/fanpt_updater.py index 4d52f71..14e0ccc 100644 --- a/pyci/fanpt/fanpt_updater.py +++ b/pyci/fanpt/fanpt_updater.py @@ -311,9 +311,8 @@ def params_updater(self): else: wfn_responses = self.responses corrections = np.sum(wfn_responses * dl.reshape(self.final_order, 1), axis=0) - active_wfn_indices = np.where(self.fanpt_container.fanci_wfn.mask[:-1])[0] - for c, active_index in zip(corrections, active_wfn_indices): - wfn_params[active_index] += c + for index, c in enumerate(corrections): + wfn_params[index] += c self.new_wfn_params = wfn_params def energy_ham_ovlp_updater(self): From 4f89f150176c4b3b73412b0c84ae9f2dff8b56b0 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Wed, 26 Jun 2024 16:35:32 -0500 Subject: [PATCH 08/11] Added tests for FanPT --- pyci/fanci/fanpt_wrapper.py | 24 +------ pyci/fanpt/base_fanpt_container.py | 2 +- pyci/test/test_fanpt.py | 105 +++++++++++++++++++++++++++++ 3 files changed, 109 insertions(+), 22 deletions(-) create mode 100644 pyci/test/test_fanpt.py diff --git a/pyci/fanci/fanpt_wrapper.py b/pyci/fanci/fanpt_wrapper.py index 9371148..0f6be5f 100644 --- a/pyci/fanci/fanpt_wrapper.py +++ b/pyci/fanci/fanpt_wrapper.py @@ -75,7 +75,7 @@ def solve_fanpt( norm_det = [(ref_sd, 1.0)] else: inorm = False - norm_det = None + norm_det = list() # Select FANPT method if energy_active: @@ -148,18 +148,8 @@ def solve_fanpt( results = fanci_wfn.optimize(fanpt_params, **solver_kwargs) params = results.x - # Output for debugging purposes - results["energy"] = fanpt_params[-1] results["residuals"] = results.fun - #output = { - # "x": params, - # "variables": { - # "steps": steps, - # "inorm": inorm, - # "norm_det": norm_det, - # "final_l": final_l, - # }, - #} + return results @@ -172,7 +162,7 @@ def update_fanci_wfn(ham, fanciwfn, norm_det, fill): nocc = fanciwfn.wfn.nocc_up return fanci_class( - ham, nocc, fanciwfn.nproj, fanciwfn.wfn + ham, nocc, fanciwfn.nproj, fanciwfn.wfn, norm_det=norm_det, fill=fill, ) @@ -190,13 +180,5 @@ def reduce_to_fock(two_int, lambda_val=0): indices = np.arange(nspatial) fock_two_int[indices[:, None, None], indices[None, :, None], indices[None, None, :], indices[None, :, None]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, None, :], indices[None, :, None]] fock_two_int[indices[:, None, None], indices[None, :, None], indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, :, None], indices[None, None, :]] - #fock_two_int[indices[None, :, None], indices[:, None, None], indices[None, None, :], indices[None, :, None]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, None, :], indices[None, :, None]] - #fock_two_int[indices[None, :, None], indices[:, None, None], indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, :, None], indices[None, None, :]] - - #occ_indices = np.arange(nelec // 2) - #fock_two_int[indices[:, None, None], occ_indices[None, :, None], indices[None, None, :], occ_indices[None, :, None]] = two_int[indices[:, None, None], occ_indices[None, :, None], indices[None, None, :], occ_indices[None, :, None]] - #fock_two_int[indices[:, None, None], occ_indices[None, :, None], occ_indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], occ_indices[None, :, None], occ_indices[None, :, None], indices[None, None, :]] - #fock_two_int[occ_indices[None, :, None], indices[:, None, None], indices[None, None, :], occ_indices[None, :, None]] = two_int[indices[:, None, None], occ_indices[None, :, None], indices[None, None, :], occ_indices[None, :, None]] - #fock_two_int[occ_indices[None, :, None], indices[:, None, None], occ_indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], occ_indices[None, :, None], occ_indices[None, :, None], indices[None, None, :]] return fock_two_int diff --git a/pyci/fanpt/base_fanpt_container.py b/pyci/fanpt/base_fanpt_container.py index 950f919..8425e97 100644 --- a/pyci/fanpt/base_fanpt_container.py +++ b/pyci/fanpt/base_fanpt_container.py @@ -195,7 +195,7 @@ def __init__( self.ref_sd = ref_sd else: raise KeyError( - "The normalization of the Slater determinant is not constrained" + "The normalization of the Slater determinant is not constrained " "in the FanCI wavefunction." ) else: diff --git a/pyci/test/test_fanpt.py b/pyci/test/test_fanpt.py new file mode 100644 index 0000000..d9f62b8 --- /dev/null +++ b/pyci/test/test_fanpt.py @@ -0,0 +1,105 @@ +# This file is part of PyCI. +# +# PyCI is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your +# option) any later version. +# +# PyCI is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License +# along with PyCI. If not, see . + +import numpy as np +import pyci +from scipy.special import comb + +import pytest +from pyci.test import datafile +from pyci.fanci import AP1roG + +from pyci.fanci.fanpt_wrapper import reduce_to_fock, solve_fanpt +from pyci.fanpt import FANPTUpdater, FANPTContainerEParam, FANPTContainerEFree + +@pytest.mark.parametrize("filename, nocc, expected", [("he_ccpvqz", 1, -2.8868091056425156), + ("be_ccpvdz", 2, -14.600556820761211), + ("li2_ccpvdz", 3, -16.862861409549044), + ("lih_sto6g", 2, -8.963531095653355), + ("h2_631gdp", 1, -1.869682842154122), + ("h2o_ccpvdz", 5, -77.96987451399201)]) +def test_fanpt_e_free(filename, nocc, expected): + nsteps = 10 + order = 1 + + # Define ham0 and ham1 + ham1 = pyci.hamiltonian(datafile("{0:s}.fcidump".format(filename))) + ham0 = pyci.hamiltonian(ham1.ecore, ham1.one_mo, reduce_to_fock(ham1.two_mo)) + + # contruct empty fci wave function class instance from # of basis functions and occupation + wfn0 = pyci.fullci_wfn(ham0.nbasis, nocc, nocc) + wfn0.add_hartreefock_det() + + # initialize sparse matrix operator (hamiltonian into wave function) + op = pyci.sparse_op(ham0, wfn0) + + # solve for the lowest eigenvalue and eigenvector + e_hf, e_vecs0 = op.solve(n=1, tol=1.0e-8) + + # Get params as the solution of the fanci wfn with ham0 (last element will be the energy of the "ideal" system). + nproj = int(comb(ham1.nbasis, nocc)) + pyci_wfn = AP1roG(ham1, nocc, nproj=nproj) + + params = np.zeros(pyci_wfn.nparam, dtype=pyci.c_double) + params[-1] = e_hf[0] + + fill = 'excitation' + fanpt_results = solve_fanpt(pyci_wfn, ham0, pyci_wfn.ham, params, + fill=fill, energy_active=False, resum=False, ref_sd=0, + final_order=order, lambda_i=0.0, lambda_f=1.0, steps=nsteps, + solver_kwargs={'mode':'lstsq', 'use_jac':True, 'xtol':1.0e-8, + 'ftol':1.0e-8, 'gtol':1.0e-5, 'max_nfev':pyci_wfn.nparam, 'verbose':2}) + + assert np.allclose(fanpt_results.x[-1], expected) + +@pytest.mark.parametrize("filename, nocc, expected", [("he_ccpvqz", 1, -2.8868091056425156), + ("be_ccpvdz", 2, -14.600556842700215), + ("li2_ccpvdz", 3, -16.86286984124269), + ("lih_sto6g", 2, -8.96353109432708), + ("h2_631gdp", 1, -1.869682842154122), + ("h2o_ccpvdz", 5, -77.96987516399848)]) +def test_fanpt_e_param(filename, nocc, expected): + nsteps = 10 + order = 1 + + # Define ham0 and ham1 + ham1 = pyci.hamiltonian(datafile("{0:s}.fcidump".format(filename))) + ham0 = pyci.hamiltonian(ham1.ecore, ham1.one_mo, reduce_to_fock(ham1.two_mo)) + + # contruct empty fci wave function class instance from # of basis functions and occupation + wfn0 = pyci.fullci_wfn(ham0.nbasis, nocc, nocc) + wfn0.add_hartreefock_det() + + # initialize sparse matrix operator (hamiltonian into wave function) + op = pyci.sparse_op(ham0, wfn0) + + # solve for the lowest eigenvalue and eigenvector + e_hf, e_vecs0 = op.solve(n=1, tol=1.0e-8) + + # Get params as the solution of the fanci wfn with ham0 (last element will be the energy of the "ideal" system). + nproj = int(comb(ham1.nbasis, nocc)) + pyci_wfn = AP1roG(ham1, nocc, nproj=nproj) + + params = np.zeros(pyci_wfn.nparam, dtype=pyci.c_double) + params[-1] = e_hf[0] + + fill = 'excitation' + fanpt_results = solve_fanpt(pyci_wfn, ham0, pyci_wfn.ham, params, + fill=fill, energy_active=True, resum=False, ref_sd=0, + final_order=order, lambda_i=0.0, lambda_f=1.0, steps=nsteps, + solver_kwargs={'mode':'lstsq', 'use_jac':True, 'xtol':1.0e-8, + 'ftol':1.0e-8, 'gtol':1.0e-5, 'max_nfev':pyci_wfn.nparam, 'verbose':2}) + + assert np.allclose(fanpt_results.x[-1], expected) From c8f7ac50892f8fce3a1870609efdb59ae077a340 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Mon, 15 Jul 2024 14:42:08 -0500 Subject: [PATCH 09/11] Support to any number of nproj Added flags to support the use of a different number of projections than the full space. --- pyci/fanpt/base_fanpt_container.py | 4 ++-- pyci/fanpt/fanpt_updater.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/pyci/fanpt/base_fanpt_container.py b/pyci/fanpt/base_fanpt_container.py index 8425e97..8d64a8e 100644 --- a/pyci/fanpt/base_fanpt_container.py +++ b/pyci/fanpt/base_fanpt_container.py @@ -171,12 +171,12 @@ def __init__( if ham_ci_op: self.ham_ci_op = ham_ci_op else: - self.ham_ci_op = pyci.sparse_op(self.ham, self.fanci_wfn.wfn, self.fanci_wfn.nproj) + self.ham_ci_op = pyci.sparse_op(self.ham, self.fanci_wfn.wfn, self.fanci_wfn.nproj, symmetric=False) if f_pot_ci_op: self.f_pot_ci_op = f_pot_ci_op else: self.f_pot = FANPTContainer.linear_comb_ham(self.ham1, self.ham0, 1.0, -1.0) - self.f_pot_ci_op = pyci.sparse_op(self.f_pot, self.fanci_wfn.wfn, self.fanci_wfn.nproj) + self.f_pot_ci_op = pyci.sparse_op(self.f_pot, self.fanci_wfn.wfn, self.fanci_wfn.nproj, symmetric=False) if ovlp_s: self.ovlp_s = ovlp_s else: diff --git a/pyci/fanpt/fanpt_updater.py b/pyci/fanpt/fanpt_updater.py index 14e0ccc..2a1ec25 100644 --- a/pyci/fanpt/fanpt_updater.py +++ b/pyci/fanpt/fanpt_updater.py @@ -335,7 +335,7 @@ def energy_ham_ovlp_updater(self): self.fanpt_container.ham1, self.fanpt_container.ham0, self.final_l, 1 - self.final_l ) new_ham_op = pyci.sparse_op( - new_ham, self.fanpt_container.fanci_wfn.wfn, self.fanpt_container.nproj + new_ham, self.fanpt_container.fanci_wfn.wfn, self.fanpt_container.nproj, symmetric=False ) new_ovlp_s = self.fanpt_container.fanci_wfn.compute_overlap(self.new_wfn_params, "S") f = np.empty(self.fanpt_container.nproj, dtype=pyci.c_double) From f89743445bc3c685d5b0845660a07b3761ff80e7 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Thu, 8 Aug 2024 12:52:19 -0400 Subject: [PATCH 10/11] Removed optimized guess FanPT guess must be provide by the user --- pyci/fanci/fanpt_wrapper.py | 54 ++++++++++++++++++++++++++----------- 1 file changed, 38 insertions(+), 16 deletions(-) diff --git a/pyci/fanci/fanpt_wrapper.py b/pyci/fanci/fanpt_wrapper.py index 0f6be5f..9eff81d 100644 --- a/pyci/fanci/fanpt_wrapper.py +++ b/pyci/fanci/fanpt_wrapper.py @@ -6,11 +6,12 @@ from .fanci import FanCI from ..fanpt import FANPTUpdater, FANPTContainerEParam, FANPTContainerEFree + def solve_fanpt( fanci_wfn, ham0, ham1, - guess_params, + params, fill, energy_active=True, resum=False, @@ -27,7 +28,7 @@ def solve_fanpt( Args: fanci_wfn : FanCI class FanCI wavefunction. - guess_params : np.ndarray + params : np.ndarray Initial guess for wave function parameters. ham0 : pyci.hamiltonian PyCI Hamiltonian of the ideal system. @@ -85,9 +86,7 @@ def solve_fanpt( if resum: if energy_active: - raise ValueError( - "The energy parameter must be inactive with the resumation option." - ) + raise ValueError("The energy parameter must be inactive with the resumation option.") nequation = fanci_wfn.nequation nparams = len(fanci_wfn.wfn_params) steps = 1 @@ -97,13 +96,11 @@ def solve_fanpt( fanci_wfn.remove_constraint(f"<\\psi_{{{ref_sd}}}|\\Psi> - v_{{{ref_sd}}}") inorm = False else: - raise ValueError( - "The necesary condition of a determined system of equations is not met." - ) + raise ValueError("The necesary condition of a determined system of equations is not met.") # Get initial guess for parameters at initial lambda value. - results = fanci_wfn.optimize(guess_params, **solver_kwargs) - params = results.x + numerical_zero = 1e-12 + params = np.where(params == 0, numerical_zero, params) # Solve FANPT equations for l in np.linspace(lambda_i, lambda_f, steps, endpoint=False): @@ -119,7 +116,7 @@ def solve_fanpt( ) final_l = l + (lambda_f - lambda_i) / steps - print(f'Solving FanPT problem at lambda={final_l}') + print(f"Solving FanPT problem at lambda={final_l}") fanpt_updater = FANPTUpdater( fanpt_container=fanpt_container, @@ -133,8 +130,8 @@ def solve_fanpt( # These params serve as initial guess to solve the fanci equations for the given lambda. fanpt_params = np.append(new_wfn_params, new_energy) - print('Frobenius Norm of parameters: {}'.format(np.linalg.norm(fanpt_params - params))) - print('Energy change: {}'.format(np.linalg.norm(fanpt_params[-1] - params[-1]))) + print("Frobenius Norm of parameters: {}".format(np.linalg.norm(fanpt_params - params))) + print("Energy change: {}".format(np.linalg.norm(fanpt_params[-1] - params[-1]))) # Initialize perturbed Hamiltonian with the current value of lambda using the static method of fanpt_container. ham = fanpt_container.linear_comb_ham(ham1, ham0, final_l, 1 - final_l) @@ -162,7 +159,12 @@ def update_fanci_wfn(ham, fanciwfn, norm_det, fill): nocc = fanciwfn.wfn.nocc_up return fanci_class( - ham, nocc, fanciwfn.nproj, fanciwfn.wfn, norm_det=norm_det, fill=fill, + ham, + nocc, + fanciwfn.nproj, + fanciwfn.wfn, + norm_det=norm_det, + fill=fill, ) @@ -178,7 +180,27 @@ def reduce_to_fock(two_int, lambda_val=0): fock_two_int = two_int * lambda_val nspatial = two_int.shape[0] indices = np.arange(nspatial) - fock_two_int[indices[:, None, None], indices[None, :, None], indices[None, None, :], indices[None, :, None]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, None, :], indices[None, :, None]] - fock_two_int[indices[:, None, None], indices[None, :, None], indices[None, :, None], indices[None, None, :]] = two_int[indices[:, None, None], indices[None, :, None], indices[None, :, None], indices[None, None, :]] + fock_two_int[ + indices[:, None, None], + indices[None, :, None], + indices[None, None, :], + indices[None, :, None], + ] = two_int[ + indices[:, None, None], + indices[None, :, None], + indices[None, None, :], + indices[None, :, None], + ] + fock_two_int[ + indices[:, None, None], + indices[None, :, None], + indices[None, :, None], + indices[None, None, :], + ] = two_int[ + indices[:, None, None], + indices[None, :, None], + indices[None, :, None], + indices[None, None, :], + ] return fock_two_int From 914a4cd2408ee7381af91599653218b573347f10 Mon Sep 17 00:00:00 2001 From: "Carlos E. V. de Moura" Date: Mon, 12 Aug 2024 11:04:56 -0400 Subject: [PATCH 11/11] Improved tests --- pyci/test/test_fanpt.py | 104 +++++++++++++++++++++++++++++----------- 1 file changed, 76 insertions(+), 28 deletions(-) diff --git a/pyci/test/test_fanpt.py b/pyci/test/test_fanpt.py index d9f62b8..7e317d8 100644 --- a/pyci/test/test_fanpt.py +++ b/pyci/test/test_fanpt.py @@ -24,12 +24,18 @@ from pyci.fanci.fanpt_wrapper import reduce_to_fock, solve_fanpt from pyci.fanpt import FANPTUpdater, FANPTContainerEParam, FANPTContainerEFree -@pytest.mark.parametrize("filename, nocc, expected", [("he_ccpvqz", 1, -2.8868091056425156), - ("be_ccpvdz", 2, -14.600556820761211), - ("li2_ccpvdz", 3, -16.862861409549044), - ("lih_sto6g", 2, -8.963531095653355), - ("h2_631gdp", 1, -1.869682842154122), - ("h2o_ccpvdz", 5, -77.96987451399201)]) + +@pytest.mark.parametrize( + "filename, nocc, expected", + [ + ("he_ccpvqz", 1, -2.8868091056425156), + ("be_ccpvdz", 2, -14.600556820761211), + ("li2_ccpvdz", 3, -16.862861409549044), + ("lih_sto6g", 2, -8.963531095653355), + ("h2_631gdp", 1, -1.869682842154122), + ("h2o_ccpvdz", 5, -77.96987451399201), + ], +) def test_fanpt_e_free(filename, nocc, expected): nsteps = 10 order = 1 @@ -53,23 +59,47 @@ def test_fanpt_e_free(filename, nocc, expected): pyci_wfn = AP1roG(ham1, nocc, nproj=nproj) params = np.zeros(pyci_wfn.nparam, dtype=pyci.c_double) - params[-1] = e_hf[0] - - fill = 'excitation' - fanpt_results = solve_fanpt(pyci_wfn, ham0, pyci_wfn.ham, params, - fill=fill, energy_active=False, resum=False, ref_sd=0, - final_order=order, lambda_i=0.0, lambda_f=1.0, steps=nsteps, - solver_kwargs={'mode':'lstsq', 'use_jac':True, 'xtol':1.0e-8, - 'ftol':1.0e-8, 'gtol':1.0e-5, 'max_nfev':pyci_wfn.nparam, 'verbose':2}) + params[-1] = e_hf[0] - ham1.ecore + + fill = "excitation" + fanpt_results = solve_fanpt( + pyci_wfn, + ham0, + pyci_wfn.ham, + params, + fill=fill, + energy_active=False, + resum=False, + ref_sd=0, + final_order=order, + lambda_i=0.0, + lambda_f=1.0, + steps=nsteps, + solver_kwargs={ + "mode": "lstsq", + "use_jac": True, + "xtol": 1.0e-8, + "ftol": 1.0e-8, + "gtol": 1.0e-5, + "max_nfev": pyci_wfn.nparam, + "verbose": 2, + }, + ) assert np.allclose(fanpt_results.x[-1], expected) -@pytest.mark.parametrize("filename, nocc, expected", [("he_ccpvqz", 1, -2.8868091056425156), - ("be_ccpvdz", 2, -14.600556842700215), - ("li2_ccpvdz", 3, -16.86286984124269), - ("lih_sto6g", 2, -8.96353109432708), - ("h2_631gdp", 1, -1.869682842154122), - ("h2o_ccpvdz", 5, -77.96987516399848)]) + +@pytest.mark.parametrize( + "filename, nocc, expected", + [ + ("he_ccpvqz", 1, -2.8868091056425156), + ("be_ccpvdz", 2, -14.600556842700215), + ("li2_ccpvdz", 3, -16.86286984124269), + ("lih_sto6g", 2, -8.96353109432708), + ("h2_631gdp", 1, -1.869682842154122), + ("h2o_ccpvdz", 5, -77.96987516399848), + ], +) def test_fanpt_e_param(filename, nocc, expected): nsteps = 10 order = 1 @@ -93,13 +123,31 @@ def test_fanpt_e_param(filename, nocc, expected): pyci_wfn = AP1roG(ham1, nocc, nproj=nproj) params = np.zeros(pyci_wfn.nparam, dtype=pyci.c_double) - params[-1] = e_hf[0] - - fill = 'excitation' - fanpt_results = solve_fanpt(pyci_wfn, ham0, pyci_wfn.ham, params, - fill=fill, energy_active=True, resum=False, ref_sd=0, - final_order=order, lambda_i=0.0, lambda_f=1.0, steps=nsteps, - solver_kwargs={'mode':'lstsq', 'use_jac':True, 'xtol':1.0e-8, - 'ftol':1.0e-8, 'gtol':1.0e-5, 'max_nfev':pyci_wfn.nparam, 'verbose':2}) + params[-1] = e_hf[0] - ham1.ecore + + fill = "excitation" + fanpt_results = solve_fanpt( + pyci_wfn, + ham0, + pyci_wfn.ham, + params, + fill=fill, + energy_active=True, + resum=False, + ref_sd=0, + final_order=order, + lambda_i=0.0, + lambda_f=1.0, + steps=nsteps, + solver_kwargs={ + "mode": "lstsq", + "use_jac": True, + "xtol": 1.0e-8, + "ftol": 1.0e-8, + "gtol": 1.0e-5, + "max_nfev": pyci_wfn.nparam, + "verbose": 2, + }, + ) assert np.allclose(fanpt_results.x[-1], expected)