(*

This code is written by Antun Balaz (antun@ipb.ac.rs) from the Scientific Computing Laboratory of the Institute of Physics Belgrade, http://www.scl.rs/

The code calculates general many-body effective potential W[x, xbar; eps] to level eps^(p-1) for time independent potential V[x]. The use of level p effective potential for calculation of short-time amplitudes leads to substantially reduces errors, proportional to eps^(p+1-Pd/2), where P is the number of particles, and d is spatial dimensionality. In Path-Integral Monte Carlo calculations of transition amplitudes, the error is proportional to 1/N^p, where N is the number of time steps used in the discretization. The effective potential is defined as

W[x, xbar, eps] = Sum[W[m, k] eps^(m-k), {m, 0, p-1}] ,

and the contractions W[m, k] (homogeneous functions in delta of order 2 k) are calculated by iterating the recursion given in Eq. (33) from the paper "Recursive Schroedinger Equation Approach to Faster Converging Path Integrals" by A. Balaz, A. Bogojevic, I. Vidanovic, and A. Pelster, Phys. Rev. E 79, 036701 (2009), available from http://arxiv.org/abs/0806.4774

This code uses MathTensor package, available from http://smc.vnet.net/mathtensor.html

The effective actions to level p=10 are available in the file EffectiveAction-ManyBody-p10.m from http://www.scl.rs/speedup/

*)

(* MathTensor initialisation. *)
Off[Syntax::com]
Off[SetDelayed::write]
Off[Set::write]
Off[General::spell]
Off[General::spell1]
<< MathTens.m
MakeMetricFlat[Euclid]
On[EvaluateODFlag]

(* Definitions. *)
DumIndCountl = 0;
DumIndCountu = 0;
DefineTensor[d, {{1}, 1}]
GenNewInd[nnn_] := tt /. {DumIndCountl = DumIndCountu + 1; DumIndCountu = DumIndCountl + nnn - 1; nn = DumIndCountl; mm = DumIndCountu; For[kk = nn, kk ≤ mm, kk++, upi[kk] = StringJoin["up", ToString[kk]]; loi[kk] = StringJoin["lo", ToString[kk]]];
For[kk = nn; upil = {}; loil = {}, kk ≤ mm, kk++, ToExpression[StringJoin["AppendTo[upil, up", ToString[kk], "]"]]; ToExpression[StringJoin["AppendTo[loil, lo", ToString[kk], "]"]]];
For[kk = nn; upils[kk - 1] =  ""; loils[kk - 1] = "", kk ≤ mm, kk++, upils[kk] = StringJoin[upils[kk - 1], If[kk > nn, ", ", ""], "up", ToString[kk]]; loils[kk] = StringJoin[loils[kk - 1], If[kk > 1, ", ", ""], "lo", ToString[kk]]]; UpLo[upil, loil]; tt -> 1}
NewDefUnique[contr_] := tempnew /. {temptemp = contr; temp1 = Union[Arglist[temptemp]]; temp2 = Length[temp1]; For[templow = {}; tempup = {}; kkk = 1, kkk ≤ temp2, kkk++, If[LowerIndexQ[temp1[[kkk]]], AppendTo[templow, temp1[[kkk]]], AppendTo[tempup, temp1[[kkk]]]]]; temp1 = Intersection[templow, Lower[tempup]]; temp2 = Length[temp1]; GenNewInd[temp2];For[mcount = 1; mcountt = DumIndCountl, mcount ≤ temp2, mcount++; mcountt++, temptemp = temptemp /. {temp1[[mcount]] -> ToExpression[loi[mcountt]], Pair[temp1[[mcount]]] -> ToExpression[upi[mcountt]]}]; tempnew -> temptemp}
NewDerivativeVec[contr_, vec_, ind_] := ttt /. {d1 = contr; d2 = Union[Arglist[d1]]; d3 = Length[d2]; For[dlow = {}; dup = {}; dk = 1, dk ≤ d3, dk++, If[LowerIndexQ[d2[[dk]]], AppendTo[dlow, d2[[dk]]], AppendTo[dup, d2[[dk]]]]]; d2 = Intersection[dlow, Lower[dup]]; d3 = Length[d2]; ttt1 = 0; For[mt = 1, mt ≤ d3, mt++, ttt1 += (D[contr, vec[d2[[mt]]]] + D[contr, vec[Raise[d2[[mt]]]]] /. {Raise[d2[[mt]]] -> Pair[ind], d2[[mt]] -> Pair[ind], OD[vec[_], _] -> 0})]; ttt1 += D[contr, vec[ind]] /.  OD[vec[_], _] -> 0; ttt -> ttt1}
NewLaplacianVec[contr_, vec_] := ttt2 /. {GenNewInd[1]; ttt3 = NewDerivativeVec[NewDerivativeVec[contr, vec, ToExpression[upi[DumIndCountl]]], vec, ToExpression[loi[DumIndCountl]]]; ttt2 -> ttt3}

(* Calculation of the effective action. *)
p = 10;
CurrentTime = FromDate[Date[]];
StartTime = CurrentTime;
W[i_, j_] = 0;
m= 0; k = 0;
W[m, k] = V;
GenNewInd[2 p + 2];
For[m = 1, m < p, m++, temp = ""; For[k1 = DumIndCountl, k1 ≤ DumIndCountl + 2 m - 1, k1++, temp = StringJoin[temp, " d[", upi[k1], "]"]]; temp = StringJoin[temp, " OD[V, ", loils[DumIndCountl + 2 m - 1], "]"]; W[m, m] = ToExpression[temp] / (2 m + 1)!]
For[m = 1, m < p, m++, For[k = m - 1, k ≥ 0, k--, GenNewInd[1]; tempnum = DumIndCountl; W[m, k] = CanAll[CanAll[(OD[W[m - 1, k], ToExpression[upi[tempnum]], ToExpression[loi[tempnum]]] + NewLaplacianVec[W[m, k + 1], d] - Sum[OD[W[l, r], ToExpression[upi[tempnum]]] NewDefUnique[OD[W[m - l - 2, k - r], ToExpression[loi[tempnum]]]],{l, 0, m - 2}, {r, Max[0, k - m + l + 2], Min[k, l]}] - Sum[NewDerivativeVec[W[l, r], d, ToExpression[upi[tempnum]]] NewDefUnique[NewDerivativeVec[W[m - l - 1, k - r + 1], d, ToExpression[loi[tempnum]]]],{l, 1, m - 2}, {r, Max[0, k - m + l + 2], Min[k, l]}]) / 8 / (m + k + 1) /. {OD[d[_], _] -> 0, OD[d[_], _, _] -> 0}]]; n1 = W[m, k]; n2 = Union[Arglist[n1]]; n3 = Length[n2]; For[nlow = {}; nup = {}; nk = 1, nk ≤ n3, nk++, If[LowerIndexQ[n2[[nk]]], AppendTo[nlow, n2[[nk]]], AppendTo[nup, n2[[nk]]]]]; n2 = Intersection[nlow, Lower[nup]]; n3 = Length[n2]; GenNewInd[n3]; For[mn = 1, mn ≤ n3, mn++, n1 = n1 /. {n2[[mn]] -> ToExpression[loi[DumIndCountl + mn - 1]], Pair[n2[[mn]]] -> ToExpression[upi[DumIndCountl + mn - 1]]}]; For[mn = 1, mn ≤ n3, mn++, n1 = n1 /. {ToExpression[loi[DumIndCountl + mn - 1]] -> ToExpression[loi[mn]], ToExpression[upi[DumIndCountl + mn - 1]] -> ToExpression[upi[mn]]}]; W[m, k] = n1]; EndTime = FromDate[Date[]]; CalculationTime = EndTime - StartTime; Print[]; Print["End of p = ", m + 1, " calculation."]; Print["Calculation time = ", CalculationTime, " seconds."]; Print["Memory used: ", N[MaxMemoryUsed[]/1024/1024], " MB."]; Print[]]

(* Exporting definitions of all W[m, k] to EffectiveAction-ManyBody-p10.m file. *)
Definition[W] >>> "EffectiveAction-ManyBody-p10.m";
(* After exporting all the definitions it is necessary replace uu -> up and ll -> lo in EffectiveAction-ManyBody-p10.m file. *)

(* Reading input from EffectiveAction-ManyBody-p10.m. Beforehand you should initialise MathTensor and load the definitions. GenNewInd should be invoked with the argument sufficient to produce enough dummy indices for a given p, at least 2 p + 2. Make sure that the replacements uu -> up and ll -> lo are made in EffectiveAction-ManyBody-p10.m. *)
p = 10;
GenNewInd[2 p + 2];
<< "EffectiveAction-ManyBody-p10.m";