Monte Carlo methods¶
Variational Monte Carlo¶
-
pyqmc.mc.initial_guess(mol, nconfig, r=1.0)[source]¶ Generate an initial guess by distributing electrons near atoms proportional to their charge.
assign electrons to atoms based on atom charges assign the minimum number first, and assign the leftover ones randomly this algorithm chooses atoms with replacement to assign leftover electrons
Parameters: - mol – A PySCF-like molecule object. Should have atom_charges(), atom_coords(), and nelec
- nconfig – How many configurations to generate.
- r – How far from the atoms to distribute the electrons
Returns: (nconfig,nelectrons,3) array of electron positions randomly distributed near the atoms.
Return type: ndarray
-
pyqmc.mc.limdrift(g, cutoff=1)[source]¶ Limit a vector to have a maximum magnitude of cutoff while maintaining direction
Parameters: - g – a [nconf,ndim] vector
- cutoff – the maximum magnitude
Returns: The vector with the cutoff applied.
-
pyqmc.mc.vmc(wf, configs, nblocks=10, nsteps_per_block=10, nsteps=None, tstep=0.5, accumulators=None, verbose=False, stepoffset=0, hdf_file=None, continue_from=None, client=None, npartitions=None)[source]¶ Run a Monte Carlo sample of a given wave function.
Parameters: - wf (a PyQMC wave-function-like class) – trial wave function for VMC
- configs (PyQMC configs object) – Initial electron coordinates
- nblocks (int) – Number of VMC blocks to run
- nsteps_per_block (int) – Number of steps to run per block
- nsteps (int) – (Deprecated) Number of steps to run, maps to nblocks = 1, nsteps_per_block = nsteps
- tstep (float) – Time step for move proposals. Only affects efficiency.
- accumulators – A dictionary of functor objects that take in (configs,wf) and return a dictionary of quantities to be averaged. np.mean(quantity,axis=0) should give the average over configurations. If None, then the coordinates will only be propagated with acceptance information.
- verbose (boolean) – Print out step information
- stepoffset (int) – If continuing a run, what to start the step numbering at.
- hdf_file (str) – Hdf_file to store vmc output.
- client – an object with submit() functions that return futures
- npartitions (int) – the number of workers to submit at a time
Returns: (df,configs) df: A list of dictionaries nstep long that contains all results from the accumulators. These are averaged across all walkers.
configs: The final coordinates from this calculation.
Return type: list of dictionaries, pyqmc.coord.Configs
Wave function optimization¶
Evaluates accumulator on the same set of configs for correlated sampling of different wave function parameters
Parameters: - wf – wave function object
- configs – (nconf, nelec, 3) array
- params – (nsteps, nparams) array list of arrays of parameters (serialized) at each step
- pgrad_acc – PGradAccumulator
Returns: a single dict with indices [parameter, values]
-
pyqmc.linemin.line_minimization(wf, coords, pgrad_acc, steprange=0.2, max_iterations=30, warmup_options=None, vmcoptions=None, lmoptions=None, update=<function sr_update>, update_kws=None, verbose=False, npts=5, hdf_file=None, client=None, npartitions=None)[source]¶ - Optimizes energy by determining gradients with stochastic reconfiguration
- and minimizing the energy along gradient directions using correlated sampling.
Parameters: - wf – initial wave function
- coords – initial configurations
- pgrad_acc – A PGradAccumulator-like object
- steprange (float) – How far to search in the line minimization
- warmup (int) – number of steps to use for vmc warmup
- max_iterations (int) – (maximum) number of steps in the gradient descent
- vmcoptions (dict) – a dictionary of options for the vmc method
- lmoptions (dict) – a dictionary of options for the lm method
- update – A function that generates a parameter change
- update_kws – Any keywords
- npts (int) – number of points to fit to in each line minimization
- verbose (boolean) – print output if True
Returns: optimized wave function
Diffusion Monte Carlo¶
-
pyqmc.dmc.branch(configs, weights)[source]¶ Perform branching on a set of walkers using the ‘stochastic comb’
Walkers are resampled with probability proportional to the weights, and the new weights are all set to be equal to the average weight.
Parameters: - configs – (nconfig,nelec,3) walker coordinates
- weights – (nconfig,) walker weights
Returns: resampled walker configurations and weights all equal to average weight
-
pyqmc.dmc.dmc_propagate(wf, configs, weights, tstep, branchcut_start, e_trial, e_est, nsteps=5, accumulators=None, ekey=('energy', 'total'))[source]¶ Propagate DMC without branching
Parameters: - wf – A Wave function-like class. recompute(), gradient(), and updateinternals() are used, as well as anything (such as laplacian() ) used by accumulators
- configs – Configs object, (nconfig, nelec, 3) - initial coordinates to start calculation.
- weights – (nconfig,) - initial weights to start calculation
- tstep – Time step for move proposals. Introduces time step error.
- nsteps – number of DMC steps to take
- accumulators – A dictionary of functor objects that take in (coords,wf) and return a dictionary of quantities to be averaged. np.mean(quantity,axis=0) should give the average over configurations. If none, a default energy accumulator will be used.
- ekey – tuple of strings; energy is needed for DMC weights. Access total energy by accumulators[ekey[0]](configs, wf)[ekey[1]
Returns: (df,coords,weights) df: A list of dictionaries nstep long that contains all results from the accumulators.
coords: The final coordinates from this calculation.
weights: The final weights from this calculation
-
pyqmc.dmc.dmc_propagate_parallel(wf, configs, weights, client, npartitions, *args, **kwargs)[source]¶ Parallelizes calls to dmc_propagate by splitting configs
If npartitions does not evenly divide nconfigs, we need to reweight the results based on the number of configs per parallel task.
The final result should be equivalent to the non-parallelized case. The average weight \(w\) and the weighted average of observables \(\langle O \rangle\) are returned. Index \(i\) refers to walker index.
\[w = \sum_i w_i / n_{\rm config} \qquad\quad \langle O \rangle = \sum_i o_{i} w_i / \sum_i w_i\]Split over parallel tasks, we need to reweight by number of walkers. The average weight \(w_p\) and weighted average of observables \(\langle O\rangle_p\) are returned from each task.
\[w_p = \sum_j^{{\rm task}\, p} w_j / n_{{\rm config}, p} \qquad\quad \langle O \rangle_p = \frac{\sum_j^{{\rm task}\, p} o_{j} w_j }{ \sum_j^{{\rm task}\, p} w_j }\]The total weight and total average (defined above) are computed from the task weights \(w_p\) and task averages \(\langle O\rangle_p\) as
\[w = \sum_p w_p n_{{\rm config}, p} / n_{\rm config}, \qquad\quad \langle O \rangle = \frac{ \sum_p \langle O\rangle_p \sum_j^{{\rm task}\, p} w_j }{ \sum_i w_i}.\]We can rewrite the weights using the equations above
\[ \begin{align}\begin{aligned}\langle O \rangle &= \frac{ \sum_p \langle O\rangle_p w_p n_{{\rm config}, p} }{ w n_{\rm config} }\\&= \sum_p \langle O\rangle_p \frac{w_p n_{{\rm config}, p}}{\sum_p w_p n_{{\rm config}, p}}\end{aligned}\end{align} \]By reweighting the task weights as \(\overline{w}_p = w_p n_{{\rm config}, p}\), we can omit the reweighting factor \(\frac{n_{{\rm config}, p}}{n_{\rm config}}\) (that we use to collect parallel vmc). Instead, we use only the reweighting factor \(\overline{w}_p / \sum_p \overline{w}_p\)
\[\langle O \rangle = \sum_p \langle O\rangle_p \frac{\overline{w}_p }{\sum_p \overline{w}_p }\]
-
pyqmc.dmc.limdrift(g, tau, acyrus=0.25)[source]¶ Use Cyrus Umrigar’s algorithm to limit the drift near nodes.
Parameters: - g – a [nconf,ndim] vector
- tau – time step
- acyrus – the maximum magnitude
Returns: The vector with the cut off applied and multiplied by tau.
-
pyqmc.dmc.propose_tmoves(wf, configs, energy_accumulator, tstep, e)[source]¶ No side effect calculation of t-moves
- Returns:
- new proposed positions probability of acceptance sum of weights
-
pyqmc.dmc.rundmc(wf, configs, weights=None, tstep=0.01, nsteps=1000, branchtime=5, stepoffset=0, branchcut_start=10, verbose=False, accumulators=None, ekey=('energy', 'total'), feedback=1.0, hdf_file=None, continue_from=None, client=None, npartitions=None, **kwargs)[source]¶ Run DMC
Parameters: - wf – A Wave function-like class. recompute(), gradient(), and updateinternals() are used, as well as anything (such as laplacian() ) used by accumulators
- configs – (nconfig, nelec, 3) - initial coordinates to start calculation.
- weights – (nconfig,) - initial weights to start calculation, defaults to uniform.
- nsteps – number of DMC steps to take
- tstep – Time step for move proposals. Introduces time step error.
- branchtime – number of steps to take between branching
- accumulators – A dictionary of functor objects that take in (coords,wf) and return a dictionary of quantities to be averaged. np.mean(quantity,axis=0) should give the average over configurations. If none, a default energy accumulator will be used.
- ekey – tuple of strings; energy is needed for DMC weights. Access total energy by accumulators[ekey[0]](configs, wf)[ekey[1]
- verbose – Print out step information
- stepoffset – If continuing a run, what to start the step numbering at.
Returns: (df,coords,weights) df: A list of dictionaries nstep long that contains all results from the accumulators.
coords: The final coordinates from this calculation.
weights: The final weights from this calculation
Orthogonal optimization¶
-
pyqmc.optimize_ortho.collect_overlap_data(wfs, configs, pgrad)[source]¶ Collect the averages assuming that configs are distributed according to
\[\rho \propto \sum_i |\Psi_i|^2\]The keys ‘overlap’ and ‘overlap_gradient’ are
overlap :
\[\langle \Psi_f | \Psi_i \rangle = \left\langle \frac{\Psi_i^* \Psi_j}{\rho} \right\rangle_{\rho}\]overlap_gradient:
\[\partial_m \langle \Psi_f | \Psi_i \rangle = \left\langle \frac{\partial_{fm} \Psi_i^* \Psi_j}{\rho} \right\rangle_{\rho}\]
Given a configs sampled from the distribution
\[\rho = \sum_i \Psi_i^2\]Compute properties for replacing the last wave function with each of parameters
For energy
\[\langle E \rangle = \left\langle \frac{H\Psi}{\Psi} \frac{|\Psi|^2}{\rho} \right\rangle\]The ratio can be computed as
\[\frac{|\Psi|^2}{\rho} = \frac{1}{\sum_i e^{2(\alpha_i - \alpha}) }\]Where we write
\[\Psi = e^{i\theta}{e^\alpha}\]We also compute
\[\langle S_i \rangle = \left\langle \frac{\Psi_i^* \Psi}{\rho} \right\rangle\]And
\[\langle N_i \rangle = \left\langle \frac{|\Psi_i|^2}{\rho} \right\rangle\]
-
pyqmc.optimize_ortho.evaluate(return_data, warmup)[source]¶ For wave functions wfs and coordinate set coords, evaluate the overlap and energy of the last wave function.
Returns a dictionary with relevant information.
-
pyqmc.optimize_ortho.optimize_orthogonal(wfs, coords, pgrad, Starget=None, forcing=None, tstep=0.1, max_iterations=30, warmup=1, warmup_options=None, Ntarget=0.5, max_step=10.0, hdf_file=None, linemin=True, step_offset=0, Ntol=0.05, weight_boundaries=0.3, sample_options=None, correlated_options=None, client=None, npartitions=None)[source]¶ Minimize
\[f(p_f) = E_f + \sum_i \lambda_{i=0}^{f-1} |S_{fi} - S_{fi}^*|^2\]Where
\[N_i = \langle \Psi_i | \Psi_i \rangle\]\[S_{fi} = \frac{\langle \Psi_f | \Psi_i \rangle}{\sqrt{N_f N_i}}\]The *’d and lambda values are respectively targets and forcings. f is the final wave function in the wave function array. We only optimize the parameters of the final wave function, so all ‘p’ values here represent a parameter in the final wave function.
Important arguments
wfs: a list of wave function objects. The last one is optimized; the rest are kept fixed and used as orthogonalization references coords: A Coord set pgrad: A Pgradient object tstep: Maximum timestep for line minimization, or timestep when line minimization is off max_iterations: Number of optimization steps to take warmup_options: a dictionary of options for warm up vmc Starget: An array-like of length len(wfs)-1, which indicates the target overlap for each reference wave function. forcing: An array-like of length len(wfs)-1, which gives the penalty (lambda) for each reference wave function hdf_file: A string that gives the filename to save the optimization data in. Arguments for experts
Other arguments should not be changed unless you know what you’re doing.Details about the implementation
The derivatives are:
\[\partial_p N_f = 2 Re \langle \partial_p \Psi_f | \Psi_f \rangle\]\[\langle \partial_p \Psi_f | \Psi_i \rangle = \int{ \frac{ \Psi_i\partial_p \Psi_f^*}{\rho} \frac{\rho}{\int \rho} }\]\[\partial_p S_{fi} = \frac{\langle \partial_p \Psi_f | \Psi_i \rangle}{\sqrt{N_f N_i}} - \frac{\langle \Psi_f | \Psi_i \rangle}{2\sqrt{N_f N_i}} \frac{\partial_p N_f}{N_f}\]In this implementation, we set
\[\rho = \sum_i |\Psi_i|^2\]Note that in the definition of N there is an arbitrary normalization of rho. The targets are set relative to the normalization of the reference wave functions.
Some implementation notes regarding the normalization:
It’s important for the normalization of the wave functions to be similar; otherwise the weights of one dominate and only one of the wave functions gets sampled. Some notes:
- One could modify the relative weights in the definition of rho, but it’s more convenient to output wave functions that are normalized with respect to each other.
- Adding a penalty to the normalization turns out to be fairly unstable.
- Moves to reduce the overlap and energy tend to change the normalization a lot (keep in mind that both the determinant and Jastrow parts can have gauge degrees of freedom). This can lead into a tailspin effect pretty quickly.
In this implementation, we handle this using three techniques:
- Most importantly, the cost function derivative is orthogonalized to the derivative of the normalization.
- In the line minimization, if the normalization deviates too far from 0.5 relative to the reference wave function, we do not consider the move. The correlated sampling is unreliable in that situation anyway.
- The wave function is renormalized if its normalization deviates too far from 0.5 relative to the first wave function.
-
pyqmc.optimize_ortho.renormalize(wfs, N)[source]¶ Normalizes the last wave function, given a current value of the normalization. Assumes that we want N to be 0.5
\[ \begin{align}\begin{aligned}b^2/(a^2 + b^2) = N\\b^2 = N a^2 /(1-N)\\f^2 b^2 = 0.5 a^2/0.5 = a^2\\f^2 = a^2/b^2 = (1-N)/N\end{aligned}\end{align} \]
-
pyqmc.optimize_ortho.sample_overlap(wfs, configs, pgrad, nblocks=10, nsteps=10, tstep=0.5)[source]¶ Sample
\[\rho(R) = \sum_i |\Psi_i(R)|^2\]pgrad is expected to be a gradient generator. returns data as follows:
overlap :
\[\left\langle \frac{\Psi_i^* \Psi_j}{\rho} \right\rangle\]overlap_gradient:
\[\left\langle \frac{\partial_{im} \Psi_i^* \Psi_j}{\rho} \right\rangle\]Note that overlap_gradient is only saved for i = f, where f is the final wave function.
In addition, any key returned by pgrad will be saved for the final wave function.