Tchipev fixes mixing rules

src/Domain.cpp

@@ -531,6 +531,7 @@ void Domain::writeCheckpointHeader(string filename,
 			}
 			unsigned int numperline=_simulation.getEnsemble()->getComponents()->size();
 			unsigned int iout=0;
+			// TODO: this should be removed, but there will be backwards compatibility issues with reading.
 			for(vector<double>::const_iterator pos=_mixcoeff.begin();pos!=_mixcoeff.end();++pos){
 				checkpointfilestream << *pos;
 				iout++;

src/Simulation.cpp

@@ -241,18 +241,16 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 	}
 
 	//The mixing coefficents have to be read in the ensemble part
-	//if this didn't happen fill the mixing coeffs with default values
+	//if this didn't happen fill the mixing coeffs with default values // now filling in happens in Ensemble part.
 	auto& dmixcoeff = _domain->getmixcoeff();
 	const std::size_t compNum = _ensemble->getComponents()->size();
 	//1 Comps: 0 coeffs; 2 Comps: 2 coeffs; 3 Comps: 6 coeffs; 4 Comps 12 coeffs
-	const std::size_t neededCoeffs = compNum*(compNum-1);
-	if(dmixcoeff.size() < neededCoeffs){
-		global_log->warning() << "Not enough mixing coefficients were given! (Filling the missing ones with 1)" << '\n';
-		global_log->warning() << "This can happen because the xml-input doesn't support these yet!" << endl;
-		unsigned int numcomponents = _simulation.getEnsemble()->getComponents()->size();
-		for (unsigned int i = dmixcoeff.size(); i < neededCoeffs; i++) {
-			dmixcoeff.push_back(1);
-		}
+	const std::size_t neededCoeffCombinations = compNum*(compNum-1);
+	const std::size_t actualCoeffCombinations = dmixcoeff.size() / 2; // two entries make one combination.
+	if(actualCoeffCombinations != neededCoeffCombinations){
+		// should never happen now
+		global_log->error() << "Mismatch in mixing coefficients vector size. Mixing coefficients should only be specified through XML." << endl;
+		Simulation::exit(1989);
 	}
 
 	/* algorithm */

src/ensemble/EnsembleBase.cpp

@@ -8,6 +8,7 @@
 #include "Domain.h"
 
 #include <vector>
+#include <array>
 
 using namespace std;
 using Log::global_log;
@@ -43,11 +44,6 @@ void Ensemble::readXML(XMLfileUnits& xmlconfig) {
 	uint32_t numMixingrules = 0;
 	numMixingrules = query.card();
 	global_log->info() << "Found " << numMixingrules << " mixing rules." << endl;
-	_mixingrules.resize(numMixingrules);
-
-	// data structure for mixing coefficients of domain class (still in use!!!)
-	std::vector<double>& dmixcoeff = global_simulation->getDomain()->getmixcoeff();
-	dmixcoeff.clear();
 
 	for(mixingruletIter = query.begin(); mixingruletIter; mixingruletIter++) {
 		xmlconfig.changecurrentnode(mixingruletIter);
@@ -65,25 +61,69 @@ void Ensemble::readXML(XMLfileUnits& xmlconfig) {
 		}
 		mixingrule->readXML(xmlconfig);
 		_mixingrules.push_back(mixingrule);
-
-		/*
-		 * Mixing coefficients
-		 *
-		 * TODO: information of mixing rules (eta, xi) is stored in Domain class and its actually in use
-		 * --> we need to decide where this information should be stored in future, in ensemble class,
-		 * in the way it is done above?
-		 *
-		 */
-		double xi, eta;
-		xmlconfig.getNodeValue("xi", xi);
-		xmlconfig.getNodeValue("eta", eta);
-		dmixcoeff.push_back(xi);
-		dmixcoeff.push_back(eta);
 	}
+
+	setVectorOfMixingCoefficientsForComp2Param();
+
 	xmlconfig.changecurrentnode(oldpath);
 	setComponentLookUpIDs();
 }
 
+void Ensemble::setVectorOfMixingCoefficientsForComp2Param() const {
+	using std::vector;
+	using std::array;
+
+	// data structure for mixing coefficients of domain class (still in use!!!)
+
+	/*
+	 * Mixing coefficients
+	 *
+	 * TODO: information of mixing rules (eta, xi) is stored in Domain class and its actually in use
+	 * --> we need to decide where this information should be stored in future, in ensemble class,
+	 * in the way it is done above?
+	 *
+	 */
+
+	int numComponents = _components.size();
+
+	std::vector<double>& dmixcoeff = global_simulation->getDomain()->getmixcoeff();
+	dmixcoeff.clear();
+
+	// two dimensional vector of Xi and Eta values
+	// NOTE: initialise non-specified ones with default values of 1.0, 1.0
+	vector<vector<array<double, 2>>> values(numComponents, vector<array<double, 2>>(numComponents,{1.0, 1.0}));
+
+	for(auto m = _mixingrules.begin(); m != _mixingrules.end(); ++m) {
+		// cast to LB mixing rule
+		LorentzBerthelotMixingRule & rule = dynamic_cast<LorentzBerthelotMixingRule &>(**m);
+		unsigned compID1 = rule.getCid1()-1;
+		unsigned compID2 = rule.getCid2()-1;
+
+		values[compID1][compID2][0] = rule.getEta();
+		values[compID1][compID2][1] = rule.getXi();
+	}
+
+	// follow the precise way of initialising shit in Comp2Param::initialize
+	// i.e. sort-of symmetrically initialised
+	for (int cid1 = 0; cid1 < numComponents; ++cid1) {
+		for (int cid2 = cid1 + 1; cid2 < numComponents; ++cid2) {
+			{
+				double eta = values[cid1][cid2][0];
+				double xi = values[cid1][cid2][1];
+				dmixcoeff.push_back(eta);
+				dmixcoeff.push_back(xi);
+			}
+			// now push symmetric value, because of how values are read
+			{
+				double eta = values[cid2][cid1][0];
+				double xi = values[cid2][cid1][1];
+				dmixcoeff.push_back(eta);
+				dmixcoeff.push_back(xi);
+			}
+		}
+	}
+}
+
 void Ensemble::setComponentLookUpIDs() {
 	// Get the maximum Component ID.
 	unsigned maxID = 0;

src/ensemble/EnsembleBase.h

@@ -123,6 +123,7 @@ class Ensemble {
 
 protected:
 
+	void setVectorOfMixingCoefficientsForComp2Param() const;
 
 	std::vector<Component> _components;
 	std::map<std::string, int> _componentnamesToIds;
@@ -135,4 +136,4 @@ class Ensemble {
 	 * and for whatever reason, Domain does not inherit from DomainBase.
 	*/
 	Domain* _simulationDomain;
-};
+};

src/io/ASCIIReader.cpp

@@ -207,14 +207,13 @@ void ASCIIReader::readPhaseSpaceHeader(Domain* domain, double timestep) {
 #endif
 
 			// Mixing coefficients
-			vector<double>& dmixcoeff = domain->getmixcoeff();
-			dmixcoeff.clear();
 			for(unsigned int i = 1; i < numcomponents; i++) {
 				for(unsigned int j = i + 1; j <= numcomponents; j++) {
 					double xi, eta;
 					_phaseSpaceHeaderFileStream >> xi >> eta;
-					dmixcoeff.push_back(xi);
-					dmixcoeff.push_back(eta);
+					global_log->warning() << "Reading mixing coefficients from phasespace file no longer supported (was producing wrong values). They will have no effect! Specify them in the .XML file." << endl;
+					//dmixcoeff.push_back(xi);
+					//dmixcoeff.push_back(eta);
 				}
 			}
 			// read in global factor \epsilon_{RF}

src/io/MPI_IOReader.cpp

@@ -234,14 +234,13 @@ void MPI_IOReader::readPhaseSpaceHeader(Domain* domain, double timestep) {
 #endif
 
 			// Mixing coefficients
-			vector<double>& dmixcoeff = domain->getmixcoeff();
-			dmixcoeff.clear();
+			// now only through XML file
 			for(unsigned int i = 1; i < numcomponents; i++) {
 				for(unsigned int j = i + 1; j <= numcomponents; j++) {
 					double xi, eta;
 					_phaseSpaceHeaderFileStream >> xi >> eta;
-					dmixcoeff.push_back(xi);
-					dmixcoeff.push_back(eta);
+					//dmixcoeff.push_back(xi);
+					//dmixcoeff.push_back(eta);
 				}
 			}
 			// read in global factor \epsilon_{RF}

src/io/TcTS.cpp

@@ -41,18 +41,7 @@ void MkTcTSGenerator::readPhaseSpaceHeader(Domain* domain, double /*timestep*/)
 unsigned long
 MkTcTSGenerator::readPhaseSpace(ParticleContainer* particleContainer, Domain* domain, DomainDecompBase*) {
 	// Mixing coefficients
-	vector<double>& dmixcoeff = domain->getmixcoeff();
-	dmixcoeff.clear();
-
-	unsigned int numcomponents = _simulation.getEnsemble()->getComponents()->size();
-	for(unsigned int i = 1; i < numcomponents; i++) {
-		for(unsigned int j = i + 1; j <= numcomponents; j++) {
-			double xi = 1., eta = 1.;
-
-			dmixcoeff.push_back(xi);
-			dmixcoeff.push_back(eta);
-		}
-	}
+	// Mixing coefficients are now only read from the XML file!
 
 	double T = _simulation.getEnsemble()->T();
 	double box[3];

src/molecules/Comp2Param.cpp

@@ -46,9 +46,9 @@ void Comp2Param::initialize(
 			++mixpos;
 			double eta = *mixpos;
 			++mixpos;
-#ifndef NDEBUG
+//#ifndef NDEBUG
 			global_log->info() << "cid+1(compi)=" << compi+1 << " <--> cid+1(compj)=" << compj+1 << ": xi=" << xi << ", eta=" << eta << endl;
-#endif
+//#endif
 			double shift6combined, sigperrc2, sigperrc6;
 			for (unsigned int centeri = 0; centeri < nci; ++centeri) {
 				const LJcenter& ljcenteri = static_cast<const LJcenter&>(components[compi].ljcenter(centeri));
@@ -73,6 +73,13 @@ void Comp2Param::initialize(
 				}
 			}
 			ParaStrm& pstrmji = m_ssparatbl(compj, compi);
+			double xi_other = *mixpos;
+			++mixpos;
+			double eta_other = *mixpos;
+			++mixpos;
+//#ifndef NDEBUG
+			global_log->info() << "cid+1(compj)=" << compj+1 << " <--> cid+1(compi)=" << compi+1 << ": xi=" << xi_other << ", eta=" << eta_other << endl;
+//#endif
 			for (unsigned int centerj = 0; centerj < ncj; ++centerj) {
 				const LJcenter& ljcenterj = static_cast<const LJcenter&>(components[compj].ljcenter(centerj));
 				epsj = ljcenterj.eps();
@@ -81,8 +88,8 @@ void Comp2Param::initialize(
 					const LJcenter& ljcenteri = static_cast<const LJcenter&>(components[compi].ljcenter(centeri));
 					epsi = ljcenteri.eps();
 					sigi = ljcenteri.sigma();
-					epsilon24 = 24. * xi * sqrt(epsi * epsj);
-					sigma2 = eta * .5 * (sigi + sigj);
+					epsilon24 = 24. * xi_other * sqrt(epsi * epsj);
+					sigma2 = eta_other * .5 * (sigi + sigj);
 					sigma2 *= sigma2;
 					sigperrc2 = sigma2 / (rcLJ * rcLJ);
 					sigperrc6 = sigperrc2 * sigperrc2 * sigperrc2;

src/molecules/mixingrules/LorentzBerthelot.h

@@ -11,6 +11,14 @@ class LorentzBerthelotMixingRule : public MixingRuleBase {
 	virtual ~LorentzBerthelotMixingRule() {}
 	void readXML(XMLfileUnits& xmlconfig);
 
+	double getEta() const {
+		return _eta;
+	}
+
+	double getXi() const {
+		return _xi;
+	}
+
 private:
 	double _eta;
 	double _xi;

src/molecules/mixingrules/MixingRuleBase.cpp

@@ -11,5 +11,7 @@ void MixingRuleBase::readXML(XMLfileUnits& xmlconfig) {
 	xmlconfig.getNodeValue("@cid2", _cid2);
 	global_log->info() << "Component id1: " << _cid1 << endl;
 	global_log->info() << "Component id2: " << _cid2 << endl;
+
+
 }
 

src/molecules/mixingrules/MixingRuleBase.h

@@ -11,9 +11,18 @@ class MixingRuleBase {
 	virtual ~MixingRuleBase(){};
 
 	virtual void readXML(XMLfileUnits& xmlconfig);
+
+	unsigned int getCid1() const {
+		return _cid1;
+	}
+
+	unsigned int getCid2() const {
+		return _cid2;
+	}
+
 private:
-	std::string _cid1;
-	std::string _cid2;
+	unsigned int _cid1;
+	unsigned int _cid2;
 };
 
 #endif /* MIXINGRULE_BASE_H_ */