Quit[]

checkIsIndex[a_]:= If[TangentM4`pmQ[a]==True|| TangentM4`Q[a]==True,Return[True],Return[False]]
replaceIteratively[expr_,x_,list_]:=Module[{n=0},expr/.x:>(n++;list[[n]])]
F[pattern_,x_,list_]:=Total@Table[replaceIteratively[pattern,x,newh],{newh,Permutations[list]}]

CDn[x_,listLegs_][Num_.h_[z1_?checkIsIndex]]:=Num  Which[ToString[h]==ToString[listLegs[[1]]],+I p1[x]h[z1],ToString[h]==ToString[listLegs[[2]]],+I p2[x]h[z1],ToString[h]==ToString[listLegs[[3]]],-I p3[x] h[z1]]

CDn2[x_,listLegs_][Num_.h_[z1_?checkIsIndex,z2_?checkIsIndex]]:=Num  Which[ToString[h]==ToString[listLegs[[1]]],+I p1[x]h[z1,z2],ToString[h]==ToString[listLegs[[2]]],+I p2[x]h[z1,z2],ToString[h]==ToString[listLegs[[3]]],-I p3[x] h[z1,z2]]

CDn3[x_,listLegs_][Num_.h_[z1_?checkIsIndex,z2_?checkIsIndex,z3_?checkIsIndex]]:=Num  Which[ToString[h]==ToString[listLegs[[1]]],+I p1[x]h[z1,z2,z3],ToString[h]==ToString[listLegs[[2]]],+I p2[x]h[z1,z2,z3],ToString[h]==ToString[listLegs[[3]]],-I p3[x] h[z1,z2,z3]]

<<xAct`xTras`

DefManifold[M4,4,{a,b,c,d,e,f,g,μ1,μ2,μ3,μ4,ν1,ν2,ν3,ν4,ρ1,ρ2,ρ3,ρ4,μ,ν,ρ,σ}]
DefMetric[-1,η[-μ1,-ν1],CD,{"|","∇"},FlatMetric->True]

DefTensor[p1[μ],M4]
DefTensor[p2[μ],M4]
DefTensor[p3[μ],M4]
DefTensor[ptot[μ],M4]
DefTensor[A[μ],M4]
DefTensor[A1[μ],M4]
DefTensor[A2[μ],M4]
DefTensor[A3[μ],M4]

(*Try with Spin 1; I give an explicit total derivative (see the variable "terms";*)

terms=CD[μ][A[-μ]A[ν]A[-ν]]//Expand
Vpermuted=Map[F[#,A,{A1,A2,A3}]&,terms]/.CD[x_]-> CDn[x,{A1,A2,A3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules//Simplification

(*It works; the term is proportional to p1+p2-p3*)

(*I try with another vertex*)

terms=CD[μ][A[-μ]CD[ρ][A[ν]]CD[-ρ][A[-ν]]]//Expand
Vpermuted=Map[F[#,A,{A1,A2,A3}]&,terms]/.CD[x_]-> CDn[x,{A1,A2,A3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules//Simplification

(*It works; the term is proportional to p1+p2-p3*)

(*Try with Spin 1 I give a linear combination of operators (like the Out[18] but with general coefficients*)

DefConstantSymbol[λ1]
DefConstantSymbol[λ2]
DefConstantSymbol[λ3]
DefConstantSymbol[λ4]
terms=({{A, ({{ }, {ν}})}}) ({{A, ({{ν}, { }})}}) (∇^μ)({{A, ({{ }, {μ}})}})+λ1({{A, ({{ }, {μ}})}}) ({{A, ({{ν}, { }})}}) (∇^μ)({{A, ({{ }, {ν}})}})+λ2({{A, ({{ }, {μ}})}}) ({{A, ({{ }, {ν}})}}) (∇^μ)({{A, ({{ν}, { }})}})
Vpermuted=Map[F[#,A,{A1,A2,A3}]&,terms]/.CD[x_]-> CDn[x,{A1,A2,A3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
(* I send p3[a_]-> p1[a]+p2[a]+ptot[a] where I will send ptot[a]-> 0 after all the computations. Then I would like to choose the coefficients λ1,λ2 in order to isolate ptot[a]*)
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical

%/.ptot[a_]->0//Simplification

(*This term Out[30] goes to zero if λ1=λ2=1 such that we recover a total derivative*)

terms=λ1({{A, ({{ }, {μ}})}}) (Subscript[∇, a])(Subscript[∇, ρ])({{A, ({{ }, {ν}})}}) (∇^μ)(∇^ρ)(∇^a)({{A, ({{ν}, { }})}})+λ2(Subscript[∇, a])(Subscript[∇, ρ])({{A, ({{ }, {ν}})}}) (∇^μ)({{A, ({{ }, {μ}})}}) (∇^ρ)(∇^a)({{A, ({{ν}, { }})}})+λ3({{A, ({{ }, {μ}})}}) (∇^μ)(Subscript[∇, a])(Subscript[∇, ρ])({{A, ({{ }, {ν}})}}) (∇^ρ)(∇^a)({{A, ({{ν}, { }})}})
Vpermuted=Map[F[#,A,{A1,A2,A3}]&,terms]/.CD[x_]-> CDn[x,{A1,A2,A3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical

%/.ptot[a_]->0//Simplification

(*Let's work now with spin-2 operators*)

DefTensor[h[μ,ν],M4,Symmetric[{μ,ν}]]
DefTensor[h1[μ,ν],M4,Symmetric[{μ,ν}]]
DefTensor[h2[μ,ν],M4,Symmetric[{μ,ν}]]
DefTensor[h3[μ,ν],M4,Symmetric[{μ,ν}]]

terms=CD[-ν][h[μ,ν]CD[-b]@h[-μ,-a]h[a,b]]//Expand
Vpermuted=Map[F[#,h,{h1,h2,h3}]&,terms]/.CD[x_]-> CDn2[x,{h1,h2,h3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical;
%/.ptot[a_]->0//Simplification

terms=λ1({{h, ({{μ, ν}, { ,  }})}}) (Subscript[∇, b])({{h, ({{ ,  }, {μ, a}})}}) (Subscript[∇, ν])({{h, ({{a, b}, { ,  }})}})+λ2({{h, ({{a, b}, { ,  }})}}) (Subscript[∇, b])({{h, ({{ ,  }, {μ, a}})}}) (Subscript[∇, ν])({{h, ({{μ, ν}, { ,  }})}})+λ3({{h, ({{a, b}, { ,  }})}}) ({{h, ({{μ, ν}, { ,  }})}}) (Subscript[∇, ν])(Subscript[∇, b])({{h, ({{ ,  }, {μ, a}})}})
Vpermuted=Map[F[#,h,{h1,h2,h3}]&,terms]/.CD[x_]-> CDn2[x,{h1,h2,h3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical;
%/.ptot[a_]->0//Simplification

(*Try with a spin-3*)

DefTensor[H[μ,ν,ρ],M4,Symmetric[{μ,ν,ρ}]]
DefTensor[H1[μ,ν,ρ],M4,Symmetric[{μ,ν,ρ}]]
DefTensor[H2[μ,ν,ρ],M4,Symmetric[{μ,ν,ρ}]]
DefTensor[H3[μ,ν,ρ],M4,Symmetric[{μ,ν,ρ}]]

terms=CD[-μ][H[μ,ν,ρ]H[-ρ,b,c]H[-b,-c,-ν]]//Expand
Vpermuted=Map[F[#,H,{H1,H2,H3}]&,terms]/.CD[x_]-> CDn3[x,{H1,H2,H3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical;
%/.ptot[a_]->0//Simplification

terms=λ1({{H, ({{μ, ν, ρ}, { ,  ,  }})}}) ({{H, ({{ , b, c}, {ρ,  ,  }})}}) (Subscript[∇, μ])({{H, ({{ ,  ,  }, {b, c, ν}})}})+λ2({{H, ({{ ,  ,  }, {b, c, ν}})}}) ({{H, ({{ , b, c}, {ρ,  ,  }})}}) (Subscript[∇, μ])({{H, ({{μ, ν, ρ}, { ,  ,  }})}})+({{H, ({{ ,  ,  }, {b, c, ν}})}}) ({{H, ({{μ, ν, ρ}, { ,  ,  }})}}) (Subscript[∇, μ])({{H, ({{ , b, c}, {ρ,  ,  }})}})
Vpermuted=Map[F[#,H,{H1,H2,H3}]&,terms]/.CD[x_]-> CDn3[x,{H1,H2,H3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical;
%/.ptot[a_]->0//Simplification

terms=CD[-μ][CD[-μ1]@H[μ,ν,ρ]CD[μ1]@H[-ρ,b,c]H[-b,-c,-ν]]//Expand
Vpermuted=Map[F[#,H,{H1,H2,H3}]&,terms]/.CD[x_]-> CDn3[x,{H1,H2,H3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical;
%/.ptot[a_]->0//Simplification

terms=λ1({{H, ({{ ,  ,  }, {b, c, ν}})}}) (Subscript[∇, μ])(∇^(μ1))({{H, ({{ , b, c}, {ρ,  ,  }})}}) (Subscript[∇, μ1])({{H, ({{μ, ν, ρ}, { ,  ,  }})}})+({{H, ({{ ,  ,  }, {b, c, ν}})}}) (Subscript[∇, μ])(Subscript[∇, μ1])({{H, ({{μ, ν, ρ}, { ,  ,  }})}}) (∇^(μ1))({{H, ({{ , b, c}, {ρ,  ,  }})}})+λ2(Subscript[∇, μ])({{H, ({{ ,  ,  }, {b, c, ν}})}}) (Subscript[∇, μ1])({{H, ({{μ, ν, ρ}, { ,  ,  }})}}) (∇^(μ1))({{H, ({{ , b, c}, {ρ,  ,  }})}})
Vpermuted=Map[F[#,H,{H1,H2,H3}]&,terms]/.CD[x_]-> CDn3[x,{H1,H2,H3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical;
%/.ptot[a_]->0//Simplification

(*However, It's very difficult to simplify the following sum of terms;*)

DefConstantSymbol[D12]
DefConstantSymbol[C15]

terms=-(D12*CD[a][H[-d,-e,-f]]*CD[-b][H[c,d,e]]*CD[-c][H[-a,b,f]])/2-(D12*H[c,d,e]*CD[-b][CD[-a][H[-d,-e,-f]]]*CD[-c][H[a,b,f]])/2-(C15*H[a,c,d]*CD[-b][H[-c,e,f]]*CD[-d][CD[-a][H[b,-e,-f]]])/2-(C15*H[-c,e,f]*CD[-b][H[-a,c,d]]*CD[-d][CD[a][H[b,-e,-f]]])/2-(C15*H[a,c,d]*H[-c,e,f]*CD[-d][CD[-b][CD[-a][H[b,-e,-f]]]])/2-(C15*H[a,c,d]*CD[-b][H[-c,e,f]]*CD[-d][CD[b][H[-a,-e,-f]]])/2-(C15*H[-c,e,f]*CD[-b][H[-a,c,d]]*CD[-d][CD[b][H[a,-e,-f]]])/2;
Vpermuted=Map[F[#,H,{H1,H2,H3}]&,terms]/.CD[x_]-> CDn3[x,{H1,H2,H3}];
FeynRules=Sum[Vpermuted[[i]],{i,1,Length[Vpermuted]}];
FeynRules/.p3[a_]-> p1[a]+p2[a]+ptot[a]//Expand//ToCanonical;
%/.ptot[a_]->0/.D12-> λ1 C15/.C15-> 1//ToCanonical