I have been developing a QTAIM package for the Wolfram Language. It won a “Staff Pick” prize! Here is a notebook that contains some of the highlights of the package’s capabilities. The package is made in honor of my friend and mentor the late Prof. Richard F. W. Bader. You can follow this project at Github.
Preamble
Preamble
In[]:=
SetDirectory[NotebookDirectory[]]
Out[]=
/Users/ecbrown/src/QTAIM.wl/QTAIM
Parallel Kernels
Parallel Kernels
In[]:=
kernels=LaunchKernels[];
In[]:=
(*CloseKernels[kernels]*)
In[]:=
Needs["QTAIM`"];(*RequirescompilationtoC,andsoevaluateone–at–a–timetopreventraceconditionin$TemporaryDirectory*)Map[ParallelEvaluate[Needs["QTAIM`"],{#}]&,kernels];
In[]:=
(*Simplecommandtotestparallelinitialization:*)(*ParallelEvaluate[$QTAIMContoursPositive]*)
MathLink
MathLink
You may compile a MathLink executable to load Fortran-formatted WFN files.
In[]:=
link=Install[$HomeDirectory<>"/src/QTAIM.wl/QTAIM/QTAIM/qtaim"];
In[]:=
(*Uninstall[link];*)(*ParallelEvaluate[Uninstall[link]];*)
Wavefunction
Wavefunction
Import Existing Wavefunction from Web Server:
In[]:=
w=Import["https://raw.githubusercontent.com/ecbrown/QTAIM.wl/master/resource/water–cas1010–augccpvtz–opt.wdx"];
Alternative ways to import wavefunctions:
From a WFX on a HTTPS Server (Water)
From a WFX on a HTTPS Server (Water)
From a WFX on a HTTPS Server (Pyridine)
From a WFX on a HTTPS Server (Pyridine)
From a WDX on a HTTPS Server (Insulin)
From a WDX on a HTTPS Server (Insulin)
From PySCF
From PySCF
In[]:=
(*session=StartExternalSession["Python"]*)session=StartExternalSession[<|"System""Python","Target""/Users/ecbrown/anaconda3/envs/pyscf/bin/python"|>]
Out[]=
ExternalSessionObject
 |  |
| ||||
In[]:=
importsInput="from pyscf import gtofrom pyscf import libfrom pyscf.x2c import x2cimport numpyfrom pyscf import scf, mcscf, symmfrom pyscf.tools import molden";
In[]:=
ExternalEvaluate[session,importsInput];
/Users/ecbrown/anaconda3/envs/pyscf/lib/python3.12/site–packages/pyscf/dft/libxc.py:771: UserWarning: Since PySCF–2.3, B3LYP (and B3P86) are changed to the VWN–RPA variant, corresponding to the original definition by Stephens et al. (issue 1480) and the same as the B3LYP functional in Gaussian. To restore the VWN5 definition, you can put the setting "B3LYP_WITH_VWN5 = True" in pyscf_conf.py
warnings.warn('Since PySCF–2.3, B3LYP (and B3P86) are changed to the VWN–RPA variant, '
In[]:=
wavefunctionInput="mol = gto.M(atom='O 0.0000000000 0.0000000000 0.1182274717; H 0.0000000000 1.4397666095 –0.9948658513; H 0.0000000000 –1.4397666095 –0.9948658513', unit='B', basis='ccpvdz', verbose=1, symmetry=1, symmetry_subgroup='c2v')mf = scf.RHF(mol).run()coeff = mf.mo_coeff[:,mf.mo_occ>0]energy = mf.mo_energy[mf.mo_occ>0]occ = mf.mo_occ[mf.mo_occ>0]mc = mcscf.CASSCF(mf, 10, 10)mc.kernel()";
In[]:=
ExternalEvaluate[session,wavefunctionInput];
In[]:=
w=ReadWavefunctionFromPySCF[session];
In[]:=
DeleteObject[session]
From a WFN File
From a WFN File
Export Examples
Export Examples
Electron Density and its Derivatives
Electron Density and its Derivatives
The electron density, its gradient, and its Hessian convenience functions are defined so that:
In[]:=
rho[w,{0.,0.,0.}]
Out[]=
10.2816
In[]:=
g[w,{0.,0.,0.}]
Out[]=
{8.90043×,–6.57391×,146.405}
–15
10
–13
10
In[]:=
H[w,{0.,0.,0.}]
Out[]=
{{–654.091,–6.5238×,1.39621×},{–6.5238×,–664.217,–9.69122×},{1.39621×,–9.69122×,2285.48}}
–14
10
–13
10
–14
10
–13
10
–13
10
–13
10
Nuclear Critical Points
Nuclear Critical Points
In[]:=
ncps=LocateNuclearCriticalPoints[w]//Chop;
In[]:=
ncps//TableForm
Out[]//TableForm=
0 | 0 | 0.222366 |
0 | 1.40105 | –0.859873 |
0 | –1.40105 | –0.859873 |
In[]:=
Graphics3D[Table[Sphere[ncps[[i]],0.1*((w["AtomicNumbers"])[[i]])],{i,1,Length[ncps]}]]
Out[]=
Bounding Box
Bounding Box
In[]:=
bbpad=3;bbxmin=Floor[Min[ncps[[All,1]]]]–bbpad;bbxmax=Ceiling[Max[ncps[[All,1]]]]+bbpad;bbymin=Floor[Min[ncps[[All,2]]]]–bbpad;bbymax=Ceiling[Max[ncps[[All,2]]]]+bbpad;bbzmin=Floor[Min[ncps[[All,3]]]]–bbpad;bbzmax=Ceiling[Max[ncps[[All,3]]]]+bbpad;bb={{bbxmin,bbxmax},{bbymin,bbymax},{bbzmin,bbzmax}};
The bounding box:
In[]:=
N@bb
Out[]=
{{–3.,3.},{–5.,5.},{–4.,4.}}
In[]:=
rhoPlot=ContourPlot[rho[w,{0.,y,z}],{y,bbymin,bbymax},{z,bbzmin,bbzmax},FrameLabel{"y","z"},Contours–>$QTAIMContoursPositive,PlotRangeAll,PlotPoints50,ColorFunction"PearlColors",ColorFunctionScalingFalse,AspectRatioAutomatic]
Out[]=
In[]:=
Graphics3D[Table[Sphere[ncps[[i]],0.1*((w["AtomicNumbers"])[[i]])],{i,1,Length[ncps]}],PlotRange–>bb]
Out[]=
Subdivision
Subdivision
Small Molecules
Small Molecules
It is beneficial to make a grid that has the nuclear critical points on the corners. This only works for small molecules!
In[]:=
xrange=Join[{–Infinity},Sort[DeleteDuplicates[Chop[ncps[[All,1]]],(EuclideanDistance[#1,#2]<10^–6)&]],{Infinity}];yrange=Join[{–Infinity},Sort[DeleteDuplicates[Chop[ncps[[All,2]]],(EuclideanDistance[#1,#2]<10^–6)&]],{Infinity}];zrange=Join[{–Infinity},Sort[DeleteDuplicates[Chop[ncps[[All,3]]],(EuclideanDistance[#1,#2]<10^–6)&]],{Infinity}];
Make sure the number of subdivisions is reasonable:
In[]:=
Length[Flatten[Table[1,{xmin,1,Length[xrange]–1},{ymin,1,Length[yrange]–1},{zmin,1,Length[zrange]–1}],2]]
Out[]=
24
In[]:=
q=AbsoluteTiming[ParallelSum[NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,xrange[[xmin]],xrange[[xmin+1]]},{y,yrange[[ymin]],yrange[[ymin+1]]},{z,zrange[[zmin]],zrange[[zmin+1]]},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity],{xmin,1,Length[xrange]–1},{ymin,1,Length[yrange]–1},{zmin,1,Length[zrange]–1},Method–>"FinestGrained"]]
Out[]=
{6.7058,9.84216}
Large Molecules
Large Molecules
For large molecules, it can be beneficial to make a grid that is on the order of the number of processors one has.
In[]:=
xrange=Subdivide[bbxmin,bbxmax,3];yrange=Subdivide[bbymin,bbymax,3];zrange=Subdivide[bbzmin,bbzmax,3];
In[]:=
Length[Flatten[Table[1,{xmin,1,Length[xrange]–1},{ymin,1,Length[yrange]–1},{zmin,1,Length[zrange]–1}],2]]
Out[]=
27
In[]:=
q=AbsoluteTiming[NIntegrate[rho[w,{x,y,z}],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{1.25324,9.95761}
Parallel Example
Parallel Example
Here, an even number of subdivisions is helpful so that the nuclear critical points are on a boundary. Then, NIntegrate will give them special treatment as it does for all boundaries.
In[]:=
xrange=Subdivide[bbxmin,bbxmax,2];yrange=Subdivide[bbymin,bbymax,2];zrange=Subdivide[bbzmin,bbzmax,2];q=AbsoluteTiming[ParallelSum[NIntegrate[rho[w,{x,y,z}],{x,xrange[[ix]],xrange[[ix+1]]},{y,xrange[[iy]],xrange[[iy+1]]},{z,xrange[[iz]],xrange[[iz+1]]},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal2,PrecisionGoalInfinity],{ix,Length[xrange]–1},{iy,Length[yrange]–1},{iz,Length[zrange]–1}]]
Out[]=
{0.141688,9.82336}
Molecular Graph
Molecular Graph
Evaluate the molecular graph, all critical points corresponding to nuclear, bond, ring, and cage critical points. This may be slower than heuristic methods, e.g. start between two NCPs but it is general, for use in cases where one may not even know where to begin. (It solves tetrahedrane, see last example.)
(It is normal to see CompiledFunction and other errors here. These will be quieted in a later version.)
In[]:=
AbsoluteTiming[cps=Quiet[LocateCriticalPoints[{gx[w,{x,y,z}],gy[w,{x,y,z}],gz[w,{x,y,z}]},{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax}]]]
Out[]=
{80.7234,{{–1.82457×,–1.40105,–0.859873},{1.92786×,1.40105,–0.859873},{1.92786×,1.40105,–0.859873},{4.62595×,1.16252,–0.680289},{2.94398×,–1.16252,–0.680289},{2.94398×,–1.16252,–0.680289}}}
–16
10
–17
10
–17
10
–17
10
–16
10
–16
10
Characterize the rank and signature, and ellipticity of the critical points with these helper functions:
In[]:=
CriticalPointRank[{x_,y_,z_}]:=Total[Map[If[#==0,0,1]&,Chop[Eigenvalues[H[w,{x,y,z}]]]]];CriticalPointSignature[{x_,y_,z_}]:=Total[Map[Which[#<0,–1,#==0,0,#>0,1]&,Chop[Eigenvalues[H[w,{x,y,z}]]]]];Ellipticity[{x_,y_,z_}]:=Module[{evals=Sort[Chop[Eigenvalues[H[w,{x,y,z}],–2]]]},(evals[[1]]/evals[[2]])–1];
In[]:=
uniqueCps=DeleteDuplicates[Chop[cps],EuclideanDistance[#1,#2]<1*^–8&]
Out[]=
{{0,–1.40105,–0.859873},{0,1.40105,–0.859873},{0,1.16252,–0.680289},{0,–1.16252,–0.680289}}
In[]:=
cpsTable=Select[Transpose[{uniqueCps[[All,1]],uniqueCps[[All,2]],uniqueCps[[All,3]],Map[CriticalPointRank,uniqueCps],Map[CriticalPointSignature,uniqueCps],Map[Ellipticity,uniqueCps]}],#[[5]]!=–3&]
Out[]=
{{0,1.16252,–0.680289,3,–1,–2.79952},{0,–1.16252,–0.680289,3,–1,–2.79952}}
In[]:=
Graphics3D[Join[Table[{Black,Sphere[cpsTable[[i,1;;3]],0.1]},{i,1,Length[cpsTable]}],Table[Sphere[ncps[[i]],0.1*((w["AtomicNumbers"])[[i]])],{i,1,Length[ncps]}]]]
Out[]=
In[]:=
bcps=Select[Transpose[{uniqueCps[[All,1]],uniqueCps[[All,2]],uniqueCps[[All,3]],Map[CriticalPointRank,uniqueCps],Map[CriticalPointSignature,uniqueCps]}],#[[5]]==–1&];
In[]:=
rcps=Select[Transpose[{uniqueCps[[All,1]],uniqueCps[[All,2]],uniqueCps[[All,3]],Map[CriticalPointRank,uniqueCps],Map[CriticalPointSignature,uniqueCps]}],#[[5]]==1&];
In[]:=
ccps=Select[Transpose[{uniqueCps[[All,1]],uniqueCps[[All,2]],uniqueCps[[All,3]],Map[CriticalPointRank,uniqueCps],Map[CriticalPointSignature,uniqueCps]}],#[[5]]==3&];
In[]:=
bcps//TableForm
Out[]//TableForm=
0 | 1.16252 | –0.680289 | 3 | –1 |
0 | –1.16252 | –0.680289 | 3 | –1 |
In[]:=
rcps//TableForm
Out[]//TableForm=
{}
In[]:=
ccps//TableForm
Out[]//TableForm=
{}
Parallel Algorithm (WIP)
Parallel Algorithm (WIP)
Unlike the case of integration, it is helpful to make an odd number of divisions, so that the planes of symmetry are not on a boundary. (Parallelization is probably not necessary or could be worse for small molecules)
In[]:=
(*xrange=Subdivide[bbxmin,bbxmax,3];yrange=Subdivide[bbymin,bbymax,3];zrange=Subdivide[bbzmin,bbzmax,3];AbsoluteTiming[cps=ParallelTable[LocateCriticalPoints[{gx[w,{x,y,z}],gy[w,{x,y,z}],gz[w,{x,y,z}]},{x,xrange[[ix]],xrange[[ix+1]]},{y,yrange[[iy]],yrange[[iy+1]]},{z,zrange[[iz]],zrange[[iz+1]]}(*,Method–>{"Newton"},Jacobian:>H[w,{x,y,z}]*)],{ix,Length[xrange]–1},{iy,Length[yrange]–1},{iz,Length[zrange]–1}]]*)
Bond Paths
Bond Paths
In[]:=
bondPaths=ParallelTable[BondPath[w,ncps,bcps[[bp,1;;3]]],{bp,1,Length[bcps]}];
In[]:=
bondPathLines3D=ParallelTable[Line[Join[bondPaths[[bp,(*bp*)1,(*dir*)2,1]],Reverse[bondPaths[[bp,(*bp*)2,(*dir*)2,1]]]][[All,2]]],{bp,1,Length[bondPaths]}];
In[]:=
graph=Show[Graphics3D[bondPathLines3D],Graphics3D[Join[Table[{Red,Sphere[bcps[[i,1;;3]],0.1]},{i,1,Length[bcps]}],Table[{Green,Sphere[rcps[[i,1;;3]],0.1]},{i,1,Length[rcps]}],Table[{Blue,Sphere[ccps[[i,1;;3]],0.1]},{i,1,Length[ccps]}],Table[Sphere[ncps[[i]],0.1*((w["AtomicNumbers"])[[i]])],{i,1,Length[ncps]}]]]]
Out[]=
In[]:=
bondPathLines=Table[Line[Join[bondPaths[[bp,(*bp*)1,(*dir*)2,1]],Reverse[bondPaths[[bp,(*bp*)2,(*dir*)2,1]]]][[All,2]][[All,{2,3}]]],{bp,1,Length[bondPaths]}];
In[]:=
rhoPlot=ContourPlot[rho[w,{0.,y,z}],{y,bbymin,bbymax},{z,bbzmin,bbzmax},FrameLabel{"y","z"},Contours–>$QTAIMContoursPositive,PlotRangeAll,PlotPoints50,ColorFunction"PearlColors",ColorFunctionScalingFalse,AspectRatioAutomatic]
Out[]=
With some trial and error, additional contours illustrate BCPs:
In[]:=
rhoPlot=ContourPlot[rho[w,{0.,y,z}],{y,bbymin,bbymax},{z,bbzmin,bbzmax},FrameLabel{"y","z"},Contours–>DeleteDuplicates[Sort[Join[$QTAIMContoursPositive,{13/40,14/40,15/40}]]],PlotRangeAll,PlotPoints50,ColorFunction"PearlColors",ColorFunctionScalingFalse,AspectRatioAutomatic]
Out[]=
In[]:=
Show[rhoPlot,Graphics[bondPathLines]]
Out[]=
Gradient of the Electron Density
Gradient of the Electron Density
In[]:=
gyz[wfn_,{y_?NumericQ,z_?NumericQ}]:=Rest[g[wfn,{0,y,z}]];streamPlot=StreamPlot[gyz[w,{y,z}],{y,bbymin,bbymax},{z,bbzmin,bbzmax},FrameLabel{"y","z"},StreamStyle–>Black,StreamColorFunctionNone,AspectRatioAutomatic]
Out[]=
In[]:=
Show[rhoPlot,streamPlot]
Out[]=
In[]:=
Show[SliceContourPlot3D[rho[w,{x,y,z}],"CenterPlanes",{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},AxesLabel{"x","y","z"},Contours–>DeleteDuplicates[Sort[Join[$QTAIMContoursPositive,{13/40,14/40,15/40}]]],PlotRange–>{Full,Full,Full,All},BoxRatiosAutomatic],graph]
Out[]=
In[]:=
Show[SliceVectorPlot3D[g[w,{x,y,z}],"CenterPlanes",{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},AxesLabel{"x","y","z"},PlotRange–>All,BoxRatiosAutomatic],graph]
Out[]=
(Negative of the) Laplacian of the Electron Density
(Negative of the) Laplacian of the Electron Density
It is fairly easy to define custom fields. For example, the Laplacian, its gradient, and its Hessian can be expressed using the ElectronDensityDerivative convenience function.
In[]:=
laprho[wfn_,{x_?NumericQ,y_?NumericQ,z_?NumericQ}]:=Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{2,0,0},{0,2,0},{0,0,2}}]]];gxl[wfn_,{x_?NumericQ,y_?NumericQ,z_?NumericQ}]:=Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{3,0,0},{1,2,0},{1,0,2}}]]];gyl[wfn_,{x_?NumericQ,y_?NumericQ,z_?NumericQ}]:=Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{2,1,0},{0,3,0},{0,1,2}}]]];gzl[wfn_,{x_?NumericQ,y_?NumericQ,z_?NumericQ}]:=Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{2,0,1},{0,2,1},{0,0,3}}]]];gl[wfn_,{x_?NumericQ,y_?NumericQ,z_?NumericQ}]:={Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{3,0,0},{1,2,0},{1,0,2}}]]],Total[First[ElectronDensityDerivative[w,{{x,y,z}},{{2,1,0},{0,3,0},{0,1,2}}]]],Total[First[ElectronDensityDerivative[w,{{x,y,z}},{{2,0,1},{0,2,1},{0,0,3}}]]]};Hl[wfn_,{x_?NumericQ,y_?NumericQ,z_?NumericQ}]:=Module[{h},h={{Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{4,0,0},{2,2,0},{2,0,2}}]]],Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{3,1,0},{1,3,0},{1,1,2}}]]],Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{3,0,1},{1,2,1},{1,0,3}}]]]},{Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{3,1,0},{1,3,0},{1,1,2}}]]],Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{2,2,0},{0,4,0},{0,2,2}}]]],Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{2,1,1},{0,3,1},{0,1,3}}]]]},{Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{3,0,1},{1,2,1},{1,0,3}}]]],Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{2,1,1},{0,3,1},{0,1,3}}]]],Total[First[ElectronDensityDerivative[wfn,{{x,y,z}},{{2,0,2},{0,2,2},{0,0,4}}]]]}}];
A plot of the Laplacian (unfortunately some “tearing” that needs to be resolved):
In[]:=
SliceContourPlot3D[–laprho[w,{x,y,z}],"CenterPlanes",{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},AxesLabel{"x","y","z"},Contours–>$QTAIMContours,PlotRange–>{Full,Full,Full,All},BoxRatiosAutomatic,ColorFunction"PearlColors",ColorFunctionScalingTrue]
Out[]=
In[]:=
ContourPlot[–laprho[w,{x,0,z}],{z,bbzmin,bbzmax},{x,bbxmin,bbxmax},FrameLabel{"z","x"},Contours–>$QTAIMContours,PlotRangeAll,ColorFunction"PearlColors",ColorFunctionScalingFalse,AspectRatioAutomatic]
Out[]=
Locate the critical points in the 3D Negative of Laplacian:
In[]:=
AbsoluteTiming[cpsl=Quiet[LocateCriticalPoints[{–gxl[w,{x,y,z}],–gyl[w,{x,y,z}],–gzl[w,{x,y,z}]},{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax}]]]
Out[]=
{4269.94,{{–1.00855,–0.576984,–0.09957},{–1.00855,0.576984,–0.09957},{–0.559035,–1.37783×,0.527661},{–0.44144,–1.10201,–1.47436},{–0.44144,1.10201,–1.47436},{–0.414831,–0.491558,–0.0637804},{–0.414831,0.491558,–0.0637804},{–0.218878,1.35982×,–0.935393},{–0.218878,1.35948×,–0.935393},{–0.0144062,–0.0000927957,0.0530215},{–1.15722×,–1.26188,–1.6019},{–1.90532×,1.84799,–1.42824},{–1.78226×,1.84799,–1.42824},{–1.77839×,1.84799,–1.42824},{–1.35247×,0.682076,0.433766},{–1.35247×,0.682076,0.433766},{–6.21382×,–1.44525×,0.0525036},{–5.64053×,0.901504,–0.473607},{–2.54044×,0.0523758,0.380154},{–1.50776×,0.0607234,0.377151},{–1.50776×,0.0607234,0.377151},{–1.50776×,0.0607234,0.377151},{–1.50776×,0.0607234,0.377151},{–1.50776×,0.0607234,0.377151},{–1.49325×,–1.43568,–0.885035},{2.43488×,1.43568,–0.885035},{6.82512×,–0.59876,–0.210374},{3.50881×,0.59876,–0.210374},{4.18046×,–0.901504,–0.473607},{4.18046×,–0.901504,–0.473607},{6.72433×,–0.682076,0.433766},{7.52043×,1.39573×,–0.518136},{7.5231×,1.396×,–0.518136},{1.33619×,1.07675,0.641846},{1.34534×,1.07675,0.641846},{1.34534×,1.07675,0.641846},{1.34534×,1.07675,0.641846},{2.20227×,1.26188,–1.6019},{2.21755×,1.26188,–1.6019},{2.62645×,–1.07675,0.641846},{7.71931×,6.86393×,1.4302},{9.39791×,–1.84799,–1.42824},{9.39791×,–1.84799,–1.42824},{9.39791×,–1.84799,–1.42824},{1.44238×,–1.06119,–1.58325},{0.218878,1.31239×,–0.935393},{0.414831,0.491558,–0.0637804},{0.414831,–0.491558,–0.0637804},{0.559035,–8.79773×,0.527661},{1.00855,0.576984,–0.09957},{1.00855,–0.576984,–0.09957}}}
–14
10
–13
10
–13
10
–13
10
–14
10
–14
10
–14
10
–14
10
–14
10
–15
10
–14
10
–16
10
–16
10
–16
10
–16
10
–16
10
–16
10
–16
10
–16
10
–17
10
–16
10
–15
10
–15
10
–15
10
–15
10
–15
10
–13
10
–15
10
–13
10
–14
10
–14
10
–14
10
–14
10
–14
10
–14
10
–14
10
–14
10
–14
10
–14
10
–14
10
–14
10
–12
10
–13
10
–15
10
In[]:=
CriticalPointRank[{x_,y_,z_}]:=Total[Map[If[#==0,0,1]&,Chop[Eigenvalues[–Hl[w,{x,y,z}]]]]];CriticalPointSignature[{x_,y_,z_}]:=Total[Map[Which[#<0,–1,#==0,0,#>0,1]&,Chop[Eigenvalues[–Hl[w,{x,y,z}]]]]];
In[]:=
uniqueCpsl=DeleteDuplicates[Chop[cpsl],EuclideanDistance[#1,#2]<1*^–8&];
In[]:=
cpslTable=Select[Transpose[{uniqueCpsl[[All,1]],uniqueCpsl[[All,2]],uniqueCpsl[[All,3]],Map[CriticalPointRank,uniqueCpsl],Map[CriticalPointSignature,uniqueCpsl]}],#[[5]]==–3&]
Out[]=
{{–0.559035,0,0.527661,3,–3},{0,–1.43568,–0.885035,3,–3},{0,1.43568,–0.885035,3,–3},{0,–0.59876,–0.210374,3,–3},{0,0.59876,–0.210374,3,–3},{0.559035,0,0.527661,3,–3}}
Lone Pairs of Electrons on the Oxygen:
In[]:=
Graphics3D[Join[Table[{Blue,Sphere[cpslTable[[i,1;;3]],0.1]},{i,1,Length[cpslTable]}],Table[Sphere[ncps[[i]],0.025*((w["AtomicNumbers"])[[i]])],{i,1,Length[ncps]}]]]
Out[]=
Better Bond Path Style Using the (Negative of the) Laplacian of the Electron Density
Better Bond Path Style Using the (Negative of the) Laplacian of the Electron Density
The sign of the Laplacian at the Bond Critical Point can be used to determine whether a bond path should be solid or dashed. The interpolating functions returned by the steepest ascent algorithm are smooth and show banana-like character in good aesthetic detail. Extract the interpolating functions for the forward/backward bond paths:
In[]:=
bondPathInterpolatingFunctions=Table[{bondPaths[[bp,(*bp*)1,(*dir*)1]],Reverse[bondPaths[[bp,(*bp*)2,(*dir*)1]]]},{bp,1,Length[bondPaths]}];
Note the sign of the negative of the Laplacian, we use these to select dashed or solid bond paths:
In[]:=
bcpNegativeLaplacians=Map[(–laprho[w,#]&),bcps[[All,1;;3]]]
Out[]=
{2.64003,2.64003}
Generate the bond paths as ParametricPlot3D functions:
In[]:=
bondPathPlots=Table[(*backward*)s0=(bondPathInterpolatingFunctions[[bcp,All,All]][[1]])[[0]][[1]][[1]][[1]];sf=(bondPathInterpolatingFunctions[[bcp,All,All]][[1]])[[0]][[1]][[1]][[2]];forwardPlot=ParametricPlot3D[(bondPathInterpolatingFunctions[[bcp,All,All]][[1]])/.{QTAIM`Private`s–>s},{s,s0,sf},PlotStyleIf[bcpNegativeLaplacians[[bcp]]>0,Directive[Black,Thick],Directive[Black,Dashed]],PlotRange–>{{–4,4},{–4,4},{–4,4}}];(*forward*)s0=(bondPathInterpolatingFunctions[[bcp,All,All]][[2]])[[0]][[1]][[1]][[1]];sf=(bondPathInterpolatingFunctions[[bcp,All,All]][[2]])[[0]][[1]][[1]][[2]];backwardPlot=ParametricPlot3D[(bondPathInterpolatingFunctions[[bcp,All,All]][[2]])/.{QTAIM`Private`s–>s},{s,s0,sf},PlotStyleIf[bcpNegativeLaplacians[[bcp]]>0,Directive[Black,Thick],Directive[Black,Dashed]],PlotRange–>{{–4,4},{–4,4},{–4,4}}];{forwardPlot,backwardPlot},{bcp,1,Length[bcps]}]
Out[]=
,
,
,
In[]:=
graph=Show[Graphics3D[Join[Table[{Red,Sphere[bcps[[i,1;;3]],0.1]},{i,1,Length[bcps]}],Table[{Green,Sphere[rcps[[i,1;;3]],0.1]},{i,1,Length[rcps]}],Table[{Blue,Sphere[ccps[[i,1;;3]],0.1]},{i,1,Length[ccps]}],Table[Sphere[ncps[[i]],0.1*((w["AtomicNumbers"])[[i]])],{i,1,Length[ncps]}]]],Flatten[bondPathPlots]]
Out[]=
In the plane of the molecule:
In[]:=
bondPathPlotLines=Table[(*backward*)s0=(bondPathInterpolatingFunctions[[bcp,All,All]][[1]])[[0]][[1]][[1]][[1]];sf=(bondPathInterpolatingFunctions[[bcp,All,All]][[1]])[[0]][[1]][[1]][[2]];forwardPlot=ParametricPlot[((bondPathInterpolatingFunctions[[bcp,All,All]][[1]])/.{QTAIM`Private`s–>s})[[{2,3}]],{s,s0,sf},PlotStyleIf[bcpNegativeLaplacians[[bcp]]>0,Directive[Black,Thick],Directive[Black,Dashed]],PlotRange–>{{–4,4},{–4,4}}];(*forward*)s0=(bondPathInterpolatingFunctions[[bcp,All,All]][[2]])[[0]][[1]][[1]][[1]];sf=(bondPathInterpolatingFunctions[[bcp,All,All]][[2]])[[0]][[1]][[1]][[2]];backwardPlot=ParametricPlot[((bondPathInterpolatingFunctions[[bcp,All,All]][[2]])/.{QTAIM`Private`s–>s})[[{2,3}]],{s,s0,sf},PlotStyleIf[bcpNegativeLaplacians[[bcp]]>0,Directive[Black,Thick],Directive[Black,Dashed]],PlotRange–>{{–4,4},{–4,4}}];{forwardPlot,backwardPlot},{bcp,1,Length[bcps]}]
Out[]=
,
,
,
In[]:=
Show[rhoPlot,bondPathPlotLines]
Out[]=
Optimizing the Beta Spheres
Optimizing the Beta Spheres
In[]:=
q=Reap[AbsoluteTiming[NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal1,PrecisionGoalInfinity,EvaluationMonitor:>Sow[{x,y,z}]]]];integrationPoints=q[[2,1]];
Determine the associated nuclear attractor for each point.
In[]:=
anas=ParallelMap[AssociatedNuclearAttractor[w,ncps,#]&,integrationPoints];
In[]:=
beta0=Table[otherPoints=Select[Transpose[{integrationPoints,anas}],(#[[2]]!=ncp)&][[All,1]];(80/100)EuclideanDistance[ncps[[ncp]],First@Nearest[otherPoints,ncps[[ncp]],1]],{ncp,1,Length[ncps]}]
Out[]=
{1.39176,0.292368,0.292368}
Comparing ODE Integration Methods and the Effects of a better Beta Sphere
Comparing ODE Integration Methods and the Effects of a better Beta Sphere
Default (Conservative) Beta Sphere Radius
Default (Conservative) Beta Sphere Radius
In[]:=
AbsoluteTiming[anasAutomatic=Map[AssociatedNuclearAttractor[w,ncps,#,Method–>Automatic]&,integrationPoints];]
Out[]=
{83.7152,Null}
In[]:=
AbsoluteTiming[anasAdams=Map[AssociatedNuclearAttractor[w,ncps,#,Method–>"Adams"]&,integrationPoints];]
Out[]=
{83.1683,Null}
In[]:=
AbsoluteTiming[anasABM=Map[AssociatedNuclearAttractor[w,ncps,#,Method–>"ABM"]&,integrationPoints];]
Out[]=
{147.062,Null}
In[]:=
AbsoluteTiming[anasBDF=Map[AssociatedNuclearAttractor[w,ncps,#,Method–>"BDF"]&,integrationPoints];]
Out[]=
{108.649,Null}
In[]:=
AbsoluteTiming[anasBS54=Map[AssociatedNuclearAttractor[w,ncps,#,Method–>"BS54"]&,integrationPoints];]
Out[]=
{279.456,Null}
In[]:=
AbsoluteTiming[anasCashKarp54=Map[AssociatedNuclearAttractor[w,ncps,#,Method–>"CashKarp45"]&,integrationPoints];]
Out[]=
{201.346,Null}
In[]:=
AbsoluteTiming[anasDOPRI54=Map[AssociatedNuclearAttractor[w,ncps,#,Method–>"DOPRI54"]&,integrationPoints];]
Out[]=
{192.267,Null}
In[]:=
AbsoluteTiming[anasFehlberg45=Map[AssociatedNuclearAttractor[w,ncps,#,Method–>"Fehlberg45"]&,integrationPoints];]
Out[]=
{182.245,Null}
With Beta Sphere
With Beta Sphere
With Beta Sphere
In[]:=
AbsoluteTiming[anasAutomaticb0=Map[AssociatedNuclearAttractor[w,ncps,#,BetaSphereRadii–>beta0,Method–>Automatic]&,integrationPoints];]
Out[]=
{24.0827,Null}
In[]:=
AbsoluteTiming[anasAdamsb0=Map[AssociatedNuclearAttractor[w,ncps,#,BetaSphereRadii–>beta0,Method–>"Adams"]&,integrationPoints];]
Out[]=
{24.2595,Null}
In[]:=
AbsoluteTiming[anasABMb0=Map[AssociatedNuclearAttractor[w,ncps,#,BetaSphereRadii–>beta0,Method–>"ABM"]&,integrationPoints];]
Out[]=
{47.9963,Null}
In[]:=
AbsoluteTiming[anasBDFb0=Map[AssociatedNuclearAttractor[w,ncps,#,BetaSphereRadii–>beta0,Method–>"BDF"]&,integrationPoints];]
Out[]=
{34.9828,Null}
In[]:=
AbsoluteTiming[anasBS54b0=Map[AssociatedNuclearAttractor[w,ncps,#,BetaSphereRadii–>beta0,Method–>"BS54"]&,integrationPoints];]
Out[]=
{216.11,Null}
In[]:=
AbsoluteTiming[anasCashKarp54b0=Map[AssociatedNuclearAttractor[w,ncps,#,BetaSphereRadii–>beta0,Method–>"CashKarp45"]&,integrationPoints];]
Out[]=
{150.191,Null}
In[]:=
AbsoluteTiming[anasDOPRI54b0=Map[AssociatedNuclearAttractor[w,ncps,#,BetaSphereRadii–>beta0,Method–>"DOPRI54"]&,integrationPoints];]
Out[]=
{153.693,Null}
In[]:=
AbsoluteTiming[anasFehlberg45b0=Map[AssociatedNuclearAttractor[w,ncps,#,BetaSphereRadii–>beta0,Method–>"Fehlberg45"]&,integrationPoints];]
Out[]=
{135.414,Null}
Plots of Atomic Basins
Plots of Atomic Basins
This package defines a function called AssociatedNuclearAttractor that returns the index of the nuclear critical point that is the terminus of the steepest ascent path from a test point. This is arguably the most “expensive” part of QTAIM: every region of space needs to be tested to ascertain to which NCP it belongs. That could be hundreds of gradient evaluations for even one nominal function evaluation!
AssociatedNuclearAttractor
AssociatedNuclearAttractor
In[]:=
AssociatedNuclearAttractor[w,ncps,{0.,0.,0.}]
Out[]=
1
In[]:=
AssociatedNuclearAttractor[w,ncps,{0.,0.,0.},BetaSphereRadii–>beta0]
Out[]=
1
In[]:=
AssociatedNuclearAttractor[w,ncps,{0.,2.,–2.},BetaSphereRadii–>beta0]
Out[]=
2
In[]:=
AssociatedNuclearAttractor[w,ncps,{0.,–2.,–2.},BetaSphereRadii–>beta0]
Out[]=
3
In[]:=
AssociatedNuclearAttractor[w,ncps,{2.,2.,2.}]
Out[]=
0
Plots
Plots
In[]:=
rp=RegionPlot[{AssociatedNuclearAttractor[w,ncps,{0.,y,z},BetaSphereRadii–>beta0]==1,AssociatedNuclearAttractor[w,ncps,{0.,y,z},BetaSphereRadii–>beta0]==2,AssociatedNuclearAttractor[w,ncps,{0.,y,z},BetaSphereRadii–>beta0]==3},{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyle{Directive[Red,Opacity[0.10]],Directive[Green,Opacity[0.10]],Directive[Blue,Opacity[0.10]]},AspectRatioAutomatic]
Out[]=
In[]:=
Show[rhoPlot,Graphics[bondPathLines],rp]
Out[]=
In[]:=
Show[rhoPlot,Graphics[bondPathLines],rp,streamPlot]
Out[]=
Zoom in on the detail of the one of the bond critical points:
In[]:=
streamPlotZoom=StreamPlot[gyz[w,{y,z}],{y,1.10,1.20},{z,–0.85,–0.75},FrameLabel{"y","z"},StreamStyle–>Black,StreamColorFunctionNone,AspectRatioAutomatic]
Out[]=
Annotate with bond critical point:
In[]:=
Show[streamPlotZoom,Graphics[{Red,Disk[bcps[[1,2;;3]],0.0025]}]]
Out[]=
One of the purposes of this package is to facilitate graphics and plots. Here, it is shown how to draw individual components so that only the necessary regions are plotted. (Otherwise, graphics get complicated and occluded.)
In[]:=
rpO=RegionPlot[AssociatedNuclearAttractor[w,ncps,{0.,y,z},BetaSphereRadii–>beta0]==1,{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Red,Opacity[0.10]],AspectRatioAutomatic]
Out[]=
In[]:=
rpH1=RegionPlot[AssociatedNuclearAttractor[w,ncps,{0.,y,z},BetaSphereRadii–>beta0]==2,{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Green,Opacity[0.10]],AspectRatioAutomatic]
Out[]=
In[]:=
rpH2=RegionPlot[AssociatedNuclearAttractor[w,ncps,{0.,y,z},BetaSphereRadii–>beta0]==3,{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Blue,Opacity[0.10]],AspectRatioAutomatic]
Out[]=
3D Basin plots can be generated in the same way. Note that increasing PlotPoints makes them smoother, with an increase in computer time.
In[]:=
rp3D=RegionPlot3D[{AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==1,AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==2,AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==3},{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyle{Directive[Red,Opacity[0.10]],Directive[Green,Opacity[0.10]],Directive[Blue,Opacity[0.10]]},BoxRatiosAutomatic]
Out[]=
In[]:=
rp3DO=RegionPlot3D[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==1,{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Red,Opacity[0.10]],BoxRatiosAutomatic]
Out[]=
In[]:=
rp3DH1=RegionPlot3D[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==2,{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Green,Opacity[0.10]],BoxRatiosAutomatic]
Out[]=
In[]:=
rp3DH2=RegionPlot3D[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==3,{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Blue,Opacity[0.10]],BoxRatiosAutomatic]
Out[]=
In[]:=
rp3DOnomesh=RegionPlot3D[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==1,{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Red,Opacity[0.10]],Mesh–>None,BoxRatiosAutomatic]
Out[]=
In[]:=
rp3DH1nomesh=RegionPlot3D[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==2,{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Green,Opacity[0.10]],Mesh–>None,BoxRatiosAutomatic]
Out[]=
In[]:=
rp3DH2nomesh=RegionPlot3D[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==3,{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},PlotStyleDirective[Blue,Opacity[0.10]],Mesh–>None,BoxRatiosAutomatic]
Out[]=
Atomic Properties with Integration over Basins in Cartesian Space
Atomic Properties with Integration over Basins in Cartesian Space
(Note that AccuracyGoals are for demo purposes only. Some analyses require increased accuracy/precision.)
It turns out beneficial to develop an integration strategy prior to integrating the irregular atomic basins. Here, we compare what Automatic gives, GlobalAdaptive, vs. LocalAdaptive. Please note these integrations are in Cartesian space and not in e.g. traditional Spherical Polar coordinates as rays shooting out, starting on a nuclear attractor. The hope is to use adaptive integration where needed in x,y, and z, and avoid strong orientation effects in spherical polar coordinates. (These pervade many plots in QTAIM.) There is a different type of mis-specification in Lebedev integration: seldom does one want perfect spherical integration either! Cartesian space representation does not suffer from unreachability problems, such as curved surfaces that rays can not reach, though a steepest ascent path exists. Let’s adapt in Cartesian space if possible, and reserve comparisons for a later date, especially adaptive spherical polar integration.
Long story short: Choose LocalAdaptive.
The whole system:
In[]:=
q=AbsoluteTiming[NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{Automatic,"SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{15.532,9.84063}
In[]:=
q=AbsoluteTiming[NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{"GlobalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{15.651,9.84063}
LocalAdaptive is much faster than GlobalAdaptive, Automatic gets it wrong every time. I wonder if this is because it can re-use expensive function evaluations! This is contrary to what I expected, since GlobalAdaptive is usually superior:
In[]:=
q=AbsoluteTiming[NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{5.95284,9.84116}
LocalAdaptive is ten times faster!
Parallel Processing:
In[]:=
xrange=Subdivide[bbxmin,bbxmax,2];yrange=Subdivide[bbymin,bbymax,2];zrange=Subdivide[bbzmin,bbzmax,2];q=AbsoluteTiming[ParallelSum[NIntegrate[rho[w,{x,y,z}],{x,xrange[[ix]],xrange[[ix+1]]},{y,xrange[[iy]],xrange[[iy+1]]},{z,xrange[[iz]],xrange[[iz+1]]},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity],{ix,Length[xrange]–1},{iy,Length[yrange]–1},{iz,Length[zrange]–1}]]
Out[]=
{0.530324,9.89058}
(Slower, but likely due to overhead for a pretty small molecule subdivided into many pieces. Should also increase accuracy since error is additive over these subdivided regions!)
Let's try a different way, employing two mirror-plane symmetry:
In[]:=
q=AbsoluteTiming[4*NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,0,bbxmax},{y,0,bbymax},{z,bbzmin,bbzmax},Method–>{Automatic,"SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{9.7147,9.83928}
In[]:=
q=AbsoluteTiming[4*NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,0,bbxmax},{y,0,bbymax},{z,bbzmin,bbzmax},Method–>{"GlobalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{9.71118,9.83928}
In[]:=
q=AbsoluteTiming[4*NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,0,bbxmax},{y,0,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{1.69383,9.83864}
Shocking time difference, same integrated charge.
Alternative definition of the envelope, when Associated Nuclear Attractor is not 0. This is probably more like what an atomic basin integration would be like, at least for exterior atoms–terminating by default at the early condition where if ρ <(x0, y0, z0) < Min[$QTAIMContourPositive] then assigned ncp index 0:
Let’s try with LocalAdaptive (skipping GlobalAdaptive this time) for the whole molecule:
In[]:=
qalt=AbsoluteTiming[NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]!=0,rho[w,{x,y,z}],0],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{2148.84,9.67444}
It’s a good thing that we spent some time investigating integration technique before tackling harder problems.
Integrate over individual atomic basins, using LocalAdaptive. It could mean orders of magnitude time difference for an expensive integration:
(Note we can help the integration to locate regions of non-zero density by placing nuclear critical point coordinates in the integration limits.)
In[]:=
qO=AbsoluteTiming[With[{ncp=1},4*NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,rho[w,{x,y,z}],0],{x,0,bbxmax},{y,0,bbymax},{z,bbzmin,ncps[[ncp,3]],bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]]
Out[]=
{763.13,8.93723}
In[]:=
qH1andH2=AbsoluteTiming[With[{ncp=2},4*NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,rho[w,{x,y,z}],0],{x,0,bbxmax},{y,0,ncps[[ncp,2]],bbymax},{z,bbzmin,ncps[[ncp,3]],bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]]
Out[]=
{287.268,0.478105}
Total atomic charge as sum of components:
In[]:=
(qO+qH1andH2)
Out[]=
{1050.4,9.41534}
Difference between whole molecule and fragments:
In[]:=
Last[qalt]–(Last[qO]+Last[qH1andH2])
Out[]=
0.259108
Integrate the Laplacian (should be zeros for each integration, by zero-flux condition). We will use the LocalAdaptive strategy, but it is really a different type of integral and we should come back and compare with GlobalAdaptive/Automatic. It is known that there are problems integrating to values that are nearly 0, and so some warnings may appear. Traditionally this is used as a check for the quality of the basin delineation, but here all grids are different because all integrands are different and adaptively sampled with regard to basin and integral value!
In[]:=
L=AbsoluteTiming[NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],laprho[w,{x,y,z}],0],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{8255.77,–0.640955}
The integrated Laplacian density should be zero if truly a quantum (sub)system bound by the zero-flux condition:
In[]:=
LO=AbsoluteTiming[With[{ncp=1},4*NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,laprho[w,{x,y,z}],0],{x,0,bbxmax},{y,0,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]]
Out[]=
{5566.73,–0.945295}
In[]:=
LH1andH2=AbsoluteTiming[With[{ncp=2},4*NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,laprho[w,{x,y,z}],0],{x,0,bbxmax},{y,0,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]]
Out[]=
{2743.49,–0.347099}
In[]:=
(LO+LH1andH2)
Out[]=
{8310.22,–1.29239}
In[]:=
Last[L]–(Last[LO]+Last[LH1andH2])
Out[]=
0.651438
Volume
Volume
In[]:=
v=AbsoluteTiming[NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],1,0],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]
Out[]=
{185.003,148.308}
In[]:=
vO=AbsoluteTiming[With[{ncp=1},4*NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,1,0],{x,0,bbxmax},{y,0,bbymax},{z,bbzmin,ncps[[ncp,3]],bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal2,PrecisionGoalInfinity]]]
Out[]=
{536.479,78.1077}
In[]:=
vH1andH2=AbsoluteTiming[With[{ncp=2},4*NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,1,0],{x,0,bbxmax},{y,0,ncps[[ncp,2]],bbymax},{z,bbzmin,ncps[[ncp,3]],bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal2,PrecisionGoalInfinity]]]
Out[]=
{263.254,15.0633}
In[]:=
vO+2*vH1andH2
Out[]=
{1062.99,108.234}
Spherical Polar Basin Delineation and Integration
Spherical Polar Basin Delineation and Integration
Interatomic Surface Radius Plot
Interatomic Surface Radius Plot
InteratomicSurface is a function that returns the {position, radius, and associated nuclear critical points} at a parametric ray shooting from a nuclear critical point along θ and ϕ.
In[]:=
(InteratomicSurface[w,ncps,1,{–3π/4,–π/8},BetaSphereRadii–>beta0])[[2]]
Out[]=
2.83735
In[]:=
ias1ContourPlot=ContourPlot[(InteratomicSurface[w,ncps,1,{theta,phi},BetaSphereRadii–>beta0])[[2]],{phi,–π,π},{theta,0,π},PlotRange–>All]
Out[]=
In[]:=
ias2ContourPlot=ContourPlot[(InteratomicSurface[w,ncps,2,{theta,phi},BetaSphereRadii–>beta0])[[2]],{phi,–π,π},{theta,0,π},PlotRange–>All]
Out[]=
In[]:=
ias3ContourPlot=ContourPlot[(InteratomicSurface[w,ncps,3,{theta,phi},BetaSphereRadii–>beta0])[[2]],{phi,–π,π},{theta,0,π},PlotRange–>All]
Out[]=
Hermite Polynomial Representation of rmax as a function of {θ, ϕ}
Hermite Polynomial Representation of rmax as a function of {θ, ϕ}
In[]:=
Options[FunctionInterpolation]
Out[]=
{AccuracyGoalAutomatic,InterpolationOrder3,InterpolationPoints11,InterpolationPrecisionAutomatic,MaxRecursion6,PrecisionGoalAutomatic}
In[]:=
iasFunctionInterpolation=ParallelTable[FunctionInterpolation[(InteratomicSurface[w,ncps,i,{theta,phi},BetaSphereRadii–>beta0])[[2]],{theta,N[0],N[π]},{phi,N[–π],N[π]}],{i,1,Length[ncps],1},Method–>"FinestGrained"];
In[]:=
ContourPlot[(iasFunctionInterpolation[[1]])[theta,phi],{phi,–π,π},{theta,0,π},PlotRange–>All]
Out[]=
In[]:=
ContourPlot[(iasFunctionInterpolation[[2]])[theta,phi],{phi,–π,π},{theta,0,π},PlotRange–>All]
Out[]=
In[]:=
ContourPlot[(iasFunctionInterpolation[[3]])[theta,phi],{phi,–π,π},{theta,0,π},PlotRange–>All]
Out[]=
Partitioning Radial and Angular Components
Partitioning Radial and Angular Components
Split the integration into radial and angular components. Define a function that determines the radial limit and then integrates an arbitrary function f:
In[]:=
Yf[w_,ncps_,ncp_,{theta_?NumericQ,phi_?NumericQ},f_]:=Module[{rmax},rmax=(InteratomicSurface[w,ncps,ncp,{theta,phi},BetaSphereRadii–>beta0,MaxIterations–>10])[[2]];NIntegrate[r^2f[w,{rCos[phi]Sin[theta],rSin[phi]Sin[theta],rCos[theta]}+ncps[[ncp]]],{r,N[0],N[rmax]},Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]];Yfi[w_,ncps_,ncp_,{theta_?NumericQ,phi_?NumericQ},f_]:=Module[{rmax},rmax=(iasFunctionInterpolation[[ncp]])[theta,phi];NIntegrate[r^2f[w,{rCos[phi]Sin[theta],rSin[phi]Sin[theta],rCos[theta]}+ncps[[ncp]]],{r,N[0],N[rmax]},Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]];
Atomic Surface
Atomic Surface
In[]:=
AbsoluteTiming[as1=Reap[With[{ncp=1},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},rho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>2,PrecisionGoalInfinity,Method–>{Automatic,"SymbolicProcessing"–>0},EvaluationMonitor:>Sow[(InteratomicSurface[w,ncps,ncp,{theta,phi},BetaSphereRadii–>beta0])[[1]]]]]];]
Out[]=
{705.985,Null}
In[]:=
With[{ncp=1},Show[ListPointPlot3D[as1[[2,1]],ColorFunctionFunction[{x,y,z},Black]],Graphics3D[Table[Line[{ncps[[ncp]],as1[[2,1,i]]}],{i,1,Length[as1[[2,1]]]}],Gray],graph,BoxRatios–>Automatic,PlotRange–>All]]
Out[]=
Electron Density
Electron Density
In[]:=
(*Experimentalf[theta_?NumericQ]:=Sin[theta]*NIntegrate[Yf[w,ncps,1,{theta,phi},rho],{phi,N[0],N[π]},Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]ParallelNInt[f,{N[–π],N[π]},2]*)
In[]:=
qO=AbsoluteTiming[With[{ncp=1},NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,rho[w,{x,y,z}],0],{x,bbxmin,ncps[[ncp,1]],bbxmax},{y,bbymin,ncps[[ncp,2]],bbymax},{z,bbzmin,ncps[[ncp,3]],bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal2,PrecisionGoalInfinity]]]
Out[]=
{106.75,8.93347}
In[]:=
AbsoluteTiming[intq1=With[{ncp=1},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},rho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>2,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
Out[]=
{89.7349,8.95469}
In[]:=
qO=AbsoluteTiming[With[{ncp=1},NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,rho[w,{x,y,z}],0],{x,bbxmin,ncps[[ncp,1]],bbxmax},{y,bbymin,ncps[[ncp,2]],bbymax},{z,bbzmin,ncps[[ncp,3]],bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal3,PrecisionGoalInfinity]]]
Out[]=
{433.645,8.93744}
In[]:=
AbsoluteTiming[intq1=With[{ncp=1},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},rho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>3,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
Out[]=
{689.079,8.92322}
In[]:=
qO=AbsoluteTiming[With[{ncp=1},NIntegrate[If[AssociatedNuclearAttractor[w,ncps,{x,y,z},BetaSphereRadii–>beta0]==ncp,rho[w,{x,y,z}],0],{x,bbxmin,ncps[[ncp,1]],bbxmax},{y,bbymin,ncps[[ncp,2]],bbymax},{z,bbzmin,ncps[[ncp,3]],bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal4,PrecisionGoalInfinity]]]
Out[]=
{3023.36,8.93723}
In[]:=
AbsoluteTiming[intq1=With[{ncp=1},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},rho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
Out[]=
{4618.82,8.92217}
In[]:=
AbsoluteTiming[intq2=With[{ncp=2},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},rho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
Out[]=
{2055.79,0.36608}
In[]:=
AbsoluteTiming[intq3=With[{ncp=3},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},rho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
Out[]=
{2051.19,0.36608}
In[]:=
intq1+intq2+intq3
Out[]=
9.65433
Laplacian of the Electron Density
Laplacian of the Electron Density
In[]:=
AbsoluteTiming[intl1=With[{ncp=1},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},laprho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
AbsoluteTiming[intl2=With[{ncp=2},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},laprho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
AbsoluteTiming[intl3=With[{ncp=3},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},laprho],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
intl1+intl2+intl3
Kinetic Energy Density G
Kinetic Energy Density G
AbsoluteTiming[intg1=With[{ncp=1},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},KineticEnergyDensityG],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
AbsoluteTiming[intg2=With[{ncp=2},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},KineticEnergyDensityG],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
AbsoluteTiming[intg3=With[{ncp=3},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},KineticEnergyDensityG],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
intg1+intg2+intg3
Kinetic Energy Density K
Kinetic Energy Density K
AbsoluteTiming[intk1=With[{ncp=1},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},KineticEnergyDensityK],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
AbsoluteTiming[intk2=With[{ncp=2},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},KineticEnergyDensityK],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
AbsoluteTiming[intk3=With[{ncp=3},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},KineticEnergyDensityK],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>4,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
intk1+intk2+intk3
Volume
Volume
In[]:=
One[x_,y_]:=1;AbsoluteTiming[vint1=With[{ncp=1},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},One],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>2,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
Out[]=
{3919.5,73.6001}
In[]:=
One[x_,y_]:=1;AbsoluteTiming[vint2=With[{ncp=2},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},One],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>2,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
Out[]=
{1604.25,12.1237}
In[]:=
One[x_,y_]:=1;AbsoluteTiming[vint3=With[{ncp=3},NIntegrate[Sin[theta]*Yf[w,ncps,ncp,{theta,phi},One],{phi,N[–π],N[π]},{theta,N[0],N[π]},AccuracyGoal–>2,PrecisionGoalInfinity,Method–>{"LocalAdaptive","SymbolicProcessing"–>0}]]]
Out[]=
{1586.66,12.1237}
In[]:=
vint1+vint2+vint3
Out[]=
97.8475
Cerjan-Miller-Baker Approach to Locating Critical Points
Cerjan-Miller-Baker Approach to Locating Critical Points
The CMB method (introduced to QTAIM by Popelier) is amenable to direct location of stationary points of a particular signature in any field so long as we can provide its gradient and Hessian. In QTAIM we need separate routines for each corresponding to nuclear (-3), bond (-1), ring (+1), and cage (+3) critical points.
In[]:=
rhoPlot
Out[]=
In[]:=
CriticalPointGradientMinus3[{g[w,{0,0,1}],H[w,{0,0,1}]}]
Out[]=
{0.,7.28116×,–2.38129}
–19
10
In[]:=
CriticalPointGradientMinus3[{g[w,{0,0,–1}],H[w,{0,0,–1}]}]
Out[]=
{0.,–7.80448×,2.3221}
–18
10
In[]:=
CriticalPointRank[{x_,y_,z_}]:=Total[Map[If[#==0,0,1]&,Chop[Eigenvalues[H[w,{x,y,z}]]]]];CriticalPointSignature[{x_,y_,z_}]:=Total[Map[Which[#<0,–1,#==0,0,#>0,1]&,Chop[Eigenvalues[H[w,{x,y,z}]]]]];
In[]:=
AbsoluteTiming[soln=NDSolveValue[{X'[s]==CriticalPointGradientMinus3[{g[w,X[s]],H[w,X[s]]}],X[0]=={ncps[[1,1]]+0.005,ncps[[1,2]]+0.005,ncps[[1,3]]–0.005}},{X[s]},{s,0,10000}(*,Method–>"BDF"*)(*,AccuracyGoal–>2,PrecisionGoalInfinity]*)]]
1
0.
1
0.
1
0.
–7
10
Out[]=
0.848964,InterpolatingFunction
[s]
|  |
|
In[]:=
ParametricPlot3D[soln[[1]],{s,0,soln[[1,0,1,1,2]]},PlotRangeAll]
Out[]=
In[]:=
soln=NDSolveValue[{X'[s]==CriticalPointGradientMinus1[{g[w,X[s]],H[w,X[s]]}],X[0]=={bcps[[1,1]]+0.005,bcps[[1,2]]+0.005,bcps[[1,3]]–0.005}},{X[s]},{s,0,10000}(*,Method–>"BDF",AccuracyGoal–>2,PrecisionGoalInfinity*)]
In[]:=
soln[[1]]
In[]:=
ParametricPlot3D[soln[[1]],{s,0,soln[[1,0,1,1,2]]},PlotRangeAll]
In[]:=
sf=soln[[1,0,1,1,2]]
In[]:=
First[soln/.{s–>soln[[1,0,1,1,2]]}]
In[]:=
FindRoot[CriticalPointGradientMinus1[{g[w,{x,y,z}],H[w,{x,y,z}]}],{x,bcps[[1,1]]+0.005},{y,bcps[[1,2]]+0.005},{z,bcps[[1,3]]–0.005},EvaluationMonitor:>Print[{x,y,z}]]
In[]:=
bcps//TableForm
Out[]//TableForm=
bcps
In[]:=
(*FindRoot[CriticalPointGradientPlus1[{g[w,{x,y,z}],H[w,{x,y,z}]}],{x,rcps[[1,1]]+0.005},{y,rcps[[1,2]]+0.005},{z,rcps[[1,3]]–0.005},EvaluationMonitor:>Print[{x,y,z}]]*)
In[]:=
(*FindRoot[CriticalPointGradientPlus3[{g[w,{x,y,z}],H[w,{x,y,z}]}],{x,ccps[[1,1]]+0.005},{y,ccps[[1,2]]+0.005},{z,ccps[[1,3]]–0.005},EvaluationMonitor:>Print[{x,y,z}]]*)
Parallel Search
Parallel Search
A strategy is to make a grid and search for each type of critical point. Let’s start with the integration of the electron density over the system:
In[]:=
q=Reap[AbsoluteTiming[NIntegrate[If[rho[w,{x,y,z}]>Min[$QTAIMContoursPositive],rho[w,{x,y,z}],0],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},Method–>{"LocalAdaptive","SymbolicProcessing"–>0},AccuracyGoal0,PrecisionGoalInfinity,EvaluationMonitor:>Sow[{x,y,z}]]]];integrationPoints=q[[2,1]];
In[]:=
Show[ListPointPlot3D[integrationPoints,AspectRatio–>Automatic],graph,BoxRatios–>Automatic,PlotRange–>All]
Out[]=
But this is probably way too many points. We would reserve this strategy for careful scans of systems of which we knew very little, and are looking for subtle/long-range bond paths, etc.
Example: Offset slightly from the first bond critical point, and then perform CJM gradient ascent to a -1 critical point:
In[]:=
soln=NDSolveValue[{X'[s]==CriticalPointGradientMinus1[{g[w,X[s]],H[w,X[s]]}],X[0]=={bcps[[1,1]]+0.005,bcps[[1,2]]+0.005,bcps[[1,3]]–0.005}},{X[s]},{s,0,10000},Method–>"BDF",MaxSteps–>100(*,AccuracyGoal–>2,PrecisionGoalInfinity*)]
(The endpoint is buried in the interpolating function object result)
This search starts at each nuclear position and performs a CJM gradient ascent:
In[]:=
AbsoluteTiming[cpsMinus3=Select[DeleteDuplicates[Select[Select[ParallelTable[Quiet[Check[soln=NDSolveValue[{X'[s]==CriticalPointGradientMinus3[{g[w,X[s]],H[w,X[s]]}],X[0]==(w["NuclearCartesianCoordinates"])[[i]]},{X[s]},{s,0,10000},Method–>"BDF",MaxSteps–>100(*,AccuracyGoal–>4,PrecisionGoalInfinity*)];First[soln/.{s–>soln[[1,0,1,1,2]]}],{}]],{i,1,w["NumberOfNuclei"],1},Method"FinestGrained"],Length[#]==3&],((bbxmin<=#[[1]]<=bbxmax)&&(bbymin<=#[[2]]<=bbymax)&&(bbzmin<=#[[3]]<=bbzmax))&],(EuclideanDistance[#1,#2]<1*^–3)&],CriticalPointCharacter[#]==3&&CriticalPointSignature[#]==–3&]];
In[]:=
cpsMinus3//Chop//TableForm
Out[]//TableForm=
{}
A table of starting values in the yz plane:
In[]:=
initialValues=Flatten[Outer[List,{0.},Subdivide[bbymin,bbymax,10],Subdivide[bbzmin,bbzmax,10]],2];
In[]:=
Show[ListPointPlot3D[initialValues,AspectRatio–>Automatic],graph,BoxRatios–>Automatic,PlotRange–>All]
Out[]=
In[]:=
AbsoluteTiming[cpsMinus3=Select[DeleteDuplicates[Select[Select[ParallelTable[Quiet[Check[soln=NDSolveValue[{X'[s]==CriticalPointGradientMinus3[{g[w,X[s]],H[w,X[s]]}],X[0]==initialValues[[i]]},{X[s]},{s,0,10000},Method–>"BDF",MaxSteps–>100(*,AccuracyGoal–>4,PrecisionGoalInfinity*)];First[soln/.{s–>soln[[1,0,1,1,2]]}],{}]],{i,1,Length[initialValues],1},Method"FinestGrained"],Length[#]==3&],((bbxmin<=#[[1]]<=bbxmax)&&(bbymin<=#[[2]]<=bbymax)&&(bbzmin<=#[[3]]<=bbzmax))&],(EuclideanDistance[#1,#2]<1*^–3)&],CriticalPointCharacter[#]==3&&CriticalPointSignature[#]==–3&];]
Out[]=
{2.34588,Null}
In[]:=
cpsMinus3//Chop//TableForm
Out[]//TableForm=
{}
Initial values on the same grid, provided they are enclosed by min/max of NCP coordinates (heuristic):
In[]:=
initialValues=Select[Flatten[Outer[List,{0.},Subdivide[bbymin,bbymax,10],Subdivide[bbzmin,bbzmax,10]],2],(Min[ncps[[All,1]]]<=#[[1]]<=Max[ncps[[All,1]]]&&Min[ncps[[All,2]]]<=#[[2]]<=Max[ncps[[All,2]]]&&Min[ncps[[All,3]]]<=#[[3]]<=Max[ncps[[All,3]]])&];
In[]:=
Show[ListPointPlot3D[initialValues,AspectRatio–>Automatic],graph,BoxRatios–>Automatic,PlotRange–>All]
Out[]=
In[]:=
AbsoluteTiming[cpsMinus1=Select[DeleteDuplicates[Select[Select[ParallelTable[Quiet[Check[soln=NDSolveValue[{X'[s]==CriticalPointGradientMinus1[{g[w,X[s]],H[w,X[s]]}],X[0]==initialValues[[i]]},{X[s]},{s,0,10000},Method–>"BDF",MaxSteps–>100(*,AccuracyGoal–>4,PrecisionGoalInfinity*)];First[soln/.{s–>soln[[1,0,1,1,2]]}],{}]],{i,1,Length[initialValues],1},Method"FinestGrained"],Length[#]==3&],((bbxmin<=#[[1]]<=bbxmax)&&(bbymin<=#[[2]]<=bbymax)&&(bbzmin<=#[[3]]<=bbzmax))&],(EuclideanDistance[#1,#2]<1*^–3)&],CriticalPointCharacter[#]==3&&CriticalPointSignature[#]==–1&];]
Out[]=
{0.079089,Null}
In[]:=
cpsMinus1//Chop//TableForm
Out[]//TableForm=
{}
In[]:=
AbsoluteTiming[cpsPlus1=Select[DeleteDuplicates[Select[Select[ParallelTable[Quiet[Check[soln=NDSolveValue[{X'[s]==CriticalPointGradientPlus1[{g[w,X[s]],H[w,X[s]]}],X[0]==initialValues[[i]]},{X[s]},{s,0,10000},Method–>"BDF",MaxSteps–>100(*,AccuracyGoal–>4,PrecisionGoalInfinity*)];First[soln/.{s–>soln[[1,0,1,1,2]]}],{}]],{i,1,Length[initialValues],1},Method"FinestGrained"],Length[#]==3&],((bbxmin<=#[[1]]<=bbxmax)&&(bbymin<=#[[2]]<=bbymax)&&(bbzmin<=#[[3]]<=bbzmax))&],(EuclideanDistance[#1,#2]<1*^–3)&],CriticalPointCharacter[#]==3&&CriticalPointSignature[#]==+1&];]
Out[]=
{0.216252,Null}
In[]:=
cpsPlus1//Chop//TableForm
Out[]//TableForm=
{}
In[]:=
AbsoluteTiming[cpsPlus3=Select[DeleteDuplicates[Select[Select[ParallelTable[Quiet[Check[soln=NDSolveValue[{X'[s]==CriticalPointGradientPlus3[{g[w,X[s]],H[w,X[s]]}],X[0]==initialValues[[i]]},{X[s]},{s,0,10000},Method–>"BDF",MaxSteps–>100(*,AccuracyGoal–>4,PrecisionGoalInfinity*)];First[soln/.{s–>soln[[1,0,1,1,2]]}],{}]],{i,1,Length[initialValues],1},Method"FinestGrained"],Length[#]==3&],((bbxmin<=#[[1]]<=bbxmax)&&(bbymin<=#[[2]]<=bbymax)&&(bbzmin<=#[[3]]<=bbzmax))&],(EuclideanDistance[#1,#2]<1*^–3)&],CriticalPointCharacter[#]==3&&CriticalPointSignature[#]==+3&];]
Out[]=
{0.216836,Null}
In[]:=
cpsPlus3//Chop//TableForm
Out[]//TableForm=
{}
In[]:=
GraphicsRow[{ncps//Chop//Sort//TableForm,cpsMinus3//Chop//Sort//TableForm}]
Out[]=
In[]:=
GraphicsRow[{bcps//Chop//Sort//TableForm,cpsMinus1//Chop//Sort//TableForm}]
Out[]=
In[]:=
GraphicsRow[{rcps//Chop//Sort//TableForm,cpsPlus1//Chop//Sort//TableForm}]
Out[]=
In[]:=
GraphicsRow[{ccps//Chop//Sort//TableForm,cpsPlus3//Chop//Sort//TableForm}]
Miscellaneous
Miscellaneous
Classic examples from traditional QTAIM and Bader’s AIM book.
Enclosing envelope at “0.001 a.u.”
Mesh:
envMesh=ContourPlot3D[rho[w,{x,y,z}],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},AxesLabel{"x","y","z"},Contours–>{0.001},PlotRange–>All,ContourStyleOpacity[0.25],ColorFunctionScalingTrue,BoxRatiosAutomatic]
Out[]=
In[]:=
Show[envMesh,graph]
Out[]=
Solid:
In[]:=
envRegion=RegionPlot3D[rho[w,{x,y,z}]>0.001,{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},AxesLabel{"x","y","z"},PlotRange–>All,PlotStyleOpacity[0.25],ColorFunctionScalingTrue,BoxRatiosAutomatic]
Out[]=
In[]:=
Show[envRegion,graph]
Out[]=
Zero Laplacian Surface:
In[]:=
ContourPlot3D[–laprho[w,{x,y,z}],{x,bbxmin,bbxmax},{y,bbymin,bbymax},{z,bbzmin,bbzmax},AxesLabel{"x","y","z"},Contours–>{0},PlotRange–>All,ColorFunction"PearlColors",ColorFunctionScalingTrue,BoxRatios–>Automatic]
Out[]=
Laplacian in a plane, lone pairs are maxima front/right.
In[]:=
Plot3D[–laprho[w,{0,y,z}],{y,bbymin,bbymax},{z,bbzmin,bbzmax},AxesLabel{"x","y","z"},PlotPoints–>50,PlotRange–>{Full,Full,{–1*^0,1*^0}},ColorFunction"PearlColors",ColorFunctionScalingTrue,BoxRatios–>Automatic]
Out[]=
Energy Decomposition over Atomic Basins
Energy Decomposition over Atomic Basins
w["Energy"]
E=Vnn+Vne+Vee+Tn+Te
Tetrahedrane
Tetrahedrane
Showcase tetrahedrane, where all possible critical point types are represented, as a demonstration of fairly comprehensive functionality:
Out[]=
Epilogue
Epilogue
In[]:=
Uninstall[link]ParallelEvaluate[Uninstall[link]]
Leave a Reply