mirror of https://github.com/QMCPACK/qmcpack.git
add better test of latdev
sample two free particles in Gaussian orbitals
This commit is contained in:
Paul Young
parent
055d680c94
commit
9089302115
|
|
@ -45,6 +45,18 @@ if (add_test)
|
|||
)
|
||||
endif()
|
||||
|
||||
set(latdev_free_python_reqs numpy;h5py)
|
||||
CHECK_PYTHON_REQS(latdev_python_reqs estimator-latdev-free add_test)
|
||||
|
||||
if (add_test)
|
||||
SIMPLE_RUN_AND_CHECK(estimator-latdev-free
|
||||
"${CMAKE_SOURCE_DIR}/tests/estimator/latdev/free"
|
||||
two.xml
|
||||
1 16
|
||||
flatdev.py
|
||||
)
|
||||
endif()
|
||||
|
||||
set(sofk_python_reqs numpy;pandas;h5py)
|
||||
CHECK_PYTHON_REQS(sofk_python_reqs estimator-sofk add_test)
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,74 @@
|
|||
#!/usr/bin/env python
|
||||
import h5py
|
||||
import numpy as np
|
||||
|
||||
def corr(trace):
|
||||
""" calculate the autocorrelation of a trace of scalar data
|
||||
|
||||
correlation time is defined as the integral of the auto-correlation
|
||||
function from t=0 to when the function first reaches 0.
|
||||
|
||||
Args:
|
||||
trace (list): should be a 1D iterable array of floating point numbers
|
||||
Return:
|
||||
float: correlation_time, the autocorrelation time of this trace of scalars
|
||||
"""
|
||||
mu = np.mean(trace)
|
||||
stddev = np.std(trace, ddof=1)
|
||||
if np.isclose(stddev, 0): # easy case
|
||||
return np.inf # infinite correlation for constant trace
|
||||
correlation_time = 0.
|
||||
for k in range(1, len(trace)):
|
||||
# calculate auto_correlation
|
||||
auto_correlation = 0.0
|
||||
num = len(trace)-k
|
||||
auto_correlation = np.dot(trace[:num]-mu, trace[k:]-mu)
|
||||
auto_correlation *= 1.0/(num*stddev**2)
|
||||
if auto_correlation > 0:
|
||||
correlation_time += auto_correlation
|
||||
else:
|
||||
break
|
||||
correlation_time = 1.0 + 2.0*correlation_time
|
||||
return correlation_time
|
||||
|
||||
def get_latdev_me(fh5, nequil):
|
||||
# get raw data
|
||||
fp = h5py.File(fh5, 'r')
|
||||
data = fp['latdev/value'][()]
|
||||
fp.close()
|
||||
|
||||
nb, ncol = data.shape
|
||||
npt = nb-nequil
|
||||
# calculate mean
|
||||
ym = data[nequil:, :].mean(axis=0)
|
||||
# calculate standard error
|
||||
yc = np.array([corr(data[nequil:, icol]) for icol in range(ncol)])
|
||||
neff = npt/yc
|
||||
ye = data[nequil:, :].std(ddof=1, axis=0)/neff**0.5
|
||||
return ym, ye
|
||||
|
||||
def check_latdev(w1, w2, ym, ye, ndim=3):
|
||||
expect = np.array([w1**2/2.]*ndim + [w2**2/2.]*ndim)
|
||||
zscore = abs(ym-expect)/ye
|
||||
fail = np.any(zscore > 3.)
|
||||
success = ~fail
|
||||
return success, zscore
|
||||
|
||||
if __name__ == '__main__':
|
||||
# calculate <x^2> <y^2> <z^2>
|
||||
nequil = 6
|
||||
fh5 = 'mt.s000.stat.h5'
|
||||
ym, ye = get_latdev_me(fh5, nequil)
|
||||
|
||||
# check answers
|
||||
w1 = 0.3
|
||||
w2 = 0.5
|
||||
success, zscore = check_latdev(w1, w2, ym, ye)
|
||||
if not success:
|
||||
print(zscore)
|
||||
|
||||
if success:
|
||||
exit(0)
|
||||
else:
|
||||
exit(1)
|
||||
# end __main__
|
||||
|
|
@ -0,0 +1,61 @@
|
|||
<simulation>
|
||||
<project id="mt" series="0"/>
|
||||
|
||||
<qmcsystem>
|
||||
<simulationcell>
|
||||
<parameter name="lattice" units="bohr">
|
||||
10 0 0
|
||||
0 10 0
|
||||
0 0 10
|
||||
</parameter>
|
||||
<parameter name="bconds">
|
||||
p p p
|
||||
</parameter>
|
||||
</simulationcell>
|
||||
|
||||
<particleset name="e">
|
||||
<group name="m5" size="2">
|
||||
<parameter name="mass"> 5.0 </parameter>
|
||||
<attrib name="position" datatype="posArray" condition="0">
|
||||
1.0001 1.0001 1.0002
|
||||
6.0001 6.0001 6.0002
|
||||
</attrib>
|
||||
</group>
|
||||
</particleset>
|
||||
<particleset name="trap_center">
|
||||
<group name="center" size="2">
|
||||
<attrib name="position" datatype="posArray" condition="0">
|
||||
1 1 1
|
||||
6 6 6
|
||||
</attrib>
|
||||
</group>
|
||||
</particleset>
|
||||
<particleset name="measure_center">
|
||||
<group name="origin" size="2">
|
||||
<attrib name="position" datatype="posArray" condition="0">
|
||||
1 1 1
|
||||
6 6 6
|
||||
</attrib>
|
||||
</group>
|
||||
</particleset>
|
||||
|
||||
<wavefunction target="e" id="psi0">
|
||||
<ionwf name="ionwf" source="trap_center" width="0.3 0.5"/>
|
||||
</wavefunction>
|
||||
|
||||
<hamiltonian name="h0" type="generic" target="e">
|
||||
<estimator name="skinetic" type="specieskinetic"/>
|
||||
<estimator type="latticedeviation" name="latdev" hdf5="yes" per_xyz="yes" target="e" tgroup="m5" source="measure_center" sgroup="origin"/>
|
||||
<!--
|
||||
<estimator name="dens" type="density" delta="0.04 0.04 0.04" x_min="0" x_max="10" y_min="0" ymax="10" z_min="0" z_max="10"/>
|
||||
-->
|
||||
</hamiltonian>
|
||||
</qmcsystem>
|
||||
|
||||
<qmc method="vmc" move="whatever">
|
||||
<parameter name="blocks"> 16 </parameter>
|
||||
<parameter name="steps"> 1024 </parameter>
|
||||
<parameter name="substeps"> 8 </parameter>
|
||||
<parameter name="timestep"> 0.3 </parameter>
|
||||
</qmc>
|
||||
</simulation>
|
||||
Loading…
Reference in New Issue