diff --git a/exploratory/citools/lasci_ominus1.py b/exploratory/citools/lasci_ominus1.py index 26049d6d..bc6f425c 100644 --- a/exploratory/citools/lasci_ominus1.py +++ b/exploratory/citools/lasci_ominus1.py @@ -33,7 +33,7 @@ def _n_m_s (dm1s, dm2s, _print_fn=print): _print_fn (',, = %f, %f, %f', n, m, ss) def kernel (fci, h1, h2, norb, nelec, norb_f=None, ci0_f=None, - tol=1e-8, gtol=1e-4, max_cycle=15000, + tol=1e-8, gtol=1e-4, max_cycle=15000, orbsym=None, wfnsym=None, ecore=0, **kwargs): if norb_f is None: norb_f = getattr (fci, 'norb_f', [norb]) @@ -183,7 +183,7 @@ def contract_ss (fci, ci, norb, nelec): return fcivec def transform_ci_for_orbital_rotation (fci, ci, norb, nelec, umat): - fcivec = np.zeros_like (ci) + fcivec = np.zeros_like (ci) for ne in product (range (norb+1), repeat=2): c = np.squeeze (fockspace.fock2hilbert (ci, norb, ne)) c = fci_addons.transform_ci_for_orbital_rotation (c, norb, ne, umat) @@ -257,7 +257,7 @@ def unpack (self, x_): x[self.var_mask] = x_[:] xconstr, x = x[:self.nconstr], x[self.nconstr:] - xcc = x[:self.uop.ngen_uniq] + xcc = x[:self.uop.ngen_uniq] x = x[self.uop.ngen_uniq:] xci = [] @@ -332,7 +332,7 @@ def hc_x (self, x, h): huc = self.contract_h2 (h, uc) uhuc = self.uop (huc, transpose=True) return c, uc, huc, uhuc, c_f - + def contract_h2 (self, h, ci, norb=None): if norb is None: norb = self.norb hci = h[0] * ci @@ -343,7 +343,7 @@ def contract_h2 (self, h, ci, norb=None): hc = direct_spin1.contract_2e (h2eff, ci_h, norb, nelec) hci += np.squeeze (fockspace.hilbert2fock (hc, norb, nelec)) return hci - + def dp_ci (self, ci_f): norb, norb_f = self.norb, self.norb_f ci = np.ones ([1,1], dtype=ci_f[0].dtype) @@ -357,7 +357,7 @@ def constr_h (self, xconstr, h): x = xconstr[0] norb, nelec = self.norb, self.nelec h = [h[0] - (x*nelec), h[1] + (x*np.eye (self.norb)), h[2]] - return h + return h def rotate_ci0 (self, xci_f): ci0, norb = self.ci_f, self.norb @@ -408,12 +408,78 @@ def get_jac_t1 (self, x, h, c=None, huc=None, uhuc=None): xconstr, xcc, xci_f = self.unpack (x) self.uop.set_uniq_amps_(xcc) if (c is None) or (uhuc is None): - c, _, _, uhuc = self.hc_x (x, h)[:4] + c, _, _, uhuc = self.hc_x (x, h)[:4] for duc, uhuc_i in zip (self.uop.gen_deriv1 (c, _full=False), self.uop.gen_partial (uhuc)): g.append (2*duc.ravel ().dot (uhuc_i.ravel ())) g = self.uop.product_rule_pack (g) return np.asarray (g) + + def get_grad_t1(self, x, h, c=None, huc=None, uhuc=None, epsilon=0.0): + """ + Compute the gradients and relevant indices based on input values. + + Parameters: + - x: array-like + Input array to unpack and set unique amplitudes. + + - h: array-like + Some form of input data to be used for computation. + + - c (optional): array-like + Precomputed value; if not provided, it will be computed using hc_x. + + - huc (optional): array-like + Precomputed value; if not provided, it will be computed using hc_x. + + - uhuc (optional): array-like + Precomputed value; if not provided, it will be computed using hc_x. + + - epsilon (optional): float, default=0.0 + Threshold value for considering a gradient. If epsilon is 0, all gradients are considered. + + Returns: + - tuple + all_g: list of all computed gradients + g: list of gradients above the epsilon threshold + gen_indices: list of indices representing a_idx and i_idx + a_idxs_lst: list of a_idx values + i_idxs_lst: list of i_idx values + len(a_idxs_lst): length of a_idx list + len(i_idxs_lst): length of i_idx list + """ + + g = [] + all_g = [] + + # Unpack and set unique amplitudes + xconstr, xcc, xci_f = self.unpack(x) + self.uop.set_uniq_amps_(xcc) + + # Compute 'c' and 'uhuc' if not provided + if (c is None) or (uhuc is None): + c, _, _, uhuc = self.hc_x(x, h)[:4] + + gen_indices = [] + a_idxs_lst = [] + i_idxs_lst = [] + # print("self.uop.init_a_idxs[i]",self.uop.init_a_idxs) + for i, (duc, uhuc_i) in enumerate(zip(self.uop.gen_deriv1(c, _full=False), self.uop.gen_partial(uhuc))): + gradient = 2 * duc.ravel().dot(uhuc_i.ravel()) + all_g.append((gradient, i)) + + # Allow all gradients if epsilon is 0, else use the abs gradient condition + if epsilon == 0.0 or abs(gradient) > epsilon: + g.append((gradient, i)) + a_idx = self.uop.a_idxs[i] + i_idx = self.uop.i_idxs[i] + + gen_indices.append((a_idx, i_idx)) + a_idxs_lst.append(a_idx) + i_idxs_lst.append(i_idx) + + return all_g, g, gen_indices, a_idxs_lst, i_idxs_lst, len(a_idxs_lst), len(i_idxs_lst) + def get_jac_ci (self, x, h, uhuc=None, uci_f=None): # "uhuc": e^-T1 H e^T1 U|ci0> # "jacci": Jacobian elements for the CI degrees of freedom @@ -542,7 +608,7 @@ def print_x (self, x, h, _print_fn=print, ci_maxlines=10, jac=None): def finalize_(self): ''' Update self.ci_f and set the corresponding part of self.x to zero ''' xconstr, xcc, xci_f = self.unpack (self.x) - self.ci_f = self.rotate_ci0 (xci_f) + self.ci_f = self.rotate_ci0 (xci_f) for xc in xci_f: xc[:] = 0.0 self.x = self.pack (xconstr, xcc, xci_f) return self.x @@ -564,7 +630,7 @@ class FCISolver (direct_spin1.FCISolver): def build_psi (self, *args, **kwargs): return LASUCCTrialState (self, *args, **kwargs) def save_psi (self, fname, psi): - psi_raw = np.concatenate ([psi.x,] + psi_raw = np.concatenate ([psi.x,] + [c.ravel () for c in psi.ci_f]) np.save (fname, psi_raw) def load_psi (self, fname, norb, nelec, norb_f=None, **kwargs): @@ -584,5 +650,3 @@ def get_uop (self, norb, norb_f): for i,j in zip (np.cumsum (norb_f)-norb_f, np.cumsum(norb_f)): freeze_mask[i:j,i:j] = True return get_uccs_op (norb, freeze_mask=freeze_mask) - - diff --git a/my_pyscf/grad/lasscf.py b/my_pyscf/grad/lasscf.py index 9927e32d..0a14a411 100644 --- a/my_pyscf/grad/lasscf.py +++ b/my_pyscf/grad/lasscf.py @@ -72,11 +72,7 @@ def grad_elec(las_grad, mo_coeff=None, ci=None, atmlst=None, verbose=None): mo_occ = mo_coeff[:,:nocc] mo_core = mo_coeff[:,:ncore] mo_cas = mo_coeff[:,ncore:nocc] - #print ("SV nocc, ncas, ncore = ", nocc, ncas, ncore) - #print ("mo_cas shape =", mo_cas.shape) - #print ("SV entering grad_elec in lasscf") lasdm1 = las.make_casdm1() - #print ("SV lasdm1 shape = ",lasdm1.shape) lasdm2 = las.make_casdm2() #casdm1, casdm2 = mc.fcisolver.make_rdm12(ci, ncas, nelecas) diff --git a/my_pyscf/mcscf/lasci.py b/my_pyscf/mcscf/lasci.py index e8e37d5d..ede696b1 100644 --- a/my_pyscf/mcscf/lasci.py +++ b/my_pyscf/mcscf/lasci.py @@ -1179,13 +1179,6 @@ def states_make_casdm2s_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, **k casdm2s.append (np.stack ([dm2aa, dm2ab, dm2bb], axis=1)) return casdm2s - #def make_casdm3_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, casdm3fr=None, **kwargs): - # if casdm3fr is None: casdm3fr = self.states_make_casdm3_sub (ci=ci, ncas_sub=ncas_sub, - # nelecas_sub=nelecas_sub, **kwargs) - # for dm3, box in zip(casdm3fr, self.fciboxes): - # casdm3_sub = np.einsum('rijklmn,r->ijklmn', dm3, box.weights) - # return casdm3_sub - def states_make_casdm3_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, **kwargs): ''' Spin-separated 3-RDMs in the MO basis for each subspace in sequence, currently this @@ -1201,6 +1194,35 @@ def states_make_casdm3_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, **kw casdm3.append (dm3_no) return casdm3 + def make_casdm3s_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, **kwargs): + '''Spin-separated 3-RDMs in the MO basis for each subspace in sequence using the PySCF make_rdm123s''' + if ci is None: ci = self.ci + if ncas_sub is None: ncas_sub = self.ncas_sub + if nelecas_sub is None: nelecas_sub = self.nelecas_sub + if ci is None: + return [np.zeros ((self.nroots,2,ncas,ncas,ncas,ncas,ncas,ncas)) for ncas in ncas_sub] + casdm3s = [] + for fcibox, ci_i, ncas, nelecas in zip (self.fciboxes, ci, ncas_sub, nelecas_sub): + dm3aaa, dm3aab, dm3abb, dm3bbb = fci.direct_spin1.make_rdm123s (np.array(ci_i), ncas, tuple(nelecas))[-1] + casdm3s.append (np.stack ([dm3aaa, dm3aab, dm3abb, dm3bbb], axis=0)) + return casdm3s + + #SV casdm3s_sub from external fake_rdm.py file + def make_casdm3s_sub_1 (self, ncas_sub=None, nelecas_sub=None, dm3s_sub_spinless=None, **kwargs): + '''Spin-separated 3-RDMs in the MO basis for each subspace in sequence''' + if ncas_sub is None: ncas_sub = self.ncas_sub + if nelecas_sub is None: nelecas_sub = self.nelecas_sub + dm3s_sub = dm3s_sub_spinless + casdm3s = [] + + for dm3, ncas in zip (dm3s_sub, ncas_sub): + dm3aaa = dm3[:ncas,:ncas,:ncas,:ncas,:ncas,:ncas] + dm3aab = dm3[:ncas,:ncas,:ncas,:ncas,ncas:,ncas:] + dm3abb = dm3[:ncas,:ncas,ncas:,ncas:,ncas:,ncas:] + dm3bbb = dm3[ncas:,ncas:,ncas:,ncas:,ncas:,ncas:] + casdm3s.append (np.stack ([dm3aaa, dm3aab, dm3abb, dm3bbb], axis=0)) + return casdm3s + def states_make_rdm1s (self, mo_coeff=None, ci=None, ncas_sub=None, nelecas_sub=None, casdm1rs=None, casdm1frs=None, **kwargs): ''' Get spin-separated, rootspace-separated, locally state-averaged 1-RDMs in the AO basis. @@ -1539,7 +1561,7 @@ def make_casdm2s (self, ci=None, ncas_sub=None, nelecas_sub=None, j = ncas_cum[isub+1] for spin in [0,1,2]: casdm2s[spin][i:j, i:j, i:j, i:j] = dm2[spin] - + # Off-diagonal for (isub1, dm1rs1), (isub2, dm1rs2) in combinations (enumerate (casdm1frs), 2): i = ncas_cum[isub1] @@ -1548,18 +1570,32 @@ def make_casdm2s (self, ci=None, ncas_sub=None, nelecas_sub=None, l = ncas_cum[isub2+1] dma1r, dmb1r = dm1rs1[:,0], dm1rs1[:,1] dma2r, dmb2r = dm1rs2[:,0], dm1rs2[:,1] - dm1r = dma1r + dmb1r - dm2r = dma2r + dmb2r + # Coulomb slice + casdm2s[0][i:j, i:j, k:l, k:l] = lib.einsum ('r,rij,rkl->ijkl', weights, dma1r, dma2r) + casdm2s[1][i:j, i:j, k:l, k:l] = lib.einsum ('r,rij,rkl->ijkl', weights, dma1r, dmb2r) # ab = ba + casdm2s[2][i:j, i:j, k:l, k:l] = lib.einsum ('r,rij,rkl->ijkl', weights, dmb1r, dmb2r) for spin in [0,1,2]: - # Coulomb slice - casdm2s[spin][i:j,i:j,k:l,k:l] = lib.einsum ('r,rij,rkl->ijkl', weights,dm1r,dm2r) - casdm2s[spin][k:l,k:l,i:j,i:j] = casdm2s[spin][i:j,i:j,k:l,k:l].transpose (2,3,0,1) - # Exchange slice - d2exc = (lib.einsum ('rij,rkl->rilkj', dma1r, dma2r) - + lib.einsum ('rij,rkl->rilkj', dmb1r, dmb2r)) - casdm2s[spin][i:j,k:l,k:l,i:j] -= np.tensordot (weights, d2exc, axes=1) - casdm2s[spin][k:l,i:j,i:j,k:l] = casdm2s[spin][i:j,k:l,k:l,i:j].transpose (1,0,3,2) - +#<<<<<<< HEAD + casdm2s[spin][k:l, k:l, i:j, i:j] = casdm2s[spin][i:j, i:j, k:l, k:l].transpose (2,3,0,1) + # Exchange slice + d2exc_aa = lib.einsum ('rij,rkl->rilkj', dma1r, dma2r) + d2exc_bb = lib.einsum ('rij,rkl->rilkj', dmb1r, dmb2r) + casdm2s[0][i:j, k:l, k:l, i:j] -= np.tensordot (weights, d2exc_aa, axes=1) + casdm2s[2][i:j, k:l, k:l, i:j] -= np.tensordot (weights, d2exc_bb, axes=1) + for spin in [0,2]: + casdm2s[spin][k:l, i:j, i:j, k:l] = casdm2s[spin][i:j, k:l, k:l, i:j].transpose (1,0,3,2) + +#======= +# # Coulomb slice +# casdm2s[spin][i:j,i:j,k:l,k:l] = lib.einsum ('r,rij,rkl->ijkl', weights,dm1r,dm2r) +# casdm2s[spin][k:l,k:l,i:j,i:j] = casdm2s[spin][i:j,i:j,k:l,k:l].transpose (2,3,0,1) +# # Exchange slice +# d2exc = (lib.einsum ('rij,rkl->rilkj', dma1r, dma2r) +# + lib.einsum ('rij,rkl->rilkj', dmb1r, dmb2r)) +# casdm2s[spin][i:j,k:l,k:l,i:j] -= np.tensordot (weights, d2exc, axes=1) +# casdm2s[spin][k:l,i:j,i:j,k:l] = casdm2s[spin][i:j,k:l,k:l,i:j].transpose (1,0,3,2) +# +#>>>>>>> upstream/master return casdm2s @@ -1657,6 +1693,7 @@ def make_casdm3 (self, ci=None, ncas_sub=None, nelecas_sub=None, ncas_cum = np.cumsum ([0] + ncas_sub.tolist ()) weights = self.weights casdm3 = np.zeros ((ncas,ncas,ncas,ncas,ncas,ncas)) + # Diagonal for isub, dm3 in enumerate (casdm3f): i = ncas_cum[isub] @@ -1700,13 +1737,21 @@ def _transpose_dm3_(i0, i1, i2, i3, j0, j1, j2, j3, k0, k1, k2, k3): casdm3[i:j, i:j, k:l, m:n, m:n, k:l] -= np.tensordot(weights, d2sigma, axes=1) _transpose_dm3_(i,j, i,j, k,l, m,n, m,n, k,l) +#<<<<<<< HEAD + # Term 3- 29f +#======= # Term 3- 29c +#>>>>>>> upstream/master d3sigma = (lib.einsum('r,rij,rkl,rmn->rinkjml',weights,dma1r,dma2r,dma3r) +lib.einsum('r,rij,rkl,rmn->rinkjml',weights,dmb1r,dmb2r,dmb3r)) casdm3[i:j, m:n, k:l, i:j, m:n, k:l] += np.tensordot (weights, d3sigma, axes=1) _transpose_dm3_(i,j, m,n, k,l, i,j, m,n, k,l) +#<<<<<<< HEAD + # Term 4- 29c +#======= # Term 4- 29d +#>>>>>>> upstream/master d4sigma = (lib.einsum('r,rij,rkl,rmn->rilkjmn',weights,dma1r,dma2r,dm3r) +lib.einsum('r,rij,rkl,rmn->rilkjmn',weights,dmb1r,dmb2r,dm3r)) casdm3[i:j, k:l, k:l, i:j, m:n, m:n] -= np.tensordot(weights, d4sigma, axes=1) @@ -1718,7 +1763,11 @@ def _transpose_dm3_(i0, i1, i2, i3, j0, j1, j2, j3, k0, k1, k2, k3): casdm3[i:j, m:n, k:l, k:l, m:n, i:j] -= np.tensordot(weights, d5sigma, axes=1) _transpose_dm3_(i,j, m,n, k,l, k,l, m,n, i,j) +#<<<<<<< HEAD + # Term 6- 29d +#======= # Term 6- 29f +#>>>>>>> upstream/master d6sigma = (lib.einsum('r,rij,rkl,rmn->rilknmj',weights,dma1r,dma2r,dma3r) +lib.einsum('r,rij,rkl,rmn->rilknmj',weights,dmb1r,dmb2r,dmb3r)) casdm3[i:j, k:l, k:l, m:n, m:n, i:j] += np.tensordot(weights, d6sigma, axes=1) @@ -1760,6 +1809,253 @@ def _transpose_dm3_(i0, i1, i2, i3, j0, j1, j2, j3, k0, k1, k2, k3): return casdm3 + def make_casdm3s (self, ci=None, ncas_sub=None, nelecas_sub=None, + casdm3r=None, casdm3f=None, casdm2frs=None, casdm3fr=None, casdm2r=None, casdm2f=None, + casdm1frs=None, casdm2fr=None, casdm3fs= None, **kwargs): + ''' Make the full-dimensional casdm3 spanning the collective active space ''' + if casdm3r is not None: + return np.einsum ('rijklmn,r->ijklmn', casdm3r, self.weights) + if ci is None: ci = self.ci + if ncas_sub is None: ncas_sub = self.ncas_sub + if nelecas_sub is None: nelecas_sub = self.nelecas_sub + if casdm1frs is None: casdm1frs = self.states_make_casdm1s_sub (ci=ci, + ncas_sub=ncas_sub, nelecas_sub=nelecas_sub) + if casdm2frs is None: casdm2frs = self.states_make_casdm2s_sub (ci=ci, + ncas_sub=ncas_sub, nelecas_sub=nelecas_sub) + if casdm2f is None: casdm2f = self.make_casdm2_sub (ci=ci, + ncas_sub=ncas_sub, nelecas_sub=nelecas_sub, casdm2fr=casdm2fr) + if casdm3f is None: casdm3f = self.states_make_casdm3_sub (ci=ci, + ncas_sub=ncas_sub, nelecas_sub=nelecas_sub) + if casdm3fs is None: casdm3fs = self.make_casdm3s_sub (ci=ci, + ncas_sub=ncas_sub, nelecas_sub=nelecas_sub) + + ncas = sum (ncas_sub) + ncas_cum = np.cumsum ([0] + ncas_sub.tolist ()) + weights = self.weights + casdm3s = np.zeros ((4,ncas,ncas,ncas,ncas,ncas,ncas)) + + # Diagonal + for isub, dm3 in enumerate (casdm3fs): + i = ncas_cum[isub] + j = ncas_cum[isub+1] + for spin in range(4): + casdm3s[spin][i:j, i:j, i:j, i:j, i:j, i:j] = dm3[spin] + + def _transpose_dm3_(dmt, i0, i1, i2, i3, j0, j1, j2, j3, k0, k1, k2, k3): + d = dmt[i0:i1, i2:i3, j0:j1, j2:j3, k0:k1, k2:k3] + dmt[i0:i1, i2:i3, k0:k1, k2:k3, j0:j1, j2:j3] = d.transpose (0,1,4,5,2,3) + dmt[k0:k1, k2:k3, i0:i1, i2:i3, j0:j1, j2:j3] = d.transpose (4,5,0,1,2,3) + dmt[k0:k1, k2:k3, j0:j1, j2:j3, i0:i1, i2:i3] = d.transpose (4,5,2,3,0,1) + dmt[j0:j1, j2:j3, k0:k1, k2:k3, i0:i1, i2:i3] = d.transpose (2,3,4,5,0,1) + dmt[j0:j1, j2:j3, i0:i1, i2:i3, k0:k1, k2:k3] = d.transpose (2,3,0,1,4,5) + + + # Off-diagonal + # First 6 terms - combs of f1+f2+f3: dm1rs1=f1, dm1rs2=f2, dm1rs3=f3 + from itertools import permutations + for (isub1,dm1rs1),(isub2,dm1rs2),(isub3,dm1rs3) in combinations (enumerate(casdm1frs), 3): + i = ncas_cum[isub1] + j = ncas_cum[isub1+1] + k = ncas_cum[isub2] + l = ncas_cum[isub2+1] + m = ncas_cum[isub3] + n = ncas_cum[isub3+1] + dma1r, dmb1r = dm1rs1[:,0], dm1rs1[:,1] + dma2r, dmb2r = dm1rs2[:,0], dm1rs2[:,1] + dma3r, dmb3r = dm1rs3[:,0], dm1rs3[:,1] + dm1r = dma1r + dmb1r + dm2r = dma2r + dmb2r + dm3r = dma3r + dmb3r + + # Term 1- 29a + d1sigma0 = lib.einsum('r,rij,rkl,rmn->rijklmn',weights,dma1r,dma2r,dma3r) + casdm3s[0][i:j, i:j, k:l, k:l, m:n, m:n] += np.tensordot(weights, d1sigma0, axes=1) + _transpose_dm3_(casdm3s[0],i,j, i,j, k,l, k,l, m,n, m,n) + + d1sigma11 = lib.einsum('r,rij,rkl,rmn->rijklmn',weights,dma1r,dma2r,dmb3r) + d1sigma12 = lib.einsum('r,rij,rkl,rmn->rijmnkl',weights,dma1r,dmb2r,dma3r) + d1sigma13 = lib.einsum('r,rij,rkl,rmn->rmnijkl',weights,dma1r,dmb2r,dma3r) + d1sigma14 = lib.einsum('r,rij,rkl,rmn->rmnklij',weights,dmb1r,dma2r,dma3r) + d1sigma15 = lib.einsum('r,rij,rkl,rmn->rklmnij',weights,dmb1r,dma2r,dma3r) + d1sigma16 = lib.einsum('r,rij,rkl,rmn->rklijmn',weights,dma1r,dma2r,dmb3r) + + d1sigma21 = lib.einsum('r,rij,rkl,rmn->rijklmn',weights,dma1r,dmb2r,dmb3r) + d1sigma22 = lib.einsum('r,rij,rkl,rmn->rijmnkl',weights,dma1r,dmb2r,dmb3r) + d1sigma23 = lib.einsum('r,rij,rkl,rmn->rmnijkl',weights,dmb1r,dmb2r,dma3r) + d1sigma24 = lib.einsum('r,rij,rkl,rmn->rmnklij',weights,dmb1r,dmb2r,dma3r) + d1sigma25 = lib.einsum('r,rij,rkl,rmn->rklmnij',weights,dmb1r,dma2r,dmb3r) + d1sigma26 = lib.einsum('r,rij,rkl,rmn->rklijmn',weights,dmb1r,dma2r,dmb3r) + + d1sigma3 = lib.einsum('r,rij,rkl,rmn->rijklmn',weights,dmb1r,dmb2r,dmb3r) + casdm3s[3][i:j, i:j, k:l, k:l, m:n, m:n] += np.tensordot(weights, d1sigma3, axes=1) + _transpose_dm3_(casdm3s[3],i,j, i,j, k,l, k,l, m,n, m,n) + + casdm3s[1][i:j, i:j, k:l, k:l, m:n, m:n] += np.tensordot(weights, d1sigma11, axes=1) + casdm3s[1][i:j, i:j, m:n, m:n, k:l, k:l] += np.tensordot(weights, d1sigma12, axes=1) + casdm3s[1][m:n, m:n, i:j, i:j, k:l, k:l] += np.tensordot(weights, d1sigma13, axes=1) + casdm3s[1][m:n, m:n, k:l, k:l, i:j, i:j] += np.tensordot(weights, d1sigma14, axes=1) + casdm3s[1][k:l, k:l, m:n, m:n, i:j, i:j] += np.tensordot(weights, d1sigma15, axes=1) + casdm3s[1][k:l, k:l, i:j, i:j, m:n, m:n] += np.tensordot(weights, d1sigma16, axes=1) + + casdm3s[2][i:j, i:j, k:l, k:l, m:n, m:n] += np.tensordot(weights, d1sigma21, axes=1) + casdm3s[2][i:j, i:j, m:n, m:n, k:l, k:l] += np.tensordot(weights, d1sigma22, axes=1) + casdm3s[2][m:n, m:n, i:j, i:j, k:l, k:l] += np.tensordot(weights, d1sigma23, axes=1) + casdm3s[2][m:n, m:n, k:l, k:l, i:j, i:j] += np.tensordot(weights, d1sigma24, axes=1) + casdm3s[2][k:l, k:l, m:n, m:n, i:j, i:j] += np.tensordot(weights, d1sigma25, axes=1) + casdm3s[2][k:l, k:l, i:j, i:j, m:n, m:n] += np.tensordot(weights, d1sigma26, axes=1) + + # Term 2- 29b + d2sigma0 = lib.einsum('r,rij,rkl,rmn->rijknml',weights,dma1r,dma2r,dma3r) + casdm3s[0][i:j, i:j, k:l, m:n, m:n, k:l] -= np.tensordot(weights, d2sigma0, axes=1) + _transpose_dm3_(casdm3s[0], i,j, i,j, k,l, m,n, m,n, k,l) + + d2sigma14 = lib.einsum('r,rij,rkl,rmn->rmlknij',weights,dmb1r,dma2r,dma3r) + d2sigma15 = lib.einsum('r,rij,rkl,rmn->rknmlij',weights,dmb1r,dma2r,dma3r) + casdm3s[1][m:n, k:l, k:l, m:n, i:j, i:j] -= np.tensordot(weights, d2sigma14, axes=1) + casdm3s[1][k:l, m:n, m:n, k:l, i:j, i:j] -= np.tensordot(weights, d2sigma15, axes=1) + + d2sigma21 = lib.einsum('r,rij,rkl,rmn->rijknml',weights,dma1r,dmb2r,dmb3r) + d2sigma22 = lib.einsum('r,rij,rkl,rmn->rijmlkn',weights,dma1r,dmb2r,dmb3r) + casdm3s[2][i:j, i:j, k:l, m:n, m:n, k:l] -= np.tensordot(weights, d2sigma21, axes=1) + casdm3s[2][i:j, i:j, m:n, k:l, k:l, m:n] -= np.tensordot(weights, d2sigma22, axes=1) + + d2sigma3 = lib.einsum('r,rij,rkl,rmn->rijknml',weights,dmb1r,dmb2r,dmb3r) + casdm3s[3][i:j, i:j, k:l, m:n, m:n, k:l] -= np.tensordot(weights, d2sigma3, axes=1) + _transpose_dm3_(casdm3s[3], i,j, i,j, k,l, m,n, m,n, k,l) + + # Term 3- 29c + d3sigma0 = lib.einsum('r,rij,rkl,rmn->rilkjmn',weights,dma1r,dma2r,dma3r) + casdm3s[0][i:j, k:l, k:l, i:j, m:n, m:n] -= np.tensordot (weights, d3sigma0, axes=1) + _transpose_dm3_(casdm3s[0], i,j, k,l, k,l, i,j, m,n, m,n) + + d3sigma11 = lib.einsum('r,rij,rkl,rmn->rilkjmn',weights,dma1r,dma2r,dmb3r) + d3sigma16 = lib.einsum('r,rij,rkl,rmn->rkjilmn',weights,dma1r,dma2r,dmb3r) + casdm3s[1][i:j, k:l, k:l, i:j, m:n, m:n] -= np.tensordot (weights, d3sigma11, axes=1) + casdm3s[1][k:l, i:j, i:j, k:l, m:n, m:n] -= np.tensordot (weights, d3sigma16, axes=1) + + d3sigma23 = lib.einsum('r,rij,rkl,rmn->rmnilkj',weights,dmb1r,dmb2r,dma3r) + d3sigma24 = lib.einsum('r,rij,rkl,rmn->rmnkjil',weights,dmb1r,dmb2r,dma3r) + casdm3s[2][m:n, m:n, i:j, k:l, k:l, i:j] -= np.tensordot (weights, d3sigma23, axes=1) + casdm3s[2][m:n, m:n, k:l, i:j, i:j, k:l] -= np.tensordot (weights, d3sigma24, axes=1) + + d3sigma3 = lib.einsum('r,rij,rkl,rmn->rilkjmn',weights,dmb1r,dmb2r,dmb3r) + casdm3s[3][i:j, k:l, k:l, i:j, m:n, m:n] -= np.tensordot (weights, d3sigma3, axes=1) + _transpose_dm3_(casdm3s[3], i,j, k,l, k,l, i,j, m,n, m,n) + + # Term 4- 29d + d4sigma0 = lib.einsum('r,rij,rkl,rmn->rilknmj',weights,dma1r,dma2r,dma3r) + casdm3s[0][i:j, k:l, k:l, m:n, m:n, i:j] += np.tensordot(weights, d4sigma0, axes=1) + _transpose_dm3_(casdm3s[0], i,j, k,l, k,l, m,n, m,n, i,j) + + d4sigma3 = lib.einsum('r,rij,rkl,rmn->rilknmj',weights,dmb1r,dmb2r,dmb3r) + casdm3s[3][i:j, k:l, k:l, m:n, m:n, i:j] += np.tensordot(weights, d4sigma3, axes=1) + _transpose_dm3_(casdm3s[3], i,j, k,l, k,l, m,n, m,n, i,j) + + # Term 5- 29e + d5sigma0 = lib.einsum('r,rij,rkl,rmn->rinklmj',weights,dma1r,dma2r,dma3r) + casdm3s[0][i:j, m:n, k:l, k:l, m:n, i:j] -= np.tensordot(weights, d5sigma0, axes=1) + _transpose_dm3_(casdm3s[0], i,j, m,n, k,l, k,l, m,n, i,j) + + d5sigma12 = lib.einsum('r,rij,rkl,rmn->rinmjkl',weights,dma1r,dma2r,dmb3r) + d5sigma13 = lib.einsum('r,rij,rkl,rmn->rmjinkl',weights,dma1r,dma2r,dmb3r) + casdm3s[1][i:j, m:n, m:n, i:j, k:l, k:l] -= np.tensordot(weights, d5sigma12, axes=1) + casdm3s[1][m:n, i:j, i:j, m:n, k:l, k:l] -= np.tensordot(weights, d5sigma13, axes=1) + + d5sigma25 = lib.einsum('r,rij,rkl,rmn->rklmjin',weights,dmb1r,dma2r,dmb3r) + d5sigma26 = lib.einsum('r,rij,rkl,rmn->rklinmj',weights,dmb1r,dma2r,dmb3r) + casdm3s[2][k:l, k:l, m:n, i:j, i:j, m:n] -= np.tensordot(weights, d5sigma25, axes=1) + casdm3s[2][k:l, k:l, i:j, m:n, m:n, i:j] -= np.tensordot(weights, d5sigma26, axes=1) + + d5sigma3 = lib.einsum('r,rij,rkl,rmn->rinklmj',weights,dmb1r,dmb2r,dmb3r) + casdm3s[3][i:j, m:n, k:l, k:l, m:n, i:j] -= np.tensordot(weights, d5sigma3, axes=1) + _transpose_dm3_(casdm3s[3], i,j, m,n, k,l, k,l, m,n, i,j) + + # Term 6- 29f + d6sigma0 = lib.einsum('r,rij,rkl,rmn->rinkjml',weights,dma1r,dma2r,dma3r) + casdm3s[0][i:j, m:n, k:l, i:j, m:n, k:l] += np.tensordot(weights, d6sigma0, axes=1) + _transpose_dm3_(casdm3s[0], i,j, m,n, k,l, i,j, m,n, k,l) + + d6sigma3 = lib.einsum('r,rij,rkl,rmn->rinkjml',weights,dmb1r,dmb2r,dmb3r) + casdm3s[3][i:j, m:n, k:l, i:j, m:n, k:l] += np.tensordot(weights, d6sigma3, axes=1) + _transpose_dm3_(casdm3s[3], i,j, m,n, k,l, i,j, m,n, k,l) + + + #Last 2 terms - combs of f1+f1+f2: dm1rs1=f1, dm2rs2=f2 + + for isub1, isub2 in permutations (range (self.nfrags), 2): + dm2rs2 = casdm2frs[isub1] + dm1rs1 = casdm1frs[isub2] + + i = ncas_cum[isub1] + j = ncas_cum[isub1+1] + k = ncas_cum[isub2] + l = ncas_cum[isub2+1] + + dma1r, dmb1r = dm1rs1[:,0], dm1rs1[:,1] + dmaa2r, dmab2r, dmbb2r = dm2rs2[:,0], dm2rs2[:,1], dm2rs2[:,2] + + # Term 7 - 29g + d4sigma0 = lib.einsum('r,rmn,rijkl->rijklmn',weights,dma1r,dmaa2r) + + d4sigma11 = lib.einsum('r,rmn,rijkl->rijklmn',weights,dmb1r,dmaa2r) + d4sigma12 = lib.einsum('r,rmn,rijkl->rijmnkl',weights,dma1r,dmab2r) + d4sigma13 = lib.einsum('r,rmn,rijkl->rmnijkl',weights,dma1r,dmab2r) + d4sigma14 = lib.einsum('r,rmn,rklij->rmnklij',weights,dma1r,dmab2r) + d4sigma15 = lib.einsum('r,rmn,rklij->rklmnij',weights,dma1r,dmab2r) + d4sigma16 = lib.einsum('r,rmn,rklij->rklijmn',weights,dmb1r,dmaa2r) + + d4sigma21 = lib.einsum('r,rmn,rijkl->rijklmn',weights,dmb1r,dmab2r) + d4sigma22 = lib.einsum('r,rmn,rijkl->rijmnkl',weights,dmb1r,dmab2r) + d4sigma23 = lib.einsum('r,rmn,rijkl->rmnijkl',weights,dma1r,dmbb2r) + d4sigma24 = lib.einsum('r,rmn,rklij->rmnklij',weights,dma1r,dmbb2r) + d4sigma25 = lib.einsum('r,rmn,rklij->rklmnij',weights,dmb1r,dmab2r) + d4sigma26 = lib.einsum('r,rmn,rklij->rklijmn',weights,dmb1r,dmab2r) + + d4sigma3 = lib.einsum('r,rmn,rijkl->rijklmn',weights,dmb1r,dmbb2r) + + casdm3s[0][i:j, i:j, i:j, i:j, k:l, k:l] += np.tensordot (weights, d4sigma0, axes=1) + d0 = casdm3s[0][i:j, i:j, i:j, i:j, k:l, k:l] + casdm3s[0][i:j, i:j, k:l, k:l, i:j, i:j] = d0.transpose(0,1,4,5,2,3) + casdm3s[0][k:l, k:l, i:j, i:j, i:j, i:j] = d0.transpose(4,5,0,1,2,3) + + casdm3s[1][i:j, i:j, i:j, i:j, k:l, k:l] += np.tensordot (weights, d4sigma11, axes=1) + casdm3s[1][i:j, i:j, k:l, k:l, i:j, i:j] += np.tensordot (weights, d4sigma12, axes=1) + casdm3s[1][k:l, k:l, i:j, i:j, i:j, i:j] += np.tensordot (weights, d4sigma13, axes=1) + + casdm3s[2][i:j, i:j, i:j, i:j, k:l, k:l] += np.tensordot (weights, d4sigma21, axes=1) + casdm3s[2][i:j, i:j, k:l, k:l, i:j, i:j] += np.tensordot (weights, d4sigma22, axes=1) + casdm3s[2][k:l, k:l, i:j, i:j, i:j, i:j] += np.tensordot (weights, d4sigma23, axes=1) + + casdm3s[3][i:j, i:j, i:j, i:j, k:l, k:l] += np.tensordot (weights, d4sigma3, axes=1) + d3 = casdm3s[3][i:j, i:j, i:j, i:j, k:l, k:l] + casdm3s[3][i:j, i:j, k:l, k:l, i:j, i:j] = d3.transpose(0,1,4,5,2,3) + casdm3s[3][k:l, k:l, i:j, i:j, i:j, i:j] = d3.transpose(4,5,0,1,2,3) + + # Term 8 - 29h + d5sigma0 = lib.einsum('r,rmn,rijkl->rijknml',weights,dma1r,dmaa2r) + + d5sigma14 = lib.einsum('r,rmn,rklij->rmlknij',weights,dma1r,dmab2r) + d5sigma15 = lib.einsum('r,rmn,rklij->rknmlij',weights,dma1r,dmab2r) + + d5sigma21 = lib.einsum('r,rmn,rijkl->rijknml',weights,dmb1r,dmab2r) + d5sigma22 = lib.einsum('r,rmn,rijkl->rijmlkn',weights,dmb1r,dmab2r) + + d5sigma3 = lib.einsum('r,rmn,rijkl->rijknml',weights,dmb1r,dmbb2r) + + casdm3s[0][i:j, i:j, i:j, k:l, k:l, i:j] -= np.tensordot (weights, d5sigma0, axes=1) + _transpose_dm3_(casdm3s[0],i,j, i,j, i,j, k,l, k,l, i,j) + + casdm3s[1][k:l, i:j, i:j, k:l, i:j, i:j] -= np.tensordot (weights, d5sigma14, axes=1) + casdm3s[1][i:j, k:l, k:l, i:j, i:j, i:j] -= np.tensordot (weights, d5sigma15, axes=1) + + casdm3s[2][i:j, i:j, i:j, k:l, k:l, i:j] -= np.tensordot (weights, d5sigma21, axes=1) + casdm3s[2][i:j, i:j, k:l, i:j, i:j, k:l] -= np.tensordot (weights, d5sigma22, axes=1) + + casdm3s[3][i:j, i:j, i:j, k:l, k:l, i:j] -= np.tensordot (weights, d5sigma3, axes=1) + _transpose_dm3_(casdm3s[3],i,j, i,j, i,j, k,l, k,l, i,j) + + return casdm3s + def get_veff (self, mol=None, dm1s=None, hermi=1, spin_sep=False, **kwargs): ''' Returns a spin-summed veff! If dm1s isn't provided, builds from self.mo_coeff, self.ci etc. ''' diff --git a/tests/lasscf/test_lasscf_rdm123s.py b/tests/lasscf/test_lasscf_rdm123s.py new file mode 100644 index 00000000..0b6f6f04 --- /dev/null +++ b/tests/lasscf/test_lasscf_rdm123s.py @@ -0,0 +1,242 @@ +#!/usr/bin/env python +# +# Author: Shreya Verma + +# The tests are performed for the following purposes: +# 1. Check accuracy of spin-summed LASSCF 2-RDMs against spin-summed 2-RDMs created by LAS ci vectors using the PySCF rdm generator- only for H8 with 4 fragments +# 2. Check accuracy of spin-separated LASSCF 2-RDMs against spin-separated 2-RDMs created by LAS ci vectors using the PySCF rdm generator- only for H8 with 4 fragments +# 3. Check accuracy of spin-summed LASSCF 3-RDMs against spin-summed 3-RDMs created by LAS ci vectors using the PySCF rdm generator- only for H8 with 4 fragments +# 4. Check accuracy of spin-separated LASSCF 3-RDMs against spin-separated 3-RDMs created by LAS ci vectors using the PySCF rdm generator- only for H8 with 4 fragments +# 5. Check accuracy of spin-separated LASSCF 3-RDMs by tracing them down to 2-RDMs +# 6. Check all of the above for different sub-spins of LASSCF fragments + +import unittest + +from pyscf import scf, gto, df, lib, fci +from pyscf import mcscf +from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF +from mrh.my_pyscf.mcscf.addons import las2cas_civec +import numpy as np + +def generate_geom (natom = 8, dist = 1.0): + """Generate H chain of natoms each with dist distance linearly""" + coords = [] + for i in range(natom): + coords.append(f"H 0.0 0.0 {i*dist}") + return "\n".join(coords) + +def generate_ci_rdm123s (ci,nelec): + nelec_sum = nelec[0][0]+nelec[0][1] + norb = nelec_sum + rdm1_fci, rdm2_fci, rdm3_fci = fci.rdm.make_dm123 ('FCI3pdm_kern_sf',ci,ci, norb, nelec_sum) + rdm1_no, rdm2_no, rdm3_no = fci.rdm.reorder_dm123 (rdm1_fci, rdm2_fci, rdm3_fci) #normal ordered RDMs + + rdm1s, rdm2s, rdm3s = fci.direct_spin1.make_rdm123s(np.array(ci), norb, nelec[0]) + + return rdm1_no, rdm2_no, rdm3_no, rdm1s, rdm2s, rdm3s + +def mult_frags(nelesub, norbsub, charge=None, spin_sub=None, frag_atom_list=None, density_fit=False): + """Used for checking systems with more than three LAS fragments to see if the 3-rdms are generated properly, here for H8 with 4 fragments.""" + xyz = generate_geom (natom = 8, dist = 1.0) + mol = gto.M(atom=xyz, basis='sto-3g', charge=charge, verbose=0) + mf = scf.RHF(mol).run() + + if spin_sub is None: + spin_sub = (1,1,1,1) + else: + spin_sub = spin_sub + + if frag_atom_list is None: frag_atom_list = ((0, 1), (2, 3), (4, 5), (6,7)) + + las = LASSCF (mf, nelesub, norbsub, spin_sub=spin_sub) + mo_loc = las.localize_init_guess (frag_atom_list, mf.mo_coeff) + las.kernel(mo_loc) + + rdm1_las = las.make_casdm1() + rdm1s_las = las.make_casdm1s() + rdm2_las = las.make_casdm2() + rdm2s_las = las.make_casdm2s() + rdm3_las = las.make_casdm3() + rdm3s_las = las.make_casdm3s() + + ci, nelec = las2cas_civec(las) + rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s = generate_ci_rdm123s(ci,nelec) # normal ordered RDMs + + return rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s + +class KnownValues(unittest.TestCase): + + def trace_down_test(self, rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s): + """ Test for tracing down spin-sep 3-RDMs down to 2-RDMs""" + nelec = 8 + fac = 1. / (nelec//2-1) + rdm2s_1a1 = np.einsum('ijkk->ij',rdm2s_las[0])*fac + rdm2s_1a2 = np.einsum('kkij->ij',rdm2s_las[0])*fac + rdm2s_1a3 = np.einsum('ijkk->ij',rdm2s_las[1])*(1.0/(nelec//2)) + rdm2s_1b1 = np.einsum('ijkk->ij',rdm2s_las[2])*fac + rdm2s_1b2 = np.einsum('kkij->ij',rdm2s_las[2])*fac + rdm2s_1b3 = np.einsum('ijkk->ij',rdm2s_las[1])*(1.0/(nelec//2)) + + self.assertAlmostEqual (lib.fp (rdm2s_1a1), lib.fp (rdm1s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1a2), lib.fp (rdm1s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1a3), lib.fp (rdm1s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1b1), lib.fp (rdm1s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1b2), lib.fp (rdm1s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1b3), lib.fp (rdm1s_las[1]), 12) + + # Tracing down ss 3-RDM to 2-RDM + nelec = 8 + + rdm3s_2a1 = np.einsum('ijklmm->ijkl',rdm3s_las[0])*(1.0/(nelec//2-2)) + rdm3s_2a2 = np.einsum('ijmmkl->ijkl',rdm3s_las[0])*(1.0/(nelec//2-2)) + rdm3s_2a3 = np.einsum('mmijkl->ijkl',rdm3s_las[0])*(1.0/(nelec//2-2)) + rdm3s_2a4 = np.einsum('ijklmm->ijkl',rdm3s_las[1])*(1.0/(nelec//2)) + + rdm3s_2b1 = np.einsum('ijklmm->ijkl',rdm3s_las[3])*(1.0/(nelec//2-2)) + rdm3s_2b2 = np.einsum('ijmmkl->ijkl',rdm3s_las[3])*(1.0/(nelec//2-2)) + rdm3s_2b3 = np.einsum('mmijkl->ijkl',rdm3s_las[3])*(1.0/(nelec//2-2)) + rdm3s_2b4 = np.einsum('mmijkl->ijkl',rdm3s_las[2])*(1.0/(nelec//2)) + + rdm3s_2ab1 = np.einsum('ijmmkl->ijkl',rdm3s_las[1])*(1.0/(nelec//2-1)) + rdm3s_2ab2 = np.einsum('mmijkl->ijkl',rdm3s_las[1])*(1.0/(nelec//2-1)) + rdm3s_2ab3 = np.einsum('ijklmm->ijkl',rdm3s_las[2])*(1.0/(nelec//2-1)) + rdm3s_2ab4 = np.einsum('ijmmkl->ijkl',rdm3s_las[2])*(1.0/(nelec//2-1)) + + self.assertAlmostEqual (lib.fp (rdm3s_2a1), lib.fp (rdm2s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2a2), lib.fp (rdm2s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2a3), lib.fp (rdm2s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2a4), lib.fp (rdm2s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2b1), lib.fp (rdm2s_las[2]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2b2), lib.fp (rdm2s_las[2]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2b3), lib.fp (rdm2s_las[2]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2b4), lib.fp (rdm2s_las[2]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2ab1), lib.fp (rdm2s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2ab2), lib.fp (rdm2s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2ab3), lib.fp (rdm2s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2ab4), lib.fp (rdm2s_las[1]), 12) + + def test_rdm3_h8_sto3g_spin1(self): + """Spin cas 1""" + spin_sub = (1,1,1,1) + frag_atom_list = ((0, 1), (2, 3), (4, 5), (6,7)) + nelesub = (2,2,2,2) + norbsub = (2,2,2,2) + rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s = mult_frags(nelesub, norbsub, charge=None, spin_sub=spin_sub, density_fit=False) + + # Direct tests with ci generated RDMs + self.assertAlmostEqual (lib.fp (rdm1_las), lib.fp (rdm1), 12) + self.assertAlmostEqual (lib.fp (rdm1s_las), lib.fp (rdm1s), 12) + self.assertAlmostEqual (lib.fp (rdm2_las), lib.fp (rdm2), 12) + self.assertAlmostEqual (lib.fp (rdm2s_las), lib.fp (rdm2s), 12) + self.assertAlmostEqual (lib.fp (rdm3_las), lib.fp (rdm3), 12) + self.assertAlmostEqual (lib.fp (rdm3s_las), lib.fp (rdm3s), 12) + + # Tracing down tests + nelec = 8 + fac = 1. / (nelec//2-1) + rdm2s_1a1 = np.einsum('ijkk->ij',rdm2s_las[0])*fac + rdm2s_1a2 = np.einsum('kkij->ij',rdm2s_las[0])*fac + rdm2s_1a3 = np.einsum('ijkk->ij',rdm2s_las[1])*(1.0/(nelec//2)) + rdm2s_1b1 = np.einsum('ijkk->ij',rdm2s_las[2])*fac + rdm2s_1b2 = np.einsum('kkij->ij',rdm2s_las[2])*fac + rdm2s_1b3 = np.einsum('ijkk->ij',rdm2s_las[1])*(1.0/(nelec//2)) + + self.assertAlmostEqual (lib.fp (rdm2s_1a1), lib.fp (rdm1s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1a2), lib.fp (rdm1s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1a3), lib.fp (rdm1s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1b1), lib.fp (rdm1s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1b2), lib.fp (rdm1s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm2s_1b3), lib.fp (rdm1s_las[1]), 12) + + # Tracing down ss 3-RDM to 2-RDM + nelec = 8 + + rdm3s_2a1 = np.einsum('ijklmm->ijkl',rdm3s_las[0])*(1.0/(nelec//2-2)) + rdm3s_2a2 = np.einsum('ijmmkl->ijkl',rdm3s_las[0])*(1.0/(nelec//2-2)) + rdm3s_2a3 = np.einsum('mmijkl->ijkl',rdm3s_las[0])*(1.0/(nelec//2-2)) + rdm3s_2a4 = np.einsum('ijklmm->ijkl',rdm3s_las[1])*(1.0/(nelec//2)) + + rdm3s_2b1 = np.einsum('ijklmm->ijkl',rdm3s_las[3])*(1.0/(nelec//2-2)) + rdm3s_2b2 = np.einsum('ijmmkl->ijkl',rdm3s_las[3])*(1.0/(nelec//2-2)) + rdm3s_2b3 = np.einsum('mmijkl->ijkl',rdm3s_las[3])*(1.0/(nelec//2-2)) + rdm3s_2b4 = np.einsum('mmijkl->ijkl',rdm3s_las[2])*(1.0/(nelec//2)) + + rdm3s_2ab1 = np.einsum('ijmmkl->ijkl',rdm3s_las[1])*(1.0/(nelec//2-1)) + rdm3s_2ab2 = np.einsum('mmijkl->ijkl',rdm3s_las[1])*(1.0/(nelec//2-1)) + rdm3s_2ab3 = np.einsum('ijklmm->ijkl',rdm3s_las[2])*(1.0/(nelec//2-1)) + rdm3s_2ab4 = np.einsum('ijmmkl->ijkl',rdm3s_las[2])*(1.0/(nelec//2-1)) + + self.assertAlmostEqual (lib.fp (rdm3s_2a1), lib.fp (rdm2s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2a2), lib.fp (rdm2s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2a3), lib.fp (rdm2s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2a4), lib.fp (rdm2s_las[0]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2b1), lib.fp (rdm2s_las[2]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2b2), lib.fp (rdm2s_las[2]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2b3), lib.fp (rdm2s_las[2]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2b4), lib.fp (rdm2s_las[2]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2ab1), lib.fp (rdm2s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2ab2), lib.fp (rdm2s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2ab3), lib.fp (rdm2s_las[1]), 12) + self.assertAlmostEqual (lib.fp (rdm3s_2ab4), lib.fp (rdm2s_las[1]), 12) + + def test_rdm3_h8_sto3g_spin2(self): + """Spin cas 2""" + spin_sub = (1,3,1,3) + frag_atom_list = ((0, 1), (2, 3), (4, 5), (6,7)) + nelesub = (2,2,2,2) + norbsub = (2,2,2,2) + + rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s = mult_frags(nelesub, norbsub, charge=None, spin_sub=spin_sub, density_fit=False) + + # Direct tests with ci generated RDMs + self.assertAlmostEqual (lib.fp (rdm1_las), lib.fp (rdm1), 12) + self.assertAlmostEqual (lib.fp (rdm1s_las), lib.fp (rdm1s), 12) + self.assertAlmostEqual (lib.fp (rdm2_las), lib.fp (rdm2), 12) + self.assertAlmostEqual (lib.fp (rdm2s_las), lib.fp (rdm2s), 12) + self.assertAlmostEqual (lib.fp (rdm3_las), lib.fp (rdm3), 12) + self.assertAlmostEqual (lib.fp (rdm3s_las), lib.fp (rdm3s), 12) + # Tracing down tests + self.trace_down_test(rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s) + + + def test_rdm3_h8_sto3g_spin3(self): + """Spin cas 3""" + spin_sub = (1,3,3,3) + frag_atom_list = ((0, 1), (2, 3), (4, 5), (6,7)) + nelesub = (2,2,2,2) + norbsub = (2,2,2,2) + + rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s = mult_frags(nelesub, norbsub, charge=None, spin_sub=spin_sub, density_fit=False) + + # Direct tests with ci generated RDMs + self.assertAlmostEqual (lib.fp (rdm1_las), lib.fp (rdm1), 12) + self.assertAlmostEqual (lib.fp (rdm1s_las), lib.fp (rdm1s), 12) + self.assertAlmostEqual (lib.fp (rdm2_las), lib.fp (rdm2), 12) + self.assertAlmostEqual (lib.fp (rdm2s_las), lib.fp (rdm2s), 12) + self.assertAlmostEqual (lib.fp (rdm3_las), lib.fp (rdm3), 12) + self.assertAlmostEqual (lib.fp (rdm3s_las), lib.fp (rdm3s), 12) + # Tracing down tests + self.trace_down_test(rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s) + + def test_rdm3_h8_sto3g_spin4(self): + """Spin cas 4""" + spin_sub = (3,3,3,3) + frag_atom_list = ((0, 1), (2, 3), (4, 5), (6,7)) + nelesub = (2,2,2,2) + norbsub = (2,2,2,2) + + rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s = mult_frags(nelesub, norbsub, charge=None, spin_sub=spin_sub, density_fit=False) + + # Direct tests with ci generated RDMs + self.assertAlmostEqual (lib.fp (rdm1_las), lib.fp (rdm1), 12) + self.assertAlmostEqual (lib.fp (rdm1s_las), lib.fp (rdm1s), 12) + self.assertAlmostEqual (lib.fp (rdm2_las), lib.fp (rdm2), 12) + self.assertAlmostEqual (lib.fp (rdm2s_las), lib.fp (rdm2s), 12) + self.assertAlmostEqual (lib.fp (rdm3_las), lib.fp (rdm3), 12) + self.assertAlmostEqual (lib.fp (rdm3s_las), lib.fp (rdm3s), 12) + # Tracing down tests + self.trace_down_test(rdm1_las, rdm2_las, rdm3_las, rdm1s_las, rdm2s_las, rdm3s_las, rdm1, rdm2, rdm3, rdm1s, rdm2s, rdm3s) + +if __name__ == "__main__": + print("Full Tests for LASSCF 3-RDMs") + unittest.main()