From 46f04e5967fe6633570820a120470dfabf3aa1e5 Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Thu, 9 Feb 2023 10:03:34 -0500 Subject: [PATCH 1/8] updated --- pycc/ccwfn.py | 10 +- pycc/lccwfn.py | 10 +- pycc/utils.py | 142 ++++++++++++++++++++++ pycc/utils_copy.py | 290 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 445 insertions(+), 7 deletions(-) create mode 100644 pycc/utils_copy.py diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py index 347e26a..53ce27f 100644 --- a/pycc/ccwfn.py +++ b/pycc/ccwfn.py @@ -10,11 +10,11 @@ import time import numpy as np import torch -from .utils import helper_diis, cc_contract -from .hamiltonian import Hamiltonian -from .local import Local -from .cctriples import t_tjl, t3c_ijk -from .lccwfn import lccwfn +from utils import helper_diis, cc_contract +from hamiltonian import Hamiltonian +from local import Local +from cctriples import t_tjl, t3c_ijk +from lccwfn import lccwfn class ccwfn(object): """ diff --git a/pycc/lccwfn.py b/pycc/lccwfn.py index 8a31551..bc1b826 100644 --- a/pycc/lccwfn.py +++ b/pycc/lccwfn.py @@ -2,7 +2,7 @@ #from timer import Timer import numpy as np from opt_einsum import contract - +from utils import helper_ldiis class lccwfn(object): """ @@ -80,7 +80,7 @@ def __init__(self, o, v, no, nv, H, local, model, eref, Local): self.t1_ii = t1_ii self.t2_ij = t2_ij - def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100): + def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=2, start_diis=1): """ Parameters ---------- @@ -111,6 +111,8 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100): #self.r2_t = Timer("r2") #self.energy_t = Timer("energy") + ldiis = helper_ldiis(self.no,self.t1_ii, self.t2_ij, max_diis) + elcc = self.lcc_energy(self.Local.Fov_ij,self.Local.Loovv_ij,self.t1_ii, self.t2_ij) print("CC Iter %3d: lCC Ecorr = %.15f dE = % .5E MP2" % (0,elcc,-elcc)) @@ -148,6 +150,10 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100): #print(Timer.timers) return elcc + ldiis.add_error_vector(self.no,self.t1_ii,self.t2_ij) + if niter >= start_diis: + self.t1_ii, self.t2_ij = ldiis.extrapolate(self.no,self.t1_ii, self.t2_ij) + def local_residuals(self, t1_ii, t2_ij): """ Constructing the two- and four-index intermediates diff --git a/pycc/utils.py b/pycc/utils.py index 807c4c1..1a30b1c 100644 --- a/pycc/utils.py +++ b/pycc/utils.py @@ -38,6 +38,7 @@ def add_error_vector(self, t1, t2): # Add DIIS vectors self.diis_vals_t1.append(t1.copy()) self.diis_vals_t2.append(t2.copy()) + # Add new error vectors error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() @@ -138,6 +139,147 @@ def extrapolate(self, t1, t2): return t1, t2 +class helper_ldiis(object): + def __init__(self, no,t1_ii, t2_ij, max_diis): + oldt1 = [] + oldt2 = [] + diis_vals_t1 = [] + diis_vals_t2 = [] + for i in range(no): + oldt1.append(t1_ii[i].copy()) + for j in range(no): + ij = i*no + j + oldt2.append(t2_ij[ij].copy()) + + self.oldt1 = oldt1 + self.oldt2 = oldt2 + diis_vals_t1.append(self.oldt1.copy()) + diis_vals_t2.append(self.oldt2.copy()) + self.diis_vals_t1 = diis_vals_t1 + self.diis_vals_t2 = diis_vals_t2 + #print("diis t2 initial", self.diis_vals_t2) + self.diis_errors_t1 = [] + self.diis_errors_t2 = [] + self.diis_size = 0 + self.max_diis = max_diis + + def add_error_vector(self, no,t1_ii, t2_ij): + # Add DIIS vectors + self.addt1 = [] + self.addt2 = [] + + for i in range(no): + ii = i*no + i + + self.addt1.append(t1_ii[i]) + + for j in range(no): + ij = i*no + j + + self.addt2.append(t2_ij[ij]) + #if ij == 24: + #print("t2 compare", ij, t2_ij[ij]) + + self.diis_vals_t1.append(self.addt1.copy()) + self.diis_vals_t2.append(self.addt2.copy()) + #for n1, e1 in enumerate(self.diis_vals_t2): + #print("n1", n1, "e1", e1[24]) + + #print("oldt2", self.oldt2[ij]) + + # Add new error vectors + error_t1 = [] + error_t2 = [] + for i in range(no): + error_t1.append((self.diis_vals_t1[-1][i] - self.oldt1[i]).ravel()) + for j in range(no): + ij = i*no + j + + error_t2.append((self.diis_vals_t2[-1][ij] - self.oldt2[ij]).ravel()) + #if ij == 24: + #print("diis t2 last", ij, self.diis_vals_t2[-1][ij]) + #print("oldt2 last", ij, self.oldt2[ij]) + #print(error_t2[24]) + self.diis_errors_t1.append(error_t1) + self.diis_errors_t2.append(error_t2) + self.oldt1 = self.addt1.copy() + self.oldt2 = self.addt2.copy() + + def extrapolate(self, no,t1_ii, t2_ij): + + if (self.max_diis == 0): + return t1_ii, t2_ij + + # Limit size of DIIS vector + if (len(self.diis_errors_t1) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors_t1[0] + del self.diis_errors_t2[0] + + self.diis_size_t1 = len(self.diis_errors_t1) + self.diis_size_t2 = len(self.diis_errors_t2) + B_t1 = np.ones((self.diis_size_t1 + 1, self.diis_size_t1 + 1)) * -1 + B_t1[-1, -1] = 0 + B_t2 = np.ones((self.diis_size_t2 + 1, self.diis_size_t2 + 1)) * -1 + B_t2[-1, -1] = 0 + + t1_ii = np.zeros_like(self.oldt1) + t2_ij = np.zeros_like(self.oldt2) + for i in range(no): + + for n1, e1 in enumerate(self.diis_errors_t1): + B_t1[n1, n1] = np.dot(e1[i], e1[i]) + for n2, e2 in enumerate(self.diis_errors_t1): + if n1 >= n2: + continue + B_t1[n1, n2] = np.dot(e1[i], e2[i]) + B_t1[n2, n1] = B_t1[n1, n2] + B_t1[:-1, :-1] /= np.abs(B_t1[:-1, :-1]).max() + + # Build residual vector + resid = np.zeros(self.diis_size_t1 + 1) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B_t1, resid) + + # Calculate new amplitudes + for num in range(self.diis_size_t1): + t1_ii[i] += ci[num] * self.diis_vals_t1[num + 1][i] + + # Save extrapolated amplitudes to old_t amplitudes + for j in range(no): + ij = i*no + j + + for n1, e1 in enumerate(self.diis_errors_t2): + B_t2[n1, n1] = np.dot(e1[ij], e1[ij]) + for n2, e2 in enumerate(self.diis_errors_t2): + if n1 >= n2: + continue + B_t2[n1, n2] = np.dot(e1[ij], e2[ij]) + B_t2[n2, n1] = B_t2[n1, n2] + + B_t2[:-1, :-1] /= np.abs(B_t2[:-1, :-1]).max() + + # Build residual vector + resid = np.zeros(self.diis_size_t2 + 1) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B_t2, resid) + #print("ci", ci, "ci.shape", ci.shape) + + # Calculate new amplitudes + for num in range(self.diis_size_t2): + t2_ij[ij] += ci[num] * self.diis_vals_t2[num + 1][ij] + #print("old", self.oldt2[24]) + #print("new", t2_ij[24]) + self.oldt1 = t1_ii.copy() + self.oldt2 = t2_ij.copy() + + return t1_ii, t2_ij + class cc_contract(object): """ A wrapper for opt_einsum.contract with tensors stored on CPU and/or GPU. diff --git a/pycc/utils_copy.py b/pycc/utils_copy.py new file mode 100644 index 0000000..e904ca2 --- /dev/null +++ b/pycc/utils_copy.py @@ -0,0 +1,290 @@ +import numpy as np +import torch +import opt_einsum + + +class helper_diis(object): + def __init__(self, t1, t2, max_diis, precision='DP'): + if isinstance(t1, torch.Tensor): + self.device0 = torch.device('cpu') + self.device1 = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + self.diis_vals_t1 = [t1.clone()] + self.diis_vals_t2 = [t2.clone()] + else: + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + print(self.oldt1.shape) + print(self.oldt1) + print(self.oldt2.shape) + self.diis_vals_t1 = [t1.copy()] + self.diis_vals_t2 = [t2.copy()] + + self.diis_errors = [] + self.diis_size = 0 + self.max_diis = max_diis + self.precision = precision + + def add_error_vector(self, t1, t2): + if isinstance(t1, torch.Tensor): + # Add DIIS vectors + self.diis_vals_t1.append(t1.clone()) + self.diis_vals_t2.append(t2.clone()) + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(torch.cat((error_t1, error_t2))) + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + else: + # Add DIIS vectors + self.diis_vals_t1.append(t1.copy()) + self.diis_vals_t2.append(t2.copy()) + + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(np.concatenate((error_t1, error_t2))) + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + + def extrapolate(self, t1, t2): + + if (self.max_diis == 0): + return t1, t2 + + # Limit size of DIIS vector + if (len(self.diis_errors) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors[0] + + self.diis_size = len(self.diis_errors) + + if isinstance(t1, torch.Tensor): + # Build error matrix B + if self.precision == 'DP': + B = torch.ones((self.diis_size + 1, self.diis_size + 1), dtype=torch.float64, device=self.device1) * -1 + elif self.precision == 'SP': + B = torch.ones((self.diis_size + 1, self.diis_size + 1), dtype=torch.float32, device=self.device1) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = torch.dot(e1, e1) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = torch.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= torch.abs(B[:-1, :-1]).max() + + # Build residual vector + if self.precision == 'DP': + resid = torch.zeros((self.diis_size + 1), dtype=torch.float64, device=self.device1) + elif self.precision == 'SP': + resid = torch.zeros((self.diis_size + 1), dtype=torch.float32, device=self.device1) + resid[-1] = -1 + + # Solve pulay equations + ci = torch.linalg.solve(B, resid) + + # Calculate new amplitudes + t1 = torch.zeros_like(self.oldt1) + t2 = torch.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += torch.real(ci[num] * self.diis_vals_t1[num + 1]) + t2 += torch.real(ci[num] * self.diis_vals_t2[num + 1]) + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + + else: + # Build error matrix B + if self.precision == 'DP': + B = np.ones((self.diis_size + 1, self.diis_size + 1)) * -1 + elif self.precision == 'SP': + B = np.ones((self.diis_size + 1, self.diis_size + 1), dtype=np.float32) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = np.dot(e1, e1) + print("diis_errors",self.diis_errors[n1]) + print("B", B.shape) + print("n1","e1",n1,e1.shape) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = np.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= np.abs(B[:-1, :-1]).max() + + # Build residual vector + if self.precision == 'DP': + resid = np.zeros(self.diis_size + 1) + elif self.precision == 'SP': + resid = np.zeros((self.diis_size + 1), dtype=np.float32) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B, resid) + + # Calculate new amplitudes + t1 = np.zeros_like(self.oldt1) + t2 = np.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += ci[num] * self.diis_vals_t1[num + 1] + t2 += ci[num] * self.diis_vals_t2[num + 1] + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + + return t1, t2 + +class helper_ldiis(object): + def __init__(self, no,t1_ii, t2_ij, max_diis): + oldt1 = [] + oldt2 = [] + for i in range(no): + + oldt1.append(np.array(t1_ii[i],dtype="object")) + for j in range(no): + ij = i*no + j + + oldt2.append(np.array(t2_ij[ij],dtype="object")) + + self.oldt1 = oldt1 + self.oldt2 = oldt2 + print(self.oldt1) + self.diis_vals_t1 = [np.array(self.oldt1,dtype="object")] + self.diis_vals_t2 = [np.array(self.oldt2,dtype="object")] + self.diis_errors = [] + self.diis_size = 0 + self.max_diis = max_diis + + def add_error_vector(self, no,t1_ii, t2_ij): + # Add DIIS vectors + self.addt1 = [] + self.addt2 = [] + for i in range(no): + ii = i*no + i + + self.addt1.append(np.array(t1_ii[i],dtype="object")) + for j in range(no): + ij = i*no + j + + self.addt2.append(np.array(t2_ij[ij],dtype="object")) + + self.diis_vals_t1.append(np.array(self.addt1,dtype="object")) + self.diis_vals_t2.append(np.array(self.addt2,dtype="object")) + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(np.concatenate((error_t1, error_t2))) + self.oldt1 = self.addt1.copy() + self.oldt2 = self.addt2.copy() + + def extrapolate(self, t1_ii, t2_ij): + + if (self.max_diis == 0): + return t1_ii, t2_ij + + # Limit size of DIIS vector + if (len(self.diis_errors) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors[0] + + self.diis_size = len(self.diis_errors) + + B = np.ones((self.diis_size + 1, self.diis_size + 1)) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + print("B",B.shape ) + print("n1","e1", n1, e1.shape) + print("diis_errors",self.diis_errors[n1]) + B[n1, n1] = np.dot(e1, e1) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = np.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= np.abs(B[:-1, :-1]).max() + + # Build residual vector + resid = np.zeros(self.diis_size + 1) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B, resid) + + # Calculate new amplitudes + t1_ii = np.zeros_like(self.oldt1) + t2_ij = np.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1_ii += ci[num] * self.diis_vals_t1[num + 1] + t2_ij += ci[num] * self.diis_vals_t2[num + 1] + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1_ii.copy() + self.oldt2 = t2_ij.copy() + + t1_ii = list(np.float_(t1_ii))#float(t1_ii.tolist()) + t2_ij = list(np.float_(t2_ij))#float(t2_ij.tolist()) + + return t1_ii, t2_ij + +class cc_contract(object): + """ + A wrapper for opt_einsum.contract with tensors stored on CPU and/or GPU. + """ + def __init__(self, device='CPU'): + """ + Parameters + ---------- + device: string + initiated in ccwfn object, default: 'CPU' + + Returns + ------- + None + """ + self.device = device + if self.device == 'GPU': + # torch.device is an object representing the device on which torch.Tensor is or will be allocated. + self.device0 = torch.device('cpu') + self.device1 = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') + + def __call__(self, subscripts, *operands): + """ + Parameters + ---------- + subscripts: string + specify the subscripts for summation (same format as numpy.einsum) + *operands: list of array_like + the arrays/tensors for the operation + + Returns + ------- + An ndarray/torch.Tensor that is calculated based on Einstein summation convention. + """ + if self.device == 'CPU': + return opt_einsum.contract(subscripts, *operands) + elif self.device == 'GPU': + # Check the type and allocation of the tensors + # Transfer the copy from CPU to GPU if needed (for ERI) + input_list = list(operands) + for i in range(len(input_list)): + if (not input_list[i].is_cuda): + input_list[i] = input_list[i].to(self.device1) + #print(len(input_list), type(input_list[0]), type(input_list[1])) + output = opt_einsum.contract(subscripts, *input_list) + del input_list + return output + From 06bb5ce3c22a354c36631495709daee859594c96 Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Mon, 13 Mar 2023 11:11:10 -0400 Subject: [PATCH 2/8] in progress to add diis to local cc --- pycc/ccwfn.py | 3 +- pycc/lccwfn.py | 845 ++++++++++++++++++++++++++++++-------------- pycc/local.py | 233 ++++++------ pycc/utils.py | 17 +- pycc/utils_copy.py | 11 +- pycc/utils_copy1.py | 331 +++++++++++++++++ 6 files changed, 1056 insertions(+), 384 deletions(-) create mode 100644 pycc/utils_copy1.py diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py index 53ce27f..66dda29 100644 --- a/pycc/ccwfn.py +++ b/pycc/ccwfn.py @@ -151,7 +151,8 @@ def __init__(self, scf_wfn, **kwargs): if local is not None: self.Local = Local(local, self.C, self.nfzc, self.no, self.nv, self.H, self.local_cutoff,self.it2_opt) if filter is not True: - self.Local._trans_integrals(self.o, self.v) + self.Local.trans_integrals(self.o, self.v) + self.Local.overlaps(self.Local.QL) self.lccwfn = lccwfn(self.o, self.v,self.no, self.nv, self.H, self.local, self.model, self.eref, self.Local) # denominators diff --git a/pycc/lccwfn.py b/pycc/lccwfn.py index bc1b826..dcd2ed8 100644 --- a/pycc/lccwfn.py +++ b/pycc/lccwfn.py @@ -1,89 +1,87 @@ -import time +import time #from timer import Timer import numpy as np -from opt_einsum import contract +from opt_einsum import contract from utils import helper_ldiis class lccwfn(object): """ - A PNO-CC object. + A PNO-CC object. Attributes ---------- o : Numpy slice - occupied orbital subspace + occupied orbital subspace v : Numpy slice virtual orbital subspace no: int number of (active) occupied orbitals - nv: int + nv: int number of virtual orbitals H: PyCC Hamiltonian object Hamiltonian object with integrals, Fock matrices, and orbital coefficients - Local: PyCC Local object - Local object with transformation matrices, local Fock matrices and integrals, etc. - + Local: PyCC Local object + Local object with transformation matrices, local Fock matrices and integrals, etc. + Parameters ---------- - local: string + local: string type of local calculation ("PAO", "PNO", etc.) model: string type of CC calculation ("CCD","CCSD", etc.) eref: float - reference energy (typically HF/SCF energy) + reference energy (typically HF/SCF energy) Methods ------- - solve_lcc() + solve_lcc() Solves the local CC T amplitude equations local_residuals() - Computes the local T1 and T2 residuals for a given set of amplitudes and Fock operator - + Computes the local T1 and T2 residuals for a given set of amplitudes and Fock operator + Notes ----- - To do: + To do: (1) need DIIS extrapolation (2) time table for each intermediate? - (3) generate and store overlap terms prior to the calculation of the residuals - (4) remove redundant transformed integrals """ - - def __init__(self, o, v, no, nv, H, local, model, eref, Local): + + def __init__(self, o, v, no, nv, H, local, model, eref, Local): self.o = o self.v = v self.no = no self.nv = nv self.H = H self.local = local - self.model = model + self.model = model self.eref = eref - self.Local = Local + self.Local = Local self.QL = self.Local.QL self.dim = self.Local.dim - self.eps = self.Local.eps + self.eps = self.Local.eps - self.t1 = np.zeros((self.no, self.nv)) - t1_ii = [] - t2_ij = [] + t1 = [] + t2 = [] for i in range(self.no): ii = i*self.no + i - - t1_ii.append(self.QL[ii].T @ self.t1[i]) - + + t1.append(np.zeros((self.Local.dim[ii]))) + for j in range(self.no): ij = i*self.no + j - - t2_ij.append(-1* self.Local.ERIoovv_ij[ij][i,j] / (self.eps[ij].reshape(1,-1) + self.eps[ij].reshape(-1,1) - - self.H.F[i,i] - self.H.F[j,j])) - self.t1_ii = t1_ii - self.t2_ij = t2_ij + t2.append(-1* self.Local.ERIoovv[ij][i,j] / (self.eps[ij].reshape(1,-1) + self.eps[ij].reshape(-1,1) + - self.H.F[i,i] - self.H.F[j,j])) + + self.t1 = t1 + self.t2 = t2 - def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=2, start_diis=1): + def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis=1): """ Parameters ---------- + e_conv : float convergence condition for correlation energy (default if 1e-7) r_conv : float @@ -93,8 +91,8 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=2, start_dii Returns ------- - elcc: float - lCC correlation energy + elcc: float + lCC correlation energy """ lcc_tstart = time.time() @@ -110,33 +108,41 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=2, start_dii #self.r1_t = Timer("r1") #self.r2_t = Timer("r2") #self.energy_t = Timer("energy") - - ldiis = helper_ldiis(self.no,self.t1_ii, self.t2_ij, max_diis) - elcc = self.lcc_energy(self.Local.Fov_ij,self.Local.Loovv_ij,self.t1_ii, self.t2_ij) + ldiis = helper_ldiis(self.no, self.t1, self.t2, max_diis) + + elcc = self.lcc_energy(self.Local.Fov,self.Local.Loovv,self.t1, self.t2) print("CC Iter %3d: lCC Ecorr = %.15f dE = % .5E MP2" % (0,elcc,-elcc)) for niter in range(1, maxiter+1): elcc_last = elcc - r1_ii, r2_ij = self.local_residuals(self.t1_ii, self.t2_ij) + r1, r2 = self.local_residuals(self.t1, self.t2) rms = 0 rms_t1 = 0 rms_t2 = 0 for i in range(self.no): + ii = i*self.no + i + + #need to change to reshape + for a in range(self.Local.dim[ii]): + self.t1[i][a] += r1[i][a]/(self.H.F[i,i] - self.Local.eps[ii][a]) + + rms_t1 += contract('Z,Z->',r1[i],r1[i]) + for j in range(self.no): ij = i*self.no + j - self.t2_ij[ij] -= r2_ij[ij]/(self.eps[ij].reshape(1,-1) + self.eps[ij].reshape(-1,1) + self.t2[ij] -= r2[ij]/(self.eps[ij].reshape(1,-1) + self.eps[ij].reshape(-1,1) - self.H.F[i,i] - self.H.F[j,j]) - rms_t2 += contract('ZY,ZY->',r2_ij[ij],r2_ij[ij]) + rms_t2 += contract('ZY,ZY->',r2[ij],r2[ij]) rms = np.sqrt(rms_t2) - elcc = self.lcc_energy(self.Local.Fov_ij,self.Local.Loovv_ij,self.t1_ii, self.t2_ij) + elcc = self.lcc_energy(self.Local.Fov,self.Local.Loovv,self.t1, self.t2) ediff = elcc - elcc_last print("lCC Iter %3d: lCC Ecorr = %.15f dE = % .5E rms = % .5E" % (niter, elcc, ediff, rms)) @@ -150,27 +156,17 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=2, start_dii #print(Timer.timers) return elcc - ldiis.add_error_vector(self.no,self.t1_ii,self.t2_ij) + ldiis.add_error_vector(self.no, self.t1,self.t2) if niter >= start_diis: - self.t1_ii, self.t2_ij = ldiis.extrapolate(self.no,self.t1_ii, self.t2_ij) + self.t1, self.t2 = ldiis.extrapolate(self.no, self.t1, self.t2) - def local_residuals(self, t1_ii, t2_ij): + def local_residuals(self, t1, t2): """ Constructing the two- and four-index intermediates Then evaluating the singles and doubles residuals, storing them in a list of length occ (single) and length occ*occ (doubles) - Naming scheme same as integrals where _ij is attach as a suffix to the intermediate name - ... special naming scheme for those involving two different pair space ij and im -> _ijm To do ------ - There are some arguments that isn't really being used and will be updated once things are good to go - Listed here are intermediates with corresponding arguments that aren't needed since it requires "on the fly" generation - Fae_ij Lovvv_ij - Fme_ij (Fme_im) Loovv_ij - Wmbej_ijim ERIovvv_ij, ERIoovv_ij, Loovv_ij - Wmbje_ijim (Wmbie_ijmj) ERIovvv_ij, ERIoovv_ij - r1_ii Lovvo_ij, ERIovvo_ij - r2_ij ERIovvo_ij, ERIovov_ij, ERIvvvo_ij """ o = self.o v = self.v @@ -178,311 +174,612 @@ def local_residuals(self, t1_ii, t2_ij): L = self.H.L ERI = self.H.ERI - Fae_ij = [] - Fme_ij = [] - - Wmbej_ijim = [] - Wmbje_ijim = [] - Wmbie_ijmj = [] - Zmbij_ij = [] - - r1_ii = [] - r2_ij = [] - - Fae_ij = self.build_lFae(Fae_ij, self.Local.Fvv_ij, self.Local.Fov_ij, - self.Local.Lovvv_ij, self.Local.Loovv_ij, t1_ii, t2_ij) - lFmi = self.build_lFmi(o, F, self.Local.Fov_ij, self.Local.Looov_ij, self.Local.Loovv_ij, t1_ii, t2_ij) - Fme_ij = self.build_lFme(Fme_ij, self.Local.Fov_ij, self.Local.Loovv_ij, t1_ii) - lWmnij = self.build_lWmnij(o, ERI, self.Local.ERIooov_ij, self.Local.ERIoovo_ij, - self.Local.ERIoovv_ij, t1_ii, t2_ij) - Zmbij = self.build_lZmbij(Zmbij_ij, self.Local.ERIovvv_ij, t1_ii, t2_ij) - Wmbej_ijim = self.build_lWmbej(Wmbej_ijim, self.Local.ERIoovv_ij, self.Local.ERIovvo_ij, - self.Local.ERIovvv_ij,self.Local.ERIoovo_ij, self.Local.Loovv_ij, t1_ii, t2_ij) - Wmbje_ijim, Wmbie_ijmj = self.build_lWmbje(Wmbje_ijim, Wmbie_ijmj, self.Local.ERIovov_ij, - self.Local.ERIovvv_ij,self.Local.ERIoovv_ij, self.Local.ERIooov_ij, t1_ii, t2_ij) - - r1_ii = self.lr_T1(r1_ii, self.Local.Fov_ij , self.Local.ERIovvv_ij, self.Local.Lovvo_ij, self.Local.Loovo_ij, - t1_ii, t2_ij,Fae_ij, Fme_ij, lFmi) - r2_ij = self.lr_T2(r2_ij, self.Local.ERIoovv_ij, self.Local.ERIvvvv_ij, self.Local.ERIovvo_ij, self.Local.ERIovoo_ij, - self.Local.ERIvvvo_ij,self.Local.ERIovov_ij, t1_ii, t2_ij, Fae_ij ,lFmi,Fme_ij, lWmnij, Zmbij_ij, Wmbej_ijim, Wmbje_ijim, Wmbie_ijmj) - - return r1_ii, r2_ij - - def build_lFae(self, Fae_ij, Fvv_ij,Fov_ij, Lovvv_ij, Loovv_ij, t1_ii, t2_ij): + Fae = [] + Fme = [] + Wmbej = [] + Wmbje = [] + Wmbie = [] + Zmbij = [] + r1 = [] + r2 = [] + + Fae = self.build_Fae(Fae, L, self.Local.Fvv, self.Local.Fov, self.Local.Sijmm, self.Local.Sijmn, t1, t2) + Fmi = self.build_Fmi(o, F, L, self.Local.Fov, self.Local.Looov, self.Local.Loovv, t1, t2) + Fme = self.build_Fme(Fme, L, self.Local.Fov, t1) + Wmnij = self.build_Wmnij(o, ERI, self.Local.ERIooov, self.Local.ERIoovo, self.Local.ERIoovv, t1, t2) + Zmbij = self.build_Zmbij(Zmbij, ERI, self.Local.ERIovvv, t1, t2) + Wmbej = self.build_Wmbej(Wmbej, ERI, L, self.Local.ERIoovo, self.Local.Sijnn, self.Local.Sijnj, self.Local.Sijjn, t1, t2) + Wmbje, Wmbie = self.build_Wmbje(Wmbje, Wmbie, ERI, self.Local.ERIooov, self.Local.Sijnn, self.Local.Sijin, self.Local.Sijjn, t1, t2) + + r1 = self.r_T1(r1, self.Local.Fov , ERI, L, self.Local.Loovo, self.Local.Sijmm, self.Local.Sijim, self.Local.Sijmn, + t1, t2, Fae, Fmi, Fme) + r2 = self.r_T2(r2, ERI, self.Local.ERIoovv, self.Local.ERIvvvv, self.Local.ERIovoo, self.Local.Sijmm, self.Local.Sijim, + self.Local.Sijmj, self.Local.Sijnn, self.Local.Sijmn, t1, t2, Fae ,Fmi, Fme, Wmnij, Zmbij, Wmbej, Wmbje, Wmbie) + + return r1, r2 + + def build_Fae(self, Fae_ij, L, Fvv, Fov, Sijmm, Sijmn, t1, t2): #self.fae_t.start() o = self.o v = self.v QL = self.QL - for ij in range(self.no*self.no): - i = ij // self.no - j = ij % self.no + if self.model == 'CCD': + for ij in range(self.no*self.no): + i = ij // self.no + j = ij % self.no - Fae = Fvv_ij[ij].copy() + Fae = Fvv[ij].copy() - Fae_3 = np.zeros_like(Fae) - for m in range(self.no): - for n in range(self.no): - mn = m *self.no +n - nn = n*self.no + n + for m in range(self.no): + for n in range(self.no): + mn = m *self.no +n + ijmn = ij*(self.no**2) + mn + + tmp = Sijmn[ijmn] @ t2[mn] + tmp1 = QL[ij].T @ L[m,n,v,v] + tmp1 = tmp1 @ QL[mn] + Fae -= tmp @ tmp1.T + + Fae_ij.append(Fae) + else: + for ij in range(self.no*self.no): + i = ij // self.no + j = ij % self.no - Sijmn = QL[ij].T @ QL[mn] + Fae = Fvv[ij].copy() - tmp2 = Sijmn @ t2_ij[mn] - tmp3_0 = QL[ij].T @ self.H.L[m,n,v,v] - tmp3_1 = tmp3_0 @ QL[mn] - Fae_3 -= tmp2 @ tmp3_1.T + for m in range(self.no): + mm = m*self.no + m + ijm = ij*(self.no) + m - Fae_ij.append(Fae + Fae_3) + tmp = Sijmm[ijm] @ t1[m] + Fae -= 0.5* contract('e,a->ae',Fov[ij][m],tmp) + + tmp1 = contract('abc,aA->Abc',L[m,v,v,v], QL[ij]) + tmp1 = contract('Abc,bB->ABc',tmp1, QL[mm]) + tmp1 = contract('ABc,cC->ABC',tmp1, QL[ij]) + Fae += contract('F,aFe->ae',t1[m],tmp1) + + for n in range(self.no): + mn = m *self.no +n + nn = n*self.no + n + ijmn = ij*(self.no**2) + mn + + tmp2 = Sijmn[ijmn] @ t2[mn] + tmp3_0 = QL[ij].T @ L[m,n,v,v] + tmp3_1 = tmp3_0 @ QL[mn] + Fae -= tmp2 @ tmp3_1.T + + tmp4 = tmp3_0 @ QL[nn] + Fae -= 0.5 *contract('a,F,eF->ae', tmp, t1[n], tmp4) + + Fae_ij.append(Fae) #self.fae_t.stop() return Fae_ij - def build_lFmi(self, o, F, Fov_ij, Looov_ij, Loovv_ij, t1_ii, t2_ij): + def build_Fmi(self, o, F, L, Fov, Looov, Loovv, t1, t2): #self.fmi_t.start() v = self.v QL = self.QL Fmi = F[o,o].copy() - Fmi_3 = np.zeros_like(Fmi) - for j in range(self.no): - for n in range(self.no): - jn = j*self.no + n + if self.model == 'CCD': + for j in range(self.no): + for n in range(self.no): + jn = j*self.no + n + + Fmi[:,j] += contract('EF,mEF->m',t2[jn],Loovv[jn][:,n,:,:]) + else: + for j in range(self.no): + jj = j*self.no +j + for n in range(self.no): + jn = j*self.no + n + nn = n*self.no + n - Fmi_3[:,j] += contract('EF,mEF->m',t2_ij[jn],Loovv_ij[jn][:,n,:,:]) + Fmi[:,j] += 0.5 * contract('e,me->m', t1[j], Fov[jj]) + Fmi[:,j] += contract('e,me->m',t1[n],Looov[nn][:,n,j]) + Fmi[:,j] += contract('EF,mEF->m',t2[jn],Loovv[jn][:,n,:,:]) + + tmp = contract('mab,aA->mAb', L[o,n,v,v],QL[jj]) + tmp = contract('mAb,bB->mAB', tmp, QL[nn]) + Fmi[:,j] += 0.5 * contract('E,F,mEF->m',t1[j], t1[n], tmp) - Fmi_tot = Fmi + Fmi_3 #self.fmi_t.stop() - return Fmi_tot + return Fmi - def build_lFme(self, Fme_ij, Fov_ij, Loovv_ij, t1_ii): - """ - Implemented up to CCSD but debugging at the CCD level to narrow down terms - """ + def build_Fme(self, Fme_ij, L, Fov, t1): #self.fme_t.start() + QL = self.QL + v = self.v + + if self.model == 'CCD': + return + else: + for ij in range(self.no*self.no): + + Fme = Fov[ij].copy() + + for m in range(self.no): + for n in range(self.no): + nn = n*self.no + n + + tmp = QL[ij].T @ L[m,n,v,v] + tmp = tmp @ QL[nn] + Fme[m] += t1[n] @ tmp.T + + Fme_ij.append(Fme) + #self.fme_t.stop() - return + return Fme_ij - def build_lWmnij(self, o, ERI, ERIooov_ij, ERIoovo_ij, ERIoovv_ij, t1_ii, t2_ij): - """ - Implemented up to CCSD but debugging at the CCD level to narrow down terms - """ + def build_Wmnij(self, o, ERI, ERIooov, ERIoovo, ERIoovv, t1, t2): #self.wmnij_t.start() + v = self.v + QL = self.Local.QL Wmnij = ERI[o,o,o,o].copy() - Wmnij_3 = np.zeros_like(Wmnij) - for i in range(self.no): - for j in range(self.no): - ij = i*self.no + j + if self.model == 'CCD': + for i in range(self.no): + for j in range(self.no): + ij = i*self.no + j - Wmnij_3[:,:,i,j] += contract('ef,mnef->mn',t2_ij[ij], ERIoovv_ij[ij]) + Wmnij[:,:,i,j] += contract('ef,mnef->mn',t2[ij], ERIoovv[ij]) + else: + for i in range(self.no): + for j in range(self.no): + ij = i*self.no + j + ii = i*self.no + i + jj = j*self.no + j + + Wmnij[:,:,i,j] += contract('E,mnE->mn', t1[j], ERIooov[jj][:,:,i,:]) + Wmnij[:,:,i,j] += contract('E,mnE->mn', t1[i], ERIoovo[ii][:,:,:,j]) + Wmnij[:,:,i,j] += contract('ef,mnef->mn',t2[ij], ERIoovv[ij]) + + tmp = contract('aA,mnab->mnAb', QL[ii], ERI[o,o,v,v]) + tmp = contract('bB,mnAb->mnAB', QL[jj], tmp) + Wmnij[:,:,i,j] += contract('e,f,mnef->mn', t1[i], t1[j], tmp) - Wmnij_tot = Wmnij + Wmnij_3 #self.wmnij_t.stop() - return Wmnij_tot + return Wmnij - def build_lZmbij(self, Zmbij_ij, ERIovvv_ij, t1_ii, t2_ij): - """ - Implemented up to CCSD but debugging at the CCD level to narrow down terms - """ + def build_Zmbij(self, Zmbij_ij, ERI, ERIovvv, t1, t2): #self.zmbij_t.start() - #self.zmbij_t.stop() - return - - def build_lWmbej(self, Wmbej_ijim, ERIoovv_ij, ERIovvo_ij, ERIovvv_ij, ERIoovo_ij, Loovv_ij, t1_ii, t2_ij): - #self.wmbej_t.start() - v = self.v o = self.o + v = self.v QL = self.QL - dim = self.dim - for ij in range(self.no*self.no): - i = ij // self.no - j = ij % self.no - jj = j*self.no + j + if self.model == 'CCD': + return + else: + for ij in range(self.no*self.no): + i = ij // self.no + j = ij % self.no + ii = i*self.no + i + jj = j*self.no + j - for m in range(self.no): - im = i*self.no + m + Zmbij = np.zeros((self.no,self.Local.dim[ij])) - Wmbej = np.zeros((dim[ij],dim[im])) + Zmbij = contract('mbef,ef->mb', ERIovvv[ij], t2[ij]) - tmp = QL[ij].T @ self.H.ERI[m,v,v,j] - Wmbej = tmp @ QL[im] + tmp = contract('iabc,aA->iAbc',ERI[o,v,v,v], QL[ij]) + tmp = contract('iAbc,bB->iABc',tmp, QL[ii]) + tmp = contract('iABc,cC->iABC',tmp, QL[jj]) + Zmbij += contract('e,f,mbef->mb',t1[i], t1[j], tmp) - Wmbej_1 = np.zeros_like(Wmbej) - Wmbej_2 = np.zeros_like(Wmbej) - Wmbej_3 = np.zeros_like(Wmbej) - Wmbej_4 = np.zeros_like(Wmbej) - for n in range(self.no): - nn = n*self.no + n - jn = j*self.no + n - nj = n*self.no + j - - Sijjn = QL[ij].T @ QL[jn] + Zmbij_ij.append(Zmbij) - tmp5 = 0.5 * t2_ij[jn] @ Sijjn.T - tmp6_1 = QL[im].T @ self.H.ERI[m,n,v,v] - tmp6_2 = tmp6_1 @ QL[jn] - Wmbej_3 -= tmp5.T @ tmp6_2.T - - Sijnj = QL[ij].T @ QL[nj] + #self.zmbij_t.stop() + return Zmbij_ij - tmp7 = t2_ij[nj] @ Sijnj.T - tmp8_1 = QL[im].T @ self.H.L[m,n,v,v] - tmp8_2 = tmp8_1 @ QL[nj] - Wmbej_4 += 0.5 * tmp7.T @ tmp8_2.T + def build_Wmbej(self, Wmbej_ijim, ERI, L, ERIoovo, Sijnn, Sijnj, Sijjn, t1, t2): + #self.wmbej_t.start() + v = self.v + o = self.o + QL = self.QL + dim = self.dim - Wmbej_ijim.append(Wmbej + Wmbej_3 + Wmbej_4) + if self.model == 'CCD': + for ij in range(self.no*self.no): + i = ij // self.no + j = ij % self.no + for m in range(self.no): + im = i*self.no + m + + Wmbej = np.zeros((dim[ij],dim[im])) + + tmp = QL[ij].T @ ERI[m,v,v,j] + Wmbej = tmp @ QL[im] + + for n in range(self.no): + jn = j*self.no + n + nj = n*self.no + j + ijn = ij*self.no + n + + tmp = 0.5 * t2[jn] @ Sijjn[ijn].T + tmp1 = QL[im].T @ ERI[m,n,v,v] + tmp1 = tmp1 @ QL[jn] + Wmbej -= tmp.T @ tmp1.T + + tmp2 = t2[nj] @ Sijnj[ijn].T + tmp3 = QL[im].T @ L[m,n,v,v] + tmp3 = tmp3 @ QL[nj] + Wmbej += 0.5 * tmp2.T @ tmp3.T + + Wmbej_ijim.append(Wmbej) + else: + for ij in range(self.no*self.no): + i = ij // self.no + j = ij % self.no + jj = j*self.no + j + for m in range(self.no): + im = i*self.no + m + + Wmbej = np.zeros((dim[ij],dim[im])) + + tmp = QL[ij].T @ self.H.ERI[m,v,v,j] + Wmbej = tmp @ QL[im] + + tmp = contract('abc,aA->Abc',ERI[m,v,v,v], QL[ij]) + tmp = contract('Abc,bB->ABc',tmp, QL[im]) + tmp = contract('ABc,cC->ABC',tmp, QL[jj]) + Wmbej += contract('F,beF->be', t1[j], tmp) + + for n in range(self.no): + nn = n*self.no + n + jn = j*self.no + n + nj = n*self.no + j + ijn = ij*(self.no) + n + + tmp1 = Sijnn[ijn] @ t1[n] + Wmbej -= contract('b,e->be', tmp1, ERIoovo[im][m,n,:,j]) + + tmp2 = 0.5 * t2[jn] @ Sijjn[ijn].T + tmp3_0 = QL[im].T @ ERI[m,n,v,v] + tmp3 = tmp3_0 @ QL[jn] + Wmbej -= tmp2.T @ tmp3.T + + tmp4 = tmp3_0 @ QL[jj] + Wmbej -= contract('f,b,ef-> be',t1[j],tmp1,tmp4) + + tmp5 = t2[nj] @ Sijnj[ijn].T + tmp6 = QL[im].T @ L[m,n,v,v] + tmp6 = tmp6 @ QL[nj] + Wmbej += 0.5 * tmp5.T @ tmp6.T + + Wmbej_ijim.append(Wmbej) #self.wmbej_t.stop() return Wmbej_ijim - def build_lWmbje(self, Wmbje_ijim,Wmbie_ijmj,ERIovov_ij, ERIovvv_ij, ERIoovv_ij, ERIooov_ij, t1_ii, t2_ij): + def build_Wmbje(self, Wmbje_ijim, Wmbie_ijmj, ERI, ERIooov, Sijnn, Sijin, Sijjn, t1, t2): #self.wmbje_t.start() o = self.o v = self.v QL = self.QL dim = self.dim - for ij in range(self.no*self.no): - i = ij // self.no - j = ij % self.no - ii = i*self.no + i - jj = j*self.no + j - - for m in range(self.no): - im = i*self.no + m - mj = m*self.no + j - - Wmbje = np.zeros(dim[ij],dim[im]) - Wmbie = np.zeros(dim[ij],dim[mj]) - - tmp_im = QL[ij].T @ self.H.ERI[m,v,j,v] - tmp_mj = QL[ij].T @ self.H.ERI[m,v,i,v] - Wmbje = -1.0 * tmp_im @ QL[im] - Wmbie = -1.0 * tmp_mj @ QL[mj] - - Wmbje_3 = np.zeros_like(Wmbje) - Wmbie_3 = np.zeros_like(Wmbie) - for n in range(self.no): - nn = n*self.no + n - jn = j*self.no + n - _in = i*self.no + n - - Sijjn = QL[ij].T @ QL[jn] - Sijin = QL[ij].T @ QL[_in] - - tmp5 = 0.5* t2_ij[jn] @ Sijjn.T - tmp6_1 = QL[jn].T @ self.H.ERI[m,n,v,v] - tmp6_2 = tmp6_1 @ QL[im] - Wmbje_3 += tmp5.T @ tmp6_2 - - tmp5_mj = 0.5 * t2_ij[_in] @ Sijin.T - tmp6_1mj = QL[_in].T @ self.H.ERI[m,n,v,v] - tmp6_2mj = tmp6_1mj @ QL[mj] - Wmbie_3 += tmp5_mj.T @ tmp6_2mj - - Wmbje_ijim.append(Wmbje + Wmbje_3) - Wmbie_ijmj.append(Wmbie + Wmbie_3) + if self.model == 'CCD': + for ij in range(self.no*self.no): + i = ij // self.no + j = ij % self.no + + for m in range(self.no): + im = i*self.no + m + mj = m*self.no + j + + Wmbje = np.zeros(dim[ij],dim[im]) + Wmbie = np.zeros(dim[ij],dim[mj]) + + tmp = QL[ij].T @ ERI[m,v,j,v] + tmp_mj = QL[ij].T @ ERI[m,v,i,v] + Wmbje = -1.0 * tmp @ QL[im] + Wmbie = -1.0 * tmp_mj @ QL[mj] + + for n in range(self.no): + jn = j*self.no + n + _in = i*self.no + n + ijn = ij*self.no + n + + tmp1 = 0.5* t2[jn] @ Sijjn[ijn].T + tmp2 = QL[jn].T @ ERI[m,n,v,v] + tmp2 = tmp2 @ QL[im] + Wmbje += tmp1.T @ tmp2 + + tmp1_mj = 0.5 * t2[_in] @ Sijin[ijn].T + tmp2_mj = QL[_in].T @ ERI[m,n,v,v] + tmp2_mj = tmp2_mj @ QL[mj] + Wmbie += tmp1_mj.T @ tmp2_mj + + Wmbje_ijim.append(Wmbje) + Wmbie_ijmj.append(Wmbie) + else: + for ij in range(self.no*self.no): + i = ij // self.no + j = ij % self.no + ii = i*self.no + i + jj = j*self.no + j + + for m in range(self.no): + im = i*self.no + m + mj = m*self.no + j + + Wmbje = np.zeros(dim[ij],dim[im]) + Wmbie = np.zeros(dim[ij],dim[mj]) + + tmp = QL[ij].T @ ERI[m,v,j,v] + tmp_mj = QL[ij].T @ ERI[m,v,i,v] + Wmbje = -1.0 * tmp @ QL[im] + Wmbie = -1.0 * tmp_mj @ QL[mj] + + tmp1_0 = contract('abc,aA->Abc',ERI[m,v,v,v], QL[ij]) + tmp1 = contract('Abc,bB->ABc',tmp1_0, QL[jj]) + tmp1 = contract('ABc,cC->ABC',tmp1, QL[im]) + Wmbje -= contract('F,bFe->be', t1[j], tmp1) + + tmp1_mj = contract('Abc,bB->ABc',tmp1_0, QL[ii]) + tmp1_mj = contract('ABc,cC->ABC',tmp1_mj, QL[mj]) + Wmbie -= contract('F,bFe->be', t1[i], tmp1_mj) + + for n in range(self.no): + nn = n*self.no + n + jn = j*self.no + n + _in = i*self.no + n + ijn = ij*self.no + n + + tmp2 = Sijnn[ijn] @ t1[n] + Wmbje += contract('b,e->be',tmp2,ERIooov[im][m,n,j]) + + Wmbie += contract('b,e->be',tmp2,ERIooov[mj][m,n,i]) + + tmp3 = 0.5 * t2[jn] @ Sijjn[ijn].T + tmp4 = QL[jn].T @ ERI[m,n,v,v] + tmp4 = tmp4 @ QL[im] + Wmbje += tmp3.T @ tmp4 + + tmp5 = QL[jj].T @ ERI[m,n,v,v] + tmp5= tmp5 @ QL[im] + Wmbje += contract('f,b,fe->be',t1[j], tmp2, tmp5) + + tmp2_mj = 0.5 * t2[_in] @ Sijin[ijn].T + tmp3_mj = QL[_in].T @ ERI[m,n,v,v] + tmp3_mj = tmp3_mj @ QL[mj] + Wmbie += tmp2_mj.T @ tmp3_mj + + tmp4_mj = QL[ii].T @ ERI[m,n,v,v] + tmp4_mj = tmp4_mj @ QL[mj] + Wmbie += contract('f,b,fe->be',t1[i], tmp2, tmp4_mj) + + Wmbje_ijim.append(Wmbje) + Wmbie_ijmj.append(Wmbie) #self.wmbje_t.stop() return Wmbje_ijim, Wmbie_ijmj - def lr_T1(self, r1_ii, Fov_ij , ERIovvv_ij, Lovvo_ij, Loovo_ij, t1_ii, t2_ij, Fae_ij , Fme_ij, lFmi): - """ - Implemented up to CCSD but debugging at the CCD level to narrow down terms - """ + def r_T1(self, r1_ii, Fov , ERI, L, Loovo, Sijmm, Sijim, Sijmn, t1, t2, Fae, Fmi, Fme): #self.r1_t.start() - for i in range(self.no): - r1_ii.append(np.zeros_like(t1_ii[i])) - #self.r1_t.stop() - return r1_ii - - def lr_T2(self,r2_ij,ERIoovv_ij, ERIvvvv_ij, ERIovvo_ij, ERIovoo_ij, ERIvvvo_ij, ERIovov_ij, t1_ii, - t2_ij, Fae_ij,lFmi,Fme_ij, lWmnij, Zmbij_ij, Wmbej_ijim, Wmbje_ijim, Wmbie_ijmj): - #self.r2_t.start() v = self.v QL = self.QL - dim = self.dim - nr2_ij = [] - for ij in range(self.no*self.no): - i = ij //self.no - j = ij % self.no - ii = i*self.no + i - jj = j*self.no + j + if self.model == 'CCD': + for i in range(self.no): + r1_ii.append(np.zeros_like(t1[i])) - r2 = np.zeros(dim[ij],dim[ij]) - r2_1one = np.zeros_like(r2) - r2_4 = np.zeros_like(r2) + else: + for i in range(self.no): + ii = i*self.no + i - r2 = 0.5 * ERIoovv_ij[ij][i,j] - r2_1one = t2_ij[ij] @ Fae_ij[ij].T - r2_4 = 0.5 * contract('ef,abef->ab',t2_ij[ij],ERIvvvv_ij[ij]) + r1 = np.zeros(self.Local.dim[ii]) - r2_2one = np.zeros_like(r2) - r2_3 = np.zeros_like(r2) - r2_6 = np.zeros_like(r2) - r2_7 = np.zeros_like(r2) - r2_8 = np.zeros_like(r2) - for m in range(self.no): - mm = m *self.no + m - im = i*self.no + m - mj = m*self.no+ j - ijm = ij*self.no + m + r1 = Fov[ii][i].copy() + r1 += contract('e,ae->a', t1[i], Fae[ii]) - Sijmm = QL[ij].T @ QL[mm] - Sijim = QL[ij].T @ QL[im] - Sijmj = QL[ij].T @ QL[mj] + for m in range(self.no): + mm = m*self.no + m + im = i*self.no + m + mi = m*self.no + i + iim = ii*(self.no) + m - tmp = Sijmm @ t1_ii[m] + tmp = Sijmm[iim] @ t1[m] + r1 -= tmp * Fmi[m,i] - tmp2_1 = Sijim @ t2_ij[im] - tmp2_2 = tmp2_1 @ Sijim.T - r2_2one -= tmp2_2 * lFmi[m,j] + tmp1 = Sijim[iim] @ (2*t2[im] - t2[im].swapaxes(0,1)) + r1 += contract('aE,E->a',tmp1, Fme[im][m]) - tmp5 = Sijim @ (t2_ij[im] - t2_ij[im].swapaxes(0,1)) - r2_6 += tmp5 @ Wmbej_ijim[ijm].T + tmp2 = contract('abc,aA->Abc',ERI[m,v,v,v], QL[ii]) + tmp2 = contract('Abc,bB->ABc',tmp2, QL[mi]) + tmp2 = contract('ABc,cC->ABC',tmp2, QL[mi]) + r1 += contract('EF,aEF->a', (2.0*t2[mi] - t2[mi].swapaxes(0,1)), tmp2) - tmp6 = Sijim @ t2_ij[im] - tmp7 = Wmbej_ijim[ijm] + Wmbje_ijim[ijm] - r2_7 += tmp6 @ tmp7.T + for n in range(self.no): + nn = n*self.no + n - tmp8 = Sijmj @ t2_ij[mj] - tmp9 = Wmbie_ijmj[ijm] - r2_8 += tmp8 @ tmp9.T + tmp3 = contract('ab,aA->Ab', L[n,v,v,i],QL[ii]) + tmp3 = contract('Ab,bB->AB', tmp3, QL[nn]) + r1 += contract('F,aF->a', t1[n], tmp3) - for n in range(self.no): - mn = m*self.no + n + for mn in range(self.no*self.no): + m = mn // self.no + n = mn % self.no + imn = i*(self.no**2) + mn + iimn =ii*(self.no**2) + mn + + tmp4 = Sijmn[iimn] @ t2[mn] + r1 -= contract('aE,E->a',tmp4,Loovo[mn][n,m,:,i]) + + r1_ii.append(r1) + #self.r1_t.stop() + return r1_ii - Sijmn = QL[ij].T @ QL[mn] + def r_T2(self,r2_ij, ERI, ERIoovv, ERIvvvv, ERIovoo, Sijmm, Sijim, Sijmj, Sijnn, Sijmn, t1, t2, Fae ,Fmi, Fme, Wmnij, Zmbij, Wmbej, Wmbje, Wmbie): + #self.r2_t.start() + v = self.v + QL = self.QL + dim = self.dim + + nr2 = [] + if self.model == 'CCD': + for ij in range(self.no*self.no): + i = ij //self.no + j = ij % self.no + ii = i*self.no + i + jj = j*self.no + j + + r2 = np.zeros(dim[ij],dim[ij]) + + r2 = 0.5 * ERIoovv[ij][i,j] + r2 += t2[ij] @ Fae[ij].T + r2 += 0.5 * contract('ef,abef->ab',t2[ij],ERIvvvv[ij]) + + for m in range(self.no): + mm = m *self.no + m + im = i*self.no + m + mj = m*self.no+ j + ijm = ij*self.no + m + + tmp = Sijim[ijm] @ t2[im] + tmp = tmp @ Sijim[ijm].T + r2 -= tmp * Fmi[m,j] + + tmp1 = Sijim[ijm] @ (t2[im] - t2[im].swapaxes(0,1)) + r2 += tmp1 @ Wmbej[ijm].T + + tmp2 = Sijim[ijm] @ t2[im] + tmp3 = Wmbej[ijm] + Wmbje[ijm] + r2 += tmp2 @ tmp3.T + + tmp4 = Sijmj[ijm] @ t2[mj] + r2 += tmp4 @ Wmbie[ijm].T + + for n in range(self.no): + mn = m*self.no + n + ijmn = ij*(self.no**2) + mn + + tmp5 = Sijmn[ijmn] @ t2[mn] + tmp5 = tmp5 @ Sijmn[ijmn].T + r2 += 0.5 * tmp5 * Wmnij[m,n,i,j] + + nr2.append(r2) + else: + for ij in range(self.no*self.no): + i = ij //self.no + j = ij % self.no + ii = i*self.no + i + jj = j*self.no + j + + r2 = np.zeros(dim[ij],dim[ij]) + + r2 = 0.5 * ERIoovv[ij][i,j].copy() + r2 += t2[ij] @ Fae[ij].T + r2 += 0.5 * contract('ef,abef->ab', t2[ij],ERIvvvv[ij]) + + tmp = contract('abcd,aA->Abcd',self.H.ERI[v,v,v,v], QL[ij]) + tmp = contract('Abcd,bB->ABcd', tmp, QL[ij]) + tmp = contract('ABcd,cC->ABCd', tmp, QL[ii]) + tmp = contract('ABCd,dD->ABCD', tmp, QL[jj]) + r2 += 0.5 *contract('e,f,abef->ab',t1[i],t1[j], tmp) + + tmp1 = contract('abc,aA->Abc',ERI[v,v,v,j], QL[ij]) + tmp1 = contract('Abc,bB->ABc',tmp1, QL[ij]) + tmp1 = contract('ABc,cC->ABC',tmp1, QL[ii]) + r2 += contract('E,abE->ab', t1[i], tmp1) + + for m in range(self.no): + mm = m *self.no + m + im = i*self.no + m + mj = m*self.no + j + ijm = ij*self.no + m + + tmp2_0 = Sijmm[ijm] @ t1[m] + tmp2 = contract('b,e->be', tmp2_0, Fme[ij][m]) + r2 -= 0.5 * t2[ij] @ tmp2.T + + tmp3 = Sijim[ijm] @ t2[im] + tmp3 = tmp3 @ Sijim[ijm].T + r2 -= tmp3 * Fmi[m,j] + + tmp4 = contract('E,E->',t1[j], Fme[jj][m]) + r2 -= 0.5 * tmp3 * tmp4 + + r2 -= contract('a,b->ab',tmp2_0,Zmbij[ij][m]) + + tmp5 = Sijim[ijm] @ (t2[im] - t2[im].swapaxes(0,1)) + r2 += tmp5 @ Wmbej[ijm].T + + tmp6 = Sijim[ijm] @ t2[im] + tmp7 = Wmbej[ijm] + Wmbje[ijm] + r2 += tmp6 @ tmp7.T + + tmp8 = Sijmj[ijm] @ t2[mj] + r2 += tmp8 @ Wmbie[ijm].T + + tmp9 = QL[ij].T @ ERI[m,v,v,j] + tmp9 = tmp9 @ QL[ii] + tmp10 = contract ('E,a->Ea', t1[i], tmp2_0) + r2 -= tmp10.T @ tmp9.T + + tmp11 = QL[ij].T @ ERI[m,v,j,v] + tmp11 = tmp11 @ QL[ii] + r2 -= tmp11 @ tmp10 + + r2 -= contract('a,b->ab',tmp2_0, ERIovoo[ij][m,:,i,j]) - tmp4_1 = Sijmn @ t2_ij[mn] - tmp4_2 = tmp4_1 @ Sijmn.T - r2_3 += 0.5 * tmp4_2 * lWmnij[m,n,i,j] + for n in range(self.no): + mn = m*self.no + n + nn = n*self.no + n + ijmn = ij*(self.no**2) + mn + ijn = ij*self.no + n + + tmp12 = Sijmn[ijmn] @ t2[mn] + tmp12 = tmp12 @ Sijmn[ijmn].T + tmp13 = Sijnn[ijn] @ t1[n] + r2 += 0.5 * tmp12 * Wmnij[m,n,i,j] + r2 += 0.5 * contract('a,b->ab',tmp2_0, tmp13) * Wmnij[m,n,i,j] - nr2_ij.append(r2 + r2_1one + r2_2one + r2_3 + r2_4 + r2_6 + r2_7 + r2_8) + nr2.append(r2) for i in range(self.no): for j in range(self.no): ij = i*self.no + j ji = j*self.no + i - - r2_ij.append(nr2_ij[ij].copy() + nr2_ij[ji].copy().transpose()) - + + r2_ij.append(nr2[ij].copy() + nr2[ji].copy().transpose()) + #self.r2_t.stop() - return r2_ij + return r2_ij - def lcc_energy(self,Fov_ij,Loovv_ij,t1_ii,t2_ij): + def lcc_energy(self, Fov, Loovv, t1, t2): #self.energy_t.start() + QL = self.QL + v = self.v + ecc_ii = 0 ecc_ij = 0 ecc = 0 - for i in range(self.no): - for j in range(self.no): - ij = i*self.no + j + if self.model == 'CCD': + for i in range(self.no): + for j in range(self.no): + ij = i*self.no + j + + ecc_ij = contract('ab,ab->',t2[ij],Loovv[ij][i,j]) + ecc += ecc_ij + else: + for i in range(self.no): + ii = i*self.no + i - ecc_ij = contract('ab,ab->',t2_ij[ij],Loovv_ij[ij][i,j]) - ecc += ecc_ij + ecc_ii = 2.0 *contract ('a,a->',Fov[ii][i], t1[i]) + ecc += ecc_ii + for j in range(self.no): + ij = i*self.no + j + ii = i*self.no + i + jj = j*self.no + j + + ecc_ij = contract('ab,ab->',t2[ij],Loovv[ij][i,j]) + tmp2 = QL[ii].T @ self.H.L[i,j,v,v] + tmp2 = tmp2 @ QL[jj] + ecc_ij += contract('a,b,ab->',t1[i], t1[j], tmp2) + ecc += ecc_ij #self.energy_t.stop() return ecc + diff --git a/pycc/local.py b/pycc/local.py index af4c5e5..fcd82e5 100644 --- a/pycc/local.py +++ b/pycc/local.py @@ -49,7 +49,7 @@ class Local(object): _build_PAO(): build PAO orbital rotation tensors _build_PNO(): build PNO orbital rotation tensors _build_PNOpp(): build PNO++ orbital rotation tensors - _trans_integral(): transform Fock matrix and ERI from the MO basis to a local basis + trans_integrals(): transform Fock matrix and ERI from the MO basis to a local basis Notes ----- @@ -267,7 +267,7 @@ def _build_PAO(self): evals = np.delete(evals,toss) for c in range(Xt.shape[1]): Xt[:,c] = Xt[:,c] / evals[c]**0.5 - dim.append(Xt.shape[1]) + dim.append(np.dtype('int64').type(Xt.shape[1])) # just below Eq 51, redundant PAO Fock Ft = contract('pq,pr,rs->qs',Rt,self.H.F_ao,Rt) @@ -295,11 +295,7 @@ def _build_PAO(self): self.L = L # transform between PAO and semicanonical PAO spaces self.dim = dim # dimension of PAO space self.eps = eps # semicananonical PAO energies - - #print("Now doing PAO_mp2") - #self._local_MP2_loop() - #self._sim_MP2_loop() - + def _build_PNO(self): """ Perform MP2 loop in non-canonical MO basis, then construct pair density based on t2 amplitudes @@ -346,17 +342,12 @@ def _build_PNO(self): #temporary way to generate make sure the phase factor of Q_ij and L_ij matches with Q_ji and L_ji for i in range(self.no): - for j in range(self.no): - if i < j: - ij = i*self.no + j - ji = j*self.no + i + for j in range(0,i): + ij = i*self.no + j + ji = j*self.no + i - self.Q[ji] = self.Q[ij] - self.L[ji] = self.L[ij] - - #print("Now doing PNO mp2") - #self._local_MP2_loop() - #self._sim_MP2_loop() + self.Q[ji] = self.Q[ij] + self.L[ji] = self.L[ij] def _build_PNOpp(self): """ @@ -400,10 +391,6 @@ def _build_PNOpp(self): self.L = L # transform between local and semicanonical local spaces self.eps = eps # semicananonical local energies self.dim = dim # dimension of local space - - #print("Now doing PNO++ mp2") - #self._local_MP2_loop() - #self._sim_MP2_loop() def _pert_pairdensity(self,t2): ''' @@ -809,14 +796,13 @@ def filter_res(self, r1, r2): return t1, t2 - def _trans_integrals(self, o, v): + def trans_integrals(self, o, v): """ Transforming all the necessary integrals to the semi-canonical PNO basis, stored in a list with length occ*occ - Naming scheme will have the tensor name with _ij such as Fov_ij - + Notes ----- - There are some integral transformation that may not be needed, double check and remove + contraction notations: i,j,a,b are MOs; A,B,C,D virtual semicanonical PNO """ trans_intstart = time.time() @@ -824,31 +810,24 @@ def _trans_integrals(self, o, v): #Initializing the transformation matrices Q = self.Q L = self.L - - #contraction notation i,j,a,b typically MO; A,B,C,D virtual semicanonical PNO; - - #Initializing transformation and integral lists + QL = [] - QLT = [] - - Fov_ij = [] - Fvv_ij = [] - - ERIoovo_ij = [] - ERIooov_ij = [] - ERIovvv_ij = [] - ERIvvvv_ij = [] - ERIoovv_ij = [] - ERIovvo_ij = [] - ERIvvvo_ij = [] - ERIovov_ij = [] - ERIovoo_ij = [] - - Loovv_ij = [] - Lovvv_ij = [] - Looov_ij = [] - Loovo_ij = [] - Lovvo_ij = [] + Fov = [] + Fvv = [] + ERIoovo = [] + ERIooov = [] + ERIovvv = [] + ERIvvvv = [] + ERIoovv = [] + ERIovvo = [] + ERIvvvo = [] + ERIovov = [] + ERIovoo = [] + Loovv = [] + Lovvv = [] + Looov = [] + Loovo = [] + Lovvo = [] for ij in range(self.no*self.no): i = ij // self.no @@ -856,72 +835,126 @@ def _trans_integrals(self, o, v): QL.append(Q[ij] @ L[ij]) - Fov_ij.append(self.H.F[o,v] @ QL[ij]) - - Fvv_ij.append(L[ij].T @ Q[ij].T @ self.H.F[v,v] @ QL[ij]) + Fov.append(self.H.F[o,v] @ QL[ij]) - ERIoovo_ij.append(contract('ijak,aA->ijAk', self.H.ERI[o,o,v,o],QL[ij])) + Fvv.append(L[ij].T @ Q[ij].T @ self.H.F[v,v] @ QL[ij]) - ERIooov_ij.append(contract('ijka,aA->ijkA', self.H.ERI[o,o,o,v],QL[ij])) + ERIoovo.append(contract('ijak,aA->ijAk', self.H.ERI[o,o,v,o],QL[ij])) - tmp = contract('ijab,aA->ijAb',self.H.ERI[o,o,v,v], QL[ij]) - ERIoovv_ij.append(contract('ijAb,bB->ijAB',tmp,QL[ij])) + ERIooov.append(ERIoovo[ij].swapaxes(0,1).swapaxes(2,3)) - tmp1 = contract('iabc,aA->iAbc',self.H.ERI[o,v,v,v], QL[ij]) - tmp2 = contract('iAbc,bB->iABc',tmp1, QL[ij]) - ERIovvv_ij.append(contract('iABc,cC->iABC',tmp2, QL[ij])) + ERIovoo.append(ERIooov[ij].swapaxes(0,2).swapaxes(1,3)) - tmp3 = contract('abcd,aA->Abcd',self.H.ERI[v,v,v,v], QL[ij]) - tmp4 = contract('Abcd,bB->ABcd',tmp3, QL[ij]) - tmp5 = contract('ABcd,cC->ABCd',tmp4, QL[ij]) - ERIvvvv_ij.append(contract('ABCd,dD->ABCD',tmp5, QL[ij])) + tmp = contract('ijab,aA->ijAb',self.H.ERI[o,o,v,v], QL[ij]) + ERIoovv.append(contract('ijAb,bB->ijAB',tmp,QL[ij])) - tmp6 = contract('iabj,aA->iAbj',self.H.ERI[o,v,v,o], QL[ij]) - ERIovvo_ij.append(contract('iAbj,bB->iABj',tmp6,QL[ij])) + ERIovvo.append(ERIoovv[ij].swapaxes(1,3)) - tmp7 = contract('abci,aA->Abci',self.H.ERI[v,v,v,o], QL[ij]) - tmp8 = contract('Abci,bB->ABci',tmp7, QL[ij]) - ERIvvvo_ij.append(contract('ABci,cC->ABCi',tmp8, QL[ij])) + tmp1 = contract('iajb,aA->iAjb',self.H.ERI[o,v,o,v], QL[ij]) + ERIovov.append(contract('iAjb,bB->iAjB',tmp1, QL[ij])) - tmp9 = contract('iajb,aA->iAjb',self.H.ERI[o,v,o,v], QL[ij]) - ERIovov_ij.append(contract('iAjb,bB->iAjB', tmp9, QL[ij])) + tmp2 = contract('iabc,aA->iAbc',self.H.ERI[o,v,v,v], QL[ij]) + tmp2 = contract('iAbc,bB->iABc',tmp2, QL[ij]) + ERIovvv.append(contract('iABc,cC->iABC',tmp2, QL[ij])) - ERIovoo_ij.append(contract('iajk,aA->iAjk', self.H.ERI[o,v,o,o], QL[ij])) + tmp3 = ERIovvv[ij].swapaxes(0,1).swapaxes(2,3) + ERIvvvo.append(tmp3.swapaxes(1,3)) - Loovo_ij.append(contract('ijak,aA->ijAk', self.H.L[o,o,v,o],QL[ij])) + tmp4 = contract('abcd,aA->Abcd',self.H.ERI[v,v,v,v], QL[ij]) + tmp4 = contract('Abcd,bB->ABcd',tmp4, QL[ij]) + tmp4 = contract('ABcd,cC->ABCd',tmp4, QL[ij]) + ERIvvvv.append(contract('ABCd,dD->ABCD',tmp4, QL[ij])) + + Loovo.append(contract('ijak,aA->ijAk', self.H.L[o,o,v,o],QL[ij])) - tmp10 = contract('ijab,aA->ijAb',self.H.L[o,o,v,v], QL[ij]) - Loovv_ij.append(contract('ijAb,bB->ijAB',tmp10,QL[ij])) + Looov.append(Loovo[ij].swapaxes(0,1).swapaxes(2,3)) - tmp11 = contract('iabc,aA->iAbc',self.H.L[o,v,v,v], QL[ij]) - tmp12 = contract('iAbc,bB->iABc',tmp11, QL[ij]) - Lovvv_ij.append(contract('iABc,cC->iABC',tmp12, QL[ij])) + tmp5 = contract('ijab,aA->ijAb',self.H.L[o,o,v,v], QL[ij]) + Loovv.append(contract('ijAb,bB->ijAB',tmp5,QL[ij])) - Looov_ij.append(contract('ijka,aA->ijkA',self.H.L[o,o,o,v], QL[ij])) + Lovvo.append(Loovv[ij].swapaxes(1,3)) - tmp13 = contract('iabj,aA->iAbj',self.H.L[o,v,v,o], QL[ij]) - Lovvo_ij.append(contract('iAbj,bB->iABj',tmp13,QL[ij])) + tmp6 = contract('iabc,aA->iAbc',self.H.L[o,v,v,v], QL[ij]) + tmp6 = contract('iAbc,bB->iABc',tmp6, QL[ij]) + Lovvv.append(contract('iABc,cC->iABC',tmp6, QL[ij])) - #Storing the list to this class self.QL = QL - self.QLT = QLT - self.Fov_ij = Fov_ij - self.Fvv_ij = Fvv_ij - - self.ERIoovo_ij = ERIoovo_ij - self.ERIooov_ij = ERIooov_ij - self.ERIovvv_ij = ERIovvv_ij - self.ERIvvvv_ij = ERIvvvv_ij - self.ERIoovv_ij = ERIoovv_ij - self.ERIovvo_ij = ERIovvo_ij - self.ERIvvvo_ij = ERIvvvo_ij - self.ERIovov_ij = ERIovov_ij - self.ERIovoo_ij = ERIovoo_ij - - self.Loovv_ij = Loovv_ij - self.Lovvv_ij = Lovvv_ij - self.Looov_ij = Looov_ij - self.Loovo_ij = Loovo_ij - self.Lovvo_ij = Lovvo_ij + self.Fov = Fov + self.Fvv = Fvv + self.ERIoovo = ERIoovo + self.ERIooov = ERIooov + self.ERIovvv = ERIovvv + self.ERIvvvv = ERIvvvv + self.ERIoovv = ERIoovv + self.ERIovvo = ERIovvo + self.ERIvvvo = ERIvvvo + self.ERIovov = ERIovov + self.ERIovoo = ERIovoo + self.Loovv = Loovv + self.Lovvv = Lovvv + self.Looov = Looov + self.Loovo = Loovo + self.Lovvo = Lovvo print("Integrals transformed in %.3f seconds." % (time.time() - trans_intstart)) + + def overlaps(self, QL): + """ + Generating and storing overlap terms + + Notes + ----- + Length: two unique index = no*no, three unique index = no*no*no, four unique index = no*no*no*no + Computational scaling is length*(nv*nv) + Storage size is length*(pno*pno) where pno dimension varies based on cutoff + + To do + ----- + Compare the timings for the use of stored overlap terms versus "on the fly" overlap terms + """ + no = self.no + + Sijmm = [] + Sijim = [] + Sijmj = [] + Sijnn = [] + Sijin = [] + Sijnj = [] + Sijjn = [] + Sijmn = [] + + for i in range(no): + for j in range(no): + ij = i*no + j + for m in range(no): + mm = m*no + m + im = i*no + m + mj = m*no + j + + Sijmm.append(QL[ij].T @ QL[mm]) + Sijim.append(QL[ij].T @ QL[im]) + Sijmj.append(QL[ij].T @ QL[mj]) + + for n in range(no): + nn = n*no + n + _in = i*no + n + nj = n*no + j + jn = j*no + n + ijn = ij*no + n + + Sijnn.append(QL[ij].T @ QL[nn]) + Sijin.append(QL[ij].T @ QL[_in]) + Sijnj.append(QL[ij].T @ QL[nj]) + Sijjn.append(QL[ij].T @ QL[jn]) + + for mn in range(no*no): + Sijmn.append(QL[ij].T @ QL[mn]) + + self.Sijmj = Sijmj + self.Sijmm = Sijmm + self.Sijim = Sijim + self.Sijnn = Sijnn + self.Sijin = Sijin + self.Sijnj = Sijnj + self.Sijjn = Sijjn + self.Sijmn = Sijmn diff --git a/pycc/utils.py b/pycc/utils.py index 1a30b1c..f99f36f 100644 --- a/pycc/utils.py +++ b/pycc/utils.py @@ -112,6 +112,9 @@ def extrapolate(self, t1, t2): if n1 >= n2: continue B[n1, n2] = np.dot(e1, e2) + print("B",B) + print("e1", n1,e1) + print("e2",n2,e2) B[n2, n1] = B[n1, n2] B[:-1, :-1] /= np.abs(B[:-1, :-1]).max() @@ -194,8 +197,11 @@ def add_error_vector(self, no,t1_ii, t2_ij): error_t1.append((self.diis_vals_t1[-1][i] - self.oldt1[i]).ravel()) for j in range(no): ij = i*no + j - + error_t2.append((self.diis_vals_t2[-1][ij] - self.oldt2[ij]).ravel()) + if ij ==2: + print(error_t2[ij]) + print("len of error_t2", len(error_t2[ij]), self.oldt2[ij].shape) #if ij == 24: #print("diis t2 last", ij, self.diis_vals_t2[-1][ij]) #print("oldt2 last", ij, self.oldt2[ij]) @@ -236,7 +242,7 @@ def extrapolate(self, no,t1_ii, t2_ij): B_t1[n1, n2] = np.dot(e1[i], e2[i]) B_t1[n2, n1] = B_t1[n1, n2] B_t1[:-1, :-1] /= np.abs(B_t1[:-1, :-1]).max() - + print("B_t1", B_t1) # Build residual vector resid = np.zeros(self.diis_size_t1 + 1) resid[-1] = -1 @@ -258,10 +264,13 @@ def extrapolate(self, no,t1_ii, t2_ij): if n1 >= n2: continue B_t2[n1, n2] = np.dot(e1[ij], e2[ij]) + print("B_t2",B_t2) + print("e1", n1, e1[ij]) + print("e2", n2, e2[ij]) B_t2[n2, n1] = B_t2[n1, n2] - + print("B_t2 before", B_t2) B_t2[:-1, :-1] /= np.abs(B_t2[:-1, :-1]).max() - + print("B_t2", B_t2) # Build residual vector resid = np.zeros(self.diis_size_t2 + 1) resid[-1] = -1 diff --git a/pycc/utils_copy.py b/pycc/utils_copy.py index e904ca2..663a093 100644 --- a/pycc/utils_copy.py +++ b/pycc/utils_copy.py @@ -146,16 +146,17 @@ def extrapolate(self, t1, t2): return t1, t2 class helper_ldiis(object): - def __init__(self, no,t1_ii, t2_ij, max_diis): + def __init__(se +lf, no,t1_ii, t2_ij, max_diis): oldt1 = [] - oldt2 = [] + oldt2 = []w for i in range(no): - oldt1.append(np.array(t1_ii[i],dtype="object")) - for j in range(no): + oldt1.append(t1_ii[i].copy(),dtype="object")) + for j in range(no):: ij = i*no + j - oldt2.append(np.array(t2_ij[ij],dtype="object")) + oldt2.append(t2_ij[ij].copy()) self.oldt1 = oldt1 self.oldt2 = oldt2 diff --git a/pycc/utils_copy1.py b/pycc/utils_copy1.py new file mode 100644 index 0000000..be09457 --- /dev/null +++ b/pycc/utils_copy1.py @@ -0,0 +1,331 @@ +import numpy as np +import torch +import opt_einsum + + +class helper_diis(object): + def __init__(self, t1, t2, max_diis, precision='DP'): + if isinstance(t1, torch.Tensor): + self.device0 = torch.device('cpu') + self.device1 = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + self.diis_vals_t1 = [t1.clone()] + self.diis_vals_t2 = [t2.clone()] + else: + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + self.diis_vals_t1 = [t1.copy()] + self.diis_vals_t2 = [t2.copy()] + + self.diis_errors = [] + self.diis_size = 0 + self.max_diis = max_diis + self.precision = precision + + def add_error_vector(self, t1, t2): + if isinstance(t1, torch.Tensor): + # Add DIIS vectors + self.diis_vals_t1.append(t1.clone()) + self.diis_vals_t2.append(t2.clone()) + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(torch.cat((error_t1, error_t2))) + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + else: + # Add DIIS vectors + self.diis_vals_t1.append(t1.copy()) + self.diis_vals_t2.append(t2.copy()) + + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(np.concatenate((error_t1, error_t2))) + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + + def extrapolate(self, t1, t2): + + if (self.max_diis == 0): + return t1, t2 + + # Limit size of DIIS vector + if (len(self.diis_errors) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors[0] + + self.diis_size = len(self.diis_errors) + + if isinstance(t1, torch.Tensor): + # Build error matrix B + if self.precision == 'DP': + B = torch.ones((self.diis_size + 1, self.diis_size + 1), dtype=torch.float64, device=self.device1) * -1 + elif self.precision == 'SP': + B = torch.ones((self.diis_size + 1, self.diis_size + 1), dtype=torch.float32, device=self.device1) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = torch.dot(e1, e1) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = torch.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= torch.abs(B[:-1, :-1]).max() + + # Build residual vector + if self.precision == 'DP': + resid = torch.zeros((self.diis_size + 1), dtype=torch.float64, device=self.device1) + elif self.precision == 'SP': + resid = torch.zeros((self.diis_size + 1), dtype=torch.float32, device=self.device1) + resid[-1] = -1 + + # Solve pulay equations + ci = torch.linalg.solve(B, resid) + + # Calculate new amplitudes + t1 = torch.zeros_like(self.oldt1) + t2 = torch.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += torch.real(ci[num] * self.diis_vals_t1[num + 1]) + t2 += torch.real(ci[num] * self.diis_vals_t2[num + 1]) + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + + else: + # Build error matrix B + if self.precision == 'DP': + B = np.ones((self.diis_size + 1, self.diis_size + 1)) * -1 + elif self.precision == 'SP': + B = np.ones((self.diis_size + 1, self.diis_size + 1), dtype=np.float32) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = np.dot(e1, e1) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = np.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= np.abs(B[:-1, :-1]).max() + + # Build residual vector + if self.precision == 'DP': + resid = np.zeros(self.diis_size + 1) + elif self.precision == 'SP': + resid = np.zeros((self.diis_size + 1), dtype=np.float32) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B, resid) + + # Calculate new amplitudes + t1 = np.zeros_like(self.oldt1) + t2 = np.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += ci[num] * self.diis_vals_t1[num + 1] + t2 += ci[num] * self.diis_vals_t2[num + 1] + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + + return t1, t2 + +class helper_ldiis(object): + def __init__(self, no,t1_ii, t2_ij, max_diis): + oldt1 = [] + oldt2 = [] + diis_vals_t1 = [] + diis_vals_t2 = [] + for i in range(no): + oldt1.append(t1_ii[i].copy()) + for j in range(no): + ij = i*no + j + oldt2.append(t2_ij[ij].copy()) + + self.oldt1 = oldt1 + self.oldt2 = oldt2 + diis_vals_t1.append(self.oldt1.copy()) + diis_vals_t2.append(self.oldt2.copy()) + self.diis_vals_t1 = diis_vals_t1 + self.diis_vals_t2 = diis_vals_t2 + #print("diis t2 initial", self.diis_vals_t2) + self.diis_errors_t1 = [] + self.diis_errors_t2 = [] + self.diis_size = 0 + self.max_diis = max_diis + + def add_error_vector(self, no,t1_ii, t2_ij): + # Add DIIS vectors + self.addt1 = [] + self.addt2 = [] + + for i in range(no): + ii = i*no + i + + self.addt1.append(t1_ii[i]) + + for j in range(no): + ij = i*no + j + + self.addt2.append(t2_ij[ij]) + #if ij == 24: + #print("t2 compare", ij, t2_ij[ij]) + + self.diis_vals_t1.append(self.addt1.copy(),dtype="object")) + self.diis_vals_t2.append(np.array(self.addt2.copy(),dtype="object")) + #for n1, e1 in enumerate(self.diis_vals_t2): + #print("n1", n1, "e1", e1[24]) + + #print("oldt2", self.oldt2[ij]) + + # Add new error vectors + error_t1 = [] + error_t2 = [] + for i in range(no): + error_t1.append((self.diis_vals_t1[-1][i] - self.oldt1[i]).ravel()) + for j in range(no): + ij = i*no + j + + error_t2.append((self.diis_vals_t2[-1][ij] - self.oldt2[ij]).ravel()) + #if ij == 24: + #print("diis t2 last", ij, self.diis_vals_t2[-1][ij]) + #print("oldt2 last", ij, self.oldt2[ij]) + #print(error_t2[24]) + self.diis_errors_t1.append(error_t1) + self.diis_errors_t2.append(error_t2) + self.oldt1 = self.addt1.copy() + self.oldt2 = self.addt2.copy() + + def extrapolate(self, no,t1_ii, t2_ij): + + if (self.max_diis == 0): + return t1_ii, t2_ij + + # Limit size of DIIS vector + if (len(self.diis_errors_t1) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors_t1[0] + del self.diis_errors_t2[0] + + self.diis_size_t1 = len(self.diis_errors_t1) + self.diis_size_t2 = len(self.diis_errors_t2) + B_t1 = np.ones((self.diis_size_t1 + 1, self.diis_size_t1 + 1)) * -1 + B_t1[-1, -1] = 0 + B_t2 = np.ones((self.diis_size_t2 + 1, self.diis_size_t2 + 1)) * -1 + B_t2[-1, -1] = 0 + + t1_ii = np.zeros_like(self.oldt1) + t2_ij = np.zeros_like(self.oldt2) + for i in range(no): + + for n1, e1 in enumerate(self.diis_errors_t1): + B_t1[n1, n1] = np.dot(e1[i], e1[i]) + for n2, e2 in enumerate(self.diis_errors_t1): + if n1 >= n2: + continue + B_t1[n1, n2] = np.dot(e1[i], e2[i]) + B_t1[n2, n1] = B_t1[n1, n2] + B_t1[:-1, :-1] /= np.abs(B_t1[:-1, :-1]).max() + + # Build residual vector + resid = np.zeros(self.diis_size_t1 + 1) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B_t1, resid) + + # Calculate new amplitudes + for num in range(self.diis_size_t1): + t1_ii[i] = t1_ii[i] + np.array(ci[num] * self.diis_vals_t1[num + 1][i]) + + # Save extrapolated amplitudes to old_t amplitudes + for j in range(no): + ij = i*no + j + + for n1, e1 in enumerate(self.diis_errors_t2): + B_t2[n1, n1] = np.dot(e1[ij], e1[ij]) + for n2, e2 in enumerate(self.diis_errors_t2): + if n1 >= n2: + continue + B_t2[n1, n2] = np.dot(e1[ij], e2[ij]) + B_t2[n2, n1] = B_t2[n1, n2] + + B_t2[:-1, :-1] /= np.abs(B_t2[:-1, :-1]).max() + + # Build residual vector + resid = np.zeros(self.diis_size_t2 + 1) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B_t2, resid) + #print("ci", ci, "ci.shape", ci.shape) + + # Calculate new amplitudes + for num in range(self.diis_size_t2): + t2_ij[ij] = t2_ij[ij] + np.array(ci[num] * self.diis_vals_t2[num + 1][ij]) + #print("old", self.oldt2[24]) + #print("new", t2_ij[24]) + print(type(t1_ii.copy())) + self.oldt1 = t1_ii.copy().tolist() + self.oldt2 = t2_ij.copy().tolist() + print("this",type(t1_ii)) + return t1_ii, t2_ij + +class cc_contract(object): + """ + A wrapper for opt_einsum.contract with tensors stored on CPU and/or GPU. + """ + def __init__(self, device='CPU'): + """ + Parameters + ---------- + device: string + initiated in ccwfn object, default: 'CPU' + + Returns + ------- + None + """ + self.device = device + if self.device == 'GPU': + # torch.device is an object representing the device on which torch.Tensor is or will be allocated. + self.device0 = torch.device('cpu') + self.device1 = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') + + def __call__(self, subscripts, *operands): + """ + Parameters + ---------- + subscripts: string + specify the subscripts for summation (same format as numpy.einsum) + *operands: list of array_like + the arrays/tensors for the operation + + Returns + ------- + An ndarray/torch.Tensor that is calculated based on Einstein summation convention. + """ + if self.device == 'CPU': + return opt_einsum.contract(subscripts, *operands) + elif self.device == 'GPU': + # Check the type and allocation of the tensors + # Transfer the copy from CPU to GPU if needed (for ERI) + input_list = list(operands) + for i in range(len(input_list)): + if (not input_list[i].is_cuda): + input_list[i] = input_list[i].to(self.device1) + #print(len(input_list), type(input_list[0]), type(input_list[1])) + output = opt_einsum.contract(subscripts, *input_list) + del input_list + return output + From 47228ad1bea2ba72e407c5f1ab00fb7fa36a48b9 Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Fri, 17 Mar 2023 13:12:01 -0400 Subject: [PATCH 3/8] cleaned up a copy of utils.py --- pycc/utils_copy.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pycc/utils_copy.py b/pycc/utils_copy.py index 663a093..4cc1917 100644 --- a/pycc/utils_copy.py +++ b/pycc/utils_copy.py @@ -146,8 +146,7 @@ def extrapolate(self, t1, t2): return t1, t2 class helper_ldiis(object): - def __init__(se -lf, no,t1_ii, t2_ij, max_diis): + def __init__(self, no,t1_ii, t2_ij, max_diis): oldt1 = [] oldt2 = []w for i in range(no): From 2e443de1b78f3e734ad7fa72ca2f74e9115de71d Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Mon, 20 Mar 2023 09:41:56 -0400 Subject: [PATCH 4/8] working code in progress --- pycc/ccwfn.py | 2 +- pycc/lccwfn.py | 2 +- pycc/test_032_pnoccd.py | 39 +++++ pycc/utils_correct.py | 349 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 390 insertions(+), 2 deletions(-) create mode 100644 pycc/test_032_pnoccd.py create mode 100644 pycc/utils_correct.py diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py index 66dda29..c583830 100644 --- a/pycc/ccwfn.py +++ b/pycc/ccwfn.py @@ -10,7 +10,7 @@ import time import numpy as np import torch -from utils import helper_diis, cc_contract +from utils_correct import helper_diis, cc_contract from hamiltonian import Hamiltonian from local import Local from cctriples import t_tjl, t3c_ijk diff --git a/pycc/lccwfn.py b/pycc/lccwfn.py index dcd2ed8..7a64d10 100644 --- a/pycc/lccwfn.py +++ b/pycc/lccwfn.py @@ -2,7 +2,7 @@ #from timer import Timer import numpy as np from opt_einsum import contract -from utils import helper_ldiis +from utils_correct import helper_ldiis class lccwfn(object): """ diff --git a/pycc/test_032_pnoccd.py b/pycc/test_032_pnoccd.py new file mode 100644 index 0000000..3141f6b --- /dev/null +++ b/pycc/test_032_pnoccd.py @@ -0,0 +1,39 @@ +""" +Test basic PNO-CCD energy +""" +# Import package, test suite, and other packages as needed +import psi4 +from ccwfn import ccwfn +from lccwfn import lccwfn +from data.molecules import * + +psi4.set_memory('2 GB') +psi4.core.set_output_file('output.dat', False) +psi4.set_options({'basis': 'cc-pvdz', + 'scf_type': 'pk', + 'mp2_type': 'conv', + 'freeze_core': 'false', + 'e_convergence': 1e-13, + 'd_convergence': 1e-13, + 'r_convergence': 1e-13, + 'diis': 1}) +mol = psi4.geometry(moldict["H2O"]) +rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + +maxiter = 100 +e_conv = 1e-5 +r_conv = 1e-5 + +#simulation code of pno-ccd +ccd_sim = ccwfn(rhf_wfn, model='CCSD',local='PNO', local_cutoff=1e-7,it2_opt=False,filter=True) +eccd_sim = ccd_sim.solve_cc(e_conv, r_conv, maxiter) + +#pno-ccd +lccd = ccwfn(rhf_wfn,model='CCSD', local='PNO', local_cutoff=1e-7,it2_opt=False) +elccd = lccd.lccwfn.solve_lcc(e_conv, r_conv, maxiter) + +print("eccd_sim",eccd_sim) +print("elccd", elccd) +assert(abs(eccd_sim - elccd) < 1e-5) + + diff --git a/pycc/utils_correct.py b/pycc/utils_correct.py new file mode 100644 index 0000000..faa2f3a --- /dev/null +++ b/pycc/utils_correct.py @@ -0,0 +1,349 @@ +import numpy as np +import torch +import opt_einsum + + +class helper_diis(object): + def __init__(self, t1, t2, max_diis, precision='DP'): + if isinstance(t1, torch.Tensor): + self.device0 = torch.device('cpu') + self.device1 = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + self.diis_vals_t1 = [t1.clone()] + self.diis_vals_t2 = [t2.clone()] + else: + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + self.diis_vals_t1 = [t1.copy()] + self.diis_vals_t2 = [t2.copy()] + + self.diis_errors = [] + self.diis_size = 0 + self.max_diis = max_diis + self.precision = precision + + def add_error_vector(self, t1, t2): + if isinstance(t1, torch.Tensor): + # Add DIIS vectors + self.diis_vals_t1.append(t1.clone()) + self.diis_vals_t2.append(t2.clone()) + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(torch.cat((error_t1, error_t2))) + print("length of error", len(self.diis_errors)) + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + else: + # Add DIIS vectors + self.diis_vals_t1.append(t1.copy()) + self.diis_vals_t2.append(t2.copy()) + + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(np.concatenate((error_t1, error_t2))) + #print("length of error", len(self.diis_errors)) + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + + def extrapolate(self, t1, t2): + + if (self.max_diis == 0): + return t1, t2 + + # Limit size of DIIS vector + if (len(self.diis_errors) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors[0] + + self.diis_size = len(self.diis_errors) + + if isinstance(t1, torch.Tensor): + # Build error matrix B + if self.precision == 'DP': + B = torch.ones((self.diis_size + 1, self.diis_size + 1), dtype=torch.float64, device=self.device1) * -1 + elif self.precision == 'SP': + B = torch.ones((self.diis_size + 1, self.diis_size + 1), dtype=torch.float32, device=self.device1) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = torch.dot(e1, e1) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = torch.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= torch.abs(B[:-1, :-1]).max() + + # Build residual vector + if self.precision == 'DP': + resid = torch.zeros((self.diis_size + 1), dtype=torch.float64, device=self.device1) + elif self.precision == 'SP': + resid = torch.zeros((self.diis_size + 1), dtype=torch.float32, device=self.device1) + resid[-1] = -1 + + # Solve pulay equations + ci = torch.linalg.solve(B, resid) + + # Calculate new amplitudes + t1 = torch.zeros_like(self.oldt1) + t2 = torch.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += torch.real(ci[num] * self.diis_vals_t1[num + 1]) + t2 += torch.real(ci[num] * self.diis_vals_t2[num + 1]) + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + + else: + # Build error matrix B + if self.precision == 'DP': + B = np.ones((self.diis_size + 1, self.diis_size + 1)) * -1 + elif self.precision == 'SP': + B = np.ones((self.diis_size + 1, self.diis_size + 1), dtype=np.float32) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = np.dot(e1, e1) + #print("B",B) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = np.dot(e1, e2) + #print("B",B) + #print("e1", n1,e1) + #print("e2",n2,e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= np.abs(B[:-1, :-1]).max() + + # Build residual vector + if self.precision == 'DP': + resid = np.zeros(self.diis_size + 1) + elif self.precision == 'SP': + resid = np.zeros((self.diis_size + 1), dtype=np.float32) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B, resid) + + # Calculate new amplitudes + t1 = np.zeros_like(self.oldt1) + t2 = np.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += ci[num] * self.diis_vals_t1[num + 1] + t2 += ci[num] * self.diis_vals_t2[num + 1] + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + + return t1, t2 + +class helper_ldiis(object): + def __init__(self, no,t1_ii, t2_ij, max_diis): + oldt1 = [] + oldt2 = [] + diis_vals_t1 = [] + diis_vals_t2 = [] + for i in range(no): + oldt1.append(t1_ii[i].copy()) + for j in range(no): + ij = i*no + j + oldt2.append(t2_ij[ij].copy()) + + self.oldt1 = oldt1.copy() + self.oldt2 = oldt2.copy() + diis_vals_t1.append(self.oldt1.copy()) + diis_vals_t2.append(self.oldt2.copy()) + self.diis_vals_t1 = diis_vals_t1 + self.diis_vals_t2 = diis_vals_t2 + #print("diis t2 initial", self.diis_vals_t2) + self.diis_errors_t1 = [] + self.diis_errors_t2 = [] + self.diis_size = 0 + self.max_diis = max_diis + + def add_error_vector(self, no,t1_ii, t2_ij): + # Add DIIS vectors + self.addt1 = [] + self.addt2 = [] + + for i in range(no): + ii = i*no + i + + self.addt1.append(t1_ii[i].copy()) + + for j in range(no): + ij = i*no + j + + self.addt2.append(t2_ij[ij].copy()) + #if ij == 2: + #print("new t2 part 1", ij, t2_ij[ij]) + + self.diis_vals_t1.append(self.addt1.copy()) + self.diis_vals_t2.append(self.addt2.copy()) + #print("new t2", self.diis_vals_t2[-1][2]) + #for n1, e1 in enumerate(self.diis_vals_t2): + #print("n1", n1, "e1", e1[24]) + + #print("oldt2", self.oldt2[2]) + + # Add new error vectors + error_t1 = [] + error_t2 = [] + for i in range(no): + #error_t1.append((self.diis_vals_t1[-1][i] - self.oldt1[i]).ravel()) + error_t1.append((self.diis_vals_t1[-1][i] - self.diis_vals_t1[len(self.diis_vals_t1) -2][i]).ravel()) + for j in range(no): + ij = i*no + j + + #error_t2.append((self.diis_vals_t2[-1][ij] - self.oldt2[ij]).ravel()) + error_t2.append((self.diis_vals_t2[-1][ij] - self.diis_vals_t2[len(self.diis_vals_t2) -2][ij]).ravel()) + #if ij ==2: + #print("new t2", self.diis_vals_t2[-1][ij]) + #print("oldt2",self.oldt2[ij]) + #print("this should be the old t2", self.diis_vals_t2[len(self.diis_vals_t2) -2][ij]) + #print("diff", error_t2[ij]) + #print("len of error_t2", len(error_t2[ij]), self.oldt2[ij].shape) + #if ij == 24: + #print("diis t2 last", ij, self.diis_vals_t2[-1][ij]) + #print("oldt2 last", ij, self.oldt2[ij]) + #print(error_t2[24]) + self.diis_errors_t1.append(error_t1) + #print("length of error", len(self.diis_errors_t1)) + self.diis_errors_t2.append(error_t2) + self.oldt1 = self.diis_vals_t1[-1].copy() + self.oldt2 = self.diis_vals_t2[-1].copy() + #print("this should be old t2", self.oldt2[2]) + + def extrapolate(self, no,t1_ii, t2_ij): + + if (self.max_diis == 0): + return t1_ii, t2_ij + + # Limit size of DIIS vector + if (len(self.diis_errors_t1) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors_t1[0] + del self.diis_errors_t2[0] + + self.diis_size_t1 = len(self.diis_errors_t1) + self.diis_size_t2 = len(self.diis_errors_t2) + B_t1 = np.ones((self.diis_size_t1 + 1, self.diis_size_t1 + 1)) * -1 + B_t1[-1, -1] = 0 + B_t2 = np.ones((self.diis_size_t2 + 1, self.diis_size_t2 + 1)) * -1 + B_t2[-1, -1] = 0 + + t1_ii = np.zeros_like(self.oldt1) + t2_ij = np.zeros_like(self.oldt2) + for i in range(no): + + for n1, e1 in enumerate(self.diis_errors_t1): + B_t1[n1, n1] = np.dot(e1[i], e1[i]) + for n2, e2 in enumerate(self.diis_errors_t1): + if n1 >= n2: + continue + B_t1[n1, n2] = np.dot(e1[i], e2[i]) + B_t1[n2, n1] = B_t1[n1, n2] + B_t1[:-1, :-1] /= np.abs(B_t1[:-1, :-1]).max() + #print("B_t1", B_t1) + # Build residual vector + resid = np.zeros(self.diis_size_t1 + 1) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B_t1, resid) + + # Calculate new amplitudes + for num in range(self.diis_size_t1): + t1_ii[i] += ci[num] * self.diis_vals_t1[num + 1][i] + + # Save extrapolated amplitudes to old_t amplitudes + for j in range(no): + ij = i*no + j + + for n1, e1 in enumerate(self.diis_errors_t2): + B_t2[n1, n1] = np.dot(e1[ij], e1[ij]) + for n2, e2 in enumerate(self.diis_errors_t2): + if n1 >= n2: + continue + B_t2[n1, n2] = np.dot(e1[ij], e2[ij]) + #print("B_t2",B_t2) + #print("e1", n1, e1[ij]) + #print("e2", n2, e2[ij]) + B_t2[n2, n1] = B_t2[n1, n2] + #print("B_t2 before", B_t2) + B_t2[:-1, :-1] /= np.abs(B_t2[:-1, :-1]).max() + #print("B_t2", B_t2) + # Build residual vector + resid = np.zeros(self.diis_size_t2 + 1) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B_t2, resid) + #print("ci", ci, "ci.shape", ci.shape) + + # Calculate new amplitudes + for num in range(self.diis_size_t2): + t2_ij[ij] += ci[num] * self.diis_vals_t2[num + 1][ij] + #print("old", self.oldt2[24]) + #print("new", t2_ij[24]) + self.oldt1 = t1_ii.copy() + self.oldt2 = t2_ij.copy() + + return t1_ii, t2_ij + +class cc_contract(object): + """ + A wrapper for opt_einsum.contract with tensors stored on CPU and/or GPU. + """ + def __init__(self, device='CPU'): + """ + Parameters + ---------- + device: string + initiated in ccwfn object, default: 'CPU' + + Returns + ------- + None + """ + self.device = device + if self.device == 'GPU': + # torch.device is an object representing the device on which torch.Tensor is or will be allocated. + self.device0 = torch.device('cpu') + self.device1 = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') + + def __call__(self, subscripts, *operands): + """ + Parameters + ---------- + subscripts: string + specify the subscripts for summation (same format as numpy.einsum) + *operands: list of array_like + the arrays/tensors for the operation + + Returns + ------- + An ndarray/torch.Tensor that is calculated based on Einstein summation convention. + """ + if self.device == 'CPU': + return opt_einsum.contract(subscripts, *operands) + elif self.device == 'GPU': + # Check the type and allocation of the tensors + # Transfer the copy from CPU to GPU if needed (for ERI) + input_list = list(operands) + for i in range(len(input_list)): + if (not input_list[i].is_cuda): + input_list[i] = input_list[i].to(self.device1) + #print(len(input_list), type(input_list[0]), type(input_list[1])) + output = opt_einsum.contract(subscripts, *input_list) + del input_list + return output From 6fc0bf34de02f9d141dce574c2b0f598f519a052 Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Tue, 21 Mar 2023 11:40:00 -0400 Subject: [PATCH 5/8] added two ways (in progress) to do DIIS --- pycc/lccwfn.py | 4 +- pycc/test_032_pnoccd.py | 6 +- pycc/utils_test.py | 326 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 331 insertions(+), 5 deletions(-) create mode 100644 pycc/utils_test.py diff --git a/pycc/lccwfn.py b/pycc/lccwfn.py index 7a64d10..e66cc0d 100644 --- a/pycc/lccwfn.py +++ b/pycc/lccwfn.py @@ -2,7 +2,7 @@ #from timer import Timer import numpy as np from opt_einsum import contract -from utils_correct import helper_ldiis +from utils_test import helper_ldiis class lccwfn(object): """ @@ -109,7 +109,7 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis #self.r2_t = Timer("r2") #self.energy_t = Timer("energy") - ldiis = helper_ldiis(self.no, self.t1, self.t2, max_diis) + ldiis = helper_ldiis(self.no, self.t1, self.t2, max_diis, self.Local) elcc = self.lcc_energy(self.Local.Fov,self.Local.Loovv,self.t1, self.t2) print("CC Iter %3d: lCC Ecorr = %.15f dE = % .5E MP2" % (0,elcc,-elcc)) diff --git a/pycc/test_032_pnoccd.py b/pycc/test_032_pnoccd.py index 3141f6b..2fcf5b6 100644 --- a/pycc/test_032_pnoccd.py +++ b/pycc/test_032_pnoccd.py @@ -9,7 +9,7 @@ psi4.set_memory('2 GB') psi4.core.set_output_file('output.dat', False) -psi4.set_options({'basis': 'cc-pvdz', +psi4.set_options({'basis': 'aug-cc-pvdz', 'scf_type': 'pk', 'mp2_type': 'conv', 'freeze_core': 'false', @@ -21,8 +21,8 @@ rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) maxiter = 100 -e_conv = 1e-5 -r_conv = 1e-5 +e_conv = 1e-12 +r_conv = 1e-12 #simulation code of pno-ccd ccd_sim = ccwfn(rhf_wfn, model='CCSD',local='PNO', local_cutoff=1e-7,it2_opt=False,filter=True) diff --git a/pycc/utils_test.py b/pycc/utils_test.py new file mode 100644 index 0000000..30c778a --- /dev/null +++ b/pycc/utils_test.py @@ -0,0 +1,326 @@ +import numpy as np +import torch +import opt_einsum + + +class helper_diis(object): + def __init__(self, t1, t2, max_diis, precision='DP'): + if isinstance(t1, torch.Tensor): + self.device0 = torch.device('cpu') + self.device1 = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + self.diis_vals_t1 = [t1.clone()] + self.diis_vals_t2 = [t2.clone()] + else: + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + print(self.oldt1.shape) + print(self.oldt1) + print(self.oldt2.shape) + self.diis_vals_t1 = [t1.copy()] + self.diis_vals_t2 = [t2.copy()] + + self.diis_errors = [] + self.diis_size = 0 + self.max_diis = max_diis + self.precision = precision + + def add_error_vector(self, t1, t2): + if isinstance(t1, torch.Tensor): + # Add DIIS vectors + self.diis_vals_t1.append(t1.clone()) + self.diis_vals_t2.append(t2.clone()) + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(torch.cat((error_t1, error_t2))) + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + else: + # Add DIIS vectors + self.diis_vals_t1.append(t1.copy()) + self.diis_vals_t2.append(t2.copy()) + + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(np.concatenate((error_t1, error_t2))) + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + + def extrapolate(self, t1, t2): + + if (self.max_diis == 0): + return t1, t2 + + # Limit size of DIIS vector + if (len(self.diis_errors) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors[0] + + self.diis_size = len(self.diis_errors) + + if isinstance(t1, torch.Tensor): + # Build error matrix B + if self.precision == 'DP': + B = torch.ones((self.diis_size + 1, self.diis_size + 1), dtype=torch.float64, device=self.device1) * -1 + elif self.precision == 'SP': + B = torch.ones((self.diis_size + 1, self.diis_size + 1), dtype=torch.float32, device=self.device1) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = torch.dot(e1, e1) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = torch.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= torch.abs(B[:-1, :-1]).max() + + # Build residual vector + if self.precision == 'DP': + resid = torch.zeros((self.diis_size + 1), dtype=torch.float64, device=self.device1) + elif self.precision == 'SP': + resid = torch.zeros((self.diis_size + 1), dtype=torch.float32, device=self.device1) + resid[-1] = -1 + + # Solve pulay equations + ci = torch.linalg.solve(B, resid) + + # Calculate new amplitudes + t1 = torch.zeros_like(self.oldt1) + t2 = torch.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += torch.real(ci[num] * self.diis_vals_t1[num + 1]) + t2 += torch.real(ci[num] * self.diis_vals_t2[num + 1]) + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1.clone() + self.oldt2 = t2.clone() + + else: + # Build error matrix B + if self.precision == 'DP': + B = np.ones((self.diis_size + 1, self.diis_size + 1)) * -1 + elif self.precision == 'SP': + B = np.ones((self.diis_size + 1, self.diis_size + 1), dtype=np.float32) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = np.dot(e1, e1) + print("diis_errors",self.diis_errors[n1]) + print("B", B.shape) + print("n1","e1",n1,e1.shape) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = np.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= np.abs(B[:-1, :-1]).max() + + # Build residual vector + if self.precision == 'DP': + resid = np.zeros(self.diis_size + 1) + elif self.precision == 'SP': + resid = np.zeros((self.diis_size + 1), dtype=np.float32) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B, resid) + + # Calculate new amplitudes + t1 = np.zeros_like(self.oldt1) + t2 = np.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += ci[num] * self.diis_vals_t1[num + 1] + t2 += ci[num] * self.diis_vals_t2[num + 1] + + # Save extrapolated amplitudes to old_t amplitudes + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + + return t1, t2 + +class helper_ldiis(object): + def __init__(self, no,t1_ii, t2_ij, max_diis, Local): + oldt1 = [] + oldt2 = [] + + + self.Local = Local + self.QL = self.Local.QL + self.dim = self.Local.dim + + self.tmp = 13 + + t1 = np.zeros((no,self.dim[self.tmp])) + t2 = np.zeros((no,no,self.dim[self.tmp], self.dim[self.tmp])) + + for i in range(no): + ii = i*no + i + + S = self.QL[self.tmp].T @ self.QL[ii] + t1[i] = t1_ii[i].copy() @ S.T + + for j in range(no): + ij = i*no + j + + S = self.QL[self.tmp].T @ self.QL[ij] + t2[i,j] = S @ t2_ij[ij].copy() @ S.T + + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + #print("t2 init", self.oldt2[0,0]) + self.diis_vals_t1 = [t1.copy()] + self.diis_vals_t2 = [t2.copy()] + #print("t2 init", self.diis_vals_t2[-1][0,0]) + self.diis_errors = [] + self.diis_size = 0 + self.max_diis = max_diis + + def add_error_vector(self, no,t1_ii, t2_ij): + # Add DIIS vectors + + t1 = np.zeros((no,self.dim[self.tmp])) + t2 = np.zeros((no,no,self.dim[self.tmp], self.dim[self.tmp])) + + for i in range(no): + ii = i*no + i + + S = self.QL[self.tmp].T @ self.QL[ii] + t1[i] = t1_ii[i].copy() @ S.T + + for j in range(no): + ij = i*no + j + + S = self.QL[self.tmp].T @ self.QL[ij] + t2[i,j] = S @ t2_ij[ij].copy() @ S.T + + self.diis_vals_t1.append(t1.copy()) + self.diis_vals_t2.append(t2.copy()) + # Add new error vectors + error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() + error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() + self.diis_errors.append(np.concatenate((error_t1, error_t2))) + #print("previous t2_add", self.oldt2[0,0]) + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + #print("incoming t2_add",self.oldt2[0,0]) + + def extrapolate(self, no, t1_ii, t2_ij): + + if (self.max_diis == 0): + return t1_ii, t2_ij + + # Limit size of DIIS vector + if (len(self.diis_errors) > self.max_diis): + del self.diis_vals_t1[0] + del self.diis_vals_t2[0] + del self.diis_errors[0] + + self.diis_size = len(self.diis_errors) + + B = np.ones((self.diis_size + 1, self.diis_size + 1)) * -1 + B[-1, -1] = 0 + + for n1, e1 in enumerate(self.diis_errors): + B[n1, n1] = np.dot(e1, e1) + for n2, e2 in enumerate(self.diis_errors): + if n1 >= n2: + continue + B[n1, n2] = np.dot(e1, e2) + B[n2, n1] = B[n1, n2] + + B[:-1, :-1] /= np.abs(B[:-1, :-1]).max() + + #print("B matrix", B) + + # Build residual vector + resid = np.zeros(self.diis_size + 1) + resid[-1] = -1 + + # Solve pulay equations + ci = np.linalg.solve(B, resid) + #print("ci", ci) + # Calculate new amplitudes + t1 = np.zeros_like(self.oldt1) + t2 = np.zeros_like(self.oldt2) + for num in range(self.diis_size): + t1 += ci[num] * self.diis_vals_t1[num + 1] + t2 += ci[num] * self.diis_vals_t2[num + 1] + + # Save extrapolated amplitudes to old_t amplitudes + #print("incoming t2_extra", self.oldt2[0,0]) + self.oldt1 = t1.copy() + self.oldt2 = t2.copy() + #print("updated t2_extra", t2[0,0]) + #print("should updated t2_extra", self.oldt2[0,0]) + t1_ii = [] + t2_ij = [] + + for i in range(no): + ii = i*no + i + + S = self.QL[self.tmp].T @ self.QL[ii] + t1_ii.append(t1[i] @ S) + + for j in range(no): + ij = i*no + j + + S = self.QL[self.tmp].T @ self.QL[ij] + t2_ij.append(S.T @ t2[i,j] @ S) + + return t1_ii, t2_ij + +class cc_contract(object): + """ + A wrapper for opt_einsum.contract with tensors stored on CPU and/or GPU. + """ + def __init__(self, device='CPU'): + """ + Parameters + ---------- + device: string + initiated in ccwfn object, default: 'CPU' + + Returns + ------- + None + """ + self.device = device + if self.device == 'GPU': + # torch.device is an object representing the device on which torch.Tensor is or will be allocated. + self.device0 = torch.device('cpu') + self.device1 = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') + + def __call__(self, subscripts, *operands): + """ + Parameters + ---------- + subscripts: string + specify the subscripts for summation (same format as numpy.einsum) + *operands: list of array_like + the arrays/tensors for the operation + + Returns + ------- + An ndarray/torch.Tensor that is calculated based on Einstein summation convention. + """ + if self.device == 'CPU': + return opt_einsum.contract(subscripts, *operands) + elif self.device == 'GPU': + # Check the type and allocation of the tensors + # Transfer the copy from CPU to GPU if needed (for ERI) + input_list = list(operands) + for i in range(len(input_list)): + if (not input_list[i].is_cuda): + input_list[i] = input_list[i].to(self.device1) + #print(len(input_list), type(input_list[0]), type(input_list[1])) + output = opt_einsum.contract(subscripts, *input_list) + del input_list + return output + From 917f149089afc1cf8ab585ac87038d32b820267c Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Tue, 21 Mar 2023 13:46:52 -0400 Subject: [PATCH 6/8] updated utils --- pycc/utils.py | 151 ---------------------- pycc/{utils_correct.py => utils_Bij.py} | 45 ++----- pycc/{utils_test.py => utils_Bproject.py} | 15 +-- 3 files changed, 19 insertions(+), 192 deletions(-) rename pycc/{utils_correct.py => utils_Bij.py} (88%) rename pycc/{utils_test.py => utils_Bproject.py} (97%) diff --git a/pycc/utils.py b/pycc/utils.py index f99f36f..807c4c1 100644 --- a/pycc/utils.py +++ b/pycc/utils.py @@ -38,7 +38,6 @@ def add_error_vector(self, t1, t2): # Add DIIS vectors self.diis_vals_t1.append(t1.copy()) self.diis_vals_t2.append(t2.copy()) - # Add new error vectors error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() @@ -112,9 +111,6 @@ def extrapolate(self, t1, t2): if n1 >= n2: continue B[n1, n2] = np.dot(e1, e2) - print("B",B) - print("e1", n1,e1) - print("e2",n2,e2) B[n2, n1] = B[n1, n2] B[:-1, :-1] /= np.abs(B[:-1, :-1]).max() @@ -142,153 +138,6 @@ def extrapolate(self, t1, t2): return t1, t2 -class helper_ldiis(object): - def __init__(self, no,t1_ii, t2_ij, max_diis): - oldt1 = [] - oldt2 = [] - diis_vals_t1 = [] - diis_vals_t2 = [] - for i in range(no): - oldt1.append(t1_ii[i].copy()) - for j in range(no): - ij = i*no + j - oldt2.append(t2_ij[ij].copy()) - - self.oldt1 = oldt1 - self.oldt2 = oldt2 - diis_vals_t1.append(self.oldt1.copy()) - diis_vals_t2.append(self.oldt2.copy()) - self.diis_vals_t1 = diis_vals_t1 - self.diis_vals_t2 = diis_vals_t2 - #print("diis t2 initial", self.diis_vals_t2) - self.diis_errors_t1 = [] - self.diis_errors_t2 = [] - self.diis_size = 0 - self.max_diis = max_diis - - def add_error_vector(self, no,t1_ii, t2_ij): - # Add DIIS vectors - self.addt1 = [] - self.addt2 = [] - - for i in range(no): - ii = i*no + i - - self.addt1.append(t1_ii[i]) - - for j in range(no): - ij = i*no + j - - self.addt2.append(t2_ij[ij]) - #if ij == 24: - #print("t2 compare", ij, t2_ij[ij]) - - self.diis_vals_t1.append(self.addt1.copy()) - self.diis_vals_t2.append(self.addt2.copy()) - #for n1, e1 in enumerate(self.diis_vals_t2): - #print("n1", n1, "e1", e1[24]) - - #print("oldt2", self.oldt2[ij]) - - # Add new error vectors - error_t1 = [] - error_t2 = [] - for i in range(no): - error_t1.append((self.diis_vals_t1[-1][i] - self.oldt1[i]).ravel()) - for j in range(no): - ij = i*no + j - - error_t2.append((self.diis_vals_t2[-1][ij] - self.oldt2[ij]).ravel()) - if ij ==2: - print(error_t2[ij]) - print("len of error_t2", len(error_t2[ij]), self.oldt2[ij].shape) - #if ij == 24: - #print("diis t2 last", ij, self.diis_vals_t2[-1][ij]) - #print("oldt2 last", ij, self.oldt2[ij]) - #print(error_t2[24]) - self.diis_errors_t1.append(error_t1) - self.diis_errors_t2.append(error_t2) - self.oldt1 = self.addt1.copy() - self.oldt2 = self.addt2.copy() - - def extrapolate(self, no,t1_ii, t2_ij): - - if (self.max_diis == 0): - return t1_ii, t2_ij - - # Limit size of DIIS vector - if (len(self.diis_errors_t1) > self.max_diis): - del self.diis_vals_t1[0] - del self.diis_vals_t2[0] - del self.diis_errors_t1[0] - del self.diis_errors_t2[0] - - self.diis_size_t1 = len(self.diis_errors_t1) - self.diis_size_t2 = len(self.diis_errors_t2) - B_t1 = np.ones((self.diis_size_t1 + 1, self.diis_size_t1 + 1)) * -1 - B_t1[-1, -1] = 0 - B_t2 = np.ones((self.diis_size_t2 + 1, self.diis_size_t2 + 1)) * -1 - B_t2[-1, -1] = 0 - - t1_ii = np.zeros_like(self.oldt1) - t2_ij = np.zeros_like(self.oldt2) - for i in range(no): - - for n1, e1 in enumerate(self.diis_errors_t1): - B_t1[n1, n1] = np.dot(e1[i], e1[i]) - for n2, e2 in enumerate(self.diis_errors_t1): - if n1 >= n2: - continue - B_t1[n1, n2] = np.dot(e1[i], e2[i]) - B_t1[n2, n1] = B_t1[n1, n2] - B_t1[:-1, :-1] /= np.abs(B_t1[:-1, :-1]).max() - print("B_t1", B_t1) - # Build residual vector - resid = np.zeros(self.diis_size_t1 + 1) - resid[-1] = -1 - - # Solve pulay equations - ci = np.linalg.solve(B_t1, resid) - - # Calculate new amplitudes - for num in range(self.diis_size_t1): - t1_ii[i] += ci[num] * self.diis_vals_t1[num + 1][i] - - # Save extrapolated amplitudes to old_t amplitudes - for j in range(no): - ij = i*no + j - - for n1, e1 in enumerate(self.diis_errors_t2): - B_t2[n1, n1] = np.dot(e1[ij], e1[ij]) - for n2, e2 in enumerate(self.diis_errors_t2): - if n1 >= n2: - continue - B_t2[n1, n2] = np.dot(e1[ij], e2[ij]) - print("B_t2",B_t2) - print("e1", n1, e1[ij]) - print("e2", n2, e2[ij]) - B_t2[n2, n1] = B_t2[n1, n2] - print("B_t2 before", B_t2) - B_t2[:-1, :-1] /= np.abs(B_t2[:-1, :-1]).max() - print("B_t2", B_t2) - # Build residual vector - resid = np.zeros(self.diis_size_t2 + 1) - resid[-1] = -1 - - # Solve pulay equations - ci = np.linalg.solve(B_t2, resid) - #print("ci", ci, "ci.shape", ci.shape) - - # Calculate new amplitudes - for num in range(self.diis_size_t2): - t2_ij[ij] += ci[num] * self.diis_vals_t2[num + 1][ij] - #print("old", self.oldt2[24]) - #print("new", t2_ij[24]) - self.oldt1 = t1_ii.copy() - self.oldt2 = t2_ij.copy() - - return t1_ii, t2_ij - class cc_contract(object): """ A wrapper for opt_einsum.contract with tensors stored on CPU and/or GPU. diff --git a/pycc/utils_correct.py b/pycc/utils_Bij.py similarity index 88% rename from pycc/utils_correct.py rename to pycc/utils_Bij.py index faa2f3a..1257bf9 100644 --- a/pycc/utils_correct.py +++ b/pycc/utils_Bij.py @@ -151,10 +151,14 @@ def __init__(self, no,t1_ii, t2_ij, max_diis): oldt2 = [] diis_vals_t1 = [] diis_vals_t2 = [] + for i in range(no): + oldt1.append(t1_ii[i].copy()) + for j in range(no): ij = i*no + j + oldt2.append(t2_ij[ij].copy()) self.oldt1 = oldt1.copy() @@ -163,7 +167,6 @@ def __init__(self, no,t1_ii, t2_ij, max_diis): diis_vals_t2.append(self.oldt2.copy()) self.diis_vals_t1 = diis_vals_t1 self.diis_vals_t2 = diis_vals_t2 - #print("diis t2 initial", self.diis_vals_t2) self.diis_errors_t1 = [] self.diis_errors_t2 = [] self.diis_size = 0 @@ -183,44 +186,25 @@ def add_error_vector(self, no,t1_ii, t2_ij): ij = i*no + j self.addt2.append(t2_ij[ij].copy()) - #if ij == 2: - #print("new t2 part 1", ij, t2_ij[ij]) self.diis_vals_t1.append(self.addt1.copy()) self.diis_vals_t2.append(self.addt2.copy()) - #print("new t2", self.diis_vals_t2[-1][2]) - #for n1, e1 in enumerate(self.diis_vals_t2): - #print("n1", n1, "e1", e1[24]) - - #print("oldt2", self.oldt2[2]) # Add new error vectors error_t1 = [] error_t2 = [] for i in range(no): - #error_t1.append((self.diis_vals_t1[-1][i] - self.oldt1[i]).ravel()) error_t1.append((self.diis_vals_t1[-1][i] - self.diis_vals_t1[len(self.diis_vals_t1) -2][i]).ravel()) + for j in range(no): ij = i*no + j - #error_t2.append((self.diis_vals_t2[-1][ij] - self.oldt2[ij]).ravel()) error_t2.append((self.diis_vals_t2[-1][ij] - self.diis_vals_t2[len(self.diis_vals_t2) -2][ij]).ravel()) - #if ij ==2: - #print("new t2", self.diis_vals_t2[-1][ij]) - #print("oldt2",self.oldt2[ij]) - #print("this should be the old t2", self.diis_vals_t2[len(self.diis_vals_t2) -2][ij]) - #print("diff", error_t2[ij]) - #print("len of error_t2", len(error_t2[ij]), self.oldt2[ij].shape) - #if ij == 24: - #print("diis t2 last", ij, self.diis_vals_t2[-1][ij]) - #print("oldt2 last", ij, self.oldt2[ij]) - #print(error_t2[24]) + self.diis_errors_t1.append(error_t1) - #print("length of error", len(self.diis_errors_t1)) self.diis_errors_t2.append(error_t2) self.oldt1 = self.diis_vals_t1[-1].copy() self.oldt2 = self.diis_vals_t2[-1].copy() - #print("this should be old t2", self.oldt2[2]) def extrapolate(self, no,t1_ii, t2_ij): @@ -252,12 +236,13 @@ def extrapolate(self, no,t1_ii, t2_ij): continue B_t1[n1, n2] = np.dot(e1[i], e2[i]) B_t1[n2, n1] = B_t1[n1, n2] + B_t1[:-1, :-1] /= np.abs(B_t1[:-1, :-1]).max() - #print("B_t1", B_t1) # Build residual vector resid = np.zeros(self.diis_size_t1 + 1) resid[-1] = -1 + #Error arise due to singular matrix # Solve pulay equations ci = np.linalg.solve(B_t1, resid) @@ -275,26 +260,22 @@ def extrapolate(self, no,t1_ii, t2_ij): if n1 >= n2: continue B_t2[n1, n2] = np.dot(e1[ij], e2[ij]) - #print("B_t2",B_t2) - #print("e1", n1, e1[ij]) - #print("e2", n2, e2[ij]) B_t2[n2, n1] = B_t2[n1, n2] - #print("B_t2 before", B_t2) + B_t2[:-1, :-1] /= np.abs(B_t2[:-1, :-1]).max() - #print("B_t2", B_t2) + # Build residual vector resid = np.zeros(self.diis_size_t2 + 1) resid[-1] = -1 - + + #Error arise due to singular matrix # Solve pulay equations ci = np.linalg.solve(B_t2, resid) - #print("ci", ci, "ci.shape", ci.shape) # Calculate new amplitudes for num in range(self.diis_size_t2): t2_ij[ij] += ci[num] * self.diis_vals_t2[num + 1][ij] - #print("old", self.oldt2[24]) - #print("new", t2_ij[24]) + self.oldt1 = t1_ii.copy() self.oldt2 = t2_ij.copy() diff --git a/pycc/utils_test.py b/pycc/utils_Bproject.py similarity index 97% rename from pycc/utils_test.py rename to pycc/utils_Bproject.py index 30c778a..9fe4e96 100644 --- a/pycc/utils_test.py +++ b/pycc/utils_Bproject.py @@ -150,16 +150,17 @@ def __init__(self, no,t1_ii, t2_ij, max_diis, Local): oldt1 = [] oldt2 = [] - self.Local = Local self.QL = self.Local.QL self.dim = self.Local.dim + #PNO space to project to self.tmp = 13 + #Project all pno space to one particular space (tmp) t1 = np.zeros((no,self.dim[self.tmp])) t2 = np.zeros((no,no,self.dim[self.tmp], self.dim[self.tmp])) - + for i in range(no): ii = i*no + i @@ -174,10 +175,8 @@ def __init__(self, no,t1_ii, t2_ij, max_diis, Local): self.oldt1 = t1.copy() self.oldt2 = t2.copy() - #print("t2 init", self.oldt2[0,0]) self.diis_vals_t1 = [t1.copy()] self.diis_vals_t2 = [t2.copy()] - #print("t2 init", self.diis_vals_t2[-1][0,0]) self.diis_errors = [] self.diis_size = 0 self.max_diis = max_diis @@ -185,6 +184,7 @@ def __init__(self, no,t1_ii, t2_ij, max_diis, Local): def add_error_vector(self, no,t1_ii, t2_ij): # Add DIIS vectors + #Porject all pno space to one particular space (tmp) t1 = np.zeros((no,self.dim[self.tmp])) t2 = np.zeros((no,no,self.dim[self.tmp], self.dim[self.tmp])) @@ -206,10 +206,8 @@ def add_error_vector(self, no,t1_ii, t2_ij): error_t1 = (self.diis_vals_t1[-1] - self.oldt1).ravel() error_t2 = (self.diis_vals_t2[-1] - self.oldt2).ravel() self.diis_errors.append(np.concatenate((error_t1, error_t2))) - #print("previous t2_add", self.oldt2[0,0]) self.oldt1 = t1.copy() self.oldt2 = t2.copy() - #print("incoming t2_add",self.oldt2[0,0]) def extrapolate(self, no, t1_ii, t2_ij): @@ -245,7 +243,7 @@ def extrapolate(self, no, t1_ii, t2_ij): # Solve pulay equations ci = np.linalg.solve(B, resid) - #print("ci", ci) + # Calculate new amplitudes t1 = np.zeros_like(self.oldt1) t2 = np.zeros_like(self.oldt2) @@ -257,11 +255,10 @@ def extrapolate(self, no, t1_ii, t2_ij): #print("incoming t2_extra", self.oldt2[0,0]) self.oldt1 = t1.copy() self.oldt2 = t2.copy() - #print("updated t2_extra", t2[0,0]) - #print("should updated t2_extra", self.oldt2[0,0]) t1_ii = [] t2_ij = [] + #reverting back to their respective pno space for i in range(no): ii = i*no + i From 03bd8bea49a5b2a9a1c3a9114796d6cc8c903e14 Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Tue, 21 Mar 2023 13:48:11 -0400 Subject: [PATCH 7/8] updated utils --- pycc/ccwfn.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py index c583830..bf3f250 100644 --- a/pycc/ccwfn.py +++ b/pycc/ccwfn.py @@ -10,11 +10,11 @@ import time import numpy as np import torch -from utils_correct import helper_diis, cc_contract -from hamiltonian import Hamiltonian -from local import Local -from cctriples import t_tjl, t3c_ijk -from lccwfn import lccwfn +from .utils_Bij import helper_diis, cc_contract #in progress +from .hamiltonian import Hamiltonian +from .local import Local +from .cctriples import t_tjl, t3c_ijk +from .lccwfn import lccwfn class ccwfn(object): """ From 3266c5bfe880e9d26fee69f7d7a46a64c50d99cd Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Tue, 21 Mar 2023 14:05:18 -0400 Subject: [PATCH 8/8] merged main to lcc_diis --- pycc/lccwfn.py | 115 ++----------------------------------------------- 1 file changed, 3 insertions(+), 112 deletions(-) diff --git a/pycc/lccwfn.py b/pycc/lccwfn.py index 3e686fe..ff7cb5a 100644 --- a/pycc/lccwfn.py +++ b/pycc/lccwfn.py @@ -73,17 +73,10 @@ def __init__(self, o, v, no, nv, H, local, model, eref, Local): t2.append(-1* self.Local.ERIoovv[ij][i,j] / (self.eps[ij].reshape(1,-1) + self.eps[ij].reshape(-1,1) - self.H.F[i,i] - self.H.F[j,j])) -<<<<<<< HEAD self.t1 = t1 self.t2 = t2 -======= - - self.t1 = t1 - self.t2 = t2 - ->>>>>>> main def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis=1): """ Parameters @@ -116,11 +109,8 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis #self.r2_t = Timer("r2") #self.energy_t = Timer("energy") -<<<<<<< HEAD ldiis = helper_ldiis(self.no, self.t1, self.t2, max_diis, self.Local) -======= - #ldiis = helper_ldiis(self.t1, self.t2, max_diis) ->>>>>>> main + elcc = self.lcc_energy(self.Local.Fov,self.Local.Loovv,self.t1, self.t2) print("CC Iter %3d: lCC Ecorr = %.15f dE = % .5E MP2" % (0,elcc,-elcc)) @@ -167,15 +157,9 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis #print(Timer.timers) return elcc -<<<<<<< HEAD - ldiis.add_error_vector(self.no, self.t1,self.t2) + ldiis.add_error_vector(self.t1,self.t2) if niter >= start_diis: - self.t1, self.t2 = ldiis.extrapolate(self.no, self.t1, self.t2) -======= - #ldiis.add_error_vector(self.t1,self.t2) - #if niter >= start_diis: - #self.t1, self.t2 = ldiis.extrapolate(self.t1, self.t2) ->>>>>>> main + self.t1, self.t2 = ldiis.extrapolate(self.t1, self.t2) def local_residuals(self, t1, t2): """ @@ -208,15 +192,9 @@ def local_residuals(self, t1, t2): Wmbej = self.build_Wmbej(Wmbej, ERI, L, self.Local.ERIoovo, self.Local.Sijnn, self.Local.Sijnj, self.Local.Sijjn, t1, t2) Wmbje, Wmbie = self.build_Wmbje(Wmbje, Wmbie, ERI, self.Local.ERIooov, self.Local.Sijnn, self.Local.Sijin, self.Local.Sijjn, t1, t2) -<<<<<<< HEAD - r1 = self.r_T1(r1, self.Local.Fov , ERI, L, self.Local.Loovo, self.Local.Sijmm, self.Local.Sijim, self.Local.Sijmn, - t1, t2, Fae, Fmi, Fme) - r2 = self.r_T2(r2, ERI, self.Local.ERIoovv, self.Local.ERIvvvv, self.Local.ERIovoo, self.Local.Sijmm, self.Local.Sijim, -======= r1 = self.r_T1(r1, self.Local.Fov , ERI, L, self.Local.Loovo, self.Local.Sijmm, self.Local.Sijim, self.Local.Sijmn, t1, t2, Fae, Fmi, Fme) r2 = self.r_T2(r2, ERI, self.Local.ERIoovv, self.Local.ERIvvvv, self.Local.ERIovoo, self.Local.Sijmm, self.Local.Sijim, ->>>>>>> main self.Local.Sijmj, self.Local.Sijnn, self.Local.Sijmn, t1, t2, Fae ,Fmi, Fme, Wmnij, Zmbij, Wmbej, Wmbje, Wmbie) return r1, r2 @@ -238,7 +216,6 @@ def build_Fae(self, Fae_ij, L, Fvv, Fov, Sijmm, Sijmn, t1, t2): for n in range(self.no): mn = m *self.no +n ijmn = ij*(self.no**2) + mn -<<<<<<< HEAD tmp = Sijmn[ijmn] @ t2[mn] tmp1 = QL[ij].T @ L[m,n,v,v] @@ -278,47 +255,6 @@ def build_Fae(self, Fae_ij, L, Fvv, Fov, Sijmm, Sijmn, t1, t2): tmp4 = tmp3_0 @ QL[nn] Fae -= 0.5 *contract('a,F,eF->ae', tmp, t1[n], tmp4) -======= - - tmp = Sijmn[ijmn] @ t2[mn] - tmp1 = QL[ij].T @ L[m,n,v,v] - tmp1 = tmp1 @ QL[mn] - Fae -= tmp @ tmp1.T - - Fae_ij.append(Fae) - else: - for ij in range(self.no*self.no): - i = ij // self.no - j = ij % self.no - - Fae = Fvv[ij].copy() - - for m in range(self.no): - mm = m*self.no + m - ijm = ij*(self.no) + m - - tmp = Sijmm[ijm] @ t1[m] - Fae -= 0.5* contract('e,a->ae',Fov[ij][m],tmp) - - tmp1 = contract('abc,aA->Abc',L[m,v,v,v], QL[ij]) - tmp1 = contract('Abc,bB->ABc',tmp1, QL[mm]) - tmp1 = contract('ABc,cC->ABC',tmp1, QL[ij]) - Fae += contract('F,aFe->ae',t1[m],tmp1) - - for n in range(self.no): - mn = m *self.no +n - nn = n*self.no + n - ijmn = ij*(self.no**2) + mn - - tmp2 = Sijmn[ijmn] @ t2[mn] - tmp3_0 = QL[ij].T @ L[m,n,v,v] - tmp3_1 = tmp3_0 @ QL[mn] - Fae -= tmp2 @ tmp3_1.T - - tmp4 = tmp3_0 @ QL[nn] - Fae -= 0.5 *contract('a,F,eF->ae', tmp, t1[n], tmp4) - ->>>>>>> main Fae_ij.append(Fae) #self.fae_t.stop() return Fae_ij @@ -334,17 +270,6 @@ def build_Fmi(self, o, F, L, Fov, Looov, Loovv, t1, t2): for j in range(self.no): for n in range(self.no): jn = j*self.no + n -<<<<<<< HEAD - - Fmi[:,j] += contract('EF,mEF->m',t2[jn],Loovv[jn][:,n,:,:]) - else: - for j in range(self.no): - jj = j*self.no +j - for n in range(self.no): - jn = j*self.no + n - nn = n*self.no + n - -======= Fmi[:,j] += contract('EF,mEF->m',t2[jn],Loovv[jn][:,n,:,:]) else: @@ -354,7 +279,6 @@ def build_Fmi(self, o, F, L, Fov, Looov, Loovv, t1, t2): jn = j*self.no + n nn = n*self.no + n ->>>>>>> main Fmi[:,j] += 0.5 * contract('e,me->m', t1[j], Fov[jj]) Fmi[:,j] += contract('e,me->m',t1[n],Looov[nn][:,n,j]) Fmi[:,j] += contract('EF,mEF->m',t2[jn],Loovv[jn][:,n,:,:]) @@ -402,21 +326,6 @@ def build_Wmnij(self, o, ERI, ERIooov, ERIoovo, ERIoovv, t1, t2): for i in range(self.no): for j in range(self.no): ij = i*self.no + j -<<<<<<< HEAD - - Wmnij[:,:,i,j] += contract('ef,mnef->mn',t2[ij], ERIoovv[ij]) - else: - for i in range(self.no): - for j in range(self.no): - ij = i*self.no + j - ii = i*self.no + i - jj = j*self.no + j - - Wmnij[:,:,i,j] += contract('E,mnE->mn', t1[j], ERIooov[jj][:,:,i,:]) - Wmnij[:,:,i,j] += contract('E,mnE->mn', t1[i], ERIoovo[ii][:,:,:,j]) - Wmnij[:,:,i,j] += contract('ef,mnef->mn',t2[ij], ERIoovv[ij]) - -======= Wmnij[:,:,i,j] += contract('ef,mnef->mn',t2[ij], ERIoovv[ij]) else: @@ -430,7 +339,6 @@ def build_Wmnij(self, o, ERI, ERIooov, ERIoovo, ERIoovv, t1, t2): Wmnij[:,:,i,j] += contract('E,mnE->mn', t1[i], ERIoovo[ii][:,:,:,j]) Wmnij[:,:,i,j] += contract('ef,mnef->mn',t2[ij], ERIoovv[ij]) ->>>>>>> main tmp = contract('aA,mnab->mnAb', QL[ii], ERI[o,o,v,v]) tmp = contract('bB,mnAb->mnAB', QL[jj], tmp) Wmnij[:,:,i,j] += contract('e,f,mnef->mn', t1[i], t1[j], tmp) @@ -443,22 +351,6 @@ def build_Zmbij(self, Zmbij_ij, ERI, ERIovvv, t1, t2): o = self.o v = self.v QL = self.QL -<<<<<<< HEAD - - if self.model == 'CCD': - return - else: - for ij in range(self.no*self.no): - i = ij // self.no - j = ij % self.no - ii = i*self.no + i - jj = j*self.no + j - - Zmbij = np.zeros((self.no,self.Local.dim[ij])) - - Zmbij = contract('mbef,ef->mb', ERIovvv[ij], t2[ij]) - -======= if self.model == 'CCD': return @@ -473,7 +365,6 @@ def build_Zmbij(self, Zmbij_ij, ERI, ERIovvv, t1, t2): Zmbij = contract('mbef,ef->mb', ERIovvv[ij], t2[ij]) ->>>>>>> main tmp = contract('iabc,aA->iAbc',ERI[o,v,v,v], QL[ij]) tmp = contract('iAbc,bB->iABc',tmp, QL[ii]) tmp = contract('iABc,cC->iABC',tmp, QL[jj])