# Maple V program to make algebraic expression corresponding to the sum of
# all diagrams contributing to the coupled-cluster equation for a
# cluster of dimension with an interaction of dimension ,
# with clusters through cluster size .
# Frank E. Harris
# Quantum Theory Project, University of Florida
# P. O. Box 118435, Gainesville FL 32611
# and
# Department of Physics
# University of Utah
# Salt Lake City, UT 84112
# Copyright (c) Frank E. Harris, 1998.
# This program has been tested on Maple V, Releases 3 and 5; it will not
# run on versions of Maple V earlier than Release 3.
# Clusters and interaction fragments are antisymmetrized, with exact
# definitions as given in Harris, Monkhorst, Freeman, "Algebraic and
# Diagrammatic Methods in Many-Fermion Theory" (Oxford, New York, 1992).
# The lists OCC and VIR contain the index names of, respectively,
# orbitals that are occupied or unoccupied (virtual) in the reference
# state. The first orbitals in each list are (in order) the
# orbitals corresponding to open lines. These lists must be long
# enough to provide all labels needed for all diagrams, and may be
# changed to desired symbols by the user.
OCC := [a,b,c,d,e,f];
VIR := [p,q,r,s,t,u];
# As an example, if this program is run with nfree=2 and the above
# lists OCC and VIR, the cluster whose evaluation is given by the program
# output will be t(p,q,-a,-b), where the minus signs indicate incoming
# lines, with p and a on one vertex, q and b on the other. All t clusters,
# as well as interaction fragments--denoted h(...)--, will have connections
# with the first + and the first - symbol on the same vertex, etc.
# If a set of summation indices is enclosed in parentheses, only distinct
# index sets are to be used in the summation (i.e., (cd) with the above
# specified OCC would mean a double sum over occupied orbitals with c
# preceding d in the orbital set).
# In this version of the program, the output is not in a form permitting
# its use as input to further programs (mainly due to the fact that the
# sums are "inert" and the summation index specifications nonstandard).
# For use in typical molecular systems, one would run the program (for
# each relevant nonzero ) for interactions of dimensions 1 and 2
# and add the dim-1 and dim-2 expressions together. The output for
# nfree=0 generated in this way will not include the reference energy.
# For use with restricted Hartree-Fock orbitals, one would
# (1) Run for interaction of dimension d=2 only;
# (2) For all t, drop all terms in which h has + and - the same occupied
# orbital (these correspond to diagrams with "bubbles")
# (3) For the energy (nfree=0), drop all terms with bubbles except the
# "double-bubble"; retain it with its sign changed.
CC_Diagrams := proc(nfree,d,maxt) local TT;
TT := make_T(nfree,d,maxt); # Obtains structures of all diagrams
# in terms of types (p vs h) and numbers
# of lines
TT := Free_T(TT,nfree); # Assigns open lines in all possible
# ways and causes identical structures
# to be combined
TT := Bound_T(TT,nfree); # Assigns closed-line indices, determines
# diagram sign, and sums the closed-line
# indices
end;
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# For each possible sum of cluster dimensions, find all partitionings
# of this sum into individual cluster dimensions, determine the number
# of closed lines, and write each possible diagram as a symbolic
# function T whose exact definition is given in connection with
# procedure T_prod (next below).
make_T := proc(nfree,d,maxt) local ttot,nc,t1,t2,TT,P;
TT := 0;
t1 := max(nfree-d,0);
t2 := nfree+d;
for ttot from t1 to t2 do
nc := d + ttot - nfree;
P := part0(ttot,maxt); # part0 is initial partitioning (as few
# and as large clusters as permitted
while P <> 0 do TT := TT + T_prod(nc,d,op(P));
P := part(P) od # part gets next partitioning (zero if
# there is no next)
od; TT end;
# Makes sum of all products of h and t-clusters with specified numbers of
# closed lines.
# Input: n=total number of closed lines, d=dimension of H, t1,...,tx
# are cluster dimensions
# Output: sum of functions T(n,d,ch,t1,u1,l1,..,tx,ux,lx),
# where n,d,t1,..,tx are same as input, ch is number of bubbles
# in h, u1,..,ux are numbers of closed particle lines connecting to t1,
# ...,tx and l1,..,lx are corresponding numbers of closed hole lines.
# Combines terms of identical structure. Coefficient of each T is
# number of identical terms in the product, divided by scaling to
# allow for multiple clusters of same dimension (this produces the
# correct diagram coefficients).
T_prod := proc() local TT,n,d,i;
if nargs<2 then ERROR(`bad input`) fi;
n := args[1];
d := args[2];
if not type(n,integer) or n<0 then ERROR(`bad input`) fi;
if not (type(d,integer) and d>0) then ERROR(`bad input`) fi;
TT := T_sumH(n,d); # Makes sum of T for h only (for all
# numbers of closed lines within limit)
for i from 3 to nargs do
TT := TT &p T_sum(n,args[i]) od; # Appends to each T contributions of all
# t (for all numbers of closed lines
# within limit)
TT := Z_delete(TT); # Removes all T with incorrect total
# numbers of closed lines
TT := T_order(TT); # sorts cluster ordering so that
# terms of identical structure combine
TT/T_scale(args) end; # imposes scale factor arising from
# the coupled-cluster expansion
# Removes all concatenated product terms with incorrect total numbers of
# closed lines
Z_delete := proc(x) local TT,i,n,d,m,m1,m2;
if type(x,`+`) then map(Z_delete,x)
elif type(x,function) and op(0,x)=T then
n := op(1,x); d := op(2,x); m := op(3,x); m1 := m; m2 := m;
for i from 5 by 3 to nops(x) do
m1 := m1 + op(i,x);
m2 := m2 + op(i+1,x) od;
if m1>d or m2>d then 0
elif m1+m2-m <> n then 0
else x fi
else 0 fi end;
# Sorts terms with clusters of same size to standard ordering so that
# terms of identical structure will combine
T_order := proc(x) local TT,i,j,i1,i2,j1,j2;
if x=0 then 0
elif type(x,`+`) then map(T_order,x)
elif type(x,function) and op(0,x)=T then
TT := x;
for i from 4 by 3 to nops(x)-3 do
for j from i+3 by 3 to nops(x) do
if op(i,TT)=op(j,TT) and (op(j+1,TT)>op(i+1,TT) or
op(j+1,TT)=op(i+1,TT) and op(j+2,TT)>op(i+2,TT))
then j1 := op(j+1,TT); j2 := op(j+2,TT);
i1 := op(i+1,TT); i2 := op(i+2,TT);
TT := subsop(i+1=j1,i+2=j2,j+1=i1,j+2=i2,TT) fi od od;
TT
else 'procname(args)' fi end;
# Calculates product of nt!, where the nt are the numbers of clusters
# of each size. Proper scaling obtained by dividing by T_scale.
T_scale := proc() local i,j,nL,L,R;
if nargs<3 then RETURN(1) fi;
L := [args[3..nargs]];
L := sort(L); nL := nops(L); R := 1;
i := 1; while i 0 then
TT := 0;
PO := perm0(nfree); # Initial permutation (ascending
# order)
while PO <> 0 do
PV := perm0(nfree); # Initial permutation
while PV <> 0 do
TT := TT + Build_T(x,PO,PV,nfree); # Does index insertion
PV := perm(PV,nfree) od; # Get next permutation
# (zero if done)
PO := perm(PO,nfree) od; # Get next permutation or zero
TT
else TT := Build_T(x,[NULL],[NULL],0) fi # case no open lines
elif type(x,constant) then x
else 'procname(args)' fi end;
# Does actual insertion of open-line indices with scale factor if
# multiple lines of same type from same t or h fragment
Build_T := proc(T,PO,PV,nfree) local i,j,n,uf,ul,nfocc,nfvir,TT,TX,S;
TX := 1;
nfocc:=0; nfvir:=0;
for i from 4 by 3 to nops(T) do
n := op(i,T); uf := n-op(i+1,T); ul := n-op(i+2,T);
S := [];
for j from 1 to uf do nfvir:=nfvir+1; S := [op(S),VIR[PV[nfvir]]] od;
S := sort(S); TT:= [op(S),seq(0,j=uf+1..n)];
S := [];
for j from 1 to ul do nfocc:=nfocc+1; S := [op(S),-OCC[PO[nfocc]]] od;
S := sort(S); TT := t(op(TT),op(S),seq(0,j=ul+1..n));
TX := TX*TT/(uf!*ul!)
od;
n := op(2,T);
S := [];
for i from nfvir+1 to nfree do S := [op(S),VIR[PV[i]]] od;
S := sort(S); TT := [op(S),seq(0,j=nfree-nfvir+1..n)];
S := [];
for i from nfocc+1 to nfree do S := [op(S),-OCC[PO[i]]] od;
S := sort(S); TT := h(op(TT),op(S),seq(0,j=nfree-nfocc+1..n));
TX*TT/((nfree-nfvir)!*(nfree-nfocc)!) end;
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Assigns closed lines, diagram sign, and applies summations
# Keeps track of closed-line indices used and groups indices of same
# type in same fragment within parentheses for summation
# Global variable SX is used to carry index information from subprocedure
# that assigns closed lines to h
# Determines sign by building opeerator product as a list P, then
# calling procedure T_sign
Bound_T := proc(x,nfree) global SX;
local P,Q,TT,SI,nocc,nvir,i,j,xx,n,nmax,ii,xxx;
if type(x,`+`) then map(Bound_T,x,nfree)
elif type(x,`*`) then TT := 1; nocc:=nfree; nvir:=nfree; P := []; SI := NULL;
for i from 1 to nops(x)-1 do xx:=op(i,x);
# if numeric, simply keep factor
# if t, write in indices
# if of form t^p, process t p times
if not type(xx,numeric) then
if type(xx,`^`) then nmax:=op(2,xx); xxx:=op(1,xx)
else nmax:=1; xxx:=xx fi;
for ii from 1 to nmax do xx:=xxx;
n:=nops(xx)/2; SX := NULL;
for j from 1 to n do if op(j,xx)=0 then nvir:=nvir+1;
SX:=cat(SX,VIR[nvir]); xx := subsop(j=VIR[nvir],xx) fi od;
if length(SX) > 1 then SI := cat(SI,`(`,op(SX),`)`)
else SI := cat(SI,SX) fi; SX := NULL;
for j from 1 to n do if op(n+j,xx)=0 then nocc:=nocc+1;
SX:=cat(SX,OCC[nocc]); xx := subsop(n+j=-OCC[nocc],xx) fi od;
Q := [op(xx)]; for j from 1 to n/2 do
Q := subsop(n+j=op(2*n-j+1,Q),2*n-j+1=op(n+j,Q),Q) od;
P := [op(P),op(Q)];
if length(SX) > 1 then SI := cat(SI,`(`,op(SX),`)`)
else SI := cat(SI,SX) fi;
TT := xx*TT od;
else TT := xx*TT
fi od;
xx := fill_h(op(nops(x),x),nfree,nocc,nvir); n:= nops(xx)/2;
Q := [op(xx)]; for i from 1 to n/2 do
Q := subsop(n+i=op(2*n-i+1,Q),2*n-i+1=op(n+i,Q),Q) od;
P := [op(Q),op(P)];
SI := cat(SI,SX);
if SI=NULL then T_sign(nfree,P)*xx*TT else
T_sign(nfree,P)*Sum(xx*TT,SI) fi
elif type(x,function) then TT := fill_h(x,nfree,nfree,nfree);
P := [op(TT)]; n:=nops(TT)/2; for i from 1 to n/2 do
P := subsop(n+i=op(2*n-i+1,P),2*n-i+1=op(n+i,P),P) od;
if SX=NULL then T_sign(nfree,P)*TT else
T_sign(nfree,P)*Sum(TT,SX) fi
else 'procname(args)' fi end;
# Fills in closed-line indices of h and puts associated indices in SX.
fill_h := proc(x,nfree,xnocc,xnvir) global SX; local n,xx,nocc,i,nvir;
n := nops(x)/2; xx:=x; SX := NULL;
nocc:=nfree; for i from 1 to n do if op(i,xx) = 0 then nocc:=nocc+1;
xx := subsop(i=OCC[nocc],xx); fi od;
nocc:=xnocc; nvir:=nfree; for i from 1 to n do if op(n+i,xx) = 0 then
if nvir 1 then SX := cat(`(`,SX,`)`) fi;
xx end;
# P is a list of creation/annihilation operators of a term of T, in
# order. + if creation, - if annihilation.
# Occupied orbitals are from list OCC, virtuals from list VIR
# Other input: nfree is excitation degree of T term; permutes the
# open-line indices to their order in OCC and VIR as part of
# the sign determination
# Output is sign of the term, presented as +1 or -1.
T_sign := proc(nfree,P) local PP,S,i,j,i1,j1;
PP := P;
S := 1;
do
if nops(PP)=0 then RETURN(S) fi;
i := 1;
while nfree > 0 and i < nops(PP) and
(member(PP[i], {op(1..nfree,VIR)}) or member(-PP[i],{op(1..nfree,OCC)}))
do i := i + 1 od;
if i = nops(PP) then
for i from 1 to nfree do
if PP[i] <> VIR[i] then
for j from i+1 while PP[j] <> VIR[i] do od;
i1:=PP[i]; j1:=PP[j]; PP := subsop(i=j1,j=i1,PP); S := -S fi od;
for i from 1 to nfree do
if PP[nfree+i] <> -OCC[nfree-i+1] then
for j from i+1 while PP[nfree+j] <> -OCC[nfree-i+1] do od;
i1:=PP[nfree+i]; j1:=PP[nfree+j];
PP := subsop(nfree+i=j1,nfree+j=i1,PP); S := -S fi od;
RETURN(S) fi;
for j from i+1 while PP[j] + PP[i] <> 0 do od;
if member(PP[i],OCC) or member (PP[j],VIR) then S := S * (-1)^(j-i-1)
else S := S * (-1)^(j-i) fi;
PP := subsop(i=NULL,j=NULL,PP);
od end;
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Utility programs--general
# Initial partition of clusters of combined dimension ttot,
# maximum individual cluster size maxt.
part0 := proc(ttot,maxt);
if ttot=0 then []
elif ttot <= maxt then [ttot]
else [maxt,op(part0(ttot-maxt,maxt))] fi end;
# Converts partition P into next allowable partition of same
# combined dimension
part := proc(P) local i,k,n,R;
n := nops(P); if n=0 then RETURN(0) fi;
k := 0;
for i from n by -1 to 1 while P[i]=1 do k := k + 1 od;
if k=n then 0 else R := part0(k+1,P[n-k]-1);
if n-k=1 then [P[1]-1,op(R)]
else [op(1..n-k-1,P),P[n-k]-1,op(R)] fi fi end;
# Sets up list containing first n integers in ascending order
perm0 := proc(n) local i; [seq(i,i=1..n)] end;
# On each call moves to the next permutation of PP; if called when
# PP already in last permutation (descending order), returns zero
perm := proc(PP,n) local i,j,P;
P := PP;
for i from n by -1 to 2 while P[i]