(*just a schematic representation of the conical intersection idea*)

E1=0.12(x+2)^2+0.1y^2 (*reactants' energy paraboloid *)
E2=0.12(x-2)^2+0.1y^2 (*products' energy paraboloid *)
Plot3D[{E1,E2},{x,-12,12},{y,-11,11},RegionFunction->Function[{x,y,z},z≤10],BoxRatios->Automatic,AxesLabel->{"x","y","Energy"},PlotPoints->21,PlotStyle->Directive[Opacity[0.90]]](*two paraboloids intersect*)
Show[
Plot3D[E2,{x,-11,11},{y,-11,11},RegionFunction->Function[{x,y,z},z≤10],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.90], RGBColor[1,0.5,0]],AxesLabel->{"x","y","Energy"}],
Plot3D[E1,{x,-11,11},{y,-11,11},RegionFunction->Function[{x,y,z},z≤10],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.90], RGBColor[0.5,1,0]],AxesLabel->{"x","y","Energy"}]
](*two paraboloids intersect; they have different colors*)
H12=Exp[-0.1x^2] (*the closer the paraboloids the stronger the coupling is; they are close for x=0*)
Delta1=2Sqrt[x^2+H12^2*Sqrt[x^2+0.01y^2]] (*the energy difference E1-E2 is proportional to x; the second term says that it is large for E1-E2=0 and the corresponding contribution increases when we go off  the origin *)
Leftparab=If[x>0,E1-Delta1,E1+Delta1]
Rightparab=If[x>0,E2+Delta1,E2-Delta1]

Show[
Plot3D[E2-Delta1,{x,-15,0},{y,-11,11},PlotRange->{{-15,10},{-11,11}},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.90], RGBColor[1,0.5,0]],AxesLabel->{"x","y","Energy"}],
Plot3D[E1+Delta1,{x,-15,0},{y,-11,11},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.70], RGBColor[0.5,1,0]]],
Plot3D[E1-Delta1,{x,0,10},{y,-11,11},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.90], RGBColor[0.5,1,0]],AxesLabel->{"x","y","Energy"}],
Plot3D[E2+Delta1,{x,0,10},{y,-11,11},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.90], RGBColor[1,0.5,0]]]
](*THIS IS THE MAIN PLOT FOR THE CONICAL INTERSECTION*)
E1b=E1*Exp[-0.12*x^2]-6Exp[-0.1(1.5(x+6)^2+0.5y^2)](*a minimum in the excited state*)
Show[
Plot3D[E2-Delta1,{x,-15,0},{y,-11,11},PlotRange->{{-15,10},{-11,11}},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.90], RGBColor[1,0.5,0],Specularity[Red,40]],AxesLabel->{"x","y","Energy"}],
Plot3D[E1b+Delta1,{x,-15,0},{y,-11,11},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.90], RGBColor[0.5,1,0],Specularity[Red,40]]],
Plot3D[E1-Delta1,{x,0,10},{y,-11,11},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.90], RGBColor[0.5,1,0],Specularity[Red,40]],AxesLabel->{"x","y","Energy"}],
Plot3D[E2+Delta1,{x,0,10},{y,-11,11},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.90], RGBColor[1,0.5,0],Specularity[Red,40]]]
](*THE CONICAL INTERSECTION WITH A MINIMUM IN THE EXCITED STATE*)
Plot3D[{E1-Delta1,E2+Delta1},{x,-20,20},{y,-11,11},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21,PlotStyle->Directive[Opacity[0.90]],AxesLabel->{"x","y","Energy"}](*conical intersection model, just an interesting view, in principle not related to the previous pictures*)
F1=0.12(x+2)^2+0.1y^2 (*reactants' energy paraboloid *)
F2=0.24(x-2)^2+0.05y^2 (*products' energy paraboloid *)
Show[
Plot3D[F2,{x,-8,8},{y,-11,11},RegionFunction->Function[{x,y,z},z≤10],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.90], RGBColor[1,0.5,0]],AxesLabel->{"x","y","Energy"}],
Plot3D[F1,{x,-8,8},{y,-11,11},RegionFunction->Function[{x,y,z},z≤10],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.90], RGBColor[0.5,1,0]],AxesLabel->{"x","y","Energy"}]
](*two paraboloids intersect; they have different colors*)
E1b=E1*Exp[-0.12*x^2]-4Exp[-0.1(1.5(x+6)^2+0.5y^2)](*a minimum in the excited state*)
Show[
Plot3D[E2-Delta1,{x,-15.,0},{y,-5.0,5.5},PlotRange->{{-15.,5},{-4.0,5.5}},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.80], RGBColor[0.1,1.,0],Specularity[Yellow,50]],AxesLabel->{"x","y","Energy"}],
Plot3D[E1b+0.5*Delta1,{x,-10.,0},{y,-5.0,5.5},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.80], RGBColor[1.,0.1,0],Specularity[Yellow,50]]],
Plot3D[E1-Delta1+Exp[-0.2(x-5)^2-y^2],{x,0,10},{y,-5.0,5.5},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.80], RGBColor[1.,0.1,0],Specularity[Yellow,50]],AxesLabel->{"x","y","Energy"}],
Plot3D[E2+0.5*Delta1,{x,0,10},{y,-5.0,5.5},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.80], RGBColor[0.1,1.,0],Specularity[Yellow,50]]],Graphics3D[Sphere[{-6.5,0,-3.8},0.5]],Graphics3D[Sphere[{-6.5,0,4},0.5]] ,Graphics3D[Sphere[{-5.2,0,2.2},0.5]],Graphics3D[Sphere[{-3.0,0,2.8},0.5]],Graphics3D[Sphere[{0,0,1.},0.5]],Graphics3D[Sphere[{4,2,-2.7},0.5]],Graphics3D[Sphere[{4,-2,-2.7},0.5]]
](*THE CONICAL INTERSECTION WITH A MINIMUM IN THE EXCITED STATE, MORE FLAT*)
E1b=E1*Exp[-0.12*x^2]-4.5Exp[-0.1(1.5(x+6)^2+0.5y^2)](*a minimum in the excited state*)
l=Line[{{-6.5,0,-3.8},{-6.5,0,4}}];
Show[
Plot3D[E2-Delta1,{x,-15.,0},{y,-5.0,5.5},PlotRange->{{-15.,5},{-5.5,5.5}},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.80], RGBColor[0.2,1.,0]],AxesLabel->{"x","y","Energy"}],
Plot3D[E1b+0.5*Delta1,{x,-10.,0},{y,-5.0,5.5},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.80], RGBColor[1.,0.2,0]]],
Plot3D[E1-Delta1+Exp[-0.2(x-5)^2-y^2],{x,0,10},{y,-5.0,5.5},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.80], RGBColor[1.,0.2,0]],AxesLabel->{"x","y","Energy"}],
Plot3D[E2+0.5*Delta1,{x,0,10},{y,-5.0,5.5},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.80], RGBColor[0.2,1.,0]]],Graphics3D[{Sphere[{-6.5,0,-3.8},0.5],Text[1,{-6.3,0,-3.3}]}],Graphics3D[{Sphere[{-6.5,0,4},0.5],Text[2,{-6.5,0,4}]}],Graphics3D[{Sphere[{-5.2,0,1.7},0.5],Text[3,{-5.2,0,1.7}]}],Graphics3D[{Sphere[{-3,0,2.5},0.5],Text[4,{-3,0,2.5}]}],Graphics3D[{Sphere[{0,0,1},0.5],Text[5,{0,0.2,0.6}]}],Graphics3D[{Sphere[{4,2,-3.2},0.5],Text[6b,{4,2,-3.0}]}],Graphics3D[{Sphere[{4,-2,-2.7},0.5],Text[6a,{4,-2,-2.7}]}],Graphics3D[{Thickness[0.015],Blue,l}],
Graphics3D[{Blue,Thickness[0.008],Line[{{-6.5,0.,3.5},{-6.1,0.,2.3}}]}],
Graphics3D[{Blue,Thickness[0.008],Line[{{-6.5,0.,3.5},{-6.9,0.,2.3}}]}],
Graphics3D[{Black,Text[StyleForm["PS",FontSize->30,FontWeight->"Bold"],{1.,0.,-0.3}]}],
Graphics3D[{Text[StyleForm["adiabat. stan podstawowy",FontColor->Black,FontSize->30,FontWeight->"Bold"],{-5.5,0.,-7.}]}],Graphics3D[{Text[StyleForm["adiabat. stan wzbudzony",FontColor->Black,FontSize->30,FontWeight->"Bold"],{-3.5,0.,8.5}]}],Graphics3D[{Text[StyleForm["obraz produktów",FontColor->Darker[Red],FontSize->20,FontWeight->"Bold"],{-3.5,0.,7.}]}],
Graphics3D[{Text[StyleForm["produkty",FontColor->Darker[Red],FontSize->20,FontWeight->"Bold"],{3.5,2.,0.}]}],Graphics3D[{Text[StyleForm["obraz substratów",FontColor->Darker[Green],FontSize->20,FontWeight->"Bold"],{3.5,0.,9.}]}],
Graphics3D[{Text[StyleForm["substraty",FontColor->Darker[Green],FontSize->20,FontWeight->"Bold"],{-18.,2.,-5.}]}],
Graphics3D[{Blue,Text[StyleForm["hν",FontSize->30,FontWeight->"Bold"],{-8.,0.,-0.3}]}],
Graphics3D[{Blue,Text[StyleForm["FC",FontSize->30,FontWeight->"Bold"],{-7.,0.,5.}]}],
Graphics3D[{Darker[Green],Text[StyleForm["S",FontSize->30,FontWeight->"Bold"],{-6.5,0.,-5.5}]}],
Graphics3D[{Darker[Red],Text[StyleForm["(Subscript[P, a])",FontSize->30,FontWeight->"Bold"],{4,-2,-4.0}]}],
Graphics3D[{Darker[Red],Text[StyleForm["(Subscript[P, b])",FontSize->30,FontWeight->"Bold"],{5,3,-4.0}]}]
](*THE CONICAL INTERSECTION WITH A MINIMUM IN THE EXCITED STATE, DEEPER*)
Show[
Plot3D[E2-Delta1,{x,-15.,0},{y,-5.0,15},PlotRange->{{-15.,5},{-5.5,15}},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.80], RGBColor[0.2,1.,0]],AxesLabel->{"x","y","Energy"}],
Plot3D[E1b+0.5*Delta1,{x,-10.,0},{y,-5.0,5},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->21 ,PlotStyle->Directive[Opacity[0.80], RGBColor[1.,0.2,0]]],
Plot3D[E1-Delta1+Exp[-0.2(x-5)^2-y^2],{x,0,10},{y,-5.0,15},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.80], RGBColor[1.,0.2,0]],AxesLabel->{"x","y","Energy"}],
Plot3D[E2+0.5*Delta1,{x,0,10},{y,-5.0,5.5},RegionFunction->Function[{x,y,z},z≤8],BoxRatios->Automatic,PlotPoints->61 ,PlotStyle->Directive[Opacity[0.80], RGBColor[0.2,1.,0]]],Graphics3D[{Sphere[{-10,0,0.},0.5],Text[1,{-6.3,0,-3.3}]}]
](*bez kul*)