#  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 <nfree> with an interaction of dimension <d>,
#  with clusters through cluster size <maxt>.

#  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 <nfree> 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 <nfree>) 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<nL do
      j := i;  while j<nL and op(j+1,L)=op(i,L) do j:=j+1 od;
      R := R*(j-i+1)!; i:=j+1
                      od;
R end;

# ------------------------------------------------------------------------

#  Further programs used in subprocedures of make_T

# Sum of T for all diagrams with <= n bubbles within H of dimension d

# Output is sum of T(n,d,# of bubbles)

T_sumH := proc(n,d) local TT,m,i; options remember;
TT := 0;
m := min(n,d);
for i from 0 to m do TT := TT + T(n,d,i) od end;

#  Sum of T for all diagrams with <= n closed lines connecting to  a
#  cluster of dimension d.  Closed lines are between upper--outgoing
#  particle (or lower--incoming hole) index of cluster and lower--in-
#  coming (or upper--outgoing) index of H.

#  Output is sum of T(d,u,l), where u and l are numbers of closed-line
#  upper and lower indices of the cluster.

T_sum := proc(n,d) local TT,m1,m2,i,j; options remember;
m1 := min(n,d);
if m1=0 then RETURN(0) fi;
TT := 0;
for i from 0 to m1 do m2 := min(n-i,d); for j from 0 to m2 do
   if not (i=0 and j=0) then TT := TT + T(d,i,j) fi od od; end;

#  Expands products of cluster sums, making concatenated terms (without
#  requiring conditions on total numbers of closed lines or identifying
#  identical terms)

`&p` := proc(x,y) local TT,i;
if x=0 or y=0 then 0
elif type(x,`+`) then TT := 0;
   for i from 1 to nops(x) do TT := TT + op(i,x) &p y od
elif type(y,`+`) then TT := 0;
   for i from 1 to nops(y) do TT := TT + x &p op(i,y) od
elif type(x,function) and type(y,function) and op(0,x)=T and op(0,y)=T then
   T(op(x),op(y))
else 'procname(args)' fi end;

# ------------------------------------------------------------------------
# ------------------------------------------------------------------------

#  Insert open-line indices in all permutations.  Terms of identical
#  structure will combine.

Free_T := proc(x,nfree) local TT,PO,PV;
if type(x,{`+`,`*`}) then map(Free_T,x,nfree)
elif type(x,function) and op(0,x)=T then
   if nfree <> 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<xnvir then nvir:=nvir+1; 
                                  xx:=subsop(n+i=-VIR[nvir],xx)
                    else nocc:=nocc+1; xx:=subsop(n+i=-OCC[nocc],xx);
                         SX := cat(SX,OCC[nocc]) fi fi od;
if length(SX) > 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]<P[i-1] do od;
if i=1 then RETURN(0) fi;
for j from n by -1 while P[j]<P[i-1] do od;
P := subsop(i-1=P[j],j=P[i-1],P);
for j from i to (i+n-1)/2 do P := subsop(j=P[n-j+i],n-j+i=P[j],P) od;
P; end;

# ------------------------------------------------------------------------
# -----------------------end of program-----------------------------------