diff --git a/README.md b/README.md
index 1bc3a8ffc..82e928244 100644
--- a/README.md
+++ b/README.md
@@ -54,6 +54,9 @@ import numpy as np
pos = np.load( 'tests/data/positions.npy')
vel = np.load( 'tests/data/velocities.npy')
box = np.load( 'tests/data/cells.npy')
+#atomic types: load an array with length equal to the number of atoms
+#in this example I set all the atoms to type 0 but the last 16,
+#wich are 8 of type 1 and 8 of type 2
types = np.zeros(pos.shape[1], dtype = np.int32)
types[-16:-8]=1
types[-8:]=2
@@ -62,13 +65,23 @@ print('first cell matrix is {}'.format(box[0]))
#create trajectory object
import pyanalisi
-analisi_traj = pyanalisi.Trajectory(pos, vel, types, box,True, False)
+analisi_traj = pyanalisi.Trajectory(
+ pos, #shape (ntimesteps, natoms, 3)
+ vel, #same shape as pos
+ types, #shape (natoms)
+ box, #shape (ntimesteps, 3,3) or (ntimesteps, 6)
+ True, # if True, box for a single frame is 3,3,
+ # if False is 6
+ False # don't apply pbc
+ )
#do the calculation that you want
+#see MSD section
msd=pyanalisi.MeanSquareDisplacement(analisi_traj,10,4000,4,True,False,False)
msd.reset(3000)
msd.calculate(0)
+#the calculation object can be used as an array
result=np.array(msd,copy=False)
@@ -81,17 +94,22 @@ result=np.array(msd,copy=False)
Heat transport coefficient calculation: correlation functions and gk integral for a multicomponent fluid example
```
import numpy as np
+#read first line, that is an header, and the column formatted data file
with open(filepath_tests + '/data/gk_integral.dat', 'r') as flog:
headers = flog.readline().split()
log = np.loadtxt(filepath_tests + '/data/gk_integral.dat', skiprows=1)
import pyanalisi
+#wrapper for the pyanalisi interface
traj = pyanalisi.ReadLog(log, headers)
+#see Green Kubo section
gk = pyanalisi.GreenKubo(analisi_log,'',
1,['c_flux[1]','c_vcm[1][1]'],
False, 2000, 2,False,0,4,False,1,100)
gk.reset(analisi_log.getNtimesteps()-2000)
gk.calculate(0)
+
+#the calculation object can be used as a python array
result = np.array(gk,copy=False)
```
@@ -105,19 +123,21 @@ If you use this program and you like it, spread the word around and give it cred
This is a framework for computing averages on molecular dynamics trajectories and block averages.
An mpi parallelization over the blocks is implemented in a very abstract and general way, so it is easy to implement a new calculation that can be parallelized with no effort using MPI. The code has two interfaces: a command line interface that is parallelized with MPI and threads, and a python interface that is parallelized on threads only. With such many parallelization levels more accurate averages can be performed, as described in details in the next sections.
+Every MPI processes uses the same amount of memory by loading the part of the trajectory that it is used to do the calculation in every block. For this reason it is suggested to use one MPI process per machine, and parallelize by using threads only inside the single machine.
+
Features:
- python interface (reads numpy array)
- command line interface (reads binary lammps-like files or time series in column format)
- multithreaded ( defaults to number of threads specifies by the shell variable OMP_NUM_THREADS )
- - command line interface has MPI too (for super-heavy calculations)
+ - command line interface has MPI too (for super-heavy multi-machine calculations)
- command line calculates variance of every quantity and every function (in python you can do it by yourselves with numpy.var )
- jupyter notebook example of the python interface
Calculations:
- g of r ( and time too!)
-- vibrational spectrum (this is nothing special)
+- vibrational spectrum with FFTW3
- histogram of number of neighbours
- mean square displacement
- green kubo integral of currents
@@ -132,7 +152,7 @@ Dependencies:
- C++17 compatible compiler (GCC 7+, clang 6+, intel 19.0.1+, [source](https://en.cppreference.com/w/cpp/compiler_support) )
- cmake
-- linux or macOS (mmap - maybe this dependency can be removed with some additional work)
+- linux or macOS (mmap - maybe this dependency can be removed with some additional work, and it is not used in the python interface)
- FFTW3 (included in the package)
- Eigen3 (included in the package)
- Boost (included in the package)
@@ -140,6 +160,14 @@ Dependencies:
- libxdrfile (for gromacs file conversion -- included in the package)
- python (optional)
+Documentation build:
+
+- r-rsvg
+- pandoc
+- rsvg-convert
+- texlive
+- readme2tex (pip)
+
## MPI build (why not?)
```
mkdir build
@@ -319,6 +347,13 @@ else:
```
+To access the atomic lammps ids and the atomic lamms types, you can use the following two functions:
+```
+analisi_traj.get_lammps_id()
+analisi_traj.get_lammps_type()
+```
+don't call them too often since every time a new numpy array of ints is created
+
### Creating a time series object
Before analyzing a time series, for example to calculate the integral of the autocorrelation function, it is necessary to create a time series object. This object will hold the data that a different function can analyze. It is created with the following:
@@ -439,11 +474,17 @@ where 
+
-where 
is the type of atom.
+where is the type of atom and 
is one of the spatial direction x,y,z.
The diffusivity can be computed as half of the zero value of D.
+The columns of the output are ordered like the following:
+
+
+
+where 
is the number of atomic species and 
the number of timesteps. Note that in the command line versione each column but the first is followed by its variance
+
### command line version
The options that you can use for this calculation are simply:
@@ -452,12 +493,15 @@ The options that you can use for this calculation are simply:
the code will generate a file with a number of line equal to the number of frequencies, and with `ntypes_atoms*2+1` columns where:
-- the first column rappresent the index of the frequncies
-- then there are the block average and variance of the spectrum for each atomic type
+
+ - the first column rappresent the index of the frequencies
+ - then there are the block average and variance of the spectrum for each atomic type
The code is trasparent ot units of measuares of the quantities. If a user wants the diffusivity in the correct units ( e.g. metal) must porcede in the following way:
-- the first column can be multiplied by `1/(nstep*dt)` to obtain the frequencies in multiples of Hz
-- the other columns can be multiplied by `dt/nstep`;
+
+ - the first column can be multiplied by `1/(nstep*dt)` to obtain the frequencies in multiples of Hz
+ - the other columns can be multiplied by `dt/nstep`;
+
where `nstep` is the total number of step of the block used to compute VDOS, `dt` is the time difference of two consecutive molecular dynamis steps.
### python interface
diff --git a/README.pdf b/README.pdf
index 7466f9d60..5e2c41fb8 100644
Binary files a/README.pdf and b/README.pdf differ
diff --git a/README_.md b/README_.md
index 8f8cc3139..66dbbd50b 100644
--- a/README_.md
+++ b/README_.md
@@ -160,6 +160,14 @@ Dependencies:
- libxdrfile (for gromacs file conversion -- included in the package)
- python (optional)
+Documentation build:
+
+- r-rsvg
+- pandoc
+- rsvg-convert
+- texlive
+- readme2tex (pip)
+
## MPI build (why not?)
```
mkdir build
@@ -339,6 +347,13 @@ else:
```
+To access the atomic lammps ids and the atomic lamms types, you can use the following two functions:
+```
+analisi_traj.get_lammps_id()
+analisi_traj.get_lammps_type()
+```
+don't call them too often since every time a new numpy array of ints is created
+
### Creating a time series object
Before analyzing a time series, for example to calculate the integral of the autocorrelation function, it is necessary to create a time series object. This object will hold the data that a different function can analyze. It is created with the following:
@@ -530,12 +545,24 @@ where $\langle \cdot \rangle$ is an average operator, and we do an additional av
The implemented equation is:
$$
-D_{\alpha}(\omega) = \frac{1}{3N_{\alpha}} \sum_{n}^{N_{\alpha}} \int_{- \infty}^{ \infty} \langle \mathbf{v}_{n}(0)\cdot \mathbf{v}_{n} (t) \rangle e^{i \omega t} dt
+D^{k}_{\alpha}(\omega) = \frac{1}{3N_{\alpha}} \sum_{n}^{N_{\alpha}} \int_{- \infty}^{ \infty} \langle v^{k}_{n}(0) v^{k}_{n} (t) \rangle e^{i \omega t} dt
$$
-where $\alpha$ is the type of atom.
+where $\alpha$ is the type of atom and $k$ is one of the spatial direction x,y,z.
The diffusivity can be computed as half of the zero value of D.
+The columns of the output are ordered like the following:
+
+$$
+\begin{aligned}
+&\{0&,&D_{0}^{x}(\omega_0),D_{0}^{y}(\omega_0),D_{0}^{z}(\omega_0),& \dots ,&D_{M}^{x}(\omega_0),D_{M}^{y}(\omega_0),D_{M}^{z}(\omega_0)&\}\\
+&\{\;\vdots&,&\dots, & \dots,&\vdots &\}\\
+&\{N&,&D_{N}^{x}(\omega_N),D_{N}^{y}(\omega_N),D_{N}^{z}(\omega_N),&\dots ,&D_{M}^{x}(\omega_N),D_{M}^{y}(\omega_N),D_{M}^{z}(\omega_N)&\}\\
+\end{aligned}
+$$
+
+where $M$ is the number of atomic species and $N$ the number of timesteps. Note that in the command line versione each column but the first is followed by its variance
+
### command line version
The options that you can use for this calculation are simply:
@@ -544,12 +571,15 @@ The options that you can use for this calculation are simply:
the code will generate a file with a number of line equal to the number of frequencies, and with `ntypes_atoms*2+1` columns where:
-- the first column rappresent the index of the frequncies
-- then there are the block average and variance of the spectrum for each atomic type
+
+ - the first column rappresent the index of the frequencies
+ - then there are the block average and variance of the spectrum for each atomic type
The code is trasparent ot units of measuares of the quantities. If a user wants the diffusivity in the correct units ( e.g. metal) must porcede in the following way:
-- the first column can be multiplied by `1/(nstep*dt)` to obtain the frequencies in multiples of Hz
-- the other columns can be multiplied by `dt/nstep`;
+
+ - the first column can be multiplied by `1/(nstep*dt)` to obtain the frequencies in multiples of Hz
+ - the other columns can be multiplied by `dt/nstep`;
+
where `nstep` is the total number of step of the block used to compute VDOS, `dt` is the time difference of two consecutive molecular dynamis steps.
### python interface
diff --git a/build_pdf.sh b/build_pdf.sh
index ef8853a34..3fa8cb676 100755
--- a/build_pdf.sh
+++ b/build_pdf.sh
@@ -1,2 +1,2 @@
-pandoc README_.md -V geometry:a4paper -o README.pdf
+pandoc README_.md -V geometry:a4paper -o README.pdf --toc
python -m readme2tex --output README.md README_.md --nocdn
diff --git a/config.h.in b/config.h.in
index cb66801ab..9edd1e3a1 100644
--- a/config.h.in
+++ b/config.h.in
@@ -53,7 +53,7 @@ const static char * _info_msg=
#ifdef XDR_FILE
"\nWith gromacs XDR file conversion support"
#endif
-"\nv0.3.2"
+"\nv0.3.3"
;
diff --git a/lib/include/traiettoria.h b/lib/include/traiettoria.h
index cb8cd7e8b..6a976c800 100644
--- a/lib/include/traiettoria.h
+++ b/lib/include/traiettoria.h
@@ -109,8 +109,8 @@ class Traiettoria : public TraiettoriaBase
Traiettoria::Errori imposta_inizio_accesso(const size_t & timesteps);
int64_t get_timestep_lammps(size_t timestep);
void index_all();
- // void set_calculate_center_of_mass(bool);
- // bool get_calculate_center_of_mass();
+ int * get_lammps_id();
+ int *get_lammps_type();
private:
std::mapid_map;
size_t * timesteps; // puntatori (offset rispetto all'inizio) all'inizio di ogni timesteps
diff --git a/lib/src/traiettoria.cpp b/lib/src/traiettoria.cpp
index f762d6c90..518a5fbf0 100644
--- a/lib/src/traiettoria.cpp
+++ b/lib/src/traiettoria.cpp
@@ -148,7 +148,7 @@ void Traiettoria::init_buffer_tipi() {
Atomo::read_id_tipo(pezzi[ichunk].atomi+iatomo*sizeof(double)*NDOUBLE_ATOMO,id__,type__);
int id=round(id__);
if (id<0) {
- throw std::runtime_error("Error: found a negativa atomic id in the trajectory");
+ throw std::runtime_error("Error: found a negative atomic id in the trajectory");
}
if (id_map.find(id)==id_map.end()){
id_map[id]=id_;
@@ -178,17 +178,44 @@ void Traiettoria::init_buffer_tipi() {
for (int ichunk=0;ichunk
#include "traiettoria.h"
#include "spettrovibrazionale.h"
#include "readlog.h"
@@ -242,6 +242,35 @@ PYBIND11_MODULE(pyanalisi,m) {
g.get_shape(), /* Buffer dimensions */
g.get_stride());
})
+ .def("get_lammps_id", [](Traiettoria & t) {
+ int * foo=t.get_lammps_id();
+ pybind11::capsule free_when_done(foo, [](void *f) {
+ int *foo = reinterpret_cast(f);
+ std::cerr << "freeing memory @ " << f << "\n";
+ delete[] foo;
+ });
+ size_t natoms=t.get_natoms();
+ return pybind11::array_t(
+ {natoms}, //shape
+ {sizeof(int)},
+ foo,
+ free_when_done
+ );})
+ .def("get_lammps_type", [](Traiettoria & t) {
+ int * foo=t.get_lammps_type();
+ pybind11::capsule free_when_done(foo, [](void *f) {
+ int *foo = reinterpret_cast(f);
+ std::cerr << "freeing memory @ " << f << "\n";
+ delete[] foo;
+ });
+ size_t natoms=t.get_natoms();
+ return pybind11::array_t(
+ {natoms}, //shape
+ {sizeof(int)},
+ foo,
+ free_when_done
+ );})
+
;
using RL = ReadLog;
diff --git a/svgs/10ccf1728c4af71363a03a251dd4f874.svg b/svgs/10ccf1728c4af71363a03a251dd4f874.svg
new file mode 100644
index 000000000..2cf8b8e60
--- /dev/null
+++ b/svgs/10ccf1728c4af71363a03a251dd4f874.svg
@@ -0,0 +1,146 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/svgs/2dffde5e18cde45e1794c82325eef7bb.svg b/svgs/2dffde5e18cde45e1794c82325eef7bb.svg
new file mode 100644
index 000000000..363bb8f57
--- /dev/null
+++ b/svgs/2dffde5e18cde45e1794c82325eef7bb.svg
@@ -0,0 +1,73 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/svgs/39796bb2cab9694a0ae799160b754d39.svg b/svgs/39796bb2cab9694a0ae799160b754d39.svg
new file mode 100644
index 000000000..f3c8ffa3d
--- /dev/null
+++ b/svgs/39796bb2cab9694a0ae799160b754d39.svg
@@ -0,0 +1,154 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/svgs/63bb9849783d01d91403bc9a5fea12a2.svg b/svgs/63bb9849783d01d91403bc9a5fea12a2.svg
new file mode 100644
index 000000000..853335889
--- /dev/null
+++ b/svgs/63bb9849783d01d91403bc9a5fea12a2.svg
@@ -0,0 +1,9 @@
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/svgs/6be1aa73ae5f695ef047e31f66db0a2a.svg b/svgs/6be1aa73ae5f695ef047e31f66db0a2a.svg
new file mode 100644
index 000000000..f510e3f51
--- /dev/null
+++ b/svgs/6be1aa73ae5f695ef047e31f66db0a2a.svg
@@ -0,0 +1,146 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/svgs/7ded6b2994b505ccc139d8bee7e8cef1.svg b/svgs/7ded6b2994b505ccc139d8bee7e8cef1.svg
new file mode 100644
index 000000000..11c806cda
--- /dev/null
+++ b/svgs/7ded6b2994b505ccc139d8bee7e8cef1.svg
@@ -0,0 +1,150 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/svgs/a016ef07bbe1fd04e814fe2391916e01.svg b/svgs/a016ef07bbe1fd04e814fe2391916e01.svg
new file mode 100644
index 000000000..91ce5fd59
--- /dev/null
+++ b/svgs/a016ef07bbe1fd04e814fe2391916e01.svg
@@ -0,0 +1,73 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file