Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Return correct shapes from Correlator #3848

Merged
merged 11 commits into from
Aug 11, 2020
4 changes: 2 additions & 2 deletions doc/tutorials/01-lennard_jones/01-lennard_jones.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -549,9 +549,9 @@
"source": [
"fig3 = plt.figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')\n",
"fig3.set_tight_layout(False)\n",
"lag_time = msd[:, 0]\n",
"lag_time = msd_corr.lag_times()\n",
"for i in range(0, N_PART, 30):\n",
" msd_particle_i = msd[:, 2+i*3] + msd[:, 3+i*3] + msd[:, 4+i*3]\n",
" msd_particle_i = msd[:, i, 0] + msd[:, i, 1] + msd[:, i, 2]\n",
" plt.plot(lag_time, msd_particle_i,\n",
" 'o-', linewidth=2, label=\"particle id =\" + str(i))\n",
"plt.xlabel(r'Lag time $\\tau$ [$\\delta t$]', fontsize=20)\n",
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,7 @@
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import logging\n",
"import sys\n",
"\n",
Expand Down Expand Up @@ -138,7 +139,9 @@
" for i in range(LOOPS):\n",
" system.integrator.run(STEPS)\n",
" correlator.finalize()\n",
" msd_results.append(correlator.result())\n",
" msd_results.append(np.column_stack((correlator.lag_times(),\n",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

" correlator.sample_sizes(),\n",
" correlator.result().reshape([-1, 3]))))\n",
"\n",
"logging.info(\"Sampling finished.\")"
]
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -147,8 +147,8 @@
"correlator.finalize()\n",
"corrdata = correlator.result()\n",
"\n",
"tau = corrdata[:, 0]\n",
"msd = corrdata[:, 2] + corrdata[:, 3] + corrdata[:, 4]\n",
"tau = correlator.lag_times()\n",
"msd = corrdata[:, 0] + corrdata[:, 1] + corrdata[:, 2]\n",
"\n",
"rh = system.analysis.calc_rh(chain_start=0, number_of_chains=1, chain_length=N_MONOMERS)[0]"
]
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -119,8 +119,8 @@
corrdata = correlator.result()

rh_results[index] = np.average(rhs)
tau = corrdata[1:, 0]
msd = corrdata[1:, 2] + corrdata[1:, 3] + corrdata[1:, 4]
tau = correlator.lag_times()[1:]
msd = corrdata[1:, 0] + corrdata[1:, 1] + corrdata[1:, 2]
np.save('msd_{}'.format(N), np.c_[tau, msd])

np.save('rh.npy', rh_results)
5 changes: 4 additions & 1 deletion doc/tutorials/06-active_matter/06-active_matter.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -385,7 +385,10 @@
"# Finalize the correlator and write to disk\n",
"system.auto_update_accumulators.remove(msd)\n",
"msd.finalize()\n",
"numpy.savetxt(\"output.dat\", msd.result())\n",
"numpy.savetxt(\"output.dat\",\n",
" numpy.column_stack((msd.lag_times(),\n",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

" msd.sample_sizes(),\n",
" msd.result().reshape([-1, 3]))))\n",
"```\n",
"\n",
"Here, the observable `pos_id` is set to the particle positions of the\n",
Expand Down
15 changes: 12 additions & 3 deletions doc/tutorials/06-active_matter/solutions/enhanced_diffusion.py
Original file line number Diff line number Diff line change
Expand Up @@ -119,12 +119,21 @@
# Finalize the accumulators and write to disk
system.auto_update_accumulators.remove(msd)
msd.finalize()
np.savetxt("{}/msd_{}_{}.dat".format(outdir, vel, run), msd.result())
np.savetxt("{}/msd_{}_{}.dat".format(outdir, vel, run),
np.column_stack((msd.lag_times(),
msd.sample_sizes(),
msd.result().reshape([-1, 3]))))

system.auto_update_accumulators.remove(vacf)
vacf.finalize()
np.savetxt("{}/vacf_{}_{}.dat".format(outdir, vel, run), vacf.result())
np.savetxt("{}/vacf_{}_{}.dat".format(outdir, vel, run),
np.column_stack((vacf.lag_times(),
vacf.sample_sizes(),
vacf.result().reshape([-1, 1]))))

system.auto_update_accumulators.remove(avacf)
avacf.finalize()
np.savetxt("{}/avacf_{}_{}.dat".format(outdir, vel, run), avacf.result())
np.savetxt("{}/avacf_{}_{}.dat".format(outdir, vel, run),
np.column_stack((avacf.lag_times(),
avacf.sample_sizes(),
avacf.result().reshape([-1, 1]))))
15 changes: 9 additions & 6 deletions samples/diffusion_coefficient.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,27 +55,30 @@
c_pos.finalize()
c_vel.finalize()

np.savetxt("msd.dat", c_pos.result())
np.savetxt("vacf.dat", c_vel.result())
msd = np.column_stack((c_pos.lag_times(),
c_pos.sample_sizes(),
c_pos.result().reshape([-1, 3])))
vacf = np.column_stack((c_vel.lag_times(),
c_vel.sample_sizes(),
c_vel.result().reshape([-1, 1])))
np.savetxt("msd.dat", msd)
np.savetxt("vacf.dat", vacf)

# Integral of vacf via Green-Kubo
# D= 1/3 int_0^infty <v(t_0)v(t_0+t)> dt

vacf = c_vel.result()
# Integrate with trapezoidal rule
I = np.trapz(vacf[:, 2], vacf[:, 0])
ratio = 1. / 3. * I / (kT / gamma)
print("Ratio of measured and expected diffusion coefficients from Green-Kubo:",
ratio)

# Check MSD
msd = c_pos.result()


def expected_msd(x):
return 2. * kT / gamma * x


# Check MSD
print("Ratio of expected and measured msd")
print("#time ratio_x ratio_y ratio_z")
for i in range(4, msd.shape[0], 4):
Expand Down
8 changes: 6 additions & 2 deletions samples/observables_correlators.py
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,10 @@

# Finalize the correlation calculation and write the results to a file
c.finalize()
np.savetxt("res.dat", c.result())
np.savetxt("res.dat", np.column_stack((c.lag_times(),
c.sample_sizes(),
c.result().reshape([-1, 3]))))
fcs.finalize()
np.savetxt("fcs.dat", fcs.result())
np.savetxt("fcs.dat", np.column_stack((fcs.lag_times(),
fcs.sample_sizes(),
fcs.result().reshape([-1, 1]))))
5 changes: 5 additions & 0 deletions src/core/accumulators/AccumulatorBase.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,9 @@
#ifndef CORE_ACCUMULATORS_ACCUMULATORBASE
#define CORE_ACCUMULATORS_ACCUMULATORBASE

#include <cstddef>
#include <vector>

namespace Accumulators {

class AccumulatorBase {
Expand All @@ -28,6 +31,8 @@ class AccumulatorBase {
virtual ~AccumulatorBase() = default;

virtual void update() = 0;
/** Dimensions needed to reshape the flat array returned by the accumulator */
virtual std::vector<size_t> shape() const = 0;

private:
// Number of timesteps between automatic updates.
Expand Down
Loading