Skip to content

Add the additional triples term when the field is on. #80

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
zw0930 wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add the additional triples term when the field is on. #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
108 changes: 92 additions & 16 deletions pycc/ccdensity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
import time
import numpy as np
import torch
from .cctriples import t3c_ijk, t3c_abc, l3_ijk, l3_abc, t3c_bc, l3_bc, t3_pert_ijk, t3_pert_bc
from .cctriples import t3c_ijk, t3c_abc, l3_ijk, l3_abc, t3c_bc, l3_bc, t3_pert_ijk, t3_pert_bc, t3_ijkabc_alt

class ccdensity(object):
"""
Expand Down Expand Up @@ -172,6 +172,8 @@ def compute_onepdm(self, t1, t2, l1, l2, real_time=False):
F = self.ccwfn.H.F
ERI = self.ccwfn.H.ERI
L = self.ccwfn.H.L

contract = self.contract

if isinstance(t1, torch.Tensor):
if self.ccwfn.precision == 'DP':
Expand All @@ -195,16 +197,37 @@ def compute_onepdm(self, t1, t2, l1, l2, real_time=False):
Wovoo = self.ccwfn.build_cc3_Wmbij(o, v, ERI, t1, Woooo)
Wvovv = self.ccwfn.build_cc3_Wamef(o, v, ERI, t1)
Wooov = self.ccwfn.build_cc3_Wmnie(o, v, ERI, t1)

if real_time is True:
t3_full = t3_ijkabc_alt(o, v, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
"""
for i in range(no):
for j in range(no):
for k in range(no):
for a in range(nv):
for b in range(nv):
for c in range(nv):
t3_full[i,j,k,a,b,c] = t3_ijkabc(o, v, i, j, k, a, b, c, t2, Wvvvo, Wovoo, F, contract)
"""
opdm[o,v] += self.build_cc3_Dov_pert(o, v, no, nv, F, L, t1, t2, l1, l2, t3_full, Wvvvo, Wovoo, Fov, Wvovv, Wooov, real_time=real_time)
# Density matrix blocks in contractions with T1-transformed dipole integrals
if isinstance(t1, torch.Tensor):
opdm_cc3 = torch.zeros_like(opdm)
else:
opdm_cc3 = np.zeros_like(opdm)
opdm_cc3[o,o] += self.build_cc3_Doo_pert(o, v, no, nv, F, L, t2, l1, l2, t3_full, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=real_time)
opdm_cc3[v,v] += self.build_cc3_Dvv_pert(o, v, no, nv, F, L, t2, l1, l2, t3_full, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=real_time)

opdm[o,v] += self.build_cc3_Dov(o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov, real_time=real_time)

# Density matrix blocks in contractions with T1-transformed dipole integrals
if isinstance(t1, torch.Tensor):
opdm_cc3 = torch.zeros_like(opdm)
else:
opdm_cc3 = np.zeros_like(opdm)
opdm_cc3[o,o] += self.build_cc3_Doo(o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov)
opdm_cc3[v,v] += self.build_cc3_Dvv(o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov)
opdm[o,v] += self.build_cc3_Dov(o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov)

# Density matrix blocks in contractions with T1-transformed dipole integrals
if isinstance(t1, torch.Tensor):
opdm_cc3 = torch.zeros_like(opdm)
else:
opdm_cc3 = np.zeros_like(opdm)
opdm_cc3[o,o] += self.build_cc3_Doo(o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov)
opdm_cc3[v,v] += self.build_cc3_Dvv(o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov)

return (opdm, opdm_cc3)

Expand Down Expand Up @@ -275,6 +298,60 @@ def build_Dov(self, t1, t2, l1, l2): # complete

# CC3 contributions to the one electron densities
def build_cc3_Dov(self, o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov, real_time=False):
contract = self.contract
if isinstance(t1, torch.Tensor):
Dov = torch.zeros_like(t1)
Zlmdi = torch.zeros_like(t2[:,:,:,:no])
else:
Dov = np.zeros_like(t1)
Zlmdi = np.zeros_like(t2[:,:,:,:no])
for i in range(no):
for j in range(no):
for k in range(no):
l3 = l3_ijk(i, j, k, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract)
# Intermediate for Dov_2
Zlmdi[i,j] += contract('def,ife->di', l3, t2[k])
# Dov_1
t3 = t3c_ijk(o, v, i, j, k, t2, Wvvvo, Wovoo, F, contract)
Dov[i] += contract('abc,bc->a', t3 - t3.swapaxes(0,1), l2[j,k])
# Dov_2
Dov -= contract('lmdi, lmda->ia', Zlmdi, t2)

return Dov

def build_cc3_Doo(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=False):
contract = self.contract
if isinstance(l1, torch.Tensor):
Doo = torch.zeros_like(l1[:,:no])
else:
Doo = np.zeros_like(l1[:,:no])
for b in range(nv):
for c in range(nv):
t3 = t3c_bc(o, v, b, c, t2, Wvvvo, Wovoo, F, contract)
l3 = l3_bc(b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract)
Doo -= 0.5 * contract('lmia,lmja->ij', t3, l3)

return Doo

def build_cc3_Dvv(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=False):
contract = self.contract
if isinstance(l1, torch.Tensor):
Dvv = torch.zeros_like(l1)
Dvv = torch.nn.functional.pad(Dvv, (0,0,0,nv-no))
else:
Dvv = np.zeros_like(l1)
Dvv = np.pad(Dvv, ((0,nv-no), (0,0)))
for i in range(no):
for j in range(no):
for k in range(no):
t3 = t3c_ijk(o, v, i, j, k, t2, Wvvvo, Wovoo, F, contract)
l3 = l3_ijk(i, j, k, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract)
Dvv += 0.5 * contract('bdc,adc->ab', t3, l3)

return Dvv

# CC3 contributions to the one electron densities when a perturbation is present
def build_cc3_Dov_pert(self, o, v, no, nv, F, L, t1, t2, l1, l2, t3_full, Wvvvo, Wovoo, Fov, Wvovv, Wooov, real_time=False):
contract = self.contract
if isinstance(t1, torch.Tensor):
Dov = torch.zeros_like(t1)
Expand All @@ -295,14 +372,14 @@ def build_cc3_Dov(self, o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, W
V = F - self.ccwfn.H.F.clone()
else:
V = F - self.ccwfn.H.F.copy()
t3 -= t3_pert_ijk(o, v, i, j, k, t2, V, F, contract)
t3 -= t3_pert_ijk(o, v, i, j, k, t2, t3_full, V, F, contract)
Dov[i] += contract('abc,bc->a', t3 - t3.swapaxes(0,1), l2[j,k])
# Dov_2
Dov -= contract('lmdi, lmda->ia', Zlmdi, t2)

return Dov

def build_cc3_Doo(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=False):
def build_cc3_Doo_pert(self, o, v, no, nv, F, L, t2, l1, l2, t3_full, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=False):
contract = self.contract
if isinstance(l1, torch.Tensor):
Doo = torch.zeros_like(l1[:,:no])
Expand All @@ -316,13 +393,13 @@ def build_cc3_Doo(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv
V = F - self.ccwfn.H.F.clone()
else:
V = F - self.ccwfn.H.F.copy()
t3 -= t3_pert_bc(o, v, b, c, t2, V, F, contract)
t3 -= t3_pert_bc(o, v, b, c, t2, t3_full, V, F, contract)
l3 = l3_bc(b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract)
Doo -= 0.5 * contract('lmia,lmja->ij', t3, l3)

return Doo

def build_cc3_Dvv(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=False):
def build_cc3_Dvv_pert(self, o, v, no, nv, F, L, t2, l1, l2, t3_full, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=False):
contract = self.contract
if isinstance(l1, torch.Tensor):
Dvv = torch.zeros_like(l1)
Expand All @@ -339,7 +416,7 @@ def build_cc3_Dvv(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv
V = F - self.ccwfn.H.F.clone()
else:
V = F - self.ccwfn.H.F.copy()
t3 -= t3_pert_ijk(o, v, i, j, k, t2, V, F, contract)
t3 -= t3_pert_ijk(o, v, i, j, k, t2, t3_full, V, F, contract)
l3 = l3_ijk(i, j, k, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract)
Dvv += 0.5 * contract('bdc,adc->ab', t3, l3)

Expand Down Expand Up @@ -610,5 +687,4 @@ def build_Mvv(self, no, nv, ints, t1):
Mvv = Mvv - contract('ie,ia->ae', ints[:no,-nv:], t1)

return Mvv



92 changes: 66 additions & 26 deletions pycc/cclambda.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
from opt_einsum import contract
from .utils import helper_diis
import torch
from .cctriples import t3c_ijk, l3_ijk, l3_ijk_alt, t3_pert_ijk
from .cctriples import t3c_ijk, l3_ijk, l3_ijk_alt, t3_pert_ijk, t3_ijkabc_alt


class cclambda(object):
Expand Down Expand Up @@ -342,20 +342,51 @@ def residuals(self, F, t1, t2, l1, l2):
Zmndi = np.zeros_like(t2[:,:,:,:no])
Zmdfa = np.zeros_like(t2)
Zmdfa = np.pad(Zmdfa, ((0,0), (0,nv-no), (0,0), (0,0)))
for m in range(no):
for n in range(no):
for l in range(no):
t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
if self.ccwfn.real_time is True:
if isinstance(t1, torch.Tensor):
V = F - self.ccwfn.H.F.clone()
else:
V = F - self.ccwfn.H.F.copy()
t3_lmn -= t3_pert_ijk(o, v, l, m, n, t2, V, F, contract)
Zmndi[m,n] += contract('def,ief->di', t3_lmn, ERI[o,l,v,v])
Zmndi[m,n] -= contract('fed,ief->di', t3_lmn, L[o,l,v,v])
Zmdfa[m] += contract('def,ea->dfa', t3_lmn, ERI[n,l,v,v])
Zmdfa[m] -= contract('dfe,ea->dfa', t3_lmn, L[n,l,v,v])
"""
if self.ccwfn.real_time is True:
t3_full = np.zeros((no,no,no,nv,nv,nv))
for i in range(no):
for j in range(no):
for k in range(no):
for a in range(nv):
for b in range(nv):
for c in range(nv):
t3_full[i,j,k,a,b,c] = t3_ijkabc(o, v, i, j, k, a, b, c, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
"""
if self.ccwfn.real_time is True:
if isinstance(t1, torch.Tensor):
V = F - self.ccwfn.H.F.clone()
else:
V = F - self.ccwfn.H.F.copy()
if abs(V[0,0]) >= 1E-6:
t3_full = t3_ijkabc_alt(o, v, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
for m in range(no):
for n in range(no):
for l in range(no):
t3_lmn = t3_full[l,m,n]
t3_lmn -= t3_pert_ijk(o, v, l, m, n, t2, t3_full, V, F, contract)
Zmndi[m,n] += contract('def,ief->di', t3_lmn, ERI[o,l,v,v])
Zmndi[m,n] -= contract('fed,ief->di', t3_lmn, L[o,l,v,v])
Zmdfa[m] += contract('def,ea->dfa', t3_lmn, ERI[n,l,v,v])
Zmdfa[m] -= contract('dfe,ea->dfa', t3_lmn, L[n,l,v,v])
else:
for m in range(no):
for n in range(no):
for l in range(no):
t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
Zmndi[m,n] += contract('def,ief->di', t3_lmn, ERI[o,l,v,v])
Zmndi[m,n] -= contract('fed,ief->di', t3_lmn, L[o,l,v,v])
Zmdfa[m] += contract('def,ea->dfa', t3_lmn, ERI[n,l,v,v])
Zmdfa[m] -= contract('dfe,ea->dfa', t3_lmn, L[n,l,v,v])
else:
for m in range(no):
for n in range(no):
for l in range(no):
t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
Zmndi[m,n] += contract('def,ief->di', t3_lmn, ERI[o,l,v,v])
Zmndi[m,n] -= contract('fed,ief->di', t3_lmn, L[o,l,v,v])
Zmdfa[m] += contract('def,ea->dfa', t3_lmn, ERI[n,l,v,v])
Zmdfa[m] -= contract('dfe,ea->dfa', t3_lmn, L[n,l,v,v])
if isinstance(t1, torch.Tensor):
Y1 = torch.zeros_like(l1)
Y2 = torch.zeros_like(l2)
Expand Down Expand Up @@ -389,17 +420,26 @@ def residuals(self, F, t1, t2, l1, l2):
Zjlid_1 = np.zeros_like(l2[:,:,:no,:])
Zjlid_2 = np.zeros_like(l2[:,:,:no,:])
# t3l1
for l in range(no):
for m in range(no):
for n in range(no):
t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
if self.ccwfn.real_time is True:
if isinstance(t1, torch.Tensor):
V = F - self.ccwfn.H.F.clone()
else:
V = F - self.ccwfn.H.F.copy()
t3_lmn -= t3_pert_ijk(o, v, l, m, n, t2, V, F, contract)
Znf[n] += contract('de,def->f', l2[l,m], (t3_lmn - t3_lmn.swapaxes(0,2)))
if self.ccwfn.real_time is True:
if abs(V[0,0]) >= 1E-6:
for l in range(no):
for m in range(no):
for n in range(no):
t3_lmn = t3_full[l,m,n]
t3_lmn -= t3_pert_ijk(o, v, l, m, n, t2, t3_full, V, F, contract)
Znf[n] += contract('de,def->f', l2[l,m], (t3_lmn - t3_lmn.swapaxes(0,2)))
else:
for l in range(no):
for m in range(no):
for n in range(no):
t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
Znf[n] += contract('de,def->f', l2[l,m], (t3_lmn - t3_lmn.swapaxes(0,2)))
else:
for l in range(no):
for m in range(no):
for n in range(no):
t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True)
Znf[n] += contract('de,def->f', l2[l,m], (t3_lmn - t3_lmn.swapaxes(0,2)))
for m in range(no):
Y1 += contract('idf,dfa->ia', l2[:,m], Zmdfa[m])
Y1 += contract('iaf,f->ia', L[o,m,v,v], Znf[m])
Expand Down
Loading