(*Remark: Not all the presented options and parameters used for the calculation are necessary. There are probably more efficient methods to implement such derivation in WM14. Therefore, we present here what worked for us, without aiming to optimize the derivations.*)
Clear["Global`*"];
(*Derivation of the rectangle-rectangle vdW interaction*)

(*Assumptions*)
SetSystemOptions["SimplificationOptions"->{"AssumptionsMaxNonlinearVariables"->25,"AssumptionsMaxVariables"->25}];

(*The dimensions of rectangles are 2*bx/2 hx and 2 by/ 2 hy. However, we first derive the law for infinite half-strip, dimensions are infinity/2hx-infinity/2hy*)
$Assumptions=hx>hy&&bx>0 && by>0&&hx>0 && hy>0&& q2>0&&  q1 ∈ Reals && y1>=q2&&0>=x1 &&Element[m,Integers]&&Ly>0 ;

(*function that helps integration*)
integrateD[x_,v_,l_,u_]:=Module[{int=x,var=v,low=l,up=u},
indef=Simplify[Integrate[int,var],TimeConstraint->1];
upper=Simplify[Limit[indef,var->up,Direction->"FromBelow"],TimeConstraint->1];
lower=Simplify[Limit[indef,var->low,Direction->"FromAbove"],TimeConstraint->1];
res=Simplify[upper-lower,TimeConstraint->1];
Return[res]]

(*Point-pair vdW potential.*)
dist =1/((y1-x1)^2+(y2-x2)^2+(q1^2))^3;

(*Let us derive an interaction potential between two infinite half-strips, HST-HST*)

(*line2-point*)
defX2=integrateD[dist,x2,-hx,hx]
(*rectangle (infinity/2hx)-point*)
defX2X1=integrateD[defX2,x1,-Infinity,0]
(*rectangle(infinity/2hx)-line2*)
defX2X1Y2=integrateD[defX2X1,y2,-hy,hy]
(*rectangle (infinity/2hx)-rectangle(infinity/2hy)*)
defX2X1Y2Y1=integrateD[defX2X1Y2,y1,q2,Infinity]

(*These are the results.*)

defX2=(1/8) ((((hx+y2) (3 (hx^2)+5 (q1^2)+5 (x1^2)-10 x1 y1+5 (y1^2)+6 hx y2+3 (y2^2)))/(((((q1^2)+(((x1-y1))^2)))^2) ((((hx^2)+(q1^2)+(x1^2)-2 x1 y1+(y1^2)+2 hx y2+(y2^2)))^2)))-((3 ArcTan[((-hx-y2)/(Sqrt[(q1^2)+(((x1-y1))^2)]))])/((((q1^2)+(((x1-y1))^2)))^(5/2))))+(1/8) ((((hx-y2) (3 (hx^2)+5 (q1^2)+5 (x1^2)-10 x1 y1+5 (y1^2)-6 hx y2+3 (y2^2)))/(((((q1^2)+(((x1-y1))^2)))^2) ((((hx^2)+(q1^2)+(x1^2)-2 x1 y1+(y1^2)-2 hx y2+(y2^2)))^2)))+((3 ArcTan[((hx-y2)/(Sqrt[(q1^2)+(((x1-y1))^2)]))])/((((q1^2)+(((x1-y1))^2)))^(5/2))));

defX2X1=-((π (-((hx-y2) ((((1/((hx^2)+(q1^2)-2 hx y2+(y2^2)))))^(3/2)) (2 (hx^2)+3 (q1^2)-4 hx y2+2 (y2^2)))-(hx+y2) ((((1/((hx^2)+(q1^2)+2 hx y2+(y2^2)))))^(3/2)) (2 (hx^2)+3 (q1^2)+4 hx y2+2 (y2^2))))/(16 (q1^4)))+(1/8) (-(y1/((q1^2) ((q1^2)+(y1^2)) (hx-y2)))-(y1/((q1^2) ((q1^2)+(y1^2)) (hx+y2)))+(y1/((hx-y2) ((hx^2)+(q1^2)-2 hx y2+(y2^2)) ((hx^2)+(q1^2)+(y1^2)-2 hx y2+(y2^2))))+(y1/((hx+y2) ((hx^2)+(q1^2)+2 hx y2+(y2^2)) ((hx^2)+(q1^2)+(y1^2)+2 hx y2+(y2^2))))+((y1 (3 (q1^2)+2 (y1^2)) ArcTan[((-hx-y2)/(Sqrt[(q1^2)+(y1^2)]))])/((q1^4) ((((q1^2)+(y1^2)))^(3/2))))-((y1 (3 (q1^2)+2 (y1^2)) ArcTan[((hx-y2)/(Sqrt[(q1^2)+(y1^2)]))])/((q1^4) ((((q1^2)+(y1^2)))^(3/2))))-(((hx-y2) (2 (hx^2)+3 (q1^2)-4 hx y2+2 (y2^2)) ArcTan[(y1/(Sqrt[(hx^2)+(q1^2)-2 hx y2+(y2^2)]))])/((q1^4) ((((hx^2)+(q1^2)-2 hx y2+(y2^2)))^(3/2))))-(((hx+y2) (2 (hx^2)+3 (q1^2)+4 hx y2+2 (y2^2)) ArcTan[(y1/(Sqrt[(hx^2)+(q1^2)+2 hx y2+(y2^2)]))])/((q1^4) ((((hx^2)+(q1^2)+2 hx y2+(y2^2)))^(3/2)))));
defX2X1Y2=(1/(8 (q1^4)))(-((π (2 (hx^2)-4 hx hy+2 (hy^2)+(q1^2)))/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))+((π (2 (hx^2)+4 hx hy+2 (hy^2)+(q1^2)))/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))+((2 (2 (hx^2)-4 hx hy+2 (hy^2)+(q1^2)) ArcTan[(y1/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))])/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))-((2 (2 (hx^2)+4 hx hy+2 (hy^2)+(q1^2)) ArcTan[(y1/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))])/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))+((hx y1 (3 (q1^2)+2 (y1^2)) ArcTan[((-hx-hy)/(Sqrt[(q1^2)+(y1^2)]))])/((((q1^2)+(y1^2)))^(3/2)))+((hy y1 (3 (q1^2)+2 (y1^2)) ArcTan[((-hx-hy)/(Sqrt[(q1^2)+(y1^2)]))])/((((q1^2)+(y1^2)))^(3/2)))+((hx y1 (3 (q1^2)+2 (y1^2)) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(y1^2)]))])/((((q1^2)+(y1^2)))^(3/2)))-((hy y1 (3 (q1^2)+2 (y1^2)) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(y1^2)]))])/((((q1^2)+(y1^2)))^(3/2)))-((hx y1 (3 (q1^2)+2 (y1^2)) ArcTan[((-hx+hy)/(Sqrt[(q1^2)+(y1^2)]))])/((((q1^2)+(y1^2)))^(3/2)))+((hy y1 (3 (q1^2)+2 (y1^2)) ArcTan[((-hx+hy)/(Sqrt[(q1^2)+(y1^2)]))])/((((q1^2)+(y1^2)))^(3/2)))-((hx y1 (3 (q1^2)+2 (y1^2)) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(y1^2)]))])/((((q1^2)+(y1^2)))^(3/2)))-((hy y1 (3 (q1^2)+2 (y1^2)) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(y1^2)]))])/((((q1^2)+(y1^2)))^(3/2))));
defX2X1Y2Y1=(1/(8 (q1^4)))(((2 (hx^2) π q2)/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))-((4 hx hy π q2)/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))+((2 (hy^2) π q2)/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))+((π (q1^2) q2)/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))-((2 (hx^2) π q2)/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))-((4 hx hy π q2)/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))-((2 (hy^2) π q2)/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))-((π (q1^2) q2)/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))-((2 (2 (hx^2)-4 hx hy+2 (hy^2)+(q1^2)) q2 ArcTan[(q2/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))])/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q1^2)]))+((2 (2 (hx^2)+4 hx hy+2 (hy^2)+(q1^2)) q2 ArcTan[(q2/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))])/(Sqrt[(hx^2)+2 hx hy+(hy^2)+(q1^2)]))-((2 hx (q1^2) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))+((2 hy (q1^2) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))-((4 hx (q2^2) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))+((4 hy (q2^2) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))+((2 hx (q1^2) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))+((2 hy (q1^2) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))+((4 hx (q2^2) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))+((4 hy (q2^2) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)])));

FullSimplify[defX2X1Y2Y1]

(*The result after simplification. HST-HST: infinite half-strip-infinite half-strip*)
infHSTinfHST[q2_]=(1/(8 (q1^4)))(π ((q1^2) (-(1/(Sqrt[(((hx-hy))^2)+(q1^2)]))+(1/(Sqrt[(((hx+hy))^2)+(q1^2)])))+2 ((Sqrt[(((hx-hy))^2)+(q1^2)])-(Sqrt[(((hx+hy))^2)+(q1^2)]))) q2-((2 (2 (((hx-hy))^2)+(q1^2)) q2 ArcTan[(q2/(Sqrt[(((hx-hy))^2)+(q1^2)]))])/(Sqrt[(((hx-hy))^2)+(q1^2)]))+((2 (2 (((hx+hy))^2)+(q1^2)) q2 ArcTan[(q2/(Sqrt[(((hx+hy))^2)+(q1^2)]))])/(Sqrt[(((hx+hy))^2)+(q1^2)]))+((2 ((q1^2)+2 (q2^2)) ((-hx+hy) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(q2^2)]))]+(hx+hy) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(q2^2)]))]))/(Sqrt[(q1^2)+(q2^2)])));

(*This function has two parts, one is linear in q2, and the other is not. When we apply RIHBI to find R-R law, the part that is linear w.r.t. q2 cancels out*)
(*e[q2_]=π ((q1^2) (-(1/(Sqrt[((k[-1])^2)+(q1^2)]))+(1/(Sqrt[((k[1])^2)+(q1^2)])))+2 ((Sqrt[((k[-1])^2)+(q1^2)])-(Sqrt[((k[1])^2)+(q1^2)]))) q2;*)

k[s_]=hx+s*hy;

v[a_,s_]=((s*a*2 (2 (((k[s]))^2)+(q1^2))ArcTan[(a/(Sqrt[(((k[s]))^2)+(q1^2)]))])/(Sqrt[(((k[s]))^2)+(q1^2)]));

l[a_,s_]=((2*(k[s]) ((q1^2)+2 (a^2)) ArcTan[((k[s])/(Sqrt[(q1^2)+(a^2)]))])/(Sqrt[(q1^2)+(a^2)]));

j[a_]=(1/(8 (q1^4)))*(v[a,1]+v[a,-1]-l[a,-1]+l[a,1]);

(*test*)
FullSimplify[j[q2]-j[2 bx+q2]-j[2 by+q2]+j[2 (bx+by)+q2]-(infHSTinfHST[q2]+infHSTinfHST[q2+2 bx+2 by]-infHSTinfHST[q2+2 bx]-infHSTinfHST[q2+2 by])]

(*So, the final expression for the R-R interaction can be written as*)
recRecDef=FullSimplify[j[q2]-j[2 bx+q2]-j[2 by+q2]+j[2 (bx+by)+q2],TimeConstraint->1,Trig->False]

recRec[q2_,bx_,by_,hx_,hy_,q1_]=(1/(8 (q1^4)))(-((2 (2 (((hx-hy))^2)+(q1^2)) q2 ArcTan[(q2/(Sqrt[(((hx-hy))^2)+(q1^2)]))])/(Sqrt[(((hx-hy))^2)+(q1^2)]))+((2 (2 (((hx+hy))^2)+(q1^2)) q2 ArcTan[(q2/(Sqrt[(((hx+hy))^2)+(q1^2)]))])/(Sqrt[(((hx+hy))^2)+(q1^2)]))+((2 (2 (((hx-hy))^2)+(q1^2)) (2 bx+q2) ArcTan[((2 bx+q2)/(Sqrt[(((hx-hy))^2)+(q1^2)]))])/(Sqrt[(((hx-hy))^2)+(q1^2)]))+((2 (2 (((hx-hy))^2)+(q1^2)) (2 by+q2) ArcTan[((2 by+q2)/(Sqrt[(((hx-hy))^2)+(q1^2)]))])/(Sqrt[(((hx-hy))^2)+(q1^2)]))-((2 (2 (((hx-hy))^2)+(q1^2)) (2 (bx+by)+q2) ArcTan[((2 (bx+by)+q2)/(Sqrt[(((hx-hy))^2)+(q1^2)]))])/(Sqrt[(((hx-hy))^2)+(q1^2)]))-((2 (hx-hy) ((q1^2)+2 (q2^2)) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))+((2 (hx+hy) ((q1^2)+2 (q2^2)) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(q2^2)]))])/(Sqrt[(q1^2)+(q2^2)]))+(-2 (2 (((hx+hy))^2)+(q1^2)) (2 bx+q2) (Sqrt[(q1^2)+(((2 bx+q2))^2)]) ArcTan[((2 bx+q2)/(Sqrt[(((hx+hy))^2)+(q1^2)]))]+2 (Sqrt[(((hx+hy))^2)+(q1^2)]) ((q1^2)+2 (((2 bx+q2))^2)) ((hx-hy) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(((2 bx+q2))^2)]))]-(hx+hy) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(((2 bx+q2))^2)]))]))/((Sqrt[((((hx+hy))^2)+(q1^2)) ((q1^2)+(((2 bx+q2))^2))]))+(-2 (2 (((hx+hy))^2)+(q1^2)) (2 by+q2) (Sqrt[(q1^2)+(((2 by+q2))^2)]) ArcTan[((2 by+q2)/(Sqrt[(((hx+hy))^2)+(q1^2)]))]+2 (Sqrt[(((hx+hy))^2)+(q1^2)]) ((q1^2)+2 (((2 by+q2))^2)) ((hx-hy) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(((2 by+q2))^2)]))]-(hx+hy) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(((2 by+q2))^2)]))]))/((Sqrt[((((hx+hy))^2)+(q1^2)) ((q1^2)+(((2 by+q2))^2))]))+(1/(Sqrt[((((hx+hy))^2)+(q1^2)) ((q1^2)+(((2 (bx+by)+q2))^2))]))2 ((2 (((hx+hy))^2)+(q1^2)) (2 (bx+by)+q2) (Sqrt[(q1^2)+(((2 (bx+by)+q2))^2)]) ArcTan[((2 (bx+by)+q2)/(Sqrt[(((hx+hy))^2)+(q1^2)]))]-(Sqrt[(((hx+hy))^2)+(q1^2)]) ((q1^2)+2 (((2 (bx+by)+q2))^2)) ((hx-hy) ArcTan[((hx-hy)/(Sqrt[(q1^2)+(((2 (bx+by)+q2))^2)]))]-(hx+hy) ArcTan[((hx+hy)/(Sqrt[(q1^2)+(((2 (bx+by)+q2))^2)]))])));

(*A numerical check - everything is correct. The numerical integration struggles with small separations.*)
r={q2->1/10.,bx->1,by->2,hx->1,hy->2,q1->1/10};
recRec@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2+(q1^2))^3/.r,{x1,-bx,bx}/.r,{x2,-hx,hx}/.r,{y2,-hy,hy}/.r,{y1,bx+q2,bx+q2+2*by}/.r]

(*Again, there is a numerical issue for q1->0*)
r={q2->1/1000,bx->1,by->2,hx->1,hy->1/2,q1->1/1000000}
1.*recRec@@r[[All,2]]
N[recRec@@r[[All,2]],16]

(*The limit case for q1->0 does not work. The force shows similar behavior.*)
Limit[recRec[q2,bx,by,hx,hy,q1],q1->0,Direction->"FromAbove"]

(*Test the scaling w.r.t. q2, for different values of parameter c=q1/q2, which is tangens of the angle between q1 and q2. The exact law shows that the factor c does not influence results for c<0.1 (all graphs overlap)*)

tab={1/100,1/10,1,10,100};
f = -D[Log[recRec[q2,1,1,1,1,c*q2]/.{q2->Exp[q2]}],q2];
as1=Table[{q2,Re[N[f,100]]},{q2,-8,8,1/4},{c,tab}];
ListLinePlot[Transpose@as1,AxesLabel->{"Log[q2]","Scaling factor of (Subscript[Π, R - R]) - m=6"},PlotLegends->{tab}]

(*Now, let us integrate the HST-HST law along q1, and find the potential between an infinite half-strip (∞/2hx) and a an infinite half-plate with one finite in-plane dimension (∞/2hy/2Ly), HST-(Subscript[HP, Ly])*)
defQ1=FullSimplify[integrateD[infHSTinfHST[q2][[1;;2]]*infHSTinfHST[q2][[3,2;;4]],q1,-Ly,Ly]]

(*The same result follows if we only integrate the function j[q2],*) 
(*FullSimplify[integrateD[j[q2],q1,-Ly,Ly]]*)

(*The result of integration, HR-(Subscript[HRP, Ly]) law*)
infHSTinfHPL[q2_]=(1/(6 (q2^2)))((((2 (((hx-hy))^2)-(Ly^2)) (Sqrt[(((hx-hy))^2)+(Ly^2)]) (q2^3) ArcTan[(q2/(Sqrt[(((hx-hy))^2)+(Ly^2)]))])/((((hx-hy))^2) (Ly^3)))+(((-2 (((hx+hy))^2)+(Ly^2)) (Sqrt[(((hx+hy))^2)+(Ly^2)]) (q2^3) ArcTan[(q2/(Sqrt[(((hx+hy))^2)+(Ly^2)]))])/((((hx+hy))^2) (Ly^3)))-((((((((hx-hy))^2)+(q2^2)))^(3/2)) ArcTan[(Ly/(Sqrt[(((hx-hy))^2)+(q2^2)]))])/(((hx-hy))^2))+((((((((hx+hy))^2)+(q2^2)))^(3/2)) ArcTan[(Ly/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(((hx+hy))^2))-(((hx-hy) ((Ly^2)-2 (q2^2)) (Sqrt[(Ly^2)+(q2^2)]) ArcTan[((hx-hy)/(Sqrt[(Ly^2)+(q2^2)]))])/(Ly^3))+(((hx+hy) ((Ly^2)-2 (q2^2)) (Sqrt[(Ly^2)+(q2^2)]) ArcTan[((hx+hy)/(Sqrt[(Ly^2)+(q2^2)]))])/(Ly^3)));

(*Now, we apply RIHBI to HST-(Subscript[HP, Ly]), and obtain R-(Subscript[RP, Ly]) (rectangle (2bx/2hx)-finite length rectangular prism (2by/2hy/2Ly))*)
k[s_]=hx+s hy;
ff[a_,s_]=-s (((2 ((k[s])^2)-(Ly^2)) (Sqrt[((k[s])^2)+(Ly^2)]) a ArcTan[(a/(Sqrt[((k[s])^2)+(Ly^2)]))])/(6*((k[s])^2) (Ly^3)));
g[a_,s_]=s (((((((k[s])^2)+(a^2)))^(3/2)) ArcTan[(Ly/(Sqrt[((k[s])^2)+(a^2)]))])/(6 (a^2)((k[s])^2)));
p[a_,s_]=s ((k[s] ((Ly^2)-2 (a^2))(Sqrt[(Ly^2)+(a^2)]) ArcTan[((k[s])/(Sqrt[(Ly^2)+(a^2)]))])/(6 (a^2) (Ly^3)));
j[a_]=ff[a,1]+ff[a,-1]+g[a,1]+g[a,-1]+p[a,1]+p[a,-1];

FullSimplify[j[q2]+j[q2+2*bx+2*by]-j[q2+2*bx]-j[q2+2*by]-(infHSTinfHPL[q2]+infHSTinfHPL[q2+2*bx+2*by]-infHSTinfHPL[q2+2*bx]-infHSTinfHPL[q2+2*by])]

(*The final expression for R-(Subscript[RP, Ly])*)
recRPL[q2_,bx_,by_,hx_,hy_,Ly_]=infHSTinfHPL[q2]+infHSTinfHPL[q2+2*bx+2*by]-infHSTinfHPL[q2+2*bx]-infHSTinfHPL[q2+2*by];

(*A numerical check - everything is correct. Numerical integration struggles with small separations. A singularity for hx->hy is present*)
r={q2->0.1,bx->1,by->2,hx->1.0001,hy->1,Ly->0.5};
recRPL@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2+(q1^2))^3,{x1,-bx,bx}/.r,{x2,-hx,hx}/.r,{y2,-hy,hy}/.r,{y1,bx+q2,bx+q2+2*by}/.r,{q1,-Ly,Ly}/.r]

(*Now, we can derive special cases. First, a vdW potential between a rectangle and a rectangular prism of infinite length. It is more elegant to first derive an interaction between an infinite half-strip and an infinite half-plate, HST-(Subscript[HP, ∞]), and then to apply RIHBI*)
defQ1infX=Limit[infHSTinfHPL[q2],Ly->Infinity]

(*This form is preferable because we will have to fix the final expression*)
HSTinfHP[q2_]=Total[(1/(Denominator[defQ1infX]))*Level[Numerator[defQ1infX],1]]

(*Now we apply RIHBI to HST-(Subscript[HP, ∞]), and obtain the R-RP law. Therefore, we can derive the R-RP law starting from the R-R law. The obtained expression is the same as the one derived by applying RIIPI to R-(Subscript[R, IP]) in Notebook 5*)
recInfRP=HSTinfHP[q2]+HSTinfHP[q2+2*bx+2*by]-HSTinfHP[q2+2*bx]-HSTinfHP[q2+2*by]

(*This expression has singularity for hx->hy that can be fixed*)

limit1=Table[Limit[recInfRP[[i]],hy->hx],{i,1,Length[recInfRP]}]
illTerms=Flatten[Position[Abs[limit1],Infinity]]

illSimplified=FullSimplify[Expand[recInfRP[[illTerms]]]]

(*One part of this expression is well-defined, the other needs fixing*)
ill=FullSimplify[illSimplified[[1;;3]] illSimplified[[4,1;;4]]]
well=FullSimplify[illSimplified[[1;;3]] illSimplified[[4,5;;7]]]

(*How to solve this? In essence, find(Sqrt[a])-(Sqrt[b]), and multiply with (Sqrt[a])+(Sqrt[b]), to get a-b*)
n=Numerator[ill]
d=Denominator[ill]
f1=π (-(Sqrt[(((hx-hy))^2)+(q2^2)])-(Sqrt[(((hx-hy))^2)+(((2 bx+q2))^2)])+(Sqrt[(((hx-hy))^2)+(((2 by+q2))^2)])+(Sqrt[(((hx-hy))^2)+(((2 (bx+by)+q2))^2)]));
fix1=FullSimplify[n*f1]
f2=(-2 (π^2) (4 bx (bx+by+q2)+(Sqrt[((((hx-hy))^2)+(q2^2)) ((((hx-hy))^2)+(((2 by+q2))^2))])+(Sqrt[((((hx-hy))^2)+(((2 bx+q2))^2)) ((((hx-hy))^2)+(((2 (bx+by)+q2))^2))])));
fix2=FullSimplify[fix1*f2]
f3=(-32 bx (π^4) (bx+by+q2) ((((hx-hy))^2)+2 by q2+(q2^2)+(Sqrt[((((hx-hy))^2)+(q2^2)) ((((hx-hy))^2)+(((2 by+q2))^2))])));
fix3=FullSimplify[fix2*f3]

fixF=f1*f2*f3;

fixed=FullSimplify[((FullSimplify[n*fixF])/(FullSimplify[d*fixF]))]

recInfRPfixed=fixed+well+recInfRP[[Complement[Range[1,Length[recInfRP]],illTerms]]]

(*the result*)
recInfRPfixed=((π ((((((hx+hy))^2)+(q2^2)))^(3/2)))/(12 (((hx+hy))^2) (q2^2)))-((π ((((((hx+hy))^2)+(((2 bx+q2))^2)))^(3/2)))/(12 (((hx+hy))^2) (((2 bx+q2))^2)))-((π ((((((hx+hy))^2)+(((2 by+q2))^2)))^(3/2)))/(12 (((hx+hy))^2) (((2 by+q2))^2)))+((π ((((((hx+hy))^2)+(((2 bx+2 by+q2))^2)))^(3/2)))/(12 (((hx+hy))^2) (((2 bx+2 by+q2))^2)))+(1/12) π (-((Sqrt[(((hx-hy))^2)+(q2^2)])/(q2^2))+((Sqrt[(((hx-hy))^2)+(((2 bx+q2))^2)])/(((2 bx+q2))^2))+((Sqrt[(((hx-hy))^2)+(((2 by+q2))^2)])/(((2 by+q2))^2))-((Sqrt[(((hx-hy))^2)+(((2 (bx+by)+q2))^2)])/(((2 (bx+by)+q2))^2)))+(16 bx (by^2) π (bx+by+q2))/(3 ((((hx-hy))^2)+2 by q2+(q2^2)+(Sqrt[((((hx-hy))^2)+(q2^2)) ((((hx-hy))^2)+(((2 by+q2))^2))])) ((Sqrt[(((hx-hy))^2)+(q2^2)])+(Sqrt[(((hx-hy))^2)+(((2 bx+q2))^2)])-(Sqrt[(((hx-hy))^2)+(((2 by+q2))^2)])-(Sqrt[(((hx-hy))^2)+(((2 (bx+by)+q2))^2)])) (4 bx (bx+by+q2)+(Sqrt[((((hx-hy))^2)+(q2^2)) ((((hx-hy))^2)+(((2 by+q2))^2))])+(Sqrt[((((hx-hy))^2)+(((2 bx+q2))^2)) ((((hx-hy))^2)+(((2 (bx+by)+q2))^2))])));

(*represent it as a compact function*)
fff[a_]=((π ((((((hx+hy))^2)+(a^2)))^(3/2)))/(12 (((hx+hy))^2) (a^2)));
gg[a_]=(1/12) π ((Sqrt[(((hx-hy))^2)+(a^2)])/(a^2));

w[a_]=(Sqrt[(((hx-hy))^2)+(a^2)]);
qQ[a_]:=((((((((a^2)+(((hx+hy))^2)))^(3/2))/(((hx+hy))^2))-w[a]) π)/(12 (a^2)));

tT=(16 bx (by^2) π (bx+by+q2))/(3 ((((hx-hy))^2)+2 by q2+(q2^2)+w[q2]*w[2 by+q2]) (w[q2]+w[2 bx+q2]-w[2 by+q2]-w[2 bx + 2 by + q2]) (4 bx (bx+by+q2)+w[q2]*w[2 by+q2]+w[2 bx+q2]*w[2(bx+by)+q2]));



FullSimplify[fff[q2]-fff[2 bx+q2]-fff[2 by+q2]+fff[2 (bx+by)+q2]-gg[q2]+gg[2 bx+q2]+gg[2 by+q2]-gg[2 (bx+by)+q2]-recInfRPfixed[[1;;5]]]

FullSimplify[qQ[q2]-qQ[2 bx+q2]-qQ[2 by+q2]+qQ[2 (bx+by)+q2]+tT-recInfRPfixed]

(*final result R-RP*)
recRP[q2_,bx_,by_,hx_,hy_]=qQ[q2]-qQ[2 bx+q2]-qQ[2 by+q2]+qQ[2 (bx+by)+q2]+tT;

(*Numerical check of the final expression, everything is correct.*)
r={q2->0.1,bx->1,by->2,hx->1,hy->1+1/100};
recInfRP/.r
recRP@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2+(q1^2))^3/.r,{x1,-bx,bx}/.r,{x2,-hx,hx}/.r,{y2,-hy,hy}/.r,{y1,bx+q2,bx+q2+2*by}/.r,{q1,-10,10}]

(*Scaling of rectangle-infinite rectangular prism, R-RP. It is 2 for small and 5 for large separations*)
f = -D[Log[recRP[q2,bx,by,hx,hy]/.{q2->Exp[q2]}],q2]/.{hx->1,hy->2,bx->1,by->1};
Plot[N[f,50],{q2,-10,10},WorkingPrecision->50]

(*The final test, let us compare this R-RP law, with the R-RP law derived by applying RIIPI to R-(Subscript[R, IP]) law, by setting hx->hy. The expression are equal*)
recRPfromRIIPI=((π (-6 (hy^2) q2+(((4 (hy^2)+(q2^2)))^(3/2))))/(48 (hy^2) (q2^2)))-((π (-6 (hy^2) (2 bx+q2)+(((4 (hy^2)+(((2 bx+q2))^2)))^(3/2))))/(48 (hy^2) (((2 bx+q2))^2)))-((π (-6 (hy^2) (2 by+q2)+(((4 (hy^2)+(((2 by+q2))^2)))^(3/2))))/(48 (hy^2) (((2 by+q2))^2)))+((π (-6 (hy^2) (2 bx+2 by+q2)+(((4 (hy^2)+(((2 bx+2 by+q2))^2)))^(3/2))))/(48 (hy^2) (((2 bx+2 by+q2))^2)));
FullSimplify[recRP[q2,bx,by,hy,hy]-recRPfromRIIPI]

(*additional special cases*)

(*Rectangle-infinite plate law. Thickness of the plate is 2 by.*)
recInfPlate=Limit[recRP[q2,bx,by,hx,hy],hy->Infinity]

(*The result is the same as the one obtained from the general rectSecInfPlate23 law for m=6, see Notebook 1.*)
((4 hx π ((q2^(4-m))-(((2 bx+q2))^(4-m))-(((2 by+q2))^(4-m))+(((2 (bx+by)+q2))^(4-m))))/((-4+m) (-3+m) (-2+m)))/.m->6

(*Plate-infinite plate law, per unit area. The thickness of the plates are 2by and 2bx. The result is the same as in Notebook 1.*)
plateInfPlate=FullSimplify[recInfPlate/.hx->1/2]

(*Rectangle-infinite half-space law. *)
recHS=FullSimplify[Limit[recInfPlate,by->Infinity]]

(*This result represents the special case of the general expression in Notebook 1.*)
Simplify[Together[((4 hx π ((q2^(4-m))-(((2 bx+q2))^(4-m))))/((-4+m) (-3+m) (-2+m)))/.m->6]]

(*Plate-infinite half-space law, per unit area of a plate. *)
plateHS=FullSimplify[recHS/.hx->1/2]

(*Again, this result represents the special case of the general expression in Notebook 1.*)
Simplify[Together[((2 π ((q2^(4-m))-(((2 bx+q2))^(4-m))))/((-4+m) (-3+m) (-2+m)))/.m->6]]

(*half-space - half-space law*)
hsHS=FullSimplify[Limit[plateHS,bx->Infinity]]

(*Again, this is the special case of the general expression in Notebook 1.*)
((2 π (q2^(4-m)))/((-4+m) (-3+m) (-2+m)))/.m->6

(*Additionally, in order to check our derivation approach, let us find an interaction between two rectangular sheets (rectangles) as in DeRocco/Parsegian. The difference w.r.t. R-R law is that rectangles are oriented such that they share the normal through the center, perpendicular to their planes. We integrate similarly as for the R-R law, but skip integration over x1 and y1. We integrate along three coordinates (x2,y2,q1=y3) to obtain line2-sheet23 law*)
$Assumptions=hx>hy&&bx>0 && by>0&&hx>0 && hy>0&& q2>0&&  q1 ∈ Reals && bx+q2+2*by>=y1>=bx+q2 &&Abs[x1]<=bx &&Abs[x2]<=hx &&Abs[x3]<=Lx &&Abs[y2]<=hy&&Abs[y3]<=Ly&&Element[x3,Reals]&&Element[m,Integers]&&Ly>0&&Lx>0 &&Element[dh,Reals]&&hh>0&&a1>0&&a2>0;
dist =1/((y1-x1)^2+(y2-x2)^2+(((y3-x3))^2))^3;
(*line2-point law*)
defX2=integrateD[dist/.y1-x1->q2,x2,-hx,hx]
(*line2-line2 law*)
defX2Y2=integrateD[defX2,y2,-hy,hy]
(*line2-sheet23 law*)
defX2Y2Y3=integrateD[defX2Y2,y3,-Ly,Ly]
(*sheet23-sheet23 law, later, first let us simplify previous integral*)
(*defX2Y2Y3X3=integrateD[defX2Y2Y3,x3,-Lx,Lx]*)

(*Some simplification*)
lev=Level[Collect[defX2Y2Y3,ArcTan[_]],1];
FullSimplify[Total[lev[[1;;8]]]];
fs1=FullSimplify[lev[[9;;Length[lev]]],Trig->False];

(*Manual simplification, arctan is an odd function, Arctan[-x]=-Arctan[x]*)
defX2Y2Y3final=Total[FullSimplify[Level[Collect[Total[fs1/.{ArcTan[((-hx+hy)/(Sqrt[a_]))]->-ArcTan[((hx-hy)/(Sqrt[a]))],ArcTan[((-Ly+x3)/(a_))]->-ArcTan[((Ly-x3)/a)],ArcTan[((-Ly+x3)/(Sqrt[a_]))]->-ArcTan[((Ly-x3)/(Sqrt[a]))]}],ArcTan[_]],1],Trig->False]]

defX2Y2Y3final=-(((hx-hy) (3 (q2^2)+2 (((Ly-x3))^2)) (Ly-x3) ArcTan[((hx-hy)/(Sqrt[(q2^2)+(((Ly-x3))^2)]))])/(4 (q2^4) ((((q2^2)+(((Ly-x3))^2)))^(3/2))))+(((hx+hy) (3 (q2^2)+2 (((Ly-x3))^2)) (Ly-x3) ArcTan[((hx+hy)/(Sqrt[(q2^2)+(((Ly-x3))^2)]))])/(4 (q2^4) ((((q2^2)+(((Ly-x3))^2)))^(3/2))))-(((2 (((hx-hy))^2)+(q2^2)) ArcTan[((Ly-x3)/(Sqrt[(((hx-hy))^2)+(q2^2)]))])/(4 (q2^4) (Sqrt[(((hx-hy))^2)+(q2^2)])))+(((2 (((hx+hy))^2)+(q2^2)) ArcTan[((Ly-x3)/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(4 (q2^4) (Sqrt[(((hx+hy))^2)+(q2^2)])))-(((2 (((hx-hy))^2)+(q2^2)) ArcTan[((Ly+x3)/(Sqrt[(((hx-hy))^2)+(q2^2)]))])/(4 (q2^4) (Sqrt[(((hx-hy))^2)+(q2^2)])))+(((2 (((hx+hy))^2)+(q2^2)) ArcTan[((Ly+x3)/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(4 (q2^4) (Sqrt[(((hx+hy))^2)+(q2^2)])))+(((-hx+hy) (Ly+x3) (3 (q2^2)+2 (((Ly+x3))^2)) ArcTan[((hx-hy)/(Sqrt[(q2^2)+(((Ly+x3))^2)]))])/(4 (q2^4) ((((q2^2)+(((Ly+x3))^2)))^(3/2))))+(((hx+hy) (Ly+x3) (3 (q2^2)+2 (((Ly+x3))^2)) ArcTan[((hx+hy)/(Sqrt[(q2^2)+(((Ly+x3))^2)]))])/(4 (q2^4) ((((q2^2)+(((Ly+x3))^2)))^(3/2))));

(*Now, let us integrate once more to obtain the sheet-sheet law*)
defX2Y2Y3X3=integrateD[defX2Y2Y3final,x3,-Lx,Lx];

(*the result*)
defX2Y2Y3X3=(1/(4 (q2^4)))(-((2 Lx (2 (((hx-hy))^2)+(q2^2)) ArcTan[((-Lx+Ly)/(Sqrt[(((hx-hy))^2)+(q2^2)]))])/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))-((2 Lx (2 (((hx-hy))^2)+(q2^2)) ArcTan[((Lx+Ly)/(Sqrt[(((hx-hy))^2)+(q2^2)]))])/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))+((2 Lx (2 (((hx+hy))^2)+(q2^2)) ArcTan[((-Lx+Ly)/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(Sqrt[(((hx+hy))^2)+(q2^2)]))+((2 Lx (2 (((hx+hy))^2)+(q2^2)) ArcTan[((Lx+Ly)/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(Sqrt[(((hx+hy))^2)+(q2^2)]))+((2 (hx-hy) (2 (Lx^2)-4 Lx Ly+2 (Ly^2)+(q2^2)) ArcTan[((hx-hy)/(Sqrt[(((Lx-Ly))^2)+(q2^2)]))])/(Sqrt[(Lx^2)-2 Lx Ly+(Ly^2)+(q2^2)]))-((2 (hx+hy) (2 (Lx^2)-4 Lx Ly+2 (Ly^2)+(q2^2)) ArcTan[((hx+hy)/(Sqrt[(((Lx-Ly))^2)+(q2^2)]))])/(Sqrt[(Lx^2)-2 Lx Ly+(Ly^2)+(q2^2)]))-((2 (hx-hy) (2 (Lx^2)+4 Lx Ly+2 (Ly^2)+(q2^2)) ArcTan[((hx-hy)/(Sqrt[(((Lx+Ly))^2)+(q2^2)]))])/(Sqrt[(((Lx+Ly))^2)+(q2^2)]))+((2 (hx+hy) (2 (Lx^2)+4 Lx Ly+2 (Ly^2)+(q2^2)) ArcTan[((hx+hy)/(Sqrt[(((Lx+Ly))^2)+(q2^2)]))])/(Sqrt[(((Lx+Ly))^2)+(q2^2)]))-(2 (((hx-hy))^2)+(q2^2)) (-((Ly ArcTan[((-Lx+Ly)/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))])/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))+(1/2) Log[(hx^2)-2 hx hy+(hy^2)+(((Lx-Ly))^2)+(q2^2)])+(2 (hx^2)-4 hx hy+2 (hy^2)+(q2^2)) Log[(hx^2)-2 hx hy+(hy^2)+(((Lx-Ly))^2)+(q2^2)]+(2 (((hx+hy))^2)+(q2^2)) (-((Ly ArcTan[((-Lx+Ly)/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(Sqrt[(((hx+hy))^2)+(q2^2)]))+(1/2) Log[(((hx+hy))^2)+(((Lx-Ly))^2)+(q2^2)])-(2 (hx^2)+4 hx hy+2 (hy^2)+(q2^2)) Log[(((hx+hy))^2)+(((Lx-Ly))^2)+(q2^2)]-(2 (((hx-hy))^2)+(q2^2)) (((Ly ArcTan[((Lx-Ly)/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))])/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))+(1/2) Log[(hx^2)-2 hx hy+(hy^2)+(Lx^2)-2 Lx Ly+(Ly^2)+(q2^2)])+(2 (((hx+hy))^2)+(q2^2)) (((Ly ArcTan[((Lx-Ly)/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(Sqrt[(((hx+hy))^2)+(q2^2)]))+(1/2) Log[(((hx+hy))^2)+(Lx^2)-2 Lx Ly+(Ly^2)+(q2^2)])+(2 (((hx-hy))^2)+(q2^2)) (((Ly ArcTan[((-Lx-Ly)/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))])/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))+(1/2) Log[(hx^2)-2 hx hy+(hy^2)+(((Lx+Ly))^2)+(q2^2)])+(2 (((hx-hy))^2)+(q2^2)) (-((Ly ArcTan[((Lx+Ly)/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))])/(Sqrt[(hx^2)-2 hx hy+(hy^2)+(q2^2)]))+(1/2) Log[(hx^2)-2 hx hy+(hy^2)+(((Lx+Ly))^2)+(q2^2)])-(2 (hx^2)-4 hx hy+2 (hy^2)+(q2^2)) Log[(hx^2)-2 hx hy+(hy^2)+(((Lx+Ly))^2)+(q2^2)]-(2 (((hx+hy))^2)+(q2^2)) (((Ly ArcTan[((-Lx-Ly)/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(Sqrt[(((hx+hy))^2)+(q2^2)]))+(1/2) Log[(((hx+hy))^2)+(((Lx+Ly))^2)+(q2^2)])-(2 (((hx+hy))^2)+(q2^2)) (-((Ly ArcTan[((Lx+Ly)/(Sqrt[(((hx+hy))^2)+(q2^2)]))])/(Sqrt[(((hx+hy))^2)+(q2^2)]))+(1/2) Log[(((hx+hy))^2)+(((Lx+Ly))^2)+(q2^2)])+(2 (hx^2)+4 hx hy+2 (hy^2)+(q2^2)) Log[(((hx+hy))^2)+(((Lx+Ly))^2)+(q2^2)]);

(*Additional simplifications*)
sheetSheet=Simplify[Total[FullSimplify[Level[Collect[Total[FullSimplify[Level[Collect[Apart[Simplify[defX2Y2Y3X3]],{ArcTan[_],Log[_]}],1]]]/.{ArcTan[((-Lx+Ly)/(Sqrt[a_]))]->-ArcTan[((Lx-Ly)/(Sqrt[a]))]},{ArcTan[_],Log[_]}],1]]],Trig->False]

kk[s_]=hx+s*hy;
ll[b_]=Lx+b*Ly;

h[s_,b_]=s*b*(((ll[b] (2 ((kk[s])^2)+(q2^2)) ArcTan[((ll[b])/(Sqrt[((kk[s])^2)+(q2^2)]))])/(Sqrt[((kk[s])^2)+(q2^2)]))+((kk[s] (2 ((ll[b])^2)+(q2^2)) ArcTan[((kk[s])/(Sqrt[((ll[b])^2)+(q2^2)]))])/(Sqrt[((ll[b])^2)+(q2^2)])));


lt={h[-1,-1],h[1,1],h[1,-1],h[-1,1]};



FullSimplify[sheetSheet[[1;;2]]*Total[lt]-sheetSheet]
(*FullSimplify[finalDefY1X1Y2X2[[1;;2]]*Total[lt1[[1;;Length[lt1]]]]-finalDefY1X1Y2X2]
FullSimplify[finalDefY1X1Y2X2[[1;;2]]*Total[lt2[[1;;Length[lt2]]]]-finalDefY1X1Y2X2]*)

sheetSheetF[q2_,hx_,hy_,Lx_,Ly_]=sheetSheet[[1;;2]]*Total[lt];

(*Now, let us check the special case, an interaction between two rectangle sheets of the same dimensions*)
sheetSheetSpec=Limit[Limit[sheetSheetF[q2,hx,hy,Lx,Ly],hx->hy],Lx->Ly]

(*The result*)
sheetSheetSpec=(1/(2 (q2^4)))(-2 hy q2 ArcTan[((2 hy)/q2)]-2 Ly q2 ArcTan[((2 Ly)/q2)]+((2 Ly (8 (hy^2)+(q2^2)) ArcTan[((2 Ly)/(Sqrt[4 (hy^2)+(q2^2)]))])/(Sqrt[4 (hy^2)+(q2^2)]))+((2 hy (8 (Ly^2)+(q2^2)) ArcTan[((2 hy)/(Sqrt[4 (Ly^2)+(q2^2)]))])/(Sqrt[4 (Ly^2)+(q2^2)])));

(*From De Roco and Hoover 1960 and Parsegian book, vdW interaction for the Hamaker limit between two equal rectangular sheets, dimensions a/b, the gap is L. In our notation, dimensions are (2hy/2Ly=2hx/2Lx) and the gap is q2. The final expressions are the same, and aligned with numerical integration. However, our sheetSheet expression is general, valid for sheets with arbitrary dimensions*)
thinSheetsDeRocco=((b*(L^2)+2*(a^2)*b)/(2*(L^4)*((((L^2)+(a^2)))^(1/2))))ArcTan[(b/((((L^2)+(a^2)))^(1/2)))]+((a*(L^2)+2*a*(b^2))/(2*(L^4)*((((L^2)+(b^2)))^(1/2))))ArcTan[(a/((((L^2)+(b^2)))^(1/2)))]-(b/(2*(L^3)))*ArcTan[(b/L)]-(a/(2*(L^3)))*ArcTan[(a/L)];
rs={b->2*Ly,a->2*hy,L->q2};
FullSimplify[(thinSheetsDeRocco/.rs)-sheetSheetSpec]

r={q2->1/1.,hx->2,hy->2,Ly->1,Lx->1};
thinSheetsDeRocco/.rs/.r
sheetSheetSpec/.r
NIntegrate[1/(q2^2+(y2-x2)^2+(((y3-x3))^2))^3/.r,{x2,-hx,hx}/.r,{y2,-hy,hy}/.r,{y3,-Ly,Ly}/.r,{x3,-Lx,Lx}/.r]

(*If we let hy->Infinity and Ly->Infinity, it is a plane figure/sheet (2bx/2Ly)-infinite plane potential for m=6, see Notebook 1*)
Limit[sheetSheetF[q2,hx,hy,Lx,Ly](*/.{hx->1/2,Lx->1/2}*),{hy->Infinity,Ly->Infinity}]
((8 hx Lx π (((-x1+y1))^(2-m)))/(-2+m))/.{m->6,-x1+y1->q2}