diff --git a/.github/workflows/check-version.yml b/.github/workflows/check-version.yml index dda4592..df2f23c 100644 --- a/.github/workflows/check-version.yml +++ b/.github/workflows/check-version.yml @@ -23,11 +23,13 @@ jobs: shell: bash run: | set -euo pipefail + if ! git tag --list | grep -q .; then echo "No tags found; skipping check." echo "skip=true" >> "$GITHUB_OUTPUT" exit 0 fi + tag="$(git describe --tags --abbrev=0)" echo "Latest tag: $tag" echo "tag=$tag" >> "$GITHUB_OUTPUT" diff --git a/Makefile b/Makefile index 0d87982..e98a759 100644 --- a/Makefile +++ b/Makefile @@ -122,9 +122,6 @@ XDR_OBJS := $(patsubst $(XDR_SRC_DIR)/%.c,$(BUILD_DIR)/xdr/%.o,$(XDR_CSRCS)) # Some distros require -lm for math symbols used by xdrfile LDLIBS += -lm -# Ensure dir exists -$(shell mkdir -p $(BUILD_DIR)/xdr) - # --- Source discovery (flat, only src/*.cpp) --------------------------------- SRCS := $(wildcard $(SRC_DIR)/*.cpp) @@ -138,9 +135,6 @@ SRCS := $(filter-out $(EXCLUDE_SRCS),$(SRCS)) OBJS := $(patsubst $(SRC_DIR)/%.cpp,$(BUILD_DIR)/%.o,$(SRCS)) DEPS := $(OBJS:.o=.d) -# Ensure build dirs exist -$(shell mkdir -p $(BUILD_DIR)) - # --- Targets ----------------------------------------------------------------- .PHONY: all clean distclean rebuild format lint help all: $(APP) @@ -204,11 +198,13 @@ $(APP): $(OBJS) $(XDR_OBJS) $(VORO_OBJS) | qvectors voro # Compile C++: strictly src/.cpp → build/.o $(BUILD_DIR)/%.o: $(SRC_DIR)/%.cpp | $(VERSION_H) @echo " CXX $<" + @mkdir -p $(dir $@) $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -o $@ # Compile C from xdrfile: third_party/xdrfile-1.1b/src/*.c → build/xdr/*.o $(BUILD_DIR)/xdr/%.o: $(XDR_SRC_DIR)/%.c @echo " CC $<" + @mkdir -p $(dir $@) $(CC) $(XDR_CFLAGS) -c $< -o $@ ## Compile voro++: third_party/voro++-0.4.6/src/*.cc → build/voro/*.o (isolated includes) diff --git a/Manual/isf.md b/Manual/isf.md new file mode 100644 index 0000000..512d362 --- /dev/null +++ b/Manual/isf.md @@ -0,0 +1,31 @@ +

isf

+

function

+ +Calculates the full intermediate scattering function of the system via the Fourier transform of the density. + +

Syntax

+ +``` +structure_factor + + (if symmetry = asymmetric only) +``` + +`` is either `auto` or `manual`. The `manual` functionality is not currently documented; employ auto. + +`` is either `symmetric` or `asymmetric`. +If `symmetric`, the analysis calculates the structure factor between the set of particles specified in target and itself. +If `asymmetric`, then a second target must be specified, and the partial structure factor describing correlations of `` with only `` is calculated. +`` and `` are the indices of the limits on time spacings to be calculated. +Options for ```` are ``xyz``, ``xy``, ``xz``, ``yz``, ``x``, ``y``, and ``z``. +This chooses which dimensions in k-space to include in the calculation of the intermediate scattering function. +``xyz`` computes the full radial three dimensional isf, ``xy``, ``yz``, and ``xz`` calculate two-dimensional in-plane radial isf's, and ``x``, ``y``, and ``z`` compute one-dimensional isf's. +```` determines the longest distance which will be decomposed into inverse space. +If a distance of 0 is given, the full box size is used. +Any deviation from '0' will in general produce incorrect results for the structure factor, especially at low k. + +`` determines what times to loop over. +If `timescheme` is -1, loop over all times. +If `timescheme` is zero or positive, only use one time per block, with the value setting the time index offset from the beginning of the block. +In most cases this should be set to either 0 or -1, with -1 giving improved statistical strength at the cost of much longer compute times. +The level of improvement in statistical strength will depend on the timescale for structural decorrelation in comparison to the length of a trajectory block. diff --git a/Manual/mean_displacement.md b/Manual/mean_displacement.md index c511fab..74c590d 100644 --- a/Manual/mean_displacement.md +++ b/Manual/mean_displacement.md @@ -1,5 +1,7 @@

mean_displacement

+[![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md) + Calculates the mean (not mean square) displacement vector of all particles as a function of time. When used with trajectory_bin_lists, particularly useful for extracting gradients in displacement of the system. Syntax is as follows. _mean\_displacement \_ diff --git a/Manual/msd.md b/Manual/msd.md index 100428e..f11ff52 100644 --- a/Manual/msd.md +++ b/Manual/msd.md @@ -1,4 +1,7 @@

msd

+ +[![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md) +

Function

Calculates mean square displacement, as follows. @@ -14,3 +17,5 @@ _msd \_ _\_ + + diff --git a/Manual/msd_2d.md b/Manual/msd_2d.md index 8d23fac..5a3a994 100644 --- a/Manual/msd_2d.md +++ b/Manual/msd_2d.md @@ -1,5 +1,7 @@

msd_2d

+[![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md) +

Function

Calculates mean square displacement in two dimensions diff --git a/Manual/multithreading.md b/Manual/multithreading.md new file mode 100644 index 0000000..26ab830 --- /dev/null +++ b/Manual/multithreading.md @@ -0,0 +1,61 @@ +

AMDAT Multithreading

+ +## Overview +High-performance computing workflows can be computationally expensive, and molecular dynamics analysis is no exception. To accelerate selected analyses, AMDAT supports multithreading with OpenMP. The analyses that currently support multithreading are listed below. + +--- +| Analysis | Description | Maximum Speedup Tested (Compared to Serial) | +|----------|-------------|------------------------| +| [`md`](mean_displacement.md) | Calculate Mean Displacement | 7x | +| [`msd`](msd.md) | Calculate Mean Square Displacement | 10x | +| [`msd_2d`](msd_2d.md) | Calculate Mean Square Displacement for two-dimensional analysis | 9x | +| [`isfs`](isfs.md) | Calculate Incoherent Scattering Function | 20x | +| [`vhs`](vhs.md) | Calculate self part of Van Hove Function | 12x | +| [`vhd`](vhd.md) | Calculate distinct part of Van Hove Function | 20x | +| `baf` | Calculate Bond Autocorrelation Function | 15x | +| [`ngp`](ngp.md) | Calculate Non Gaussian Parameter | 30x | +| [`composition`](composition.md) | Calculate Composition | 7x | +| [`rdf`](rdf.md) | Calculate Radial Distribution Function | 21x | +| [`rnf`](rnf.md) | Calculate Radial Count | 23x | + + +## Inconsistency in performance between analysis + +AMDAT is currently implemented with a highly object-oriented design and an array-of-structures (AoS) memory layout. Hence, data are organized in adjacent memory blocks, making them more prone to Cache Coherency issue. In parallel system, this conflict results in false-sharing, meaning different threads try to have writing access to one common cache line, despite aiming at independent memory location. To fix performance bottleneck, AMDAT uses padding to separate each memory location further from each other, reducing Cache Coherency conflicts. This approach, however, trades spatial locality from serial implementation for mitigated memory conflict in parallel tasks, making perfect speed-up realistically unachievable. Therefore, speed-ups are inconsistent between analyses due to the differences in how we store the data. + +## System Requirements + +To use multithreading effectively in AMDAT, make sure the following requirements are met: + +- Your system has a multi-core CPU +- The OpenMP runtime is available in your environment (This is already handled by Conda) + +## Analysis Requirements + +Parallelism is introduced to reduce the runtime of time-based loops. +However, to avoid oversubscription and ensure efficient resource utilization, +the optimal number of threads is allocated can be given by: + +$$N_{\text{threads}} = \min\left(N_{\text{iterations}},\;N_{\text{CPUs}} \right).$$ + +where: +- $$\ N_{\text{CPUs}} \$$ is the maximum number of hardware threads (logical CPUs) available on the node. +- $$\ N_{\text{iterations}} \$$ is the number of independent time-based iterations. In AMDAT's input files, it is highlighted as shown below: + +> system +> xyz_log +> ./testfiles/traj.xyz ./testfiles/log.lammps +> exponential **`70`** 77 1.05 0 0 .01 +> polymer 50 +> +> ... + +## Usage + +Running AMDAT is very simple with conda build, because OpenMP is already preinstalled. The flag `-n` will allocate the desired number of parallel regions for your analysis. + +### Example with msd + +```bash +./AMDAT -n 8 -i ./testfiles/msd.in + diff --git a/Manual/ngp.md b/Manual/ngp.md index ccef5f9..9010109 100644 --- a/Manual/ngp.md +++ b/Manual/ngp.md @@ -1,5 +1,7 @@

ngp

+[![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md) +

Function

Calculates non-Gaussian parameter of the mean squared displacement. diff --git a/Manual/overview.md b/Manual/overview.md index 1408aba..3b08b4a 100644 --- a/Manual/overview.md +++ b/Manual/overview.md @@ -6,7 +6,7 @@ AMDAT can read multiple file formats produced by the LAMMPS molecular dynamics s A central feature of AMDAT is its ability to read in and analyze blocked, exponentially-spaced trajectories. This is essential for analysis for dynamics over many orders of magnitude in time. The handling of this is currently documented in the [input_files](input_file.md) documentation. Both the documentation and the underlying functionality for this are expected to be expanded in the near future. -AMDAT includes no built-in tools for visualization of either trajectories or properties of trajectories. All results are output in a text-based form and must be viewed in another software (exe python, julia, excel, matlab,...). AMDAT does not include the ability to aggregate statistics acquired based on multiple simulations; this requires use of post-analysis software, such as python, julia, matlab, igor, excel, or similar. AMDAT presently does not include multithreading, although this capability is planned for introduction in a future update. +AMDAT includes no built-in tools for visualization of either trajectories or properties of trajectories. All results are output in a text-based form and must be viewed in another software (exe python, julia, excel, matlab,...). AMDAT does not include the ability to aggregate statistics acquired based on multiple simulations; this requires use of post-analysis software, such as python, julia, matlab, igor, excel, or similar. AMDAT includes CPU multithreading for select analysis function; more information is provided below or at this link: [multithreading](multithreading.md).

Making AMDAT

@@ -38,7 +38,7 @@ AMDAT is run from the command line, specifying the filename of a user [input fil `-n <# OF PROCESSORS>` -Presently, no analysis methods allow for multithreading, so this option does nothing. Eventually, this will permit AMDAT to run some analyses on `<# OF PROCESSORS>` processing cores. Defaults to 1 if not given.. +This sets AMDAT to run certain analyses (those for which multithreading is enabled) on `<# OF PROCESSORS>` processing cores. Defaults to 1 if not given. Analysis tools with multithreading capabilities are identified by this flag in their documentation page: [![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md). Additional information and a full list of analysis tools with multithreading capabilities can be found at [multithreading](multithreading.md). `-c ` diff --git a/Manual/radial_count.md b/Manual/radial_count.md index 8513142..dbcfd8b 100644 --- a/Manual/radial_count.md +++ b/Manual/radial_count.md @@ -1,5 +1,8 @@

radial_count

+[![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md) + +--- Computes a non-normalized radial distribution function - essentially the mean number of particles in shells a distance $r$ to $r+\Delta r$ from a central particle. Syntax is as follows. diff --git a/Manual/rdf.md b/Manual/rdf.md index 921d7f6..e4a54f1 100644 --- a/Manual/rdf.md +++ b/Manual/rdf.md @@ -1,4 +1,7 @@

rdf

+ +[![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md) +

Function

diff --git a/Manual/trajectory_lists.md b/Manual/trajectory_lists.md index 9615c99..b58c2fb 100644 --- a/Manual/trajectory_lists.md +++ b/Manual/trajectory_lists.md @@ -50,7 +50,7 @@ A large number of analysis tools target `trajectory_list` objects. | [`displacement_dist`](displacement_dist.md) | Computes the displacement of each particle at some time separation and reports a distribution | Tab demarcated data file. | | [`displacement_map`](displacement_map.md) | Generates maps of particle displacement magnitude and writes to pdb files with value beta specifying displacement.| PDB file | | [`incremental_mean_displacement`](incremental_mean_displacement.md) | Calculates mean displacement of particles as a function of time. | Tab demarcated data file. | -| [`isf_list`](isf_list.md) | Calculate full intermediate scattering function | Tab demarkated data file | +| [`isf`](isf.md) | Calculate full intermediate scattering function | Tab demarkated data file | | [`isfs`](isfs.md) | Calculates the self-part of the intermediate scattering function. | Tab demarkated data file | | [`mean_closest_distance`](mean_closest_distance.md) | Calculates the average distance between particles in two distinct `trajectory_list`s. | Tab demarkated data file | | [`mean_displacement`](mean_displacement.md) | Calculates the mean (not mean square) displacement vector of all particles as a function of time. | Tab demarcated data file | @@ -81,7 +81,6 @@ A large number of analysis tools target `trajectory_list` objects. | `find_fast_fixedthreshold` | Finds fast particles based on displacements higher than some value. Likely to become deprecated in favor of valuelist approach. | new `trajectory_list` | | [`find_between`](find_between.md) | Finds particles with displacements in some value range. Likely to become deprecated in favor of valuelist approach. | new `trajectory_list` | | `gyration_radius` | Deprecated; to be removed. | - | -| `isf` | Deprecated; to be removed. | - | | `radial_debye_waller` | Computes dwf as a function of distance from (0,0,0). Planned for deprecation.| | `vector_autocorrelation_function` | Computes autocorrelation function for a specified set of vectors. Essentially replaced by [`raf`](raf.md) [`multibody`](multibodies.md) functionality and planned for deprecation. | Tab-demarcated data file. | diff --git a/Manual/vhd.md b/Manual/vhd.md index 6e1ed68..8061b03 100644 --- a/Manual/vhd.md +++ b/Manual/vhd.md @@ -1,5 +1,7 @@

vhd

+[![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md) +

Function

Calculates distinct part of the Van Hove function. The distinct van Hove is by construction asymmetric, computed between two distinct sets of atoms $a$ and $b$ ($a$ and $b$ may be the same). When asymmetric, $a$ is the central particle in the $G(r,\Delta t)$ calculation and $b$ is the set of particles for which densities around the central particles are computed. diff --git a/Manual/vhs.md b/Manual/vhs.md index c908289..e6e837b 100644 --- a/Manual/vhs.md +++ b/Manual/vhs.md @@ -1,5 +1,7 @@

vhs

+[![multi](https://img.shields.io/badge/Multithreaded-blue)](multithreading.md) + Calculates self-part of the Van Hove function. ``` diff --git a/README.md b/README.md index 2210361..7873f55 100644 --- a/README.md +++ b/README.md @@ -15,12 +15,12 @@ It provides a wide variety of high-performance analyses integral to molecular mo ## Why AMDAT (at a glance) * Fast, in-memory engine – load once, analyze many time delays without re-reading files. -* Blocked exponential time spacing – efficient long-timescale dynamics across orders of magnitude. -* Modular data abstractions – trajectory / neighbor / multibody / value lists compose into rich workflows. -* Validated analyses – - * static and time-resolved structure factors, - * radial distribution functions, - * mean-square displacements, +* Blocked [exponential time spacing](/Manual/time_scheme.md) – efficient long-timescale dynamics across orders of magnitude. +* Modular data abstractions – [trajectory](Manual/trajectory_lists.md) / [neighbor](Manual/neighborlist.md) / [multibody](Manual/multibodies.md) / [value](Manual/value_list.md) lists compose into rich workflows. +* Validated analyses such as - + * [static](Manual/structure_factor.md) and [time-resolved](Manual/isf_list.md) structure factors, + * [radial distribution functions](Manual/rdf.md), + * [mean-square displacements](Manual/msd.md), * neighbor correlations, * clustering, and more. * Plain-text outputs – easy post-processing in Python/Matlab/Excel/VMD/OVITO. @@ -162,6 +162,9 @@ If you use AMDAT, please cite: > Kawak, Drayer, and Simmons, "AMDAT: An Open-Source Molecular Dynamics Analysis Toolkit for Supercooled Liquids, Glass-Forming Materials, and Complex Fluids". arXiv, 2026. DOI: [10.48550/arXiv.2602.05865](https://doi.org/10.48550/arXiv.2602.05865) +and +> https://arxiv.org/abs/2602.05865 + *Also see [CITATION.cff](CITATION.cff).* ## Authors & Maintainers diff --git a/VERSION b/VERSION index 7dea76e..9084fa2 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.0.1 +1.1.0 diff --git a/src/analysis.h b/src/analysis.h index 53311cf..c1b5117 100644 --- a/src/analysis.h +++ b/src/analysis.h @@ -88,6 +88,7 @@ class Analysis virtual void listkernel(Trajectory*, int, int, int){cout<<"Error: Trajectory list targets not fully implemented for this analysis method.\n";}; //added by Michael? virtual void listkernel2(Trajectory*, Trajectory*, int, int, int){cout<<"Error: Trajectory list targets not fully implemented for this analysis method.\n";}; //listkernel for use only when two (nested) trajectory loops are needed. virtual void listkernel(Trajectory*){cout<<"Error: Trajectory list targets not fully implemented for this analysis method.\n";}; + virtual void listkernel(Trajectory*,int){cout<<"Error: Trajectory list targets not fully implemented for this analysis method.\n";}; /*System loop methods*/ diff --git a/src/analysis_onetime.cpp b/src/analysis_onetime.cpp index 2aa3a29..da3e3b8 100644 --- a/src/analysis_onetime.cpp +++ b/src/analysis_onetime.cpp @@ -29,6 +29,7 @@ void Analysis_Onetime::analyze(Trajectory_List * t_list) if(time_scheme==-1) { + #pragma omp parallel for schedule(static) if(this->isThreadSafe()) for (timeii=0; timeiishow_n_timesteps();timeii++) { timekernel(timeii); @@ -40,6 +41,7 @@ void Analysis_Onetime::analyze(Trajectory_List * t_list) } else { + #pragma omp parallel for schedule(static) if(this->isThreadSafe()) for (timeii=time_scheme; timeiishow_n_exponentials();timeii+=system->show_n_exponential_steps()) { timekernel(timeii); @@ -62,6 +64,7 @@ void Analysis_Onetime::analyze(Trajectory_List * t_list,Trajectory_List* t_list2 if(time_scheme==-1) { + #pragma omp parallel for schedule(static) if(this->isThreadSafe()) for (timeii=0; timeiishow_n_timesteps();timeii++) { timekernel2(timeii); @@ -73,6 +76,7 @@ void Analysis_Onetime::analyze(Trajectory_List * t_list,Trajectory_List* t_list2 } else { + #pragma omp parallel for schedule(static) if(this->isThreadSafe()) for (timeii=time_scheme; timeiishow_n_exponentials();timeii+=system->show_n_exponential_steps()) { timekernel2(timeii); diff --git a/src/bond_autocorrelation_function.cpp b/src/bond_autocorrelation_function.cpp index 9d87b0a..f346fa0 100644 --- a/src/bond_autocorrelation_function.cpp +++ b/src/bond_autocorrelation_function.cpp @@ -22,8 +22,8 @@ Bond_Autocorrelation_Function::Bond_Autocorrelation_Function() n_times = 0; //allocate memory for mean square displacement data - baf = new float [n_times]; - weighting = new int [n_times]; + baf = new float * [n_times]; + weighting = new int * [n_times]; atomcount = 0; dimensions.set(1,1,1); @@ -43,15 +43,20 @@ Bond_Autocorrelation_Function::Bond_Autocorrelation_Function(const Bond_Autocorr n_times = copy.n_times; atomcount = copy.atomcount; - baf = new float [n_times]; - weighting = new int [n_times]; + baf = new float * [n_times]; + weighting = new int * [n_times]; + + for(int i = 0; i < n_times; i++){ + baf[i] = new float[16]; + weighting[i] = new int[16]; + } timetable = system->displacement_times(); for(timeii=0;timeiishow_n_timegaps(); //allocate memory for mean square displacement data - baf = new float [n_times]; - weighting = new int [n_times]; + baf = new float * [n_times]; + weighting = new int * [n_times]; timetable = system->displacement_times(); for(timeii=0;timeiishow_n_timegaps(); //allocate memory for mean square displacement data - baf = new float [n_times]; - weighting = new int [n_times]; + baf = new float * [n_times]; + weighting = new int * [n_times]; + + for(int i = 0; i < n_times; i++){ + baf[i] = new float[16]; + weighting[i] = new int[16]; + } timetable = system->displacement_times(); for(timeii=0;timeiishow_n_timegaps(); //allocate memory for mean square displacement data - baf = new float [n_times]; - weighting = new int [n_times]; + baf = new float * [n_times]; + weighting = new int * [n_times]; + + for(int i = 0; i < n_times; i++){ + baf[i] = new float[16]; + weighting[i] = new int[16]; + } timetable = system->displacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiishow_n_multibodies(thisii); + weighting[timegapii][0]+=multibody_list->show_n_multibodies(thisii); multibody_list->listloop(this,timegapii, thisii, nextii); } @@ -339,7 +392,7 @@ void Bond_Autocorrelation_Function::listkernel(Multibody* current_multibody, int float dotproduct = (this->*vprep)((*current_multibody)(1)->show_unwrapped(thisii)-(*current_multibody)(0)->show_unwrapped(thisii))&(this->*vprep)((*current_multibody)(1)->show_unwrapped(nextii)-(*current_multibody)(0)->show_unwrapped(nextii)); //compute dot product between unit vectors at initial and later times - baf[timegapii]+=(this->*legendre_p)(dotproduct);//increment baf by chosen legendre polynomial of dot product above + baf[timegapii][0]+=(this->*legendre_p)(dotproduct);//increment baf by chosen legendre polynomial of dot product above } @@ -349,7 +402,7 @@ void Bond_Autocorrelation_Function::postprocess_list() for(int timeii=0;timeiishow_n_trajectories(timeii); - current_time = timeii; - current_density[current_time] = current_total_atoms/volume; + int current_total_atoms=trajectory_list->show_n_trajectories(timeii); + int current_time = timeii; + current_density[current_time][0] = current_total_atoms/volume; - trajectory_list->listloop(this,current_time); + trajectory_list->listloop(this,0,current_time,0); + #pragma omp atomic total_atoms += current_total_atoms; - average_density+=current_density[current_time]; + #pragma omp atomic + average_density+=current_density[current_time][0]; for (int typeii=0;typeiishow_type()-1; - #pragma omp atomic current_comp[current_time][traj_type]++; } diff --git a/src/composition.h b/src/composition.h index 48eacee..3bd7900 100644 --- a/src/composition.h +++ b/src/composition.h @@ -23,12 +23,12 @@ class Composition: public Analysis_Onetime int n_atomtypes; int n_molecules; int n_times; - int current_time; - int current_total_atoms; + //int current_time; + //int current_total_atoms; int total_atoms; float volume; - float * current_density; + float ** current_density; float average_density; float * time_average_comp; float** current_comp; @@ -46,12 +46,13 @@ class Composition: public Analysis_Onetime Analysis_Type what_are_you(){Analysis_Type type = composition; return type;}; - void listkernel(Trajectory *); + //void listkernel(Trajectory *); + void listkernel(Trajectory *,int,int,int); void timekernel(int); void postprocess_list(); void write(string); void write(ofstream&); - bool isThreadSafe(){return false;}; + bool isThreadSafe(){return true;}; diff --git a/src/composition_timedependent.cpp b/src/composition_timedependent.cpp index 78a119a..40f1b76 100644 --- a/src/composition_timedependent.cpp +++ b/src/composition_timedependent.cpp @@ -43,8 +43,8 @@ void Composition_TimeDependent::write(string filename) for (int timeii=0; timeiishow_time(timeii); - output << "\t" << current_density[timeii]; - output << "\t" << current_density[timeii]*volume; + output << "\t" << current_density[timeii][0]; + output << "\t" << current_density[timeii][0]*volume; for(int typeii=0;typeiishow_time(timeii); - output << "\t" << current_density[timeii]; + output << "\t" << current_density[timeii][0]; for (int typeii=0;typeii (ndf,runline,filename); // pass run_analysis template the analysis type 'Neighbor_Decorrelation_Function' - finish = time(NULL); cout << "\nCalculated neighbor_decorrelation_function in " << finish-start<<" seconds."< (ngpar,runline,filename); // pass run_analysis template the analysis type //ngpar.write(filename); cout << "\n Peak time index of non-Gaussian parameter is " << ngpar.max() << "."; - } @@ -5132,6 +5130,7 @@ void Control::orientational_correlation() start = time(NULL); oc.analyze(multibodylist); // pass run_analysis template the analysis type 'Mean_Square_Displacement' finish = time(NULL); + cerr << duration; cout << "\nCalculated bond autocorrelation function in " << finish-start<<" seconds."< //#include #include +#include #include #include #include @@ -87,9 +88,18 @@ class Control Space_Time_Correlation_Function vht; Gaussian_Comparison * gaussian_comparison; //array of Gaussian comparison objects int n_gaussian_comparisons; - time_t start; //timer start - time_t finish; //timer stop - + //Timer + time_t start, finish; + + //High-resolution timer (For algorithm benchmarking) + using clock = std::chrono::steady_clock; + using timepoint = clock::time_point; + timepoint high_start; //timer start + timepoint high_finish; //timer stop + float duration; + float time_count(timepoint start, timepoint finish){ + return std::chrono::duration(finish - start).count(); + } /*Arrays to store analysis results with a name given by the user, for later recall and use in other analysis techniques*/ //To be added @@ -291,6 +301,7 @@ Analysis_type Control::run_analysis(Analysis_type analyzer, string setline, stri Tokenize runtokenize(setline); n_setargs = runtokenize.count(); //n_setargs = tokenize(setline, setargs); + //high_start = clock::now(); if ( n_setargs==0 ) { cout << "Error: No atom set command found."; @@ -385,7 +396,9 @@ Analysis_type Control::run_analysis(Analysis_type analyzer, string setline, stri //analyzer.analyze(setline); //analyzer.write(filename); } - + //high_finish = clock::now(); + //duration = time_count(high_start, high_finish); + //cerr << duration; return analyzer; } diff --git a/src/incoherent_scattering_function.cpp b/src/incoherent_scattering_function.cpp index 42f4a1f..90d2dc5 100644 --- a/src/incoherent_scattering_function.cpp +++ b/src/incoherent_scattering_function.cpp @@ -227,9 +227,9 @@ void Incoherent_Scattering_Function::analyze(Trajectory_List * t_list) void Incoherent_Scattering_Function::list_displacementkernel(int timegapii,int thisii, int nextii) { - currenttime=thisii; - nexttime=nextii; - timegap=timegapii; + int currenttime=thisii; + int nexttime=nextii; + int timegap=timegapii; // //cout<show_n_trajectories(currenttime)<<"\t"<show_n_trajectories(currenttime))/float(timegap_weighting[timegap]); // //n_atoms[timegap]+=float(trajectory_list->show_n_trajectories(currenttime)); @@ -293,7 +293,6 @@ void Incoherent_Scattering_Function::listkernel(Trajectory* current_trajectory, coordinate2 = current_trajectory->show_unwrapped(nextii); tempcorrelation += double(cos(vectorlist[wavevectorii]&(coordinate2-coordinate1))) / double(vectorcount); } - #pragma omp atomic correlation[timegapii][wavenumberii] += float(tempcorrelation); } } diff --git a/src/incoherent_scattering_function.h b/src/incoherent_scattering_function.h index a32cdf1..cf124a2 100644 --- a/src/incoherent_scattering_function.h +++ b/src/incoherent_scattering_function.h @@ -34,7 +34,7 @@ class Incoherent_Scattering_Function:public Correlation_2D void bin_hook(Trajectory_List * t_list, int timegapii, int thisii, int nextii); void postprocess_bins(); - //bool isThreadSafe(){return true;}; + bool isThreadSafe(){return true;}; }; } diff --git a/src/intermediate_scattering_function.h b/src/intermediate_scattering_function.h index 132592d..4231aac 100644 --- a/src/intermediate_scattering_function.h +++ b/src/intermediate_scattering_function.h @@ -32,7 +32,7 @@ class Intermediate_Scattering_Function:public Correlation_2D //void displacementkernel(int timegapii, int thisii, int nextii, int species_index, int molecule_index, int atom_type, int atom_index); void displacementkernel(int timegapii, int thisii, int nextii, Trajectory * traj); - + bool isThreadSafe(){return true;} }; } diff --git a/src/mean_displacement.cpp b/src/mean_displacement.cpp index cc7f2b8..a499a79 100644 --- a/src/mean_displacement.cpp +++ b/src/mean_displacement.cpp @@ -12,14 +12,15 @@ #include "static_trajectory_list.h" using namespace std; +const int PAD = 16; Mean_Displacement::Mean_Displacement() { n_times = 0; //allocate memory for mean square displacement data - md = new Coordinate [n_times]; - weighting = new int [n_times]; + md = new Coordinate * [n_times]; + weighting = new int * [n_times]; atomcount = 0; } @@ -35,15 +36,21 @@ Mean_Displacement::Mean_Displacement(const Mean_Displacement & copy) n_times = copy.n_times; atomcount = copy.atomcount; - md = new Coordinate [n_times]; - weighting = new int [n_times]; + md = new Coordinate * [n_times]; + weighting = new int * [n_times]; + + for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiishow_n_timegaps(); //allocate memory for mean square displacement data - md = new Coordinate [n_times]; - weighting = new int [n_times]; + md = new Coordinate * [n_times]; + weighting = new int * [n_times]; + + for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiitimegap_weighting(); + int timeii; + weighting_temp = system->timegap_weighting(); + for(timeii=0;timeiishow_n_trajectories(currenttime); + weighting[timegapii][0]+=trajectory_list->show_n_trajectories(currenttime); //weighting[timegapii]+=(trajectory_list[0]).show_n_trajectories(currenttime); - (trajectory_list[0]).listloop(this,currenttime); + (trajectory_list[0]).listloop(this,currenttimegap,currenttime,nexttime); } -void Mean_Displacement::listkernel(Trajectory* current_trajectory) +void Mean_Displacement::listkernel(Trajectory* current_trajectory, int currenttimegap, int currenttime, int nexttime) { - md[currenttimegap]+=current_trajectory->show_unwrapped(nexttime)-current_trajectory->show_unwrapped(currenttime); + Coordinate c_next = current_trajectory->show_unwrapped(nexttime); + Coordinate c_this = current_trajectory->show_unwrapped(currenttime); + md[currenttimegap][0]+=c_next-c_this; } @@ -173,7 +211,7 @@ void Mean_Displacement::postprocess_list() for(int timeii=0;timeii using namespace std; +const int PAD = 16; Mean_Square_Displacement::Mean_Square_Displacement() { n_times = 0; //allocate memory for mean square displacement data - msd = new float [n_times]; - weighting = new float [n_times]; + msd = new float * [n_times]; + weighting = new float * [n_times]; atomcount = 0; } @@ -36,15 +37,20 @@ Mean_Square_Displacement::Mean_Square_Displacement(const Mean_Square_Displacemen n_times = copy.n_times; atomcount = copy.atomcount; - msd = new float [n_times]; - weighting = new float [n_times]; + msd = new float * [n_times]; + weighting = new float * [n_times]; + + for(int i = 0; i < n_times; i++){ + msd[i] = new float[PAD]; + weighting[i] = new float[PAD]; + } timetable = system->displacement_times(); for(timeii=0;timeiishow_n_timegaps(); //allocate memory for mean square displacement data - msd = new float [n_times]; - weighting = new float [n_times]; + msd = new float * [n_times]; + weighting = new float * [n_times]; + + for(int i = 0; i < n_times; i++){ + msd[i] = new float[PAD]; + weighting[i] = new float[PAD]; + } + timetable = system->displacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiishow_n_timegaps(); //allocate memory for mean square displacement data + for(int i = 0; i < n_times; i++){ + delete [] msd[i]; + delete [] weighting[i]; + } delete [] msd; delete [] weighting; - msd = new float [n_times]; - weighting = new float [n_times]; + msd = new float * [n_times]; + weighting = new float * [n_times]; + + for(int i = 0; i < n_times; i++){ + msd[i] = new float[PAD]; + weighting[i] = new float[PAD]; + } timetable = system->displacement_times(); for(timeii=0;timeiishow_n_trajectories(currenttime); // //weighting[timegapii]+=(trajectory_list[0]).show_n_trajectories(currenttime); // (trajectory_list[0]).listloop(this,currenttime); - #pragma omp atomic - weighting[timegapii]+=trajectory_list->show_n_trajectories(thisii); + weighting[timegapii][0]+=trajectory_list->show_n_trajectories(thisii); (trajectory_list[0]).listloop(this,timegapii, thisii, nextii); } @@ -163,8 +189,7 @@ void Mean_Square_Displacement::list_displacementkernel(int timegapii,int thisii, void Mean_Square_Displacement::listkernel(Trajectory* current_trajectory, int timegapii,int thisii, int nextii) { - #pragma omp atomic - msd[timegapii]+=current_trajectory->distance(thisii,nextii)*current_trajectory->distance(thisii,nextii); + msd[timegapii][0]+=current_trajectory->distance(thisii,nextii)*current_trajectory->distance(thisii,nextii); } @@ -174,7 +199,7 @@ void Mean_Square_Displacement::postprocess_list() for(int timeii=0;timeiidisplacement_times(); for (int timeii=0;timeiidisplacement_times(); for (int timeii=0;timeiishow_n_timegaps(); //allocate memory for mean square displacement data - msd = new float [n_times]; - weighting =system->timegap_weighting(); + msd = new float * [n_times]; + weighting = new int * [n_times]; + + for(int i = 0; i < n_times; i++){ + msd[i] = new float[PAD]; + weighting[i] = new int[PAD]; + } + + weighting_temp = system->timegap_weighting(); timetable = system->displacement_times(); + for(timeii=0;timeiidisplacement_list(this); postprocess_list(); @@ -134,16 +160,16 @@ void Mean_Square_Displacement_2D::analyze(Trajectory_List * t_list) void Mean_Square_Displacement_2D::list_displacementkernel(int timegapii,int thisii, int nextii) { - currenttime=thisii; - nexttime=nextii; - currenttimegap=timegapii; - weighting[timegapii]+=trajectory_list[0].show_n_trajectories(currenttime); - (trajectory_list[0]).listloop(this,currenttime); + int currenttimegap=timegapii; + int currenttime=thisii; + int nexttime=nextii; + weighting[timegapii][0]+=trajectory_list[0].show_n_trajectories(currenttime); + (trajectory_list[0]).listloop(this,currenttimegap,currenttime,nexttime); } -void Mean_Square_Displacement_2D::listkernel(Trajectory* current_trajectory) +void Mean_Square_Displacement_2D::listkernel(Trajectory* current_trajectory, int currenttimegap,int currenttime, int nexttime) { - msd[currenttimegap]+=(current_trajectory->*distancefun)(currenttime,nexttime)*(current_trajectory->*distancefun)(currenttime,nexttime); + msd[currenttimegap][0]+=(current_trajectory->*distancefun)(currenttime,nexttime)*(current_trajectory->*distancefun)(currenttime,nexttime); } void Mean_Square_Displacement_2D::postprocess_list() @@ -151,7 +177,7 @@ void Mean_Square_Displacement_2D::postprocess_list() int timeii; for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiishow_n_timegaps(); //allocate memory for mean square displacement data - ndf = new float [n_times]; - weighting = new float [n_times]; + ndf = new float * [n_times]; + weighting = new float * [n_times]; timetable = system->displacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiidisplacement_times(); for(timeii=0;timeiishow_n_trajectories(thisii); + weighting[timegapii][0]+=trajectory_list->show_n_trajectories(thisii); (trajectory_list[0]).listloop(this,timegapii, thisii, nextii); } @@ -169,7 +177,7 @@ void Neighbor_Decorrelation_Function::listkernel(Trajectory* current_trajectory, trajID=current_trajectory->show_trajectory_ID(); - ndf[timegapii]+=n_list->n_persistent_neighbors(trajID,thisii,nextii); + ndf[timegapii][0]+=n_list->n_persistent_neighbors(trajID,thisii,nextii); } @@ -179,7 +187,7 @@ void Neighbor_Decorrelation_Function::postprocess_list() for(int timeii=0;timeii displacement_times(); for(int timeii=0;timeii displacement_times(); for(int timeii=0;timeii displacement_times(); for(timeii=0;timeiishow_n_trajectories(thisii); + int currenttime=thisii; + int nexttime=nextii; + int currenttimegap=timegapii; + weighting[timegapii][0]+=trajectory_list->show_n_trajectories(thisii); (trajectory_list[0]).listloop(this,timegapii,thisii,nextii); } @@ -124,8 +142,8 @@ void Non_Gaussian_Parameter::listkernel(Trajectory* current_trajectory) void Non_Gaussian_Parameter::listkernel(Trajectory* current_trajectory, int timegapii,int thisii, int nextii) { - #pragma omp atomic - ngp[timegapii]+=pow(current_trajectory->distance(thisii,nextii),4); + double dist = current_trajectory->distance(thisii,nextii); + ngp[timegapii][0]+=dist*dist*dist*dist; } @@ -135,8 +153,9 @@ void Non_Gaussian_Parameter::postprocess_list() for(timeii=0;timeiishow(timeii)),2.0)); - ngp[timeii] -= 1.0; + double timestepmsd = msd->show(timeii); + ngp[timeii][0] *= (3.0/(float(weighting[timeii][0])))/(5*timestepmsd*timestepmsd); + ngp[timeii][0] -= 1.0; } } @@ -153,7 +172,7 @@ void Non_Gaussian_Parameter::write(string filename)const output << "Non-Gaussian parameter data created by AMDAT v." << amdat::build::SEMVER << "\n"; for(timeii=0;timeiimaxvalue) + if(ngp[timeii][0]>maxvalue) { - maxvalue = ngp[timeii]; + maxvalue = ngp[timeii][0]; maxtime = timeii; } } diff --git a/src/non_gaussian_parameter.h b/src/non_gaussian_parameter.h index 0dee15f..3ff4815 100644 --- a/src/non_gaussian_parameter.h +++ b/src/non_gaussian_parameter.h @@ -16,15 +16,15 @@ namespace std{ class Non_Gaussian_Parameter: public Analysis { int n_times; - float * ngp; - long int * weighting; + float ** ngp; + long int ** weighting; Mean_Square_Displacement const * msd; float * timetable; int atomcount; float * n_atoms; //calculation variables - int currenttime, nexttime, currenttimegap; + //int currenttime, nexttime, currenttimegap; public: Non_Gaussian_Parameter(); @@ -47,7 +47,7 @@ class Non_Gaussian_Parameter: public Analysis void bin_hook(Trajectory_List * t_list, int timegapii, int thisii, int nextii); void postprocess_bins(); - //bool isThreadSafe(){return true;}; + bool isThreadSafe(){return true;}; }; } diff --git a/src/orientational_correlation.cpp b/src/orientational_correlation.cpp index c0cd4af..f6a24d7 100644 --- a/src/orientational_correlation.cpp +++ b/src/orientational_correlation.cpp @@ -21,8 +21,8 @@ Orientational_Correlation::Orientational_Correlation() { n_times=0; - correlation = new float [n_times]; - weighting = new float [n_times]; + correlation = new float * [n_times]; + weighting = new float * [n_times]; overall_correlation = 0; weighting = 0; correlated_vector.set(0,0,0); @@ -37,12 +37,14 @@ Orientational_Correlation::Orientational_Correlation(const Orientational_Correla n_times = copy.n_times; multibody_list = copy.multibody_list; - correlation = new float [n_times]; - weighting = new float [n_times]; + correlation = new float * [n_times]; + weighting = new float * [n_times]; for(timeii=0;timeiishow_n_timesteps(); - correlation = new float [n_times]; - weighting = new float [n_times]; + correlation = new float * [n_times]; + weighting = new float * [n_times]; for(timeii=0;timeiishow_n_timesteps(); - correlation = new float [n_times]; - weighting = new float [n_times]; + correlation = new float * [n_times]; + weighting = new float * [n_times]; for(timeii=0;timeiiisThreadSafe()) for(int timeii=0; timeiishow_n_timesteps(); timeii++) { - weighting[timeii]+=multibody_list->show_n_multibodies(timeii); + weighting[timeii][0]+=multibody_list->show_n_multibodies(timeii); multibody_list->listloop(this,0, timeii, 0); } postprocess_list(); @@ -129,7 +138,7 @@ void Orientational_Correlation::analyze(Multibody_List * t_list) void Orientational_Correlation::listkernel(Multibody* current_multibody, int timegapii,int thisii, int nextii) { float dotproduct= (((*current_multibody)(1)->show_unwrapped(thisii)-(*current_multibody)(0)->show_unwrapped(thisii)).unit_vector())&((correlated_vector).unit_vector()); //compute dot product between unit vectors at initial and later times - correlation[thisii]+=0.5*(3.0*dotproduct*dotproduct - 1.0); //increment baf by second legendre polynomial of dot product above + correlation[thisii][0]+=0.5*(3.0*dotproduct*dotproduct - 1.0); //increment baf by second legendre polynomial of dot product above } @@ -138,9 +147,9 @@ void Orientational_Correlation::postprocess_list() float cumulative_weighting=0; for(int timeii=0;timeiishow_time(timeii)<<"\t"<show_time(timeii)<<"\t"<show_time(timeii)<<"\t"<show_time(timeii)<<"\t"<show_n_trajectories(system_time(timeii)); - n_atoms_j[timeii]=trajectory_list2->show_n_trajectories(system_time(timeii)); + n_atoms_i[timeii][0]=trajectory_list->show_n_trajectories(system_time(timeii)); + n_atoms_j[timeii][0]=trajectory_list2->show_n_trajectories(system_time(timeii)); trajectory_list->listloop(this,0, timeii, 0); } @@ -266,14 +270,16 @@ void Radial_Count::bin(int timestep, float distance) { boxsize = system->size(system_time(timeii)); boxvolume=boxsize.show_x()*boxsize.show_y()*boxsize.show_z(); - rhoj=float(n_atoms_j[timeii])/boxvolume; - n_i_total+=n_atoms_i[timeii]; + rhoj=float(n_atoms_j[timeii][0])/boxvolume; + n_i_total+=n_atoms_i[timeii][0]; for(binii=0;biniishow_n_trajectories(system_time(timeii)); - n_atoms_j[timeii]=trajectory_list2->show_n_trajectories(system_time(timeii)); + n_atoms_i[timeii][0]=trajectory_list->show_n_trajectories(system_time(timeii)); + n_atoms_j[timeii][0]=trajectory_list2->show_n_trajectories(system_time(timeii)); trajectory_list->listloop(this,0, timeii, 0); } @@ -266,14 +292,17 @@ void Radial_Distribution_Function::bin(int timestep, float distance) { boxsize = system->size(system_time(timeii)); boxvolume=boxsize.show_x()*boxsize.show_y()*boxsize.show_z(); - rhoj=float(n_atoms_j[timeii])/boxvolume; - n_i_total+=n_atoms_i[timeii]; + rhoj=float(n_atoms_j[timeii][0])/boxvolume; + n_i_total+=n_atoms_i[timeii][0]; for(binii=0;biniisize(system_time(timeii)); boxvolume=boxsize.show_x()*boxsize.show_y()*boxsize.show_z(); - mean_rho+=float(n_atoms_j[timeii])/boxvolume/n_times; + mean_rho+=float(n_atoms_j[timeii][0])/boxvolume/n_times; } min_k = 2.*PI/max_distance; diff --git a/src/radial_distribution_function.h b/src/radial_distribution_function.h index eaecd87..754ea4c 100644 --- a/src/radial_distribution_function.h +++ b/src/radial_distribution_function.h @@ -21,8 +21,8 @@ class Radial_Distribution_Function: public Analysis_Onetime int n_times; float ** time_rdf; float * mean_rdf; - int * n_atoms_i; - int * n_atoms_j; + int ** n_atoms_i; + int ** n_atoms_j; public: @@ -53,7 +53,7 @@ class Radial_Distribution_Function: public Analysis_Onetime -// bool isThreadSafe(){return true;}; + bool isThreadSafe(){return true;}; }; } diff --git a/src/space_time_correlation_function.cpp b/src/space_time_correlation_function.cpp index a6b0eb1..55a17b3 100644 --- a/src/space_time_correlation_function.cpp +++ b/src/space_time_correlation_function.cpp @@ -36,6 +36,7 @@ void Space_Time_Correlation_Function::clear_memory() for(ii=0;iiisThreadSafe()) // TODO: Test if we can use the old loop + #pragma omp parallel for schedule(static) if(analysis->isThreadSafe()) // TODO: Test if we can use the old loop for(int timegapii=0;timegapiilist_displacementkernel(timegapii,thisii,nextii); //#pragma omp atomic displacement_count++; @@ -4136,7 +4134,7 @@ void System::displacement_list(Analysis* analysis, bool fullblock)const { - //#pragma omp parallel for schedule(dynamic) if(analysis->isThreadSafe()) // This makes this loop execute in parallel, splitting by time values. + #pragma omp parallel for schedule(dynamic) if(analysis->isThreadSafe()) // This makes this loop execute in parallel, splitting by time values. for(int timegapii=n_exponential_steps; timegapiiisThreadSafe()) // TODO: Test if we can use the old loop for(int timegapii=0;timegapiilist_displacementkernel(timegapii,thisii,nextii); - #pragma omp atomic displacement_count++; abort = (displacement_count == displacement_limit && displacement_limit != 0); - #pragma omp flush(abort) } //if(displacement_count == displacement_limit) break; // cout << thisii << "\t" << nextii << "\n"; @@ -4252,10 +4246,8 @@ void System::displacement_list(Multibody_Analysis* analysis, bool fullblock)cons int thisii = n_exponential_steps*blockii+expii+int(frt); int nextii = thisii + n_exponential_steps*block_timegapii; analysis->list_displacementkernel(timegapii,thisii,nextii); - #pragma omp atomic displacement_count++; abort = (displacement_count == displacement_limit && displacement_limit != 0); - #pragma omp flush(abort) // analysis->list_displacementkernel(timegapii,thisii,nextii); // displacement_count++; // // if (displacement_count == displacement_limit) break; diff --git a/src/van_hove_distinct.cpp b/src/van_hove_distinct.cpp index 7ef0d75..dc0a5f4 100644 --- a/src/van_hove_distinct.cpp +++ b/src/van_hove_distinct.cpp @@ -11,7 +11,7 @@ #include "generated/version.h" using namespace std; - +const int PAD = 16; /*----------------------------------------------------------------------------------------------*/ @@ -26,8 +26,9 @@ Van_Hove_Distinct::Van_Hove_Distinct() bin_size = 0; timetable = 0; correlation = new float * [1]; - weighting = new int [0]; - correlation[0]=new float [0]; + correlation[0]= new float [0]; + weighting = new int * [1]; + weighting[0] = new int[0]; } @@ -56,10 +57,11 @@ Van_Hove_Distinct::Van_Hove_Distinct(System*sys, Trajectory_List_Bins binnedtraj // currentlists = new Trajectory_List [2]; correlation = new float * [n_times]; - weighting = new int [n_times]; + weighting = new int * [n_times]; for(timeii=0;timeiishow_n_trajectories(thisii); + weighting[timegapii][0]+=trajectory_list->show_n_trajectories(thisii); // current_list1->listloop(this,timegapii, thisii, nextii); trajectory_list->listloop(this,timegapii, thisii, nextii); } @@ -202,12 +207,12 @@ void Van_Hove_Distinct::listkernel2(Trajectory* traj1, Trajectory* traj2, int ti float distance; if(traj1!=traj2) { - (traj2->show_coordinate(nextii)-(traj1->show_coordinate(thisii))).length_unwrapped(system->size()); //calculate shortest distance between two coordinates, taking into account periodic boundaries + distance = (traj2->show_coordinate(nextii)-(traj1->show_coordinate(thisii))).length_unwrapped(system->size()); //calculate shortest distance between two coordinates, taking into account periodic boundaries bin(timegapii,distance); } else { - weighting[timegapii]--; + //weighting[timegapii][0]--; } } diff --git a/src/van_hove_distinct.h b/src/van_hove_distinct.h index 74449d9..4885da2 100644 --- a/src/van_hove_distinct.h +++ b/src/van_hove_distinct.h @@ -36,7 +36,7 @@ class Van_Hove_Distinct: public Space_Time_Correlation_Function void list_displacementkernel(int,int,int); void listkernel(Trajectory*, int, int, int); void listkernel2(Trajectory*, Trajectory*, int, int, int); - + bool isThreadSafe(){return true;} }; } diff --git a/src/van_hove_self.cpp b/src/van_hove_self.cpp index 6fdb980..69f668e 100644 --- a/src/van_hove_self.cpp +++ b/src/van_hove_self.cpp @@ -11,7 +11,6 @@ using namespace std; - Van_Hove_Self::Van_Hove_Self(System* sys, int bin_count, float value_max) { initialize(sys, bin_count, value_max); @@ -29,7 +28,8 @@ Van_Hove_Self::Van_Hove_Self() timetable=0; correlation = new float * [1]; - weighting = new int [0]; + weighting = new int * [1]; + weighting[0] = new int [0]; correlation[0]=new float [0]; } @@ -62,14 +62,15 @@ void Van_Hove_Self::initialize(System* sys, int bin_count, float value_max) //allocate memory for van hove self-correlation function and weighting and initialize to zero correlation = new float * [n_times]; - weighting = new int [n_times]; + weighting = new int * [n_times]; for(timeii=0;timeiishow_unwrapped(nextii)-current_trajectory->show_unwrapped(thisii)).length(); + Coordinate this_c = current_trajectory->show_unwrapped(nextii); + Coordinate next_c = current_trajectory->show_unwrapped(thisii); + float distance = (next_c - this_c).length(); bin(timegapii,distance); - weighting[timegapii]++; - - + weighting[timegapii][0]++; } diff --git a/src/van_hove_self.h b/src/van_hove_self.h index efe4714..e56f774 100644 --- a/src/van_hove_self.h +++ b/src/van_hove_self.h @@ -26,7 +26,7 @@ class Van_Hove_Self: public Space_Time_Correlation_Function void analyze(Trajectory_List * t_list); void list_displacementkernel(int,int,int); void listkernel (Trajectory* current_trajectory, int timegapii, int thisii, int nextii); - + bool isThreadSafe(){return true;} }; }