Fix velocity assignment in CubicGridGenerator (and others)

src/io/CubicGridGeneratorInternal.cpp

@@ -17,7 +17,6 @@
 #include "parallel/DomainDecompBase.h"
 #include "particleContainer/ParticleContainer.h"
 #include "utils/Logger.h"
-#include "utils/Random.h"
 #include "utils/mardyn_assert.h"
 
 #include <cmath>
@@ -86,8 +85,7 @@ unsigned long CubicGridGeneratorInternal::readPhaseSpace(ParticleContainer *part
 		global_simulation->getEnsemble()->getComponents()->at(1).updateMassInertia();
 	}
 
-	unsigned long int id = particleContainer->initCubicGrid(
-		numMoleculesPerDim, simBoxLength, static_cast<size_t>(domainDecomp->getRank()) * mardyn_get_max_threads());
+	unsigned long int id = particleContainer->initCubicGrid(numMoleculesPerDim, simBoxLength);
 
 	Log::global_log->info() << "Finished reading molecules: 100%" << std::endl;
 

src/io/IOHelpers.cpp

@@ -1,7 +1,9 @@
 #include "IOHelpers.h"
 
+#include "Simulation.h"
 #include "parallel/DomainDecompBase.h"
 #include "particleContainer/ParticleContainer.h"
+#include "utils/generator/EqualVelocityAssigner.h"
 
 void IOHelpers::removeMomentum(ParticleContainer* particleContainer, const std::vector<Component>& components,
 							   DomainDecompBase* domainDecomp) {
@@ -38,9 +40,9 @@ void IOHelpers::removeMomentum(ParticleContainer* particleContainer, const std::
 	Log::global_log->info() << "momentumsum prior to removal: " << momentum_sum[0] << " " << momentum_sum[1] << " "
 							<< momentum_sum[2] << std::endl;
 	Log::global_log->info() << "mass_sum: " << mass_sum << std::endl;
-	double v_sub0 = momentum_sum[0] / mass_sum;
-	double v_sub1 = momentum_sum[1] / mass_sum;
-	double v_sub2 = momentum_sum[2] / mass_sum;
+	const double v_sub0 = momentum_sum[0] / mass_sum;
+	const double v_sub1 = momentum_sum[1] / mass_sum;
+	const double v_sub2 = momentum_sum[2] / mass_sum;
 #if defined(_OPENMP)
 #pragma omp parallel
 #endif
@@ -69,38 +71,56 @@ void IOHelpers::removeMomentum(ParticleContainer* particleContainer, const std::
 	}
 	Log::global_log->info() << "momentumsum after removal: " << momentum_sum[0] << " " << momentum_sum[1] << " "
 							<< momentum_sum[2] << std::endl;
-}
 
-void IOHelpers::initializeVelocityAccordingToTemperature(ParticleContainer* particleContainer,
-														 DomainDecompBase* /*domainDecomp*/, double temperature) {
-	Random rng;
-	for (auto iter = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); iter.isValid(); ++iter) {
-		auto velocity = getRandomVelocityAccordingToTemperature(temperature, rng);
-		for (unsigned short d = 0; d < 3; ++d) {
-			iter->setv(d, velocity[d]);
+	// After subtraction, the temperature was effectively changed
+	// Therefore, the velocities have to be scaled again to match target temperature
+	const double temperature_target = global_simulation->getEnsemble()->T();
+	double ekin_sum = 0.0;
+	unsigned long numDOFs = 0;
+
+#if defined(_OPENMP)
+#pragma omp parallel reduction(+ : ekin_sum, numDOFs)
+#endif
+	{
+		for (auto molecule = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); molecule.isValid();
+			 ++molecule) {
+			ekin_sum += molecule->U_kin();
+			numDOFs += (3+molecule->component()->getRotationalDegreesOfFreedom());
 		}
 	}
-}
 
-std::array<double, 3> IOHelpers::getRandomVelocityAccordingToTemperature(double temperature, Random& RNG) {
-	std::array<double, 3> ret{};
-	double dotprod_v = 0;
-	// Doing this multiple times ensures that we have equally distributed directions for the velocities!
-	do {
-		// Velocity
-		for (int dim = 0; dim < 3; dim++) {
-			ret[dim] = RNG.uniformRandInRange(-0.5f, 0.5f);
-		}
-		dotprod_v = ret[0] * ret[0] + ret[1] * ret[1] + ret[2] * ret[2];
-	} while (dotprod_v < 0.0625);
+	domainDecomp->collCommInit(2);
+	domainDecomp->collCommAppendDouble(ekin_sum);
+	domainDecomp->collCommAppendUnsLong(numDOFs);
+	domainDecomp->collCommAllreduceSum();
+	ekin_sum = domainDecomp->collCommGetDouble();
+	numDOFs = domainDecomp->collCommGetUnsLong();
+	domainDecomp->collCommFinalize();
 
-	// Velocity Correction
-	double vCorr = sqrt(3. * temperature / dotprod_v);
-	for (unsigned int i = 0; i < ret.size(); i++) {
-		ret[i] *= vCorr;
+	const double temperature_current = 2*ekin_sum/numDOFs;
+	const double scaleFactor = std::sqrt(temperature_target/temperature_current);
+	Log::global_log->info() << "current temperature: " << temperature_current
+							<< "scale velocities by: " << scaleFactor << std::endl;
+#if defined(_OPENMP)
+#pragma omp parallel
+#endif
+	{
+		for (auto mol = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); mol.isValid();
+			 ++mol) {
+			for (int d = 0; d < 3; d++) {
+				// Scale velocities to match target temperature
+				mol->setv(d,mol->v(d)*scaleFactor);
+			}
+		}
 	}
+}
 
-	return ret;
+void IOHelpers::initializeVelocityAccordingToTemperature(ParticleContainer* particleContainer,
+														 DomainDecompBase* /*domainDecomp*/, double temperature) {
+	EqualVelocityAssigner eqVeloAssigner(temperature);
+	for (auto iter = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); iter.isValid(); ++iter) {
+		eqVeloAssigner.assignVelocity(&(*iter));
+	}
 }
 
 unsigned long IOHelpers::makeParticleIdsUniqueAndGetTotalNumParticles(ParticleContainer* particleContainer,

src/io/IOHelpers.h

@@ -1,7 +1,6 @@
 #pragma once
 
 #include "molecules/Component.h"
-#include "utils/Random.h"
 
 class ParticleContainer;
 class DomainDecompBase;
@@ -46,14 +45,4 @@ void initializeVelocityAccordingToTemperature(ParticleContainer* particleContain
 unsigned long makeParticleIdsUniqueAndGetTotalNumParticles(ParticleContainer* particleContainer,
 														   DomainDecompBase* domainDecomp);
 
-/**
- * Generates a velocity vector according to the given temperature.
- * The generated velocity vector lies on the surface of a sphere with radius sqrt(3. * temperature).
- * If this function is called multiple times, the generated velocity vectors will be equally distributed on the surface
- * of the sphere.
- * @param temperature
- * @param RNG
- * @return
- */
-std::array<double, 3> getRandomVelocityAccordingToTemperature(double temperature, Random& RNG);
 }  // namespace IOHelpers

src/io/PerCellGenerator.cpp

@@ -2,6 +2,7 @@
 
 #include <iostream>
 #include <random>
+#include <array>
 
 #include "Domain.h"
 #include "IOHelpers.h"

src/particleContainer/AutoPasContainer.cpp

@@ -14,6 +14,7 @@
 #include "autopas/utils/StringUtils.h"
 #include "autopas/utils/logging/Logger.h"
 #include "parallel/DomainDecompBase.h"
+#include "utils/generator/EqualVelocityAssigner.h"
 
 // Declare the main AutoPas class and the iteratePairwise() methods with all used functors as extern template
 // instantiation. They are instantiated in the respective cpp file inside the templateInstantiations folder.
@@ -674,33 +675,10 @@ std::variant<ParticleIterator, SingleCellIterator<ParticleCell>> AutoPasContaine
 }
 
 unsigned long AutoPasContainer::initCubicGrid(std::array<unsigned long, 3> numMoleculesPerDimension,
-											  std::array<double, 3> simBoxLength, size_t seed_offset) {
+											  std::array<double, 3> simBoxLength) {
 
-	// Stolen from ParticleCellBase.cpp
-	auto getRandomVelocity = [](auto temperature, Random &RNG) {
-		using T = vcp_real_calc;
-		std::array<T,3> ret{};
-
-		// Velocity
-		for (int dim = 0; dim < 3; dim++) {
-			ret[dim] = RNG.uniformRandInRange(-0.5f, 0.5f);
-		}
-		T dotprod_v = 0;
-		for (unsigned int i = 0; i < ret.size(); i++) {
-			dotprod_v += ret[i] * ret[i];
-		}
-		// Velocity Correction
-		const T three = static_cast<T>(3.0);
-		T vCorr = sqrt(three * temperature / dotprod_v);
-		for (unsigned int i = 0; i < ret.size(); i++) {
-			ret[i] *= vCorr;
-		}
-
-		return ret;
-	};
-
-	Random myRNG{static_cast<int>(seed_offset) + mardyn_get_thread_num()};
 	vcp_real_calc T = global_simulation->getEnsemble()->T();
+	EqualVelocityAssigner eqVeloAssigner(T);
 
 	const std::array<double, 3> spacing = autopas::utils::ArrayMath::div(simBoxLength,
 																		 autopas::utils::ArrayUtils::static_cast_copy_array<double>(
@@ -715,10 +693,10 @@ unsigned long AutoPasContainer::initCubicGrid(std::array<unsigned long, 3> numMo
 				const double posY = offset[1] + spacing[1] * y;
 				for (int x = 0; x < numMoleculesPerDimension[0]; ++x) {
 					const double posX = offset[0] + spacing[0] * x;
-					std::array<vcp_real_calc, 3> v = getRandomVelocity(T, myRNG);
+					// Init molecule with zero velocity and use the EqualVelocityAssigner in the next step
 					Molecule m(numMolecules++, &(global_simulation->getEnsemble()->getComponents()->at(0)), posX, posY,
-							   posZ, v[0], v[1],
-							   v[2]);
+							   posZ, 0.0, 0.0, 0.0);
+					eqVeloAssigner.assignVelocity(&m);
 					addParticle(m, true, false, false);
 				}
 			}

src/particleContainer/AutoPasContainer.h

@@ -111,7 +111,7 @@ class AutoPasContainer : public ParticleContainer {
 	getMoleculeAtPosition(const double pos[3]) override;
 
 	unsigned long initCubicGrid(std::array<unsigned long, 3> numMoleculesPerDimension,
-								std::array<double, 3> simBoxLength, size_t seed_offset) override;
+								std::array<double, 3> simBoxLength) override;
 
 	double *getCellLength() override;
 

src/particleContainer/LinkedCells.cpp

@@ -13,7 +13,6 @@
 #include "particleContainer/adapter/ParticlePairs2PotForceAdapter.h"
 #include "particleContainer/handlerInterfaces/ParticlePairsHandler.h"
 #include "utils/Logger.h"
-#include "utils/Random.h"
 #include "utils/mardyn_assert.h"
 #include "utils/GetChunkSize.h"
 
@@ -939,7 +938,7 @@ void LinkedCells::getCellIndicesOfRegion(const double startRegion[3], const doub
 	endIndex = getCellIndexOfPoint(endRegion);
 }
 
-unsigned long LinkedCells::initCubicGrid(std::array<unsigned long, 3> numMoleculesPerDimension, std::array<double, 3> simBoxLength, size_t seed_offset) {
+unsigned long LinkedCells::initCubicGrid(std::array<unsigned long, 3> numMoleculesPerDimension, std::array<double, 3> simBoxLength) {
 	const unsigned long numCells = _cells.size();
 
 	std::vector<unsigned long> numMoleculesPerThread;
@@ -953,15 +952,13 @@ unsigned long LinkedCells::initCubicGrid(std::array<unsigned long, 3> numMolecul
 		const int myID = mardyn_get_thread_num();
 		const unsigned long myStart = numCells * myID / numThreads;
 		const unsigned long myEnd = numCells * (myID + 1) / numThreads;
-		const int seed = seed_offset + myID;
 
 		unsigned long numMoleculesByThisThread = 0;
-		Random threadPrivateRNG{seed};
 
 		// manual "static" scheduling important, because later this thread needs to traverse the same cells
 		for (unsigned long cellIndex = myStart; cellIndex < myEnd; ++cellIndex) {
 			ParticleCell & cell = _cells[cellIndex];
-			numMoleculesByThisThread += cell.initCubicGrid(numMoleculesPerDimension, simBoxLength, threadPrivateRNG);
+			numMoleculesByThisThread += cell.initCubicGrid(numMoleculesPerDimension, simBoxLength);
 		}
 
 		// prefix sum of numMoleculesByThisThread

src/particleContainer/LinkedCells.h

@@ -260,7 +260,7 @@ class LinkedCells : public ParticleContainer {
 	bool requiresForceExchange() const override; // new
 
 	unsigned long initCubicGrid(std::array<unsigned long, 3> numMoleculesPerDimension,
-								std::array<double, 3> simBoxLength, size_t seed_offset) override;
+								std::array<double, 3> simBoxLength) override;
 
 	std::vector<unsigned long> getParticleCellStatistics() override;
 

src/particleContainer/ParticleCellBase.cpp

@@ -7,8 +7,8 @@
 
 #include "ParticleCellBase.h"
 #include "ensemble/EnsembleBase.h"
-#include "utils/Random.h"
 #include "Simulation.h"
+#include "utils/generator/EqualVelocityAssigner.h"
 
 CellBorderAndFlagManager ParticleCellBase::_cellBorderAndFlagManager;
 
@@ -29,28 +29,6 @@ bool ParticleCellBase::deleteMoleculeByID(unsigned long molid) {
 	return found;
 }
 
-template <typename T>
-std::array<T, 3> getRandomVelocity(T temperature, Random& RNG) {
-	std::array<T,3> ret{};
-
-	// Velocity
-	for (int dim = 0; dim < 3; dim++) {
-		ret[dim] = RNG.uniformRandInRange(-0.5f, 0.5f);
-	}
-	T dotprod_v = 0;
-	for (unsigned int i = 0; i < ret.size(); i++) {
-		dotprod_v += ret[i] * ret[i];
-	}
-	// Velocity Correction
-	const T three = static_cast<T>(3.0);
-	T vCorr = sqrt(three * temperature / dotprod_v);
-	for (unsigned int i = 0; i < ret.size(); i++) {
-		ret[i] *= vCorr;
-	}
-
-	return ret;
-}
-
 template <typename T>
 bool PositionIsInBox1D(const T l, const T u, const T r) {
 #ifdef __INTEL_COMPILER
@@ -71,7 +49,7 @@ bool PositionIsInBox3D(const T l[3], const T u[3], const T r[3]) {
 }
 
 unsigned long ParticleCellBase::initCubicGrid(const std::array<unsigned long, 3> &numMoleculesPerDimension,
-											  const std::array<double, 3> &simBoxLength, Random &RNG) {
+											  const std::array<double, 3> &simBoxLength) {
 
 	std::array<double, 3> spacing{};
 	std::array<double, 3> origin1{}; // origin of the first DrawableMolecule
@@ -83,6 +61,7 @@ unsigned long ParticleCellBase::initCubicGrid(const std::array<unsigned long, 3>
 	}
 
 	vcp_real_calc T = global_simulation->getEnsemble()->T();
+	EqualVelocityAssigner eqVeloAssigner(T);
 
 	double boxMin[3] = {getBoxMin(0), getBoxMin(1), getBoxMin(2)};
 	double boxMax[3] = {getBoxMax(0), getBoxMax(1), getBoxMax(2)};
@@ -141,17 +120,19 @@ unsigned long ParticleCellBase::initCubicGrid(const std::array<unsigned long, 3>
 
 				if (x1In and y1In and z1In) {
 					++numInserted;
-					std::array<vcp_real_calc,3> v = getRandomVelocity<vcp_real_calc>(T, RNG);
+					// Init molecule with zero velocity and use the EqualVelocityAssigner in the next step
 					Molecule dummy(0, &(global_simulation->getEnsemble()->getComponents()->at(0)),
-						x1, y1, z1, v[0], -v[1], v[2]);
+						x1, y1, z1, 0.0, 0.0, 0.0);
+					eqVeloAssigner.assignVelocity(&dummy);
 					buffer.push_back(dummy);
 				}
 
 				if (x2In and y2In and z2In) {
 					++numInserted;
-					std::array<vcp_real_calc,3> v = getRandomVelocity<vcp_real_calc>(T, RNG);
+					// Init molecule with zero velocity and use the EqualVelocityAssigner in the next step
 					Molecule dummy(0, &(global_simulation->getEnsemble()->getComponents()->at(0)),
-						x2, y2, z2, v[0], -v[1], v[2]);
+						x2, y2, z2, 0.0, 0.0, 0.0);
+					eqVeloAssigner.assignVelocity(&dummy);
 					buffer.push_back(dummy);
 				}
 			}

src/particleContainer/ParticleCellBase.h

@@ -16,8 +16,6 @@
 #include <quicksched.h>
 #endif
 
-class Random;
-
 /** @brief ParticleCellBase defines the interface for cells used by the LinkedCells data structure to store molecule data.
  */
 class ParticleCellBase: public Cell {
@@ -83,7 +81,7 @@ class ParticleCellBase: public Cell {
 	virtual void prefetchForForce() const {/*TODO*/}
 
 	unsigned long initCubicGrid(const std::array<unsigned long, 3> &numMoleculesPerDimension,
-								const std::array<double, 3> &simBoxLength, Random &RNG);
+								const std::array<double, 3> &simBoxLength);
 
 //protected: Do not use! use SingleCellIterator instead!
 	// multipurpose:

src/particleContainer/ParticleContainer.h

@@ -235,10 +235,9 @@ class ParticleContainer: public MemoryProfilable {
      * Generates a body-centered cubic grid.
      * @param numMoleculesPerDimension
      * @param simBoxLength
-     * @param seed_offset
      * @return
      */
-	virtual unsigned long initCubicGrid(std::array<unsigned long, 3> numMoleculesPerDimension, std::array<double, 3> simBoxLength, size_t seed_offset) = 0;
+	virtual unsigned long initCubicGrid(std::array<unsigned long, 3> numMoleculesPerDimension, std::array<double, 3> simBoxLength) = 0;
 
 	virtual double* getCellLength() = 0;
 

src/utils/generator/ReplicaFiller.cpp

@@ -116,7 +116,7 @@ class ParticleContainerToBasisWrapper : public ParticleContainer {
 	std::variant<ParticleIterator, SingleCellIterator<ParticleCell>> getMoleculeAtPosition(const double pos[3]) override { return {}; }
 
 	unsigned long initCubicGrid(std::array<unsigned long, 3> numMoleculesPerDimension,
-								std::array<double, 3> simBoxLength, size_t seed_offset) override { return 0; }
+								std::array<double, 3> simBoxLength) override { return 0; }
 
 	size_t getTotalSize() override { return _basis.numMolecules() * sizeof(Molecule); }