(*  Paper: Quapp / Bofill, J. Phys Chem B ... 2024 ... *)
(* Publication Date: April 18, 2024
   https://doi.org/10.1021/acs.jpcb.4c00468 *)
(* First we calculate the PES and Newton trajectories (NT) on it *)
(* main direction (L+d) for dashed NT, comparison with (-1,1) for black NT  *)
(* blue is the corresponding singular NT, red the NT to direction (1,0) *)


(* first molecule LSelA108A + 2GSP6 : no.4 in Table of BB *)

ph0=0.58*Pi ;
d0=0.33;
k0=266.0;
D0=201.0;
cc=0.0;
ph1=0.97*Pi ;
sig=0.17*Pi ;

PHI[L_]=2.0*ArcSin[L/5.6]
De[ph_]=D0*Exp[-(ph-ph1)^2/(2.*sig^2)]+cc
Vphi[L_]=0.5*k0*(PHI[L]-ph0)^2
Be[d_,ph_]=De[ph]*((1.0-Exp[-d/d0])^2-1.)
pes[L_,d_,F_]=Vphi[L]+Be[d,PHI[L]]-(L+d)*F

conBarLSelA=ContourPlot[pes[L,d,0],{L,4.25,5.6},{d,-0.2,0.75},ContourShading -> False, Contours ->50,PlotRange->All,PlotPoints -> 53, ContourStyle->Thickness[0.00215],FrameLabel->{"L","d"},AspectRatio->1,FormatType -> TraditionalForm,FrameStyle->FontSize->16]

(* derivations, here for traditional reasons, the second variable has no.3 *) 
 p3[L_,d_]=D[pes[L,d,0],d];
 p1[L_,d_]=D[pes[L,d,0],L];
 p11[L_,d_]=D[p1[L,d],L];
 p33[L_,d_]=D[p3[L,d],d];
 p31[L_,d_]=D[p3[L,d],L];

(* For NTs the gradient has to be parallel to the search direction. 
In 2D we can plot the NT by the zero line of the interchanged coefficients of the gradient componenets NT1 is to (1,0),   NT2 is to direction (-1,1) *)

NewT1=ContourPlot[p3[L,d],{L,4.25,5.6},{d,-0.2,0.15},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Red}]

NewT2=ContourPlot[p1[L,d]+p3[L,d],{L,4.25,5.6},{d,-0.22,0.75},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Black}]

(* In 2D also the determinant of the Hessian matrix is simple *)
DetH[L_,d_]=Simplify[p11[L,d]*p33[L,d]-p31[L,d]*p31[L,d]];
(* We search for the line Det(H)=0 on the PES  *)
deter1=ContourPlot[Evaluate[DetH[L,d],{L,4.25,5.6},{d,-0.2,0.75},ContourShading->False,PlotPoints->40,Contours->{0.0},ContourStyle->{Thickness[0.0046],Green}]]

(* NT11 to direction (1,1), that one of main interest *)

NewT11=ContourPlot[p3[L,d]-1.0*p1[L,d],{L,4.25,5.6},{d,-0.2,0.75},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Dashed,Thickness[0.0085]}]

(* next are data for the figures *)

spArr=Graphics[Arrow[{{4.537,0.65},{4.537,0.73}}]]
textLSelA=Show[Graphics[{Text[Style["R"],{5.5540,-0.05}],
Text[Style["TS"],{4.615,0.7}],
Text[Style["LSelA108H"],{5.29,0.71}],
Text[Style["+2GSP6"],{5.293,0.635}]},
PlotRange->{{L,1.5,5.5},{d,-.25,0.5}}]];
BiLSelA=Show[conBarLSelA,NewT2,deter1,NewT1,NewTNT11,spArr,textLSelA,TextStyle->FontSize->19]

Export["BildLSelA108H.pdf",BiLSelA,ImageResolution->400,ImageSize->Automatic]

(*  NEXT EXAMPLE  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*)

(*  LSelN138G + PSGL1   :  no.2 in BB tabelle *)
ph0=0.58*Pi ;
d0=0.29;
k0=220.0;
D0=210.0;
cc=0.0;
ph1=Pi ;
sig=0.12*Pi;

PHI[L_]=2.0*ArcSin[L/5.6]
De[ph_]=D0*Exp[-(ph-ph1)^2/(2.*sig^2)]+cc
Vphi[L_]=0.5*k0*(PHI[L]-ph0)^2
Be[d_,ph_]=De[ph]*((1.0-Exp[-d/d0])^2-1.)
pes[L_,d_,F_]=Vphi[L]+Be[d,PHI[L]]-(L+d)*F

conBarLSelN=ContourPlot[pes[L,d,0],{L,4.25,5.6},{d,-0.2,0.65},ContourShading -> False, Contours ->50,PlotRange->All,PlotPoints -> 53, ContourStyle->Thickness[0.00215],FrameLabel->{"L","d"},AspectRatio->1,FormatType -> TraditionalForm,FrameStyle->FontSize->16]

p3[L_,d_]=D[pes[L,d,0],d];
 p1[L_,d_]=D[pes[L,d,0],L];
 p11[L_,d_]=D[p1[L,d],L];
 p33[L_,d_]=D[p3[L,d],d];
 p31[L_,d_]=D[p3[L,d],L];
NewT1=ContourPlot[p3[L,d],{L,4.25,5.6},{d,-0.2,0.15},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Red}]

NewT2=ContourPlot[p1[L,d]+p3[L,d],{L,4.25,5.6},{d,-0.22,0.65},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Black}]

DetH[L_,d_]=Simplify[p11[L,d]*p33[L,d]-p31[L,d]*p31[L,d]];
deter1=ContourPlot[Evaluate[DetH[L,d],{L,4.25,5.6},{d,-0.2,0.75},ContourShading->False,PlotPoints->40,Contours->{0.0},ContourStyle->{Thickness[0.0046],Green}]]
NewTNT11=ContourPlot[p3[L,d]-1.0*p1[L,d],{L,4.25,5.6},{d,-0.2,0.65},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Dashed,Thickness[0.0085]}]

(* The singular NT through the VRI point must be determined by hand, it means one has to probe the coefficients for the searched intersection of the two branches of the NT , drawn in blue *)

NewTSi=ContourPlot[-0.296*p1[L,d]+1.445*p3[L,d],{L,4.25,5.6},{d,-0.2,0.55},ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Blue}]

textLSelN=Show[Graphics[{Text[Style["R"],{5.5540,-0.05}],
Text[Style["TS"],{5.05,-0.05}],
Text[Style["LSelN138G"],{5.29,0.626}],
Text[Style["+ PSGL1"],{5.273,0.57}],
Text[Style["VRI"],{4.62,0.12}],
Text[Style["Int"],{4.49,-0.05}]},
PlotRange->{{L,1.5,5.5},{d,-.25,0.5}}]];
BiLSelN=Show[conBarLSelN,NewT2,deter1,NewT1,NewTSi,NewTNT11,textLSelN,TextStyle->FontSize->19]

Export["BildLSelN138G.pdf",BiLSelN,ImageResolution->400,ImageSize->Automatic]

(*  L-selectin + 2GSP6  no.3 in BB Table *)
ph0=0.58*Pi ;
d0=0.33;
k0=266.0;
D0=217.0;
cc=0.0;
ph1=0.98*Pi ;
sig=0.09*Pi;

PHI[L_]=2.0*ArcSin[L/5.6]
De[ph_]=D0*Exp[-(ph-ph1)^2/(2.*sig^2)]+cc
Vphi[L_]=0.5*k0*(PHI[L]-ph0)^2
Be[d_,ph_]=De[ph]*((1.0-Exp[-d/d0])^2-1.)
pes[L_,d_,F_]=Vphi[L]+Be[d,PHI[L]]-(L+d)*F

conBarL2G=ContourPlot[pes[L,d,0],{L,4.25,5.6},{d,-0.2,0.75},ContourShading -> False, Contours ->50,PlotRange->All,PlotPoints -> 53, ContourStyle->Thickness[0.00215],FrameLabel->{"L","d"},AspectRatio->1,FormatType -> TraditionalForm,FrameStyle->FontSize->16]

p3[L_,d_]=D[pes[L,d,0],d];
 p1[L_,d_]=D[pes[L,d,0],L];
 p11[L_,d_]=D[p1[L,d],L];
 p33[L_,d_]=D[p3[L,d],d];
 p31[L_,d_]=D[p3[L,d],L];

NewT2=ContourPlot[p1[L,d]+p3[L,d],{L,4.25,5.6},{d,-0.22,0.75},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Black}]

NewT1=ContourPlot[p3[L,d],{L,4.25,5.6},{d,-0.2,0.15},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Red}]

DetH[L_,d_]=Simplify[p11[L,d]*p33[L,d]-p31[L,d]*p31[L,d]];
deter1=ContourPlot[Evaluate[DetH[L,d],{L,4.25,5.6},{d,-0.2,0.75},ContourShading->False,PlotPoints->40,Contours->{0.0},ContourStyle->{Thickness[0.0046],Green}]]

NewTNT11=ContourPlot[p3[L,d]-1.0*p1[L,d],{L,4.25,5.6},{d,-0.2,0.75},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Dashed,Thickness[0.0085]}]

p3[4.25,0.1273]

NewTSi=ContourPlot[-0.04040*p1[L,d]+9.133*p3[L,d],{L,3.25,5.56},{d,-0.2,0.75},ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Blue}]

spArr=Graphics[Arrow[{{4.537,0.65},{4.537,0.73}}]]
textLG2=Show[Graphics[{Text[Style["R"],{5.5540,-0.05}],
Text[Style["TS"],{4.615,0.7}],
Text[Style["L-selectin"],{5.29,0.71}],
Text[Style["+ 2GSP6"],{5.293,0.635}]},
PlotRange->{{L,1.5,5.5},{d,-.25,0.5}}]];
BiL2G=Show[conBarL2G,deter1,NewTSi,NewT1,NewTNT11,spArr,textLG2,TextStyle->FontSize->19]

Export["BildLselec2GSP6.pdf",BiL2G,ImageResolution->400,ImageSize->Automatic]

(*  L-selectin + PSGL1  vorher falsch k0=299.0  XXXXXXXXXXXXXXXXX Haupteil !!!  *)
(*  No 1 in Tabelle in PNAS ist auch gleich !!!!! D.h. L-selectin alles wie vorher *) 
ph0=0.58*Pi ;
d0=0.29;
k0=266.0;
D0=237.0;
cc=0.0;
ph1=Pi ;
sig=0.12*Pi;

(* phi(L) indirekt gegeben !! W=2.8  *)
PHI[L_]=2.0*ArcSin[L/5.6]
Plot[PHI[L],{L,4.5,5.6}]

theta0=0.58 *Pi
L0=5.6* Sin[theta0/2.]

(* phi(L) indirekt gegeben !! W=2.8  *)
PHI[L_]=2.0*ArcSin[L/5.6]
Plot[PHI[L],{L,4.5,5.6}]
De[ph_]=D0*Exp[-(ph-ph1)^2/(2.*sig^2)]+cc
Vphi[L_]=0.5*k0*(PHI[L]-ph0)^2
Be[d_,ph_]=De[ph]*((1.0-Exp[-d/d0])^2-1.)
pes[L_,d_,F_]=Vphi[L]+Be[d,PHI[L]]-(L+d)*F

Plot[De[ph],{ph,3.5,5.6}]

Plot[Vphi[L],{L,4.05,5.8}]

Plot[(1.0-Exp[-d/d0])^2-1.,{d,-0.2,1.0}]

conBarL=ContourPlot[pes[L,d,0],{L,4.25,5.6},{d,-0.2,0.55},ContourShading -> False, Contours ->50,PlotRange->All,PlotPoints -> 53, ContourStyle->Thickness[0.00215],FrameLabel->{"L","d"},AspectRatio->1,FormatType -> TraditionalForm,FrameStyle->FontSize->16]

Plot3D[pes[L,d,0],{L,4.25,5.6},{d,-0.2,0.45},PlotRange->All]

p3[L_,d_]=D[pes[L,d,0],d];
 p1[L_,d_]=D[pes[L,d,0],L];
 p11[L_,d_]=D[p1[L,d],L];
 p33[L_,d_]=D[p3[L,d],d];
 p31[L_,d_]=D[p3[L,d],L];

FindRoot[{p1[L,d],p3[L,d]},{{L,4.46},{d,0.0}}]

{L->4.461915723575546`,d->0.`}   (*  Intermediate *)

pes[4.461915723575546,0.0,0]

FindRoot[{p1[L,d],p3[L,d]},{{L,5.15},{d,0.0}}]

{L->5.150914868876448`,d->0.`}  (*  SP *)

pes[5.15091,0.0,0]

FindRoot[{p1[L,d],p3[L,d]},{{L,5.59},{d,0.0}}]

{L->5.570252566046982`,d->0.`}    (*  Min R *)

pes[5.57025256,0.0,0]

Bar0=10.956699577153884+39.23677075785412

(*  Control of kind of stat points *)

H[L_,d_]={{p11[L,d],p31[L,d]},{p31[L,d],p33[L,d]}};

Eigensystem[H[4.461915723,0.0]]

Eigensystem[H[5.15091,0.0]]

Eigensystem[H[5.57025256,0.0]]

NewT2=ContourPlot[p1[L,d]+p3[L,d],{L,4.25,5.6},{d,-0.22,0.55},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Black}]

(*  vorher NewT1 *)

NewTed=ContourPlot[p3[L,d],{L,4.25,5.6},{d,-0.2,0.15},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Thickness[0.005],Red}]  ;

DetH[L_,d_]=Simplify[p11[L,d]*p33[L,d]-p31[L,d]*p31[L,d]];
deter1=ContourPlot[Evaluate[DetH[L,d],{L,4.25,5.6},{d,-0.2,0.45},ContourShading->False,PlotPoints->40,Contours->{0.0},ContourStyle->{Thickness[0.0046],Green}]]

NewT11=ContourPlot[p3[L,d]-1.0*p1[L,d],{L,4.25,5.6},{d,-0.2,0.55},
       ContourShading->False,PlotPoints->77,Contours->{0.0},
       ContourStyle->{Dashed,Thickness[0.0085]}]

(*  VRI point ~ (4.64, 0.1675)   *)
p1[4.64, 0.1675]

p3[4.64, 0.1675]

NewTSi=ContourPlot[-0.27963771199820813*p1[L,d]+1.5131569585841813*p3[L,d],{L,4.25,5.56},{d,-0.2,0.55},ContourShading->False,PlotPoints->77,Contours->{0.0},ContourStyle->{Thickness[0.005],Blue}]

spArr=Graphics[Arrow[{{4.537,0.5},{4.537,0.553}}]]

text=Show[Graphics[{Text[Style["R"],{5.5540,-0.05}],
Text[Style["TS"],{5.15,-0.05}],
Text[Style["TS"],{4.615,0.53}],
Text[Style["L-selectin"],{5.29,0.52}],
Text[Style["+ PSGL1"],{5.273,0.45}],
Text[Style["VRI"],{4.62,0.12}],
Text[Style["Int"],{4.4,-0.05}]},
PlotRange->{{L,1.5,5.5},{d,-.25,0.5}}]];

BiL=Show[conBarL,NewT1,deter1,NewT1,NewTSi,NewTNT11,spArr,text,TextStyle->FontSize->19]

Export["BildLselectin.pdf",BiL,ImageResolution->400,ImageSize->Automatic]

(*   xxxxxxxxxxxxxxxxx XXXXXXXXXXXXXXXXXXXX  xxxxxxxxxxxxxx  *)
(*   !!! Calculation of the catch bond barriers !!! *)

p1F[L_,d_,F_]=D[pes[L,d,F],L];
p3F[L_,d_,F_]=D[pes[L,d,F],d];

Lo=5.560;
do=0.0;
Minweg=Table[{mini=FindRoot[p1F[L,d,F]==0&&p3F[L,d,F]==0,{L,Lo},{d,do}];{Lo,do}={L,d} /.mini;pes[Lo,do,F]},{F,0.0,100,5}]
(

(*  +++++++++++++++++ SP   *)

(*  TS  *)
Lo=5.160;
do=0.01;
MinSP=Table[{saddle=FindRoot[p1F[L,d,F]==0&&p3F[L,d,F]==0,{L,Lo},{d,do}];{Lo,do}={L,d} /.saddle;pes[Lo,do,F]},{F,0.0,100,5}]

barrier=Flatten[MinSP-Minweg]

BiListe=ListLinePlot[barrier,AxesLabel->{"F/5","Barrier"},AxesStyle->FontSize->16]

(* here we obtain a nice catch bond behavior *)