(* Solving the Schrodinger equation for the ground state
      in a torsion potential of S1 state of 1 - 
    phenylpyrrole molecule (JCP 109, 7185 (1998) *)
Clear["@"];
ndig = 24;
{b, v2, v4, v6, v8} = 
    SetPrecision[#, ndig] & /@ {0.4891, 1721, -132, -300, -28};
v[x_] := (phi = Sqrt[2 b]x;
      v2/2(1 - Cos[2 phi]) + v4/2(1 - Cos[4 phi]) + v6/2(1 - Cos[6 phi]) + 
        v8/2(1 - Cos[8 phi]));
en0 = -69.981428390952; (* calculated by FORTRAN program *)
xend = 0.9; (* Where the eigenfunction still does not diverge *)
{en0, xend} = SetPrecision[#, ndig] & /@ {en0, xend};
pltv = Plot[v[x], {x, -Pi/2 Sqrt[2 b], Pi/2 Sqrt[2 b]},
      PlotStyle -> {Thickness[0.01]}, DisplayFunction -> Identity];
plte = Plot[en0, {x, -xend, xend},
      PlotStyle -> {Thickness[0.005], RGBColor[0, 0, 1]}, 
      DisplayFunction -> Identity];
spsi = N[NDSolve[{-1/2psi''[x] + (v[x] - en0)psi[x] == 0, psi[0] == 1,
            psi'[0] == 0}, psi, {x, 0, -xend, xend}, 
          WorkingPrecision -> ndig][[1]]];
npsi = 2 NIntegrate[(psi[x]^2) /. spsi, {x, 0, xend}];
scpsi = 450;
pltpsi = Plot[scpsi (psi[x] /. spsi)/Sqrt[npsi], {x, -xend, xend},
      PlotStyle -> {Thickness[0.007], RGBColor[1, 0, 0]}, 
      DisplayFunction -> Identity];
rho[x1_, p1_] := 
    2/Pi NIntegrate[(psi[x1 + x2] /. spsi)(psi[x1 - x2] /. spsi)/npsi Cos[
            2 x2 p1],
        {x2, 0, 0, xend - Abs[x1]}];
scrho = 2000;
pltrho = Plot[scrho rho[x1, 0], {x1, -xend, xend},
      PlotStyle -> {Thickness[0.01], RGBColor[0, 1, 0]}, 
      DisplayFunction -> Identity];
$TextStyle = {FontFamily -> "Times", FontSize -> 18};
pltup = Show[pltv, plte, pltpsi, pltrho,
      AxesLabel -> {"\!\(q = \((2  B\()\^\(-\(1\/2\)\)\)\) \[CurlyPhi]\)", 
          "V, " <>
            ToString[scpsi] <> 
            "\[Times]\!\(\* StyleBox[\"\[Psi]\",FontColor->RGBColor[1, 0, 0]]\), " <>            
            ToString[scrho] <> 
            "\[Times]\!\(\* StyleBox[\"\[Rho]\",FontColor->RGBColor[0, 1, 0]]\)"},
      DisplayFunction -> Identity,
      PlotRange -> {{-0.95, 0.95}, {-130, 750}},
      Frame -> True,
      AspectRatio -> 1/3.9, ImageSize -> 600];
(* Finding Rho - min and Rho - max *)
{xmin00, ymin00} = {0, 5};
{dx, dy} = {0.1, 1};
{zmin0, fmin} = 
    FindMinimum[
      rho[x, y], {x, {xmin00 - dx/100, 
          xmin00 + dx/100}}, {y, {ymin00 - dy/100, ymin00 + dy/100}}];
{xmin0, ymin0} = {x, y} /. fmin;
{xmax00, ymax00} = {0, 0};
{zmax0, fmax} = 
    FindMinimum[-rho[x, y], {x, {xmax00 - dx/100, 
          xmax00 + dx/100}}, {y, {ymax00 - dy/100, ymax00 + dy/100}}];
{xmax0, ymax0} = {x, y} /. fmax;
zmax0 = -zmax0;
zscale = Max[Abs[zmin0], Abs[zmax0]];
(* Density plot, takes some time depending on value of npoints *)
{xmin, xmax} = {-0.52, 0.52};
{ymin, ymax} = {-12, 12};
npoints = 100;
z1max = 1.3;
z1white = 0.7;
cf = (z = zmin0 + #(zmax0 - zmin0);
        z1 = Abs[z/zscale] z1max;
        Hue[If[z > 0, 0, 0.65], 
          If[z1 <= 1, 0.25, 0.25(z1max - z1white z1)/(z1max - z1white)], 
          If[z1 <= 1, z1, 1]]) &;
plt1 = ContourPlot[rho[x, y]/zscale, {x, xmin, xmax}, {y, ymin, ymax},
      PlotPoints -> npoints + 1,
      Contours -> 200,
      PlotRange -> {zmin0, zmax0}/zscale,
      FrameLabel -> {"\!\(\* StyleBox[\"q\", FontSlant->\"Italic\"]\)", 
          "\!\(\* StyleBox[\"p\", FontSlant->\"Italic\"]\)"},
      ContourLines -> False,
      DisplayFunction -> Identity,
      AspectRatio -> 1/2,
      ImageSize -> 600,
      PlotRegion -> {{-0.165, 1}, {0, 1}},
      ColorFunction -> cf];
comb = GraphicsArray[{{pltup}, {pltdown}}, GraphicsSpacing -> -0.0];
gr = Show[Graphics[comb],
      ImageSize -> 900,
      DisplayFunction -> $DisplayFunction];