sample.orbital.density.ais.parallel: Sample the wave function's density using AIS in parallel

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

Sample the wave function's density using AIS in parallel

Description

Sample the wave function's density using but now in parallel

Usage

sample.orbital.density.ais.parallel(
  nqn,
  lqn,
  mqn,
  orb.func.info,
  printQ = FALSE,
  num.processes = 1
)

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 annealed importance sampling (AIS). This fancy-er kind of MCMC is a little better at catching all the nodes in the wave function. Samples each chain in a separate process.

Value

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

References

Recipe adapted from https://wiseodd.github.io/techblog/2017/12/23/annealed-importance-sampling/

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.