sample.orbital.density.mcmc: MCMC sample orbital density

In npetraco/che302r: Handy functions for CHEM 302, Physical Chemistry II

MCMC sample orbital density

Description

Sample the orbital wave function's density by simple Metropolis Markov Chain Monte Carlo

Usage

sample.orbital.density.mcmc(nqn, lqn, mqn, orb.func.info, printQ = FALSE)

Arguments

nqn	n quantum number
lqn	l quantum number
mqn	m quantum number
orb.func.info	Orbital pdfs info generated by splined.sqorb.functs()

Details

Sample the wave function's density and plot corresponding r,theta, phi (x,y,z) values. Uses MCMC. Under the hood, this is a VERY stripped down wrapper (to a wrapper: sample.density3b) to run MCMC. CHECK THAT YOU GET ALL THE NODES IN YOUR WAVE FUNCTIONS! MCMC by itself can easily miss them.

Value

orbital sample matrix (num.sample-rows by 3-columns)

References

XXXXX

Examples

library(che302r)
library(parallel)
library(doSNOW)
library(rgl)

#Pick an orbital:
n <- 1
l <- 0
m <- (0)

# Create squared orbital functions for sampling
orb.sq.funcs  <- splined.sqorb.functs(n, l, m, r.max = 320, num.knots = 1000)
rax           <- orb.sq.funcs$r.axis
R.sq          <- orb.sq.funcs$R.sq.func
plot(rax, rax^2 * R.sq(rax), typ="l")

#Sample the orbital density:
#detectCores() # Number of processors available

#density.samples <- sample.orbital.density.ais(n, l, m, printQ = T, orb.func.info = orb.sq.funcs)
density.samples <- sample.orbital.density.mcmc(n, l, m, printQ = T, orb.func.info = orb.sq.funcs)
#density.samples <- sample.orbital.density.ais.parallel(n, l, m, printQ = T, orb.func.info = orb.sq.funcs, num.processes = 2)
#density.samples <- sample.orbital.density.mcmc.parallel(n, l, m, printQ = T, orb.func.info = orb.sq.funcs, num.processes = 4)

# Sample in cartesian coordinates
x               <- density.samples[,1]
y               <- density.samples[,2]
z               <- density.samples[,3]

#Plot the sampled orbital density:
#2D
plot(x,y)
plot(x,z)
plot(y,z)

#3D
plot3d(x,y,z,type="s",radius=1.5,xlab="x",ylab="y",zlab="z")

# 3D slices
plot3d(x,y,rep(0,length(z)),type="s",radius=0.5,xlab="x",ylab="y",zlab="z")
plot3d(x,rep(0,length(y)),z,type="s",radius=0.5,xlab="x",ylab="y",zlab="z")
plot3d(rep(0,length(x)), y, z,type="s",radius=0.5,xlab="x",ylab="y",zlab="z")

npetraco/che302r documentation built on April 17, 2025, 10:34 p.m.