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/__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/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/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 new file mode 100644 index 0000000..9eff81d --- /dev/null +++ b/pyci/fanci/fanpt_wrapper.py @@ -0,0 +1,206 @@ +""" 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, + 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. + 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 = list() + + # Select FANPT method + if energy_active: + fanptcontainer = FANPTContainerEParam + else: + fanptcontainer = FANPTContainerEFree + + if resum: + if energy_active: + raise ValueError("The energy parameter must be inactive with the resumation option.") + nequation = fanci_wfn.nequation + nparams = len(fanci_wfn.wfn_params) + steps = 1 + if not inorm and (nequation == nparams): + norm_det = [(ref_sd, 1.0)] + elif inorm and (nequation - 1) == nparams: + 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. + 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): + 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) + params = results.x + + results["residuals"] = results.fun + + 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 + + return fanci_class( + ham, + nocc, + fanciwfn.nproj, + fanciwfn.wfn, + norm_det=norm_det, + fill=fill, + ) + + +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, :], + ] + + 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..8d64a8e --- /dev/null +++ b/pyci/fanpt/base_fanpt_container.py @@ -0,0 +1,296 @@ +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, + active_energy, + 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. + 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 + 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.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 + + # 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, 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, symmetric=False) + 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 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..8250060 --- /dev/null +++ b/pyci/fanpt/fanpt_cont_e_free.py @@ -0,0 +1,262 @@ +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, + active_energy=False, + 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 active_energy: + raise TypeError("The energy cannot be an active parameter.") + else: + super().__init__( + fanci_wfn, + params, + ham0, + ham1, + active_energy, + 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() + 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) + 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..44a7ea1 --- /dev/null +++ b/pyci/fanpt/fanpt_cont_e_param.py @@ -0,0 +1,241 @@ +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, + active_energy=True, + 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, + active_energy, + 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, 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 + + 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, self.nactive). + """ + 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..2a1ec25 --- /dev/null +++ b/pyci/fanpt/fanpt_updater.py @@ -0,0 +1,373 @@ +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) + for index, c in enumerate(corrections): + wfn_params[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, 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) + 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) diff --git a/pyci/test/test_fanpt.py b/pyci/test/test_fanpt.py new file mode 100644 index 0000000..7e317d8 --- /dev/null +++ b/pyci/test/test_fanpt.py @@ -0,0 +1,153 @@ +# 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] - 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), + ], +) +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] - 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) 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 = {