Skip to content

DVR

Jump to bottom
b3m2a1 edited this page Jul 3, 2018 · 2 revisions

A fully-featured discrete variable representation package is built into ChemTools. It provides a large base of built-in DVR utilities and also may be extended at will.

The heart of the package is the ChemDVRObject :

<<ChemTools`DVR`

We can load a 1D DVR on a Cartesian grid like so:

dvr = ChemDVRObject["Cartesian1DDVR"]
(*Out:*)

dvr-373887066624419482

As with large portions of the package it is implemented in an object-oriented way to allow for more flexible use and better encapsulation.

‌

Basic DVR

Test DVR functionality is built in directly. Calling a DVR object without any arguments will calculate and plot the DVR wavefunctions. We can supply arguments to select how these display and how many:

dvr[
  "PlotDisplayMode"->Show, 
  "WavefunctionSelection"->5,
  "WavefunctionShifting"->"Energy",
  "WavefunctionScaling"->.5,
  "WavefunctionClipping"->None
  ]
(*Out:*)

dvr-8369702441338047581

If we want to set arguments that apply every time we call the DVR we can pass them via the "RuntimeOptions" :

dvr["RuntimeOptions"]=
  {
    "PlotDisplayMode"->Show, 
    "WavefunctionSelection"->5,
    "WavefunctionShifting"->"Energy",
    "WavefunctionScaling"->.5,
    "WavefunctionClipping"->None
    };

Now when we call our DVR again it will use these automatically;

dvr[]
(*Out:*)

dvr-795583022599405158

Potentials

We can also use an arbitrary potential. There are a number of special potentials we can use, too:

dvr["PotentialFunction"->{"MultiWellPolynomial", "TurningPoints"->{-5, 0, 5}}]
(*Out:*)

dvr-2843937171153480684

Or we can supply a function directly:


dvr["PotentialFunction"->(10(1-E-1/15*#)^2&), PlotRange->{0, 5}]


(*Out:*)

dvr-4074255523434672647

One incredibly useful case of this is in using an InterpolatingFunction returned via Interpolation . This allows for potentials extracted from potential energy scans using the electronic structure package of one's choice to be used directly as source potentials for a DVR.

Another case where this can be useful is when the potential function is expensive. If that's the case we can use FunctionInterpolation to construct a fast approximation for with varying grids and discretizations and things. Here's an example of that.

First we'll evaluate our DVR over a wonky, slow function, expressed in terms of integrals of the AiryAi function. We'll do this with an excessive number of points, just to highlight the slowness.

{time, plot} = 
  AbsoluteTiming@
    dvr[
      "PotentialFunction"->
        Function[x, 5*NIntegrate[AiryAi[(t-10)/2], {t, -10, x}]//Quiet],
      "Points"->666, 
      PlotRange->{-2, 3}, 
      PlotRangePadding->None
      ];
time

(*Out:*)

8.30746`

plot
(*Out:*)

dvr-4769566160543927881

Now instead we'll create an interpolation over it:

pot = Quiet@FunctionInterpolation[5*NIntegrate[AiryAi[(t-10)/2], {t, -10, x}], {x, -10, 10}]
(*Out:*)

dvr-5242833578834456675

And now we can use this in our DVR without it taking forever:

{time, plot} = 
  AbsoluteTiming@
    dvr[
      "PotentialFunction"->
        pot,
      "Points"->666, 
      PlotRange->{-2, 3}, 
      PlotRangePadding->None
      ];
time
(*Out:*)

4.071804`

plot
(*Out:*)

dvr-2124589457410979137

And we see we gained a full four seconds over the slow, integration at each point. As we play with the number of points and the dimensionality of the problem the savings can be even more dramatic.

Masses

We can also play with the mass of our system:


dvr["PotentialFunction"->(10(1-E-1/15*#)^2&), PlotRange->{0, 5}, "Mass"->.25]


(*Out:*)

dvr-4586967124960494851

Probability Density

Plot the probability density instead of the wavefunction:


dvr[
  "PotentialFunction"->(10(1-E-1/15*#)^2&), 
  PlotRange->{0, 5}, "Mass"->.25,
  "PlotProbabilityDensity"->True
  ]


(*Out:*)

dvr-3741032298068959063

Range

Use a larger grid range


dvr[
  "Range"->{-20, 20},
  "PotentialFunction"->(10(1-E-1/15*#)^2&), 
  PlotRange->{0, 5}, 
  "Mass"->.25,
  "PlotProbabilityDensity"->True
  ]


(*Out:*)

dvr-4299460644312072876

Points

Use fewer points on the grid


dvr[
  "Range"->{-20, 20},
  "Points"->25,
  "PotentialFunction"->(10(1-E-1/15*#)^2&),
  PlotRange->{0, 5}, 
  "Mass"->.25,
  "PlotProbabilityDensity"->True
  ]


(*Out:*)

dvr-1512047220194205403

One thing worth noting is that the object will store any passed values for the "Points" and the "Range" . It does this to make further a bit cleaner, although this can come back to bite you if you're not careful about it.

‌

Direct Passing

If we already have the "KineticEnergy" , "PotentialEnergy" , we can pass these directly.

As an example, we'll pre-calculate the kinetic energy and potential energy matrices:

ke = dvr["KineticEnergy"]; // RepeatedTiming
(*Out:*)

{0.0482854500000000075`2.,Null}

pe = dvr["PotentialEnergy"]; // RepeatedTiming
(*Out:*)

{0.0006658179271708684`2.,Null}

MatrixPlot/@{ke, pe}//GraphicsRow
(*Out:*)

dvr-990246090326943613

Then we'll compare how quickly we can compute the wavefunctions by directly passing this or not:

wfns = dvr["Wavefunctions"]; // RepeatedTiming
(*Out:*)

{0.0850724999999999953`2.,Null}

wfns2 = dvr["Wavefunctions", "KineticEnergy"->ke, "PotentialEnergy"->pe ]; // RepeatedTiming
(*Out:*)

{0.0305115937499999997`2.,Null}

And these agree up to inherent numerical instabilities/phases:

Round[wfns[[1]], 10^-8]==Round[wfns2[[1]], 10^-8]
(*Out:*)

True

We can also notice that the speed up is almost exactly that of the time spent computing the matrices in the earlier case.

For things that only need the wavefunctions and the grid (e.g. the plotting routines) we can pass these directly, too. In this case we'll allow the "Grid" to be computed instead of passing it.

dvr[ "Wavefunctions"->wfns ]
(*Out:*)

dvr-691257566539814854

‌

DVR Properties

There are many parameters and things that may be requested of these objects, which we can see via:

dvr["Properties"]
(*Out:*)

{"Name","Range","Points","Defaults","File","Class","UUID","Grid","KineticEnergy","PotentialEnergy","GridPotentialEnergy","Wavefunctions","View","GridWavefunctions","InterpolatingWavefunctions","ExpectationValues","OperatorMatrix","FormatGrid"}

Hamiltonian

We can access the Hamiltonian directly:

MatrixPlot[dvr["Hamiltonian"]]
(*Out:*)

dvr-7263548070700961953

GridPotentialEnergy

Or see the potential energy plotted on the DVR grid:

ListLinePlot[dvr["GridPotentialEnergy"]]
(*Out:*)

dvr-3626467297131689715

Energies

We can access just the energies

dvr["Energies"]
(*Out:*)

{0.4999999999997726`,1.4999999999997726`,2.500000000000682`,3.4999999999995453`,4.499999999999545`,5.5000000000009095`,6.499999999999773`,7.499999999999773`,8.500000000000114`,9.500000000000114`}

(Note the harmonic oscillator energy pattern—a comforting result as the default for this 1D DVR is a harmonic oscillator potential)

These will track directly with ℏ

dvr["Energies", "HBar"->2]
(*Out:*)

{0.9999999999882903`,3.0000000000070486`,4.999999999983288`,6.999999999993406`,9.000000000003524`,11.00000000000125`,13.000000000020918`,14.999999999980673`,17.000000000005002`,19.00000000000989`}

Expectation Values

We can supply arbitrary functions for which to calculate the expectation value. These functions should take the grid or a gridpoint as the first argument. If the value can be computed listably it will.

The second argument is optional and will be the value of the wavefunction on the grid. Here's a way to calculate the expectation value of the coordinate:

dvr["ExpectationValues", #&, "WavefunctionSelection"->10]//Chop
(*Out:*)

{0,0,0,0,0,0,0,0,0,0}

Operator Matrix Calculations

We can also compute matrix elements for an arbitrary operator (i.e. the elements of the matrix for which the expectation values are the diagonal elements). Here's a way to calculate the operator matrix for the coordinate, which is in turn just the overlaps between the wavefunctions:

MatrixPlot[dvr["OperatorMatrix", #&, "WavefunctionSelection"->10]//Chop]
(*Out:*)

dvr-3889900728127672864

‌

Multidimensional DVR

Custom DVR classes for various 2D, 3D and a single 4D DVR are implemented. The most basic, a 2D DVR on a Cartesian grid, is accessible like so:

dvr2D = ChemDVRObject["Cartesian2DDVR"];

For good measure we'll use a sum of a 2D polynomial and a Bessel function:

dvr2D[
  "GridPotentialEnergy",
  "Points"->{30, 30},
  "PotentialFunction"->
    {"MultiWellPolynomial", "TurningPoints"->{-4, 0, 4}, "Depth"->10^-1}**
      Function[.5*BesselJ[0, #/1.5]]
  ]//ListPlot3D[#, AxesLabel->{"x", "y"}]&
(*Out:*)

dvr-2076129964559736229

dvr2D[
  "Points"->{30, 30},
  "PotentialFunction"->
    {"MultiWellPolynomial", "TurningPoints"->{-4, 0, 4}, "Depth"->10^-1}**
      Function[.5*BesselJ[0, #/1.5]],
  "PlotDisplayMode"->List,
  "ShowPotential"->True,
  "NumberOfWavefunctions"->3,
  "WavefunctionClipping"->None,
  "WavefunctionRescaling"->None,
  "WavefunctionScaling"->1
  ]//GraphicsRow[#, ImageSize->600]&
(*Out:*)

dvr-3723850491367664870

That faint gray outline is the potential drawn on top of the wavefunctions.

We could also plot these as density plots:

dvr2D[
  "Points"->{30, 30},
  "PotentialFunction"->
    {"MultiWellPolynomial", "TurningPoints"->{-4, 0, 4}, "Depth"->10^-1}**
      Function[.5*BesselJ[0, #/1.5]],
  "NumberOfWavefunctions"->3,
  "PlotDisplayMode"->List,
  ColorFunction->"Rainbow",
  "PlotMode"->"Density",
  "WavefunctionClipping"->None
  ]//GraphicsRow[#, ImageSize->600]&
(*Out:*)

dvr-3353478607919207490

Potential Optimized DVR

If degrees of freedom are weakly coupled in a system it is possible to be more efficient in how we compute the DVR by optimizing the 2D grid based on the 1D solutions. We can do this with the "PotentialOptimizationOptions" parameter passed to the DVR.

Here's a new potential that looks a lot like the previous one, but now with a term coupling the degrees of freedom:

potCoup = 
  Function[
    .5*BesselJ[0, #[[1]]/1.5]+
    .5*BesselJ[0, #[[2]]/1.5]+
      .05*Sin[(#[[1]]*#[[2]])/15]
    ];
Plot3D[potCoup[{x, y}], {x, -10, 10}, {y, -10, 10}]
(*Out:*)

dvr-6550985844468222297

We can get our DVR energies here on a refined grid:

{time, plots}=
  dvr2D[
    "Points"->{80, 80},
    "PotentialFunction"->potCoup,
    "NumberOfWavefunctions"->3,
    "PlotDisplayMode"->List,
    ColorFunction->"Rainbow",
    "PlotMode"->"Density",
    "WavefunctionClipping"->None
    ]//GraphicsRow[#1,ImageSize->600]&//AbsoluteTiming;
time
(*Out:*)

30.244599`

plots
(*Out:*)

dvr-7788038674918945049

But we can also do this much faster by pre-optimizing in each degree of freedom. To do this we specify a higher precision grid, and then pass our potential in each coordinate in the "PotentialOptimizationOptions" . In that we specify how many optimized grid points to use via the "OptimizedBasisSize" option

{time2, plots2}=
  dvr2D[
    "Points"->{125, 125},
    "PotentialFunction"->potCoup,
    "NumberOfWavefunctions"->3,
    "PotentialOptimizationOptions"->
      {
        "PotentialFunction"->
          {
            Function[.5*BesselJ[0, #/1.5]],
            Function[.5*BesselJ[0, #/1.5]]
            },
        "OptimizedBasisSize"->15
        },
    "PlotDisplayMode"->List,
    ColorFunction->"Rainbow",
    "PlotMode"->"Density",
    "WavefunctionClipping"->None
    ]//GraphicsRow[#1,ImageSize->600]&//AbsoluteTiming;
time2
(*Out:*)

0.369797`

plots2
(*Out:*)

dvr-8743390599993568298

And if we check our frequencies they match up (although the second takes much less time to compute):

#-#[[1]]&@
  dvr2D[
    "Energies",
    "Points"->{80, 80},
    "PotentialFunction"->potCoup,
    "NumberOfWavefunctions"->10
    ]
(*Out:*)

{0.`,7.972492994667846`*^-7,0.05925895275277071`,0.05925980369028139`,0.27000204785366577`,0.27001449442769854`,0.2885098059199436`,0.28861302316073534`,0.31971689444992535`,0.31982572504330165`}

#-#[[1]]&@
  dvr2D[
    "Energies",
    "Points"->{125, 125},
    "PotentialFunction"->potCoup,
    "NumberOfWavefunctions"->10,
    "PotentialOptimizationOptions"->
      {
        "PotentialFunction"->
          {
            Function[.5*BesselJ[0, #/1.5]],
            Function[.5*BesselJ[0, #/1.5]]
            },
        "OptimizedBasisSize"->15
        }
    ]
(*Out:*)

{0.`,8.062354943660921`*^-7,0.059847641045838174`,0.059848502458651254`,0.2720771145047962`,0.272089773331718`,0.2905200297948607`,0.29062477560112576`,0.3228685824490185`,0.32298081748454166`}

‌

Miscellaneous Classes

Currently these are the supported DVR classes:

(*Out:*)

{"Cartesian1DDVR","Cartesian2DDVR","Cartesian3DDVR","Cartesian4DDVR","DiskDVR","LegendreDVR","RadialDVR","RingDVR"}

Still need to write up how to use this for

  • Other options tuning

  • Custom DVRs

Clone this wiki locally