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

@@ -240,21 +240,6 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 		Simulation::exit(1);
 	}
 
-	//The mixing coefficents have to be read in the ensemble part
-	//if this didn't happen fill the mixing coeffs with default values
-	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);
-		}
-	}
-
 	/* algorithm */
 	if(xmlconfig.changecurrentnode("algorithm")) {
 		/* cutoffs */

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,9 @@ 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();
+	if(numMixingrules > 0) {
+		global_log->info() << "Mixing rules are applied symmetrically: specifying a rule for AB also sets BA." << endl;
+	}
 
 	for(mixingruletIter = query.begin(); mixingruletIter; mixingruletIter++) {
 		xmlconfig.changecurrentnode(mixingruletIter);
@@ -61,29 +60,140 @@ void Ensemble::readXML(XMLfileUnits& xmlconfig) {
 
 		} else {
 			global_log->error() << "Unknown mixing rule " << mixingruletype << endl;
+			global_log->error() << "Only rule type LB supported." << endl;
 			Simulation::exit(1);
 		}
 		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);
 	}
+
+	checkMixingRules();
+	setVectorOfMixingCoefficientsForComp2Param();
+
 	xmlconfig.changecurrentnode(oldpath);
 	setComponentLookUpIDs();
 }
 
+
+void Ensemble::checkMixingRuleCID(unsigned cid) const {
+	int internalCID = static_cast<int>(cid) - 1;
+	if (internalCID < 0) {
+		global_log->error() << "Found mixing rule for component ID " << cid << ", which is not allowed." << endl;
+		global_log->error() << "Component IDs start from 1." << endl;
+		global_log->error() << "Aborting." << endl;
+		Simulation::exit(2019);
+	}
+
+	int numComponents = _components.size();
+	if (internalCID >= numComponents) {
+		global_log->error() << "Found mixing rule for component ID " << cid << ", which is not allowed." << endl;
+		global_log->error() << "Maximal component ID: " << numComponents << "." << endl;
+		global_log->error() << "Aborting." << endl;
+		Simulation::exit(2019);
+	}
+}
+
+void Ensemble::checkMixingRules() const {
+	using std::vector;
+
+	int numMixingRules = _mixingrules.size();
+	int numComponents = _components.size();
+
+	//1 Comps: 0 coeffs; 2 Comps: 1 coeffs; 3 Comps: 3 coeffs; 4 Comps 6 coeffs : C * (C-1) / 2
+	int exactNumMixingRules = numComponents * (numComponents - 1) / 2;
+
+	if (numMixingRules < exactNumMixingRules) {
+		global_log->info() << "Found " << numMixingRules << " mixing rules, which is smaller than exact number of mixing rules: " << exactNumMixingRules << "." << endl;
+		global_log->info() << "The remaining mixing rules are populated with default values of (xi, eta) = (1.0, 1.0)." << endl;
+	} else if (numMixingRules > exactNumMixingRules) {
+		global_log->error() << "Found " << numMixingRules << " mixing rules, which is larger than exact number of mixing rules: " << exactNumMixingRules << "." << endl;
+		global_log->error() << "Aborting." << endl;
+		Simulation::exit(2019);
+	}
+
+
+	// check that:
+	// II doesn't appear,
+	// if IJ appears, then JI doesn't appear.
+	vector<vector<int>> flags(numComponents, vector<int>(numComponents, 0));
+
+	for(auto m = _mixingrules.begin(); m != _mixingrules.end(); ++m) {
+		unsigned outerCID1 = (*m)->getCid1();
+		unsigned outerCID2 = (*m)->getCid2();
+		checkMixingRuleCID(outerCID1);
+		checkMixingRuleCID(outerCID2);
+
+		unsigned compID1 = outerCID1-1;
+		unsigned compID2 = outerCID2-1;
+
+		if (compID1 == compID2) {
+			global_log->error() << "Specified rule for same component ids: " << compID1+1 << " and " << compID2+1 <<", which is not allowed." << endl;
+			global_log->error() << "Aborting." << endl;
+			Simulation::exit(2019);
+		}
+
+		flags[compID1][compID2] ++;
+	}
+
+	for (int cid1 = 0; cid1 < numComponents; ++cid1) {
+		for (int cid2 = cid1 + 1; cid2 < numComponents; ++cid2) {
+			int flag = flags[cid1][cid2];
+			if (flag > 1) {
+				global_log->error() << "Specified rule for component ids: " << cid1+1 << " and " << cid2+1 <<" more than once, which is not allowed." << endl;
+				global_log->error() << "If you specify " << cid1+1 << " and " << cid2+1 <<", then you must not specify " << cid2+1 << " and " << cid1+1 << "." << endl;
+				global_log->error() << "Aborting." << endl;
+				Simulation::exit(2019);
+			}
+		}
+	}
+}
+
+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);
+		}
+	}
+}
+
 void Ensemble::setComponentLookUpIDs() {
 	// Get the maximum Component ID.
 	unsigned maxID = 0;

src/ensemble/EnsembleBase.h

@@ -123,6 +123,15 @@ class Ensemble {
 
 protected:
 
+	/*! Check that CID's are between 1 and the number of components */
+	void checkMixingRuleCID(unsigned cid) const;
+
+	/*! We don't support unsymmetric mixing rules.
+	 * Check that I-J rule appears only once and no II combination appears.
+	 * */
+	void checkMixingRules() const;
+
+	void setVectorOfMixingCoefficientsForComp2Param() const;
 
 	std::vector<Component> _components;
 	std::map<std::string, int> _componentnamesToIds;
@@ -135,4 +144,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));

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

@@ -3,13 +3,23 @@
 #include "utils/Logger.h"
 #include "utils/xmlfileUnits.h"
 
+#include <algorithm>
+
 using namespace std;
 using Log::global_log;
 
 void MixingRuleBase::readXML(XMLfileUnits& xmlconfig) {
-	xmlconfig.getNodeValue("@cid1", _cid1);
-	xmlconfig.getNodeValue("@cid2", _cid2);
-	global_log->info() << "Component id1: " << _cid1 << endl;
-	global_log->info() << "Component id2: " << _cid2 << endl;
+	unsigned int readFirst, readSecond;
+
+	xmlconfig.getNodeValue("@cid1", readFirst);
+	xmlconfig.getNodeValue("@cid2", readSecond);
+
+	global_log->info() << "Component id1: " << readFirst << endl;
+	global_log->info() << "Component id2: " << readSecond << endl;
+
+	_cid1 = std::min(readFirst, readSecond);
+	_cid2 = std::max(readFirst, readSecond);
+
+	// sanity checked in Ensemble::checkMixingRules()
 }
 

src/molecules/mixingrules/MixingRuleBase.h

@@ -10,10 +10,22 @@ class MixingRuleBase {
 	MixingRuleBase(){};
 	virtual ~MixingRuleBase(){};
 
+	/** Automatically sets cid1 <= cid2.
+	 * Validity checks are performed elsewhere.
+	 */
 	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_ */