diff --git a/src/Domain.cpp b/src/Domain.cpp index 62af0bc566..7aa6b39dd3 100644 --- a/src/Domain.cpp +++ b/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::const_iterator pos=_mixcoeff.begin();pos!=_mixcoeff.end();++pos){ checkpointfilestream << *pos; iout++; diff --git a/src/Simulation.cpp b/src/Simulation.cpp index f41e697f4f..3f1f95ec3c 100644 --- a/src/Simulation.cpp +++ b/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 */ diff --git a/src/ensemble/EnsembleBase.cpp b/src/ensemble/EnsembleBase.cpp index bdffdd90fb..060b5caa4c 100644 --- a/src/ensemble/EnsembleBase.cpp +++ b/src/ensemble/EnsembleBase.cpp @@ -8,6 +8,7 @@ #include "Domain.h" #include +#include 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& 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(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> flags(numComponents, vector(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& 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>> values(numComponents, vector>(numComponents,{1.0, 1.0})); + + for(auto m = _mixingrules.begin(); m != _mixingrules.end(); ++m) { + // cast to LB mixing rule + LorentzBerthelotMixingRule & rule = dynamic_cast(**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; diff --git a/src/ensemble/EnsembleBase.h b/src/ensemble/EnsembleBase.h index 8fb1bec20c..f747d54c0e 100644 --- a/src/ensemble/EnsembleBase.h +++ b/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 _components; std::map _componentnamesToIds; @@ -135,4 +144,4 @@ class Ensemble { * and for whatever reason, Domain does not inherit from DomainBase. */ Domain* _simulationDomain; -}; \ No newline at end of file +}; diff --git a/src/io/ASCIIReader.cpp b/src/io/ASCIIReader.cpp index 126320b2d1..e978716764 100644 --- a/src/io/ASCIIReader.cpp +++ b/src/io/ASCIIReader.cpp @@ -207,14 +207,13 @@ void ASCIIReader::readPhaseSpaceHeader(Domain* domain, double timestep) { #endif // Mixing coefficients - vector& 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} diff --git a/src/io/MPI_IOReader.cpp b/src/io/MPI_IOReader.cpp index 2f9f5fefd0..5f29004884 100644 --- a/src/io/MPI_IOReader.cpp +++ b/src/io/MPI_IOReader.cpp @@ -234,14 +234,13 @@ void MPI_IOReader::readPhaseSpaceHeader(Domain* domain, double timestep) { #endif // Mixing coefficients - vector& 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} diff --git a/src/io/TcTS.cpp b/src/io/TcTS.cpp index 7e12344859..58bc27944a 100644 --- a/src/io/TcTS.cpp +++ b/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& 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]; diff --git a/src/molecules/Comp2Param.cpp b/src/molecules/Comp2Param.cpp index c322a2c221..389632e701 100644 --- a/src/molecules/Comp2Param.cpp +++ b/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(components[compi].ljcenter(centeri)); diff --git a/src/molecules/mixingrules/LorentzBerthelot.h b/src/molecules/mixingrules/LorentzBerthelot.h index 1ff88c0014..44775c298b 100644 --- a/src/molecules/mixingrules/LorentzBerthelot.h +++ b/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; diff --git a/src/molecules/mixingrules/MixingRuleBase.cpp b/src/molecules/mixingrules/MixingRuleBase.cpp index cbe82ea51a..23a7db7c72 100644 --- a/src/molecules/mixingrules/MixingRuleBase.cpp +++ b/src/molecules/mixingrules/MixingRuleBase.cpp @@ -3,13 +3,23 @@ #include "utils/Logger.h" #include "utils/xmlfileUnits.h" +#include + 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() } diff --git a/src/molecules/mixingrules/MixingRuleBase.h b/src/molecules/mixingrules/MixingRuleBase.h index 9445619c08..4a5c94b7d4 100644 --- a/src/molecules/mixingrules/MixingRuleBase.h +++ b/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_ */