pyvibdmc.simulation_utilities

Submodules

Classes

`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
`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
`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`

Functions

`get_atomic_num`(atms)
`get_atomic_string`(atomic_num)

Package Contents

class pyvibdmc.simulation_utilities.Constants[source]

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)[source]

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)[source]: 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)[source]: 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.simulation_utilities.get_atomic_num(atms)[source]

Parameters:: atms – A list (or single string) of atomic element symbols
Returns:: The atomic numbers of each of the atom strings you provide

pyvibdmc.simulation_utilities.get_atomic_string(atomic_num)[source]

Parameters:: atomic_num – The atomic numbers of each of the atom strings you provide
Returns:: A list of atomic element symbols

class pyvibdmc.simulation_utilities.XYZNPY[source]

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

static extract_xyz(file_name, num_atoms)[source]

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)[source]

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.simulation_utilities.Potential(potential_function, potential_directory, python_file, num_cores=1, pass_timestep=False, pot_kwargs=None)[source]

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)[source]: 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()[source]

class pyvibdmc.simulation_utilities.Potential_NoMP(potential_function, potential_directory, python_file, pass_timestep=False, ch_dir=False, pot_kwargs=None)[source]

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)[source]: 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.simulation_utilities.NN_Potential(potential_function, potential_directory, python_file, model, ch_dir=False, pot_kwargs=None, pass_timestep=False)[source]

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)[source]: Subclass of Potential_NoMP, but with an explicit argument passed for the NN model for convenience.

class pyvibdmc.simulation_utilities.Potential_Direct(potential_function, pot_kwargs=None, pass_timestep=False)[source]

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)[source]

class pyvibdmc.simulation_utilities.MolFiniteDifference[source]

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)[source]: 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)[source]

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

eq_geom

atoms

potential

masses = None

dx = 0.001

points_diag = 5

points_off_diag = 3

num_elems

dipole = None

generate_hessian()[source]

dipole_derivs()[source]

diagonalize(hessian)[source]

run()[source]

class pyvibdmc.simulation_utilities.InitialConditioner(coord, atoms, num_walkers, technique, technique_kwargs, masses=None)[source]

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)[source]

displace_along_nms(freqz, nmz, massez, ensemble)[source]

run_harm()[source]

run_permute()[source]: Must pass in a list of lists.

run()[source]

class pyvibdmc.simulation_utilities.DistIt(zs, method, eq_xyz=None, sorted_groups=None, sorted_atoms=None, full_mat=False, force_numpy=False)[source]

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()[source]

dist_matrix(atm_vec)[source]: If full matrix required, calculate it here.

atm_atm_dists(cds)[source]: Takes in coordinates, will return a vector of the atm-atm dists.

sort_atoms(d_mat)[source]

sort_groups(d_mat)[source]: 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)[source]

run(cds)[source]: Takes in cartesian coordinates, outputs the descriptor matrix in either vector form (upper triangle) or matrix form.

class pyvibdmc.simulation_utilities.Potential(potential_function, potential_directory, python_file, num_cores=1, pass_timestep=False, pot_kwargs=None)[source]

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)[source]: 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()[source]

class pyvibdmc.simulation_utilities.Potential_NoMP(potential_function, potential_directory, python_file, pass_timestep=False, ch_dir=False, pot_kwargs=None)[source]

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)[source]: 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.simulation_utilities.NN_Potential(potential_function, potential_directory, python_file, model, ch_dir=False, pot_kwargs=None, pass_timestep=False)[source]

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)[source]: Subclass of Potential_NoMP, but with an explicit argument passed for the NN model for convenience.

class pyvibdmc.simulation_utilities.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)[source]

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)[source]: Get trial wave function using multiprocessing

call_trial_no_mp(cds)[source]: 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)[source]: For when derivatives are not supplied, call finite difference function. This is still parallelized.

class pyvibdmc.simulation_utilities.ImpSampManager_NoMP(trial_function, trial_directory, python_file, chdir=False, pass_timestep=False, deriv_function=None, trial_kwargs=None, deriv_kwargs=None)[source]

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)[source]: Convenience function for trial, deriv, and sderiv so I don’t have to have triplicates of code

call_trial(cds)[source]: Call trial wave function.

call_derivs(cds)[source]: For when derivatives are not supplied, call finite difference function. Returns derivatives divided by psi already

class pyvibdmc.simulation_utilities.ImpSamp(imp_samp_manager)[source]

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)[source]: Internally returns the direct product wfn

drift(cds)[source]: 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)[source]

static local_kin(inv_masses_trip, sec_deriv)[source]

static finite_diff(cds, trial_func)[source]: Internal finite diff of needed derivatives for imp samp DMC, dpsi/dx and d2spi/dx2

class pyvibdmc.simulation_utilities.AnalyzeWfn(coordinates)[source]

xx

static dot_pdt(v1, v2)[source]: Takes two stacks of vectors and dots the pairs of vectors together

static exp_val(operator, dw)[source]

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)[source]

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)[source]

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)[source]

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)[source]: 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)[source]

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)[source]: 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)[source]

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)[source]

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.simulation_utilities.ChainRuleHelper(coords, module)[source]

cds

xp

analyzer

dpsidx(dpsi_dr, dr_dx)[source]: 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)[source]

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)[source]: 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)[source]: Calculates d2r/dx2 second derivatives for bond lengths

dcth_dx(atm_pair, cos_theta=None, dr_da=None, dr_dc=None)[source]: 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)[source]

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)[source]

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)[source]

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.