From d44e2d1aed1f95c49d9ab6ea53ac290114482e35 Mon Sep 17 00:00:00 2001 From: tchipevn Date: Fri, 12 Apr 2019 19:26:05 +0200 Subject: [PATCH 1/5] change application of mixing rules to be non-symmetric --- src/molecules/Comp2Param.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/molecules/Comp2Param.cpp b/src/molecules/Comp2Param.cpp index c322a2c221..0db592416b 100644 --- a/src/molecules/Comp2Param.cpp +++ b/src/molecules/Comp2Param.cpp @@ -73,6 +73,10 @@ void Comp2Param::initialize( } } ParaStrm& pstrmji = m_ssparatbl(compj, compi); + double xi_other = *mixpos; + ++mixpos; + double eta_other = *mixpos; + ++mixpos; for (unsigned int centerj = 0; centerj < ncj; ++centerj) { const LJcenter& ljcenterj = static_cast(components[compj].ljcenter(centerj)); epsj = ljcenterj.eps(); @@ -81,8 +85,8 @@ void Comp2Param::initialize( const LJcenter& ljcenteri = static_cast(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; From f6547c436035fbff4764bc0e6bf04916beb4ffd9 Mon Sep 17 00:00:00 2001 From: Tchipev Date: Fri, 12 Apr 2019 21:23:34 +0200 Subject: [PATCH 2/5] Fix initialization of mixing rules through XML. --- src/Simulation.cpp | 16 ++-- src/ensemble/EnsembleBase.cpp | 79 +++++++++++++++----- src/ensemble/EnsembleBase.h | 4 +- src/molecules/Comp2Param.cpp | 7 +- src/molecules/mixingrules/LorentzBerthelot.h | 8 ++ src/molecules/mixingrules/MixingRuleBase.cpp | 2 + src/molecules/mixingrules/MixingRuleBase.h | 13 +++- 7 files changed, 96 insertions(+), 33 deletions(-) diff --git a/src/Simulation.cpp b/src/Simulation.cpp index f41e697f4f..d6dc6b7212 100644 --- a/src/Simulation.cpp +++ b/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 */ diff --git a/src/ensemble/EnsembleBase.cpp b/src/ensemble/EnsembleBase.cpp index bdffdd90fb..35a2f68ee7 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,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& 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& 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); + } + // 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; @@ -98,3 +138,4 @@ void Ensemble::setComponentLookUpIDs() { centers += c->numLJcenters(); } } + diff --git a/src/ensemble/EnsembleBase.h b/src/ensemble/EnsembleBase.h index 8fb1bec20c..7c061bcb7d 100644 --- a/src/ensemble/EnsembleBase.h +++ b/src/ensemble/EnsembleBase.h @@ -4,6 +4,7 @@ #include #include +#include #include #include "DomainBase.h" @@ -123,6 +124,7 @@ class Ensemble { protected: + void setVectorOfMixingCoefficientsForComp2Param() const; std::vector _components; std::map _componentnamesToIds; @@ -135,4 +137,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/molecules/Comp2Param.cpp b/src/molecules/Comp2Param.cpp index 0db592416b..fb32ae2f49 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)); @@ -77,6 +77,9 @@ void Comp2Param::initialize( ++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(components[compj].ljcenter(centerj)); epsj = ljcenterj.eps(); 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..e41e8cfd55 100644 --- a/src/molecules/mixingrules/MixingRuleBase.cpp +++ b/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; + + } diff --git a/src/molecules/mixingrules/MixingRuleBase.h b/src/molecules/mixingrules/MixingRuleBase.h index 9445619c08..46716591dc 100644 --- a/src/molecules/mixingrules/MixingRuleBase.h +++ b/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_ */ From ae428956536aef5f39ee9de35fa4fc267172e357 Mon Sep 17 00:00:00 2001 From: Tchipev Date: Fri, 12 Apr 2019 21:36:06 +0200 Subject: [PATCH 3/5] remove modification of mixing coefficients from phasespace readers --- src/Domain.cpp | 1 + src/ensemble/EnsembleBase.cpp | 1 - src/ensemble/EnsembleBase.h | 1 - src/io/ASCIIReader.cpp | 7 +++---- src/io/MPI_IOReader.cpp | 7 +++---- src/io/TcTS.cpp | 13 +------------ 6 files changed, 8 insertions(+), 22 deletions(-) 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/ensemble/EnsembleBase.cpp b/src/ensemble/EnsembleBase.cpp index 35a2f68ee7..d3b4bd13a5 100644 --- a/src/ensemble/EnsembleBase.cpp +++ b/src/ensemble/EnsembleBase.cpp @@ -138,4 +138,3 @@ void Ensemble::setComponentLookUpIDs() { centers += c->numLJcenters(); } } - diff --git a/src/ensemble/EnsembleBase.h b/src/ensemble/EnsembleBase.h index 7c061bcb7d..83d4f7605b 100644 --- a/src/ensemble/EnsembleBase.h +++ b/src/ensemble/EnsembleBase.h @@ -4,7 +4,6 @@ #include #include -#include #include #include "DomainBase.h" 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]; From d5334e4b972bcaea91531cb583a1065adc268db7 Mon Sep 17 00:00:00 2001 From: tchipev Date: Mon, 15 Apr 2019 19:11:03 +0200 Subject: [PATCH 4/5] Roll back possibility for unsymmetric mixing rules. More sanity checks for mixing rules. --- src/Simulation.cpp | 13 ---- src/ensemble/EnsembleBase.cpp | 73 ++++++++++++++++---- src/ensemble/EnsembleBase.h | 5 ++ src/molecules/Comp2Param.cpp | 11 +-- src/molecules/mixingrules/MixingRuleBase.cpp | 16 +++-- src/molecules/mixingrules/MixingRuleBase.h | 3 + 6 files changed, 82 insertions(+), 39 deletions(-) diff --git a/src/Simulation.cpp b/src/Simulation.cpp index d6dc6b7212..3f1f95ec3c 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -240,19 +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 // 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 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 */ if(xmlconfig.changecurrentnode("algorithm")) { /* cutoffs */ diff --git a/src/ensemble/EnsembleBase.cpp b/src/ensemble/EnsembleBase.cpp index d3b4bd13a5..c03e3adb21 100644 --- a/src/ensemble/EnsembleBase.cpp +++ b/src/ensemble/EnsembleBase.cpp @@ -44,6 +44,9 @@ void Ensemble::readXML(XMLfileUnits& xmlconfig) { uint32_t numMixingrules = 0; numMixingrules = query.card(); global_log->info() << "Found " << numMixingrules << " mixing rules." << endl; + 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); @@ -57,18 +60,71 @@ 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); } + checkMixingRules(); setVectorOfMixingCoefficientsForComp2Param(); xmlconfig.changecurrentnode(oldpath); setComponentLookUpIDs(); } +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 << " which is smaller than total 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 << " which is larger than total 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 compID1 = (*m)->getCid1()-1; + unsigned compID2 = (*m)->getCid2()-1; + + if (compID1 == compID2) { + global_log->error() << "Specified rule for same component ids: " << compID1 << " and " << compID2 <<", 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 << " and " << cid2 <<" more than once, which is not allowed." << endl; + global_log->error() << "If you specify " << cid1 << " and " << cid2 <<", then you must not specify " << cid2 << " and " << cid1 << "." << endl; + global_log->error() << "Aborting." << endl; + Simulation::exit(2019); + } + } + } +} + void Ensemble::setVectorOfMixingCoefficientsForComp2Param() const { using std::vector; using std::array; @@ -107,19 +163,10 @@ void Ensemble::setVectorOfMixingCoefficientsForComp2Param() const { // 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); - } + double eta = values[cid1][cid2][0]; + double xi = values[cid1][cid2][1]; + dmixcoeff.push_back(eta); + dmixcoeff.push_back(xi); } } } diff --git a/src/ensemble/EnsembleBase.h b/src/ensemble/EnsembleBase.h index 83d4f7605b..4e241bdc5f 100644 --- a/src/ensemble/EnsembleBase.h +++ b/src/ensemble/EnsembleBase.h @@ -123,6 +123,11 @@ class Ensemble { protected: + /*! 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; diff --git a/src/molecules/Comp2Param.cpp b/src/molecules/Comp2Param.cpp index fb32ae2f49..389632e701 100644 --- a/src/molecules/Comp2Param.cpp +++ b/src/molecules/Comp2Param.cpp @@ -73,13 +73,6 @@ 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(components[compj].ljcenter(centerj)); epsj = ljcenterj.eps(); @@ -88,8 +81,8 @@ void Comp2Param::initialize( const LJcenter& ljcenteri = static_cast(components[compi].ljcenter(centeri)); epsi = ljcenteri.eps(); sigi = ljcenteri.sigma(); - epsilon24 = 24. * xi_other * sqrt(epsi * epsj); - sigma2 = eta_other * .5 * (sigi + sigj); + epsilon24 = 24. * xi * sqrt(epsi * epsj); + sigma2 = eta * .5 * (sigi + sigj); sigma2 *= sigma2; sigperrc2 = sigma2 / (rcLJ * rcLJ); sigperrc6 = sigperrc2 * sigperrc2 * sigperrc2; diff --git a/src/molecules/mixingrules/MixingRuleBase.cpp b/src/molecules/mixingrules/MixingRuleBase.cpp index e41e8cfd55..23a7db7c72 100644 --- a/src/molecules/mixingrules/MixingRuleBase.cpp +++ b/src/molecules/mixingrules/MixingRuleBase.cpp @@ -3,15 +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 46716591dc..4a5c94b7d4 100644 --- a/src/molecules/mixingrules/MixingRuleBase.h +++ b/src/molecules/mixingrules/MixingRuleBase.h @@ -10,6 +10,9 @@ class MixingRuleBase { MixingRuleBase(){}; virtual ~MixingRuleBase(){}; + /** Automatically sets cid1 <= cid2. + * Validity checks are performed elsewhere. + */ virtual void readXML(XMLfileUnits& xmlconfig); unsigned int getCid1() const { From 9cadcfc53bf3d5d4563fee50b2367ddd82e33a1a Mon Sep 17 00:00:00 2001 From: tchipev Date: Tue, 30 Apr 2019 17:04:56 +0200 Subject: [PATCH 5/5] More sanity checks for mixing rules and better error output. --- src/ensemble/EnsembleBase.cpp | 37 ++++++++++++++++++++++++++++------- src/ensemble/EnsembleBase.h | 3 +++ 2 files changed, 33 insertions(+), 7 deletions(-) diff --git a/src/ensemble/EnsembleBase.cpp b/src/ensemble/EnsembleBase.cpp index c03e3adb21..060b5caa4c 100644 --- a/src/ensemble/EnsembleBase.cpp +++ b/src/ensemble/EnsembleBase.cpp @@ -74,6 +74,25 @@ void Ensemble::readXML(XMLfileUnits& xmlconfig) { 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; @@ -84,10 +103,10 @@ void Ensemble::checkMixingRules() const { int exactNumMixingRules = numComponents * (numComponents - 1) / 2; if (numMixingRules < exactNumMixingRules) { - global_log->info() << "Found " << numMixingRules << " which is smaller than total number of mixing rules: " << exactNumMixingRules << "." << endl; + 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 << " which is larger than total number of mixing rules: " << exactNumMixingRules << "." << endl; + 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); } @@ -99,12 +118,16 @@ void Ensemble::checkMixingRules() const { 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 = (*m)->getCid1()-1; - unsigned compID2 = (*m)->getCid2()-1; + unsigned compID1 = outerCID1-1; + unsigned compID2 = outerCID2-1; if (compID1 == compID2) { - global_log->error() << "Specified rule for same component ids: " << compID1 << " and " << compID2 <<", which is not allowed." << endl; + 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); } @@ -116,8 +139,8 @@ void Ensemble::checkMixingRules() const { 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 << " and " << cid2 <<" more than once, which is not allowed." << endl; - global_log->error() << "If you specify " << cid1 << " and " << cid2 <<", then you must not specify " << cid2 << " and " << cid1 << "." << endl; + 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); } diff --git a/src/ensemble/EnsembleBase.h b/src/ensemble/EnsembleBase.h index 4e241bdc5f..f747d54c0e 100644 --- a/src/ensemble/EnsembleBase.h +++ b/src/ensemble/EnsembleBase.h @@ -123,6 +123,9 @@ 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. * */