Add plugin for extended sampling

doc/Plugins_Summary.dox

@@ -23,6 +23,7 @@ Dropaccelerator | utility for droplet simulations | Speeds up particles that are
 Dropaligner | utility for droplet simulations | aligns the center of mass of particles within a spherical region with a specified center of mass in a similar fashion as the COMaligner. Used to hold a spherical droplet in position.
 EnergyLogWriter | print, file, thermodynamic macroscopic quantities | Write the global N, U_pot, U_kin, U_kinTrans, U_kinRot, T, p to file. Duplication with ResultWriter?
 ExamplePlugin | example, sample, illustrate, usage | Print user-specifed text at a user-specified position in the code.
+ExtendedProfileSampling | sampling, sample, profiles | Sample several regular quantities, e.g. temperature and density, besides more exotic quantities, e.g. chemical potential and energy flux, in spatially resolved profiles along the y-axis.
 ExtractPhase | extract, phase, interface, vapor-liquid equilibrium (VLE), local density | Extracts the liquid phase from a vapor-liquid equilibrium (VLE) simulation by evaluating the local density surrounding individual particles.
 FixRegion | todo | Fixes particles within a cuboid region to their respective positions. Fixed particles still exert forces onto other particles but they do not move.
 FlopRateWriter | performance, FLOP/sec | Calculate the floating point operations per second at which the simulation runs. Only FLOPs in the force calculation are calculated.

examples/run-validation.sh

@@ -0,0 +1,31 @@
+#!/bin/bash
+# Testscript running validation tests according to GitHub Actions
+#
+# Copyright (c) 2023 Simon Homes
+#
+
+repoPath="$PWD/.."  # Path to root directory
+
+export OMP_NUM_THREADS=2
+
+logfileName="${repoPath}/examples/run-validation.log"
+rm -f ${logfileName}
+
+IFS=$'\n'
+for i in $(cat "${repoPath}/examples/example-list.txt" )
+do
+	# skip if comment or empty line
+	if [[ $i == \#* || -z "$i" ]]
+	then
+		continue
+	fi
+	cd ${repoPath}/examples/$(dirname $i)
+
+	# run the examples
+	printf "Running example: ${i} ... " | tee -a ${logfileName}
+	EXE=${repoPath}/build/src/MarDyn
+
+	mpirun --oversubscribe -np 4 ${EXE} $(basename $i) --steps=20 | tee -a ${logfileName}
+	printf "done\n"
+done
+

src/Simulation.cpp

@@ -552,7 +552,7 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 			{
 				delete _longRangeCorrection;
 				Log::global_log->info() << "Initializing homogeneous LRC." << std::endl;
-				_longRangeCorrection = new Homogeneous(_cutoffRadius, _LJCutoffRadius, _domain, global_simulation);
+				_longRangeCorrection = new Homogeneous(_cutoffRadius, _LJCutoffRadius, _domain, _moleculeContainer, global_simulation);
 			}
 			else if("none" == type)
 			{
@@ -568,7 +568,7 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 		} else {
 			delete _longRangeCorrection;
 			Log::global_log->info() << "Initializing default homogeneous LRC, as no LRC was defined." << std::endl;
-			_longRangeCorrection = new Homogeneous(_cutoffRadius, _LJCutoffRadius, _domain, global_simulation);
+			_longRangeCorrection = new Homogeneous(_cutoffRadius, _LJCutoffRadius, _domain, _moleculeContainer, global_simulation);
 		}
 
 		xmlconfig.changecurrentnode(".."); /* algorithm section */

src/Simulation.h

@@ -226,6 +226,9 @@ class Simulation {
 	/** Get pointer to the molecule container */
 	ParticleContainer* getMoleculeContainer() { return _moleculeContainer; }
 
+	/** Get pointer to longrange correction */
+	LongRangeCorrection* getLongRangeCorrection() { return _longRangeCorrection; }
+
 	/** Set the number of time steps to be performed in the simulation */
 	void setNumTimesteps( unsigned long steps ) { _numberOfTimesteps = steps; }
 	/** Get the number of time steps to be performed in the simulation */

src/bhfmm/cellProcessors/VectorizedChargeP2PCellProcessor.h

@@ -132,7 +132,9 @@ class VectorizedChargeP2PCellProcessor {
 		const RealCalcVec& r2_x, const RealCalcVec& r2_y, const RealCalcVec& r2_z,
 		const RealCalcVec& qjj,
 		RealCalcVec& f_x, RealCalcVec& f_y, RealCalcVec& f_z,
-		RealAccumVec& V_x, RealAccumVec& V_y, RealAccumVec& V_z,
+		RealAccumVec& V_xx, RealAccumVec& V_yy, RealAccumVec& V_zz,
+		RealAccumVec& V_xy, RealAccumVec& V_xz, RealAccumVec& V_yz,
+		RealAccumVec& V_yx, RealAccumVec& V_zx, RealAccumVec& V_zy,
 		RealAccumVec& sum_upotXpoles, RealAccumVec& sum_virial,
 		const MaskCalcVec& forceMask);
 

src/bhfmm/cellProcessors/VectorizedLJP2PCellProcessor.cpp

@@ -179,7 +179,9 @@ void VectorizedLJP2PCellProcessor :: _loopBodyLJ(
 		const RealCalcVec& m2_r_x, const RealCalcVec& m2_r_y, const RealCalcVec& m2_r_z,
 		const RealCalcVec& r2_x, const RealCalcVec& r2_y, const RealCalcVec& r2_z,
 		RealCalcVec& f_x, RealCalcVec& f_y, RealCalcVec& f_z,
-		RealAccumVec& V_x, RealAccumVec& V_y, RealAccumVec& V_z,
+		RealAccumVec& V_xx, RealAccumVec& V_yy, RealAccumVec& V_zz,
+		RealAccumVec& V_xy, RealAccumVec& V_xz, RealAccumVec& V_yz,
+		RealAccumVec& V_yx, RealAccumVec& V_zx, RealAccumVec& V_zy,
 		RealAccumVec& sum_upot6lj, RealAccumVec& sum_virial,
 		const MaskCalcVec& forceMask,
 		const RealCalcVec& eps_24, const RealCalcVec& sig2,
@@ -211,9 +213,15 @@ void VectorizedLJP2PCellProcessor :: _loopBodyLJ(
 	const RealCalcVec m_dy = m1_r_y - m2_r_y;//1FP (virial) (does not count)
 	const RealCalcVec m_dz = m1_r_z - m2_r_z;//1FP (virial) (does not count)
 
-	V_x = RealAccumVec::convertCalcToAccum(m_dx * f_x);//1FP (virial)
-	V_y = RealAccumVec::convertCalcToAccum(m_dy * f_y);//1FP (virial)
-	V_z = RealAccumVec::convertCalcToAccum(m_dz * f_z);//1FP (virial)
+	V_xx = RealAccumVec::convertCalcToAccum(m_dx * f_x);
+	V_yy = RealAccumVec::convertCalcToAccum(m_dy * f_y);
+	V_zz = RealAccumVec::convertCalcToAccum(m_dz * f_z);
+	V_xy = RealAccumVec::convertCalcToAccum(m_dx * f_y);
+	V_xz = RealAccumVec::convertCalcToAccum(m_dx * f_z);
+	V_yz = RealAccumVec::convertCalcToAccum(m_dy * f_z);
+	V_yx = RealAccumVec::convertCalcToAccum(m_dy * f_x);
+	V_zx = RealAccumVec::convertCalcToAccum(m_dz * f_x);
+	V_zy = RealAccumVec::convertCalcToAccum(m_dz * f_y);
 
 	// Check if we have to add the macroscopic values up
 	if (calculateMacroscopic) {
@@ -224,7 +232,7 @@ void VectorizedLJP2PCellProcessor :: _loopBodyLJ(
 
 		sum_upot6lj = sum_upot6lj + upot_masked_accum;//1FP (sum macro)
 
-		sum_virial = sum_virial +  V_x + V_y + V_z;//1 FP (sum macro) + 2 FP (virial)
+		sum_virial = sum_virial +  V_xx + V_yy + V_zz;//1 FP (sum macro) + 2 FP (virial)
 	}
 }
 
@@ -254,9 +262,15 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 	     vcp_real_accum * const soa1_ljc_f_x = soa1.getBeginAccum(QuantityType::FORCE, LJC, Coordinate::X);
 	     vcp_real_accum * const soa1_ljc_f_y = soa1.getBeginAccum(QuantityType::FORCE, LJC, Coordinate::Y);
 	     vcp_real_accum * const soa1_ljc_f_z = soa1.getBeginAccum(QuantityType::FORCE, LJC, Coordinate::Z);
-	     vcp_real_accum * const soa1_ljc_V_x = soa1.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::X);
-	     vcp_real_accum * const soa1_ljc_V_y = soa1.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::Y);
-	     vcp_real_accum * const soa1_ljc_V_z = soa1.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::Z);
+	     vcp_real_accum * const soa1_ljc_V_xx = soa1.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::X);
+		 vcp_real_accum * const soa1_ljc_V_yy = soa1.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::Y);
+		 vcp_real_accum * const soa1_ljc_V_zz = soa1.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::Z);
+		 vcp_real_accum * const soa1_ljc_V_xy = soa1.getBeginAccum(QuantityType::VIRIALND1, LJC, Coordinate::X);
+		 vcp_real_accum * const soa1_ljc_V_xz = soa1.getBeginAccum(QuantityType::VIRIALND1, LJC, Coordinate::Y);
+		 vcp_real_accum * const soa1_ljc_V_yz = soa1.getBeginAccum(QuantityType::VIRIALND1, LJC, Coordinate::Z);
+		 vcp_real_accum * const soa1_ljc_V_yx = soa1.getBeginAccum(QuantityType::VIRIALND2, LJC, Coordinate::X);
+		 vcp_real_accum * const soa1_ljc_V_zx = soa1.getBeginAccum(QuantityType::VIRIALND2, LJC, Coordinate::Y);
+		 vcp_real_accum * const soa1_ljc_V_zy = soa1.getBeginAccum(QuantityType::VIRIALND2, LJC, Coordinate::Z);
 	const int * const soa1_mol_ljc_num = soa1._mol_ljc_num;
 	const vcp_ljc_id_t * const soa1_ljc_id = soa1._ljc_id;
 
@@ -269,9 +283,15 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 	     vcp_real_accum * const soa2_ljc_f_x = soa2.getBeginAccum(QuantityType::FORCE, LJC, Coordinate::X);
 	     vcp_real_accum * const soa2_ljc_f_y = soa2.getBeginAccum(QuantityType::FORCE, LJC, Coordinate::Y);
 	     vcp_real_accum * const soa2_ljc_f_z = soa2.getBeginAccum(QuantityType::FORCE, LJC, Coordinate::Z);
-	     vcp_real_accum * const soa2_ljc_V_x = soa2.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::X);
-	     vcp_real_accum * const soa2_ljc_V_y = soa2.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::Y);
-	     vcp_real_accum * const soa2_ljc_V_z = soa2.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::Z);
+	     vcp_real_accum * const soa2_ljc_V_xx = soa2.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::X);
+		 vcp_real_accum * const soa2_ljc_V_yy = soa2.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::Y);
+		 vcp_real_accum * const soa2_ljc_V_zz = soa2.getBeginAccum(QuantityType::VIRIAL, LJC, Coordinate::Z);
+		 vcp_real_accum * const soa2_ljc_V_xy = soa2.getBeginAccum(QuantityType::VIRIALND1, LJC, Coordinate::X);
+		 vcp_real_accum * const soa2_ljc_V_xz = soa2.getBeginAccum(QuantityType::VIRIALND1, LJC, Coordinate::Y);
+		 vcp_real_accum * const soa2_ljc_V_yz = soa2.getBeginAccum(QuantityType::VIRIALND1, LJC, Coordinate::Z);
+		 vcp_real_accum * const soa2_ljc_V_yx = soa2.getBeginAccum(QuantityType::VIRIALND2, LJC, Coordinate::X);
+		 vcp_real_accum * const soa2_ljc_V_zx = soa2.getBeginAccum(QuantityType::VIRIALND2, LJC, Coordinate::Y);
+		 vcp_real_accum * const soa2_ljc_V_zy = soa2.getBeginAccum(QuantityType::VIRIALND2, LJC, Coordinate::Z);
 	const vcp_ljc_id_t * const soa2_ljc_id = soa2._ljc_id;
 
 	vcp_lookupOrMask_single* const soa2_ljc_dist_lookup = my_threadData._ljc_dist_lookup;
@@ -320,9 +340,15 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 				RealAccumVec sum_fy1 = RealAccumVec::zero();
 				RealAccumVec sum_fz1 = RealAccumVec::zero();
 
-				RealAccumVec sum_Vx1 = RealAccumVec::zero();
-				RealAccumVec sum_Vy1 = RealAccumVec::zero();
-				RealAccumVec sum_Vz1 = RealAccumVec::zero();
+				RealAccumVec sum_V1_xx = RealAccumVec::zero();
+				RealAccumVec sum_V1_yy = RealAccumVec::zero();
+				RealAccumVec sum_V1_zz = RealAccumVec::zero();
+				RealAccumVec sum_V1_xy = RealAccumVec::zero();
+				RealAccumVec sum_V1_xz = RealAccumVec::zero();
+				RealAccumVec sum_V1_yz = RealAccumVec::zero();
+				RealAccumVec sum_V1_yx = RealAccumVec::zero();
+				RealAccumVec sum_V1_zx = RealAccumVec::zero();
+				RealAccumVec sum_V1_zy = RealAccumVec::zero();
 
 				const RealCalcVec c_r_x1 = RealCalcVec::broadcast(soa1_ljc_r_x + i_ljc_idx);
 				const RealCalcVec c_r_y1 = RealCalcVec::broadcast(soa1_ljc_r_y + i_ljc_idx);
@@ -344,7 +370,7 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 
 						const size_t id_i = soa1_ljc_id[i_ljc_idx];
 						RealCalcVec fx, fy, fz;
-						RealAccumVec Vx, Vy, Vz;
+						RealAccumVec Vxx, Vyy, Vzz, Vxy, Vxz, Vyz, Vyx, Vzx, Vzy;
 
 						RealCalcVec eps_24;
 						RealCalcVec sig2;
@@ -357,7 +383,9 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 							m1_r_x, m1_r_y, m1_r_z, c_r_x1, c_r_y1, c_r_z1,
 							m_r_x2, m_r_y2, m_r_z2, c_r_x2, c_r_y2, c_r_z2,
 							fx, fy, fz,
-							Vx, Vy, Vz,
+							Vxx, Vyy, Vzz,
+							Vxy, Vxz, Vyz,
+							Vyx, Vzx, Vzy,
 							sum_upot6lj, sum_virial,
 							MaskGatherChooser::getForceMask(lookupORforceMask),
 							eps_24, sig2,
@@ -375,13 +403,25 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 						sum_fy1 = sum_fy1 + a_fy;
 						sum_fz1 = sum_fz1 + a_fz;
 
-						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_x, j, Vx, lookupORforceMask);
-						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_y, j, Vy, lookupORforceMask);
-						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_z, j, Vz, lookupORforceMask);
-
-						sum_Vx1 = sum_Vx1 + Vx;
-						sum_Vy1 = sum_Vy1 + Vy;
-						sum_Vz1 = sum_Vz1 + Vz;
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_xx, j, Vxx, lookupORforceMask);
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_yy, j, Vyy, lookupORforceMask);
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_zz, j, Vzz, lookupORforceMask);
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_xy, j, Vxy, lookupORforceMask);
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_xz, j, Vxz, lookupORforceMask);
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_yz, j, Vyz, lookupORforceMask);
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_yx, j, Vyx, lookupORforceMask);
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_zx, j, Vzx, lookupORforceMask);
+						vcp_simd_load_add_store<MaskGatherChooser>(soa2_ljc_V_zy, j, Vzy, lookupORforceMask);
+
+						sum_V1_xx = sum_V1_xx + Vxx;
+						sum_V1_yy = sum_V1_yy + Vyy;
+						sum_V1_zz = sum_V1_zz + Vzz;
+						sum_V1_xy = sum_V1_xy + Vxy;
+						sum_V1_xz = sum_V1_xz + Vxz;
+						sum_V1_yz = sum_V1_yz + Vyz;
+						sum_V1_yx = sum_V1_yx + Vyx;
+						sum_V1_zx = sum_V1_zx + Vzx;
+						sum_V1_zy = sum_V1_zy + Vzy;
 					}
 				}
 #if VCP_VEC_TYPE == VCP_VEC_KNL_GATHER or VCP_VEC_TYPE == VCP_VEC_AVX512F_GATHER
@@ -405,7 +445,7 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 
 						const size_t id_i = soa1_ljc_id[i_ljc_idx];
 						RealCalcVec fx, fy, fz;
-						RealAccumVec Vx, Vy, Vz;
+						RealAccumVec Vxx, Vyy, Vzz, Vxy, Vxz, Vyz, Vyx, Vzx, Vzy;
 
 						RealCalcVec eps_24;
 						RealCalcVec sig2;
@@ -418,7 +458,9 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 							m1_r_x, m1_r_y, m1_r_z, c_r_x1, c_r_y1, c_r_z1,
 							m_r_x2, m_r_y2, m_r_z2, c_r_x2, c_r_y2, c_r_z2,
 							fx, fy, fz,
-							Vx, Vy, Vz,
+							Vxx, Vyy, Vzz,
+							Vxy, Vxz, Vyz,
+							Vyx, Vzx, Vzy,
 							sum_upot6lj, sum_virial,
 							remainderM,//use remainder mask as forcemask
 							eps_24, sig2,
@@ -436,13 +478,25 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 						sum_fy1 = sum_fy1 + a_fy;
 						sum_fz1 = sum_fz1 + a_fz;
 
-						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_x, j, Vx, lookupORforceMask, remainderM);
-						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_y, j, Vy, lookupORforceMask, remainderM);
-						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_z, j, Vz, lookupORforceMask, remainderM);
-
-						sum_Vx1 = sum_Vx1 + Vx;
-						sum_Vy1 = sum_Vy1 + Vy;
-						sum_Vz1 = sum_Vz1 + Vz;
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_xx, j, Vxx, lookupORforceMask, remainderM);
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_yy, j, Vyy, lookupORforceMask, remainderM);
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_zz, j, Vzz, lookupORforceMask, remainderM);
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_xy, j, Vxy, lookupORforceMask, remainderM);
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_xz, j, Vxz, lookupORforceMask, remainderM);
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_yz, j, Vyz, lookupORforceMask, remainderM);
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_yx, j, Vyx, lookupORforceMask, remainderM);
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_zx, j, Vzx, lookupORforceMask, remainderM);
+						vcp_simd_load_add_store_masked<MaskGatherChooser>(soa2_ljc_V_zy, j, Vzy, lookupORforceMask, remainderM);
+
+						sum_V1_xx = sum_V1_xx + Vxx;
+						sum_V1_yy = sum_V1_yy + Vyy;
+						sum_V1_zz = sum_V1_zz + Vzz;
+						sum_V1_xy = sum_V1_xy + Vxy;
+						sum_V1_xz = sum_V1_xz + Vxz;
+						sum_V1_yz = sum_V1_yz + Vyz;
+						sum_V1_yx = sum_V1_yx + Vyx;
+						sum_V1_zx = sum_V1_zx + Vzx;
+						sum_V1_zy = sum_V1_zy + Vzy;
 					}
 				}
 #endif
@@ -451,9 +505,15 @@ void VectorizedLJP2PCellProcessor::_calculatePairs(CellDataSoA & soa1, CellDataS
 				hSum_Add_Store(soa1_ljc_f_y + i_ljc_idx, sum_fy1);
 				hSum_Add_Store(soa1_ljc_f_z + i_ljc_idx, sum_fz1);
 
-				hSum_Add_Store(soa1_ljc_V_x + i_ljc_idx, sum_Vx1);
-				hSum_Add_Store(soa1_ljc_V_y + i_ljc_idx, sum_Vy1);
-				hSum_Add_Store(soa1_ljc_V_z + i_ljc_idx, sum_Vz1);
+				hSum_Add_Store(soa1_ljc_V_xx + i_ljc_idx, sum_V1_xx);
+				hSum_Add_Store(soa1_ljc_V_yy + i_ljc_idx, sum_V1_yy);
+				hSum_Add_Store(soa1_ljc_V_zz + i_ljc_idx, sum_V1_zz);
+				hSum_Add_Store(soa1_ljc_V_xy + i_ljc_idx, sum_V1_xy);
+				hSum_Add_Store(soa1_ljc_V_xz + i_ljc_idx, sum_V1_xz);
+				hSum_Add_Store(soa1_ljc_V_yz + i_ljc_idx, sum_V1_yz);
+				hSum_Add_Store(soa1_ljc_V_yx + i_ljc_idx, sum_V1_yx);
+				hSum_Add_Store(soa1_ljc_V_zx + i_ljc_idx, sum_V1_zx);
+				hSum_Add_Store(soa1_ljc_V_zy + i_ljc_idx, sum_V1_zy);
 
 				i_ljc_idx++;
 			}

src/bhfmm/cellProcessors/VectorizedLJP2PCellProcessor.h

@@ -170,7 +170,9 @@ class VectorizedLJP2PCellProcessor : public CellProcessor {
 			const RealCalcVec& m2_r_x, const RealCalcVec& m2_r_y, const RealCalcVec& m2_r_z,
 			const RealCalcVec& r2_x, const RealCalcVec& r2_y, const RealCalcVec& r2_z,
 			RealCalcVec& f_x, RealCalcVec& f_y, RealCalcVec& f_z,
-			RealAccumVec& V_x, RealAccumVec& V_y, RealAccumVec& V_z,
+			RealAccumVec& V_xx, RealAccumVec& V_yy, RealAccumVec& V_zz,
+			RealAccumVec& V_xy, RealAccumVec& V_xz, RealAccumVec& V_yz,
+			RealAccumVec& V_yx, RealAccumVec& V_zx, RealAccumVec& V_zy,
 			RealAccumVec& sum_upot6lj, RealAccumVec& sum_virial,
 			const MaskCalcVec& forceMask,
 			const RealCalcVec& eps_24, const RealCalcVec& sig2,

src/ensemble/ChemicalPotential.cpp

@@ -469,6 +469,8 @@ void ChemicalPotential::grandcanonicalStep(
 	bool hasInsertion = true;
 	double ins[3];
 	unsigned nextid = 0;
+	double zeroVec[3] = {0.0, 0.0, 0.0};
+	double zeroVec9[9] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 	while (hasDeletion || hasInsertion) {
 		if (hasDeletion) {
 			auto m = this->getDeletion(moleculeContainer, minco, maxco);
@@ -491,10 +493,9 @@ void ChemicalPotential::grandcanonicalStep(
 					// m->upd_cache(); TODO what to do here? somebody deleted the method "upd_cache"!!! why???
 					// reset forces and momenta to zero
 					{
-						double zeroVec[3] = {0.0, 0.0, 0.0};
 						m->setF(zeroVec);
 						m->setM(zeroVec);
-						m->setVi(zeroVec);
+						m->setVi(zeroVec9);
 					}
 
 					this->storeMolecule(*m);
@@ -536,10 +537,9 @@ void ChemicalPotential::grandcanonicalStep(
 			// m->upd_cache(); TODO: what to do here? somebody deleted the method "upd_cache"!!! why??? -- update caches do no longer exist ;)
 			// reset forces and torques to zero
 			if (!this->isWidom()) {
-				double zeroVec[3] = { 0.0, 0.0, 0.0 };
 				mit->setF(zeroVec);
 				mit->setM(zeroVec);
-				mit->setVi(zeroVec);
+				mit->setVi(zeroVec9);
 			}
 			mit->check(nextid);
 #ifndef NDEBUG
@@ -571,17 +571,15 @@ void ChemicalPotential::grandcanonicalStep(
 #endif
 			if (accept) {
 				this->_localInsertionsMinusDeletions++;
-				double zeroVec[3] = { 0.0, 0.0, 0.0 };
-				tmp.setVi(zeroVec);
+				tmp.setVi(zeroVec9);
 				bool inBoxCheckedAlready = false, checkWhetherDuplicate = false, rebuildCaches = true;
 				moleculeContainer->addParticle(tmp, inBoxCheckedAlready, checkWhetherDuplicate, rebuildCaches);
 			} else {
 				// moleculeContainer->deleteMolecule(m->getID(), m->r(0), m->r(1), m->r(2));
 //				moleculeContainer->_cells[cellid].deleteMolecule(m->getID());
-				double zeroVec[3] = { 0.0, 0.0, 0.0 };
 				mit->setF(zeroVec);
 				mit->setM(zeroVec);
-				mit->setVi(zeroVec);
+				mit->setVi(zeroVec9);
 				mit->check(mit->getID());
 //				moleculeContainer->_particles.erase(mit);
 			}