pyvibdmc

PyVibDMC A general purpose diffusion monte carlo code for studying vibrational problems

Submodules

• pyvibdmc.analysis
• pyvibdmc.pyvibdmc
• pyvibdmc.simulation_utilities

Attributes

versions

Classes

DMC_Sim	The main object that you will use to run the DMC simulation. You will instantiate an object DMC_SIM(...)
Constants	Thanks, Mark Boyer, for this silly little class.
XYZNPY	Handler of xyz <-> npy conversion, like a personal version of openBabel.
Potential	A potential handler that is able to call python functions that
Potential_NoMP	Version of Potential where no multiprocessing is included. As such, it does not leave a worker in the potential
NN_Potential	Subclass of Potential_NoMP, where a model argument is provided for convenience
Potential_Direct	Version of Potential where the actual Python function is passed rather than
MolFiniteDifference	Helper to calculate derivatives of some value as a function of Cartesian displacements in a molecule.
HarmonicAnalysis
InitialConditioner	If given a minimum energy geometry of the system you are trying to run, it will generate a preliminary ensemble.
DistIt	Calculates the Distance matrix, coulomb matrix, or SPF matrix (delta_r / r) for a given molecule. Returns a cupy or
ImpSampManager	Imports and Wraps around the user-provided trial wfn and (optionally) the first and second derivatives.
ImpSampManager_NoMP	Version of the manager that does not use any multiprocessing. If we ever evaluate the trial wfns with GPUs
ImpSamp	Internal class that:
AnalyzeWfn
ChainRuleHelper
AnalyzeWfn
SimInfo	Object that takes in a .hdf5 file (one of the outputs of the simulation) and provides tools for analysis.
Plotter	A very basic plotting class that will use matplotlib to generate various plots using results from PyVibDMC/AnalyzeWfn/SimInfo.
MolRotator	A helper class that will rotate a stack of molecules and generate 3D rotation matrices using vectorized

Functions

dmc_restart(potential, chkpt_folder, sim_name[, ...])	TODO: Need to add in impsamp infrastructure here for restarting...
get_atomic_num(atms)
get_atomic_string(atomic_num)

Package Contents

class pyvibdmc.DMC_Sim(sim_name='DMC_Sim', output_folder='exSimResults', weighting='discrete', num_walkers=10000, num_timesteps=20000, equil_steps=2000, chkpt_every=1000, wfn_every=1000, desc_wt_steps=100, atoms=[], delta_t=1, potential=None, masses=None, start_structures=None, start_cont_wts=None, branch_every=1, log_every=1, cur_timestep=0, cont_wt_thresh=None, imp_samp=None, imp_samp_oned=False, second_impsamp_displacement=False, excited_state_imp_samp=False, adiabatic_dmc=None, fixed_node=None, DEBUG_alpha=None, DEBUG_save_desc_wt_tracker=None, DEBUG_save_training_every=None, DEBUG_save_before_bod=False, DEBUG_mass_change=None)

The main object that you will use to run the DMC simulation. You will instantiate an object DMC_SIM(…) and then use object.run() to begin the simulation.

Parameters:

• sim_name (str) – Simulation name for saving wave functions and checkpointing

• output_folder (str) – The folder where the results will be stored, including wave functions and energies. Abs. or rel. path.

• weighting (str) – Discrete or Continuous weighting DMC. Continuous means that there are fixed number of walkers

• num_walkers (int) – Number of walkers we will start the simulation with

• num_timesteps (int) – Total time steps we will be propagating the walkers. num_timesteps*delta_t = total time in A.U.

• equil_steps (int) – Time steps before we start collecting wave function

• chkpt_every (int) – How many time steps in between checkpoints (only checkpoints once we reached equil_steps.

• desc_wt_steps (int) – Number of time steps monitoring the wave function’s desc_wt_timeendants.

• atoms (list) – List of atoms for the simulation

• delta_t (int) – The length of the time step; how many atomic units of time are you going in one time step.

• potential (function) – Takes in coordinates, gives back energies

• masses (list) – For feeding in artificial masses in atomic units. If not, then the atoms param will designate masses

• start_structures (np.ndarray) – An initial structure to initialize all your walkers

• branch_every (int) – Branch the walkers every x time steps, also known as birth/death

• log_every (int) – Gather simulation information every n time steps. This includes timing of potential call, max/min of weight and energy of ensemble, and number of births/deaths or branching

• cur_timestep (int) – The current time step, should be zero unless you are restarting

• cont_wt_thresh (list of floats) – If a weight goes past this bound, branch it. If only supplied one number, it will use that for the lower bound.

• DEBUG_alpha (float) – The number that will be used instead of 1/(2*delta_t) for alpha.

• DEBUG_save_desc_wt_tracker (bool) – If true, will save the array that keeps track of births/deaths during desc_wt_time weighting. Will save as a binary .npy file (see numpy documentation)

• DEBUG_save_training_every (int) – If true, will collect coordinates and energies every n time steps

• DEBUG_mass_change – Dictionary that will scale the mass every change_every steps. The scaling factor_per_change can be an array or an int/float

“type DEBUG_mass_change: dict

atoms = []

sim_name = 'DMC_Sim'

output_folder = 'exSimResults'

num_walkers = 10000

num_timesteps = 20000

potential_info

potential

weighting = ''

desc_wt_time_steps = 100

branch_every = 1

delta_t = 1

start_structures = None

start_cont_wts = None

masses = None

equil_steps = 2000

chkpt_every = 1000

wfn_every = 1000

log_every = 1

cur_timestep = 0

cont_wt_thresh = None

impsamp_manager = None

imp1d = False

second_impsamp_displacement = False

adiabatic_dmc = None

fixed_node = None

excited_state_imp_samp = False

property vref_vs_tau

Returns the vref array, including zeros from initialization

property walkers

update_effetive_timestep()

birth_or_death()

Chooses whether or not the walker made a bad enough random walk to be removed from the simulation. For discrete weighting, this leads to removal or duplication of the walkers. For continuous, this leads to an update of the weights and a potential branching of a large weight walker to the smallest one :return: Updated Continus Weights , the “who from” array for desc_wt_timeendent weighting, walker coords, and pot vals.

birth_or_death_old()

Chooses whether or not the walker made a bad enough random walk to be removed from the simulation. For discrete weighting, this leads to removal or duplication of the walkers. For continuous, this leads to an update of the weights and a potential branching of a large weight walker to the smallest one :return: Updated Continus Weights , the “who from” array for desc_wt_timeendent weighting, walker coords, and pot vals.

move_randomly()

The random displacement of each of the coordinates of each of the walkers, done in a vectorized fashion. Displaces self._walker_coords

imp_move_randomly()

The random displacement of each of the coordinates of each of the walkers, done in a vectorized fashion. Displaces self._walker_coords

imp_move_randomly_second_type()

The random displacement of each of the coordinates of each of the walkers, done in a vectorized fashion. Displaces self._walker_coords

calc_vref()

Use the energy of all walkers to calculate self._vref with a correction for the fluctuation in the population or weight. Updates Vref

calc_desc_wts()

At the end of desc_wt_timeendent weighting, count up which walkers came from other walkers (desc_wt_timeendants)

update_sim_arrays(prop_step)

Parameters:

prop_step – the time step that we are currently on

recrossing(q_1, q_2, cds)

propagate()

The main DMC loop. 1. Move Randomly 2. Calculate the Potential Energy 3. Birth/Death 4. Update Vref Additionally, checks when the wavefunction has hit a point where it should save / do desc_wt_timeendent weighting.

run()

This function calls propagate and saves simulation results

pyvibdmc.dmc_restart(potential, chkpt_folder, sim_name, additional_timesteps=0, impsamp=None, imp_samp_oned=False, fixed_node=None)

TODO: Need to add in impsamp infrastructure here for restarting…

class pyvibdmc.Constants

Thanks, Mark Boyer, for this silly little class. Converter that handles energy, distance, and mass conversions for DMC. Can be expanded upon.

atomic_units

classmethod convert(val, unit, to_AU=True)

Parameters:

• val (np.ndarray) – The value or values that will be converted

• unit (str) – The units (not atomic units) that we will be converting to or from

• to_AU – If true, converting from non-a.u. to a.u. If false, converting to a.u. from non-a.u.

:type to_AU:boolean :return: converted values

classmethod mass(atom, to_AU=True)

Given a string that corresponds to an atomic element, output the atomic mass of that element :param atom: The string of an atomic element :type atom:str :param to_AU: If true, converting from non-a.u. to a.u. If false, converting to a.u. from non-a.u. :type to_AU:boolean :return: mass in atomic units unless user changes to_AU to False, then AMU

classmethod reduced_mass(atoms, to_AU=True)

Given a string like ‘O-H’ or ‘N-N’ , output the reduced mass of that diatomic :param atoms: A string that is composed of two atoms :type atom:str :param to_AU: If true, converting from non-a.u. to a.u. If false, converting to a.u. from non-a.u. :type to_AU:boolean :return: mass in atomic units unless user changes to_AU to False, then AMU

pyvibdmc.get_atomic_num(atms)

Parameters:

atms – A list (or single string) of atomic element symbols

Returns:

The atomic numbers of each of the atom strings you provide

pyvibdmc.get_atomic_string(atomic_num)

Parameters:

atomic_num – The atomic numbers of each of the atom strings you provide

Returns:

A list of atomic element symbols

class pyvibdmc.XYZNPY

Handler of xyz <-> npy conversion, like a personal version of openBabel.

static extract_xyz(file_name, num_atoms)

Extracts the coordinates from an xyz file and returns it as an np array of dimension nxmx3, where n = number of geometries, m = number of atoms, and 3 = cartesian coordinates

Parameters:

• file_name (str) – Name of the .xyz file, including extension

• num_atoms (int) – Number of atoms in molecule.

Returns:

nxmx3 numpy array of coordinates

static write_xyz(coords, fname, atm_strings, cmt=None)

Writes a numpy array of x,y,z coordinates to a .xyz file

Parameters:

• fname – name of xyz file

• coords – numpy array, either mx3 or nxmx3, where n = number of geometries and m = number of atoms

• atm_strings – list of strings that correspond to the atom type e.g. [“H”,”H”,”O”]

Returns:

np.ndarray

class pyvibdmc.Potential(potential_function, potential_directory, python_file, num_cores=1, pass_timestep=False, pot_kwargs=None)

A potential handler that is able to call python functions that call .so files, either generated by f2py or loaded in by ctypes. :param potential_function: The name of a python function (user specified) that will take in a n x m x 3 stack of geometries and return a 1D numpy array filled with potential values in hartrees. :type potential_function: str :param potential_dir: The absolute path to the directory that contains the .so file and .py file. If it”s a python function, then just the absolute path to your .py file. :type: str :param num_cores: Will create a pool of <num_cores> processes using Python”s multiprocessing module. This should never be larger than the number of processors on the machine this code is run. :type: int

potential_function

python_file

potential_directory

num_cores = 1

pass_timestep = False

pot_kwargs = None

property pool

Returns the potential manager’s pool so that it can be used internally with Imp Samp or with the user elsewhere

getpot(cds, timeit=False)

Uses the potential function we got to call potential :param cds: A stack of geometries (nxmx3, n=num_geoms;m=num_atoms;3=x,y,z) whose energies we need :type cds: np.ndarray :param timeit: The logger telling the potential manager whether or not to time the potential call :type timeit: bool

mp_close()

class pyvibdmc.Potential_NoMP(potential_function, potential_directory, python_file, pass_timestep=False, ch_dir=False, pot_kwargs=None)

Version of Potential where no multiprocessing is included. As such, it does not leave a worker in the potential directory of interest. If you need to cd into the directory to call a potential (for example, the potential loads in a file using a relative path or something), then use ch_dir=True. Not the default since it is computationally inefficient to cd during each potential call when not necessary.

potential_function

python_file

pass_timestep = False

potential_directory

pot_kwargs = None

ch_dir = False

getpot(cds, timeit=False)

Uses the potential function we got to call potential :param cds: A stack of geometries (nxmx3, n=num_geoms;m=num_atoms;3=x,y,z) whose energies we need :type cds: np.ndarray :param timeit: The logger telling the potential manager whether or not to time the potential call :type timeit: bool

class pyvibdmc.NN_Potential(potential_function, potential_directory, python_file, model, ch_dir=False, pot_kwargs=None, pass_timestep=False)

Bases: Potential_NoMP

Subclass of Potential_NoMP, where a model argument is provided for convenience

model

ch_dir = False

pot_kwargs = None

getpot(cds, timeit=False)

Subclass of Potential_NoMP, but with an explicit argument passed for the NN model for convenience.

class pyvibdmc.Potential_Direct(potential_function, pot_kwargs=None, pass_timestep=False)

Version of Potential where the actual Python function is passed rather than being imported from an external file. At the request of Mark.

potential_function

pot_kwargs = None

pass_timestep = False

getpot(cds, timeit=False)

class pyvibdmc.MolFiniteDifference

Helper to calculate derivatives of some value as a function of Cartesian displacements in a molecule.

weights_der1

weights_der2

static displace_molecule(eq_geom, atm_cd, dx, num_disps)

Displace atm along cd :param eq_geom: Geometry from which you will be displaced :param atm_cd: a tuple that has a paricular atom of interest to displace in a particular dimension (x,y,or,z) :param dx: The amount each geometry will be replaced :param num_disps: int of how many displacements to do, will take the form of 3 or 5 for harmonic analysis. :return: Displaced coordinates in a 3D array (n,m,3). If displaced in two directions, then still (n,m,3)

classmethod differentiate(values, dx, num_points, der)

class pyvibdmc.HarmonicAnalysis(eq_geom, atoms, potential, masses=None, dx=0.001, points_diag=5, points_off_diag=3, dipole=None)

eq_geom

atoms

potential

masses = None

dx = 0.001

points_diag = 5

points_off_diag = 3

num_elems

dipole = None

generate_hessian()

dipole_derivs()

diagonalize(hessian)

run()

class pyvibdmc.InitialConditioner(coord, atoms, num_walkers, technique, technique_kwargs, masses=None)

If given a minimum energy geometry of the system you are trying to run, it will generate a preliminary ensemble. In one instance, you can calculate the harmonic frequencies and normal modes, and then sample the harmonic 3N-6 ground state wave function along those normal modes. Will, in the future, handle other initial conditions.

coord

atoms

masses = None

num_walkers

technique

technique_kwargs

gen_disps(sigmas)

displace_along_nms(freqz, nmz, massez, ensemble)

run_harm()

run_permute()

Must pass in a list of lists.

run()

class pyvibdmc.DistIt(zs, method, eq_xyz=None, sorted_groups=None, sorted_atoms=None, full_mat=False, force_numpy=False)

Calculates the Distance matrix, coulomb matrix, or SPF matrix (delta_r / r) for a given molecule. Returns a cupy or numpy multidimensional array. If cupy is installed, this will default to using it. This means in the run(cds) function, you must pass cds as a cupy array

zs

method

eq_xyz = None

sorted_groups = None

sorted_atoms = None

full_mat = False

force_numpy = False

check_cupy()

dist_matrix(atm_vec)

If full matrix required, calculate it here.

atm_atm_dists(cds)

Takes in coordinates, will return a vector of the atm-atm dists.

sort_atoms(d_mat)

sort_groups(d_mat)

Swap 2 or more groups of atoms according to the sum of the norm of the group’s columns. Only one type of group can swapped.

get_prepped_vec(atm_atm_vec)

run(cds)

Takes in cartesian coordinates, outputs the descriptor matrix in either vector form (upper triangle) or matrix form.

class pyvibdmc.ImpSampManager(trial_function, trial_directory, python_file, pot_manager, pass_timestep=False, new_pool_num_cores=None, deriv_function=None, trial_kwargs=None, deriv_kwargs=None)

Imports and Wraps around the user-provided trial wfn and (optionally) the first and second derivatives. Parallelized using multiprocessing, which is considered the default for pyvibdmc.

trial_func

trial_dir

python_file

deriv_func = None

trial_kwargs = None

deriv_kwargs = None

pot_manager

pass_timestep = False

nomp_pool_cores = None

call_trial(cds)

Get trial wave function using multiprocessing

call_trial_no_mp(cds)

For call_derivs (finite diff), get trial wave function. Still used in the mp.pool context, just doesn’t call pool itself

call_derivs(cds)

For when derivatives are not supplied, call finite difference function. This is still parallelized.

class pyvibdmc.ImpSampManager_NoMP(trial_function, trial_directory, python_file, chdir=False, pass_timestep=False, deriv_function=None, trial_kwargs=None, deriv_kwargs=None)

Version of the manager that does not use any multiprocessing. If we ever evaluate the trial wfns with GPUs this could be useful. Could also be useful if multiprocessing is incompatible with your workflow.

trial_fuc

trial_dir

python_file

pass_timestep = False

deriv_func = None

trial_kwargs = None

deriv_kwargs = None

chdir = False

call_imp_func(func, cds, func_kwargs=None)

Convenience function for trial, deriv, and sderiv so I don’t have to have triplicates of code

call_trial(cds)

Call trial wave function.

call_derivs(cds)

For when derivatives are not supplied, call finite difference function. Returns derivatives divided by psi already

class pyvibdmc.ImpSamp(imp_samp_manager)

Internal class that: 1. Calculates local energy 2. Calculates drift terms 3. Calls trial wave function Directly interfaces with the simulation. Also a good place to put finite difference derivatives.

imp_manager

trial(cds)

Internally returns the direct product wfn

drift(cds)

Internally returns (dpsi/psi), since it’s more convenient in the workflow. Also returns second derivatives divided by psi

static metropolis(sigma_trip, trial_x, trial_y, disp_x, disp_y, D_x, D_y, dt)

static local_kin(inv_masses_trip, sec_deriv)

static finite_diff(cds, trial_func)

Internal finite diff of needed derivatives for imp samp DMC, dpsi/dx and d2spi/dx2

class pyvibdmc.AnalyzeWfn(coordinates)

xx

static dot_pdt(v1, v2)

Takes two stacks of vectors and dots the pairs of vectors together

static exp_val(operator, dw)

Calculation of an expectation value using Monte Carlo Integration

Parameters:

• operator (np.ndarray) – array of bond lengths, bond angles, potentials, dipole moments, etc.

• dw (np.ndarray) – Descendant Weights

Returns:

Expectation value (single float)

bond_length(atm1, atm2)

Computes the norm of the vector between two atoms.

Parameters:

• atm1 (int) – desired atom number 1 (python index starts at 0)

• atm2 (int) – desired atom number 2

Returns:

norm of x[atm1]-x[atm2]

bond_angle(atm1, atm_vert, atm3)

Calculate atm1 - atm_vert - atm3 angle

Parameters:

• atm1 (int) – Index of the first external atom (python index starts at 0)

• atm_vert (int) – Index of the atom at the vertex (python index starts at 0)

• atm3 (int) – Index of the second external atom (python index starts at 0)

Returns:

Angle in radians

bisecting_vector(atm1, atm_vert, atm3)

Get normalized bisecting vector of two vectors (two bond lengths). It is assumed that the two flanking atoms have the same central atom. |b|*a + |a|*b

Parameters:

• atm1 (int) – Index of the first external atom (python index starts at 0)

• atm_vert (int) – Index of the atom at the vertex (python index starts at 0)

• atm3 (int) – Index of the second external atom (python index starts at 0)

Returns:

normalized bisecting vector

dihedral(atm_1, atm_2, atm_3, atm_4)

Defines the dihedral coordinate using 3 vectors. vec_1 and vec_2 form a plane, and vec_2 and vec_3 form a plane. The dihedral angle is the calculated angle betwen these two planes. The sign is dependent on the direction of the vectors, as the cross product is taken between the two pairs of two vectors https://stackoverflow.com/questions/20305272/dihedral-torsion-angle-from-four-points-in-cartesian-coordinates-in-python https://en.wikipedia.org/w/index.php?title=Dihedral_angle&oldid=689165217#Angle_between_three_vectors

get_component(atm, xyz)

Get x, y, or z component of a vector that corresponds to a particular atom in some predetermined cooridnate system.

Parameters:

• atm – atom’s index in numpy array

• xyz (int) – Either 0 (x), 1 (y), or 2 (z)

Type:

int

Returns:

vector of x, y, or z components of a stack of vectors

static cart_to_spherical(vecc)

Takes stack of cartesian vectors (nx3) and translates them to stack of r theta phi coordinates, with the assumption that you are already in the coordinate frame you desire.

static projection_1d(attr, desc_weights, bins=25, range=None, normalize=True)

Project the probability amplitude onto a particular coordinate, 1D Histogram.

Parameters:

• attr (bool) – the coordinate to be projected onto (bond length, bond angle, etc.)

• desc_weights (np.ndarray) – the descendant weights for psi**2

• bins – Number of bins in histogram (higher number, more number of walkers needed)

• range (tuple) – Range over which attr is histogrammed

• normalize – Normalizes the wave function along coordinate

Returns:

np.ndarray of shape (bin_num-1 x 2), bin centers, amplitude at bin centers

static projection_2d(attr1, attr2, desc_weights, bins=[25, 25], range=None, normalize=True)

Project the probability amplitude onto a 2 coordinates, 2D Histogram.

Parameters:

• attr1 – coordinate 1 to be projected onto (bond length, bond angle, etc.)

• attr2 – coordinate 2 to be projected onto (bond length, bond angle, etc.)

• desc_weights (np.ndarray) – the descendant weights for psi**2

• bins – Number of bins for histogram in both directions (higher number, more number of walkers needed)

• range (list of two tuples) – Range over which attr is histogrammed

• normalize – Normalizes the wave function along both coordinates (integrated area under 2D surface is 1)

Returns:

np.ndarray of shape (bin_num-1 x 2), bin centers x, bin centers y, and then the 2D array with amplitude.

class pyvibdmc.ChainRuleHelper(coords, module)

cds

xp

analyzer

dpsidx(dpsi_dr, dr_dx)

Generic function that takes in a series of dpsi/dr matrices and dr/dx matrices and generates the dpsi/dx matrix. Assumes direct product wave function fed in appropriately.

d2psidx2(d2psi_dr2, d2r_dx2, dpsi_dr, dr_dx)

Generic function that takes in a series of d2psi_dr2 and d2r/dx matrices and generates the second derivative of

psi wrt x, doing out the chain rule explicitly.

dr_dx(atm_pair)

calculates the dr/dx first derivatives for bond lengths :param self.cds: num_walkers x num_atms x 3 array :param atm_bonds: list of bond lengths that are relevant to the trial wave function. like [0,1] :return: num_walkers x num_atoms x 3 array of dr/dx

d2r_dx2(atm_pair, dr_dx=None)

Calculates d2r/dx2 second derivatives for bond lengths

dcth_dx(atm_pair, cos_theta=None, dr_da=None, dr_dc=None)

calculates the dcos(theta)/dx first derivatives for bond lengths :param self.cds: num_walkers x num_atms x 3 array :param atm_pair: list of bond lengths that are relevant to the trial wave function. like [0,2,1]. ALWAYS put vertex of angle in middle index atm_bonds[x][1]. :param dr_das: The derivative of the bond length corresponding to atm_bonds[i][0] and atm_bonds[i][1]. i is the index of the atm_bonds list of lists :param dr_dcs: The derivative of the bond length corresponding to atm_bonds[i][1] and atm_bonds[i][2]. i is the index of the atm_bonds list of lists :return: len(atm_bonds) x num_walkers x num_atoms x 3 array of dr/dx

d2cth_dx2(atm_pair, cos_theta=None, dr_da=None, dr_dc=None, d2r_da2=None, d2r_dc2=None)

Parameters:

• self.cds – num_walkers x num_atms x 3 array

• atm_bonds – list of lists of bond lengths that are relevant to the trial wave function. like [[0,2,1]]. ALWAYS put vertex of angle in middle index atm_bonds[x][1].

• dr_das – The derivative of the bond length corresponding to atm_bonds[i][0] and atm_bonds[i][1]. i is the index of the atm_bonds list of lists

• dr_dcs – The derivative of the bond length corresponding to atm_bonds[i][1] and atm_bonds[i][2]. i is the index of the atm_bonds list of lists

• d2r_da2s – The second derivative of the bond length corresponding to atm_bonds[i][0] and atm_bonds[i][1]. i is the index of the atm_bonds list of lists

• d2r_dc2s – The second derivative of the bond length corresponding to atm_bonds[i][1] and atm_bonds[i][2]. i is the index of the atm_bonds list of lists

dth_dx(atm_pair, cos_theta=None, dcth_dx=None, dr_da=None, dr_dc=None)

Parameters:

• self.cds – num_walkers x num_atoms x 3 array

• atm_pair – list that corresponds to indices of the ABC angle, where B is the vertex

• dcth_dxs – Derivative of cosine theta wrt x. This is needed for this derivative and the second derivative.

d2th_dx2(atm_pair, cos_theta=None, dcth_dx=None, dr_da=None, dr_dc=None, d2r_da2=None, d2r_dc2=None)

Parameters:

• self.cds – num_walkers x num_atoms x 3 array

• atm_pair – list that corresponds to indices of the ABC angle, where B is the vertex

• dcth_dxs – Derivative of cosine theta wrt x. This is needed for this derivative and the second derivative.

class pyvibdmc.AnalyzeWfn(coordinates)

xx

static dot_pdt(v1, v2)

Takes two stacks of vectors and dots the pairs of vectors together

static exp_val(operator, dw)

Calculation of an expectation value using Monte Carlo Integration

Parameters:

• operator (np.ndarray) – array of bond lengths, bond angles, potentials, dipole moments, etc.

• dw (np.ndarray) – Descendant Weights

Returns:

Expectation value (single float)

bond_length(atm1, atm2)

Computes the norm of the vector between two atoms.

Parameters:

• atm1 (int) – desired atom number 1 (python index starts at 0)

• atm2 (int) – desired atom number 2

Returns:

norm of x[atm1]-x[atm2]

bond_angle(atm1, atm_vert, atm3)

Calculate atm1 - atm_vert - atm3 angle

Parameters:

• atm1 (int) – Index of the first external atom (python index starts at 0)

• atm_vert (int) – Index of the atom at the vertex (python index starts at 0)

• atm3 (int) – Index of the second external atom (python index starts at 0)

Returns:

Angle in radians

bisecting_vector(atm1, atm_vert, atm3)

Get normalized bisecting vector of two vectors (two bond lengths). It is assumed that the two flanking atoms have the same central atom. |b|*a + |a|*b

Parameters:

• atm1 (int) – Index of the first external atom (python index starts at 0)

• atm_vert (int) – Index of the atom at the vertex (python index starts at 0)

• atm3 (int) – Index of the second external atom (python index starts at 0)

Returns:

normalized bisecting vector

dihedral(atm_1, atm_2, atm_3, atm_4)

Defines the dihedral coordinate using 3 vectors. vec_1 and vec_2 form a plane, and vec_2 and vec_3 form a plane. The dihedral angle is the calculated angle betwen these two planes. The sign is dependent on the direction of the vectors, as the cross product is taken between the two pairs of two vectors https://stackoverflow.com/questions/20305272/dihedral-torsion-angle-from-four-points-in-cartesian-coordinates-in-python https://en.wikipedia.org/w/index.php?title=Dihedral_angle&oldid=689165217#Angle_between_three_vectors

get_component(atm, xyz)

Get x, y, or z component of a vector that corresponds to a particular atom in some predetermined cooridnate system.

Parameters:

• atm – atom’s index in numpy array

• xyz (int) – Either 0 (x), 1 (y), or 2 (z)

Type:

int

Returns:

vector of x, y, or z components of a stack of vectors

static cart_to_spherical(vecc)

Takes stack of cartesian vectors (nx3) and translates them to stack of r theta phi coordinates, with the assumption that you are already in the coordinate frame you desire.

static projection_1d(attr, desc_weights, bins=25, range=None, normalize=True)

Project the probability amplitude onto a particular coordinate, 1D Histogram.

Parameters:

• attr (bool) – the coordinate to be projected onto (bond length, bond angle, etc.)

• desc_weights (np.ndarray) – the descendant weights for psi**2

• bins – Number of bins in histogram (higher number, more number of walkers needed)

• range (tuple) – Range over which attr is histogrammed

• normalize – Normalizes the wave function along coordinate

Returns:

np.ndarray of shape (bin_num-1 x 2), bin centers, amplitude at bin centers

static projection_2d(attr1, attr2, desc_weights, bins=[25, 25], range=None, normalize=True)

Project the probability amplitude onto a 2 coordinates, 2D Histogram.

Parameters:

• attr1 – coordinate 1 to be projected onto (bond length, bond angle, etc.)

• attr2 – coordinate 2 to be projected onto (bond length, bond angle, etc.)

• desc_weights (np.ndarray) – the descendant weights for psi**2

• bins – Number of bins for histogram in both directions (higher number, more number of walkers needed)

• range (list of two tuples) – Range over which attr is histogrammed

• normalize – Normalizes the wave function along both coordinates (integrated area under 2D surface is 1)

Returns:

np.ndarray of shape (bin_num-1 x 2), bin centers x, bin centers y, and then the 2D array with amplitude.

class pyvibdmc.SimInfo(h5_name)

Object that takes in a .hdf5 file (one of the outputs of the simulation) and provides tools for analysis.

Parameters:

h5Name – hdf5 file

:type str

fname

static get_wfn(wfn_fl, ret_ang=False, get_parent_wts=False)

Given a .hdf5 file, return wave function and descendant weights associated with that wave function.

Parameters:

• wfn_fl – A resultant .hdf5 file from a PyVibDMC simulation

• ret_ang – boolean indicating returning the coordinates in angtstroms. Bohr is the default.

Returns:

Coordinates array in angstroms (nxmx3), descendant weights array (n).

get_wfns(time_step_list, ret_ang=False, get_parent_wts=False)

Extract the wave function (walker set) and descendant weights given a time step number or numbers :param time_step_list: a list of ints that correspond to the time steps you want the wfn from given the simulation you are working with :type time_step_list: int or list :param ret_ang: boolean indicating returning the coordinates in angtstroms. Bohr is the default. :param get_parent_wts: Return the continuous weights associated with the walkers at the beginning of descendant weighting.

get_vref(ret_cm=False)

Returns vref_vs_tau array

get_pop()

Returns population array, either ensemble size or sum of weights

get_atomic_nums()

Returns list of atoms used in the simulation (by atomic number)

get_atom_masses()

Returns masses used in the simulation in atomic units (mass electron)

get_zpe(onwards=1000, ret_cm=False)

onwards is an int that tells us where to start averaging (python indexing starts at 0)

window_avg(blocks=5, ret_cm=False)

Splits vref into blocks, calculates zpe in each of those blocks

static get_training(training_file, ret_ang=False, ret_cm=False)

If using deb_training_every argument, read the files with this. Returns walkers in angstr and engs in cm-1

class pyvibdmc.Plotter

A very basic plotting class that will use matplotlib to generate various plots using results from PyVibDMC/AnalyzeWfn/SimInfo.

params

static plt_vref_vs_tau(vref_vs_tau, save_name='vref_vs_tau.png')

Takes in the vref vs tau array from a DMC sim and plots it. Can also take in many vref_vs_tau arrays and plot them successively

Parameters:

• vref_vs_tau (str or list) – The vref_vs_tau array from the *sim_info.hdf5 file. Can be a list of these as well.

• save_name (str) – The name of the plot that will save as a .png

Returns:

static plt_pop_vs_tau(pop_vs_tau, save_name='pop_vs_tau.png')

Takes in the pop_vs_tau array from a DMC sim and plots it. Can also take in many pop_vs_tau arrays and plot them successively

Parameters:

• pop_vs_tau (str or list) – The vref_vs_tau array from the *sim_info.hdf5 file. Can be a list of these as well.

• save_name (str) – The name of the plot that will save as a .png

Returns:

static plt_hist1d(hist, xlabel, ylabel='Probability Amplitude ($\\mathrm{\\Psi^{2}}$)', title='', save_name='histogram.png')

Plots the histogram generated from AnalyzeWfn.projection_1d.

Parameters:

• hist (np.ndarray) – Output from AnalyzeWfn.projection_1d ; array of shape (bin_num-1 x 2)

• xlabel (str) – What to label the x-axis

• ylabel (str) – What to label the y-axis (units in Parenthesis is always good)

• title (str) – Title of the plot

• save_name (str) – name of the .png file that this plot will be saved to

static plt_hist2d(binsx, binsy, hist_2d, xlabel, ylabel, title='', save_name='histogram.png')

Plots the 2D-histogram generated from AnalyzeWfn.projection_2d.

Parameters:

• bins (np.ndarray) – First output from AnalyzeWfn.projection_2d; the array containing the x and y bins

• hist_2d (np.ndarray) – Second output from AnalyzeWfn.projection_2d; The matrix that contains the amplitude at each of the bins in the 2d histogram.

• xlabel (str) – What to label the x-axis

• ylabel (str) – What to label the y-axis (units in Parenthesis is always good)

• title (str) – Title of the plot

• save_name (str) – name of the .png file that this plot will be saved to

class pyvibdmc.MolRotator

A helper class that will rotate a stack of molecules and generate 3D rotation matrices using vectorized numpy operations.

static rotate_geoms(rot_mats, geoms)

Takes in a stack of rotation matrices and applies it to a stack of geometries.

static rotate_vec(rot_mats, vecc)

Takes in a stack of rotation matrices and applies it to a stack of vector

static gen_rot_mats(theta, xyz_int)

Generates the 3D rotation matrix about X, Y, or Z by theta radians

static rotate_to_xy_plane(geoms, orig, xax, xyp, return_mat=False)

Rotate geometries to XY plane, placing one atom at the origin, one on the xaxis, and one on the xyplane. Done through successive rotations about the X-axis, Z-axis then X-axis again.

static gen_eulers(xyz, XYZ)

Takes in cartesian vectors and gives you the 3 euler angles that bring xyz to XYZ based on a ‘ZYZ’ rotation

static extract_eulers(rot_mats)

From a rotation matrix, calculate the three euler angles theta,phi and Chi. This is based on a ‘ZYZ’ euler rotation

pyvibdmc.versions