Initial Conditions in DMC
Typically, a blown up equillibrium geometry is enough for the DMC to begin sampling the ground state. However, to aid in equilibration, we can use more advanced methods.
Permutations of Like Atoms
If you want to start with a ground state geometry that has a random distribution of like atoms, you can do this with the
InitialConditioner object. For example, if you have a methane molecule, you may want to
start with initial conditions so that there is not a bias in picking where “Hydrogen 1” is versus “Hydrogen 2”. This
technique will duplicate the geometry to the specified ensemble size, and for each duplication it will swap like atoms
randomly (if you start with ["C","H1","H2","H3","H4"], this code will randomly permute all Hs to get things like
["C","H1","H2","H4","H3"], ["C","H4","H3","H2","H1"], and ["C","H3","H1","H2","H4"]).
You can also use this for selective permutations, in the protonated water dimer (H5O2+), you can just permute the two hydrogen atoms on either side of the water. You will pass in an list of lists:
import pyvibdmc as pv
from pyvibdmc import Constants
from pyvibdmc import potential_manager as pm
ch4 = np.array([[0.000000000000000, 0.000000000000000, 0.000000000000000],
[0.1318851447521099, 2.088940054609643, 0.000000000000000],
[1.786540362044548, -1.386051328559878, 0.000000000000000],
[2.233806981137821, 0.3567096955165336, 0.000000000000000],
[-0.8247121421923925, -0.6295306113384560, -1.775332267901544])
ch4 = Constants.convert(ch4, "angstroms", to_AU=True) # To Bohr from angstroms
atms = ["C", "H", "H", "H", "H"]
initializer = pv.InitialConditioner(coord=ch4,
atoms=atms,
num_walkers=50000,
technique='permute_atoms',
technique_kwargs={'like_atoms': [[0],[1, 2, 3, 4]],
'ensemble': None})
new_coords = initializer.run()
h5o2 = np.array([[ 2.25486805, 0. , 0. ],
[-2.25486805, 0. , 0. ],
[ 0. , 0. , 0.12351035],
[-3.01319472, 1.48918379, -0.746598 ],
[-3.20158756, -0.45448108, 1.49813637],
[ 3.01319472, -1.48918379, -0.746598 ],
[ 3.20158756, 0.45448108, 1.49813637]])
atms = ["O", "O", "H", "H", "H", "H", "H"]
# 3 and 4 are the hydrogen atoms on one side, the 5 and 6 are on the other. They will not permute with each other,
# only the two pairs by themselves.
initializer = InitialConditioner(coord=h5o2,
atoms=atms,
num_walkers=50000,
technique='permute_atoms',
technique_kwargs={'like_atoms': [[0],[1],[2],[3,4],[5,6]],
'ensemble': None})
new_coords = initializer.run()
The ensemble option is for if you want to pass in more than just a duplicated minimum energy structure.
If you have an ensemble and you want to permute its like atoms, you can feed in a num_walkers, num_atoms, 3 array.
Sampling from the Harmonic Ground State
To sample from the separable, 3N-6 dimensional Gaussian distribution that is the harmonic ground state, you must run a (Cartesian) normal mode analysis calculation first. If you have run any electronic structure calculations before (Psi4, NWChem, Gaussian, etc.), this is equivalent to running a frequency calculation, just using fitted potential energy surface rather than an ab initio one.:
import pyvibdmc as pv
from pyvibdmc import Constants
from pyvibdmc import potential_manager as pm
dxx = 1.e-3 # bohr
water_geom = np.array([[0.9578400, 0.0000000, 0.0000000],
[-0.2399535, 0.9272970, 0.0000000],
[0.0000000, 0.0000000, 0.0000000]])
# Everything is in Atomic Units going into generating the Hessian.
pot_dir = path/to/Partridge_Schwenke_H2O/')
py_file = 'h2o_potential.py'
pot_func = 'water_pot'
partridge_schwenke = pm.Potential(potential_function=pot_func,
potential_directory=pot_dir,
python_file=py_file,
num_cores=1)
geom = Constants.convert(water_geom, "angstroms", to_AU=True) # To Bohr from angstroms
atms = ["H", "H", "O"]
harm_h2o = pv.HarmonicAnalysis(eq_geom=geom,
atoms=atms,
potential=partridge_schwenke,
dx=dxx)
freqs, normal_modes = harm_h2o.run()
# Turns of scientific notation
np.set_printoptions(suppress=True)
print(f"Freqs (cm-1): {freqs}")
The HarmonicAnalysis object can also take in the arguments points_diag=__ and points_off_diag=__. This
refers to the number of finite difference points used in the generation of the Hessian matrix. These numbers default to
5 and 3 respectively, meaning that the on-diagonal second derivatives are generated using a 5-point finite difference,
and the off-diagonal mixed derivatives use a 3 point finite difference in both dimensions. Currently, this code only
supports using 3 or 5 point finite difference for either argument.
The 3N frequencies and normal modes that are returned from the harmonic analysis include the 6 near-zero modes from
the translational and rotational degrees of freedom (this code does not support linear molecules).
From there, you will pass these frequencies and normal modes to the InitialConditioner, which will generate the
desired ensemble of walkers that we will feed into the DMC.:
# Do initial conditions based on freqs and normal modes
initializer = pv.InitialConditioner(coord=water_geom,
atoms=atms,
num_walkers=50000,
technique='harmonic_sampling',
technique_kwargs={'freqs': freqs,
'normal_modes': normal_modes,
'scaling_factor': 1.5,
'ensemble': None})
new_coords = initializer.run()
The technique_kwargs you see above are all necessary to pass in. The scaling_factor broadens the 3N-6 dimensional
Gaussian distribution by a uniform factor in all dimensions. In the case above, it is equivalent to saying the
harmonic frequencies are all divided by 1.5, which will give you a broader distribution that the
walkers will sample from. This technique is described in more detail
in this paper.
The ensemble argument is present so that you can pass in a whole ensemble that will be displaced along those normal
modes randomly if desired. If left as None, then it will simply duplicate the minimum energy geometry you supplied,
and you can ignore the next code block in the tutorial.
If you feed in a num_walkers, num_atoms, 3 array, you can combine this with the permute_atoms method above;
start by randomly displacing along the harmonic ground state, then permuting like atoms:
from pyvibdmc import InitialConditioner, HarmonicAnalysis
from pyvibdmc import Constants
from pyvibdmc import potential_manager as pm
ch4 = np.array([[0.000000000000000, 0.000000000000000, 0.000000000000000],
[0.1318851447521099, 2.088940054609643, 0.000000000000000],
[1.786540362044548, -1.386051328559878, 0.000000000000000],
[2.233806981137821, 0.3567096955165336, 0.000000000000000],
[-0.8247121421923925, -0.6295306113384560, -1.775332267901544])
ch4 = Constants.convert(ch4, "angstroms", to_AU=True) # To Bohr from angstroms
atms = ["C", "H", "H", "H", "H"]
# Run harmonic analysis
freqs, normal_modes = HarmonicAnalysis(...)
# Then, push the freqs, normal modes, and ensemble to the InitialConditioner
initializer = InitialConditioner(coord=ch4,
atoms=atms,
num_walkers=50000,
technique='harmonic_sampling',
technique_kwargs={'freqs': freqs,
'normal_modes': normal_modes,
'scaling_factor': 1.5},
'ensemble': None)
harm_coords = initializer.run()
# Finally, then permute like atoms for each walker that are now spread along the harmonic ground state.
initializer = InitialConditioner(coord=ch4,
atoms=atms,
num_walkers=50000,
technique='permute_atoms',
technique_kwargs={'like_atoms': [[0],[1, 2, 3, 4]],
'ensemble': harm_coords})
new_coords = initializer.run()
Now, the harmonically-sampled-then-permuted new_coords are passed to the DMC_Sim object and used during the DMC run:
import pyvibdmc as pv
myDMC = pv.DMC_Sim(sim_name=f"conditioner_{sim_num}",
output_folder="initial_conditions_tutorial",
weighting='discrete', #or 'continuous'. 'continuous' keeps the ensemble size constant.
num_walkers=50000, #number of geometries exploring the potential surface
num_timesteps=10000, #how long the simulation will go. (num_timesteps * delta_t atomic units of time)
equil_steps=500, #how long before we start collecting wave functions
chkpt_every=500, #checkpoint the simulation every "chkpt_every" time steps
wfn_every=1000, #collect a wave function every "wfn_every" time steps
desc_wt_steps=100, #number of time steps you allow for descendant weighting per wave function
atoms=['H','H','O'],
delta_t=10, #the size of the time step in atomic units
potential=water_pot,
start_structures=new_coords,
masses=None #can put in artificial masses, otherwise it auto-pulls values from the atoms string
)
Note, you should only run the harmonic calculation THEN permute, not the other way around. This is because this code produces the eigenvectors of the Hessian that only correspond to the atom ordering of the non-permuted molecular system. You can, of course, do either individually or use neither technique before a DMC run.