diff --git a/AUTHORS b/AUTHORS index c7bc327..bb770ec 100644 --- a/AUTHORS +++ b/AUTHORS @@ -12,6 +12,10 @@ Hauke Rabbel, Germany Martin Wengenmayr, Germany Marco Werner, Germany +Other contributors: +-------------------- +Cornelia StrĂ¼big, Germany: FeatureWall + These authors keep copyright of their contributions. They just grant a license to everyone to use it as detailed in LICENSE.txt. diff --git a/CMakeLists.txt b/CMakeLists.txt index 6b5fecb..11c092e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,9 +32,9 @@ PROJECT (LeMonADE) SET (APPLICATION_NAME "LeMonADE") SET (APPLICATION_CODENAME "${PROJECT_NAME}") SET (APPLICATION_COPYRIGHT_YEARS "2015") -SET (APPLICATION_VERSION_MAJOR 1) -SET (APPLICATION_VERSION_MINOR 1) -SET (APPLICATION_VERSION_PATCH 3) +SET (APPLICATION_VERSION_MAJOR 2) +SET (APPLICATION_VERSION_MINOR 0) +SET (APPLICATION_VERSION_PATCH 0) SET (APPLICATION_VERSION_TYPE SNAPSHOT) SET (APPLICATION_VERSION_STRING "${APPLICATION_VERSION_MAJOR}.${APPLICATION_VERSION_MINOR}.${APPLICATION_VERSION_PATCH}-${APPLICATION_VERSION_TYPE}") SET (APPLICATION_ID "${APPLICATION_VENDOR_ID}.${PROJECT_NAME}") diff --git a/include/LeMonADE/feature/FeatureMoleculesIO.h b/include/LeMonADE/feature/FeatureMoleculesIO.h index 8ffa6d0..2ecc667 100644 --- a/include/LeMonADE/feature/FeatureMoleculesIO.h +++ b/include/LeMonADE/feature/FeatureMoleculesIO.h @@ -28,6 +28,8 @@ along with LeMonADE. If not, see . #ifndef LEMONADE_FEATURE_FEATUREMOLECULESIO_H #define LEMONADE_FEATURE_FEATUREMOLECULESIO_H +#include + #include #include #include diff --git a/include/LeMonADE/feature/FeatureNNInteractionBcc.h b/include/LeMonADE/feature/FeatureNNInteractionBcc.h new file mode 100644 index 0000000..3c5b118 --- /dev/null +++ b/include/LeMonADE/feature/FeatureNNInteractionBcc.h @@ -0,0 +1,602 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ + + +#ifndef FEATURE_NN_INTERACTION_BCC_H +#define FEATURE_NN_INTERACTION_BCC_H + +/** + * @file + * @date 2016/06/18 + * @author Hauke Rabbel + * @brief Definition and implementation of class template FeatureNNInteractionBcc + * + * @todo MoveAddMonomerBcc is used here, which might be obsolete. +**/ + +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * @class FeatureNNInteractionBcc + * @brief Provides interaction of monomers on distances d<=sqrt(12)a for bccBFM + * + * @tparam FeatureLatticeType Underlying lattice feature, e.g. FeatureLattice or + * FeatureLatticePowerOfTwo (template template parameter) + * + * @details + * The interaction energy can be set for pairs of monomer-types A,B, where + * the type can be any integer between 1 and 255. + * The feature automatically adds FeatureExcludedVolumeBcc + * to the system. Given an energy E in kT between two types, the interaction potential + * as a function of the distance d is: + * - inf d < sqrt(3) (implicitly through excluded volume) + * - E sqrt(3)< d <= sqrt(12) + * - 0 d>sqrt(12) + * . + * Usage: In the feature list defining Ingredients use this feature as + * FeatureNNInteractionBcc (arbitrary lattices), or as + * FeatureNNInteractionBcc (2**n lattices) + * The feature adds the bfm-file command !nn_interaction A B E + * for monomers of types A B with interaction energy of E in kT. +**/ + +template class FeatureLatticeType> +class FeatureNNInteractionBcc:public Feature +{ + +private: + //! Type for the underlying lattice, used as template parameter for FeatureLatticeType<...> + typedef uint8_t lattice_value_type; + + //! Interaction energies between monomer types. Max. type=255 given by max(uint8_t)=255 + double interactionTable[256][256]; + + //! Lookup table for exp(-interactionTable[a][b]) + double probabilityLookup[256][256]; + + //! Returns this feature's factor for the acceptance probability for the given Monte Carlo move + template + double calculateAcceptanceProbability(const IngredientsType& ingredients, + const MoveLocalBcc& move) const; + + //! Occupies the lattice with the attribute tags of all monomers + template + void fillLattice(IngredientsType& ingredients); + + //! Access to array probabilityLookup with extra checks in Debug mode + double getProbabilityFactor(int32_t typeA,int32_t typeB) const; + + +public: + + FeatureNNInteractionBcc(); + ~FeatureNNInteractionBcc(){} + + + //This feature adds interaction energies, so it requires FeatureBoltzmann + typedef LOKI_TYPELIST_1(FeatureBoltzmann) required_features_back; + + //FeatureExcludedVolumeSc needs to be in front, because FeatureNNInteractionBcc + //re-initializes the lattice and overwrites what FeatureExcludedVolumeSc has written. + //FeatureAttributes needs to be in front, because when a monomer is added to the system + //by a MoveAddScMonomer, its attribute has to be set before it is written to the lattice. + typedef LOKI_TYPELIST_2( + FeatureAttributes, + FeatureExcludedVolumeBcc >) + required_features_front; + + + //! check for all Monte Carlo moves without special check functions (always true) + template + bool checkMove(const IngredientsType& ingredients,const MoveBase& move) const; + + //! check move for bcc-BFM local move. always throws std::runtime_error + template + bool checkMove(const IngredientsType& ingredients,const MoveLocalSc& move) const; + + //! check move for bcc-BFM local move + template + bool checkMove(const IngredientsType& ingredients,MoveLocalBcc& move) const; + + //! apply function for all Monte Carlo moves without special apply functions (does nothing) + template + void applyMove(const IngredientsType& ing, const MoveBase& move){} + + //! apply function for sc-BFM local move (always throws std::runtime_error) + template + void applyMove(const IngredientsType& ing, const MoveLocalSc& move); + + //! apply function for adding a monomer in bcc-BFM + template + void applyMove(IngredientsType& ing, const MoveAddMonomerBcc& move); + + //note: apply function for bcc-BFM local move is not necessary, because + //job of moving lattice entries is done by the underlying FeatureLatticeType + + //! guarantees that the lattice is properly occupied with monomer attributes + template + void synchronize(IngredientsType& ingredients); + + //!adds interaction energy between two types of monomers + void setNNInteraction(int32_t typeA,int32_t typeB,double energy); + + //!returns the interaction energy between two types of monomers + double getNNInteraction(int32_t typeA,int32_t typeB) const; + + //!export bfm-file read command !nn_interaction + template + void exportRead(FileImport & fileReader); + + //!export bfm-file write command !nn_interaction + template + void exportWrite(AnalyzerWriteBfmFile & fileWriter) const; + +}; + + +////////////////// IMPLEMENTATION OF MEMBERS ////////////////////////////////// + +/** + * @brief Constructor + **/ +template class LatticeClassType> +FeatureNNInteractionBcc::FeatureNNInteractionBcc() +{ + //initialize the energy and probability lookups with default values + for(size_t n=0;n<256;n++) + { + for(size_t m=0;m<256;m++) + { + interactionTable[m][n]=0.0; + probabilityLookup[m][n]=1.0; + } + } +} + +/** + * @details Returns true for all moves other than the ones that have specialized versions of this function. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move other than moves with specialized functions. + * @return true (always) + **/ +template class LatticeClassType> +template +bool FeatureNNInteractionBcc::checkMove(const IngredientsType& ingredients, + const MoveBase& move) const +{ + return true; +} + +/** + * @details calculates the factor for the acceptance probability of the move + * arising from the contact interactions and adds it to the move. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move of type MoveLocalSc + * @return true (always) + **/ +template class LatticeClassType> +template +bool FeatureNNInteractionBcc::checkMove(const IngredientsType& ingredients, + MoveLocalBcc& move) const +{ + //add the probability factor coming from this feature, then return true, + //because the total probability is evaluated by FeatureBoltzmann at the end + double prob=calculateAcceptanceProbability(ingredients,move); + move.multiplyProbability(prob); + return true; +} + +/** + * @details Because moves of type MoveLocalSc must not be used with this + * feature, this function always throws an exception when called. The function + * is only implemented for savety purposes. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move of type MoveLocalBcc + * @throw std::runtime_error + * @return false always throws exception before returning + **/ +template class LatticeClassType> +template +bool FeatureNNInteractionBcc::checkMove(const IngredientsType& ingredients, + const MoveLocalSc& move) const +{ + //throw exception in case someone accidentaly uses a bcc-BFM move with this feature + std::stringstream errormessage; + errormessage<<"FeatureNNInteractionBcc::checkMove(...):\n"; + errormessage<<"attempting to use sc-BFM move, which is not allowed\n"; + throw std::runtime_error(errormessage.str()); + + return false; +} + + +/** + * @details When a new monomer is added to the system, the lattice site occupied + * by this monomer must be filled with the attribute tag. This function takes care + * of this. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move of type MoveAddMonomerBcc + * @throw std::runtime_error if attribute tag is not in range [1,255] + **/ +template class LatticeClassType> +template +void FeatureNNInteractionBcc::applyMove(IngredientsType& ing, + const MoveAddMonomerBcc& move) +{ + //get the position and attribute tag of the monomer to be inserted + VectorInt3 pos=move.getPosition(); + lattice_value_type type=lattice_value_type(move.getTag()); + + //the feature is based on a uint8_t lattice, thus the max type must not + //exceed the max value of uint8_t (255) + if(int32_t(type) != move.getTag() || int32_t(type)==0) + { + std::stringstream errormessage; + errormessage<<"FeatureNNInteractionBcc::applyMove(MoveAddMonomerBcc)\n"; + errormessage<<"Trying to add monomer with type "<maxType=255\n"; + throw std::runtime_error(errormessage.str()); + } + + //update lattice + ing.setLatticeEntry(pos,type); +} + +/** + * @details Because moves of type MoveLocalSc must not be used with this + * feature, this function always throws an exception when called. The function + * is only implemented for savety purposes. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move of type MoveLocalSc + * @throw std::runtime_error + **/ +template class LatticeClassType> +template +void FeatureNNInteractionBcc::applyMove(const IngredientsType& ing, + const MoveLocalSc& move) +{ + //throw exception in case someone accidentaly uses a sc-BFM move with this feature + std::stringstream errormessage; + errormessage<<"FeatureNNInteractionBcc::applyMove(...):\n"; + errormessage<<"attempting to use sc-BFM move, which is not allowed\n"; + throw std::runtime_error(errormessage.str()); + +} + +/** + * @tparam IngredientsType The type of the system including all features + * @param [in] ingredients A reference to the IngredientsType - mainly the system + **/ +template class LatticeClassType> +template +void FeatureNNInteractionBcc::synchronize(IngredientsType& ingredients) +{ + + //refill the lattice with attribute tags + //caution: this overwrites, what is currently written on the lattice + fillLattice(ingredients); + +} + + +/** + * @details If not compiled with DEBUG flag this function only returns the content + * of the lookup table probabilityLookup. If compiled with DEBUG flag it checks + * that the attribute tags typeA, typeB are within the allowed range. + * @param typeA monomer attribute type in range [1,255] + * @param typeB monomer attribute type in range [1,255] + * @throw std::runtime_error In debug mode, if types are not in range [1,255] + **/ +template class LatticeClassType> +inline double FeatureNNInteractionBcc::getProbabilityFactor(int32_t typeA, + int32_t typeB) const +{ +#ifdef DEBUG + //extra checks only in debug mode, because this is very frequently called + //and this costs performance + if(typeA<0 || typeA>255 || typeB<0 || typeB>255){ + std::stringstream errormessage; + errormessage<<"***FeatureNNInteractionBcc::getInteraction(typeA,typeB)***\n"; + errormessage<<"probability undefined between types "< class LatticeClassType> +template +void FeatureNNInteractionBcc::exportRead(FileImport< IngredientsType >& fileReader) +{ + typedef FeatureNNInteractionBcc my_type; + fileReader.registerRead("!nn_interaction",new ReadNNInteraction(*this)); +} + + +/** + * The function is called by the Ingredients class when an object of type Ingredients + * is associated with an object of type AnalyzerWriteBfmFile. The export of the Writes is thus + * taken care automatically when it becomes necessary.\n + * Registered Write-Out Commands: + * - !nn_interaction + * + * @tparam IngredientsType The type of the system including all features + * @param fileWriter File writer for the bfm-file. + */ +template class LatticeClassType> +template +void FeatureNNInteractionBcc::exportWrite(AnalyzerWriteBfmFile< IngredientsType >& fileWriter) const +{ + typedef FeatureNNInteractionBcc my_type; + fileWriter.registerWrite("!nn_interaction",new WriteNNInteraction(*this)); +} + +/** + * @details occupies the lattice with the attribute tags of the monomers + * as this is required to determine the contact interactions in this feature. + * An additional check is performed asserting that the tags are in the range [1,255] + * + * @tparam IngredientsType The type of the system including all features + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @throw std::runtime_error In case a monomer has attribute tag not in [1,255] + **/ +template class LatticeClassType> +template +void FeatureNNInteractionBcc::fillLattice(IngredientsType& ingredients) +{ + const typename IngredientsType::molecules_type& molecules=ingredients.getMolecules(); + + for(size_t n=0;n class LatticeClassType> +template +double FeatureNNInteractionBcc::calculateAcceptanceProbability( + const IngredientsType& ingredients, + const MoveLocalBcc& move) const +{ + + VectorInt3 oldPos=ingredients.getMolecules()[move.getIndex()]; + VectorInt3 direction=move.getDir(); + int32_t monoType=ingredients.getMolecules()[move.getIndex()].getAttributeTag(); + + //get three directions that define which lattice sites have to be checked + //these directions are v1=(2*deltaX,0,0), v2=(0,2*deltaY,0), v3=(0,0,2*deltaZ) + //the idea is that the relative positions v1,v1+v2,v1+v3,... do not have to + //be checked here, as they cannot be occupied due to excluded volume restrictions + VectorInt3 v1(2*direction.getX(),0,0); + VectorInt3 v2(0,2*direction.getY(),0); + VectorInt3 v3(0,0,2*direction.getZ()); + + + //the probability is calculated by going through all possible lattice sites + //at which the contacts may have changed. At every site the type of the + //monomer sitting there is retrieved from the lattice. the additional + //factor for the probability (exp(-deltaE/kT)) is retrieved from the + //lookup using getProbabilityFactor. For new contacts this factor is multiplied + //with the probability, for contacts taken away the probability is devided. + VectorInt3 actual=oldPos; + double prob_div=1.0; + + actual-=v2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + + //now check back side (contacts taken away) + //again, not the complete contact shell has to be checked. the sites + //-v1,-v1-v2,-v1-v3,... relative to the new position cannot be occupied + //because they were previously blocked by the excluded volume of the monomer + //to be moved. + double prob=1.0; + actual=oldPos+direction; + + actual+=v2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=v1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=v3; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + + prob/=prob_div; + return prob; + +} + +/** + * @param typeA monomer attribute tag in range [1,255] + * @param typeB monomer attribute tag in range [1,255] + * @param interaction energy between typeA and typeB + * @throw std::runtime_error In case typeA or typeB exceed range [1,255] + **/ +template class LatticeClassType> +void FeatureNNInteractionBcc::setNNInteraction(int32_t typeA, + int32_t typeB, + double energy) +{ + if(0 class LatticeClassType> +double FeatureNNInteractionBcc::getNNInteraction(int32_t typeA, + int32_t typeB) const +{ + + if(0. + +--------------------------------------------------------------------------------*/ + +#ifndef FEATURE_NN_INTERACTION_READ_WRITE_H +#define FEATURE_NN_INTERACTION_READ_WRITE_H + +/** + * @file + * @date 2016/06/18 + * @author Hauke Rabbel + * @brief Def. and impl. of class templates ReadNNInteraction and WriteNNInteraction +**/ + +#include +#include +#include + +/** + * @class ReadNNInteraction + * @brief Handles BFM-file read command !nn_interaction + * @tparam IngredientsType Ingredients class storing all system information. +**/ +template < class IngredientsType> +class ReadNNInteraction: public ReadToDestination +{ +public: + ReadNNInteraction(IngredientsType& i):ReadToDestination(i){} + virtual ~ReadNNInteraction(){} + virtual void execute(); +}; + +/** + * @class WriteNNInteraction + * @brief Handles BFM-file write command !nn_interaction + * @tparam IngredientsType Ingredients class storing all system information. +**/ +template +class WriteNNInteraction:public AbstractWrite +{ +public: + + //constructor sets the headerOnly tag, such that the interaction + //is written only once at the beginning of the output file. + WriteNNInteraction(const IngredientsType& i) + :AbstractWrite(i){this->setHeaderOnly(true);} + + virtual ~WriteNNInteraction(){} + + virtual void writeStream(std::ostream& strm); +}; + +/////////////MEMBER IMPLEMENTATIONS //////////////////////////////////////////// + +/** + * @brief Executes the reading routine to extract \b !nn_interaction. + * + * @throw fail to read monomer types or interaction energy. + **/ +template +void ReadNNInteraction::execute() +{ + IngredientsType& ingredients=this->getDestination(); + std::istream& file=this->getInputStream(); + + int32_t typeA,typeB; + double interactionConstant; + + //set stream to throw exception on fail + file.exceptions(file.exceptions() | std::ifstream::failbit); + + try + { + file>>typeA; + file>>typeB; + file>>interactionConstant; + } + catch(std::ifstream::failure e) + { + std::stringstream errormessage; + errormessage<<"ReadNNInteraction::execute().\n"; + errormessage<<"Could not read interaction from file\n"; + errormessage<<"Previous error: "< +void WriteNNInteraction::writeStream(std::ostream& stream) +{ + int32_t nSpecies=255; //number must fit into 8 bit (uint8_t based lattice) + stream<<"## nearest neighbor interactions between types in kT (default 0.0kT)\n"; + + for(int32_t typeA=1;typeA<=nSpecies;typeA++) + { + for(int32_t typeB=1;typeB<=typeA;typeB++) + { + if(this->getSource().getNNInteraction(typeA,typeB)!=0.0) + { + stream<<"!nn_interaction "<getSource().getNNInteraction(typeB,typeA)<<"\n"; + } + + } + + } + stream<<"\n\n"; + +} + + + + +#endif // FEATURE_NN_INTERACTION_READ_WRITE_H diff --git a/include/LeMonADE/feature/FeatureNNInteractionSc.h b/include/LeMonADE/feature/FeatureNNInteractionSc.h new file mode 100644 index 0000000..2bfcecc --- /dev/null +++ b/include/LeMonADE/feature/FeatureNNInteractionSc.h @@ -0,0 +1,599 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ + + +#ifndef FEATURE_NN_INTERACTION_H +#define FEATURE_NN_INTERACTION_H + +/** + * @file + * @date 2016/06/18 + * @author Hauke Rabbel + * @brief Definition and implementation of class template FeatureNNInteractionSc +**/ + +#include +#include +#include +#include +#include +#include +#include + +/** + * @class FeatureNNInteractionSc + * @brief Provides interaction of monomers on distances d<=sqrt(6) for standard BFM + * + * @tparam FeatureLatticeType Underlying lattice feature, e.g. FeatureLattice or + * FeatureLatticePowerOfTwo (template template parameter) + * + * @details + * The interaction energy can be set for pairs of monomer-types A,B, where + * the type can be any integer between 1 and 255. + * The feature automatically adds FeatureExcludedVolumeSc + * to the system. Given an energy E in kT between two types, the interaction potential + * as a function of the distance d is: + * - inf d<2 (implicitly through excluded volume) + * - 4*E d=2 + * - 2*E d=sqrt(5) + * - 1*E d=sqrt(6) + * - 0 d>sqrt(6) + * . + * Usage: In the feature list defining Ingredients use this feature as + * FeatureNNInteractionSc (arbitrary lattices), or as + * FeatureNNInteractionSc (2**n lattices) + * The feature adds the bfm-file command !nn_interaction A B E + * for monomers of types A B with interaction energy of E in kT. +**/ + +template class FeatureLatticeType> +class FeatureNNInteractionSc:public Feature +{ + +private: + //! Type for the underlying lattice, used as template parameter for FeatureLatticeType<...> + typedef uint8_t lattice_value_type; + + //! Interaction energies between monomer types. Max. type=255 given by max(uint8_t)=255 + double interactionTable[256][256]; + + //! Lookup table for exp(-interactionTable[a][b]) + double probabilityLookup[256][256]; + + //! Returns this feature's factor for the acceptance probability for the given Monte Carlo move + template + double calculateAcceptanceProbability(const IngredientsType& ingredients, + const MoveLocalSc& move) const; + + //! Occupies the lattice with the attribute tags of all monomers + template + void fillLattice(IngredientsType& ingredients); + + //! Access to array probabilityLookup with extra checks in Debug mode + double getProbabilityFactor(int32_t typeA,int32_t typeB) const; + + +public: + + FeatureNNInteractionSc(); + ~FeatureNNInteractionSc(){} + + + //This feature adds interaction energies, so it requires FeatureBoltzmann + typedef LOKI_TYPELIST_1(FeatureBoltzmann) required_features_back; + + //FeatureExcludedVolumeSc needs to be in front, because FeatureNNInteractionSc + //re-initializes the lattice and overwrites what FeatureExcludedVolumeSc has written. + //FeatureAttributes needs to be in front, because when a monomer is added to the system + //by a MoveAddMonomerSc, its attribute has to be set before it is written to the lattice. + typedef LOKI_TYPELIST_2( + FeatureAttributes, + FeatureExcludedVolumeSc >) + required_features_front; + + + //! check for all Monte Carlo moves without special check functions (always true) + template + bool checkMove(const IngredientsType& ingredients,const MoveBase& move) const; + + //! check for standard sc-BFM local move + template + bool checkMove(const IngredientsType& ingredients,MoveLocalSc& move) const; + + //! check move for bcc-BFM local move. always throws std::runtime_error + template + bool checkMove(const IngredientsType& ingredients,const MoveLocalBcc& move) const; + + //! apply function for all Monte Carlo moves without special apply functions (does nothing) + template + void applyMove(const IngredientsType& ing, const MoveBase& move){} + + //! apply function for bcc-BFM local move (always throws std::runtime_error) + template + void applyMove(const IngredientsType& ing, const MoveLocalBcc& move); + + //! apply function for adding a monomer in sc-BFM + template + void applyMove(IngredientsType& ing, const MoveAddMonomerSc& move); + + //note: apply function for sc-BFM local move is not necessary, because + //job of moving lattice entries is done by the underlying FeatureLatticeType + + //! guarantees that the lattice is properly occupied with monomer attributes + template + void synchronize(IngredientsType& ingredients); + + //!adds interaction energy between two types of monomers + void setNNInteraction(int32_t typeA,int32_t typeB,double energy); + + //!returns the interaction energy between two types of monomers + double getNNInteraction(int32_t typeA,int32_t typeB) const; + + //!export bfm-file read command !nn_interaction + template + void exportRead(FileImport & fileReader); + + //!export bfm-file write command !nn_interaction + template + void exportWrite(AnalyzerWriteBfmFile & fileWriter) const; + +}; + + +////////////////// IMPLEMENTATION OF MEMBERS ////////////////////////////////// + +/** + * @brief Constructor + **/ +template class LatticeClassType> +FeatureNNInteractionSc::FeatureNNInteractionSc() +{ + //initialize the energy and probability lookups with default values + for(size_t n=0;n<256;n++) + { + for(size_t m=0;m<256;m++) + { + interactionTable[m][n]=0.0; + probabilityLookup[m][n]=1.0; + } + } +} + +/** + * @details Returns true for all moves other than the ones that have specialized versions of this function. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move other than moves with specialized functions. + * @return true (always) + **/ +template class LatticeClassType> +template +bool FeatureNNInteractionSc::checkMove(const IngredientsType& ingredients, + const MoveBase& move) const +{ + return true; +} + +/** + * @details calculates the factor for the acceptance probability of the move + * arising from the contact interactions and adds it to the move. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move of type MoveLocalSc + * @return true (always) + **/ +template class LatticeClassType> +template +bool FeatureNNInteractionSc::checkMove(const IngredientsType& ingredients, + MoveLocalSc& move) const +{ + //add the probability factor coming from this feature, then return true, + //because the total probability is evaluated by FeatureBoltzmann at the end + double prob=calculateAcceptanceProbability(ingredients,move); + move.multiplyProbability(prob); + return true; +} + +/** + * @details Because moves of type MoveLocalBcc must not be used with this + * feature, this function always throws an exception when called. The function + * is only implemented for savety purposes. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move of type MoveLocalBcc + * @throw std::runtime_error + * @return false always throws exception before returning + **/ +template class LatticeClassType> +template +bool FeatureNNInteractionSc::checkMove(const IngredientsType& ingredients, + const MoveLocalBcc& move) const +{ + //throw exception in case someone accidentaly uses a bcc-BFM move with this feature + std::stringstream errormessage; + errormessage<<"FeatureNNInteractionSc::checkMove(...):\n"; + errormessage<<"attempting to use bcc-BFM move, which is not allowed\n"; + throw std::runtime_error(errormessage.str()); + + return false; +} + + +/** + * @details When a new monomer is added to the system, the lattice sites occupied + * by this monomer must be filled with the attribute tag. This function takes care + * of this. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move of type MoveAddMonomerSc + * @throw std::runtime_error if attribute tag is not in range [1,255] + **/ +template class LatticeClassType> +template +void FeatureNNInteractionSc::applyMove(IngredientsType& ing, + const MoveAddMonomerSc& move) +{ + //get the position and attribute tag of the monomer to be inserted + VectorInt3 pos=move.getPosition(); + VectorInt3 dx(1,0,0); + VectorInt3 dy(0,1,0); + VectorInt3 dz(0,0,1); + lattice_value_type type=lattice_value_type(move.getTag()); + + //the feature is based on a uint8_t lattice, thus the max type must not + //exceed the max value of uint8_t (255) + if(int32_t(type) != move.getTag() || int32_t(type)==0) + { + std::stringstream errormessage; + errormessage<<"FeatureNNInteractionSc::applyMove(MoveAddMonomerSc)\n"; + errormessage<<"Trying to add monomer with type "<maxType=255\n"; + throw std::runtime_error(errormessage.str()); + } + + //update lattice + ing.setLatticeEntry(pos,type); + ing.setLatticeEntry(pos+dx,type); + ing.setLatticeEntry(pos+dy,type); + ing.setLatticeEntry(pos+dx+dy,type); + ing.setLatticeEntry(pos+dz,type); + ing.setLatticeEntry(pos+dz+dx,type); + ing.setLatticeEntry(pos+dz+dy,type); + ing.setLatticeEntry(pos+dz+dx+dy,type); +} + +/** + * @details Because moves of type MoveLocalBcc must not be used with this + * feature, this function always throws an exception when called. The function + * is only implemented for savety purposes. + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move Monte Carlo move of type MoveLocalBcc + * @throw std::runtime_error + **/ +template class LatticeClassType> +template +void FeatureNNInteractionSc::applyMove(const IngredientsType& ing, + const MoveLocalBcc& move) +{ + //throw exception in case someone accidentaly uses a bcc-BFM move with this feature + std::stringstream errormessage; + errormessage<<"FeatureNNInteractionSc::applyMove(...):\n"; + errormessage<<"attempting to use bcc-BFM move, which is not allowed\n"; + throw std::runtime_error(errormessage.str()); + +} + +/** + * @tparam IngredientsType The type of the system including all features + * @param [in] ingredients A reference to the IngredientsType - mainly the system + **/ +template class LatticeClassType> +template +void FeatureNNInteractionSc::synchronize(IngredientsType& ingredients) +{ + + //refill the lattice with attribute tags + //caution: this overwrites, what is currently written on the lattice + fillLattice(ingredients); + +} + + +/** + * @details If not compiled with DEBUG flag this function only returns the content + * of the lookup table probabilityLookup. If compiled with DEBUG flag it checks + * that the attribute tags typeA, typeB are within the allowed range. + * @param typeA monomer attribute type in range [1,255] + * @param typeB monomer attribute type in range [1,255] + * @throw std::runtime_error In debug mode, if types are not in range [1,255] + **/ +template class LatticeClassType> +inline double FeatureNNInteractionSc::getProbabilityFactor(int32_t typeA, + int32_t typeB) const +{ +#ifdef DEBUG + //extra checks only in debug mode, because this is very frequently called + //and this costs performance + if(typeA<0 || typeA>255 || typeB<0 || typeB>255){ + std::stringstream errormessage; + errormessage<<"***FeatureNaNInteractionSc::getInteraction(typeA,typeB)***\n"; + errormessage<<"probability undefined between types "< class LatticeClassType> +template +void FeatureNNInteractionSc::exportRead(FileImport< IngredientsType >& fileReader) +{ + typedef FeatureNNInteractionSc my_type; + fileReader.registerRead("!nn_interaction",new ReadNNInteraction(*this)); +} + + +/** + * The function is called by the Ingredients class when an object of type Ingredients + * is associated with an object of type AnalyzerWriteBfmFile. The export of the Writes is thus + * taken care automatically when it becomes necessary.\n + * Registered Write-Out Commands: + * - !contact_interaction + * + * @tparam IngredientsType The type of the system including all features + * @param fileWriter File writer for the bfm-file. + */ +template class LatticeClassType> +template +void FeatureNNInteractionSc::exportWrite(AnalyzerWriteBfmFile< IngredientsType >& fileWriter) const +{ + typedef FeatureNNInteractionSc my_type; + fileWriter.registerWrite("!nn_interaction",new WriteNNInteraction(*this)); +} + +/** + * @details occupies the lattice with the attribute tags of the monomers + * as this is required to determine the contact interactions in this feature. + * An additional check is performed asserting that the tags are in the range [1,255] + * + * @tparam IngredientsType The type of the system including all features + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @throw std::runtime_error In case a monomer has attribute tag not in [1,255] + **/ +template class LatticeClassType> +template +void FeatureNNInteractionSc::fillLattice(IngredientsType& ingredients) +{ + const typename IngredientsType::molecules_type& molecules=ingredients.getMolecules(); + + for(size_t n=0;n class LatticeClassType> +template +double FeatureNNInteractionSc::calculateAcceptanceProbability( + const IngredientsType& ingredients, + const MoveLocalSc& move) const +{ + + VectorInt3 oldPos=ingredients.getMolecules()[move.getIndex()]; + VectorInt3 direction=move.getDir(); + + double prob=1.0; + int32_t monoType=ingredients.getMolecules()[move.getIndex()].getAttributeTag(); + + /*get two directions perpendicular to vector directon of the move*/ + VectorInt3 perp1,perp2; + /* first perpendicular direction is either (0 1 0) or (1 0 0)*/ + int32_t x1=((direction.getX()==0) ? 1 : 0); + int32_t y1=((direction.getX()!=0) ? 1 : 0); + perp1.setX(x1); + perp1.setY(y1); + perp1.setZ(0); + + /* second perpendicular direction is either (0 0 1) or (0 1 0)*/ + int32_t y2=((direction.getZ()==0) ? 0 : 1); + int32_t z2=((direction.getZ()!=0) ? 0 : 1); + perp2.setX(0); + perp2.setY(y2); + perp2.setZ(z2); + + //the probability is calculated by going through all possible lattice sites + //at which the contacts may have changed. At every site the type of the + //monomer sitting there is retrieved from the lattice. the additional + //factor for the probability (exp(-deltaE/kT)) is retrieved from the + //lookup using getProbabilityFactor. For new contacts this factor is multiplied + //with the probability, for contacts taken away the probability is devided. + VectorInt3 actual=oldPos; + + //first check front,i.e newly acquired contacts + if(direction.getX()>0 || direction.getY()>0 || direction.getZ()>0) actual+=direction; + actual+=direction; + + actual-=perp1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=perp2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual=actual+perp2+perp1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=perp1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual=actual+perp1-perp2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=perp2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual=actual-perp1-perp2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=perp1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual=actual+perp2+direction; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=perp2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=perp1; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=perp2; + prob*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + + //now check back side (contacts taken away) + double prob_div=1.0; + actual=oldPos; + if(direction.getX()<0 || direction.getY()<0 || direction.getZ()<0) actual-=direction; + actual-=perp1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=perp2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual=actual+perp2+perp1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=perp1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual=actual+perp1-perp2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=perp2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual=actual-perp1-perp2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=perp1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual=actual+perp2-direction; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=perp2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual+=perp1; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + actual-=perp2; + prob_div*=getProbabilityFactor(monoType,int32_t(ingredients.getLatticeEntry(actual))); + + prob/=prob_div; + return prob; + +} + +/** + * @param typeA monomer attribute tag in range [1,255] + * @param typeB monomer attribute tag in range [1,255] + * @param interaction energy between typeA and typeB + * @throw std::runtime_error In case typeA or typeB exceed range [1,255] + **/ +template class LatticeClassType> +void FeatureNNInteractionSc::setNNInteraction(int32_t typeA, + int32_t typeB, + double energy) +{ + if(0 class LatticeClassType> +double FeatureNNInteractionSc::getNNInteraction(int32_t typeA, + int32_t typeB) const +{ + + if(0. + +--------------------------------------------------------------------------------*/ + +#ifndef LEMONADE_FEATURE_WALL_H +#define LEMONADE_FEATURE_WALL_H + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + + +/** + * @file + * @brief Enable arbitrary walls in the simulation box by the FeatureWall + * @details For the \b sc-BFM this feature can hold an arbitrary + * number of walls with normal vectors (1,0,0), (0,1,0) and (0,0,1) being the + * unit vectors along the principal lattice directions. + * + * @todo Enable this feature for the \b bcc-BFM. + * */ + +/** + * @class Wall + * @brief class providing a single wall with definitions, setter and getter functions + * */ +class Wall +{ + +public: + //! standard constructor creating an empty wall + Wall(): base(), normal() {} + + //! getter function for the base vector of the wall + const VectorInt3 getBase() const { + return base; + } + + //! setter function for the base vector of the wall + VectorInt3 setBase(uint32_t baseX_, uint32_t baseY_, uint32_t baseZ_) { + base.setAllCoordinates(baseX_,baseY_,baseZ_); + } + + //! getter function for the normal vector of the wall + const VectorInt3 getNormal() const { + return normal; + } + + //! setter function for the normal vector of the wall + VectorInt3 setNormal(uint32_t norX_, uint32_t norY_, uint32_t norZ_) { + VectorInt3 test(norX_, norY_, norZ_); + if(test==VectorInt3(1,0,0) || test==VectorInt3(0,1,0) || test==VectorInt3(0,0,1) ){ + normal.setAllCoordinates(norX_,norY_,norZ_); + }else{ + throw std::runtime_error("normal vector should be P(1,0,0)"); + } + } + +private: + + //! base vector of the wall: arbitrary position somewhere on the wall to provide a well defined wall + VectorInt3 base; + + //! normal vector of the wall: vector perpendiculat to the walls surface + VectorInt3 normal; + +}; + + +/** + * @class FeatureWall + * @brief Feature holding a vector of walls. + * */ +class FeatureWall: public Feature +{ +public: + + //! standard constructor + FeatureWall() {} + + //! standard destructor + virtual ~FeatureWall(){} + + //! read walls from bfm file + template + void exportRead(FileImport& fileReader); + + //! write walls to bfm file + template + void exportWrite(AnalyzerWriteBfmFile& fileWriter) const; + + //! check move function for local sc move + template + bool checkMove(const IngredientsType& ingredients,MoveLocalSc& move); + + //! check move function for add sc move + template + bool checkMove(const IngredientsType& ingredients,MoveAddMonomerSc& addmove); + + //! implemantation of synchronize + template + void synchronize(const IngredientsType& ingredients); + + //! getter function for the walls container + std::vector getWalls() const{ + return walls; + } + + /** + * @brief add function to add a wall to the walls container + * @throw monomer occupies a position on the walls + **/ + void addWall(Wall wall){ + if(wall.getNormal().getLength()!=0.0) + walls.push_back(wall); + else + throw std::runtime_error("wall is not well defined: normal vector has length 0"); + } + + //! empty walls container + void clearAllWalls(){ + walls.clear(); + } + + +private: + //! walls container + std::vector walls; + +}; + + +/*****************************************************************/ +/** + * @class WriteWall + * + * @brief Handles BFM-File-Write \b #!wall + * @tparam [in] IngredientsType Ingredients class storing all system information. + **/ +template +class WriteWall: public AbstractWrite +{ + public: + + WriteWall(const IngredientsType& src):AbstractWrite(src){this->setHeaderOnly(true);} + + void writeStream(std::ostream& strm){ + const IngredientsType& ingredients=(this->getSource()); + + for (size_t i=0; i +class ReadWall : public ReadToDestination < IngredientsType > +{ + public: + + ReadWall(IngredientsType& destination):ReadToDestination< IngredientsType > (destination){} + + void execute(); +}; + + +/** + * @details checking if move is going to touch one of the walls + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] move local sc move + */ +template +bool FeatureWall::checkMove(const IngredientsType& ingredients, MoveLocalSc& move) +{ + uint32_t counter(0); + + for (size_t i = 0; i < ingredients.getWalls().size(); i++) { //check all walls in the system + + uint8_t direction; + VectorInt3 unitX(1,0,0); + VectorInt3 unitY(0,1,0); + VectorInt3 unitZ(0,0,1); + + if (ingredients.getWalls()[i].getNormal() == unitX) { direction = 0; } //points direction of normal and sets variable direction to look only on this part of the new position vector later + if (ingredients.getWalls()[i].getNormal() == unitY) { direction = 1; } + if (ingredients.getWalls()[i].getNormal() == unitZ) { direction = 2; } + + if ((ingredients.getMolecules()[move.getIndex()] + move.getDir()).getCoordinate(direction) != ingredients.getWalls()[i].getBase().getCoordinate(direction)-1) { //if new position has NOT the same value in "direction" as the wall -1, permit the move + counter++; + } + } + + if (counter == ingredients.getWalls().size()) {return true;} + + return false; +} + +/** + * @details checking if added monomer is going to touch one of the walls + * + * @param [in] ingredients A reference to the IngredientsType - mainly the system + * @param [in] addmove move to add sc monomer + */ +template +bool FeatureWall::checkMove(const IngredientsType& ingredients, MoveAddMonomerSc& addmove) +{ + uint32_t counter(0); + + for (size_t i = 0; i < ingredients.getWalls().size(); i++) { //check all walls in the system + + uint8_t direction; + VectorInt3 unitX(1,0,0); + VectorInt3 unitY(0,1,0); + VectorInt3 unitZ(0,0,1); + + if (ingredients.getWalls()[i].getNormal() == unitX) { direction = 0; } //points direction of normal and sets variable direction to look only on this part of the new position vector later + if (ingredients.getWalls()[i].getNormal() == unitY) { direction = 1; } + if (ingredients.getWalls()[i].getNormal() == unitZ) { direction = 2; } + + if (addmove.getPosition().getCoordinate(direction) != ingredients.getWalls()[i].getBase().getCoordinate(direction)-1) { //if new molecule has NOT the same value in "direction" as the wall -1, permit the move + counter++; + } + } + + if (counter == ingredients.getWalls().size()) {return true;} + + return false; +} + +/** + * @brief Synchronize this feature with the system given as argument + * + * @details checking all monomer positions to be not in conflict with one of the walls + * + * @throw monomer occupies a position on the walls + * @param [in] ingredients a reference to the IngredientsType - mainly the system + **/ +template +void FeatureWall::synchronize(const IngredientsType& ingredients) +{ + uint8_t direction; + VectorInt3 unitX(1,0,0); + VectorInt3 unitY(0,1,0); + VectorInt3 unitZ(0,0,1); + + for (size_t i=0; i walls could not be read. + * @tparam IngredientsType Features used in the system. See Ingredients. + **/ +template < class IngredientsType > +void ReadWall::execute() +{ + std::cout << "reading #!wall..."; + + + uint32_t baseX, baseY, baseZ, normalX, normalY, normalZ; + + this->getInputStream() >> baseX >> baseY >> baseZ >> normalX >> normalY >> normalZ; + + if(!this->getInputStream().fail()) + { + Wall wall; + wall.setBase(baseX,baseY,baseZ); + wall.setNormal(normalX,normalY,normalZ); + this->getDestination().addWall(wall); + + std::cout << "Base: " << baseX << " " << baseY << " " << baseZ << ", Normal: " << normalX << " " << normalY << " " << normalZ << std::endl; + } + else + throw std::runtime_error("ReadWall::execute()\n Could not read wall"); +} + + +/** + * @details The function is called by the Ingredients class when an object of type Ingredients + * is associated with an object of type FileImport. The export of the Reads is thus + * taken care automatically when it becomes necessary.\n + * Registered Read-In Commands: + * * #!wall: bx by bz nx ny nz + * + * @param fileReader File importer for the bfm-file + * @tparam IngredientsType Features used in the system. See Ingredients. + **/ +template < class IngredientsType > +void FeatureWall::exportRead(FileImport& fileReader) +{ + IngredientsType& destination=fileReader.getDestination(); + fileReader.registerRead("#!wall:", new ReadWall< IngredientsType > (destination)); +} + +/** + * The function is called by the Ingredients class when an object of type Ingredients + * is associated with an object of type AnalyzerWriteBfmFile. The export of the Writes is thus + * taken care automatically when it becomes necessary.\n + * Registered Write-Out Commands: + * * #!wall: bx by bz nx ny nz + * + * @param fileWriter File writer for the bfm-file. + * @tparam IngredientsType Features used in the system. See Ingredients. + */ +template < class IngredientsType > +void FeatureWall::exportWrite(AnalyzerWriteBfmFile& fileWriter) const +{ + const IngredientsType& source=fileWriter.getIngredients_(); + fileWriter.registerWrite("#!wall:", new WriteWall (source)); +} + +#endif //LEMONADE_FEATURE_WALL_H diff --git a/include/LeMonADE/updater/UpdaterAddLinearChains.h b/include/LeMonADE/updater/UpdaterAddLinearChains.h new file mode 100644 index 0000000..38bd7f2 --- /dev/null +++ b/include/LeMonADE/updater/UpdaterAddLinearChains.h @@ -0,0 +1,173 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ + +#ifndef LEMONADE_UPDATER_SETUP_LINEAR_CHAINS +#define LEMONADE_UPDATER_SETUP_LINEAR_CHAINS +/** + * @file + * + * @class UpdaterAddLinearChains + * + * @brief Updater to create a solution of monomdisperse linear chains. + * + * @details This is a simple implementation of a system setup starting from an empty ingredients + * or a system with some monomers inside. This updater requires FeatureAttributes. + * Two tags are added to the monomers in alternating manner, usually needed for GPU computing. + * + * @tparam IngredientsType + * + * @param ingredients_ The system, holding eigther an empty simulation box for system setup + * or a prefilled ingredients where the linear chains shall be added + * @param NChain_ number of chains that are added to ingredients + * @param NMonoPerChain_ number of monomer is each chain + * @param type1_ attribute tag of "even" monomers + * @param type2_ attribute tag of "odd" monomers + **/ + +#include +#include + +template +class UpdaterAddLinearChains: public UpdaterAbstractCreate +{ + typedef UpdaterAbstractCreate BaseClass; + +public: + UpdaterAddLinearChains(IngredientsType& ingredients_, uint32_t NChain_, uint32_t NMonoPerChain_, int32_t type1_=1, int32_t type2_=2); + + virtual void initialize(); + virtual bool execute(); + virtual void cleanup(); + +private: + // provide access to functions of UpdaterAbstractCreate used in this updater + using BaseClass::ingredients; + using BaseClass::addMonomerToParent; + using BaseClass::addSingleMonomer; + using BaseClass::linearizeSystem; + + //! number of monomers in a chain + uint32_t NMonoPerChain; + + //!number of linear chains in the box + uint32_t NChain; + + //! lattice occupation density + double density; + + //! bool for execution + bool wasExecuted; + + //! attribute tag of even monomers + int32_t type1; + + //!getAttributeTag of odd monomers + int32_t type2; + +}; + +/** +* @brief Constructor handling the new systems paramters +* +* @param ingredients_ a reference to the IngredientsType - mainly the system +* @param NChain_ number of chains to be added in the system instead of solvent +* @param NMonoPerChain_ number of monomers in one chain +*/ +template < class IngredientsType > +UpdaterAddLinearChains::UpdaterAddLinearChains(IngredientsType& ingredients_, uint32_t NChain_, uint32_t NMonoPerChain_, int32_t type1_, int32_t type2_): +BaseClass(ingredients_), NChain(NChain_), NMonoPerChain(NMonoPerChain_), density(0), wasExecuted(false), +type1(type1_), type2(type2_) +{} + +/** +* @brief initialise function, calculate the target density to compare with at the end. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template < class IngredientsType > +void UpdaterAddLinearChains::initialize(){ + std::cout << "initialize UpdaterAddLinearChains" << std::endl; + + // get the target density from the sum of existing monomers and the new added chains + density=(double)( ingredients.getMolecules().size() + NMonoPerChain*NChain ) * 8 /(double)( ingredients.getBoxX()*ingredients.getBoxY()*ingredients.getBoxZ() ); + + std::cout << "add "< +bool UpdaterAddLinearChains::execute(){ + if(wasExecuted) + return true; + + std::cout << "execute UpdaterAddLinearChains" << std::endl; + + //loop over chains and chain monomers and build it up + for(uint32_t i=0;i<(NChain);i++){ + for(uint32_t j=0;j<(NMonoPerChain);j++){ + if(j==0) + addSingleMonomer(type1); + else{ + if(ingredients.getMolecules()[ingredients.getMolecules().size()-1].getAttributeTag() == type1) + addMonomerToParent(ingredients.getMolecules().size()-1,type2); + else + addMonomerToParent(ingredients.getMolecules().size()-1,type1); + } + } + } + + ingredients.synchronize(); + double lattice_volume(ingredients.getBoxX()*ingredients.getBoxY()*ingredients.getBoxZ()); + if(density =! ( (double)(ingredients.getMolecules().size()*8) / lattice_volume ) ){ + std::cout << density << " " <<( (ingredients.getMolecules().size()*8) / lattice_volume)< +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +class NNInteractionBccTest: public ::testing::Test{ + /* suppress cout output for better readability -->un-/comment here:*/ +public: + //redirect cout output + virtual void SetUp(){ + originalBuffer=cout.rdbuf(); + cout.rdbuf(tempStream.rdbuf()); + }; + //restore original output + virtual void TearDown(){ + cout.rdbuf(originalBuffer); + }; +private: + std::streambuf* originalBuffer; + std::ostringstream tempStream; + /* ** */ +}; + + +TEST_F(NNInteractionBccTest,CheckApplyBccMovePowerOfTwoLattice) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureNNInteractionBcc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients1; + + + + //prepare myIngredients1 + myIngredients1.setBoxX(64); + myIngredients1.setBoxY(64); + myIngredients1.setBoxZ(64); + myIngredients1.setPeriodicX(1); + myIngredients1.setPeriodicY(1); + myIngredients1.setPeriodicZ(1); + + //set interaction between types 1,2 + myIngredients1.setNNInteraction(1,2,0.8); + double epsilon0=0.8; + + //add two monomers. the test then proceeds like this: monomer 1 is moved + //to all positions around monomer 0 that give contributions to the + //acceptance probability. the probability is checked using a testmove + //and compared with the expected value. in the testmove, monomer 0 is + //simply moved one step in positive x-direction, and then back to its + //original position. + //the reverse move is checked using + //a different testmove (backmove) and also compared with the expectation. + //the monomer is then moved to a new position using updatemove. + //all used moves are of type MoveLocalBcc + typename Ing1::molecules_type& molecules1=myIngredients1.modifyMolecules(); + molecules1.resize(2); + molecules1[0].setAllCoordinates(10,10,10); + molecules1[0].setAttributeTag(1); + molecules1[1].setAttributeTag(2); + + molecules1[1].setAllCoordinates(13,7,7); + myIngredients1.synchronize(myIngredients1); + + //define moves in all possible directions for later use + MoveLocalBcc DXDYDZ,DXDYmDZ,DXmDYDZ,DXmDYmDZ,mDXDYDZ,mDXDYmDZ,mDXmDYDZ,mDXmDYmDZ; + + DXDYDZ.init(myIngredients1); + while((DXDYDZ.getDir()[0]!=1 ) ||(DXDYDZ.getDir()[1]!=1 ) ||(DXDYDZ.getDir()[2]!=1 ) || (DXDYDZ.getIndex()!=1)){ DXDYDZ.init(myIngredients1);} + + DXDYmDZ.init(myIngredients1); + while((DXDYmDZ.getDir()[0]!=1 ) ||(DXDYmDZ.getDir()[1]!=1 ) ||(DXDYmDZ.getDir()[2]!=-1 ) || (DXDYmDZ.getIndex()!=1)){ DXDYmDZ.init(myIngredients1);} + + DXmDYDZ.init(myIngredients1); + while((DXmDYDZ.getDir()[0]!=1 ) ||(DXmDYDZ.getDir()[1]!=-1 ) ||(DXmDYDZ.getDir()[2]!=1 ) || (DXmDYDZ.getIndex()!=1)){ DXmDYDZ.init(myIngredients1);} + + DXmDYmDZ.init(myIngredients1); + while((DXmDYmDZ.getDir()[0]!=1 ) ||(DXmDYmDZ.getDir()[1]!=-1 ) ||(DXmDYmDZ.getDir()[2]!=-1 ) || (DXmDYmDZ.getIndex()!=1)){ DXmDYmDZ.init(myIngredients1);} + + mDXDYDZ.init(myIngredients1); + while((mDXDYDZ.getDir()[0]!=-1 ) ||(mDXDYDZ.getDir()[1]!=1 ) ||(mDXDYDZ.getDir()[2]!=1 ) || (mDXDYDZ.getIndex()!=1)){ mDXDYDZ.init(myIngredients1);} + + mDXDYmDZ.init(myIngredients1); + while((mDXDYmDZ.getDir()[0]!=-1 ) ||(mDXDYmDZ.getDir()[1]!=1 ) ||(mDXDYmDZ.getDir()[2]!=-1 ) || (mDXDYmDZ.getIndex()!=1)){ mDXDYmDZ.init(myIngredients1);} + + mDXmDYDZ.init(myIngredients1); + while((mDXmDYDZ.getDir()[0]!=-1 ) ||(mDXmDYDZ.getDir()[1]!=-1 ) ||(mDXmDYDZ.getDir()[2]!=1 ) || (mDXmDYDZ.getIndex()!=1)){ mDXmDYDZ.init(myIngredients1);} + + mDXmDYmDZ.init(myIngredients1); + while((mDXmDYmDZ.getDir()[0]!=-1 ) ||(mDXmDYmDZ.getDir()[1]!=-1 ) ||(mDXmDYmDZ.getDir()[2]!=-1 ) || (mDXmDYmDZ.getIndex()!=1)){ mDXmDYmDZ.init(myIngredients1);} + + //now start checking moves of monomer 1 around the position of monomer 0 + + ///--- first move upwards in a zig-zag line from (3,-3,-3) to (3,-1,3) relative to monomer 0 + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),8); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),9); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),10); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),11); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),12); + + //last step to the next row + DXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYDZ.getProbability(),exp(epsilon0)); + DXDYDZ.apply(myIngredients1); + DXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),13); + + //----- now move downwards in zig-zag from (3,-1,3) to (3,1,-3) + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(-epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),12); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),11); + + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(-epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),10); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),9); + + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(-epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),8); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),7); + + //------ now move upwards in zig-zag to (1,1,3) + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),8); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),9); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),10); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),11); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),12); + + //last step to the next row + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),11); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),13); + + //--- now again zig-za downwards to (-1,3,-3) + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(-epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),12); + + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),13); + EXPECT_EQ(molecules1[1].getZ(),11); + + DXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYmDZ.getProbability(),exp(-epsilon0)); + DXmDYmDZ.apply(myIngredients1); + DXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),10); + + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),13); + EXPECT_EQ(molecules1[1].getZ(),9); + + DXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYmDZ.getProbability(),exp(-epsilon0)); + DXmDYmDZ.apply(myIngredients1); + DXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),8); + + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),13); + EXPECT_EQ(molecules1[1].getZ(),7); + + //---move up again to (-3,1,3) relative to monomer 0 + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(-epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),8); + + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),9); + + DXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYDZ.getProbability(),exp(-epsilon0)); + DXDYDZ.apply(myIngredients1); + DXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),10); + + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),11); + + DXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYDZ.getProbability(),exp(-epsilon0)); + DXDYDZ.apply(myIngredients1); + DXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),12); + + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),13); + + //---move down to (-2,0,-2), which has a contact, then move to (-1,1,-3), which doesn't + DXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYmDZ.getProbability(),exp(-epsilon0)); + DXmDYmDZ.apply(myIngredients1); + DXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),12); + + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),11); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(-epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),10); + + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),9); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(-epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),8); + + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),7); + + //-- now check the position below monomer 0, i.e. (0,0,-2) + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(-epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),8); + + //--now check the position (0,-2,-2), go there in a first step + DXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYmDZ.getProbability(),exp(epsilon0)); + DXmDYmDZ.apply(myIngredients1); + DXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),11); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),7); + + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(-epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),8); + + //--now check (-2,-2,-2), go there in a first step + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),7); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),8); + + //--check -2,-2,0, go there in a first step + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),9); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(-epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),10); + + //-- check (0,-2,0), go there first + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),11); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(-epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),10); + + //--now check (0,-2,2), go there in one step first + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),11); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),11); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),12); + + //-- now the last unchecked position with a contact is (0,0,2) + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),13); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(-epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),12); +} + +//same test but with any lattice + +TEST_F(NNInteractionBccTest,CheckApplyBccMoveAnyLattice) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureNNInteractionBcc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients1; + + + + //prepare myIngredients1 + myIngredients1.setBoxX(65); + myIngredients1.setBoxY(65); + myIngredients1.setBoxZ(65); + myIngredients1.setPeriodicX(1); + myIngredients1.setPeriodicY(1); + myIngredients1.setPeriodicZ(1); + + //set interaction between types 1,2 + myIngredients1.setNNInteraction(1,2,0.8); + double epsilon0=0.8; + + //add two monomers. the test then proceeds like this: monomer 1 is moved + //to all positions around monomer 0 that give contributions to the + //acceptance probability. the probability is checked using a testmove + //and compared with the expected value. in the testmove, monomer 0 is + //simply moved one step in positive x-direction, and then back to its + //original position. + //the reverse move is checked using + //a different testmove (backmove) and also compared with the expectation. + //the monomer is then moved to a new position using updatemove. + //all used moves are of type MoveLocalBcc + typename Ing1::molecules_type& molecules1=myIngredients1.modifyMolecules(); + molecules1.resize(2); + molecules1[0].setAllCoordinates(10,10,10); + molecules1[0].setAttributeTag(1); + molecules1[1].setAttributeTag(2); + + molecules1[1].setAllCoordinates(13,7,7); + myIngredients1.synchronize(myIngredients1); + + //define moves in all possible directions for later use + MoveLocalBcc DXDYDZ,DXDYmDZ,DXmDYDZ,DXmDYmDZ,mDXDYDZ,mDXDYmDZ,mDXmDYDZ,mDXmDYmDZ; + + DXDYDZ.init(myIngredients1); + while((DXDYDZ.getDir()[0]!=1 ) ||(DXDYDZ.getDir()[1]!=1 ) ||(DXDYDZ.getDir()[2]!=1 ) || (DXDYDZ.getIndex()!=1)){ DXDYDZ.init(myIngredients1);} + + DXDYmDZ.init(myIngredients1); + while((DXDYmDZ.getDir()[0]!=1 ) ||(DXDYmDZ.getDir()[1]!=1 ) ||(DXDYmDZ.getDir()[2]!=-1 ) || (DXDYmDZ.getIndex()!=1)){ DXDYmDZ.init(myIngredients1);} + + DXmDYDZ.init(myIngredients1); + while((DXmDYDZ.getDir()[0]!=1 ) ||(DXmDYDZ.getDir()[1]!=-1 ) ||(DXmDYDZ.getDir()[2]!=1 ) || (DXmDYDZ.getIndex()!=1)){ DXmDYDZ.init(myIngredients1);} + + DXmDYmDZ.init(myIngredients1); + while((DXmDYmDZ.getDir()[0]!=1 ) ||(DXmDYmDZ.getDir()[1]!=-1 ) ||(DXmDYmDZ.getDir()[2]!=-1 ) || (DXmDYmDZ.getIndex()!=1)){ DXmDYmDZ.init(myIngredients1);} + + mDXDYDZ.init(myIngredients1); + while((mDXDYDZ.getDir()[0]!=-1 ) ||(mDXDYDZ.getDir()[1]!=1 ) ||(mDXDYDZ.getDir()[2]!=1 ) || (mDXDYDZ.getIndex()!=1)){ mDXDYDZ.init(myIngredients1);} + + mDXDYmDZ.init(myIngredients1); + while((mDXDYmDZ.getDir()[0]!=-1 ) ||(mDXDYmDZ.getDir()[1]!=1 ) ||(mDXDYmDZ.getDir()[2]!=-1 ) || (mDXDYmDZ.getIndex()!=1)){ mDXDYmDZ.init(myIngredients1);} + + mDXmDYDZ.init(myIngredients1); + while((mDXmDYDZ.getDir()[0]!=-1 ) ||(mDXmDYDZ.getDir()[1]!=-1 ) ||(mDXmDYDZ.getDir()[2]!=1 ) || (mDXmDYDZ.getIndex()!=1)){ mDXmDYDZ.init(myIngredients1);} + + mDXmDYmDZ.init(myIngredients1); + while((mDXmDYmDZ.getDir()[0]!=-1 ) ||(mDXmDYmDZ.getDir()[1]!=-1 ) ||(mDXmDYmDZ.getDir()[2]!=-1 ) || (mDXmDYmDZ.getIndex()!=1)){ mDXmDYmDZ.init(myIngredients1);} + + //now start checking moves of monomer 1 around the position of monomer 0 + + ///--- first move upwards in a zig-zag line from (3,-3,-3) to (3,-1,3) relative to monomer 0 + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),8); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),9); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),10); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),11); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),12); + + //last step to the next row + DXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYDZ.getProbability(),exp(epsilon0)); + DXDYDZ.apply(myIngredients1); + DXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),13); + + //----- now move downwards in zig-zag from (3,-1,3) to (3,1,-3) + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(-epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),12); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),11); + + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(-epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),10); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),9); + + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(-epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),8); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),7); + + //------ now move upwards in zig-zag to (1,1,3) + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),8); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),9); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),10); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),13); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),11); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),12); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),12); + + //last step to the next row + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),11); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),13); + + //--- now again zig-za downwards to (-1,3,-3) + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(-epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),12); + + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),13); + EXPECT_EQ(molecules1[1].getZ(),11); + + DXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYmDZ.getProbability(),exp(-epsilon0)); + DXmDYmDZ.apply(myIngredients1); + DXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),10); + + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),13); + EXPECT_EQ(molecules1[1].getZ(),9); + + DXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYmDZ.getProbability(),exp(-epsilon0)); + DXmDYmDZ.apply(myIngredients1); + DXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),8); + + mDXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYmDZ.getProbability(),exp(epsilon0)); + mDXDYmDZ.apply(myIngredients1); + mDXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),13); + EXPECT_EQ(molecules1[1].getZ(),7); + + //---move up again to (-3,1,3) relative to monomer 0 + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(-epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),8); + + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),9); + + DXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYDZ.getProbability(),exp(-epsilon0)); + DXDYDZ.apply(myIngredients1); + DXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),10); + + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),11); + + DXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYDZ.getProbability(),exp(-epsilon0)); + DXDYDZ.apply(myIngredients1); + DXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),12); + EXPECT_EQ(molecules1[1].getZ(),12); + + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),13); + + //---move down to (-2,0,-2), which has a contact, then move to (-1,1,-3), which doesn't + DXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYmDZ.getProbability(),exp(-epsilon0)); + DXmDYmDZ.apply(myIngredients1); + DXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),12); + + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),11); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(-epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),10); + + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),9); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(-epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),8); + + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),11); + EXPECT_EQ(molecules1[1].getZ(),7); + + //-- now check the position below monomer 0, i.e. (0,0,-2) + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(-epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),8); + + //--now check the position (0,-2,-2), go there in a first step + DXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYmDZ.getProbability(),exp(epsilon0)); + DXmDYmDZ.apply(myIngredients1); + DXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),11); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),7); + + mDXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYDZ.getProbability(),exp(-epsilon0)); + mDXmDYDZ.apply(myIngredients1); + mDXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),8); + + //--now check (-2,-2,-2), go there in a first step + mDXmDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXmDYmDZ.getProbability(),exp(epsilon0)); + mDXmDYmDZ.apply(myIngredients1); + mDXmDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),7); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),8); + + //--check -2,-2,0, go there in a first step + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),7); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),9); + + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(-epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),8); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),10); + + //-- check (0,-2,0), go there first + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),11); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(-epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),10); + + //--now check (0,-2,2), go there in one step first + DXmDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXmDYDZ.getProbability(),exp(epsilon0)); + DXmDYDZ.apply(myIngredients1); + DXmDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),11); + EXPECT_EQ(molecules1[1].getY(),7); + EXPECT_EQ(molecules1[1].getZ(),11); + + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(-epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),8); + EXPECT_EQ(molecules1[1].getZ(),12); + + //-- now the last unchecked position with a contact is (0,0,2) + mDXDYDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(mDXDYDZ.getProbability(),exp(epsilon0)); + mDXDYDZ.apply(myIngredients1); + mDXDYDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),9); + EXPECT_EQ(molecules1[1].getY(),9); + EXPECT_EQ(molecules1[1].getZ(),13); + + DXDYmDZ.check(myIngredients1); + EXPECT_DOUBLE_EQ(DXDYmDZ.getProbability(),exp(-epsilon0)); + DXDYmDZ.apply(myIngredients1); + DXDYmDZ.resetProbability(); + EXPECT_EQ(molecules1[1].getX(),10); + EXPECT_EQ(molecules1[1].getY(),10); + EXPECT_EQ(molecules1[1].getZ(),12); +} + + + + + +TEST_F(NNInteractionBccTest,getSetInteraction) +{ + typedef LOKI_TYPELIST_2(FeatureBondset<>,FeatureNNInteractionBcc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients; + + for(int32_t i=1;i<=255;i++) + { + for(int32_t j=1;j<=255;j++) + { + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(i,j),0.0); + } + } + EXPECT_THROW(myIngredients.getNNInteraction(-1,5),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(1,-5),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(0,1),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(1,0),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(0,256),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(256,0),std::runtime_error); + + //set interaction between types 1,2 + myIngredients.setNNInteraction(1,2,0.8); + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(1,2),0.8); + + myIngredients.setNNInteraction(1,2,-0.8); + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(1,2),-0.8); + + + myIngredients.setNNInteraction(1,2,0.0); + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(1,2),0.0); +} + +TEST_F(NNInteractionBccTest,Synchronize) +{ + typedef LOKI_TYPELIST_2(FeatureBondset<>,FeatureNNInteractionBcc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients1; + //prepare myIngredients1 + myIngredients1.setBoxX(32); + myIngredients1.setBoxY(32); + myIngredients1.setBoxZ(32); + myIngredients1.setPeriodicX(1); + myIngredients1.setPeriodicY(1); + myIngredients1.setPeriodicZ(1); + + typename Ing1::molecules_type& molecules1=myIngredients1.modifyMolecules(); + molecules1.resize(2); + molecules1[0].setAllCoordinates(9,9,9); + molecules1[1].setAllCoordinates(1,1,1); + molecules1[0].setAttributeTag(1); + molecules1[1].setAttributeTag(2); + + myIngredients1.synchronize(myIngredients1); + + VectorInt3 pos; + //go through lattice and count contacts at each position + for(int32_t x=0;x<32;x++) + { + for(int32_t y=0;y<32;y++) + { + for(int32_t z=0;z<32;z++) + { + //get the lattice entry + pos.setAllCoordinates(x,y,z); + if(pos==VectorInt3(9,9,9)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(1,1,1)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else + EXPECT_EQ(0,myIngredients1.getLatticeEntry(pos)); + + + } + } + } + + molecules1.resize(3); + molecules1[2].setAllCoordinates(1,1,1); + molecules1[2].setAttributeTag(2); + EXPECT_THROW(myIngredients1.synchronize(myIngredients1),std::runtime_error); + + +} + +TEST_F(NNInteractionBccTest,ReadWrite) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureNNInteractionBcc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients; + + + //set interaction between types 1,2 + myIngredients.setNNInteraction(1,2,0.8); + myIngredients.setNNInteraction(1,3,-0.8); + myIngredients.setNNInteraction(2,3,0.0); + myIngredients.setNNInteraction(9,4,10.0); + + //prepare myIngredients + myIngredients.setBoxX(32); + myIngredients.setBoxY(32); + myIngredients.setBoxZ(32); + myIngredients.setPeriodicX(1); + myIngredients.setPeriodicY(1); + myIngredients.setPeriodicZ(1); + + typename Ing1::molecules_type& molecules1=myIngredients.modifyMolecules(); + molecules1.resize(2); + molecules1[0].setAllCoordinates(9,9,9); + molecules1[1].setAllCoordinates(1,1,1); + molecules1[0].setAttributeTag(1); + molecules1[1].setAttributeTag(2); + + myIngredients.synchronize(myIngredients); + + AnalyzerWriteBfmFile outfile("./interactionRW.test",myIngredients,AnalyzerWriteBfmFile::NEWFILE); + outfile.initialize(); + outfile.execute(); + + Ing1 myIngredients2; + UpdaterReadBfmFile infile("./interactionRW.test",myIngredients2,UpdaterReadBfmFile::READ_LAST_CONFIG_FAST); + infile.initialize(); + + for(int32_t i=1;i<=255;i++) + { + for(int32_t j=1;j<=255;j++) + { + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(i,j),myIngredients2.getNNInteraction(i,j)); + } + } + + outfile.closeFile(); + infile.closeFile(); + //remove the temporary file + EXPECT_EQ(0,remove("./interactionRW.test")); + + +} + + +TEST_F(NNInteractionBccTest,ApplyMoveAddMonomerBcc) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureNNInteractionBcc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients1; + + //prepare myIngredients1 + myIngredients1.setBoxX(65); + myIngredients1.setBoxY(65); + myIngredients1.setBoxZ(65); + myIngredients1.setPeriodicX(1); + myIngredients1.setPeriodicY(1); + myIngredients1.setPeriodicZ(1); + + myIngredients1.synchronize(); + + typename Ing1::molecules_type& molecules1=myIngredients1.modifyMolecules(); + + MoveAddMonomerBcc addMonomer; + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(256); + EXPECT_THROW(addMonomer.apply(myIngredients1),std::runtime_error); + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(0); + EXPECT_THROW(addMonomer.apply(myIngredients1),std::runtime_error); + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(-1); + EXPECT_THROW(addMonomer.apply(myIngredients1),std::runtime_error); + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(2); + addMonomer.apply(myIngredients1); + + EXPECT_EQ(2,myIngredients1.getLatticeEntry(5,6,7)); + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(255); + addMonomer.apply(myIngredients1); + + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(5,6,7))); + +} diff --git a/tests/feature/TestFeatureNNInteractionSc.cpp b/tests/feature/TestFeatureNNInteractionSc.cpp new file mode 100644 index 0000000..c633ba0 --- /dev/null +++ b/tests/feature/TestFeatureNNInteractionSc.cpp @@ -0,0 +1,2370 @@ +#include + +#include "gtest/gtest.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +class NNInteractionScTest: public ::testing::Test{ + /* suppress cout output for better readability -->un-/comment here:*/ +public: + //redirect cout output + virtual void SetUp(){ + originalBuffer=cout.rdbuf(); + cout.rdbuf(tempStream.rdbuf()); + }; + //restore original output + virtual void TearDown(){ + cout.rdbuf(originalBuffer); + }; +private: + std::streambuf* originalBuffer; + std::ostringstream tempStream; + /* ** */ +}; + + +TEST_F(NNInteractionScTest,CheckApplyScMovePowerOfTwoLattice) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureNNInteractionSc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients1; + + + + //prepare myIngredients1 + myIngredients1.setBoxX(64); + myIngredients1.setBoxY(64); + myIngredients1.setBoxZ(64); + myIngredients1.setPeriodicX(1); + myIngredients1.setPeriodicY(1); + myIngredients1.setPeriodicZ(1); + + //set interaction between types 1,2 + myIngredients1.setNNInteraction(1,2,0.8); + double epsilon0=0.8; + + //add two monomers. the test then proceeds like this: monomer 1 is moved + //to all positions around monomer 0 that give contributions to the + //acceptance probability. the probability is checked using a testmove + //and compared with the expected value. in the testmove, monomer 0 is + //simply moved one step in positive x-direction, and then back to its + //original position. + //the reverse move is checked using + //a different testmove (backmove) and also compared with the expectation. + //the monomer is then moved to a new position using updatemove. + //all used moves are of type MoveLocalSc + typename Ing1::molecules_type& molecules1=myIngredients1.modifyMolecules(); + molecules1.resize(2); + molecules1[0].setAllCoordinates(9,10,10); + molecules1[0].setAttributeTag(1); + molecules1[1].setAttributeTag(2); + + + //now check EV for move in positive x-direction + MoveLocalSc testmove,backmove,updateMove; + + molecules1[1].setAllCoordinates(12,9,10); + myIngredients1.synchronize(myIngredients1); + + //cout<<"-------------> checking single Interaction same hight <--------------------"< checking single Interaction above <--------------------"< checking single Interaction below <--------------------"< checking single Interaction twice below <--------------------"< checking single Interaction twice above <--------------------"<,FeatureNNInteractionSc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients1; + + + + //prepare myIngredients1 + myIngredients1.setBoxX(65); + myIngredients1.setBoxY(65); + myIngredients1.setBoxZ(65); + myIngredients1.setPeriodicX(1); + myIngredients1.setPeriodicY(1); + myIngredients1.setPeriodicZ(1); + + //set interaction between types 1,2 + myIngredients1.setNNInteraction(1,2,0.8); + double epsilon0=0.8; + + //add two monomers. the test then proceeds like this: monomer 1 is moved + //to all positions around monomer 0 that give contributions to the + //acceptance probability. the probability is checked using a testmove + //and compared with the expected value. in the testmove, monomer 0 is + //simply moved one step in positive x-direction, and then back to its + //original position. + //the reverse move is checked using + //a different testmove (backmove) and also compared with the expectation. + //the monomer is then moved to a new position using updatemove. + //all used moves are of type MoveLocalSc + typename Ing1::molecules_type& molecules1=myIngredients1.modifyMolecules(); + molecules1.resize(2); + molecules1[0].setAllCoordinates(9,10,10); + molecules1[0].setAttributeTag(1); + molecules1[1].setAttributeTag(2); + + + //now check EV for move in positive x-direction + MoveLocalSc testmove,backmove,updateMove; + + molecules1[1].setAllCoordinates(12,9,10); + myIngredients1.synchronize(myIngredients1); + + //cout<<"-------------> checking single Interaction same hight <--------------------"< checking single Interaction above <--------------------"< checking single Interaction below <--------------------"< checking single Interaction twice below <--------------------"< checking single Interaction twice above <--------------------"<,FeatureNNInteractionSc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients; + + for(int32_t i=1;i<=255;i++) + { + for(int32_t j=1;j<=255;j++) + { + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(i,j),0.0); + } + } + EXPECT_THROW(myIngredients.getNNInteraction(-1,5),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(1,-5),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(0,1),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(1,0),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(0,256),std::runtime_error); + EXPECT_THROW(myIngredients.getNNInteraction(256,0),std::runtime_error); + + //set interaction between types 1,2 + myIngredients.setNNInteraction(1,2,0.8); + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(1,2),0.8); + + myIngredients.setNNInteraction(1,2,-0.8); + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(1,2),-0.8); + + + myIngredients.setNNInteraction(1,2,0.0); + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(1,2),0.0); +} + +TEST_F(NNInteractionScTest,Synchronize) +{ + typedef LOKI_TYPELIST_2(FeatureBondset<>,FeatureNNInteractionSc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients1; + //prepare myIngredients1 + myIngredients1.setBoxX(32); + myIngredients1.setBoxY(32); + myIngredients1.setBoxZ(32); + myIngredients1.setPeriodicX(1); + myIngredients1.setPeriodicY(1); + myIngredients1.setPeriodicZ(1); + + typename Ing1::molecules_type& molecules1=myIngredients1.modifyMolecules(); + molecules1.resize(2); + molecules1[0].setAllCoordinates(9,10,10); + molecules1[1].setAllCoordinates(1,1,1); + molecules1[0].setAttributeTag(1); + molecules1[1].setAttributeTag(2); + + myIngredients1.synchronize(myIngredients1); + + VectorInt3 pos; + //go through lattice and count contacts at each position + for(int32_t x=0;x<32;x++) + { + for(int32_t y=0;y<32;y++) + { + for(int32_t z=0;z<32;z++) + { + //get the lattice entry + pos.setAllCoordinates(x,y,z); + if(pos==VectorInt3(9,10,10)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(9,11,10)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(9,10,11)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(9,11,11)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(10,10,10)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(10,11,10)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(10,10,11)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(10,11,11)) + EXPECT_EQ(1,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(1,1,1)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(1,2,1)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(1,1,2)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(1,2,2)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(2,1,1)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(2,2,1)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(2,1,2)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else if(pos==VectorInt3(2,2,2)) + EXPECT_EQ(2,myIngredients1.getLatticeEntry(pos)); + else + EXPECT_EQ(0,myIngredients1.getLatticeEntry(pos)); + + + } + } + } + + molecules1.resize(3); + molecules1[2].setAllCoordinates(1,1,1); + molecules1[2].setAttributeTag(2); + EXPECT_THROW(myIngredients1.synchronize(myIngredients1),std::runtime_error); + + +} + +TEST_F(NNInteractionScTest,ReadWrite) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureNNInteractionSc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients; + + + //set interaction between types 1,2 + myIngredients.setNNInteraction(1,2,0.8); + myIngredients.setNNInteraction(1,3,-0.8); + myIngredients.setNNInteraction(2,3,0.0); + myIngredients.setNNInteraction(9,4,10.0); + + //prepare myIngredients + myIngredients.setBoxX(32); + myIngredients.setBoxY(32); + myIngredients.setBoxZ(32); + myIngredients.setPeriodicX(1); + myIngredients.setPeriodicY(1); + myIngredients.setPeriodicZ(1); + + typename Ing1::molecules_type& molecules1=myIngredients.modifyMolecules(); + molecules1.resize(2); + molecules1[0].setAllCoordinates(9,10,10); + molecules1[1].setAllCoordinates(1,1,1); + molecules1[0].setAttributeTag(1); + molecules1[1].setAttributeTag(2); + + myIngredients.synchronize(myIngredients); + + AnalyzerWriteBfmFile outfile("./interactionRW.test",myIngredients,AnalyzerWriteBfmFile::NEWFILE); + outfile.initialize(); + outfile.execute(); + + Ing1 myIngredients2; + UpdaterReadBfmFile infile("./interactionRW.test",myIngredients2,UpdaterReadBfmFile::READ_LAST_CONFIG_FAST); + infile.initialize(); + + for(int32_t i=1;i<=255;i++) + { + for(int32_t j=1;j<=255;j++) + { + EXPECT_DOUBLE_EQ(myIngredients.getNNInteraction(i,j),myIngredients2.getNNInteraction(i,j)); + } + } + + outfile.closeFile(); + infile.closeFile(); + //remove the temporary file + EXPECT_EQ(0,remove("./interactionRW.test")); + + +} + + +TEST_F(NNInteractionScTest,ApplyMoveAddMonomerSc) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureNNInteractionSc) Features1; + typedef ConfigureSystem Config1; + typedef Ingredients Ing1; + Ing1 myIngredients1; + + //prepare myIngredients1 + myIngredients1.setBoxX(65); + myIngredients1.setBoxY(65); + myIngredients1.setBoxZ(65); + myIngredients1.setPeriodicX(1); + myIngredients1.setPeriodicY(1); + myIngredients1.setPeriodicZ(1); + + myIngredients1.synchronize(); + + typename Ing1::molecules_type& molecules1=myIngredients1.modifyMolecules(); + + MoveAddMonomerSc addMonomer; + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(256); + EXPECT_THROW(addMonomer.apply(myIngredients1),std::runtime_error); + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(0); + EXPECT_THROW(addMonomer.apply(myIngredients1),std::runtime_error); + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(-1); + EXPECT_THROW(addMonomer.apply(myIngredients1),std::runtime_error); + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(2); + addMonomer.apply(myIngredients1); + + EXPECT_EQ(2,myIngredients1.getLatticeEntry(5,6,7)); + EXPECT_EQ(2,myIngredients1.getLatticeEntry(6,6,7)); + EXPECT_EQ(2,myIngredients1.getLatticeEntry(5,7,7)); + EXPECT_EQ(2,myIngredients1.getLatticeEntry(6,7,7)); + EXPECT_EQ(2,myIngredients1.getLatticeEntry(5,6,8)); + EXPECT_EQ(2,myIngredients1.getLatticeEntry(6,6,8)); + EXPECT_EQ(2,myIngredients1.getLatticeEntry(5,7,8)); + EXPECT_EQ(2,myIngredients1.getLatticeEntry(6,7,8)); + + addMonomer.init(myIngredients1); + addMonomer.setPosition(5,6,7); + addMonomer.setTag(255); + addMonomer.apply(myIngredients1); + + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(5,6,7))); + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(6,6,7))); + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(5,7,7))); + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(6,7,7))); + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(5,6,8))); + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(6,6,8))); + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(5,7,8))); + EXPECT_EQ(255,int32_t(myIngredients1.getLatticeEntry(6,7,8))); + + +} diff --git a/tests/feature/TestFeatureWall.cpp b/tests/feature/TestFeatureWall.cpp new file mode 100644 index 0000000..933d676 --- /dev/null +++ b/tests/feature/TestFeatureWall.cpp @@ -0,0 +1,316 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ + +#include +#include + +#include "gtest/gtest.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +class TestFeatureWall: public ::testing::Test{ +public: + //redirect cout output + virtual void SetUp(){ + originalBuffer=std::cout.rdbuf(); + std::cout.rdbuf(tempStream.rdbuf()); + }; + + //restore original output + virtual void TearDown(){ + std::cout.rdbuf(originalBuffer); + }; + +private: + std::streambuf* originalBuffer; + std::ostringstream tempStream; +}; + + + + +TEST(TestFeatureWall,Moves) +{ + + typedef LOKI_TYPELIST_2(FeatureMoleculesIO,FeatureWall) Features; + + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + IngredientsType ingredients; + + //prepare ingredients + ingredients.setBoxX(32); + ingredients.setBoxY(32); + ingredients.setBoxZ(32); + ingredients.setPeriodicX(1); + ingredients.setPeriodicY(1); + ingredients.setPeriodicZ(1); + + //one move of every type + MoveLocalSc scmove; + MoveAddMonomerSc addmove; + + ingredients.modifyMolecules().resize(3); + ingredients.modifyMolecules()[0].setAllCoordinates(0,0,0); + ingredients.modifyMolecules()[1].setAllCoordinates(2,0,0); + ingredients.modifyMolecules()[2].setAllCoordinates(2,2,0); + + EXPECT_NO_THROW(ingredients.synchronize(ingredients)); + + + + //************ scmove *************** + + //build wall in y-z-plane at x=2 -> x=1 should be forbidden + Wall wall1; + wall1.setBase(2,0,0); + // check forbidden nomal vector: must be a unit vector! P(1,0,0) + EXPECT_ANY_THROW(wall1.setNormal(2,0,0)); + EXPECT_ANY_THROW(wall1.setNormal(0,1,1)); + + //add correct normal vector + wall1.setNormal(1,0,0); + ingredients.addWall(wall1); +/* + std::cout << "getWalls(): " << ingredients.getWalls()[0].getBase().getX() << " " << ingredients.getWalls()[0].getBase().getY() << " " << ingredients.getWalls()[0].getBase().getZ() << " " << ingredients.getWalls()[0].getNormal().getX() << " " << ingredients.getWalls()[0].getNormal().getY() << " " << ingredients.getWalls()[0].getNormal().getZ() << std::endl; +*/ + + //check monomer 0 movement to backside of wall + while((scmove.getDir().getX()!=1) || (scmove.getIndex()!=0)) scmove.init(ingredients); //init(ingredients) searches randomly one monomer in ingredients and one direction for the move -> while you have not the right monomer and the prefered direction, keep on searching randomly + EXPECT_FALSE(ingredients.checkMove(ingredients,scmove)); + + //check monomer 1 movement to frontside of wall + while((scmove.getDir().getX()!=-1) || (scmove.getIndex()!=1)) scmove.init(ingredients); + EXPECT_FALSE(ingredients.checkMove(ingredients,scmove)); + + //check monomer 2 movement to frontside of wall + while((scmove.getDir().getX()!=-1) || (scmove.getIndex()!=2)) scmove.init(ingredients); + EXPECT_FALSE(ingredients.checkMove(ingredients,scmove)); + + //check possible monomer 0 movements + while((scmove.getDir().getX()==1) || (scmove.getIndex()!=0)) scmove.init(ingredients); //search randomly for any direction except direction to the wall (x=1) for monomer 0 + EXPECT_TRUE(ingredients.checkMove(ingredients,scmove)); + + //check possible monomer 1 movements + while((scmove.getDir().getX()==-1) || (scmove.getIndex()!=1)) scmove.init(ingredients); + EXPECT_TRUE(ingredients.checkMove(ingredients,scmove)); + + //check possible monomer 2 movements + while((scmove.getDir().getX()==-1) || (scmove.getIndex()!=2)) scmove.init(ingredients); + EXPECT_TRUE(ingredients.checkMove(ingredients,scmove)); + + //build additional wall in x-z-plane at y=2 -> x=1 and y=1 should be forbidden now + Wall wall2; + wall2.setBase(0,2,0); + wall2.setNormal(0,1,0); + ingredients.addWall(wall2); +/* + for(size_t i=0; i< ingredients.getWalls().size(); i++){ + std::cout << "getWalls()[" << i << "]: " << ingredients.getWalls()[i].getBase().getX() << " " << ingredients.getWalls()[i].getBase().getY() << " " << ingredients.getWalls()[i].getBase().getZ() << " " << ingredients.getWalls()[i].getNormal().getX() << " " << ingredients.getWalls()[i].getNormal().getY() << " " << ingredients.getWalls()[i].getNormal().getZ() << std::endl; + } +*/ + //check monomer 0 movement to backside of both walls + while((scmove.getDir().getY()!=1) || (scmove.getIndex()!=0)) scmove.init(ingredients); + EXPECT_FALSE(ingredients.checkMove(ingredients,scmove)); + + //check monomer 1 movement to frontside of wall1 and backside of wall2 + while((scmove.getDir().getY()!=1) || (scmove.getIndex()!=1)) scmove.init(ingredients); + EXPECT_FALSE(ingredients.checkMove(ingredients,scmove)); + + //check monomer 2 movement to frontside of both walls + while((scmove.getDir().getY()!=-1) || (scmove.getIndex()!=2)) scmove.init(ingredients); + EXPECT_FALSE(ingredients.checkMove(ingredients,scmove)); + + //check possible monomer 0 movements + while((scmove.getDir().getX()==1) || (scmove.getDir().getY()==1) || (scmove.getIndex()!=0)) scmove.init(ingredients); + EXPECT_TRUE(ingredients.checkMove(ingredients,scmove)); + + //check possible monomer 1 movements + while((scmove.getDir().getX()==-1) || (scmove.getDir().getY()==1) || (scmove.getIndex()!=1)) scmove.init(ingredients); + EXPECT_TRUE(ingredients.checkMove(ingredients,scmove)); + + //check possible monomer 2 movements + while((scmove.getDir().getX()==-1) || (scmove.getDir().getY()==-1) || (scmove.getIndex()!=2)) scmove.init(ingredients); + EXPECT_TRUE(ingredients.checkMove(ingredients,scmove)); + + + + ingredients.clearAllWalls(); + + //add former wall1 again + ingredients.addWall(wall1); + + //build second wall in y-z-plane at x=4 -> x=3 should be forbidden + Wall wall3; + wall3.setBase(4,0,0); + wall3.setNormal(1,0,0); + ingredients.addWall(wall3); + + //check monomer 1 movement to backside of wall3 (wall1 was already checked above) + while((scmove.getDir().getX()!=1) || (scmove.getIndex()!=1)) scmove.init(ingredients); + EXPECT_FALSE(ingredients.checkMove(ingredients,scmove)); + + //check monomer 2 movement to backside of wall3 + while((scmove.getDir().getX()!=1) || (scmove.getIndex()!=2)) scmove.init(ingredients); + EXPECT_FALSE(ingredients.checkMove(ingredients,scmove)); + + //check possible monomer 1 movements + while((scmove.getDir().getX()==-1) || (scmove.getDir().getX()==1) || (scmove.getIndex()!=1)) scmove.init(ingredients); + EXPECT_TRUE(ingredients.checkMove(ingredients,scmove)); + + //check possible monomer 2 movements + while((scmove.getDir().getX()==-1) || (scmove.getDir().getX()==1) || (scmove.getIndex()!=2)) scmove.init(ingredients); + EXPECT_TRUE(ingredients.checkMove(ingredients,scmove)); + + + + + //************ addmove *************** + + ingredients.clearAllWalls(); + + //add former walls again + ingredients.addWall(wall1); + ingredients.addWall(wall2); + ingredients.addWall(wall3); + + //create random system and check it by synchronize + std::cout << "create random config by addmove and check by snychronize"< Config; + typedef Ingredients IngredientsType; + IngredientsType ingredients1; //to write bfm-file + + typedef LOKI_TYPELIST_2(FeatureMoleculesIO,FeatureWall) Features; + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + IngredientsType ingredients2; //to read bfm-file + + //prepare ingredients1 to write bfm-file + ingredients1.setBoxX(32); + ingredients1.setBoxY(32); + ingredients1.setBoxZ(32); + ingredients1.setPeriodicX(1); + ingredients1.setPeriodicY(1); + ingredients1.setPeriodicZ(1); + + Wall wall4; + wall4.setBase(5,0,0); + wall4.setNormal(1,0,0); + ingredients1.addWall(wall4); + + Wall wall5; + wall5.setBase(0,4,0); + wall5.setNormal(0,1,0); + ingredients1.addWall(wall5); + + EXPECT_NO_THROW(ingredients1.synchronize()); + + TaskManager taskmanager; + taskmanager.addAnalyzer(new AnalyzerWriteBfmFile("test_wall.bfm",ingredients1,AnalyzerWriteBfmFile::NEWFILE)); + + taskmanager.initialize(); + taskmanager.run(); + taskmanager.cleanup(); + + TaskManager taskmanager1; + taskmanager1.addUpdater(new UpdaterReadBfmFile("test_wall.bfm",ingredients2,UpdaterReadBfmFile::READ_LAST_CONFIG_SAVE)); + + taskmanager1.initialize(); + taskmanager1.run(); + taskmanager1.cleanup(); + + remove("test_wall.bfm"); + + EXPECT_TRUE(ingredients1.getWalls()[0].getBase().getX()==ingredients2.getWalls()[0].getBase().getX()); + EXPECT_TRUE(ingredients1.getWalls()[1].getBase().getX()==ingredients2.getWalls()[1].getBase().getX()); + + EXPECT_NO_THROW(ingredients1.synchronize()); + EXPECT_NO_THROW(ingredients2.synchronize()); + +} diff --git a/tests/updater/TestUpdaterSetupLinearChains.cpp b/tests/updater/TestUpdaterSetupLinearChains.cpp new file mode 100644 index 0000000..ecafd35 --- /dev/null +++ b/tests/updater/TestUpdaterSetupLinearChains.cpp @@ -0,0 +1,130 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ + +/*****************************************************************************/ +/** + * @file + * @brief Tests for UpdaterAddLinearChains + * */ +/*****************************************************************************/ + +#include "gtest/gtest.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace std; + + + + +/************************************************************************/ +//define test fixtures for the different tests their purpose is to set up +//the tests to suppress cout's output such that is does not display on the +//standard output during the tests. this makes google test's output more readeable +/************************************************************************/ + +class TestUpdaterAddLinearChains: public ::testing::Test{ +public: + + //redirect cout output + virtual void SetUp(){ + originalBuffer=cout.rdbuf(); + cout.rdbuf(tempStream.rdbuf()); + }; + + //restore original output + virtual void TearDown(){ + cout.rdbuf(originalBuffer); + }; + +private: + std::streambuf* originalBuffer; + std::ostringstream tempStream; + +}; + +TEST_F(TestUpdaterAddLinearChains, Constructor) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureExcludedVolumeSc< FeatureLatticePowerOfTwo >,FeatureAttributes) Features; + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + + IngredientsType ingredients; + + //Constructor call + EXPECT_NO_THROW(UpdaterAddLinearChains Timmy(ingredients,1,2)); +} + +TEST_F(TestUpdaterAddLinearChains, TestUpdater) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureExcludedVolumeSc< FeatureLatticePowerOfTwo >,FeatureAttributes) Features; + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + + IngredientsType ingredients; + ingredients.setBoxX(16); + ingredients.setBoxY(16); + ingredients.setBoxZ(16); + ingredients.setPeriodicX(true); + ingredients.setPeriodicY(true); + ingredients.setPeriodicZ(true); + ingredients.modifyBondset().addBFMclassicBondset(); + EXPECT_NO_THROW(ingredients.synchronize()); + + UpdaterAddLinearChains Tommy(ingredients, 4, 16); + + // first execution + EXPECT_TRUE(Tommy.execute()); + EXPECT_NO_THROW(ingredients.synchronize()); + EXPECT_EQ((4*16),ingredients.getMolecules().size()); + // check (default) monomer taggs + EXPECT_EQ(1,ingredients.getMolecules()[0].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[1].getAttributeTag()); + EXPECT_EQ(1,ingredients.getMolecules()[2].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[15].getAttributeTag()); + + // check connectivity + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(0)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(1)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(14)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(15)); + + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(16)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(17)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(30)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(31)); + +} +