From 67151b03e8cd26e7084e6742e2d2d703caff180e Mon Sep 17 00:00:00 2001 From: Ankus7890 Date: Wed, 21 Aug 2019 22:49:59 +0200 Subject: [PATCH 1/9] lipidbilayerrelated --- .../analyzer/AnalyzerPermeabilityCalculator.h | 598 ++++++++++++++ .../LeMonADE/analyzer/AnalyzerPoreFinder.h | 682 +++++++++++++++ .../feature/FeatureBendingPotential.h | 515 ++++++++++++ .../FeatureBendingPotentialReadWrite.h | 145 ++++ .../LeMonADE/updater/UpdaterLipidsCreator.h | 780 ++++++++++++++++++ .../updater/UpdaterLipidsRandomInitializer.h | 251 ++++++ .../TestAnalyzerPermeabilityCalculator.cpp | 373 +++++++++ tests/analyzer/TestAnalyzerPoreFinder.cpp | 459 +++++++++++ tests/feature/TestFeatureBendingPotential.cpp | 408 +++++++++ tests/updater/TestUpdaterLipidsCreator.cpp | 358 ++++++++ .../TestUpdaterLipidsRandomInitializer.cpp | 270 ++++++ tests/updater/TesterHelperFunctionsLipids.h | 76 ++ 12 files changed, 4915 insertions(+) create mode 100644 include/LeMonADE/analyzer/AnalyzerPermeabilityCalculator.h create mode 100644 include/LeMonADE/analyzer/AnalyzerPoreFinder.h create mode 100644 include/LeMonADE/feature/FeatureBendingPotential.h create mode 100644 include/LeMonADE/feature/FeatureBendingPotentialReadWrite.h create mode 100644 include/LeMonADE/updater/UpdaterLipidsCreator.h create mode 100644 include/LeMonADE/updater/UpdaterLipidsRandomInitializer.h create mode 100644 tests/analyzer/TestAnalyzerPermeabilityCalculator.cpp create mode 100644 tests/analyzer/TestAnalyzerPoreFinder.cpp create mode 100644 tests/feature/TestFeatureBendingPotential.cpp create mode 100644 tests/updater/TestUpdaterLipidsCreator.cpp create mode 100644 tests/updater/TestUpdaterLipidsRandomInitializer.cpp create mode 100644 tests/updater/TesterHelperFunctionsLipids.h diff --git a/include/LeMonADE/analyzer/AnalyzerPermeabilityCalculator.h b/include/LeMonADE/analyzer/AnalyzerPermeabilityCalculator.h new file mode 100644 index 0000000..dcc48db --- /dev/null +++ b/include/LeMonADE/analyzer/AnalyzerPermeabilityCalculator.h @@ -0,0 +1,598 @@ +#ifndef AnalyzerPermeabilityCalculator_H +#define AnalyzerPermeabilityCalculator_H + +#include +#include + +/** + * @file + * + * @class AnalyzerPermeabilityCalculator + * + * @brief Analyzer to calculate permeability across the bilayer which has the normal along z axis on LatticePowerOfTwo. + * + * @details This analyzer calculates the average permeability across the bilayer(normal along z axis) + * of objects (like copolymers and nano particles (see UpdaterLipidsCreator)) which are setup in as + * std::vector > and solvent monomers. The analyzer calculates the solvent + * and object permeability based on bufferDumpSize choosen by user and dumps it in file with default name + * Permeability.dat.
+ * + * \b Important:
+ * (i)It is reccomended to create solvent and objects monomer with tags greater than and equal to 5(see UpdaterLipidsCreator). + * If one needs to change this number from 5 to some other tag, it can be done by changing \SOLVENTAG below. + * Thus, to use this analyzer, the general scheme for tags should be \b { Lipid monomers(heads and tails)< SOLVENTAG \b } + * and \b { Non-Lipids monomers including Solvent >= SOLVENTAG \b }.
+ * + * (ii)In the case of the bilayers with pores, it is essential to find the position of pores are in xy plane. + * The calculation of permeability is difficult in the region of pore, as there is no interface in this region. + * To solve this problem, this analyzer uses method similar to S. Podogin et. al. ACS Nano, 6:10555–10561, 2012.
+ * The bilayer plane is divided into bins and permeability is calculated in seperately in these bins. If a bin + * has less than 30% of average lipid monomer in other bins, it is considered to be a pore state. This threshold + * of 30% can be controlled by \b PORE_THRESHOLD below.
+ * + * (iii) It is reccomended to simulate the system for 100 MCS(simInterval_=100)(or track permeable objects after 100 MCS) + * before executing this analyzer to get proper results.
+ * + * @tparam IngredientsType + * + * @param ingredients_ the system holding the simulation box with bilayer, solvent and objects. + * @param groups_ subgroup of molecules_type having groups of objects which translocate through bilayer. + * + * @param binSize_ size of bin in which the bilayer plane is divided and permeability + * is calculated in each of the bins. + * @param factorSigma_ this is the factor which is multiplied to spread of lipid monomers along z-axis(standard deviation). + * This was choosen to be 3 in S. Podogin et. al. ACS Nano, 6:10555–10561, 2012 and 5 in Werner et. al., Soft Matter + * 8:11714,2012. Please see these articles for further information. + * @param outputFile_ the user defined output file name for dumping permeability values on disk. + * @param simInterval_ the interval of MCS for which system is simulated before the analyzer is executed. It basically + * tells interval(MCS) after which objects and solvent molecules are tracked. + * @param bufferDumpSize_ the buffer size after which to print the permeability values into the file. This value + * also automatically sets interval of MCS after which permeability is calculated. For example, for default permeability + * value is calculated after every simInterval_*bufferDumpSize_=10000 MCS. + * @param midplaneUpdaterInterval_ the interval of MCS after which midplane values should be updated. For example, for default + * values, it would be updated every simInterval_*midplaneUpdaterInterval_=15000MCS. + **/ + +#define SOLVENTAG 5 +#define PORE_THRESHOLD 0.3 + +using namespace std; +template +class AnalyzerPermeabilityCalculator{ + +protected: + //!Typedefs for various underlying container holding the monomers + //!and monomer subgroups. + typedef typename IngredientsType::molecules_type molecules_type; + typedef std::vector< MonomerGroup< molecules_type> > group_type; + typedef MonomerGroup< molecules_type> group_type_ind; + const IngredientsType& ingredients; + const molecules_type& molecules; + const group_type& groups; + + //!More explainantion of these variables could be found below. + int32_t binSize,factorSigma,bufferDumpSize,simInterval,midplaneUpdaterInterval; + std::vector groupsCOM; + int32_t boxXm_1,boxYm_1,boxZm_1,binBoxX,binBoxY,boxZ; + + //!Maps for storing indices of particle which are inside the boundaries(close to the bilayer). + std::map particlesInsideBoundaries; + std::map particlesInsideBoundariesOld; + + //!Vector for storing the indices of monomers which are solvent. + std::vector < int32_t > solventIndices; + + //!Vectors and variables used by midplaneUpdater() function. + std::vector > midplane; + std::vector > counterMidplane; + std::vector > poreFlag; + std::vector neigbours; + int32_t counterSolvent,counterObjects,counterExecute; + float sigma; + + //!Variables used by dumpPermeabilityPerMcs() function. + std::string outputFile; + bool isFirstFileDump; + +public: + AnalyzerPermeabilityCalculator(const IngredientsType& ingredients_,const group_type& groups_,int32_t binSize_=8, int32_t factorSigma_=3,std::string outputFile_="Permeability.dat", int32_t bufferDumpSize_=100,int32_t midplaneUpdaterInterval_=150,int32_t simInterval_=100); + + //!General updater functions + void execute(); + void initialize(){}; + void cleanup(){}; + + //!Function for updating permeability values for objects and solvent monomers. + void permeabilityUpdate(); + + //!Function for updating position of midplane and width of bilayer. + void midplaneUpdater(); + + //!Function used by midplaneUpdater() to find the nearest neigbouring bin's midplane value. + void findNearestMidplaneValue(int32_t x, int32_t y); + + //!Function to initiate the array with initial maps of particles inside/outside boundaries. + void initParticleArrays(); + + //!Function to find distances using minimum image convention. + int32_t reduceDistanceInPeriodicSpace(int32_t distance, int32_t period); + + //!Function to find the center of mass of subgroup of monomers(generally objects). + VectorDouble3 centerOfMass(const group_type_ind& m); + + //!Function for dumping the permeability values in the file. + void dumpPermeabilityPerMcs(); + +}; + +/** +* @brief Constructor +**/ +template +AnalyzerPermeabilityCalculator::AnalyzerPermeabilityCalculator(const IngredientsType& ingredients_,const group_type& groups_, int32_t binSize_, int32_t factorSigma_,std::string outputFile_ ,int32_t bufferDumpSize_, int32_t midplaneUpdaterInterval_, int32_t simInterval_) +:ingredients(ingredients_) +,molecules(ingredients_.getMolecules()) +,groups(groups_) +,binSize(binSize_) +,factorSigma(factorSigma_) +,outputFile(outputFile_) +,bufferDumpSize(bufferDumpSize_) +,midplaneUpdaterInterval(midplaneUpdaterInterval_) +,simInterval(simInterval_) +{ + + //The list of neigbours on 2-dimensional surface used by + //function findNearestMidplaneValue() to find a neigbour + //of a bin which doesn't have a pore. + + + neigbours.push_back(VectorInt3(-1,0,0)); + neigbours.push_back(VectorInt3(0,-1,0)); + neigbours.push_back(VectorInt3(1,0,0)); + neigbours.push_back(VectorInt3(0,1,0)); + neigbours.push_back(VectorInt3(-1,1,0)); + neigbours.push_back(VectorInt3(1,-1,0)); + neigbours.push_back(VectorInt3(1,1,0)); + neigbours.push_back(VectorInt3(-1,-1,0)); + + //The variables for number of bins in x and y direction + //depending on periodic box length in those directions + //respectively. + + binBoxX=ingredients.getBoxX()/binSize; + binBoxY=ingredients.getBoxY()/binSize; + + //Counters for various utilities below. + + counterExecute=0; + counterObjects=0; + counterSolvent=0; + + //Default values for box_length and boxXm_1 + //for folding back and minimum image convention + //utilities. + + boxZ=ingredients.getBoxZ(); + + boxXm_1=ingredients.getBoxX()-1; + boxYm_1=ingredients.getBoxY()-1; + boxZm_1=ingredients.getBoxZ()-1; + + //Initiate arrays of the size as number of bins in x and + //y direction. + + midplane.resize(binBoxX, std::vector(binBoxY, boxZ/2.0)); + counterMidplane.resize(binBoxX, std::vector(binBoxY, 0)); + poreFlag.resize(binBoxX, std::vector(binBoxY, 0)); + sigma=5.0; + isFirstFileDump=1; + + //Create an array of solvent indices so that it can be used to + //directly access the position of solvent monomer using + //ingredients::molecules_type + + for(int32_t i=0;i +void AnalyzerPermeabilityCalculator::execute(){ + + //Increase the counter. + counterExecute++; + + int32_t x,y,z,distance,f2s,sign; + + //Loop to go over all the particles(objects and solvent) + //to see if they are in the boundaries. + + for(size_t i=0;i +void AnalyzerPermeabilityCalculator::permeabilityUpdate(){ + + int32_t x,y,upDownIndicator; + + //Loop to go over all the particles(objects and solvent) + //which were present in last execution and check if they + //still present, if not then, check if they went from + //opposite direction they came in and count a translocation + //event. If a particle is still present, then store the + //direction it came inside from. + + for(int32_t i=0;i +void AnalyzerPermeabilityCalculator::midplaneUpdater(){ + + //Reset all values to zero + + for(int32_t x1=0;x1=SOLVENTAG) continue; + + x=(int32_t(ingredients.getMolecules()[i].getX())&boxXm_1)/binSize; + y=(int32_t(ingredients.getMolecules()[i].getY())&boxYm_1)/binSize; + + z=ingredients.getMolecules()[i].getZ(); + + midplane[x][y] += (reduceDistanceInPeriodicSpace(z-z_Ref,boxZ)); + counterMidplane[x][y]++; + } + + float maxValue=0; + + for(int32_t x1=0;x1counterMidplane[x1][y1])?maxValue:counterMidplane[x1][y1]; + } + + //Find pore flags. + for(int32_t x1=0;x1=SOLVENTAG) continue; + + x=(int32_t(ingredients.getMolecules()[i].getX())&boxXm_1)/binSize; + y=(int32_t(ingredients.getMolecules()[i].getY())&boxYm_1)/binSize; + z=ingredients.getMolecules()[i].getZ(); + + if(!poreFlag[x][y]){ + distance=(z-midplane[x][y]); + distance=reduceDistanceInPeriodicSpace(distance,boxZ); + sigma += distance*distance; + counterSigma++; + } + } + + sigma=sqrt(sigma/counterSigma); +}; + +/** +* @details Finds a niegbouring bin without any pore and sets a midplane of bin with given +* coordinates to that neigbouring bin. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +* +* @param x1 x coordinate of the bin with a pore. +* @param y1 y coordinate of the bin with a pore. +*/ + +template +void AnalyzerPermeabilityCalculator::findNearestMidplaneValue(int32_t x1, int32_t y1){ + + int32_t neighX,neighY=0; + + for(int nL=0;nL +void AnalyzerPermeabilityCalculator::initParticleArrays(){ + + int32_t x,y,distance,sign,f2s; + + for(size_t i=0;i +void AnalyzerPermeabilityCalculator::dumpPermeabilityPerMcs(){ + + std::vector< std::vector< double > > resultsTimeseries(3,std::vector< double >(1)); + resultsTimeseries[0][0] = ingredients.getMolecules().getAge(); + resultsTimeseries[1][0] = counterObjects/(bufferDumpSize*simInterval); + resultsTimeseries[2][0] = counterSolvent/(bufferDumpSize*simInterval); + + //if it is written for the first time, include comment in the output file + if(isFirstFileDump){ + + std::stringstream commentTimeSeries; + commentTimeSeries<<"Created by AnalyzerPermeabilityCalculator\n"; + commentTimeSeries<<"file contains time series of average permeability averaged per \n"; + commentTimeSeries<<"format: mcs\t P_objects\t P_solvent\n"; + + ResultFormattingTools::writeResultFile( + outputFile, + ingredients, + resultsTimeseries, + commentTimeSeries.str() + ); + + isFirstFileDump=false; + } + //otherwise just append the new data + else{ + ResultFormattingTools::appendToResultFile(outputFile, + resultsTimeseries); + } + + counterObjects=0; + counterSolvent=0; + + } + +/** +* @details Calculates the center of mass of objects given a subgroup of +* monomers of the object. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ + +template < class IngredientsType > +VectorDouble3 AnalyzerPermeabilityCalculator::centerOfMass(const group_type_ind& m) +{ + VectorDouble3 CoM_sum; + + for ( uint i = 0; i < m.size(); ++i) + { + CoM_sum.setX( CoM_sum.getX() + m[i].getX() ); + CoM_sum.setY( CoM_sum.getY() + m[i].getY() ); + CoM_sum.setZ( CoM_sum.getZ() + m[i].getZ() ); + } + + float inv_N = 1.0 / double ( m.size() ); + + VectorDouble3 CoM ( + float ( CoM_sum.getX() ) * inv_N, + float ( CoM_sum.getY() ) * inv_N, + float ( CoM_sum.getZ() ) * inv_N); + + return CoM; + +}; + +/** +* @details Finds the distances in the minimum image convention give a period. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ + +template +int32_t AnalyzerPermeabilityCalculator::reduceDistanceInPeriodicSpace(int32_t distance, int32_t period) +{ + if(distance > period/2) + {while(-(distance-period) period/2) + {while((distance+period) +#include + +using namespace std; + +/** + * @file AnalyzerPoreFinder.h + * + * @class VectorInt2 + * + * @brief This class is created as substitute to VectorInt3 in in class Vector3D.h, as the + * most calculations required in class AnalyzerPoreFinder only requires a two dimensional space. + **/ + + +class VectorInt2{ +public: + VectorInt2(int32_t x_=0,int32_t y_=0): + x(x_), + y(y_){} + + void setX(int32_t x_){x=x_;}; + void setY(int32_t y_){y=y_;}; + int32_t getX(){return x;}; + int32_t getY(){return y;}; + + float getLength(){return sqrt(x*x+y*y);}; + + void setAllCoordinates(int32_t x_,int32_t y_){x=x_; y=y_;}; + + VectorInt2 operator+(const VectorInt2& b) { + VectorInt2 vec; + vec.x = this->x + b.x; + vec.y = this->y + b.y; + return vec; + }; + + VectorInt2 operator/(int32_t divisor) { + VectorInt2 vec; + vec.x = round(this->x, divisor); + vec.y = round(this->y, divisor); + return vec; + }; + + friend ostream& operator<<(ostream& out, const VectorInt2& b){ + out< + * + * Also important to set all lipid monomers below SOLVENTAG and all non lipid monomers including solvent to SOLVENTAG.
+ * Please see AnalyzerPermeabilityCalculator and UpdaterLipidsCreator for further details.
+ * + * This class produces a file with time series of pore sizes, averaged position of pore points and copolymers + * around the centriod of pore in files named as "PoreSize.dat","PoreCoordinates.dat" and "PolymerCoordinates.dat" + * respectively, with an option to add suffix to file name using an option in the constructor. + * + * @tparam IngredientsType + * + * @param groups_ subgroup of molecules_type having groups of objects which translocate through bilayer. + * @param ingredients_ the system holding the simulation box with bilayer, solvent and objects. + * + * @param filesuffix_ option to add suffix to filename which are output at the end. + * + **/ + +#define CLUSTER_SIZE_THRESHOLD 20 +#define POLYMER_DIST_THRESHOLD 10 +#define SOLVENTAG 5 +#define TAILTAG 2 + +template +class AnalyzerPoreFinder: public AbstractAnalyzer +{ +protected: + //!Typedefs for various underlying container holding the monomers + //!and monomer subgroups. + typedef typename IngredientsType::molecules_type molecules_type; + typedef std::vector< MonomerGroup< molecules_type> > group_type; + IngredientsType& ingredients; + const molecules_type& molecules; + const group_type& groups; + typedef MonomerGroup< molecules_type> group_type_ind; + + //!Box variables. + int32_t boxX,boxY,boxZ,boxXm_1,boxYm_1,boxZm_1; + + //!vectors storing VectorInt2 objects + std::vector coordinatesOfPore; + std::vector centriod; + std::vector neigbours; + std::vector clusterList; + + //!vectors various arrays used by functions below. + + std::vector > tailOccupied; + std::vector > coordinatesOfPoreStore; + std::vector > coordinatesOfPolymers; + std::vector sizeOfCluster; + std::vector< std::vector< int32_t > > poreSizeStore; + + //!string variables for name of output files. + std::string outputFilePoreSize,outputFilePore,outputFilePolymer; + + int32_t counterForPore,counterForPolymer; + int32_t binforradius; + int32_t referenceI; + +public: + AnalyzerPoreFinder(IngredientsType& i, const group_type& g,std::string fileSuffix_=""); + virtual ~AnalyzerPoreFinder(){} + + //!Function to fill lattice with tail monomer types. + void LatticeFiller(); + + //!Function to store lattice points which are occupied by tail monomers. + void tailOccupiedFinder(); + + //!Function for cluster analysis of points free of any tail monomers. + void clusterAnalysis(); + + //!Function used by clusterAnalysis() to find a lipid site which is free of lipids and not been checked already. + VectorInt2 lookForLipidFreeSites(); + + //!Function used by clusterAnalysis() addNeigbours of an last unchecked values in the cluster list. + void addNeigbours(int32_t lastUnChecked); + + //!Function used by clusterAnalysis() to store various values to the memory for later dump. + void storeValues(); + + //!Function finding the polymer coordinates wrt the centriod of pore. + void polymerCoordinatesFinder(); + + //!Function to find the center of mass of subgroup of monomers(generally objects). + VectorDouble3 centerOfMass(const group_type_ind& m); + + //!Function to find distances using minimum image convention. + int32_t reduceDistanceInPeriodicSpace( int32_t distance, int32_t period); + + //!General updater functions + bool execute(); + void cleanup(); + void initialize(){}; + +}; + +/** +* @brief Constructor +**/ +template +AnalyzerPoreFinder::AnalyzerPoreFinder(IngredientsType& i, const group_type& g, std::string fileSuffix_) +:ingredients(i) +,molecules(i.getMolecules()) +,groups(g) +{ + //Initate the box variables and counter for execute function + // and polymers as zero. + + counterForPore=0; + counterForPolymer=0; + boxX=ingredients.getBoxX(); + boxY=ingredients.getBoxY(); + boxZ=ingredients.getBoxZ(); + boxXm_1=ingredients.getBoxX()-1; + boxYm_1=ingredients.getBoxY()-1; + boxZm_1=ingredients.getBoxZ()-1; + + //Add suffix if needed to added for output files. + + outputFilePoreSize="PoreSize"+fileSuffix_+".dat"; + outputFilePore="PoreCoordinates"+fileSuffix_+".dat"; + outputFilePolymer="PolymerCoordinates"+fileSuffix_+".dat"; + + //Finding the index of the lipid monomer on the bilayer for reference to + //calculate value of midplane. This is mainly done for case when bilayer is + //at the of periodic box in z direction(z=0 or z=boxZ-1). + + for(int32_t i=0;i(boxY)); + coordinatesOfPoreStore.resize(boxX,std::vector(boxY,0)); + coordinatesOfPolymers.resize(boxX,std::vector(boxY,0)); + + //Initate pore size array with three rows. + //for MCS, poresize, and third one radius of gyration. + poreSizeStore.resize(3); + + //Initate nearest neigbours on 2d surface on simple cubic lattice. + VectorInt2 V(1,0); + neigbours.push_back(V); + V.setAllCoordinates(-1,0); + neigbours.push_back(V); + V.setAllCoordinates(0,1); + neigbours.push_back(V); + V.setAllCoordinates(0,-1); + neigbours.push_back(V); +} + +/** +* @details Execution of analyzer where it finds the point filled with tail monomers, +* then analyses the cluster of points free of tail monomers. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ + +template +bool AnalyzerPoreFinder::execute(){ + + //Fill the lattice + //and find all the sites filled with tail monomers in x-y plane. + LatticeFiller(); + tailOccupiedFinder(); + + //Perform cluster analysis. + clusterAnalysis(); + + //If no pore found, return true as there is no need to find polymer coordinates + // in closed pore. + if(coordinatesOfPore.size()==0) + return true; + + //If there is a pore find where centerOfMass of copolymers + //lie wrt centriod of the pore + polymerCoordinatesFinder(); + + return true; +} + +/** +* @details Fill lattice if there is a tail monomer with its type or zero otherwise +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ + +template +void AnalyzerPoreFinder::LatticeFiller(){ + + for(int32_t x=0;x +void AnalyzerPoreFinder::tailOccupiedFinder(){ + + for(int32_t x=0;x +void AnalyzerPoreFinder::clusterAnalysis(){ + //Restart the arrays with zero size + coordinatesOfPore.resize(0); + centriod.resize(0); + sizeOfCluster.resize(0); + sizeOfCluster.push_back(0); + + + // Outside do-while loop looks for points on lattice in x-y plane + // which are free of lipids, once one such point is found, the inside + // loop finds the cluster associated with this point. At the end, all + // the clusters with the size bigger than CLUSTER_SIZE_THRESHOLD are saved + // in coordinatesofPore for further processing. + + do{ + //Restart the clusterList + clusterList.resize(0); + //initial coordinate of point which is tail free + VectorInt2 iniCoordinate=lookForLipidFreeSites(); + + //If lookForLipidFreeSites returns an impossible coordinate + //this means all the lattice points have been looked upon. + // thus break the loop. + + if(iniCoordinate.getX()==(boxX+1)) break; + + //Otherwise add this clusterList + clusterList.push_back(iniCoordinate); + int32_t lastUnChecked=0; + + do{ + //check the for neigbours of the + //lastUnChecked cluster points which are tail free + addNeigbours(lastUnChecked); + //increase the counter of unchecked lattice sites. + lastUnChecked++; + //If the lastUnChecked has yet gone through + // all the clusterList, keep do-ing + }while(lastUnChecked!=clusterList.size()); + + //if clusterlist size is larger than threshold then save. + if(clusterList.size()>CLUSTER_SIZE_THRESHOLD){ + sizeOfCluster.push_back(clusterList.size()+coordinatesOfPore.size()); + coordinatesOfPore.insert(coordinatesOfPore.end(), clusterList.begin(), clusterList.end()); + } + + }while(true); + + //Once the cluster is found, + // store the values. + storeValues(); +} + +/** +* @details Adds neigbours of a lattice point which are free of tails. +* This function also sets tail occupied tag to true as an indicator +* that the lattice site has been checked before. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ + +template +void AnalyzerPoreFinder::addNeigbours(int32_t lastUnChecked){ + + VectorInt2 neigbourToPoint; + int32_t x,y; + VectorInt2 pointTobeChecked=clusterList[lastUnChecked]; + + for(int32_t i=0;i +VectorInt2 AnalyzerPoreFinder::lookForLipidFreeSites(){ + + for(int32_t x=0;x +void AnalyzerPoreFinder::storeValues(){ + + //Centriods calculation + //Similar to midplane calculation here we need a reference point on + //x-y plane to cover for boundary cases. + + for(int32_t i=1;i +void AnalyzerPoreFinder::polymerCoordinatesFinder(){ + /**find the midplane */ + int32_t midplane=0,lipidcounter=0; + for(int32_t i=0;i=SOLVENTAG) continue; + midplane += (reduceDistanceInPeriodicSpace(ingredients.getMolecules()[i].getZ()-ingredients.getMolecules()[referenceI].getZ(),boxZ)); + lipidcounter++; + } + + midplane /= lipidcounter; + midplane = (midplane+ingredients.getMolecules()[referenceI].getZ())&boxZm_1; + //store polymer COM wrt the centriod. + for(int32_t CN=0;CN +VectorDouble3 AnalyzerPoreFinder::centerOfMass(const group_type_ind& m) +{ + VectorDouble3 CoM_sum; + + for ( uint i = 0; i < m.size(); ++i) + { + CoM_sum.setX( CoM_sum.getX() + m[i].getX() ); + CoM_sum.setY( CoM_sum.getY() + m[i].getY() ); + CoM_sum.setZ( CoM_sum.getZ() + m[i].getZ() ); + } + + float inv_N = 1.0 / double ( m.size() ); + + VectorDouble3 CoM ( + float ( CoM_sum.getX() ) * inv_N, + float ( CoM_sum.getY() ) * inv_N, + float ( CoM_sum.getZ() ) * inv_N); + + return CoM; + +}; + +/** +* @details Finds the distances in the minimum image convention give a period. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ + + +template +int32_t AnalyzerPoreFinder::reduceDistanceInPeriodicSpace(int32_t distance, int32_t period) +{ + + if(distance > period/2){ + while(-(distance-period) period/2){ + while((distance+period) +void AnalyzerPoreFinder::cleanup(){ + + std::stringstream comments; + comments<<"Created by AnalyzerPoreFinder\n"; + comments<<"file contains time series of weighted average pore size and radius of gyration of pore(second cumulant of average pore size) \n"; + comments<<"format: mcs\t PoreSize\t Radius_Gyration\n"; + + ResultFormattingTools::writeResultFile( + outputFilePoreSize, + ingredients, + poreSizeStore, + comments.str() + ); + + comments.str(""); + comments<<"Created by AnalyzerPoreFinder\n"; + comments<<"file contains averaged pore points from centriod of a pore\n"; + comments<<"format: x\t y\t Average_value\n"; + std::vector< std::vector< float > > coordinatesOfPoreFormat; + coordinatesOfPoreFormat.resize(3); + + for(int32_t x=0;x > coordinatesOfPolyFormat; + coordinatesOfPolyFormat.resize(3); + for(int32_t x=0;x. + +--------------------------------------------------------------------------------*/ + + +#ifndef FEATURE_BENDING_POTENTIAL_H +#define FEATURE_BENDING_POTENTIAL_H + +/** + * @file + * @date 2019/08/11 + * @author Ankush Checkervarty + * @brief Definition and implementation of class template FeatureBendingPotential +**/ + +#include +#include +#include +#include +#include +#include + + + +namespace BendingPotentials +{ + /*This namespace is created so different functional form could be created or used + *for kind of bending potential which is supposed to be implemented. Here for example, + *a simple harmonic potential has been shown. + * Functional form:
+ * + * U(\theta)=\theta**2. + * + */ + + float simpleHarmonic(VectorInt3& bondVec1, VectorInt3& bondVec2){ + + float magnitude,cosine; + + magnitude =bondVec1*bondVec1; + magnitude *= bondVec2*bondVec2; + magnitude =sqrt(magnitude); + + cosine= (bondVec1*bondVec2)/magnitude; + + return acos(cosine)*acos(cosine); + } +} + + +/** + * @class FeatureBendingPotential + * @brief Provides implementation of bending potential on the BFM polymer chains. + * + * @details + * Different potential strength can be set for a bonded monomer type by using + * setBendingPotential(type,energy). This feature only works for a linear chain + * of a monomer type, thus, the monomer to be moved could be attached to either + * one or two of the same monomer type. + * The linear chain could bonded to chains other monomer type. For example, + * in a lipid(see UpdaterLipidsCreator). Even though tail monomers are of different + * type than head monomer and head monomers are bonded to tail monomers. One can + * set different (or None at default) bending potential strength for head and tail monomers. + * \b IMPORTANT:
+ * (i) It is important to create solvent(or non bonded monomers) after the bonded objects. + * This reduces memory requirements. This feature throws runtime_error if opposite is done.
+ * (ii)The type can be any integer between 1 and 10.
+ * The feature adds the bfm-file command !bending_potential A energy + * for monomers of type A with bending potential of E in kT + **/ + +class FeatureBendingPotential:public Feature +{ + +private: + + //! total number of non solvent monomers in the system + int32_t numNonSolvent; + + //! bending potential strength of the type of the monomers. type range [1:10] + double bpStrengthTable[11]; + + //! Lookup table for exp(-bpStrengthTable[a]) + double probabilityLookup[11][180][180]; + + //! Lookup to tag chains ends and their neigbours. + std::vector chainsEnds; + + //! Returns this feature's factor for the acceptance probability for the given Monte Carlo move and type. + template + double calculateAcceptanceProbability(const IngredientsType& ingredients, + const MoveLocalSc& move, + int32_t monoType) const; + + //! Used to fill in the values in chainsEnds. + void tagNiegbsandEnds(int32_t index, int32_t sameAttNeigbIndex); + +public: + + FeatureBendingPotential(); + ~FeatureBendingPotential(){} + + + //This feature adds interaction energies, so it requires FeatureBoltzmann + typedef LOKI_TYPELIST_1(FeatureBoltzmann) required_features_back; + + //FeatureAttributes needs to be in front to decide which type will + //have which bending potential strength. + //FeatureBondset needs to be in front. As it required to fill probabilityLookup table. + typedef LOKI_TYPELIST_2( + FeatureAttributes, + FeatureBondset<>) + 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; + + //! This is simple function to convert bond vector to integer ID in one to one mapping. + uint32_t bondVectorToIndex(const VectorInt3& bondVector) const; + + //!adds bending potential strength to type. + template + void setBendingPotential(IngredientsType& ingredients, int32_t type,double energy); + + //!sets up the chain ends tags in chainsEnds. + template + void setChainEnds(IngredientsType& ingredients); + + //!returns the bending potential strength for type. + double getBendingPotential(int32_t type) const; + + //!export bfm-file read command !bending_potential + template + void exportRead(FileImport & fileReader); + + //!export bfm-file write command !bending_potential + template + void exportWrite(AnalyzerWriteBfmFile & fileWriter) const; + +}; + + +////////////////// IMPLEMENTATION OF MEMBERS ////////////////////////////////// + +/** + * @brief Constructor + **/ + +FeatureBendingPotential::FeatureBendingPotential() +{ + //initialize the bpStrengthTable and probability lookups with default values + for(size_t n=0;n<11;n++) + { + bpStrengthTable[n]=0; + for(size_t m=0;m<180;m++) + for(size_t o=0;o<180;o++) + probabilityLookup[n][m][o]=1.0; + } +} + +/** + * @details calculates the factor for the acceptance probability of the move + * arising from the bending potential 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 +bool FeatureBendingPotential::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. + + //if monoType has zero strength then returns true without making any change in + //probability + int32_t monoType=ingredients.getMolecules()[move.getIndex()].getAttributeTag(); + if(!bpStrengthTable[monoType]) return true; + + double prob=calculateAcceptanceProbability(ingredients,move,monoType); + move.multiplyProbability(prob); + return true; +} + + +/** + * @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: + * - !bending_potential + * . + * + * @tparam IngredientsType The type of the system including all features + * @param fileReader File importer for the bfm-file + **/ +template +void FeatureBendingPotential::exportRead(FileImport< IngredientsType >& fileReader) +{ +// typedef FeatureBendingPotential my_type; + fileReader.registerRead("!bending_potential",new ReadBendingPotential(fileReader.getDestination())); +} + + +/** + * 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: + * - !bending_potential + * + * @tparam IngredientsType The type of the system including all features + * @param fileWriter File writer for the bfm-file. + */ +template +void FeatureBendingPotential::exportWrite(AnalyzerWriteBfmFile< IngredientsType >& fileWriter) const +{ +// typedef FeatureBendingPotential my_type; + fileWriter.registerWrite("!bending_potential",new WriteBendingPotential(fileWriter.getIngredients_())); +} + + +/** + * @details The function calculates the factor for the acceptance probability + * for the local move given as argument. The calculation is based on whether the + * choosen monomers is at middle of the chain, last monomer away from start of + * the chain monomer, one monomer away from end of the chain monomer or any end monomers. + * + * @tparam IngredientsType The type of the system including all features. + * @param [in] ingredients A reference to the IngredientsType - mainly the system. + * @param [in] move reference to the local move for which the calculation is performed. + * @param [in] monoType attribute of the monomer moved. + * + * @return acceptance probability factor for the move arising from bending potential interactions. + **/ + +template +double FeatureBendingPotential::calculateAcceptanceProbability( + const IngredientsType& ingredients, + const MoveLocalSc& move, + int32_t monoType) const +{ + + int32_t index=move.getIndex(); + + VectorInt3 presentPos=ingredients.getMolecules()[move.getIndex()]; + VectorInt3 direction=move.getDir(); + VectorInt3 bondvector1,bondvector2; + + double prob=1.0; + + VectorInt3 futurePos = presentPos + direction; + /**distance from the start of the chain and end of the chain*/ + int32_t distfromSOC=3&chainsEnds[index]; + int32_t distfromEOC=(3<<2)&chainsEnds[index]; + distfromEOC = distfromEOC>>2; + + float probBeforeMove,probAfterMove; + + /**Generally, a move effects all the bending angles of two monomer in postive + * direction and two monomer negative direction to moved monomer. However, special + * care has to be taken if the monomer is end(or near) of the chain or start of chain.*/ + + //check if it is the end monomer. + if(distfromSOC&&distfromEOC){ + VectorInt3 presentPosm_1=ingredients.getMolecules()[index-1]; + VectorInt3 presentPosp_1=ingredients.getMolecules()[index+1]; + + bondvector1=presentPosp_1-futurePos; + bondvector2=futurePos-presentPosm_1; + + probAfterMove=probabilityLookup[monoType][bondVectorToIndex(bondvector1)][bondVectorToIndex(bondvector2)]; + + bondvector1=presentPosp_1-presentPos; + bondvector2=presentPos-presentPosm_1; + + probBeforeMove=probabilityLookup[monoType][bondVectorToIndex(bondvector1)][bondVectorToIndex(bondvector2)]; + + prob *= probAfterMove/probBeforeMove; + } + + //check end monomer lies in negative direction. + if(distfromSOC>1){ + + VectorInt3 presentPosm_1=ingredients.getMolecules()[index-1]; + VectorInt3 presentPosm_2=ingredients.getMolecules()[index-2]; + + bondvector1=futurePos-presentPosm_1; + bondvector2=presentPosm_1-presentPosm_2; + + probAfterMove=probabilityLookup[monoType][bondVectorToIndex(bondvector1)][bondVectorToIndex(bondvector2)]; + + bondvector1=presentPos-presentPosm_1; + + probBeforeMove=probabilityLookup[monoType][bondVectorToIndex(bondvector1)][bondVectorToIndex(bondvector2)]; + + prob *= probAfterMove/probBeforeMove; +} + + //check end monomer lies in positive direction. + if(distfromEOC>1){ + + VectorInt3 presentPosp_1=ingredients.getMolecules()[index+1]; + VectorInt3 presentPosp_2=ingredients.getMolecules()[index+2]; + + bondvector1=presentPosp_2-presentPosp_1; + bondvector2=presentPosp_1-futurePos; + + probAfterMove=probabilityLookup[monoType][bondVectorToIndex(bondvector1)][bondVectorToIndex(bondvector2)]; + + bondvector2=presentPosp_1-presentPos; + + probBeforeMove=probabilityLookup[monoType][bondVectorToIndex(bondvector1)][bondVectorToIndex(bondvector2)]; + + prob *= probAfterMove/probBeforeMove; + } + + return prob; + +} +/** + * @details This function fills in tags in chainsEnds. Using bit shift operations it fills + * first 2 bits with distance from start of the chain and second 2 bits distance from end + * of the chain. + * + * @tparam IngredientsType The type of the system including all features. + * @param [in] ingredients A reference to the IngredientsType - mainly the system. + * + **/ + +template +void FeatureBendingPotential::setChainEnds(IngredientsType& ingredients){ + //First find out all the bonded monomers to initiate the size of chain ends. + //Throw error if solvent monomers are created before objects. + numNonSolvent=0; + bool isSolventFirst=true; + for(size_t i=0;i1 would act same. + chainsEnds.resize(numNonSolvent,2+(2<<2)); + + for(size_t i=0;i1){ + + int32_t numSameConnection=0; + int32_t sameAttNeigbIndex=0; + + for(int32_t neigb=0;neigbindex){ + int32_t tempoOldVal = chainsEnds[index]&(3<<2); + chainsEnds[index] = (0&3); + chainsEnds[index] += tempoOldVal; + + tempoOldVal = chainsEnds[index+1]&(3<<2); + chainsEnds[index+1] = (1&3); + chainsEnds[index+1] += tempoOldVal; + } + else{ + int32_t tempoOldVal = chainsEnds[index]&3; + chainsEnds[index] = (0&3)<<2; + chainsEnds[index] += tempoOldVal; + + tempoOldVal = chainsEnds[index-1]&3; + chainsEnds[index-1] = (1&3)<<2; + chainsEnds[index-1] += tempoOldVal; + } +} +/** + * @param type monomer attribute tag in range [1,10]. + * @param energy magnitude of bending potential energy. + * @throw std::runtime_error In case type exceed range [1,255]. + **/ +template +void FeatureBendingPotential::setBendingPotential(IngredientsType& ingredients, + int32_t type, + double energy) +{ + VectorInt3 bondVec1,bondVec2; + if(0::const_iterator it,it2; + //go over all the bondvector pair possible + for(it=ingredients.getBondset().begin();it!=ingredients.getBondset().end();it++) + for(it2=ingredients.getBondset().begin();it2!=ingredients.getBondset().end();it2++){ + bondVec1=it->second; + bondVec2=it2->second; + //find the value from bending potential form + float potentialVal=BendingPotentials::simpleHarmonic(bondVec1,bondVec2); + probabilityLookup[type][bondVectorToIndex(bondVec1)][bondVectorToIndex(bondVec2)]=exp(-potentialVal*energy); + } + + std::cout<<"set bending potential for types "; + std::cout<. + +--------------------------------------------------------------------------------*/ + +#ifndef FEATURE_BENDING_POTENTIAL_READ_WRITE_H +#define FEATURE_BENDING_POTENTIAL_READ_WRITE_H + +/** + * @file + * @date 2019/08/11 + * @author Ankush Checkervarty + * @brief Def. and impl. of class templates ReadBendingPotential and WriteBendingPotential +**/ + +#include +#include +#include + +/** + * @class ReadBendingPotential + * @brief Handles BFM-file read command !bending_potential + * @tparam IngredientsType Ingredients class storing all system information. +**/ +template < class IngredientsType> +class ReadBendingPotential: public ReadToDestination +{ +public: + ReadBendingPotential(IngredientsType& i):ReadToDestination(i){ +// i.modifyBondset().addBFMclassicBondset(); + } + virtual ~ReadBendingPotential(){} + virtual void execute(); +}; + +/** + * @class WriteBendingPotential + * @brief Handles BFM-file write command !bending_potential + * @tparam IngredientsType Ingredients class storing all system information. +**/ +template +class WriteBendingPotential:public AbstractWrite +{ +public: + + //constructor sets the headerOnly tag, such that the potential + //is written only once at the beginning of the output file. + WriteBendingPotential(const IngredientsType& i) + :AbstractWrite(i){this->setHeaderOnly(true);} + + virtual ~WriteBendingPotential(){} + + virtual void writeStream(std::ostream& strm); +}; + +/////////////MEMBER IMPLEMENTATIONS //////////////////////////////////////////// + +/** + * @brief Executes the reading routine to extract \b !bending_potential. + * + * @throw fail to read monomer types or bending potential strength. + **/ +template +void ReadBendingPotential::execute() +{ + IngredientsType& ingredients=this->getDestination(); + std::istream& file=this->getInputStream(); + + + int32_t type; + float bpStrength; + + //set stream to throw exception on fail + file.exceptions(file.exceptions() | std::ifstream::failbit); + + try + { + file>>type; + file>>bpStrength; + } + catch(std::ifstream::failure e) + { + std::stringstream errormessage; + errormessage<<"ReadBendingPotential::execute().\n"; + errormessage<<"Could not read interaction from file\n"; + errormessage<<"Previous error: "< +void WriteBendingPotential::writeStream(std::ostream& stream) +{ + int32_t nSpecies=10; //number must fit into 8 bit (uint8_t based lattice) + stream<<"## bending potential magnitude at type in kT (default 0.0kT)\n"; + + for(int32_t type=1;type<=nSpecies;type++) + { + if(this->getSource().getBendingPotential(type)!=0.0) + { + stream<<"!bending_potential "<getSource().getBendingPotential(type)<<"\n"; + } + + } + stream<<"\n\n"; + +} + + + + +#endif // FEATURE_NN_INTERACTION_READ_WRITE_H diff --git a/include/LeMonADE/updater/UpdaterLipidsCreator.h b/include/LeMonADE/updater/UpdaterLipidsCreator.h new file mode 100644 index 0000000..44aa8cc --- /dev/null +++ b/include/LeMonADE/updater/UpdaterLipidsCreator.h @@ -0,0 +1,780 @@ +#ifndef UpdaterLipidsCreator_H +#define UpdaterLipidsCreator_H + +#include +#include +#include + +/** + * @file + * + * @class UpdaterLipidsCreator + * + * @brief Updater to create bilayer in the middle of the box along z axis with solvent and/or nano particles + * and facially amphiphillic copolymers at random positions. + * + * @details This is a simple implementation of a system setup starting from an empty ingredients. + * This updater requires FeatureAttributes. There are choices of two types of head structure (i) Ring shaped + * (ii) linear polymer with 3 monomers as was used in Werner et. al., Soft Matter 8:11714,2012. + * There is a possibility to create two types with different tail lengths specified by n1_ and n2_.
+ * \b Attributes:
+ * \b Lipids_shorter_length: Head monomers Tag 1 and Tail monomers Tag 2.
+ * \b Lipids_longer_length: Head monomers Tag 3 and Tail monomers Tag 4.
+ * \b Solvent_monomers: Tag 5.
+ * \b Nano_particles: Tag 7.
+ * \b Facially_Amphiphlic_Copolymers: backbone monomers Tag 7 and side group monomers Tag 8.
+ * + * @tparam IngredientsType + * + * @param ingredients_ The system, holding an empty simulation box for system setup. + * + * @param nlipids_ number of lipids required in bilayer. + * @param Bilayer_Half_ if needed a half membranes(covers the 75% of xy plane). + * @param ringHead if needed a ring shaped head instead of normal three monomers. + * @param MeanLatticeOccupany_ Mean Lattice Occupany of simple cubic lattice(default 0.5). + * @param n1_ is length of tails of lipid type 1. + * @param n2_ is length of tails of lipid type 2. + * @param n_fraction_ fraction of lipid type 1. + * @param NanoType_ if Nano particles are to be included type 0 then tetrahedral nano particle + * would be created otherwise triangular shaped. + * @param NanoNum_ number of nano particles in the system. + * @param polylength_ length of facially amphiphillic copolymers to be created. + * @param numpoly_ number of facially amphiphillic copolymers to be created. + **/ + +using namespace std; + +template +class UpdaterLipidsCreator:public AbstractUpdater +{ +public: + UpdaterLipidsCreator(IngredientsType& ingredients_,int32_t nlipids_=300,double MeanLatticeOccupany_=0.5, int32_t n1_=5, int32_t n2_=0, double n_fraction_=0,bool ringHead_=false,bool Bilayer_Half_=false,int32_t NanoType_=1,int32_t NanoNum_=0,int32_t polylength_=5,int32_t numpoly_=0); + + //!General updater functions + bool execute(); + void initialize(){execute();}; + void cleanup(){}; + + + //!Function used by function lipidsBilayerCreator to find a vacant position where lipid has to be created. + VectorInt3 vacantLipidPosFinder(MoveAddMonomerSc& move,int32_t n1,int32_t upperleaflet_); + + //!Function used by function lipidsBilayerCreator to create a single lipid. + void singleLipidCreator(MoveAddMonomerSc& move, VectorInt3 pos,int32_t TotalMonLipids,int32_t upperleaflet); + + + //!Function for creating lipid bilayer in the middle of the box. + void lipidsBilayerCreator(MoveAddMonomerSc& move); + + //!Function used by the functions faciallyAmphCopolymers, nanoTriangle and nanoTetrahedron + //!to find shapes vacant for objects to be created. + VectorInt3 vacantShapeFinder(MoveAddMonomerSc& move, std::vector shape); + + //!Function for creating Facially Amphiphlic Copolymers. + void faciallyAmphCopolymers(MoveAddMonomerSc& move,int32_t polylength, int32_t numpoly); + + //!Function for creating Nano particles in the triangular shape. + void nanoTriangle(MoveAddMonomerSc& move,int32_t TotNanoNum_); + + //!Function for creating Nano particles in the tetrahedral shape. + void nanoTetrahedron(MoveAddMonomerSc& move, int32_t TotNanoN_); + + //!Function for creating ring shaped head on the lipid molecules. + void ringHeadCreator(MoveAddMonomerSc& move,VectorInt3 pos, int32_t upperleaflet_, int32_t type); + + //!Function for creating solvent when all other particles have been created. + void solventCreator(MoveAddMonomerSc& move, double solventdens_); + + //!Function to get upperleaflet lipids. + int32_t getUpperleafletLipids(){return leafletcounter;} + +protected: + //!Ingredients holding the simulation box + IngredientsType& ingredients; + + //!Number of lipids to be created + int32_t nlipids; + + //!Bool variable for whether half bilayer is needed or ring head is needed in the lipids. + bool Bilayer_Half,ringHead; + + //!Mean lattice Occupancy of the lattice. + double MeanLatticeOccupany; + + //!Mean lattice Occupancy of the lattice. + std::vector upperleafletBondvecs,lowerleafletBondvecs; + + RandomNumberGenerators randomNumbers; + + int32_t leafletcounter; + int32_t n1,n2,shorterlipid_length,longerlipid_length,shorterlipid_number,longerlipid_number; + int32_t NanoType,NanoNumber,polylength,numpoly; + double n_fraction; +}; + +/** +* @brief Constructor +**/ + +template +UpdaterLipidsCreator::UpdaterLipidsCreator(IngredientsType& ingredients_,int32_t nlipids_,double MeanLatticeOccupany_, int32_t n1_, int32_t n2_, double n_fraction_, bool ringHead_,bool Bilayer_Half_, int32_t NanoType_,int32_t NanoNum_,int32_t polylength_,int32_t numpoly_) + :ingredients(ingredients_), + nlipids(nlipids_), + MeanLatticeOccupany(MeanLatticeOccupany_), + Bilayer_Half(Bilayer_Half_), + ringHead(ringHead_), + n1(n1_), + n2(n2_), + n_fraction(n_fraction_), + NanoType(NanoType_), + NanoNumber(NanoNum_), + polylength(polylength_), + numpoly(numpoly_) +{ + + upperleafletBondvecs.push_back(ingredients.getBondset().getBondVector(21)); + upperleafletBondvecs.push_back(ingredients.getBondset().getBondVector(33)); + upperleafletBondvecs.push_back(ingredients.getBondset().getBondVector(27)); + + lowerleafletBondvecs.push_back(ingredients.getBondset().getBondVector(18)); + lowerleafletBondvecs.push_back(ingredients.getBondset().getBondVector(30)); + lowerleafletBondvecs.push_back(ingredients.getBondset().getBondVector(24)); + + randomNumbers.seedAll(); + + /** The choice of these bond vectors implies two conditions: + * + * (i) the height of lipids in z axis would be lipid_length+2. + * + * (ii) all the lipids are created 4 lattice units away from each other + * to avoid the excluded volume conflicts. This restricts the number of + * lipids that can be filled, but as it turns out maximum number of lipids + * which are filled this way already voilates the "almost flat condition", + * which is generally required by the bilayer related simulations, therefore, + * these bond vectors suffice the simulations of bilayer. + */ + + /********************************************************************************/ + + /**If both n1 and n2 are present decide what are the number and total length of shorter + * and longer lipids (assuming three head monomers). There are actually + * 4 head monomers in the case of ring head. But this case is separately handled + * function singleLipidCreator. + */ + if(n1&&n2){ + shorterlipid_length=(n1n2)?2*n1+3:2*n2+3; + shorterlipid_number=(shorterlipid_length==2*n1+3)?(nlipids*n_fraction):((1-n_fraction)*nlipids); + longerlipid_number=(longerlipid_length==2*n1+3)?(nlipids*n_fraction):((1-n_fraction)*nlipids); + } + + //if one of them is zero there is only one type of lipid, + //assign its important properties to shorterlipid_length + else{ + shorterlipid_length=(n1)?2*n1+3:2*n2+3; + longerlipid_length=shorterlipid_length; + shorterlipid_number=nlipids+1; + longerlipid_number=-1;} +} + +/** +* @brief Execution of the system creation +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template +bool UpdaterLipidsCreator::execute() +{ + + MoveAddMonomerSc move; + + //First create the Lipids + lipidsBilayerCreator(move); + + //Create facially Amphiphlic Copolymers + faciallyAmphCopolymers(move,polylength, numpoly); + + //Create Nano particles + if(NanoType) + nanoTriangle(move, NanoNumber); + else + nanoTetrahedron(move, NanoNumber); + + //Lastly create Solvent + solventCreator(move, MeanLatticeOccupany); + return true; +} + +/** +* @brief Creates a bilayer with midplane at the middle of z axis of the box. +* It distributes equal number of lipids in both leaflets of bilayer. +* Also, it creates n_fraction of n1 lipids. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ + +template +void UpdaterLipidsCreator::lipidsBilayerCreator(MoveAddMonomerSc& move) +{ + //Initialize the variables to calculate number of lipids in leaflet and number + //of types of lipids. + int32_t typeofLipid; + int32_t shorterlipidcounter=0,longerlipidcounter=0; + int32_t leaflet; + leafletcounter=0; + bool Oneleaffull=false,OneofTypeComplete=false; + VectorInt3 InitialPos; + + for(int32_t lipidnumber =0;lipidnumber +VectorInt3 UpdaterLipidsCreator::vacantLipidPosFinder(MoveAddMonomerSc& move,int32_t lipid_type,int32_t leaflet){ + + //Initializing a move to check for vacant position in lattice. + move.init(ingredients); + VectorInt3 position; + int32_t x,y,z; + + // \grow decides if lipid grows upwards or downwards (lower or upper leaflets respectively). + int32_t grow=(leaflet)?1:-1; + + //Here, we choose shorterlipid_length as if position is vacant for shorter lipids, + //it is also vacant for longer lipids. + + int32_t height_lipids=shorterlipid_length+2; + + //do-while loop performs till a vacant position is found + + do{ + + //In the case of bilayer half only 3/4 of x-axis is covered + //the factor 1/4 comes due to the choice of bond vectors. + // See the explanation in constructor. + + if(Bilayer_Half) + x=rand()%(3*(ingredients.getBoxX())/16); + else + x=rand()%((ingredients.getBoxY())/4); + + x=x*4; + + //In the case of ringHead, head monomers covers more lattice + //units on y-axis. + + if(ringHead){ + y=rand()%((ingredients.getBoxY())/6); + y=y*6;} + else{ + y=rand()%((ingredients.getBoxY())/4); + y=y*4;} + + //The point on the z axis where lipid starts growing from. + + z=(ingredients.getBoxZ()/2)+height_lipids*grow; + + position.setAllCoordinates(x,y,z); + move.setPosition(position); + }while(move.check(ingredients)==false); + + //Change z position if longer lipids have to be grown. + + if(lipid_type==shorterlipid_length) + return position; + + else{ + z=(ingredients.getBoxZ()/2)+(longerlipid_length+2)*grow; + position.setAllCoordinates(x,y,z); + return position; + } +} + +/** +* @brief Creates a single lipid on given any intial position, total monomers present in lipid and leaflet(upper or lower). +* +* @tparam IngredientsType Features used in the system. See Ingredients. +* +* @param move MoveAddMonomerSc type move used to apply or create a single lipid. +* @param Initpos intial position of lipid from where lipid starts growing. +* @param TotalMonLipids Total monomers in the lipid to be grown. +* @param leaflet upper or lower leaflet. +*/ + +template +void UpdaterLipidsCreator::singleLipidCreator(MoveAddMonomerSc& move,VectorInt3 initPos,int32_t totalMonLipids,int32_t leaflet){ + + std::vector bondVecSet; + + if(leaflet) + bondVecSet=upperleafletBondvecs; + else + bondVecSet=lowerleafletBondvecs; + + // Initialize tail length variable which length of tails which can calculated by + // assuming three head monomers. + + int32_t tail_length=(totalMonLipids-3)/2,lastMonomerAdded; + + //type_chr handles the tag to be given to the head and tail monomers + //of different kind of lipids. + + int32_t type_chr=(totalMonLipids==shorterlipid_length)?1:3; + + VectorInt3 position; + position=initPos; + VectorInt3 bondVector(0,0,0); + + for(int32_t monNum=0;monNum 0&& monNum!=(totalMonLipids-tail_length)) + ingredients.modifyMolecules().connect(lastMonomerAdded,lastMonomerAdded-1); + + // Choosing the different bondvector if it is the last head monomer. + + if(monNum==2){ + bondVector=bondVecSet[1]; + } + + //Connection to second tail and last head monomer. + + if(monNum==totalMonLipids-tail_length) + ingredients.modifyMolecules().connect(lastMonomerAdded,lastMonomerAdded-tail_length-1); + + //Change the position and bondvector at the end first tail. + + if(monNum==totalMonLipids-tail_length-1){ + position=ingredients.getMolecules()[lastMonomerAdded-tail_length]; + bondVector=bondVecSet[2]; + } + } + + else + throw std::runtime_error("Excluded volume conflict: could not add a single lipid"); + } + } + +/** +* @brief Creates a ring shape head in lipids, given the position, leaflet(upper or lower) and type. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +* +* @param move MoveAddMonomerSc type move used to apply or create a ring head. +* @param Initpos intial position of lipid; from where lipid starts growing. +* @param leaflet upper or lower leaflet. +* @param type tag assigned to the monomers. +*/ + +template +void UpdaterLipidsCreator::ringHeadCreator(MoveAddMonomerSc& move, VectorInt3 Initpos,int32_t leaflet,int32_t type){ + + move.init(ingredients); + int32_t Factor=(leaflet)?-1:1; + VectorInt3 v1(1,0,Factor*2); + VectorInt3 v2(-1,0,Factor*2); + VectorInt3 v3(0,2,Factor*2); + VectorInt3 v4(0,0,Factor*4); + move.setPosition(Initpos); + move.setTag(type); + move.apply(ingredients); + + move.setPosition(Initpos+v1); + move.setTag(type); + move.apply(ingredients); + int32_t lastMonomerAdded = ingredients.getMolecules().size()-1; + ingredients.modifyMolecules().connect(lastMonomerAdded,lastMonomerAdded-1); + + move.setPosition(Initpos+v2); + move.setTag(type); + move.apply(ingredients); + lastMonomerAdded = ingredients.getMolecules().size()-1; + ingredients.modifyMolecules().connect(lastMonomerAdded,lastMonomerAdded-2); + + move.setPosition(Initpos+v4); + move.setTag(type); + move.apply(ingredients); + lastMonomerAdded = ingredients.getMolecules().size()-1; + ingredients.modifyMolecules().connect(lastMonomerAdded,lastMonomerAdded-1); + + ingredients.modifyMolecules().connect(lastMonomerAdded,lastMonomerAdded-2); +} + +/** +* @brief Finds the vacant shape on avoiding excluded volume conflicts. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +* +* @param move MoveAddMonomerSc type move used to find a vacant shape on the lattice. +* @param shape vector of 3d coordinates of shape which has to be found on the lattice. +*/ + + +template +VectorInt3 UpdaterLipidsCreator::vacantShapeFinder(MoveAddMonomerSc& move, std::vector shape) +{ + //x,y,z are the initial coordinates of shape which would be returned later. + //Midplane and ShiftValue are found so that they do not have any conflicts with the lipids. + + int32_t x,y,z,counteroutside=0; + VectorInt3 position; + + bool occupied; + + //Do till an unoccupied position is found. + + do{ + occupied =false; + + //Random initial position + x=rand()%(ingredients.getBoxX()); + y=rand()%(ingredients.getBoxY()); + z=rand()%(ingredients.getBoxZ()); + + move.init(ingredients); + position.setAllCoordinates(x,y,z); + move.setPosition(position); + + //Find if the object can inserted at this place + for(int32_t monNum=0;monNum +void UpdaterLipidsCreator::faciallyAmphCopolymers(MoveAddMonomerSc& move,int32_t polylength, int32_t numpoly) +{ + //Create a vector of 3d coordinates to store shape of copolymer. + + std::vector shapevector; + VectorInt3 vec; + + for(int32_t polyMono=0;polyMono +void UpdaterLipidsCreator::nanoTriangle(MoveAddMonomerSc& move,int32_t totNanoNum) +{ + // Create a vector of 3d coordinates to store shape of nano particle. + + VectorInt3 v1(2,1,0); + VectorInt3 v2(1,0,2); + + std::vector shapevector; + VectorInt3 vec(0,0,0); + shapevector.push_back(vec); + shapevector.push_back(vec+v1); + shapevector.push_back(vec+v2); + + VectorInt3 Initpos; + + for(int32_t nanoNumber=0;nanoNumber +void UpdaterLipidsCreator::nanoTetrahedron(MoveAddMonomerSc& move,int32_t totNanoNum) +{ + VectorInt3 v1(1,0,2); + VectorInt3 v2(-1,0,2); + VectorInt3 v3(0,2,1); + + std::vector shapevector; + VectorInt3 vec(0,0,0); + shapevector.push_back(vec); + shapevector.push_back(vec+v1); + shapevector.push_back(vec+v2); + shapevector.push_back(vec+v3); + + VectorInt3 Initpos; + + for(int32_t nanoNumber=0;nanoNumber +void UpdaterLipidsCreator::solventCreator(MoveAddMonomerSc& move,double meanLatticeOccupany) +{ + + //Initialize the variables required to create solvent. + //Calculate the number of solvent particles needed to be created in order to maintain the given + // mean lattice occupancy. + + int32_t x,y,z; + int32_t total_monomers_added = ingredients.getMolecules().size(); + size_t numsol = (((ingredients.getBoxX())*(ingredients.getBoxY())*(ingredients.getBoxZ())*MeanLatticeOccupany)-(total_monomers_added*8))/8; + + int32_t initIndex=total_monomers_added-1; + + size_t i=0; + int32_t Midplane=ingredients.getBoxZ()/2; + + while(inumsol/2) + z=rand()%(Midplane-10) + (Midplane+10); + } + } + + else{ + if(inumsol/2) + z=rand()%(Midplane-10) + (Midplane+10);} + + VectorInt3 position(x,y,z); + move.setTag(5); + move.setPosition(position); + if(move.check(ingredients)==true){ + move.apply(ingredients); + i++; + } + int32_t lastMonomerAdded=ingredients.getMolecules().size()-1; + ingredients.setCompressedOutputIndices(initIndex+1,lastMonomerAdded); + } +}; + +#endif diff --git a/include/LeMonADE/updater/UpdaterLipidsRandomInitializer.h b/include/LeMonADE/updater/UpdaterLipidsRandomInitializer.h new file mode 100644 index 0000000..7452d86 --- /dev/null +++ b/include/LeMonADE/updater/UpdaterLipidsRandomInitializer.h @@ -0,0 +1,251 @@ +#include +#include +/** + * @file + * + * @class UpdaterLipidsRandomInitializer + * + * @brief Updater to create lipids at random position and random orientations with choice of two types of lipid + * solvent monomers. + * + * @details This is upgrade to UpdaterLipidsCreator where instead of creating the lipids in form of bilayer, + * lipids are generated at random postions and orientations.There is a possibility to create two types + * with different tail lengths specified by n1_ and n2_.
+ * \b Attributes:
+ * \b Lipids_shorter_length: Head monomers Tag 1 and Tail monomers Tag 2.
+ * \b Lipids_longer_length: Head monomers Tag 3 and Tail monomers Tag 4.
+ * \b Solvent_monomers: Tag 5.
+ * + * @tparam IngredientsType + * + * @param ingredients_ The system, holding an empty simulation box for system setup. + * + * @param nlipids_ number of lipids required in bilayer. + * @param MeanLatticeOccupany_ Mean Lattice Occupany of simple cubic lattice(default 0.5). + * @param n1_ is length of tails of lipid type 1. + * @param n2_ is length of tails of lipid type 2. + * @param n_fraction_ fraction of lipid type 1. + * + **/ + +template +class UpdaterLipidsRandomInitializer:public UpdaterLipidsCreator +{ + typedef UpdaterLipidsCreator BaseClass; + typedef UpdaterSimpleSimulator SimpSimu; + +public: + UpdaterLipidsRandomInitializer(); + //!Allowing access to constructors of UpdaterLipidsCreator. + using BaseClass::UpdaterLipidsCreator; + //!Allowing access to other function of this class. + using BaseClass::solventCreator; + using BaseClass::ingredients; + //general updater functions + void initialize(); + bool execute(); + void cleanup(){}; + + //!Function for initiating lipid molecules at random position in the box. + void randomInitializer(MoveAddMonomerSc& move); + + //!Function for creating a single lipid in random orientation. + void singleLipidCreator(MoveAddMonomerSc& move,int32_t totalMonLipids,SimpSimu& simulater); + + //!Function draw a random bond vector which are allowed by BFM creatorian. + VectorInt3 drawRandomDirection() const; + +private: + + //!Redecleration of the variables from the base class to be used in this class. + int32_t longerlipid_length,shorterlipid_length,nlipids,n_fraction,shorterlipid_number,longerlipid_number; + float MeanLatticeOccupany; + bool rI; + +}; + +/** +* @brief Default constructor +**/ +template +UpdaterLipidsRandomInitializer::UpdaterLipidsRandomInitializer(){ + std::cout<<"hello"; + exit(0); +}; + +template +void UpdaterLipidsRandomInitializer::initialize() +{ + //Use the base class variables in this class. + this->nlipids = BaseClass::nlipids; + this->shorterlipid_length = BaseClass::shorterlipid_length; + this->longerlipid_length = BaseClass::longerlipid_length; + this->shorterlipid_number = BaseClass::shorterlipid_number; + this->longerlipid_number = BaseClass::longerlipid_number; + this->n_fraction = BaseClass::n_fraction; + this->MeanLatticeOccupany = BaseClass::MeanLatticeOccupany; + + execute(); +}; + +template +bool UpdaterLipidsRandomInitializer::execute() +{ + //Move to add monomers + MoveAddMonomerSc move; + + //Initialize the lipids at random position and random orientation. + // and then create solvent. + randomInitializer(move); + solventCreator(move, rI=true); +}; + +/** +* @brief Creates lipids randomly everywher in the box. +* It creates n_fraction of n1 lipids. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ + +template +void UpdaterLipidsRandomInitializer::randomInitializer(MoveAddMonomerSc& move) +{ + //Initialize the variables to calculate number of lipids number + //of lipids type of lipid. Also Instantiate an object of simple simulater + //as the it might be required to move the system. + + int32_t typeofLipid; + SimpSimu sim(ingredients,1000); + int32_t shorterlipidcounter=0,longerlipidcounter=0; + bool OneofTypeComplete=false; + VectorInt3 InitialPos; + + for(int32_t lipidnumber =0;lipidnumber +void UpdaterLipidsRandomInitializer::singleLipidCreator(MoveAddMonomerSc& move,int32_t totalMonLipids,SimpSimu& sim){ + + //Initialize the variables to calculate rejected trials to create monomers + //and other variables to be used later. + + int32_t rejectedMonoTrials,x,y,z,type,lastMonomerAdded,monNum; + VectorInt3 position; + + int32_t tail_length=(totalMonLipids-3)/2; + int32_t type_chr=(totalMonLipids==shorterlipid_length)?1:3; + + VectorInt3 bondVector(0,0,0); + + //Random initial position. + x=rand()%(ingredients.getBoxX()); + y=rand()%(ingredients.getBoxY()); + z=rand()%(ingredients.getBoxZ()); + position.setAllCoordinates(x,y,z); + rejectedMonoTrials=0; + monNum=0; + while(monNum0 && monNum!=(totalMonLipids-tail_length)) + ingredients.modifyMolecules().connect(lastMonomerAdded,lastMonomerAdded-1); + + if(monNum==totalMonLipids-tail_length) + ingredients.modifyMolecules().connect(lastMonomerAdded,lastMonomerAdded-tail_length-1); + + + if(monNum==totalMonLipids-tail_length-1) + position=ingredients.getMolecules()[lastMonomerAdded-tail_length]; + + monNum++; + } + + else{ + rejectedMonoTrials++; + + //If on the average all the bond vectors are used, + //move the system by simply simulating for 1000 MCS. + + if(rejectedMonoTrials==108){ + sim.initialize(); + sim.execute(); + lastMonomerAdded = ingredients.getMolecules().size()-1; + position=ingredients.getMolecules()[lastMonomerAdded]; + + if(monNum==totalMonLipids-tail_length-1) + position=ingredients.getMolecules()[lastMonomerAdded-tail_length]; + + rejectedMonoTrials=0; + } + } + + bondVector=drawRandomDirection(); + } +}; + +/** +* @brief Returns a random vector from the set of bond vectorss. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template +VectorInt3 UpdaterLipidsRandomInitializer::drawRandomDirection() const +{ + VectorInt3 bondVector; + size_t randVectorId=rand()%(ingredients.getBondset().size()); + randVectorId=randVectorId%46; + randVectorId += (randVectorId<17) ? 17 : 0; + bondVector=ingredients.getBondset().getBondVector(randVectorId); + return bondVector; +}; + diff --git a/tests/analyzer/TestAnalyzerPermeabilityCalculator.cpp b/tests/analyzer/TestAnalyzerPermeabilityCalculator.cpp new file mode 100644 index 0000000..f13a4ef --- /dev/null +++ b/tests/analyzer/TestAnalyzerPermeabilityCalculator.cpp @@ -0,0 +1,373 @@ +/*-------------------------------------------------------------------------------- + * 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 the class AnalyzerPermeabilityCalculator + * + * @author Ankush Checkervarty + * @date 09.08.2019 + * */ +/*****************************************************************************/ + +#include "gtest/gtest.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + + +using namespace std; +// class with body that is needed for every test, only used by TEST_F() +/*****************************************************************************/ +/** + * @class AnalyzerPermeabilityCalculatorTest + * @brief prepare system and check functions for AnalyzerPermeabilityCalculatorTest + * */ +/*****************************************************************************/ +class AnalyzerPermeabilityCalculatorTest: public ::testing::Test{ +protected: + void setInitConfig() + { + ingredients.setBoxX(64); + ingredients.setBoxY(64); + ingredients.setBoxZ(64); + ingredients.setPeriodicX(true); + ingredients.setPeriodicY(true); + ingredients.setPeriodicZ(true); + ingredients.synchronize(ingredients); + } + + /** This function setups the system with two lipid monomers, two + * disconnected object monomers and a solvent monomer. + * First two particles are objects and 3rd one is solvent. + */ + void setConfig(int32_t midplane, int32_t width) + { + + ingredients.modifyMolecules().resize(5); + + ingredients.modifyMolecules()[0].setAllCoordinates(2,3,midplane+width); + ingredients.modifyMolecules()[0].setAttributeTag(1); + + ingredients.modifyMolecules()[1].setAllCoordinates(2,3,midplane-width); + ingredients.modifyMolecules()[1].setAttributeTag(2); + + ingredients.modifyMolecules()[2].setAllCoordinates(2,5,midplane+width); + ingredients.modifyMolecules()[2].setAttributeTag(7); + + ingredients.modifyMolecules()[3].setAllCoordinates(2,5,midplane-width); + ingredients.modifyMolecules()[3].setAttributeTag(8); + + ingredients.modifyMolecules()[4].setAllCoordinates(2,5,midplane+factor2sigma*width+1); + ingredients.modifyMolecules()[4].setAttributeTag(5); + ingredients.synchronize(ingredients); + } + + + //! set the system into one of two model configs, designed to give easily calculated Rg^2 + std::string passComments(std::ifstream& file) + { + std::string linecontent; + + for(int32_t i=0;;i++){ + std::getline(file,linecontent); + if(linecontent[0]!='#'&&linecontent.size()>0) + break;} + + return linecontent; + } + + //define system + typedef LOKI_TYPELIST_2(FeatureMoleculesIO,FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients < Config> MyIngredients; + MyIngredients ingredients; + int32_t factor2sigma; + + //! nested derived class from AnalyzerPermeabilityCalculator to test protected functions + class PmAnalyzerDerived:public AnalyzerPermeabilityCalculator + { + public: + using AnalyzerPermeabilityCalculator::AnalyzerPermeabilityCalculator; + using AnalyzerPermeabilityCalculator::execute; + using AnalyzerPermeabilityCalculator::particlesInsideBoundaries; + using AnalyzerPermeabilityCalculator::particlesInsideBoundariesOld; + using AnalyzerPermeabilityCalculator::counterSolvent; + using AnalyzerPermeabilityCalculator::counterObjects; + using AnalyzerPermeabilityCalculator::midplane; + using AnalyzerPermeabilityCalculator::sigma; + + }; + /* 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(AnalyzerPermeabilityCalculatorTest, CheckConstructorandPermeabilityupdate) +{ + setInitConfig(); + + int32_t midplane=32; + int32_t width=4; + factor2sigma=3; + + setConfig(midplane,width); + /**setConfig() produces + * Particle 0:object1 inside the boundaries from positive side. + * Particle 1:object2 inside the boundaries from negative side. + * Particle 2:solvent outside the boundaries. + */ + std::vector > objects; + hasOneOfTheseTwoTypes<7,8> hastype; + fill_connected_groups(ingredients.getMolecules(),objects,MonomerGroup((ingredients.getMolecules())),hastype); + + + //making a bin size same as box, so there is only one bin. + PmAnalyzerDerived pmAnalyzer(ingredients,objects,64); + + //checking the constructor + EXPECT_EQ(midplane,pmAnalyzer.midplane[0][0]); + EXPECT_EQ(width,pmAnalyzer.sigma); + + EXPECT_EQ(3,pmAnalyzer.particlesInsideBoundariesOld.size()); + EXPECT_EQ(1,pmAnalyzer.particlesInsideBoundariesOld[0]); + EXPECT_EQ(-1,pmAnalyzer.particlesInsideBoundariesOld[1]); + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundariesOld[2]); + + EXPECT_EQ(3,pmAnalyzer.particlesInsideBoundaries.size()); + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundaries[0]); + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundaries[1]); + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundaries[2]); + + EXPECT_EQ(0,pmAnalyzer.counterObjects); + EXPECT_EQ(0,pmAnalyzer.counterSolvent); + + //moving objects out of boundaries to opposite sides. + //and moving solvent positive side from boundary. + ingredients.modifyMolecules()[2].setZ(midplane-factor2sigma*width-1); + ingredients.modifyMolecules()[3].setZ(midplane+factor2sigma*width+1); + ingredients.modifyMolecules()[4].setZ(midplane+width); + + pmAnalyzer.execute(); + + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundariesOld[0]); + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundariesOld[1]); + EXPECT_EQ(1,pmAnalyzer.particlesInsideBoundariesOld[2]); + + //Two object translocated + EXPECT_EQ(2,pmAnalyzer.counterObjects); + EXPECT_EQ(0,pmAnalyzer.counterSolvent); + + //Objects move back, solvent is opposite side of midplane + ingredients.modifyMolecules()[2].setZ(midplane-(factor2sigma-1)*width); + ingredients.modifyMolecules()[3].setZ(midplane+(factor2sigma-1)*width); + ingredients.modifyMolecules()[4].setZ(midplane-(factor2sigma-1)*width); + + pmAnalyzer.execute(); + + EXPECT_EQ(-1,pmAnalyzer.particlesInsideBoundariesOld[0]); + EXPECT_EQ(1,pmAnalyzer.particlesInsideBoundariesOld[1]); + + //Analyzer remember where solvent entered into boundaries + EXPECT_EQ(1,pmAnalyzer.particlesInsideBoundariesOld[2]); + + EXPECT_EQ(2,pmAnalyzer.counterObjects); + EXPECT_EQ(0,pmAnalyzer.counterSolvent); + + //Objects move to other side of midplane, solvent + //moved out opposite side from where entered in. + ingredients.modifyMolecules()[2].setZ(midplane+(factor2sigma-1)*width); + ingredients.modifyMolecules()[3].setZ(midplane-(factor2sigma-1)*width); + ingredients.modifyMolecules()[4].setZ(midplane-factor2sigma*width-1); + + pmAnalyzer.execute(); + //analyzer remembers where it entered through + EXPECT_EQ(-1,pmAnalyzer.particlesInsideBoundariesOld[0]); + EXPECT_EQ(1,pmAnalyzer.particlesInsideBoundariesOld[1]); + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundariesOld[2]); + + EXPECT_EQ(2,pmAnalyzer.counterObjects); + EXPECT_EQ(1,pmAnalyzer.counterSolvent); + + //Objects move out of boundaries however one of them from + //same side that it came in. Thus only one event is counted. + ingredients.modifyMolecules()[2].setZ(midplane-factor2sigma*width-1); + ingredients.modifyMolecules()[3].setZ(midplane-factor2sigma*width-3); + ingredients.modifyMolecules()[4].setZ(midplane-factor2sigma*width-5); + + pmAnalyzer.execute(); + + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundariesOld[0]); + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundariesOld[1]); + EXPECT_EQ(0,pmAnalyzer.particlesInsideBoundariesOld[2]); + + EXPECT_EQ(3,pmAnalyzer.counterObjects); + EXPECT_EQ(1,pmAnalyzer.counterSolvent); + +} + +TEST_F(AnalyzerPermeabilityCalculatorTest, CheckFileDumpandMidplaneUpdate) +{ + setInitConfig(); + int32_t midplane=32; + int32_t width=4; + + setConfig(midplane,width); + + std::vector > objects; + hasOneOfTheseTwoTypes<7,8> hastype; + fill_connected_groups(ingredients.getMolecules(),objects,MonomerGroup((ingredients.getMolecules())),hastype); + + //making a bin size same as box, so there is only one bin. + std::string outputFilename="TestDump.dat"; + PmAnalyzerDerived pmAnalyzer1(ingredients,objects,64,3,outputFilename,1,150,1); + + pmAnalyzer1.counterObjects=1; + pmAnalyzer1.counterSolvent=1; + + pmAnalyzer1.execute(); + + std::ifstream file(outputFilename.c_str()); + + EXPECT_TRUE(file.is_open()); + + std::string linecontent = passComments(file); + EXPECT_EQ(linecontent[0],'0'); + EXPECT_EQ(linecontent[2],'1'); + EXPECT_EQ(linecontent[4],'1'); + + pmAnalyzer1.counterObjects=2; + pmAnalyzer1.counterSolvent=3; + + + ingredients.modifyMolecules().setAge(1); + pmAnalyzer1.execute(); + + std::getline(file,linecontent); + EXPECT_EQ(linecontent[0],'1'); + EXPECT_EQ(linecontent[2],'2'); + EXPECT_EQ(linecontent[4],'3'); + + remove("TestDump.dat"); + //midplane update Check + PmAnalyzerDerived pmAnalyzer2(ingredients,objects,64,3,outputFilename,100,1,1); + + EXPECT_EQ(midplane,pmAnalyzer2.midplane[0][0]); + EXPECT_EQ(width,pmAnalyzer2.sigma); + + midplane=35; + width=2; + setConfig(midplane,width); + pmAnalyzer2.execute(); + + EXPECT_EQ(midplane,pmAnalyzer2.midplane[0][0]); + EXPECT_EQ(width,pmAnalyzer2.sigma); + + midplane=28; + width=6; + setConfig(midplane,width); + pmAnalyzer2.execute(); + + EXPECT_EQ(midplane,pmAnalyzer2.midplane[0][0]); + EXPECT_EQ(width,pmAnalyzer2.sigma); + +} + +TEST_F(AnalyzerPermeabilityCalculatorTest, CheckPoreUpdate) +{ + setInitConfig(); + int32_t midplane=32; + int32_t width=4; + setConfig(midplane,width); + + std::vector > objects; + hasOneOfTheseTwoTypes<7,8> hastype; + fill_connected_groups(ingredients.getMolecules(),objects,MonomerGroup((ingredients.getMolecules())),hastype); + + PmAnalyzerDerived pmAnalyzer3(ingredients,objects,32); + + EXPECT_EQ(midplane, pmAnalyzer3.midplane[0][0]); + EXPECT_EQ(midplane, pmAnalyzer3.midplane[0][1]); + EXPECT_EQ(midplane, pmAnalyzer3.midplane[1][0]); + EXPECT_EQ(midplane, pmAnalyzer3.midplane[1][1]); + + midplane=35; + width=4; + setConfig(midplane,width); + + std::string outputFilename="TestDump.dat"; + PmAnalyzerDerived pmAnalyzer4(ingredients,objects,16,3,outputFilename,100,1,1); + + EXPECT_EQ(midplane, pmAnalyzer4.midplane[0][0]); + EXPECT_EQ(midplane, pmAnalyzer4.midplane[0][1]); + EXPECT_EQ(midplane, pmAnalyzer4.midplane[1][0]); + EXPECT_EQ(midplane, pmAnalyzer4.midplane[1][1]); + + EXPECT_NE(midplane, pmAnalyzer4.midplane[2][0]); + EXPECT_NE(midplane, pmAnalyzer4.midplane[0][2]); + EXPECT_NE(midplane, pmAnalyzer4.midplane[2][2]); + + EXPECT_EQ(midplane, pmAnalyzer4.midplane[3][0]); + EXPECT_EQ(midplane, pmAnalyzer4.midplane[3][1]); + + ingredients.modifyMolecules().addMonomer(2*16+1,0,midplane+width); + ingredients.modifyMolecules()[5].setAttributeTag(1); + ingredients.modifyMolecules().addMonomer(2*16+1,0,midplane-width); + ingredients.modifyMolecules()[6].setAttributeTag(2); + + pmAnalyzer4.execute(); + EXPECT_EQ(midplane, pmAnalyzer4.midplane[2][0]); + EXPECT_EQ(midplane, pmAnalyzer4.midplane[2][1]); + EXPECT_EQ(midplane, pmAnalyzer4.midplane[3][0]); + EXPECT_EQ(midplane, pmAnalyzer4.midplane[3][1]); + remove("TestDump.dat"); + +} diff --git a/tests/analyzer/TestAnalyzerPoreFinder.cpp b/tests/analyzer/TestAnalyzerPoreFinder.cpp new file mode 100644 index 0000000..dbe8f1b --- /dev/null +++ b/tests/analyzer/TestAnalyzerPoreFinder.cpp @@ -0,0 +1,459 @@ +/*-------------------------------------------------------------------------------- + * 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 the class AnalyzerPoreFinder + * + * @author Ankush Checkervarty + * @date 09.08.2019 + * */ +/*****************************************************************************/ + +#include "gtest/gtest.h" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + + +using namespace std; +// class with body that is needed for every test, only used by TEST_F() +/*****************************************************************************/ +/** + * @class AnalyzerPoreFinderTest + * @brief prepare system and check functions for AnalyzerPoreFinderTest + * */ +/*****************************************************************************/ +class AnalyzerPoreFinderTest: public ::testing::Test{ +protected: + + void setInitConfig() + { + ingredients.setBoxX(32); + ingredients.setBoxY(32); + ingredients.setBoxZ(32); + ingredients.setPeriodicX(true); + ingredients.setPeriodicY(true); + ingredients.setPeriodicZ(true); + ingredients.synchronize(ingredients); + } + + int32_t reduceDistanceInPeriodicSpace(int32_t distance, int32_t period) + { + if(distance > period/2){ + while(-(distance-period) period/2){ + while((distance+period)(radius)*(radius)&& + shift_xp1*shift_xp1+shift_y*shift_y>(radius)*(radius)&& + shift_x*shift_x+shift_yp1*shift_yp1>(radius)*(radius)&& + shift_xp1*shift_xp1+shift_yp1*shift_yp1>(radius)*(radius) + ){ + ingredients.modifyMolecules().addMonomer(x,y,boxZ/2); + ingredients.modifyMolecules()[size].setAttributeTag(2); + size++; + } + } + } + /** This function sets up two pores with given radii and pore centers.*/ + void setTwoPores(int32_t radius1,int32_t radius2, VectorInt2 center1, VectorInt2 center2) + { + int32_t boxX=ingredients.getBoxX(); + int32_t boxY=ingredients.getBoxY(); + int32_t boxZ=ingredients.getBoxZ(); + int32_t boxXm_1=ingredients.getBoxX()-1; + int32_t boxYm_1=ingredients.getBoxY()-1; + int32_t size=0; + int32_t poreCenterX1=center1.getX(); + int32_t poreCenterY1=center1.getY(); + int32_t poreCenterX2=center2.getX(); + int32_t poreCenterY2=center2.getY(); + + for(int32_t x=0;x(radius1)*(radius1)&& + shift_x1p1*shift_x1p1+shift_y1*shift_y1>(radius1)*(radius1)&& + shift_x1*shift_x1+shift_y1p1*shift_y1p1>(radius1)*(radius1)&& + shift_x1p1*shift_x1p1+shift_y1p1*shift_y1p1>(radius1)*(radius1)&& + + shift_x2*shift_x2+shift_y2*shift_y2>(radius2)*(radius2)&& + shift_x2p1*shift_x2p1+shift_y2*shift_y2>(radius2)*(radius2)&& + shift_x2*shift_x2+shift_y2p1*shift_y2p1>(radius2)*(radius2)&& + shift_x2p1*shift_x2p1+shift_y2p1*shift_y2p1>(radius2)*(radius2) + ){ + ingredients.modifyMolecules().addMonomer(x,y,boxZ/2); + ingredients.modifyMolecules()[size].setAttributeTag(2); + size++; + } + } + } + + + //! set the system into one of two model configs, designed to give easily calculated Rg^2 + std::string passComments(std::ifstream& file) + { + std::string linecontent; + + for(int32_t i=0;;i++){ + std::getline(file,linecontent); + if(linecontent[0]!='#'&&linecontent.size()>0) + break;} + + return linecontent; + } + + //define system + typedef LOKI_TYPELIST_3(FeatureMoleculesIO,FeatureLatticePowerOfTwo,FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients MyIngredients; + MyIngredients ingredients; + + //!nested derived class from AnalyzerPoreFinder to test protected functions + class PmAnalyzerDerived: public AnalyzerPoreFinder + { + public: + using AnalyzerPoreFinder::AnalyzerPoreFinder; + using AnalyzerPoreFinder::execute; + using AnalyzerPoreFinder::cleanup; + using AnalyzerPoreFinder::coordinatesOfPore; + using AnalyzerPoreFinder::centriod; + }; + + /* 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(AnalyzerPoreFinderTest, ClusterAnalysisOnePore) +{ + RandomNumberGenerators randomNumbers; + randomNumbers.seedAll(); + + setInitConfig(); + /**create a pore in the middle the box of radius + */ + int32_t radius=3; + int32_t poreCenterX=16; + int32_t poreCenterY=16; + VectorInt2 center(poreCenterX,poreCenterY); + + setPore(radius,center); + /**setPore(3,16,16) produces + * monomers with tag 2 all over the middle of box in z direction + * in xy plane leaving a radius of 3(pore). + */ + + std::vector > objects; + hasOneOfTheseTwoTypes<7,8> hastype; + fill_connected_groups(ingredients.getMolecules(),objects,MonomerGroup((ingredients.getMolecules())),hastype); + + PmAnalyzerDerived analPore(ingredients,objects); + + analPore.execute(); + + float pi=3.142; + int32_t area=pi*radius*radius; + + //Since this is comparing discrete lattice points and real space, + //there is some uncertainty of 3-4 lattice points + int32_t diff= area-analPore.coordinatesOfPore.size(); + EXPECT_LE(abs(diff),4); + + //Only Centriod should be there + EXPECT_EQ(1,analPore.centriod.size()); + + //This is again due to difference between discrete and real space + EXPECT_LE(abs(analPore.centriod[0].getX()-poreCenterX),1); + + EXPECT_LE(abs(analPore.centriod[0].getY()-poreCenterY),1); + + //Now let's test it for random radius and random center position + //radius can be maximum half of box size + radius=rand()%5+3; + + poreCenterX=rand()%32; + poreCenterY=rand()%32; + VectorInt2 center2(poreCenterX,poreCenterY); + + //set Molecules size to zero again. + ingredients.modifyMolecules().resize(0); + + setPore(radius,center2); + analPore.execute(); + + area=pi*radius*radius; + diff= area-analPore.coordinatesOfPore.size(); + EXPECT_LE(abs(diff),4); + + //only one Centriod should be there + EXPECT_EQ(1,analPore.centriod.size()); + + //this is again due to difference between discrete and real space + EXPECT_LE(abs(analPore.centriod[0].getX()-poreCenterX),1); + + EXPECT_LE(abs(analPore.centriod[0].getY()-poreCenterY),1); + + radius=3; + poreCenterX=16; + poreCenterY=16; + + + analPore.execute(); + //check the cleanup and files. + analPore.cleanup(); + + std::ifstream file("PoreCoordinates.dat"); + + EXPECT_TRUE(file.is_open()); + + std::string linecontent = passComments(file); + EXPECT_EQ(linecontent[0],'0'); + EXPECT_EQ(linecontent[2],'0'); + EXPECT_EQ(linecontent[4],'0'); + + for(int32_t x=0;x<16;x++) + for(int32_t y=0;y<32;y++) + std::getline(file,linecontent); + + for(int32_t y=0;y<16;y++) + std::getline(file,linecontent); + + EXPECT_EQ(linecontent.substr(0,2),"16"); + EXPECT_EQ(linecontent.substr(3,2),"16"); + EXPECT_EQ(linecontent.substr(6,1),"1"); + + std::getline(file,linecontent); + + EXPECT_EQ(linecontent.substr(0,2),"16"); + EXPECT_EQ(linecontent.substr(3,2),"17"); + EXPECT_EQ(linecontent.substr(6,1),"1"); + + file.close(); + remove("PoreCoordinates.dat"); + + } +TEST_F(AnalyzerPoreFinderTest, ClusterAnalysisTwoPoreandCopolymers) +{ + setInitConfig(); + + int32_t radius1=3; + int32_t poreCenterX1=16; + int32_t poreCenterY1=16; + VectorInt2 center1(poreCenterX1,poreCenterY1); + + int32_t radius2=4; + int32_t poreCenterX2=26; + int32_t poreCenterY2=16; + VectorInt2 center2(poreCenterX2,poreCenterY2); + + setTwoPores(radius1,radius2,center1,center2); + + //let's create some copolymers + ingredients.modifyMolecules().addMonomer(poreCenterX1,poreCenterY1-9,16+5); + int32_t last_Index=ingredients.getMolecules().size()-1; + ingredients.modifyMolecules()[last_Index].setAttributeTag(8); + + ingredients.modifyMolecules().addMonomer(poreCenterX1,poreCenterY1-3,16+5); + last_Index=ingredients.getMolecules().size()-1; + ingredients.modifyMolecules()[last_Index].setAttributeTag(8); + + ingredients.modifyMolecules().addMonomer(poreCenterX1,poreCenterY1+3,16+5); + last_Index=ingredients.getMolecules().size()-1; + ingredients.modifyMolecules()[last_Index].setAttributeTag(7); + + + std::vector > objects; + hasOneOfTheseTwoTypes<7,8> hastype; + fill_connected_groups(ingredients.getMolecules(),objects,MonomerGroup((ingredients.getMolecules())),hastype); + + PmAnalyzerDerived analPore(ingredients,objects); + + analPore.execute(); + + float pi=3.142; + int32_t area=pi*radius1*radius1 + pi*radius2*radius2; + int32_t diff= area-analPore.coordinatesOfPore.size(); + + EXPECT_LE(abs(diff),4); + + //Two Centriods should be there. + EXPECT_EQ(2,analPore.centriod.size()); + + //According to the search pore with less value of x and y. + EXPECT_LE(abs(analPore.centriod[0].getX()-poreCenterX1),1); + + EXPECT_LE(abs(analPore.centriod[0].getY()-poreCenterY1),1); + + EXPECT_LE(abs(analPore.centriod[1].getX()-poreCenterX2),1); + + EXPECT_LE(abs(analPore.centriod[1].getY()-poreCenterY2),1); + + //check the cleanup and files. + analPore.cleanup(); + + std::ifstream file("PoreCoordinates.dat"); + + EXPECT_TRUE(file.is_open()); + + std::string linecontent = passComments(file); + + //check for the values at the center. + for(int32_t x=0;x<16;x++) + for(int32_t y=0;y<32;y++) + std::getline(file,linecontent); + + for(int32_t y=0;y<16;y++) + std::getline(file,linecontent); + + EXPECT_EQ(linecontent.substr(0,2),"16"); + EXPECT_EQ(linecontent.substr(3,2),"16"); + EXPECT_EQ(linecontent.substr(6,1),"1"); + + //If there are two pores, the analyzer + //treats them like two pore events. + //Thus, value at 4 radius points should 0.5! + for(int32_t y=16;y<16+4;y++) + std::getline(file,linecontent); + + EXPECT_EQ(linecontent.substr(0,2),"16"); + EXPECT_EQ(linecontent.substr(3,2),"20"); + EXPECT_EQ(linecontent.substr(6,3),"0.5"); + + file.close(); + + //check for copolymers. + std::ifstream file1("PolymerCoordinates.dat"); + + EXPECT_TRUE(file1.is_open()); + + linecontent = passComments(file1); + + EXPECT_EQ(linecontent[0],'0'); + EXPECT_EQ(linecontent[2],'0'); + EXPECT_EQ(linecontent[4],'0'); + + //For the first copolymer. + for(int32_t x=0;x<16;x++) + for(int32_t y=0;y<32;y++) + std::getline(file1,linecontent); + + for(int32_t y=0;y<16-9;y++) + std::getline(file1,linecontent); + + EXPECT_EQ(linecontent.substr(0,2),"16"); + EXPECT_EQ(linecontent.substr(3,1),"7"); + EXPECT_EQ(linecontent.substr(5,3),"0.5"); + + //second copolymer + for(int32_t y=16-9;y<16-3;y++) + std::getline(file1,linecontent); + + EXPECT_EQ(linecontent.substr(0,2),"16"); + EXPECT_EQ(linecontent.substr(3,2),"13"); + EXPECT_EQ(linecontent.substr(6,3),"0.5"); + + //third copolymer + for(int32_t y=16-3;y<16+3;y++) + std::getline(file1,linecontent); + + EXPECT_EQ(linecontent.substr(0,2),"16"); + EXPECT_EQ(linecontent.substr(3,2),"19"); + EXPECT_EQ(linecontent.substr(6,3),"0.5"); + + remove("PoreCoordinates.dat"); + remove("PolymerCoordinates.dat"); + +} diff --git a/tests/feature/TestFeatureBendingPotential.cpp b/tests/feature/TestFeatureBendingPotential.cpp new file mode 100644 index 0000000..b2e1563 --- /dev/null +++ b/tests/feature/TestFeatureBendingPotential.cpp @@ -0,0 +1,408 @@ +#include + +#include "gtest/gtest.h" + +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +class BendingPotentialTest: public ::testing::Test{ + +public: + /** This function setups a linear chain which has only one + * bondvector=(2,0,0) between monomers. + */ + template < class IngredientsType> + void setSimpleLinearChain(IngredientsType& ingredients, int32_t length) + { + ingredients.modifyMolecules().resize(length); + + for(int32_t i=0;i + void setMixedLinearChain(IngredientsType& ingredients, int32_t length) + { + ingredients.modifyMolecules().resize(length); + + for(int32_t i=0;iun-/comment here:*/ + //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(BendingPotentialTest, CheckLocalmoveProbability) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureBendingPotential) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + + + + //prepare myIngredients. + myIngredients.setBoxX(64); + myIngredients.setBoxY(64); + myIngredients.setBoxZ(64); + myIngredients.setPeriodicX(1); + myIngredients.setPeriodicY(1); + myIngredients.setPeriodicZ(1); + myIngredients.modifyBondset().addBFMclassicBondset(); + + //setup a simple Linear chain. + setSimpleLinearChain(myIngredients,5); + + //set lipid chainsEnds + myIngredients.setChainEnds(myIngredients); + + //set interaction for type 1 + myIngredients.setBendingPotential(myIngredients,1,0.5); + double bpStrength=0.5; + + //setup a move type and direction to move. + MoveLocalSc move; + VectorInt3 dir(0,1,0); + + //first let's choose the middle most monomer + move.init(myIngredients,2,dir); + float probMiddleMost=1; + + VectorInt3 bondVector1,bondVector2; + + //Initial bondvectors at monomer position. + bondVector1.setAllCoordinates(2,0,0); + bondVector2.setAllCoordinates(2,0,0); + + float probPresentMonPos=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength); + + //Future bondvectors at monomer position after the move. + bondVector1.setAllCoordinates(2,1,0); + bondVector2.setAllCoordinates(2,-1,0); + + float probFutureMonPos=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength); + + probMiddleMost *= probFutureMonPos/probPresentMonPos; + + + //Initial bondvectors at monomer position one direction positive to the moved monomer. + bondVector1.setAllCoordinates(2,0,0); + bondVector2.setAllCoordinates(2,0,0); + + float probPresentOnePositive=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength); + + bondVector1.setAllCoordinates(2,-1,0); + bondVector2.setAllCoordinates(2,0,0); + + float probFutureOnePositive=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength); + + probMiddleMost *= probFutureOnePositive/probPresentOnePositive; + + //Since the changes in bond vectors in positive and negative direction is symmetric + //we can multiply same probability(calculated for positive) for negative direction. + + probMiddleMost *= probFutureOnePositive/probPresentOnePositive; + + move.check(myIngredients); + + //use less than instead of equal to get rid of slight mistmatch problem. + EXPECT_LE(move.getProbability()-probMiddleMost,1e-8); + + + //choose the monomer second last to end of the chain. + move.init(myIngredients,3,dir); + + float probSecondLastMon = 1; + + //It will have contribution from the monomer position change. + probSecondLastMon *= probFutureMonPos/probPresentMonPos; + + //and contribution from position one positive to monomer. + probSecondLastMon *= probFutureOnePositive/probPresentOnePositive; + + move.check(myIngredients); + EXPECT_LE(move.getProbability()-probSecondLastMon,1e-8); + + //Start of chain monomer is symmetric to second last monomer to end of chain. + move.init(myIngredients,1,dir); + + move.check(myIngredients); + EXPECT_LE(move.getProbability()-probSecondLastMon,1e-8); + + //Monomer at the start of the chain. + move.init(myIngredients,0,dir); + + float probEndMon=1; + + //It has only contribution from monomer position one positive to monomer position. + probEndMon *= probFutureOnePositive/probPresentOnePositive; + + move.check(myIngredients); + EXPECT_LE(move.getProbability()-probEndMon,1e-8); + + //Similarly symmetric for the other end monomer + move.init(myIngredients,4,dir); + + move.check(myIngredients); + EXPECT_LE(move.getProbability()-probEndMon,1e-8); +} + + +TEST_F(BendingPotentialTest, CheckforMixedChain) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureBendingPotential) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + + + + //prepare myIngredients. + myIngredients.setBoxX(64); + myIngredients.setBoxY(64); + myIngredients.setBoxZ(64); + myIngredients.setPeriodicX(1); + myIngredients.setPeriodicY(1); + myIngredients.setPeriodicZ(1); + myIngredients.modifyBondset().addBFMclassicBondset(); + + //setup a mixed Linear chain. + setMixedLinearChain(myIngredients,6); + + //set lipid chainsEnds + myIngredients.setChainEnds(myIngredients); + + //set interaction for type 1 + myIngredients.setBendingPotential(myIngredients,1,0.5); + myIngredients.setBendingPotential(myIngredients,2,0.2); + double bpStrength1=0.5; + double bpStrength2=0.2; + + //setup a move type and direction to move. + MoveLocalSc move; + VectorInt3 dir(0,1,0); + + //first middle monomer in type 1 part of the chain + //This is where bit shift entry chainsEnds is useful. + move.init(myIngredients,1,dir); + float probMiddleMost1=1; + + VectorInt3 bondVector1,bondVector2; + + //Initial bondvectors at monomer position. + bondVector1.setAllCoordinates(2,0,0); + bondVector2.setAllCoordinates(2,0,0); + + float probPresentMonPos=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength1); + + //Future bondvectors at monomer position after the move. + bondVector1.setAllCoordinates(2,1,0); + bondVector2.setAllCoordinates(2,-1,0); + + float probFutureMonPos=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength1); + + probMiddleMost1 *= probFutureMonPos/probPresentMonPos; + + //There can't be any other contributions as the first chain + //only three monomers with type1 + move.check(myIngredients); + + //use less than instead of equal to get rid of slight mistmatch problem. + EXPECT_LE(move.getProbability()-probMiddleMost1,1e-8); + + //choose the monomer end of the chain for attribute 1. + move.init(myIngredients,0,dir); + + //Initial bondvectors at monomer position one direction positive to the moved monomer. + bondVector1.setAllCoordinates(2,0,0); + bondVector2.setAllCoordinates(2,0,0); + + float probPresentOnePositive=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength1); + + bondVector1.setAllCoordinates(2,-1,0); + bondVector2.setAllCoordinates(2,0,0); + + float probFutureOnePositive=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength1); + + + float probEndMon1=1; + + //It has only contribution from monomer position one positive to monomer position. + probEndMon1 *= probFutureOnePositive/probPresentOnePositive; + + move.check(myIngredients); + EXPECT_LE(move.getProbability()-probEndMon1,1e-8); + + //Similarly symmetric for the other end monomer for type 1. + move.init(myIngredients,2,dir); + + move.check(myIngredients); + EXPECT_LE(move.getProbability()-probEndMon1,1e-8); + + /********Type2********/ + + // middle monomer in Type 2 + move.init(myIngredients,4,dir); + EXPECT_EQ(myIngredients.getMolecules()[4].getAttributeTag(),2); + float probMiddleMost2=1; + + //Initial bondvectors at monomer position. + bondVector1.setAllCoordinates(2,0,0); + bondVector2.setAllCoordinates(2,0,0); + + //multiply with the bpStrength2 this time. + probPresentMonPos=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength2); + + //Future bondvectors at monomer position after the move. + bondVector1.setAllCoordinates(2,1,0); + bondVector2.setAllCoordinates(2,-1,0); + + probFutureMonPos=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength2); + + probMiddleMost2 *= probFutureMonPos/probPresentMonPos; + + move.check(myIngredients); + + //use less than instead of equal to get rid of slight mistmatch problem. + EXPECT_LE(move.getProbability()-probMiddleMost2,1e-8); + + + //choose the monomer end of the chain for attribute 2. + move.init(myIngredients,3,dir); + + //Initial bondvectors at monomer position one direction positive to the moved monomer. + bondVector1.setAllCoordinates(2,0,0); + bondVector2.setAllCoordinates(2,0,0); + + //bpStrength2 + probPresentOnePositive=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength2); + + bondVector1.setAllCoordinates(2,-1,0); + bondVector2.setAllCoordinates(2,0,0); + + probFutureOnePositive=exp(-BendingPotentials::simpleHarmonic(bondVector1,bondVector2)*bpStrength2); + + + float probEndMon2=1; + + //It has only contribution from monomer position one positive to monomer position. + probEndMon2 *= probFutureOnePositive/probPresentOnePositive; + + move.check(myIngredients); + EXPECT_LE(move.getProbability()-probEndMon2,1e-8); + + //Similarly symmetric for the other end monomer for type 2. + move.init(myIngredients,5,dir); + + move.check(myIngredients); + EXPECT_LE(move.getProbability()-probEndMon2,1e-8); + +} + +TEST_F(BendingPotentialTest, CheckReadWriteAndbondVectorToIndex) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureBondset<>,FeatureBendingPotential) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + myIngredients.modifyBondset().addBFMclassicBondset(); + + + //set interaction between type + myIngredients.setBendingPotential(myIngredients,1,0.8); + myIngredients.setBendingPotential(myIngredients,2,-0.8); + myIngredients.setBendingPotential(myIngredients,3,0.0); + myIngredients.setBendingPotential(myIngredients,9,10.0); + + //prepare myIngredients + myIngredients.setBoxX(32); + myIngredients.setBoxY(32); + myIngredients.setBoxZ(32); + myIngredients.setPeriodicX(1); + myIngredients.setPeriodicY(1); + myIngredients.setPeriodicZ(1); + + typename Ing::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("./bpStrength.test",myIngredients,AnalyzerWriteBfmFile::NEWFILE); + outfile.initialize(); + outfile.execute(); + + Ing myIngredients2; + UpdaterReadBfmFile infile("./bpStrength.test",myIngredients2,UpdaterReadBfmFile::READ_LAST_CONFIG_FAST); + infile.initialize(); + + for(int32_t i=1;i<=10;i++) + EXPECT_LE(myIngredients.getBendingPotential(i)-myIngredients2.getBendingPotential(i),1e-7); + + + outfile.closeFile(); + infile.closeFile(); + //remove the temporary file + EXPECT_EQ(0,remove("./bpStrength.test")); + + /***testing the new bondVectorToIndex function**/ + std::map ::const_iterator it; + + std::vector store(180,0); + + for(it=myIngredients.getBondset().begin();it!=myIngredients.getBondset().end();it++){ + int32_t iD=myIngredients.bondVectorToIndex(it->second); + store[iD]++;} + + int32_t count1=0; + for(int32_t i=0;i. + +--------------------------------------------------------------------------------*/ + +/*****************************************************************************/ +/** + * @file + * @brief Tests for UpdaterLipidsCreator + * */ +/*****************************************************************************/ + + +#include "gtest/gtest.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include "TesterFunctions.h" + +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 TestUpdaterLipidsCreator: 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(TestUpdaterLipidsCreator, TestSimpleBilayerSystem) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureExcludedVolumeSc< FeatureLatticePowerOfTwo >,FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + + IngredientsType ingredients; + ingredients.setBoxX(64); + ingredients.setBoxY(64); + ingredients.setBoxZ(64); + int32_t box_volume=64*64*64; + ingredients.setPeriodicX(true); + ingredients.setPeriodicY(true); + ingredients.setPeriodicZ(true); + ingredients.modifyBondset().addBFMclassicBondset(); + EXPECT_NO_THROW(ingredients.synchronize()); + + UpdaterLipidsCreator first(ingredients,2,0.1); + + int32_t defaultLipidTail=5; + int32_t lipidsLength=2*defaultLipidTail+3; + int32_t numLipids=2; + + + TesterFunctions TF; + + // first execution + EXPECT_NO_THROW(first.initialize()); + + EXPECT_NO_THROW(ingredients.synchronize()); + + EXPECT_EQ((numLipids*lipidsLength),TF.getTotalNonsolventMonomers(ingredients)); + + EXPECT_EQ(int32_t((0.1*box_volume-numLipids*lipidsLength*8)/8),TF.getTotalsolventMonomers(ingredients)); + + EXPECT_EQ(1,first.getUpperleafletLipids()); + + // check (default) monomer tags + //first lipid + EXPECT_EQ(1,ingredients.getMolecules()[0].getAttributeTag()); + EXPECT_EQ(1,ingredients.getMolecules()[2].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[3].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[10].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[12].getAttributeTag()); + + //second lipid + EXPECT_EQ(1,ingredients.getMolecules()[13].getAttributeTag()); + EXPECT_EQ(1,ingredients.getMolecules()[15].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[17].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[20].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[25].getAttributeTag()); + + //solvent + EXPECT_EQ(5, ingredients.getMolecules()[30].getAttributeTag()); + EXPECT_EQ(5, ingredients.getMolecules()[32].getAttributeTag()); + EXPECT_EQ(5, ingredients.getMolecules()[27].getAttributeTag()); + EXPECT_EQ(5, ingredients.getMolecules()[28].getAttributeTag()); + + //check connectivity + //first lipid + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(0)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(1)); + EXPECT_EQ(3, ingredients.getMolecules().getNumLinks(2)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(5)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(7)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(9)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(12)); + + //second lipid + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(13)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(14)); + EXPECT_EQ(3, ingredients.getMolecules().getNumLinks(15)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(18)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(20)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(22)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(25)); + + //solvent + EXPECT_EQ(0, ingredients.getMolecules().getNumLinks(30)); + EXPECT_EQ(0, ingredients.getMolecules().getNumLinks(32)); + EXPECT_EQ(0, ingredients.getMolecules().getNumLinks(27)); + EXPECT_EQ(0, ingredients.getMolecules().getNumLinks(28)); + + //Orientations of lipids. + std::vector Orientations = TF.getOrientations(ingredients,defaultLipidTail); + + EXPECT_EQ(2, Orientations.size()); + + EXPECT_FALSE(Orientations[0].getX()); + EXPECT_FALSE(Orientations[0].getY()); + EXPECT_TRUE(Orientations[0].getZ()); + + EXPECT_FALSE(Orientations[1].getX()); + EXPECT_FALSE(Orientations[1].getY()); + EXPECT_TRUE(Orientations[1].getZ()); + +} + +TEST_F(TestUpdaterLipidsCreator, TestMixedMembraneAndRinghead) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureExcludedVolumeSc< FeatureLatticePowerOfTwo >,FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + + IngredientsType ingredients; + ingredients.setBoxX(64); + ingredients.setBoxY(64); + ingredients.setBoxZ(64); + int32_t box_volume=64*64*64; + ingredients.setPeriodicX(true); + ingredients.setPeriodicY(true); + ingredients.setPeriodicZ(true); + ingredients.modifyBondset().addBFMclassicBondset(); + EXPECT_NO_THROW(ingredients.synchronize()); + + UpdaterLipidsCreator first_1(ingredients,4,0.1,5,7,0.5,true); + + int32_t longerLipidsLength=2*5+4; + int32_t shorterLipidsLength=2*7+4; + int32_t longerTail=7; + int32_t shorterTail=5; + //number of longerlipids. + int32_t n1=2; + //number of shorterlipids. + int32_t n2=2; + + + TesterFunctions TF; + // first execution + EXPECT_NO_THROW(first_1.initialize()); + + EXPECT_NO_THROW(ingredients.synchronize()); + + + EXPECT_EQ((n1*longerLipidsLength + n2*shorterLipidsLength),TF.getTotalNonsolventMonomers(ingredients)); + + EXPECT_EQ(int32_t((0.1*box_volume-(n1*longerLipidsLength + n2*shorterLipidsLength)*8)/8),TF.getTotalsolventMonomers(ingredients)); + + EXPECT_EQ(2,first_1.getUpperleafletLipids()); + + //check (default) monomer tags + //shorter lipids(4 head monomers for ring) + int32_t numMonHead1=TF.getNumMonWithAttribute(ingredients,1); + EXPECT_EQ(4*n2,numMonHead1); + + int32_t numMonTail1=TF.getNumMonWithAttribute(ingredients,2); + EXPECT_EQ(shorterTail*2*n2,numMonTail1); + + //longer lipids(4 head monomers for ring) + int32_t numMonHead2=TF.getNumMonWithAttribute(ingredients,3); + EXPECT_EQ(4*n1,numMonHead2); + + int32_t numMonTail2=TF.getNumMonWithAttribute(ingredients,4); + EXPECT_EQ(longerTail*2*n1,numMonTail2); + + //check connectivity + //All the ring heads with four connections. + int32_t numHeads4links=TF.getNumMonWithLinks(ingredients,4); + EXPECT_EQ(1*(n1+n2),numHeads4links); + + //All the tail ends connections + int32_t numTailsEnds=TF.getNumMonWithLinks(ingredients,1); + EXPECT_EQ(2*(n1+n2),numTailsEnds); + + + //Orientations of shorter lipids + std::vector Orientations = TF.getOrientations(ingredients,shorterTail,true); + + EXPECT_EQ(2, Orientations.size()); + + EXPECT_FALSE(Orientations[0].getX()); + EXPECT_FALSE(Orientations[0].getY()); + EXPECT_TRUE(Orientations[0].getZ()); + + EXPECT_FALSE(Orientations[1].getX()); + EXPECT_FALSE(Orientations[1].getY()); + EXPECT_TRUE(Orientations[1].getZ()); + + //Orientations of longer lipids + Orientations = TF.getOrientations(ingredients,longerTail,true); + + EXPECT_EQ(2, Orientations.size()); + + EXPECT_FALSE(Orientations[0].getX()); + EXPECT_FALSE(Orientations[0].getY()); + EXPECT_TRUE(Orientations[0].getZ()); + + EXPECT_FALSE(Orientations[1].getX()); + EXPECT_FALSE(Orientations[1].getY()); + EXPECT_TRUE(Orientations[1].getZ()); + + +} + +TEST_F(TestUpdaterLipidsCreator, TestObjectsandSolvent) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureExcludedVolumeSc< FeatureLatticePowerOfTwo >,FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + + IngredientsType ingredients; + ingredients.setBoxX(64); + ingredients.setBoxY(64); + ingredients.setBoxZ(64); + int32_t box_volume=64*64*64; + ingredients.setPeriodicX(true); + ingredients.setPeriodicY(true); + ingredients.setPeriodicZ(true); + ingredients.modifyBondset().addBFMclassicBondset(); + EXPECT_NO_THROW(ingredients.synchronize()); + + UpdaterLipidsCreator first_2(ingredients,0,0.1,0,0,0.5,false,false,1,2,5,2); + + //Triangular Nano Particle monomers. + int32_t tNP_Monomers=3; + //Facially Amphiphlic Copolymers monomers. + int32_t fAC_Monomers=5*2; + //Number of Triangular Nano Particle monomers. + int32_t num_tNP=2; + //Number of Facially Amphiphlic Copolymers monomers. + int32_t num_fAC=2; + + TesterFunctions TF; + //first execution + EXPECT_NO_THROW(first_2.initialize()); + + EXPECT_NO_THROW(ingredients.synchronize()); + + + EXPECT_EQ((num_tNP*tNP_Monomers + num_fAC*fAC_Monomers),TF.getTotalNonsolventMonomers(ingredients)); + + EXPECT_EQ(int32_t((0.1*box_volume-(num_tNP*tNP_Monomers + num_fAC*fAC_Monomers)*8)/8),TF.getTotalsolventMonomers(ingredients)); + + //Half of Facially Amphiphlic Copolymers monomers must be tag 7 and all of nano particles must be 7. + int32_t type7monomers=TF.getNumMonWithAttribute(ingredients,7); + EXPECT_EQ(num_tNP*tNP_Monomers+num_fAC*int(0.5*fAC_Monomers),type7monomers); + + int32_t type8monomers=TF.getNumMonWithAttribute(ingredients,8); + EXPECT_EQ(num_fAC*int(0.5*fAC_Monomers),type8monomers); + + // check connectivity + // Chain ends only exist in Facially Amphiphlic Copolymers*/ + int32_t numEnds=TF.getNumMonWithLinks(ingredients,1); + EXPECT_EQ(num_fAC*5,numEnds); + + //Monomer with dual connection present in FACs and all the monomers of TNPs. + int32_t num2Links=TF.getNumMonWithLinks(ingredients,2); + EXPECT_EQ(num_fAC*2+num_tNP*3,num2Links); + + //Monomer with three connections present in FACs (Copolymers_length-2*ends_monomers). + int32_t num3Links=TF.getNumMonWithLinks(ingredients,3); + EXPECT_EQ(num_fAC*(5-2),num3Links); + + /**Generally both Triangular and tetrahedral NPs are not required, however, + * as the functions for generating tetrahedral NPs is publically available we can + * generate them here as the mean lattice occupancy(MLO) is set to be quite low. + * This procedure is not reccomended for general setups where MLO is high. + */ + + //Initiate move to create nano particles. + MoveAddMonomerSc move; + first_2.nanoTetrahedron(move,2); + + //Tetrahedral Nano Particle monomers. + int32_t ttNP_Monomers=4; + //Number of Tetrahedral Nano Particle monomers. + int32_t num_ttNP=2; + + // check (default) monomer tags + //All monomer ttNP are type 7. + type7monomers=TF.getNumMonWithAttribute(ingredients,7); + EXPECT_EQ(num_tNP*tNP_Monomers+num_fAC*int(0.5*fAC_Monomers)+ttNP_Monomers*num_ttNP,type7monomers); + + //Monomer with three connections in FACs and ttNPs. + num3Links=TF.getNumMonWithLinks(ingredients,3); + EXPECT_EQ(num_fAC*(5-2)+num_ttNP*4,num3Links); + +} diff --git a/tests/updater/TestUpdaterLipidsRandomInitializer.cpp b/tests/updater/TestUpdaterLipidsRandomInitializer.cpp new file mode 100644 index 0000000..a5cf759 --- /dev/null +++ b/tests/updater/TestUpdaterLipidsRandomInitializer.cpp @@ -0,0 +1,270 @@ +/*-------------------------------------------------------------------------------- + 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 UpdaterLipidsRandomInitializer + * */ +/*****************************************************************************/ + + +#include "gtest/gtest.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include "TesterFunctions.h" + +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 TestUpdaterLipidsRandomInitializer: 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; + +}; + +//This tests are similar to tests for UpdaterLipidsCreator, the only difference +//here is that Orientations would random instead only towards z axis. +TEST_F(TestUpdaterLipidsRandomInitializer, TestSimpleRandomisedSystem) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureExcludedVolumeSc< FeatureLatticePowerOfTwo >,FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + + IngredientsType ingredients; + ingredients.setBoxX(64); + ingredients.setBoxY(64); + ingredients.setBoxZ(64); + int32_t box_volume=64*64*64; + ingredients.setPeriodicX(true); + ingredients.setPeriodicY(true); + ingredients.setPeriodicZ(true); + ingredients.modifyBondset().addBFMclassicBondset(); + EXPECT_NO_THROW(ingredients.synchronize()); + + UpdaterLipidsRandomInitializer first(ingredients,2,0.1); + + int32_t defaultLipidTail=5; + int32_t lipidsLength=2*defaultLipidTail+3; + int32_t numLipids=2; + + + TesterFunctions TF; + + // first execution + EXPECT_NO_THROW(first.initialize()); + + EXPECT_NO_THROW(ingredients.synchronize()); + + EXPECT_EQ((numLipids*lipidsLength),TF.getTotalNonsolventMonomers(ingredients)); + + EXPECT_EQ(int32_t((0.1*box_volume-numLipids*lipidsLength*8)/8),TF.getTotalsolventMonomers(ingredients)); + + + // check (default) monomer tags + // first lipid. + EXPECT_EQ(1,ingredients.getMolecules()[0].getAttributeTag()); + EXPECT_EQ(1,ingredients.getMolecules()[2].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[3].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[10].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[12].getAttributeTag()); + + // second lipid. + EXPECT_EQ(1,ingredients.getMolecules()[13].getAttributeTag()); + EXPECT_EQ(1,ingredients.getMolecules()[15].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[17].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[20].getAttributeTag()); + EXPECT_EQ(2,ingredients.getMolecules()[25].getAttributeTag()); + + // solvent . + EXPECT_EQ(5, ingredients.getMolecules()[30].getAttributeTag()); + EXPECT_EQ(5, ingredients.getMolecules()[32].getAttributeTag()); + EXPECT_EQ(5, ingredients.getMolecules()[27].getAttributeTag()); + EXPECT_EQ(5, ingredients.getMolecules()[28].getAttributeTag()); + + // check connectivity + //first lipid. + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(0)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(1)); + EXPECT_EQ(3, ingredients.getMolecules().getNumLinks(2)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(5)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(7)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(9)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(12)); + + //second lipid. + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(13)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(14)); + EXPECT_EQ(3, ingredients.getMolecules().getNumLinks(15)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(18)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(20)); + EXPECT_EQ(2, ingredients.getMolecules().getNumLinks(22)); + EXPECT_EQ(1, ingredients.getMolecules().getNumLinks(25)); + + //solvent. + EXPECT_EQ(0, ingredients.getMolecules().getNumLinks(30)); + EXPECT_EQ(0, ingredients.getMolecules().getNumLinks(32)); + EXPECT_EQ(0, ingredients.getMolecules().getNumLinks(27)); + EXPECT_EQ(0, ingredients.getMolecules().getNumLinks(28)); + + //Orientations of lipids. + std::vector Orientations = TF.getOrientations(ingredients,defaultLipidTail); + + EXPECT_EQ(2, Orientations.size()); + + //Different from Lipids creator. + EXPECT_GE(fabs(Orientations[0].getX()),0); + EXPECT_GE(fabs(Orientations[0].getY()),0); + EXPECT_GE(fabs(Orientations[0].getZ()),0); + + EXPECT_GE(fabs(Orientations[1].getX()),0); + EXPECT_GE(fabs(Orientations[1].getY()),0); + EXPECT_GE(fabs(Orientations[1].getZ()),0); + +} + +TEST_F(TestUpdaterLipidsRandomInitializer, TestMixedlipids) +{ + typedef LOKI_TYPELIST_3(FeatureMoleculesIO, FeatureExcludedVolumeSc< FeatureLatticePowerOfTwo >,FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients IngredientsType; + + IngredientsType ingredients; + ingredients.setBoxX(64); + ingredients.setBoxY(64); + ingredients.setBoxZ(64); + int32_t box_volume=64*64*64; + ingredients.setPeriodicX(true); + ingredients.setPeriodicY(true); + ingredients.setPeriodicZ(true); + ingredients.modifyBondset().addBFMclassicBondset(); + EXPECT_NO_THROW(ingredients.synchronize()); + + UpdaterLipidsRandomInitializer first_1(ingredients,4,0.1,5,7,0.5); + + int32_t longerLipidsLength=2*5+3; + int32_t shorterLipidsLength=2*7+3; + int32_t longerTail=7; + int32_t shorterTail=5; + //number of longerlipids. + int32_t n1=2; + //number of shorterlipids. + int32_t n2=2; + + + TesterFunctions TF; + // first execution + EXPECT_NO_THROW(first_1.initialize()); + + EXPECT_NO_THROW(ingredients.synchronize()); + + + EXPECT_EQ((n1*longerLipidsLength + n2*shorterLipidsLength),TF.getTotalNonsolventMonomers(ingredients)); + + EXPECT_EQ(int32_t((0.1*box_volume-(n1*longerLipidsLength + n2*shorterLipidsLength)*8)/8),TF.getTotalsolventMonomers(ingredients)); + + //check (default) monomer tags + //shorter lipids(3 head monomers here). + int32_t numMonHead1=TF.getNumMonWithAttribute(ingredients,1); + EXPECT_EQ(3*n2,numMonHead1); + + int32_t numMonTail1=TF.getNumMonWithAttribute(ingredients,2); + EXPECT_EQ(shorterTail*2*n2,numMonTail1); + + //longer lipids(3 head monomers). + int32_t numMonHead2=TF.getNumMonWithAttribute(ingredients,3); + EXPECT_EQ(3*n1,numMonHead2); + + int32_t numMonTail2=TF.getNumMonWithAttribute(ingredients,4); + EXPECT_EQ(longerTail*2*n1,numMonTail2); + + //check connectivity + //All the heads with three connections. + int32_t numHeads3links=TF.getNumMonWithLinks(ingredients,3); + EXPECT_EQ(1*(n1+n2),numHeads3links); + + //Two tail ends connections plus one head monomer. + int32_t numTailsEnds=TF.getNumMonWithLinks(ingredients,1); + EXPECT_EQ((2+1)*(n1+n2),numTailsEnds); + + + //Orientations of shorter lipids. + std::vector Orientations = TF.getOrientations(ingredients,shorterTail,false); + + EXPECT_EQ(2, Orientations.size()); + + EXPECT_GE(fabs(Orientations[0].getX()),0); + EXPECT_GE(fabs(Orientations[0].getY()),0); + EXPECT_GE(fabs(Orientations[0].getZ()),0); + + EXPECT_GE(fabs(Orientations[1].getX()),0); + EXPECT_GE(fabs(Orientations[1].getY()),0); + EXPECT_GE(fabs(Orientations[1].getZ()),0); + + //Orientations of longer lipids. + Orientations = TF.getOrientations(ingredients,longerTail,false); + + EXPECT_EQ(2, Orientations.size()); + + EXPECT_GE(fabs(Orientations[0].getX()),0); + EXPECT_GE(fabs(Orientations[0].getY()),0); + EXPECT_GE(fabs(Orientations[0].getZ()),0); + + EXPECT_GE(fabs(Orientations[1].getX()),0); + EXPECT_GE(fabs(Orientations[1].getY()),0); + EXPECT_GE(fabs(Orientations[1].getZ()),0); +} diff --git a/tests/updater/TesterHelperFunctionsLipids.h b/tests/updater/TesterHelperFunctionsLipids.h new file mode 100644 index 0000000..baae913 --- /dev/null +++ b/tests/updater/TesterHelperFunctionsLipids.h @@ -0,0 +1,76 @@ + +template +class TesterFunctions +{ +public: + //!Function to get total number of monomers in system with a given attribute. + int32_t getNumMonWithAttribute(IngredientsType& ingredients, int32_t attribute ){ + + int32_t numMon=0; + for(int32_t i=0;i getOrientations(IngredientsType& ingredients, int32_t tail_length, bool ringHead=false){ + + std::vector Orientations; + int headlinks; + if(ringHead) + headlinks=4; + else + headlinks=3; + + for(int32_t i=0;i Date: Wed, 21 Aug 2019 23:02:52 +0200 Subject: [PATCH 2/9] lipidbilayerrelated --- tests/feature/TestFeatureBendingPotential.cpp | 10 ++++++++++ tests/updater/TestUpdaterLipidsCreator.cpp | 5 ++++- tests/updater/TestUpdaterLipidsRandomInitializer.cpp | 5 ++++- 3 files changed, 18 insertions(+), 2 deletions(-) diff --git a/tests/feature/TestFeatureBendingPotential.cpp b/tests/feature/TestFeatureBendingPotential.cpp index b2e1563..55ff43f 100644 --- a/tests/feature/TestFeatureBendingPotential.cpp +++ b/tests/feature/TestFeatureBendingPotential.cpp @@ -11,6 +11,16 @@ #include using namespace std; +/*****************************************************************************/ +/** + * @file + * @brief Tests for the class FeatureBendingPotential + * + * @author Ankush Checkervarty + * @date 01.08.2019 + * */ +/*****************************************************************************/ + class BendingPotentialTest: public ::testing::Test{ diff --git a/tests/updater/TestUpdaterLipidsCreator.cpp b/tests/updater/TestUpdaterLipidsCreator.cpp index 3894742..74e4f52 100644 --- a/tests/updater/TestUpdaterLipidsCreator.cpp +++ b/tests/updater/TestUpdaterLipidsCreator.cpp @@ -28,7 +28,10 @@ along with LeMonADE. If not, see . /*****************************************************************************/ /** * @file - * @brief Tests for UpdaterLipidsCreator + * @brief Tests for the class UpdaterLipidsCreator + * + * @author Ankush Checkervarty + * @date 01.08.2019 * */ /*****************************************************************************/ diff --git a/tests/updater/TestUpdaterLipidsRandomInitializer.cpp b/tests/updater/TestUpdaterLipidsRandomInitializer.cpp index a5cf759..e62a4f1 100644 --- a/tests/updater/TestUpdaterLipidsRandomInitializer.cpp +++ b/tests/updater/TestUpdaterLipidsRandomInitializer.cpp @@ -28,7 +28,10 @@ along with LeMonADE. If not, see . /*****************************************************************************/ /** * @file - * @brief Tests for UpdaterLipidsRandomInitializer + * @brief Tests for the class UpdaterLipidsRandomInitializer + * + * @author Ankush Checkervarty + * @date 01.08.2019 * */ /*****************************************************************************/ From eb4344af6db4ef24b095be9126337248ddd631af Mon Sep 17 00:00:00 2001 From: Ankush Check Date: Thu, 11 Mar 2021 13:34:56 +0100 Subject: [PATCH 3/9] bugfix --- CMakeLists.txt | 2 +- tests/updater/TestUpdaterLipidsCreator.cpp | 2 +- tests/updater/TestUpdaterLipidsRandomInitializer.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index eb0d436..adaa615 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -102,7 +102,7 @@ ADD_SUBDIRECTORY(projects) # # Add option for building tests # -option(LEMONADE_TESTS "Build the test" OFF) +option(LEMONADE_TESTS "Build the test" ON) if(LEMONADE_TESTS) add_subdirectory(tests) endif(LEMONADE_TESTS) diff --git a/tests/updater/TestUpdaterLipidsCreator.cpp b/tests/updater/TestUpdaterLipidsCreator.cpp index 74e4f52..99bf9ef 100644 --- a/tests/updater/TestUpdaterLipidsCreator.cpp +++ b/tests/updater/TestUpdaterLipidsCreator.cpp @@ -47,7 +47,7 @@ along with LeMonADE. If not, see . #include #include #include -#include "TesterFunctions.h" +#include "TesterHelperFunctionsLipids.h" using namespace std; diff --git a/tests/updater/TestUpdaterLipidsRandomInitializer.cpp b/tests/updater/TestUpdaterLipidsRandomInitializer.cpp index e62a4f1..3b44b4a 100644 --- a/tests/updater/TestUpdaterLipidsRandomInitializer.cpp +++ b/tests/updater/TestUpdaterLipidsRandomInitializer.cpp @@ -47,7 +47,7 @@ along with LeMonADE. If not, see . #include #include #include -#include "TesterFunctions.h" +#include "TesterHelperFunctionsLipids.h" using namespace std; From 4038946b48dd7568f4c33e3551ac7da73116fb5a Mon Sep 17 00:00:00 2001 From: bondoki Date: Thu, 20 May 2021 19:59:29 +0200 Subject: [PATCH 4/9] * change MoveConnectScReactive for selecting only reactive monomers inside shell (+-2,0,0),(+-2,+-1,0),(+-2,+-1,+-1) by random chance * improving tests for TestMoveConnectionScReactive --- .../updater/moves/MoveConnectScReactive.h | 182 +++++++++++++++--- .../updater/TestMoveConnectionScReactive.cpp | 37 +++- 2 files changed, 193 insertions(+), 26 deletions(-) diff --git a/include/LeMonADE/updater/moves/MoveConnectScReactive.h b/include/LeMonADE/updater/moves/MoveConnectScReactive.h index fa4fd96..352abf4 100644 --- a/include/LeMonADE/updater/moves/MoveConnectScReactive.h +++ b/include/LeMonADE/updater/moves/MoveConnectScReactive.h @@ -3,7 +3,7 @@ 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 D evelopment | Copyright (C) 2013-2015,2021 by o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) ooo | ---------------------------------------------------------------------------------- @@ -46,12 +46,72 @@ class MoveConnectScReactive:public MoveConnectBase { public: MoveConnectScReactive(){ - shellPositions[0]=VectorInt3( 2, 0, 0); + // old implementation until May2021 + /*shellPositions[0]=VectorInt3( 2, 0, 0); shellPositions[1]=VectorInt3(-2, 0, 0); shellPositions[2]=VectorInt3( 0, 2, 0); shellPositions[3]=VectorInt3( 0,-2, 0); shellPositions[4]=VectorInt3( 0, 0, 2); shellPositions[5]=VectorInt3( 0, 0,-2); + */ + + shellPositions[0] = VectorInt3(2,0,0); + shellPositions[1] = VectorInt3(0,0,2); + shellPositions[2] = VectorInt3(0,2,0); + shellPositions[3] = VectorInt3(-2,0,0); + shellPositions[4] = VectorInt3(0,0,-2); + shellPositions[5] = VectorInt3(0,-2,0); + + shellPositions[6] = VectorInt3(2,1,0); + shellPositions[7] = VectorInt3(1,0,2); + shellPositions[8] = VectorInt3(0,2,1); + shellPositions[9] = VectorInt3(-2,1,0); + shellPositions[10] = VectorInt3(1,0,-2); + shellPositions[11] = VectorInt3(0,-2,1); + shellPositions[12] = VectorInt3(2,-1,0); + shellPositions[13] = VectorInt3(-1,0,2); + shellPositions[14] = VectorInt3(0,2,-1); + shellPositions[15] = VectorInt3(-2,-1,0); + shellPositions[16] = VectorInt3(-1,0,-2); + shellPositions[17] = VectorInt3(0,-2,-1); + shellPositions[18] = VectorInt3(1,2,0); + shellPositions[19] = VectorInt3(2,0,1); + shellPositions[20] = VectorInt3(0,1,2); + shellPositions[21] = VectorInt3(-1,2,0); + shellPositions[22] = VectorInt3(2,0,-1); + shellPositions[23] = VectorInt3(0,-1,2); + shellPositions[24] = VectorInt3(1,-2,0); + shellPositions[25] = VectorInt3(-2,0,1); + shellPositions[26] = VectorInt3(0,1,-2); + shellPositions[27] = VectorInt3(-1,-2,0); + shellPositions[28] = VectorInt3(-2,0,-1); + shellPositions[29] = VectorInt3(0,-1,-2); + + shellPositions[30] = VectorInt3(2,1,1); + shellPositions[31] = VectorInt3(1,1,2); + shellPositions[32] = VectorInt3(1,2,1); + shellPositions[33] = VectorInt3(-2,1,1); + shellPositions[34] = VectorInt3(1,1,-2); + shellPositions[35] = VectorInt3(1,-2,1); + shellPositions[36] = VectorInt3(2,-1,1); + shellPositions[37] = VectorInt3(-1,1,2); + shellPositions[38] = VectorInt3(1,2,-1); + shellPositions[39] = VectorInt3(-2,-1,1); + shellPositions[40] = VectorInt3(-1,1,-2); + shellPositions[41] = VectorInt3(1,-2,-1); + shellPositions[42] = VectorInt3(2,1,-1); + shellPositions[43] = VectorInt3(1,-1,2); + shellPositions[44] = VectorInt3(-1,2,1); + shellPositions[45] = VectorInt3(-2,1,-1); + shellPositions[46] = VectorInt3(1,-1,-2); + shellPositions[47] = VectorInt3(-1,-2,1); + shellPositions[48] = VectorInt3(2,-1,-1); + shellPositions[49] = VectorInt3(-1,-1,2); + shellPositions[50] = VectorInt3(-1,2,-1); + shellPositions[51] = VectorInt3(-2,-1,-1); + shellPositions[52] = VectorInt3(-1,-1,-2); + shellPositions[53] = VectorInt3(-1,-2,-1); + } // overload initialise function to be able to set the moves index and direction if neccessary @@ -66,18 +126,22 @@ class MoveConnectScReactive:public MoveConnectBase private: // holds the possible move directions /** - * @brief Array that holds the 6 possible move directions + * @brief Array that holds the 54 possible move directions * * @details * * shellPositions = (dx, dy, dz) - * * shellPositions[0]= ( 2, 0, 0); - * * shellPositions[1]= (-2, 0, 0); - * * shellPositions[2]= ( 0, 2, 0); - * * shellPositions[3]= ( 0, -2, 0); - * * shellPositions[4]= ( 0, 0, 2); - * * shellPositions[5]= ( 0, 0, -2); + * * shellPositions[..]=(+-2, 0, 0); // 6 + * * shellPositions[..]=(+-2,+-1, 0); // 24 + * * shellPositions[..]=(+-2,+-1,+-1); // 24 */ - VectorInt3 shellPositions[6]; + VectorInt3 shellPositions[54]; + const size_t shellPositionsEntries = 54; + + //! number of randomly choosen direction out of the shell + const size_t numRandomSelectedDirections = 1; + + //! contains the id of the reactive partner in the randomly selected shell direction + std::map IdReactivePartnerInShell; }; @@ -98,6 +162,9 @@ template void MoveConnectScReactive::init(const IngredientsType& ing) { this->resetProbability(); + + // clear the list of reactive partners in shell + IdReactivePartnerInShell.clear(); //draw index auto nUnreactedMonomers(ing.getNUnreactedMonomers()); @@ -106,14 +173,48 @@ void MoveConnectScReactive::init(const IngredientsType& ing) this->setPartner(std::numeric_limits::max() ); } else { + // randomly select a monomer auto index(this->randomNumbers.r250_rand32() % nUnreactedMonomers); auto it (ing.getUnreactiveMonomers().cbegin()); std::advance( it, index); this->setIndex(it->first); - //draw direction - VectorInt3 randomDir(shellPositions[ this->randomNumbers.r250_rand32() % 6]); - this->setDir(randomDir); - this->setPartner(ing.getIdFromLattice(ing.getMolecules()[this->getIndex()]+randomDir) ); + + // search the vicinity for reactive monomers + for(size_t test = 0; test < numRandomSelectedDirections; test++) + { + //draw direction + VectorInt3 randomDir(shellPositions[ this->randomNumbers.r250_rand32() % shellPositionsEntries]); + + // Get the lattice value at a certain point - see FeatureConnectionSc + // return monomer index between (0; molecules.size()-1) if there's a monomer + // return uint32_t(-1)=4294967295 if place is empty + uint32_t idPartner=ing.getIdFromLattice(ing.getMolecules()[this->getIndex()]+randomDir); + + // check if selected position contains a partner + if ((std::numeric_limits::max() != idPartner ) && ing.getMolecules()[idPartner].isReactive() ) + { + IdReactivePartnerInShell[idPartner]=randomDir; + } + } + + // select a random partner and direction for connection move + if(IdReactivePartnerInShell.size() == 0) + { + // no partner at all - move will rejected anyway + this->setDir(shellPositions[0]); // it's arbitary + this->setPartner(std::numeric_limits::max() ); + } + else + { + // select a random reactive monomer from the list + auto idx(this->randomNumbers.r250_rand32() % IdReactivePartnerInShell.size()); + auto itID (IdReactivePartnerInShell.cbegin()); + std::advance( itID, idx); + + this->setDir(itID->second); + this->setPartner(itID->first); + } + } } @@ -130,6 +231,9 @@ template void MoveConnectScReactive::init(const IngredientsType& ing, uint32_t index) { this->resetProbability(); + + // clear the list of reactive partners in shell + IdReactivePartnerInShell.clear(); if ( !ing.checkCapableFormingBonds(index) ) { @@ -143,10 +247,41 @@ void MoveConnectScReactive::init(const IngredientsType& ing, uint32_t index) else throw std::runtime_error("MoveConnectScReactive::init(ing, index): index out of range!"); - //draw direction - VectorInt3 randomDir(shellPositions[ this->randomNumbers.r250_rand32() % 6]); - this->setDir(randomDir); - this->setPartner(ing.getIdFromLattice(ing.getMolecules()[this->getIndex()]+randomDir)); + // search the vicinity for reactive monomers + for(size_t test = 0; test < numRandomSelectedDirections; test++) + { + //draw direction + VectorInt3 randomDir(shellPositions[ this->randomNumbers.r250_rand32() % shellPositionsEntries]); + + // Get the lattice value at a certain point - see FeatureConnectionSc + // return monomer index between (0; molecules.size()-1) if there's a monomer + // return uint32_t(-1)=4294967295 if place is empty + uint32_t idPartner=ing.getIdFromLattice(ing.getMolecules()[this->getIndex()]+randomDir); + + // check if selected position contains a partner + if ((std::numeric_limits::max() != idPartner ) && ing.getMolecules()[idPartner].isReactive() ) + { + IdReactivePartnerInShell[idPartner]=randomDir; + } + } + + // select a random partner and direction for connection move + if(IdReactivePartnerInShell.size() == 0) + { + // no partner at all - move will rejected anyway + this->setDir(shellPositions[0]); // it's arbitary + this->setPartner(std::numeric_limits::max() ); + } + else + { + // select a random reactive monomer from the list + auto idx(this->randomNumbers.r250_rand32() % IdReactivePartnerInShell.size()); + auto itID (IdReactivePartnerInShell.cbegin()); + std::advance( itID, idx); + + this->setDir(itID->second); + this->setPartner(itID->first); + } } } @@ -210,14 +345,11 @@ void MoveConnectScReactive::init(const IngredientsType& ing, uint32_t index, Vec else throw std::runtime_error("MoveConnectScReactive::init(ing, index, bondpartner): index out of range!"); - //set direction - if(dir==shellPositions[0] || - dir==shellPositions[1] || - dir==shellPositions[2] || - dir==shellPositions[3] || - dir==shellPositions[4] || - dir==shellPositions[5] ) + //set direction if applicable + if (std::find(std::begin(shellPositions), std::end(shellPositions), dir) != std::end(shellPositions)) + { this->setDir(dir); + } else throw std::runtime_error("MoveLocalSc::init(ing, dir): direction vector out of range!"); this->setPartner(ing.getIdFromLattice(ing.getMolecules()[this->getIndex()]+dir)); diff --git a/tests/updater/TestMoveConnectionScReactive.cpp b/tests/updater/TestMoveConnectionScReactive.cpp index 2642dd3..8914a6d 100644 --- a/tests/updater/TestMoveConnectionScReactive.cpp +++ b/tests/updater/TestMoveConnectionScReactive.cpp @@ -3,7 +3,7 @@ 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 D evelopment | Copyright (C) 2013-2015,2021 by o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) ooo | ---------------------------------------------------------------------------------- @@ -98,18 +98,46 @@ TEST_F(TestMoveConnectionScReactive, checkAll) MoveConnectScReactive move; + move.init(ingredients); + EXPECT_TRUE((move.getIndex() >= 0) && (move.getIndex() <= 1)); // only 0 and 1 are reactive + EXPECT_TRUE((move.getDir()*move.getDir() >= 4) && (move.getDir()*move.getDir() <= 6)); + // will maybe fail as the direction is selected randomly hence the partner also + //EXPECT_TRUE((move.getPartner() >= 0) && (move.getPartner() <= 1)); // only 0 and 1 are reactive + EXPECT_TRUE(((move.getPartner() >= 0) && (move.getPartner() <= 1)) || (move.getPartner()==std::numeric_limits::max()) ); // only 0 and 1 are reactive + EXPECT_FALSE((move.getPartner() >= 2) && (move.getPartner() < std::numeric_limits::max()) ); // only 0 and 1 are reactive + + move.init(ingredients,0); EXPECT_EQ(0, move.getIndex()); + EXPECT_TRUE((move.getIndex() >= 0) && (move.getIndex() <= 1)); // only 0 and 1 are reactive + EXPECT_TRUE((move.getDir()*move.getDir() >= 4) && (move.getDir()*move.getDir() <= 6)); + EXPECT_TRUE(((move.getPartner() == 0)) || (move.getPartner()==std::numeric_limits::max()) ); // only 0 and 1 are reactive + EXPECT_FALSE((move.getPartner() >= 2) && (move.getPartner() < std::numeric_limits::max()) ); // only 0 and 1 are reactive + move.init(ingredients,1); EXPECT_EQ(1, move.getIndex()); + EXPECT_TRUE((move.getIndex() >= 0) && (move.getIndex() <= 1)); // only 0 and 1 are reactive + EXPECT_TRUE((move.getDir()*move.getDir() >= 4) && (move.getDir()*move.getDir() <= 6)); + EXPECT_TRUE(((move.getPartner() == 0)) || (move.getPartner()==std::numeric_limits::max()) ); // only 0 and 1 are reactive + EXPECT_FALSE((move.getPartner() >= 2) && (move.getPartner() < std::numeric_limits::max()) ); // only 0 and 1 are reactive + move.init(ingredients,1,2); EXPECT_EQ(move.getDir().getY(),2); EXPECT_FALSE(move.check(ingredients)); + EXPECT_TRUE(move.getIndex() == 1); // only 0 and 1 are reactive + EXPECT_TRUE((move.getDir()*move.getDir() == 4)); + EXPECT_TRUE((move.getPartner() == 2)); // only 0 and 1 are reactive + EXPECT_FALSE((move.getPartner() < 2) && (move.getPartner() < std::numeric_limits::max()) ); // only 0 and 1 are reactive + move.init(ingredients,0,1); EXPECT_EQ(move.getDir().getX(),2); + EXPECT_TRUE(move.getIndex() == 0); // only 0 and 1 are reactive + EXPECT_TRUE((move.getDir()*move.getDir() == 4)); + EXPECT_TRUE((move.getPartner() == 1)); // only 0 and 1 are reactive + EXPECT_FALSE((move.getPartner() < 1) && (move.getPartner() < std::numeric_limits::max()) ); // only 0 and 1 are reactive EXPECT_TRUE(move.check(ingredients)); move.apply(ingredients); @@ -118,6 +146,13 @@ TEST_F(TestMoveConnectionScReactive, checkAll) move.init(ingredients,0,1); //already occupied EXPECT_FALSE(move.check(ingredients)); + //no reactive groups left + EXPECT_EQ(ingredients.getNUnreactedMonomers(), 0); + move.init(ingredients); + EXPECT_EQ(move.getIndex(), 0 ); // default value for invalid move + EXPECT_TRUE(move.getPartner() == std::numeric_limits::max()); + EXPECT_FALSE(move.check(ingredients)); + } From ec4acd6c59dbad9b69620ef9e5d87dab7cadebfa Mon Sep 17 00:00:00 2001 From: bondoki Date: Fri, 21 May 2021 10:25:42 +0200 Subject: [PATCH 5/9] * improve test for MoveConnectionScReactive --- tests/updater/TestMoveConnectionScReactive.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/updater/TestMoveConnectionScReactive.cpp b/tests/updater/TestMoveConnectionScReactive.cpp index 8914a6d..e060df3 100644 --- a/tests/updater/TestMoveConnectionScReactive.cpp +++ b/tests/updater/TestMoveConnectionScReactive.cpp @@ -146,11 +146,12 @@ TEST_F(TestMoveConnectionScReactive, checkAll) move.init(ingredients,0,1); //already occupied EXPECT_FALSE(move.check(ingredients)); - //no reactive groups left - EXPECT_EQ(ingredients.getNUnreactedMonomers(), 0); + //one reactive group left (monomer id 1) + EXPECT_EQ(ingredients.getNUnreactedMonomers(), 1); move.init(ingredients); - EXPECT_EQ(move.getIndex(), 0 ); // default value for invalid move - EXPECT_TRUE(move.getPartner() == std::numeric_limits::max()); + EXPECT_EQ(move.getIndex(), 1 ); + EXPECT_TRUE(move.getPartner() == std::numeric_limits::max()); // default value for invalid move + EXPECT_TRUE((move.getDir()*move.getDir() == 4)); // default value for invalid move EXPECT_FALSE(move.check(ingredients)); From 48baf1fee11c76b1906d82a215922bed9280dded Mon Sep 17 00:00:00 2001 From: bondoki Date: Fri, 21 May 2021 11:02:48 +0200 Subject: [PATCH 6/9] * MoveConnectScReactive + check for possible connection inside shell and for vacant reactive monomer + change default behavior for invalid MoveConnectScReactive.init() --- include/LeMonADE/updater/moves/MoveConnectScReactive.h | 6 +++--- tests/updater/TestMoveConnectionScReactive.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/LeMonADE/updater/moves/MoveConnectScReactive.h b/include/LeMonADE/updater/moves/MoveConnectScReactive.h index 352abf4..d1359df 100644 --- a/include/LeMonADE/updater/moves/MoveConnectScReactive.h +++ b/include/LeMonADE/updater/moves/MoveConnectScReactive.h @@ -169,7 +169,7 @@ void MoveConnectScReactive::init(const IngredientsType& ing) //draw index auto nUnreactedMonomers(ing.getNUnreactedMonomers()); if( nUnreactedMonomers ==0 ) { - this->setIndex(0); + this->setIndex(std::numeric_limits::max()); this->setPartner(std::numeric_limits::max() ); } else { @@ -190,8 +190,8 @@ void MoveConnectScReactive::init(const IngredientsType& ing) // return uint32_t(-1)=4294967295 if place is empty uint32_t idPartner=ing.getIdFromLattice(ing.getMolecules()[this->getIndex()]+randomDir); - // check if selected position contains a partner - if ((std::numeric_limits::max() != idPartner ) && ing.getMolecules()[idPartner].isReactive() ) + // check if selected position contains a reactive partner with vaccant reaction site + if ((std::numeric_limits::max() != idPartner ) && (ing.getMolecules()[idPartner].isReactive()) && (ing.getMolecules().getNumLinks(idPartner) < ing.getMolecules()[idPartner].getNumMaxLinks()) ) { IdReactivePartnerInShell[idPartner]=randomDir; } diff --git a/tests/updater/TestMoveConnectionScReactive.cpp b/tests/updater/TestMoveConnectionScReactive.cpp index e060df3..9292849 100644 --- a/tests/updater/TestMoveConnectionScReactive.cpp +++ b/tests/updater/TestMoveConnectionScReactive.cpp @@ -111,7 +111,7 @@ TEST_F(TestMoveConnectionScReactive, checkAll) EXPECT_EQ(0, move.getIndex()); EXPECT_TRUE((move.getIndex() >= 0) && (move.getIndex() <= 1)); // only 0 and 1 are reactive EXPECT_TRUE((move.getDir()*move.getDir() >= 4) && (move.getDir()*move.getDir() <= 6)); - EXPECT_TRUE(((move.getPartner() == 0)) || (move.getPartner()==std::numeric_limits::max()) ); // only 0 and 1 are reactive + EXPECT_TRUE(((move.getPartner() == 1)) || (move.getPartner()==std::numeric_limits::max()) ); // only 0 and 1 are reactive EXPECT_FALSE((move.getPartner() >= 2) && (move.getPartner() < std::numeric_limits::max()) ); // only 0 and 1 are reactive From cc91c8eb991e1c5921f613d4ee63cb6353b6d526 Mon Sep 17 00:00:00 2001 From: bondoki Date: Thu, 27 May 2021 11:35:41 +0200 Subject: [PATCH 7/9] * Hotfix for ResultFormattingTools + floating point types (float, double) is printed with maximum precision to avoid rounding errors --- include/LeMonADE/utility/ResultFormattingTools.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/LeMonADE/utility/ResultFormattingTools.h b/include/LeMonADE/utility/ResultFormattingTools.h index da2a1cf..9ce2cff 100644 --- a/include/LeMonADE/utility/ResultFormattingTools.h +++ b/include/LeMonADE/utility/ResultFormattingTools.h @@ -35,7 +35,8 @@ along with LeMonADE. If not, see . #include #include #include - +#include +#include /** @@ -120,6 +121,8 @@ void ResultFormattingTools::writeTable(std::ostream& stream, for (size_t row = 0; row < results[0].size(); ++row) { for (size_t column = 0; column < results.size(); ++column) { + stream.precision(std::numeric_limits::max_digits10); + stream.setf( std::ios::fixed, std:: ios::floatfield ); stream << results[column][row] << "\t"; } stream << std::endl; From a4c2d8b0f14cf8c6ad23634b169f7fc0c597b4cf Mon Sep 17 00:00:00 2001 From: bondoki Date: Thu, 27 May 2021 11:35:41 +0200 Subject: [PATCH 8/9] * Hotfix for ResultFormattingTools + floating point types (float, double) is printed with maximum precision to avoid rounding errors + change output for writeTable() and appendToResultFile() --- include/LeMonADE/utility/ResultFormattingTools.h | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/include/LeMonADE/utility/ResultFormattingTools.h b/include/LeMonADE/utility/ResultFormattingTools.h index da2a1cf..db635fe 100644 --- a/include/LeMonADE/utility/ResultFormattingTools.h +++ b/include/LeMonADE/utility/ResultFormattingTools.h @@ -35,7 +35,8 @@ along with LeMonADE. If not, see . #include #include #include - +#include +#include /** @@ -120,6 +121,8 @@ void ResultFormattingTools::writeTable(std::ostream& stream, for (size_t row = 0; row < results[0].size(); ++row) { for (size_t column = 0; column < results.size(); ++column) { + stream.precision(std::numeric_limits::max_digits10); + stream.setf( std::ios::fixed, std:: ios::floatfield ); stream << results[column][row] << "\t"; } stream << std::endl; @@ -177,6 +180,8 @@ void ResultFormattingTools::appendToResultFile(std::string filename,ResultType& //write content for (size_t row = 0; row < results[0].size(); ++row) { for (size_t column = 0; column < results.size(); ++column) { + file.precision(std::numeric_limits::max_digits10); + file.setf( std::ios::fixed, std:: ios::floatfield ); file << results[column][row] << "\t"; } file << std::endl; From 083ae421c227c967639f3b445efa1da6ab1a51f6 Mon Sep 17 00:00:00 2001 From: bondoki Date: Mon, 31 May 2021 14:06:38 +0200 Subject: [PATCH 9/9] * changing test to fit new high precision writing of ResultFormatingTools --- .../TestAnalyzerPermeabilityCalculator.cpp | 33 +++-- tests/analyzer/TestAnalyzerPoreFinder.cpp | 124 +++++++++++++----- 2 files changed, 117 insertions(+), 40 deletions(-) diff --git a/tests/analyzer/TestAnalyzerPermeabilityCalculator.cpp b/tests/analyzer/TestAnalyzerPermeabilityCalculator.cpp index f13a4ef..925c5bd 100644 --- a/tests/analyzer/TestAnalyzerPermeabilityCalculator.cpp +++ b/tests/analyzer/TestAnalyzerPermeabilityCalculator.cpp @@ -31,7 +31,8 @@ * @brief Tests for the class AnalyzerPermeabilityCalculator * * @author Ankush Checkervarty - * @date 09.08.2019 + * @author Ron Dockhorn + * @date 31.05.2021 * */ /*****************************************************************************/ @@ -39,6 +40,7 @@ #include #include +#include #include #include @@ -280,9 +282,16 @@ TEST_F(AnalyzerPermeabilityCalculatorTest, CheckFileDumpandMidplaneUpdate) EXPECT_TRUE(file.is_open()); std::string linecontent = passComments(file); - EXPECT_EQ(linecontent[0],'0'); - EXPECT_EQ(linecontent[2],'1'); - EXPECT_EQ(linecontent[4],'1'); + std::string buf; // Have a buffer string + std::stringstream ss(linecontent); // Insert the string into a stream + std::vector tokens; // Create vector to hold our words + + while (ss >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),0.0); + EXPECT_EQ(atof(tokens[1].c_str()),1.0); + EXPECT_EQ(atof(tokens[2].c_str()),1.0); pmAnalyzer1.counterObjects=2; pmAnalyzer1.counterSolvent=3; @@ -292,9 +301,15 @@ TEST_F(AnalyzerPermeabilityCalculatorTest, CheckFileDumpandMidplaneUpdate) pmAnalyzer1.execute(); std::getline(file,linecontent); - EXPECT_EQ(linecontent[0],'1'); - EXPECT_EQ(linecontent[2],'2'); - EXPECT_EQ(linecontent[4],'3'); + std::stringstream ss2(linecontent); // Insert the string into a stream + tokens.clear(); + + while (ss2 >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),1.0); + EXPECT_EQ(atof(tokens[1].c_str()),2.0); + EXPECT_EQ(atof(tokens[2].c_str()),3.0); remove("TestDump.dat"); //midplane update Check @@ -343,7 +358,7 @@ TEST_F(AnalyzerPermeabilityCalculatorTest, CheckPoreUpdate) width=4; setConfig(midplane,width); - std::string outputFilename="TestDump.dat"; + std::string outputFilename="TestDumpA.dat"; PmAnalyzerDerived pmAnalyzer4(ingredients,objects,16,3,outputFilename,100,1,1); EXPECT_EQ(midplane, pmAnalyzer4.midplane[0][0]); @@ -368,6 +383,6 @@ TEST_F(AnalyzerPermeabilityCalculatorTest, CheckPoreUpdate) EXPECT_EQ(midplane, pmAnalyzer4.midplane[2][1]); EXPECT_EQ(midplane, pmAnalyzer4.midplane[3][0]); EXPECT_EQ(midplane, pmAnalyzer4.midplane[3][1]); - remove("TestDump.dat"); + remove("TestDumpA.dat"); } diff --git a/tests/analyzer/TestAnalyzerPoreFinder.cpp b/tests/analyzer/TestAnalyzerPoreFinder.cpp index dbe8f1b..de724e6 100644 --- a/tests/analyzer/TestAnalyzerPoreFinder.cpp +++ b/tests/analyzer/TestAnalyzerPoreFinder.cpp @@ -31,7 +31,8 @@ * @brief Tests for the class AnalyzerPoreFinder * * @author Ankush Checkervarty - * @date 09.08.2019 + * @author Ron Dockhorn + * @date 31.05.2021 * */ /*****************************************************************************/ @@ -39,7 +40,8 @@ #include #include -#include +#include +#include #include #include @@ -310,16 +312,29 @@ TEST_F(AnalyzerPoreFinderTest, ClusterAnalysisOnePore) for(int32_t y=0;y<16;y++) std::getline(file,linecontent); - - EXPECT_EQ(linecontent.substr(0,2),"16"); - EXPECT_EQ(linecontent.substr(3,2),"16"); - EXPECT_EQ(linecontent.substr(6,1),"1"); + + std::string buf; // Have a buffer string + std::stringstream ss(linecontent); // Insert the string into a stream + std::vector tokens; // Create vector to hold our words + + while (ss >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),16.0); + EXPECT_EQ(atof(tokens[1].c_str()),16.0); + EXPECT_EQ(atof(tokens[2].c_str()),1.0); std::getline(file,linecontent); - EXPECT_EQ(linecontent.substr(0,2),"16"); - EXPECT_EQ(linecontent.substr(3,2),"17"); - EXPECT_EQ(linecontent.substr(6,1),"1"); + std::stringstream ss2(linecontent); // Insert the string into a stream + tokens.clear(); + + while (ss2 >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),16.0); + EXPECT_EQ(atof(tokens[1].c_str()),17.0); + EXPECT_EQ(atof(tokens[2].c_str()),1.0); file.close(); remove("PoreCoordinates.dat"); @@ -397,10 +412,17 @@ TEST_F(AnalyzerPoreFinderTest, ClusterAnalysisTwoPoreandCopolymers) for(int32_t y=0;y<16;y++) std::getline(file,linecontent); - - EXPECT_EQ(linecontent.substr(0,2),"16"); - EXPECT_EQ(linecontent.substr(3,2),"16"); - EXPECT_EQ(linecontent.substr(6,1),"1"); + + std::string buf; // Have a buffer string + std::stringstream ss(linecontent); // Insert the string into a stream + std::vector tokens; // Create vector to hold our words + + while (ss >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),16.0); + EXPECT_EQ(atof(tokens[1].c_str()),16.0); + EXPECT_EQ(atof(tokens[2].c_str()),1.0); //If there are two pores, the analyzer //treats them like two pore events. @@ -408,9 +430,17 @@ TEST_F(AnalyzerPoreFinderTest, ClusterAnalysisTwoPoreandCopolymers) for(int32_t y=16;y<16+4;y++) std::getline(file,linecontent); - EXPECT_EQ(linecontent.substr(0,2),"16"); - EXPECT_EQ(linecontent.substr(3,2),"20"); - EXPECT_EQ(linecontent.substr(6,3),"0.5"); + ss.str(""); //clear stringstream + ss.clear(); + ss.str(linecontent); // Insert the string into a stream + tokens.clear(); + + while (ss >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),16.0); + EXPECT_EQ(atof(tokens[1].c_str()),20.0); + EXPECT_EQ(atof(tokens[2].c_str()),0.5); file.close(); @@ -421,10 +451,18 @@ TEST_F(AnalyzerPoreFinderTest, ClusterAnalysisTwoPoreandCopolymers) linecontent = passComments(file1); - EXPECT_EQ(linecontent[0],'0'); - EXPECT_EQ(linecontent[2],'0'); - EXPECT_EQ(linecontent[4],'0'); - + ss.str(""); //clear stringstream + ss.clear(); + ss.str(linecontent); // Insert the string into a stream + tokens.clear(); + + while (ss >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),0.0); + EXPECT_EQ(atof(tokens[1].c_str()),0.0); + EXPECT_EQ(atof(tokens[2].c_str()),0.0); + //For the first copolymer. for(int32_t x=0;x<16;x++) for(int32_t y=0;y<32;y++) @@ -432,26 +470,50 @@ TEST_F(AnalyzerPoreFinderTest, ClusterAnalysisTwoPoreandCopolymers) for(int32_t y=0;y<16-9;y++) std::getline(file1,linecontent); - - EXPECT_EQ(linecontent.substr(0,2),"16"); - EXPECT_EQ(linecontent.substr(3,1),"7"); - EXPECT_EQ(linecontent.substr(5,3),"0.5"); + + ss.str(""); //clear stringstream + ss.clear(); + ss.str(linecontent); // Insert the string into a stream + tokens.clear(); + + while (ss >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),16.0); + EXPECT_EQ(atof(tokens[1].c_str()),7.0); + EXPECT_EQ(atof(tokens[2].c_str()),0.5); //second copolymer for(int32_t y=16-9;y<16-3;y++) std::getline(file1,linecontent); - - EXPECT_EQ(linecontent.substr(0,2),"16"); - EXPECT_EQ(linecontent.substr(3,2),"13"); - EXPECT_EQ(linecontent.substr(6,3),"0.5"); + + ss.str(""); //clear stringstream + ss.clear(); + ss.str(linecontent); // Insert the string into a stream + tokens.clear(); + + while (ss >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),16.0); + EXPECT_EQ(atof(tokens[1].c_str()),13.0); + EXPECT_EQ(atof(tokens[2].c_str()),0.5); //third copolymer for(int32_t y=16-3;y<16+3;y++) std::getline(file1,linecontent); - EXPECT_EQ(linecontent.substr(0,2),"16"); - EXPECT_EQ(linecontent.substr(3,2),"19"); - EXPECT_EQ(linecontent.substr(6,3),"0.5"); + ss.str(""); //clear stringstream + ss.clear(); + ss.str(linecontent); // Insert the string into a stream + tokens.clear(); + + while (ss >> buf) + tokens.push_back(buf); + + EXPECT_EQ(atof(tokens[0].c_str()),16.0); + EXPECT_EQ(atof(tokens[1].c_str()),19.0); + EXPECT_EQ(atof(tokens[2].c_str()),0.5); remove("PoreCoordinates.dat"); remove("PolymerCoordinates.dat");