diff --git a/README.md b/README.md index f3dfca7..ea4bc23 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ This project needs the [LeMonADE-library](https://github.com/LeMonADE-project/Le This will run `cmake` and create the `build` directory. Now type in `make` which compiles the executables and runs the test if this option was given. -## Linear chains with cross linkers +## Linear chains with cross-links The main [program](https://github.com/LeMonADE-project/LeMonADE_PhantomModulus/blob/master/projects/ForceEquilibrium.cpp) reads in a .bfm file and a data file. It is expected that the bfm file contains the following flags ``` #!number_of_linear_chains=value @@ -23,12 +23,26 @@ The main [program](https://github.com/LeMonADE-project/LeMonADE_PhantomModulus/b 1-1:1/2 ... ``` -Moreover, it is expected that in the bfm file first all chains and afterwards all crosslinks are given. +Moreover, it is expected that in the bfm file first all chains and afterwards all cross-links are given. The data file should be created during the connection process and needs to contain: ``` Time ChainID MonID1 P1X P1Y P1Z MonID2 P2X P2Y P2Z ``` -where MonID1 and MonID2 are the cross link ID starting at 0. ChainID starts at 0 as well. +where MonID1 and MonID2 are the cross-link ID starting at 0. ChainID starts at 0 as well. For the calculations only the ChainID and MonID1 is needed. The remaining fields are arbitrary for a correct calculation of the positions of the cross links in a force equilibrium. (The format was given from a previous simulation.) + +## Write out parts of the networks +For the force equilibrium only active chains are of interest. +Therefore, there exist two programs extracting the active gel structure from the total system. + +## Tendomer +There are programs specialized on tendomers. They read the tendomer specific tag and can use a +custom force extension relation (e.g. numerical solution of the tendomer) to calculate the forces +between cross-links. + +## Ideal reference +These systems are tendomer specific. A star is created and a monomer along the chain with a distance to the center +according to the tendomer distribution of the number of elastic segments is marked as fixed in space. +In a second system a double star is simulated. \ No newline at end of file diff --git a/include/LeMonADE_PM/analyzer/AnalyzerEquilbratedPosition.h b/include/LeMonADE_PM/analyzer/AnalyzerEquilbratedPosition.h index ea38b60..339a013 100644 --- a/include/LeMonADE_PM/analyzer/AnalyzerEquilbratedPosition.h +++ b/include/LeMonADE_PM/analyzer/AnalyzerEquilbratedPosition.h @@ -75,7 +75,13 @@ template < class IngredientsType > class AnalyzerEquilbratedPosition : public Ab std::vector< std::vector > CalculateDistance(); //! just collects the id and the position for the cross links - std::vector > CollectAveragePositions(); + std::vector > CollectAveragePositions(); + + //! setter for the filenames + void setFilenames(std::string outAvPosBasename_, std::string outDistBasename_ ){ + outAvPosBasename= outAvPosBasename_; + outDistBasename=outDistBasename_; + } }; /************************************************************************* @@ -102,18 +108,18 @@ std::vector< std::vector > AnalyzerEquilbratedPosition auto neighbors(ingredients.getCrossLinkNeighborIDs(IDx)); for (size_t j=0; j < neighbors.size() ;j++){ VectorDouble3 vec(ingredients.getMolecules()[neighbors[j].ID].getVector3D()-ingredients.getMolecules()[IDx].getVector3D()-neighbors[j].jump); - VectorDouble3 vec2=LemonadeDistCalcs::MinImageVector(ingredients.getMolecules()[neighbors[j].ID].getVector3D(),ingredients.getMolecules()[IDx].getVector3D(),ingredients); - if ( vec.getLength() != vec2.getLength()) { - std::stringstream error_message; - error_message << "AnalyzerEquilbratedPosition:\n" - << "distance from jump=(" << vec << ") \n" - << "distance from IMC=(" << vec2 << ") \n" - << "jump vector=("<< neighbors[j].jump <<") \n" - << "position1=("< > AnalyzerEquilbratedPosition } //////////////////////////////////////////////////////////////////////////////// template< class IngredientsType > -std::vector< std::vector > AnalyzerEquilbratedPosition::CollectAveragePositions(){ - std::vector > AveragePosition(4, std::vector()); +std::vector< std::vector > AnalyzerEquilbratedPosition::CollectAveragePositions(){ + std::vector > AveragePosition(4, std::vector()); auto crosslinkID(ingredients.getCrosslinkIDs()); for (size_t i = 0 ; i < crosslinkID.size(); i++){ auto IDx(crosslinkID[i]); @@ -189,14 +195,14 @@ void AnalyzerEquilbratedPosition::dumpData() std::cout << "////////////////////////////////////"< > CrossLinkPositions=CollectAveragePositions() ; + std::vector< std::vector > CrossLinkPositions=CollectAveragePositions() ; std::stringstream commentAveragePosition; commentAveragePosition<<"Created by AnalyzerEquilbratedPosition\n"; commentAveragePosition<<"ID's start at 0 \n"; commentAveragePosition<<"conversion="<::dumpData() commentDistribution<<"Chain ID's start at 1 \n"; commentDistribution<<"ID1 ID2 vector length ChainID \n"; std::stringstream outDist; - outDist<< std::setprecision(2) << "C" << conversion; + outDist<< std::setprecision(3) << "C" << conversion; outDist << "_" << outDistBasename; ResultFormattingTools::writeResultFile( diff --git a/include/LeMonADE_PM/analyzer/AnalyzerMolecularWeight.h b/include/LeMonADE_PM/analyzer/AnalyzerMolecularWeight.h new file mode 100644 index 0000000..e13a979 --- /dev/null +++ b/include/LeMonADE_PM/analyzer/AnalyzerMolecularWeight.h @@ -0,0 +1,129 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ + +#ifndef LEMONADE_PM_ANALYZER_MOLECULAR_WEIGTH_H +#define LEMONADE_PM_ANALYZER_MOLECULAR_WEIGTH_H + +#include + +#include +#include +#include +#include +#include +#include + + +/************************************************************************* + * definition of AnalyzerMolecularWeight class + * ***********************************************************************/ + +/** + * @file + * + * @class AnalyzerMolecularWeight + * + * @tparam IngredientsType Ingredients class storing all system information( e.g. monomers, bonds, etc). + * + * @details Calculates the number and weight averaged molecular weight. + * + */ +template < class IngredientsType > +class AnalyzerMolecularWeight : public AnalyzerAbstractDump +{ + typedef AnalyzerAbstractDump BaseClass; +private: + //! typedef for the underlying container holding the monomers + typedef typename IngredientsType::molecules_type molecules_type; + + std::string basename; + +protected: + using BaseClass::ingredients; + using BaseClass::Data; + using BaseClass::MCSTimes; + using BaseClass::dumpTimeSeries; + +public: + //! constructor + AnalyzerMolecularWeight(const IngredientsType& ing, std::string fileSuffix_); + //! destructor. does nothing + virtual ~AnalyzerMolecularWeight(){} + //! calculate the moleculare weight distribution + virtual void initialize(); + virtual bool execute(); +}; + +/************************************************************************* + * implementation of memebers + * ***********************************************************************/ + +/** + * @param ing reference to the object holding all information of the system + * @param fileSuffix output file name. defaults to "Rg2TimeSeries.dat". + * */ +template +AnalyzerMolecularWeight::AnalyzerMolecularWeight( + const IngredientsType& ing,std::string fileSuffix_) +:BaseClass(ing,fileSuffix_) {} + +template +void AnalyzerMolecularWeight::initialize() +{ + basename=BaseClass::getOutputFilename() ; + BaseClass::setBufferSize(1); + BaseClass::setNumberOfColumns(1); + execute(); +} +template< class IngredientsType > +bool AnalyzerMolecularWeight::execute() +{ + double conversion = ingredients.getConversion(); + std::stringstream comment; + comment << "Created by AnalyzerMolecularWeight\n" + << "conversion=" << conversion << "%"; + BaseClass::setComment(comment.str()); + + //calculate the molecular weight distribution + std::vector > groups; + fill_connected_groups(ingredients.getMolecules(), groups, MonomerGroup(ingredients.getMolecules()), alwaysTrue()); + //key is the size of the molecule and the value is the number of occurence + std::map dist; + for (auto it=groups.begin(); it!=groups.end();it++ ) + dist[it->size()]++; + for ( auto& it : dist ) + { + Data[0].push_back(it.second); + MCSTimes.push_back(it.first); + } + std::stringstream output; + output << "MolecularWeight_c" << conversion << basename<< ".dat"; + BaseClass::resetIsFirstFileDump(); + BaseClass::setOutputFilename(output.str()); + dumpTimeSeries(); +} +#endif diff --git a/include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUp.h b/include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUp.h index 1f38489..2b2b3ff 100644 --- a/include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUp.h +++ b/include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUp.h @@ -71,14 +71,24 @@ class FeatureCrosslinkConnectionsLookUp : public Feature { void synchronize(IngredientsType& ingredients) { fillTables(ingredients); }; - + //! set the jump vector + void setCrossLinkNeighborJump(uint32_t CrossLinkID, uint32_t idx, VectorDouble3 vec) { + if ( idx > CrossLinkNeighbors.at(CrossLinkID).size() ){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUpTendomers::setCrossLinkNeighborJump neighbor idx " << idx <<" is to high."; + throw std::runtime_error(errormessage.str()); + } + CrossLinkNeighbors.at(CrossLinkID)[idx].jump=vec; + } //!getter function for the neighboring crosslinks std::vector getCrossLinkNeighborIDs(uint32_t CrossLinkID) const{ + #ifdef DEBUG if (CrossLinkNeighbors.find(CrossLinkID) == CrossLinkNeighbors.end()){ std::stringstream errormessage; errormessage << "FeatureCrosslinkConnectionsLookUp::getCrossLinkNeighborIDs Cross Link ID " << CrossLinkID <<" does not exist."; throw std::runtime_error(errormessage.str()); } + #endif return CrossLinkNeighbors.at(CrossLinkID); }; diff --git a/include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference.h b/include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference.h new file mode 100644 index 0000000..218c82e --- /dev/null +++ b/include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference.h @@ -0,0 +1,192 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ + + +#ifndef LENONADE_PM_FEATURE_FEATURECROSSLINKCONNECTIONLOOKUPIDEALREFERENCE_H +#define LENONADE_PM_FEATURE_FEATURECROSSLINKCONNECTIONLOOKUPIDEALREFERENCE_H + +#include +#include +#include +#include +#include +#include + + +/*****************************************************************************/ +/** + * @file + * @date 2021/04/01 + * @author Toni + * + * @class FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference + * @brief Creates a lookup for crosslink IDs and their neighbors + * @details For calculating the equilibrium position of crosslinks in a + * a phantom network, one needs to know the crosslink neighbors and the number + * of segments between them. + * + **/ +/*****************************************************************************/ +/////////////////////////////////////////////////////////////////////////////// +//DEFINITION OF THE CLASS TEMPLATE /////// +//Implementation of the members below /////// +/////////////////////////////////////////////////////////////////////////////// + +class FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference : public Feature { + +public: + //! This Feature requires a monomer_extensions. + typedef LOKI_TYPELIST_1(MonomerReactivity) monomer_extensions; + + //! check bas connect move - always true + template + bool checkMove(const IngredientsType& ingredients, const MoveBase& move) const { return true;}; + + //! synchronize lookup table + template + void synchronize(IngredientsType& ingredients) { + fillTables(ingredients); + }; + //! set the jump vector + void setCrossLinkNeighborJump(uint32_t CrossLinkID, uint32_t idx, VectorDouble3 vec) { + if ( idx > CrossLinkNeighbors.at(CrossLinkID).size() ){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUpTendomers::setCrossLinkNeighborJump neighbor idx " << idx <<" is to high."; + throw std::runtime_error(errormessage.str()); + } + CrossLinkNeighbors.at(CrossLinkID)[idx].jump=vec; + } + //!getter function for the neighboring crosslinks + std::vector getCrossLinkNeighborIDs(uint32_t CrossLinkID) const{ + #ifdef DEBUG + if (CrossLinkNeighbors.find(CrossLinkID) == CrossLinkNeighbors.end()){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference::getCrossLinkNeighborIDs Cross Link ID " << CrossLinkID <<" does not exist."; + throw std::runtime_error(errormessage.str()); + } + #endif + return CrossLinkNeighbors.at(CrossLinkID); + }; + + //!get the ID of crosslinks (determined by nConnections>3 and connected to another crosslink) + const std::vector& getCrosslinkIDs() const {return crosslinkIDs;} + +private: + //! convinience function to fill all tables + template + void fillTables(IngredientsType& ingredients); + //!key - CrossLink ID; value - neighboring cross link and the number of segments to them + std::map > CrossLinkNeighbors; + //!ID for crosslinks + std::vector crosslinkIDs; +}; +/** + *@details Create look up table + **/ +template +void FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference::fillTables(IngredientsType& ingredients){ + + const typename IngredientsType::molecules_type& molecules=ingredients.getMolecules(); + std::cout << "FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference::fillTables" < NeighborIDs; + // VectorDouble3 jumpVector(0.,0.,0.); + // uint32_t nSegments(ingredients.getNumOfMonomersPerChain()); + // for (auto i=0; i < ingredients.getMolecules().size(); i++){ + // if (ingredients.getMolecules()[i].getMovableTag() == false) + // NeighborIDs.push_back( neighborX(i, ( (i-1)%((2*nSegments+1)) )+1, jumpVector) ); + + // } + + // // for( auto i=0; i NeighborIDs1; + // // // NeighborIDs1 + // // } + // CrossLinkNeighbors[0]=NeighborIDs; + + + + for (uint32_t i = 0 ;i < molecules.size();i++){ + //find next crosslink + if( molecules.getNumLinks(i) > 2 ){ + //temporary storage for the neighboring crosslinks + std::vector NeighborIDs; + auto posX(molecules[i].getVector3D()); + for (size_t j = 0 ; j < molecules.getNumLinks(i); j++){ + uint32_t tail(i); + uint32_t head(molecules.getNeighborIdx(i,j)); + VectorDouble3 posHead(molecules[head].getVector3D()); + VectorDouble3 bond(LemonadeDistCalcs::MinImageVector( posX,posHead,ingredients)); + if(bond.getLength() > std::sqrt(10)){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUp: Wrong bond " << bond << " between " << i << " and " << head << "\n"; + throw std::runtime_error(errormessage.str()); + } + VectorDouble3 jumpVector(posHead-bond-posX); // tracks if one bond jumps across periodic images + + //direct connection of two cross links + if ( molecules.getNumLinks(head) > 2 || molecules[head].getMovableTag()==false ) { + NeighborIDs.push_back( neighborX(head, 1, jumpVector) ); + }else{ + uint32_t nSegments(1); + //cross links are connected by a chain + while( molecules.getNumLinks(head) == 2 && head != i ){ + //find next head + for (size_t k = 0 ; k < molecules.getNumLinks(head); k++){ + uint32_t NextMonomer( molecules.getNeighborIdx(head,k)); + if ( NextMonomer != tail ) { + tail=head; + head=NextMonomer; + break; + } + } + posHead=molecules[head].getVector3D(); + auto posTail=molecules[tail].getVector3D(); + bond=LemonadeDistCalcs::MinImageVector( posTail,posHead,ingredients); + jumpVector+=(posHead-bond-posTail); // tracks if one bond jumps across periodic images + nSegments++; + //a cross link has more than 2 connections + if( molecules.getNumLinks(head) > 2 || molecules[head].getMovableTag()==false ){ + NeighborIDs.push_back(neighborX(head, nSegments, jumpVector)); + break; + } + } + } + } + CrossLinkNeighbors[i]=NeighborIDs; + // if(NeighborIDs.size()>0) + crosslinkIDs.push_back(i); + } + } + + + std::cout << "FeatureCrosslinkConnectionsLookUpIdealDoubleStarReference::fillTables.done" <. + +--------------------------------------------------------------------------------*/ + + +#ifndef LENONADE_PM_FEATURE_FEATURECROSSLINKCONNECTIONLOOKUPIDEALREFERENCE_H +#define LENONADE_PM_FEATURE_FEATURECROSSLINKCONNECTIONLOOKUPIDEALREFERENCE_H + +#include +#include +#include +#include +#include +#include + + +/*****************************************************************************/ +/** + * @file + * @date 2021/04/01 + * @author Toni + * + * @class FeatureCrosslinkConnectionsLookUpIdealReference + * @brief Creates a lookup for crosslink IDs and their neighbors + * @details For calculating the equilibrium position of crosslinks in a + * a phantom network, one needs to know the crosslink neighbors and the number + * of segments between them. + * + **/ +/*****************************************************************************/ +/////////////////////////////////////////////////////////////////////////////// +//DEFINITION OF THE CLASS TEMPLATE /////// +//Implementation of the members below /////// +/////////////////////////////////////////////////////////////////////////////// + +class FeatureCrosslinkConnectionsLookUpIdealReference : public Feature { + +public: + //! This Feature requires a monomer_extensions. + typedef LOKI_TYPELIST_1(MonomerReactivity) monomer_extensions; + + //! check bas connect move - always true + template + bool checkMove(const IngredientsType& ingredients, const MoveBase& move) const { return true;}; + + //! synchronize lookup table + template + void synchronize(IngredientsType& ingredients) { + fillTables(ingredients); + }; + //! set the jump vector + void setCrossLinkNeighborJump(uint32_t CrossLinkID, uint32_t idx, VectorDouble3 vec) { + if ( idx > CrossLinkNeighbors.at(CrossLinkID).size() ){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUpTendomers::setCrossLinkNeighborJump neighbor idx " << idx <<" is to high."; + throw std::runtime_error(errormessage.str()); + } + CrossLinkNeighbors.at(CrossLinkID)[idx].jump=vec; + } + //!getter function for the neighboring crosslinks + std::vector getCrossLinkNeighborIDs(uint32_t CrossLinkID) const{ + #ifdef DEBUG + if (CrossLinkNeighbors.find(CrossLinkID) == CrossLinkNeighbors.end()){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUpIdealReference::getCrossLinkNeighborIDs Cross Link ID " << CrossLinkID <<" does not exist."; + throw std::runtime_error(errormessage.str()); + } + #endif + return CrossLinkNeighbors.at(CrossLinkID); + }; + + //!get the ID of crosslinks (determined by nConnections>3 and connected to another crosslink) + const std::vector& getCrosslinkIDs() const {return crosslinkIDs;} + +private: + //! convinience function to fill all tables + template + void fillTables(IngredientsType& ingredients); + //!key - CrossLink ID; value - neighboring cross link and the number of segments to them + std::map > CrossLinkNeighbors; + //!ID for crosslinks + std::vector crosslinkIDs; +}; +/** + *@details Create look up table + **/ +template +void FeatureCrosslinkConnectionsLookUpIdealReference::fillTables(IngredientsType& ingredients){ + + const typename IngredientsType::molecules_type& molecules=ingredients.getMolecules(); + std::cout << "FeatureCrosslinkConnectionsLookUpIdealReference::fillTables" < NeighborIDs; + VectorDouble3 jumpVector(0.,0.,0.); + uint32_t nSegments(ingredients.getNumOfMonomersPerChain()); + for (auto i=0; i < ingredients.getMolecules().size(); i++){ + if (ingredients.getMolecules()[i].getMovableTag() == false) + NeighborIDs.push_back( neighborX(i, ( (i-1)%((2*nSegments+1)) )+1, jumpVector) ); + + } + + // for( auto i=0; i NeighborIDs1; + // // NeighborIDs1 + // } + CrossLinkNeighbors[0]=NeighborIDs; + + std::cout << "FeatureCrosslinkConnectionsLookUpIdealReference::fillTables.done" <. + +--------------------------------------------------------------------------------*/ + + +#ifndef LENONADE_PM_FEATURE_FEATURECROSSLINKCONNECTIONLOOKUPTENDOMERS_H +#define LENONADE_PM_FEATURE_FEATURECROSSLINKCONNECTIONLOOKUPTENDOMERS_H + +#include +#include +#include +#include +#include +#include + + +/*****************************************************************************/ +/** + * @file + * @date 2021/04/01 + * @author Toni + * + * @class FeatureCrosslinkConnectionsLookUpTendomers + * @brief Creates a lookup for crosslink IDs and their neighbors + * @details For calculating the equilibrium position of crosslinks in a + * a phantom network, one needs to know the crosslink neighbors and the number + * of segments between them. + * + **/ +/*****************************************************************************/ +/////////////////////////////////////////////////////////////////////////////// +//DEFINITION OF THE CLASS TEMPLATE /////// +//Implementation of the members below /////// +/////////////////////////////////////////////////////////////////////////////// + +class FeatureCrosslinkConnectionsLookUpTendomers : public Feature { + +public: + //! This Feature requires a monomer_extensions. + typedef LOKI_TYPELIST_1(MonomerReactivity) monomer_extensions; + + //! check bas connect move - always true + template + bool checkMove(const IngredientsType& ingredients, const MoveBase& move) const { return true;}; + + //! synchronize lookup table + template + void synchronize(IngredientsType& ingredients) { + fillTables(ingredients); + }; + + //!getter function for the neighboring crosslinks + std::vector getCrossLinkNeighborIDs(uint32_t CrossLinkID) const{ + if (CrossLinkNeighbors.find(CrossLinkID) == CrossLinkNeighbors.end()){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUpTendomers::getCrossLinkNeighborIDs Cross Link ID " << CrossLinkID <<" does not exist."; + throw std::runtime_error(errormessage.str()); + } + return CrossLinkNeighbors.at(CrossLinkID); + }; + + //! set the jump vector + void setCrossLinkNeighborJump(uint32_t CrossLinkID, uint32_t idx, VectorDouble3 vec) { + if ( idx > CrossLinkNeighbors.at(CrossLinkID).size() ){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUpTendomers::setCrossLinkNeighborJump neighbor idx " << idx <<" is to high."; + throw std::runtime_error(errormessage.str()); + } + CrossLinkNeighbors.at(CrossLinkID)[idx].jump=vec; + } + + //!get the ID of crosslinks (determined by nConnections>3 and connected to another crosslink) + const std::vector& getCrosslinkIDs() const {return crosslinkIDs;} + +private: + //! convinience function to fill all tables + template + void fillTables(IngredientsType& ingredients); + //!key - CrossLink ID; value - neighboring cross link and the number of segments to them + std::map > CrossLinkNeighbors; + //!ID for crosslinks + std::vector crosslinkIDs; +}; +/** + *@details Create look up table + **/ +template +void FeatureCrosslinkConnectionsLookUpTendomers::fillTables(IngredientsType& ingredients){ + + const typename IngredientsType::molecules_type& molecules=ingredients.getMolecules(); + std::cout << "FeatureCrosslinkConnectionsLookUpTendomers::fillTables" < 2 ){ + //temporary storage for the neighboring crosslinks + std::vector NeighborIDs; + auto posX(molecules[i].getVector3D()); + for (size_t j = 0 ; j < molecules.getNumLinks(i); j++){ + uint32_t tail(i); + uint32_t head(molecules.getNeighborIdx(i,j)); + VectorDouble3 posHead(molecules[head].getVector3D()); + VectorDouble3 bond(LemonadeDistCalcs::MinImageVector( posX,posHead,ingredients)); + if(bond.getLength() > std::sqrt(10)){ + std::stringstream errormessage; + errormessage << "FeatureCrosslinkConnectionsLookUpTendomers: Wrong bond " << bond << " between " << i << " and " << head << "\n"; + throw std::runtime_error(errormessage.str()); + } + VectorDouble3 jumpVector(posHead-bond-posX); // tracks if one bond jumps across periodic images + //go to tendomer end + tail=(head/nMonomersPerChain % 2 == 0 ) ? head +nMonomersPerChain : head-nMonomersPerChain; + //cross links are connected by a chain + if( molecules.getNumLinks(tail) == 2 ){ + //get crosslink + for (size_t k = 0 ; k < molecules.getNumLinks(tail); k++){ + uint32_t NextMonomer( molecules.getNeighborIdx(tail,k)); + if( molecules[NextMonomer].isReactive() && molecules[NextMonomer].getNumMaxLinks() > 2 ){ + head=NextMonomer; + break; + } + } + posHead=molecules[head].getVector3D(); + auto posTail=molecules[tail].getVector3D(); + bond=LemonadeDistCalcs::MinImageVector( posTail,posHead,ingredients); + jumpVector+=(posHead-bond-posTail); // tracks if one bond jumps across periodic images + // auto vecJ(jumpVector); + // vecJ.setAllCoordinates( + // static_cast(jumpVector.getX()) % 256 , + // static_cast(jumpVector.getY()) % 256 , + // static_cast(jumpVector.getZ()) % 256 ); + + // std::cout << "Jumpvector=" << jumpVector << " for ID=" << i << " to " << head << ": " <. + +--------------------------------------------------------------------------------*/ +#ifndef LEMONADEPM_UPDATER_UPDATERABSTRACTCREATEALLBONDS_H +#define LEMONADEPM_UPDATER_UPDATERABSTRACTCREATEALLBONDS_H + + +/** + * @file + * + * @class UpdaterAbstractCreateAllBonds + * + * @brief abstract updater to create systems + * + * @details This abstract class provides the three basic functions to create systems: add a single monomer, add a connected monomer and move the system to find some free space. + * This Updater requires FeatureAttributes. + * + * @tparam IngredientsType + * + **/ + +#include +#include +#include +#include +#include +#include + +template +class UpdaterAbstractCreateAllBonds:public AbstractUpdater +{ +public: + UpdaterAbstractCreateAllBonds(IngredientsType& ingredients_):ingredients(ingredients_),lookupFilled(false){} + + virtual void initialize(); + virtual bool execute(); + virtual void cleanup(); + +protected: + IngredientsType& ingredients; + + //! function to add a standalone monomer + bool addSingleMonomer(int32_t type=1); + + //! function to add a monomer to a parent monomer + bool addMonomerToParent(uint32_t parent_id, int32_t type=1); + + //! function to add a monomer at a spezific place + bool addMonomerAtPosition(VectorInt3 pos, int32_t type=1); + + //! function to add a new monomer inbetween two others + bool addMonomerInsideConnectedPair(uint32_t indexA, uint32_t indexB, int32_t type=1); + + //! function to add a new monomer at two others + bool addMonomerAtConnectedPair(uint32_t indexA, uint32_t indexB, int32_t type=1); + + //! function to add a ring around a chain monomer + bool addRing(uint32_t parent_id, int32_t type=1, uint32_t NRingMonomers=5); + + //! function to move the whole system + void moveSystem(int32_t nsteps); + + //! function to find groups of connected monomers + void linearizeSystem(); + + //! function to get a random bondvector of length 2 + VectorInt3 randomBondvector(); + +private: + RandomNumberGenerators rng; + + //! bond table is filled with all valid bonds or not + bool lookupFilled; + //!set of vectors which are valid for the conformations + std::vector bondvectorSet; + +}; + +/** +* The initialize function handles the new systems information. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template < class IngredientsType > +void UpdaterAbstractCreateAllBonds::initialize(){ + +} + +/** +* Execution of the system creation +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template < class IngredientsType > +bool UpdaterAbstractCreateAllBonds::execute(){ + +} + +/** +* Standard clean up. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template < class IngredientsType > +void UpdaterAbstractCreateAllBonds::cleanup(){ + +} + +/******************************************************************************/ +/** + * @brief function to add a standalone monomer + * @param type attribute tag of the new monomer + * @return if position is not free, if move was applied + */ +template +bool UpdaterAbstractCreateAllBonds::addSingleMonomer(int32_t type){ + // set properties of add Monomer Move + MoveAddMonomerSc<> addmove; + addmove.init(ingredients); + addmove.setTag(type); + + int32_t counter(0); + while(counter<10000){ + VectorInt3 newPosition((rng.r250_rand32() % (ingredients.getBoxX()-1)), + (rng.r250_rand32() % (ingredients.getBoxY()-1)), + (rng.r250_rand32() % (ingredients.getBoxZ()-1))); + addmove.setPosition(newPosition); + if(addmove.check(ingredients)==true){ + addmove.apply(ingredients); + return true; + } + counter++; + } + return false; +} + +/******************************************************************************/ +/** + * @brief function to add a monomer to a parent monomer + * @param parent_id id of monomer to connect with + * @param type attribute tag of the new monomer + * @return if position is not free, if move was applied + */ +template +bool UpdaterAbstractCreateAllBonds::addMonomerToParent(uint32_t parent_id, int32_t type){ + // set properties of add Monomer Move + MoveAddMonomerSc<> addmove; + addmove.init(ingredients); + addmove.setTag(type); + + int32_t counter(0); + + while(counter<10000){ + //try at most 30 random bondvectors to find a new monomer position + for(uint i=0;i<30;i++){ + VectorInt3 bondvector(randomBondvector()); + // set position of new monomer + addmove.setPosition(ingredients.getMolecules()[parent_id]+bondvector); + + // check new position (excluded volume) + if(addmove.check(ingredients)==true){ + addmove.apply(ingredients); + ingredients.modifyMolecules().connect( parent_id, (ingredients.getMolecules().size()-1) ); + return true; + } + } + // if no position matches, we need to move the system a bit + moveSystem(2); + counter++; + } + return false; +} + +/******************************************************************************/ +/** + * @brief function to add a monomer to a specific position. if position is not free return false + * @param VectorInt3 position + * @param type attribute tag of the new monomer + * @return if position is not free, if move was applied + */ +template +bool UpdaterAbstractCreateAllBonds::addMonomerAtPosition(VectorInt3 position, int32_t type){ + MoveAddMonomerSc<> addmove; + addmove.init(ingredients); + addmove.setTag(type); + + addmove.setPosition(position); + if(addmove.check(ingredients)==true){ + addmove.apply(ingredients); + return true; + }else{ + return false; + } +} + +/******************************************************************************/ +/** + * @brief function to add a new monomer between two already existing ones instead of a bond. + * @param indexA + * @param indexB + * @param type attribute tag of the new monomer + * @return if position is not free, if move was applied + */ +template +bool UpdaterAbstractCreateAllBonds::addMonomerInsideConnectedPair(uint32_t indexA, uint32_t indexB, int32_t type){ + //first check if monomers are connected + if( ! ingredients.getMolecules().areConnected(indexA,indexB)) + return false; + + MoveAddMonomerSc<> addmove; + addmove.init(ingredients); + addmove.setTag(type); + + int32_t counter(0); + + while(counter<10000){ + //try at most 30 random bondvectors to find a new monomer position + for(uint i=0;i<30;i++){ + // get a random bondvector + VectorInt3 bondvector(randomBondvector()); + // set position of new monomer + addmove.setPosition(ingredients.getMolecules()[indexA]+bondvector); + + // check new position (excluded volume, other features) + if(addmove.check(ingredients)==true){ + //check the new bondvector bewten the new monomer and indexB + VectorInt3 checkBV(addmove.getPosition()-ingredients.getMolecules()[indexB]); + if( (checkBV.getLength() < 3) && (ingredients.getBondset().isValidStrongCheck(checkBV)) ){ + addmove.apply(ingredients); + ingredients.modifyMolecules().connect( indexA, (ingredients.getMolecules().size()-1) ); + ingredients.modifyMolecules().connect( indexB, (ingredients.getMolecules().size()-1) ); + ingredients.modifyMolecules().disconnect( indexA, indexB ); + return true; + } + } + } + // if no position matches, we need to move the system a bit + moveSystem(2); + counter++; + } + return false; +} +/******************************************************************************/ +/** + * @brief function to add a new monomer to two already existing ones such that all three monomers are connected. + * @param indexA + * @param indexB + * @param type attribute tag of the new monomer + * @return if position is not free, if move was applied + */ +template +bool UpdaterAbstractCreateAllBonds::addMonomerAtConnectedPair(uint32_t indexA, uint32_t indexB, int32_t type){ + //first check if monomers are connected + if( ! ingredients.getMolecules().areConnected(indexA,indexB)) + return false; + + MoveAddMonomerSc<> addmove; + addmove.init(ingredients); + addmove.setTag(type); + + int32_t counter(0); + + while(counter<10000){ + //try at most 30 random bondvectors to find a new monomer position + for(uint i=0;i<30;i++){ + // get a random bondvector + VectorInt3 bondvector(randomBondvector()); + // set position of new monomer + addmove.setPosition(ingredients.getMolecules()[indexA]+bondvector); + + // check new position (excluded volume, other features) + if(addmove.check(ingredients)==true){ + //check the new bondvector between the new monomer and indexB + VectorInt3 checkBV(addmove.getPosition()-ingredients.getMolecules()[indexB]); + if( (checkBV.getLength() < 3) && (ingredients.getBondset().isValidStrongCheck(checkBV)) ){ + addmove.apply(ingredients); + ingredients.modifyMolecules().connect( indexA, (ingredients.getMolecules().size()-1) ); + ingredients.modifyMolecules().connect( indexB, (ingredients.getMolecules().size()-1) ); + return true; + } + } + } + // if no position matches, we need to move the system a bit + moveSystem(2); + counter++; + } + return false; +} +/** + * @brief add a ring threaded on a chain at position of parent monomer + * + * @todo are there conformations where bonds of the ring cross the bonds of the chain? + * + * @param parent + * @param type monomer tag + * @param NRingMonomers number of monomers per ring + */ +template < class IngredientsType > +bool UpdaterAbstractCreateAllBonds::addRing(uint32_t parent, int32_t type, uint32_t NRingMonomers) + { + if(parent < ingredients.getMolecules().size()) + { + uint32_t attempts(0); + while (attempts<10000) + { + uint32_t neighborDirection(0); + uint32_t currentParent(parent); + while (true) + { + uint32_t neighborID(ingredients.getMolecules().getNeighborIdx(currentParent,neighborDirection)); + if(ingredients.getMolecules().getNumLinks(neighborID) != 2 ) + { + if(neighborDirection==1) break; + else {neighborDirection++; currentParent=parent;} + } + else + { + VectorInt3 neighborPosition(ingredients.getMolecules()[neighborID]); + VectorInt3 parentPosition(ingredients.getMolecules()[currentParent]); + VectorInt3 bond(neighborPosition-parentPosition); + uint32_t otherParentNeighbor; + uint32_t otherNeighborNeighbor; + neighborDirection ? otherParentNeighbor=(ingredients.getMolecules().getNeighborIdx(currentParent,0)) : otherParentNeighbor=(ingredients.getMolecules().getNeighborIdx(currentParent,1)); + (ingredients.getMolecules().getNeighborIdx(neighborID,0)==currentParent) ? otherNeighborNeighbor=(ingredients.getMolecules().getNeighborIdx(neighborID,1)): otherNeighborNeighbor=(ingredients.getMolecules().getNeighborIdx(neighborID,0)); + VectorInt3 neighborParentBond(parentPosition-ingredients.getMolecules()[otherParentNeighbor]); + VectorInt3 neighborNeighborBond(neighborPosition-ingredients.getMolecules()[otherNeighborNeighbor]); + if(bond.getLength() == 2 && (neighborParentBond.getLength()<3) &&(neighborNeighborBond.getLength()<3)) + { + for(uint32_t i=0; i<8; i++) + { + + VectorInt3 StartPosition=parentPosition; + //set the possible conformations + int32_t dx1, dx2; + int32_t dy1, dz1; + int32_t dy2, dz2; + VectorInt3 vec1,vec2; + if(bond.getX() == 2 || bond.getX() == -2) + { + StartPosition+=VectorInt3(bond.getX()/2,0,0); + dx1=0;dx2=0;//dy1=2;dy2=0;dz1=0;dz2=2; + if(rng.r250_drand()>0.5){dy1=2;dy2=0;dz1=0;dz2=2;} + else {dy1=0;dy2=2;dz1=2;dz2=0;} + int32_t dx3,dx4; + if(rng.r250_drand()>0.5){dx3=1;dx4=-1;} + else{dx3=-1;dx4=1;} +// rng.r250_drand()>0.5 ? dx3=1 : dx3=-1; +// rng.r250_drand()>0.5 ? dx4=1 : dx4=-1; + vec1=VectorInt3( dx2+dx3, dy2, dz2); + vec2=VectorInt3(-dx2+dx4, -dy2, -dz2); + } + else if(bond.getY() == 2 || bond.getY() == -2) + { + StartPosition+=VectorInt3(0,bond.getY()/2,0); + dy1=0;dy2=0; + if(rng.r250_drand()>0.5){dx1=2;dx2=0;dz1=0;dz2=2;} + else {dx1=0;dx2=2;dz1=2;dz2=0;} + int32_t dy3,dy4; + if(rng.r250_drand()>0.5){dy3=1;dy4=-1;} + else{dy3=-1;dy4=1;} +// rng.r250_drand()>0.5 ? dy3=1 : dy3=-1; +// rng.r250_drand()>0.5 ? dy4=1 : dy4=-1; + vec1=VectorInt3( dx2, dy2+dy3, dz2); + vec2=VectorInt3(-dx2, -dy2+dy4, -dz2); + } + else if(bond.getZ() == 2 || bond.getZ() == -2) + { + StartPosition+=VectorInt3(0,0,bond.getZ()/2); + dz1=0;dz2=0; + if(rng.r250_drand()>0.5){dx1=2;dx2=0;dy1=0;dy2=2;} + else {dx1=0;dx2=2;dy1=2;dy2=0;} + int32_t dz3,dz4; + if(rng.r250_drand()>0.5){dz3=1;dz4=-1;} + else{dz3=-1;dz4=1;} +// rng.r250_drand()>0.5 ? dz3=1 : dz3=-1; +// rng.r250_drand()>0.5 ? dz4=1 : dz4=-1; + vec1=VectorInt3( dx2, dy2, dz2+dz3); + vec2=VectorInt3(-dx2, -dy2, -dz2+dz4); + } + // set positions of monomers and two which guarantee that ring is threaded + std::vector PotentialPositions(6,StartPosition); + PotentialPositions[0]+=VectorInt3( dx1, dy1, dz1); + PotentialPositions[2]+=VectorInt3(-dx1,-dy1,-dz1); + PotentialPositions[4]+=VectorInt3( dx2, dy2, dz2);//position to check, not to use for a monomer + PotentialPositions[5]+=VectorInt3(-dx2,-dy2,-dz2);//position to check, not to use for a monomer + PotentialPositions[1]+=vec1; + PotentialPositions[3]+=vec2; + + // check if positions are occupied + bool PositionsFit(true); + for(uint32_t i=0; i addmove; + addmove.init(ingredients); + addmove.setPosition(PotentialPositions[i]); + if(addmove.check(ingredients)==false){PositionsFit=false;} + } + // add ring to system + if (PositionsFit) + { + bool AddRing(true); + for(uint32_t i=0; i3 && AddRing ) + { + for (uint32_t j=0;j +void UpdaterAbstractCreateAllBonds::moveSystem(int32_t nsteps){ + MoveLocalSc move; + for(int32_t n=0;n +void UpdaterAbstractCreateAllBonds::linearizeSystem(){ + //call ingredients copy constructor + IngredientsType newIngredients(ingredients); + + //delete all the informations in molecules + newIngredients.modifyMolecules().clear(); + + std::vector < MonomerGroup > LinearMonomerGroupsVector; + + fill_connected_groups( ingredients.getMolecules(), LinearMonomerGroupsVector, MonomerGroup(ingredients.getMolecules()), alwaysTrue() ); + + for(size_t groups=0; groups < LinearMonomerGroupsVector.size(); ++groups){ + if(groups==0) + newIngredients.modifyMolecules() = LinearMonomerGroupsVector[groups].copyGroup(); + else + newIngredients.modifyMolecules() += LinearMonomerGroupsVector[groups].copyGroup(); + } + + ingredients=newIngredients; + +} + +/******************************************************************************/ +/** + * @brief function to get a random bondvector of length 2 which is part of the bondvectorset + * @param nsteps number of MCS to move + */ +template +VectorInt3 UpdaterAbstractCreateAllBonds::randomBondvector(){ + if ( lookupFilled == false ){ + + for (auto it=ingredients.getBondset().begin(); it!= ingredients.getBondset().end();it++){ + auto bondvector(it->second); + if (ingredients.getBondset().isValidStrongCheck( bondvector ) ) + bondvectorSet.push_back(bondvector); + } + lookupFilled=true; + } + //get a random direction for the bondvector of length 2 + return bondvectorSet[ rng.r250_rand32() % bondvectorSet.size()] ; +} +#endif /* LEMONADE_UPDATER_ABSTRACT_CREATE_H */ diff --git a/include/LeMonADE_PM/updater/UpdaterAddStars.h b/include/LeMonADE_PM/updater/UpdaterAddStars.h new file mode 100644 index 0000000..d0c8ace --- /dev/null +++ b/include/LeMonADE_PM/updater/UpdaterAddStars.h @@ -0,0 +1,196 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ + +#ifndef LEMONADE_UPDATER_SETUP_STARSAB +#define LEMONADE_UPDATER_SETUP_STARSAB +/** + * @file + * + * @class UpdaterAddStars + * + * @brief Updater to create a solution of monodisperse branched stars. + * + * @details This is a simple implementation of a system setup starting from an empty ingredients + * or a system with some monomers inside. This updater requires FeatureAttributes. + * Two tags are added to the monomers in alternating manner, usually needed for GPU computing. // not realized in this first attempt + * // RS, 05 Nov 2019, but presumably 29 July 2020 + * + * @tparam IngredientsType + * + * @param ingredients_ The system, holding either an empty simulation box for system setup + * or a prefilled ingredients where the stars shall be added + * @param NStar_ number of stars that are added to ingredients + * @param NMonoPerStar_ number of monomers in each star + * @param NMonoPerBranch_ number of monomers in each branch (excluding central monomer) + * @param NBranchPerStar_ number of branches in each star + * @param type1_ attribute tag of "even" monomers + * @param type2_ attribute tag of "odd" monomers + **/ + +// #include +#include +#include + +#include + +template +class UpdaterAddStars: public UpdaterAbstractCreateAllBonds +{ + typedef UpdaterAbstractCreateAllBonds BaseClass; + +public: + UpdaterAddStars(IngredientsType& ingredients_, uint32_t NStar_, + uint32_t NMonoPerBranch_, uint32_t NBranchPerStar_); + + virtual void initialize(); + virtual bool execute(); + virtual void cleanup(); + + //! getter function for number of stars + const int32_t getNStar() const {return NStar;} + + //! getter function for number of monomers in branch + const int32_t getNMonoPerBranch() const {return NMonoPerBranch;} + + //! getter function for number of branches in stars + const int32_t getNBranchPerStar() const {return NBranchPerStar;} + + //! getter function for calculated density + const double getDensity() const {return density;} + +private: + // provide access to functions of UpdaterAbstractCreate used in this updater + using BaseClass::ingredients; + using BaseClass::addMonomerToParent; + using BaseClass::addSingleMonomer; + using BaseClass::linearizeSystem; + + //! number of stars in the box + uint32_t NStar; + + //! number of monomers in a branch + uint32_t NMonoPerBranch; + + //! number of branches in a star + uint32_t NBranchPerStar; + + //! lattice occupation density + double density; + + //! bool for execution + bool wasExecuted; +}; + +/** +* @brief Constructor handling the new systems paramters +* +* @param ingredients_ a reference to the IngredientsType - mainly the system +* @param NStar_ number of stars to be added in the system instead of solvent +* @param NBranchPerStar_ number of branches in each star +* @param NMonoPerBranch_ number of monomers in each branch (excluding central monomer) +*/ + + +template < class IngredientsType > +UpdaterAddStars::UpdaterAddStars( + IngredientsType& ingredients_, + uint32_t NStar_, + uint32_t NMonoPerBranch_, + uint32_t NBranchPerStar_ + ): + BaseClass(ingredients_), NStar(NStar_), NMonoPerBranch(NMonoPerBranch_), NBranchPerStar(NBranchPerStar_), + density(0.0), wasExecuted(false) + {} + +/** +* @brief initialise function, calculate the target density to compare with at the end. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template < class IngredientsType > +void UpdaterAddStars::initialize(){ + std::cout << "initialize UpdaterAddStars" << std::endl; + + // get the target density from the sum of existing monomers and the new added chains + density=(double)( ingredients.getMolecules().size() + (NMonoPerBranch*NBranchPerStar+1) *NStar ) * 8 /(double)( ingredients.getBoxX()*ingredients.getBoxY()*ingredients.getBoxZ() ); + + std::cout << "add "< +bool UpdaterAddStars::execute(){ + if(wasExecuted) return true; + std::cout << "execute UpdaterAddStars" << std::endl; + + for(uint32_t i=0;i 0.0000000001 ){ + std::cout << density << " " <<( (ingredients.getMolecules().size()*8) / lattice_volume)< +#include +#include + +#include + +template +class UpdaterAddTMDoubleStars: public UpdaterAbstractCreateAllBonds +{ + typedef UpdaterAbstractCreateAllBonds BaseClass; + +public: + UpdaterAddTMDoubleStars(IngredientsType& ingredients_, uint32_t NStar_, + uint32_t NMonoPerBranch_, uint32_t nRings_, uint32_t NBranchPerStar_); + + virtual void initialize(); + virtual bool execute(); + virtual void cleanup(); + + //! getter function for number of stars + const int32_t getNStar() const {return NStar;} + + //! getter function for number of monomers in branch + const int32_t getNMonoPerBranch() const {return NMonoPerBranch;} + + //! getter function for number of branches in stars + const int32_t getNBranchPerStar() const {return NBranchPerStar;} + + //! getter function for calculated density + const double getDensity() const {return density;} + +private: + // provide access to functions of UpdaterAbstractCreate used in this updater + using BaseClass::ingredients; + using BaseClass::addMonomerToParent; + using BaseClass::addSingleMonomer; + using BaseClass::linearizeSystem; + double prob_q(double n, double N, double m) { + return (m+1.)/(N-n-1.); + } + //! number of rings + uint32_t nRings; + //! number of stars in the box + uint32_t NStar; + + //! number of monomers in a branch + uint32_t NMonoPerBranch; + + //! number of branches in a star + uint32_t NBranchPerStar; + + //! lattice occupation density + double density; + + //! bool for execution + bool wasExecuted; + + // adds a chain of with nMonomers monomers to the parentID + void createChain(uint parentID, uint32_t nMonomers); + // add a chain to the system at a random positions + void createChain(uint32_t nMonomers); + + //! + const uint32_t steps; + //! + std::vector invCPF; + // random number generator (globally seeded) + RandomNumberGenerators rng; +}; + +/** +* @brief Constructor handling the new systems paramters +* +* @param ingredients_ a reference to the IngredientsType - mainly the system +* @param NStar_ number of double stars to be added in the system instead of solvent +* @param NBranchPerStar_ number of branches in each star +* @param NMonoPerBranch_ number of monomers in each branch (excluding central monomer) +*/ + + +template < class IngredientsType > +UpdaterAddTMDoubleStars::UpdaterAddTMDoubleStars( + IngredientsType& ingredients_, + uint32_t NStar_, + uint32_t NMonoPerBranch_, + uint32_t nRings_, + uint32_t NBranchPerStar_ + ): + BaseClass(ingredients_), NStar(NStar_), NMonoPerBranch(NMonoPerBranch_), nRings(nRings_), NBranchPerStar(NBranchPerStar_), + density(0.0), wasExecuted(false),steps(50000) + { + invCPF.resize(steps-1,0); + } + +/** +* @brief initialise function, calculate the target density to compare with at the end. +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template < class IngredientsType > +void UpdaterAddTMDoubleStars::initialize(){ + std::cout << "initialize UpdaterAddTMDoubleStars" << std::endl; + + + //probability distribution function + std::cout << "Calculate PDF" << std::endl; + std::vector PF((NMonoPerBranch-nRings),0); + double sum(0.); + for (auto i = 1; i < (NMonoPerBranch- nRings); i++ ){ + PF[i]=prob_q(i,NMonoPerBranch,nRings)*(1.-sum); + // std::cout << "probabilityDisti "< convPF((NMonoPerBranch- nRings)*2,0); + for (auto i=0; i < 2*(NMonoPerBranch-nRings);i++){ + for(auto j=0; j < i ; j++){ + convPF[i]+=PF[j]*PF[i-j]; + } + // std::cout << "ConvprobabilityDisti "< CPF((NMonoPerBranch- nRings)*2,0); + for (auto i=1; i < 2*(NMonoPerBranch-nRings);i++){ + CPF[i]+=CPF[i-1]+convPF[i]; + } + //calculate the inverse cummulative convoluted probability distribution + std::cout << "Calculate inverse cummulative convoluted PDF" << std::endl; + for (auto i =1 ; i < steps-1; i++) { + uint32_t n=0 ; + auto prob=static_cast (i) / static_cast(steps ); + while ( prob > CPF [ n ] ){ n++;} + invCPF[i]=(n-1) + static_cast(round( (n-(n-1)) *( prob- CPF[n-1] )/(CPF[n]-CPF[n-1]) )); + // std::cout << "invCPF: " << prob << " " << invCPF[i] <<" "<< n-1 <<" "< +void UpdaterAddTMDoubleStars::createChain(uint32_t parentID, uint32_t nMonomers){ + auto ID(parentID); + for (auto i=0; i< nMonomers ; i++){ + addMonomerToParent(parentID,1); + parentID=(ingredients.getMolecules().size()-1); + } +} + +template < class IngredientsType > +void UpdaterAddTMDoubleStars::createChain(uint32_t nMonomers){ + addSingleMonomer(1); + createChain(ingredients.getMolecules().size()-1,nMonomers-1); +} +/** +* @brief Execution of the system creation +* +* @tparam IngredientsType Features used in the system. See Ingredients. +*/ +template < class IngredientsType > +bool UpdaterAddTMDoubleStars::execute(){ + if(wasExecuted) return true; + std::cout << "execute UpdaterAddTMDoubleStars" << std::endl; + //initial number of monomer in the system + auto nMonomers(ingredients.getMolecules().size()); + //keeps track on the monomers which are added during the procedure + auto nAddedMonomers(0); + for(uint32_t i=0;i0 ) + start=(ingredients.getMolecules().size()-1); + std::cout << "Start: " << start < +bool UpdaterAffineDeformation::execute(){ + std::cout << "UpdaterAffineDeformation::initialize():"<< std::endl; + //adjusting the box size is not neccessary, because it is used only once in the FeatureCrosslinkConnections* + //there the jump vectors are calculated + + //The jump vector needs to be adjusted according the deformation labmda! + auto CrossLinkIDs(ing.getCrosslinkIDs()); + for (uint32_t i =0 ; i Neighbors(ing.getCrossLinkNeighborIDs(ID) ); + int32_t number_of_neighbors(Neighbors.size()); + if (number_of_neighbors > 0) { + for (size_t j = 0; j < number_of_neighbors; j++){ + if (i < 20 ) + std::cout << "ID= "<< i << " initial jump " << Neighbors[j].jump << " "; + ing.setCrossLinkNeighborJump(ID,j,deform(Neighbors[j].jump)); + if (i < 20 ) + std::cout << "ID= "<< i << " deformed jump " << ing.getCrossLinkNeighborIDs(ID)[j].jump << "\n"; + } + } + } + + //The initial posistions of all monomers needs to be adjusted according the deformation labmda! + for(size_t i = 0; i< ing.getMolecules().size();i++){ + if (i < 20 ) + std::cout << "ID="<< i << "initial position " << ing.getMolecules()[i].getVector3D() << " "; + ing.modifyMolecules()[i].modifyVector3D()=deform(ing.getMolecules()[i].getVector3D()); + if (i < 20 ) + std::cout << " deformed position " << ing.getMolecules()[i].getVector3D() << "\n"; + } + std::cout << "UpdaterAffineDeformation::initialize():done.\n"; +} +#endif /*LEMONADE_PM_UPDATER_UPDATERAFFINEDEFORMATION_H*/ \ No newline at end of file diff --git a/include/LeMonADE_PM/updater/UpdaterForceBalancedPosition.h b/include/LeMonADE_PM/updater/UpdaterForceBalancedPosition.h index 223ded4..bdd70e0 100644 --- a/include/LeMonADE_PM/updater/UpdaterForceBalancedPosition.h +++ b/include/LeMonADE_PM/updater/UpdaterForceBalancedPosition.h @@ -40,10 +40,10 @@ class UpdaterForceBalancedPosition:public AbstractUpdater { public: //! constructor for UpdaterForceBalancedPosition - UpdaterForceBalancedPosition(IngredientsType& ing_, double threshold_ ): - ing(ing_),threshold(threshold_){}; + UpdaterForceBalancedPosition(IngredientsType& ing_, double threshold_ , double decreaseFactor_=1.0): + ing(ing_),threshold(threshold_),decreaseFactor(decreaseFactor_){}; - virtual void initialize(){}// init(move);}; + virtual void initialize(){}; bool execute(); virtual void cleanup(){}; @@ -55,24 +55,15 @@ class UpdaterForceBalancedPosition:public AbstractUpdater //! threshold for the certainty double threshold; - + //! move to equilibrate the cross links by force equilibrium moveType move; //! random number generator RandomNumberGenerators rng; - - // //! input filename for the MoveNonLinearForceEquilibrium - // std::string filename; - // //! relaxation parameter - // double relaxationParameter; - // //!initialize of the move : MoveForceEquilibrium - // void init ( MoveForceEquilibrium& move ) {std::cout << "Nothing to initialize for this move type.\n";} - // //!initialize of the move :MoveNonLinearForceEquilibrium - // void init (MoveNonLinearForceEquilibrium& move ) { - // move.setRelaxationParameter(relaxationParameter); - // move.setFilename(filename); - // }; + + //! + double decreaseFactor; }; template @@ -96,13 +87,15 @@ bool UpdaterForceBalancedPosition::execute(){ NSuccessfulMoves++; } } - if( NSuccessfulMoves>0 ){ - avShift/=(NSuccessfulMoves); - }else - avShift=threshold*1.1; ing.modifyMolecules().setAge(ing.getMolecules().getAge()+1); - if (ing.getMolecules().getAge() %1000 == 0 ) + if (ing.getMolecules().getAge() %1000 == 0 ){ std::cout << "MCS: " << ing.getMolecules().getAge() << " and average shift: " << avShift << std::endl; + + } + if (ing.getMolecules().getAge() %100 == 0 ){ + setRelaxationParameter(move.getRelaxationParameter()*decreaseFactor); + // threshold*=decreaseFactor; + } } std::cout << "Finish equilibration with average shift per cross link < " << avShift << " after " << ing.getMolecules().getAge()-StartMCS <. --------------------------------------------------------------------------------*/ +#ifndef LEMONADE_PM_UPDATER_UPDATERREADCROSSLINKCONNECTIONS_H +#define LEMONADE_PM_UPDATER_UPDATERREADCROSSLINKCONNECTIONS_H #include #include #include @@ -90,7 +92,7 @@ class UpdaterReadCrosslinkConnections : public AbstractUpdater uint32_t nExecutions; //! connects the cross link to the chain - void ConnectCrossLinkToChain(uint32_t MonID, uint32_t chainID); + bool ConnectCrossLinkToChain(uint32_t MonID, uint32_t chainID); //!bond Table std::map,std::vector > bondTable; @@ -98,18 +100,31 @@ class UpdaterReadCrosslinkConnections : public AbstractUpdater template -void UpdaterReadCrosslinkConnections::ConnectCrossLinkToChain(uint32_t MonID, uint32_t chainID){ +bool UpdaterReadCrosslinkConnections::ConnectCrossLinkToChain(uint32_t MonID, uint32_t chainID){ + // std::pair key(MonID,chainID-1); + // if (bondTable.find(key) != bondTable.end()){ + // if ( !ing.getMolecules().areConnected( MonID,bondTable.at(key)[0] ) ) + // ing.modifyMolecules().connect(MonID,bondTable.at(key)[0] ); + // else + // ing.modifyMolecules().connect(MonID,bondTable.at(key)[1] ); + // }else{ + // std::stringstream errormessage; + // errormessage << "There was no such a connection in the bfm file for monomer ID= " << MonID << " with chainID=" << chainID-1<< "\n"; + // throw std::runtime_error(errormessage.str()); + // } std::pair key(MonID,chainID-1); if (bondTable.find(key) != bondTable.end()){ - if ( !ing.getMolecules().areConnected( MonID,bondTable.at(key)[0] ) ) + if ( !ing.getMolecules().areConnected( MonID,bondTable.at(key)[0] ) ){ ing.modifyMolecules().connect(MonID,bondTable.at(key)[0] ); - else + return true ; + } + //chain can only have one neighbor (without this statement a more or less random partner would be connected to the structure !!) + if ( bondTable.at(key).size()>1 ) { ing.modifyMolecules().connect(MonID,bondTable.at(key)[1] ); - }else{ - std::stringstream errormessage; - errormessage << "There was no such a connection in the bfm file for monomer ID= " << MonID << " with chainID=" << chainID-1<< "\n"; - throw std::runtime_error(errormessage.str()); + return true ; + } } + return false; } /** @@ -151,18 +166,6 @@ bool UpdaterReadCrosslinkConnections::execute(){ stream.open(input); if (stream.fail()) throw std::runtime_error(std::string("error opening input file ") + input + std::string("\n")); - bool findStartofData(false); - //set head to beginning of the data block - while (!findStartofData){ - std::string line; - getline(stream, line); - if (line.empty()) - findStartofData = true; - else if (line.at(0) == '#') //go to next line - continue; - else - throw std::runtime_error("Wrong input format!"); - } //current conversion auto conversion = minConversion + static_cast(nExecutions) * stepwidth; std::cout << "Current conversion is " <::execute(){ while (NewConnections < ReadNLines && stream.good()){ std::string line; getline(stream, line); - if (line.empty()) + if (line.empty() || line.at(0) == '#') break; std::stringstream ss; uint32_t Time, ChainID, MonID1, P1X, P1Y, P1Z, MonID2, P2X, P2Y, P2Z; ss << line; ss >> Time >> ChainID >> MonID1 >> P1X >> P1Y >> P1Z >> MonID2 >> P2X >> P2Y >> P2Z; - // if (NewConnections == 0 ) - // std::cout << Time << " " << ChainID<< " " - // << MonID1<< " " << P1X << " " << P1Y << " " << P1Z << " " - // << MonID2<< " " << P2X << " " << P2Y << " " << P2Z << "\n"; + #ifdef DEBUG + std::cout << ss.str() << std::endl; + #endif //DEBUG// ing.modifyMolecules().setAge(Time); - ConnectCrossLinkToChain(MonID1, ChainID); + // ConnectCrossLinkToChain(MonID1, ChainID); + if ( ConnectCrossLinkToChain(MonID1, ChainID+1) ) { + }else if( ConnectCrossLinkToChain(MonID2, ChainID+1) ) { + }else { + std::stringstream errormessage; + errormessage << "There was no such a connection in the bfm file for monomer ID= " << MonID1 <<" with ID=" << MonID2 << " with chainID=" << ChainID<< "\n"; + throw std::runtime_error(errormessage.str()); + } //update positions : I think this is not needed // ing.modifyMolecules()[MonID1].modifyVector3D().setAllCoordinates(P1X, P1Y, P1Z); // if (MonID2 > 0) @@ -208,3 +217,5 @@ bool UpdaterReadCrosslinkConnections::execute(){ return true; } } + +#endif \ No newline at end of file diff --git a/include/LeMonADE_PM/updater/UpdaterReadCrosslinkConnectionsTendomer.h b/include/LeMonADE_PM/updater/UpdaterReadCrosslinkConnectionsTendomer.h new file mode 100644 index 0000000..016d43d --- /dev/null +++ b/include/LeMonADE_PM/updater/UpdaterReadCrosslinkConnectionsTendomer.h @@ -0,0 +1,202 @@ +/*-------------------------------------------------------------------------------- + ooo L attice-based | + o\.|./o e xtensible | LeMonADE: An Open Source Implementation of the + o\.\|/./o Mon te-Carlo | Bond-Fluctuation-Model for Polymers +oo---0---oo A lgorithm and | + o/./|\.\o D evelopment | Copyright (C) 2013-2015 by + o/.|.\o E nvironment | LeMonADE Principal Developers (see AUTHORS) + ooo | +---------------------------------------------------------------------------------- + +This file is part of LeMonADE. + +LeMonADE is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LeMonADE is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LeMonADE. If not, see . + +--------------------------------------------------------------------------------*/ +#ifndef LEMONADE_PM_UPDATER_UPDATERREADCROSSLINKCONNECTIONSTENDOMER_H +#define LEMONADE_PM_UPDATER_UPDATERREADCROSSLINKCONNECTIONSTENDOMER_H +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * @class UpdaterReadCrosslinkConnectionsTendomer + * @brief reads in connections for a system made of linear chains and crosslinker + * + * @details + * -the input file needs to have a format like: + * #Time, ChainID, MonID1, P1X, P1Y, P1Z, MonID2, P2X, P2Y, P2Z + * -the chains are before the crosslinks in the bfm file + * -the chain length must be at least 1 + * + */ + +template +class UpdaterReadCrosslinkConnectionsTendomer : public AbstractUpdater +{ +public: + UpdaterReadCrosslinkConnectionsTendomer( + IngredientsType& ing_, + const std::string input_, + const double stepwidth_, + const double minConversion_): + ing(ing_), + input(input_), + stepwidth(stepwidth_), + minConversion(minConversion_), + nExecutions(0){}; + virtual void initialize(); + virtual bool execute(); + virtual void cleanup(){}; + +private: + //! container storing system information about monomers + IngredientsType& ing; + + //! container storing system information about monomers + IngredientsType initialIng; + + //! input of the cross link positions and so on.. + const std::string input; + + //!conversion step used in the execute function to increase the conversion + const double stepwidth; + + //! minimal conversion to be read in before further analyse + double minConversion; + + //!number of maximum connections for a cross link + uint32_t NMaxConnection; + + //!numbe of monomers per chain + uint32_t NMonomerPerChain; + + //!number of executions; + uint32_t nExecutions; + + //! connects the cross link to the chain + bool ConnectCrossLinkToChain(uint32_t MonID, uint32_t chainID); + + //!bond Table + std::map,std::vector > bondTable; +}; +template +bool UpdaterReadCrosslinkConnectionsTendomer::ConnectCrossLinkToChain(uint32_t MonID, uint32_t chainID){ + std::pair key(MonID,chainID-1); + if (bondTable.find(key) != bondTable.end()){ + if ( !ing.getMolecules().areConnected( MonID,bondTable.at(key)[0] ) ){ + ing.modifyMolecules().connect(MonID,bondTable.at(key)[0] ); + return true ; + } + //chain can only have one neighbor (without this statement a more or less random partner would be connected to the structure !!) + if ( bondTable.at(key).size()>1 ) { + ing.modifyMolecules().connect(MonID,bondTable.at(key)[1] ); + return true ; + } + } + return false; +} +/** + * @brief read in connections up to the minimum conversion given + * */ +template +void UpdaterReadCrosslinkConnectionsTendomer::initialize(){ + //assume a stochiometric mixture + NMaxConnection=4*ing.getNumCrossLinkers(); + NMonomerPerChain = 2*ing.getNumMonomersPerChain(); + std::cout << "Number of maximum connection: " << NMaxConnection << std::endl; + //erase bonds between reactive monomers + for (uint32_t i =0 ; i < ing.getMolecules().size(); i++) + if (ing.getMolecules()[i].isReactive() ) + for(size_t j=0; j < ing.getMolecules().getNumLinks(i);j++){ + uint32_t neighbor(ing.getMolecules().getNeighborIdx(i,j)); + if (ing.getMolecules()[neighbor].isReactive()){ + ing.modifyMolecules().disconnect(i, neighbor ); + uint32_t chainMonomer(std::min(i,neighbor) ); + uint32_t crosslink(std::max(i,neighbor)); + uint32_t chainID( (chainMonomer-chainMonomer%NMonomerPerChain)/NMonomerPerChain); + bondTable[std::pair(crosslink,chainID) ].push_back(chainMonomer) ; + } + } + std::cout << "Erase " << bondTable.size() << " bonds." < +bool UpdaterReadCrosslinkConnectionsTendomer::execute(){ + //reset the ingredients container to the inital one + ing = initialIng; + //stream reading input file + std::ifstream stream; + stream.open(input); + if (stream.fail()) + throw std::runtime_error(std::string("error opening input file ") + input + std::string("\n")); + //current conversion + auto conversion = minConversion + static_cast(nExecutions) * stepwidth; + std::cout << "Current conversion is " <> Time >> createBreak>> ChainID >> nSegments >>MonID1 >> P1X >> P1Y >> P1Z >> MonID2 >> P2X >> P2Y >> P2Z; + #ifdef DEBUG + std::cout << ss.str() << std::endl; + #endif //DEBUG// + ing.modifyMolecules().setAge(Time); + if ( ConnectCrossLinkToChain(MonID1, ChainID+1) ) { + }else if( ConnectCrossLinkToChain(MonID2, ChainID+1) ) { + }else { + std::stringstream errormessage; + errormessage << "There was no such a connection in the bfm file for monomer ID= " << MonID1 <<" of with ID=" << MonID2 << " with chainID=" << ChainID<< "\n"; + throw std::runtime_error(errormessage.str()); + } + NewConnections++; + } + std::cout << "Read and add " << NewConnections << "/" << NMaxConnection + << " to the system at time " + << ing.getMolecules().getAge() << std::endl; + ing.synchronize(); + nExecutions++; + std::cout << "UpdaterReadCrosslinkConnectionsTendomer::execute " << nExecutions << " times.\n"; + //close the filestream and return false if the file has ended and thus the updater has nothing more to do + if (stream.eof()) { + stream.close(); + return false; + }else{ + stream.close(); + return true; + } +} + +#endif \ No newline at end of file diff --git a/include/LeMonADE_PM/updater/moves/MoveForceEquilibrium.h b/include/LeMonADE_PM/updater/moves/MoveForceEquilibrium.h index 4909893..47bd8c0 100644 --- a/include/LeMonADE_PM/updater/moves/MoveForceEquilibrium.h +++ b/include/LeMonADE_PM/updater/moves/MoveForceEquilibrium.h @@ -47,7 +47,9 @@ along with LeMonADE. If not, see . class MoveForceEquilibrium:public MoveForceEquilibriumBase{ public: - MoveForceEquilibrium():bondlength(2.68){}; + MoveForceEquilibrium():bondlength2(2.68*2.68){ + std::cout << "Use the MoveForceEquilibrium\n"; + }; // overload initialise function to be able to set the moves index and direction if neccessary template void init(const IngredientsType& ing); @@ -57,19 +59,28 @@ class MoveForceEquilibrium:public MoveForceEquilibriumBase template bool check(IngredientsType& ing); template< class IngredientsType> void apply(IngredientsType& ing); + void setFilename(std::string filename_){} + //! get the filename for the force extension data + std::string const getFilename(){} + + //! set the relaxation parameter for the cross link + void setRelaxationParameter(double relaxationChain_){} + //! get the relaxation parameter for the cross link + double getRelaxationParameter(){} private: //average square bond length - const double bondlength; + const double bondlength2; + // //Gaussina force extension relation //f=R*3/(N*b^2) VectorDouble3 FE(VectorDouble3 extensionVector, double nSegs){ - return extensionVector*3./((nSegs)*bondlength*bondlength); + return extensionVector*3./((nSegs)*bondlength2); } //Gaussian extension force relation //R=-f/3*N*b^2 VectorDouble3 EF(VectorDouble3 force, double nSegs){ - return force/(3.)*(nSegs)*bondlength*bondlength; + return force/(3.)*(nSegs)*bondlength2; } //calculate the shift for the cross link @@ -85,6 +96,7 @@ class MoveForceEquilibrium:public MoveForceEquilibriumBase VectorDouble3 vec(ing.getMolecules()[Neighbors[i].ID].getVector3D()-Position-Neighbors[i].jump); avNSegments+=1./Neighbors[i].segDistance; force+=FE(vec,Neighbors[i].segDistance); + // std::cout <<"MoveLFE " << ing.getMolecules()[Neighbors[i].ID].getVector3D() << "\t" << Neighbors[i].jump<< "\t" << Position << "\t" << vec << "\t" << force << std::endl; } shift=EF(force,1./avNSegments); } diff --git a/include/LeMonADE_PM/updater/moves/MoveForceEquilibriumBase.h b/include/LeMonADE_PM/updater/moves/MoveForceEquilibriumBase.h index fe42bcc..24a325b 100644 --- a/include/LeMonADE_PM/updater/moves/MoveForceEquilibriumBase.h +++ b/include/LeMonADE_PM/updater/moves/MoveForceEquilibriumBase.h @@ -88,6 +88,17 @@ class MoveForceEquilibriumBase:public MoveBase { ShiftVector.setAllCoordinates(dx,dy,dz); } +public: + void setFilename(std::string filename_){static_cast(this)->setFilename(filename_);} + //! get the filename for the force extension data + std::string const getFilename(){static_cast(this)->getFilename();} + + //! set the relaxation parameter for the cross link + void setRelaxationParameter(double relaxationChain_){static_cast(this)->setRelaxationParameter(relaxationChain_);} + //! get the relaxation parameter for the cross link + double getRelaxationParameter(){static_cast(this)->getRelaxationParameter();} + + //! Random Number Generator (RNG) RandomNumberGenerators randomNumbers; diff --git a/include/LeMonADE_PM/updater/moves/MoveNonLinearForceEquilibrium.h b/include/LeMonADE_PM/updater/moves/MoveNonLinearForceEquilibrium.h index 973fc96..8f9d9ed 100644 --- a/include/LeMonADE_PM/updater/moves/MoveNonLinearForceEquilibrium.h +++ b/include/LeMonADE_PM/updater/moves/MoveNonLinearForceEquilibrium.h @@ -29,6 +29,7 @@ along with LeMonADE. If not, see . #define LEMONADE_PM_UPDATER_MOVES_MOVENONLINEARFORCEEQUILIBRIUM_H #include #include +#include #include #include #include @@ -50,7 +51,7 @@ class MoveNonLinearForceEquilibrium:public MoveForceEquilibriumBase(round(extensionVector.getLength())) ); - if ( max_extension < length ){ + if (extensionVector == VectorDouble3(0.,0.,0.) ) + return VectorDouble3(0.,0.,0.); + double length( extensionVector.getLength() ); + if ( max_extension < length ){ + return EFGauss(extensionVector); + } + #ifdef DEBUG + if ( max_extension < length ){ std::stringstream errormessage; errormessage << "The length of the extension vector is greater than the maximum given in the file: \n" << "length is " << length @@ -92,8 +97,25 @@ class MoveNonLinearForceEquilibrium:public MoveForceEquilibriumBase(round(length/accuracy))] << "\t" << EFGauss(extensionVector).getLength() << "\t" + << extensionVector << "\t" + << EFGauss(extensionVector) << "\t"<< extensionVector.normalize()*(force_extension[static_cast(round(length/accuracy))])<<"\t" + << "\n"; + #endif + if (length == 0) + return VectorDouble3(0.,0.,0.); + auto x(length/accuracy); + auto down(static_cast(floor(x))); + auto up (down+1); + auto amplitude=force_extension[down] + (force_extension[up]-force_extension[down])/static_cast(up-down)*(x-down)/static_cast(up-down); + return extensionVector.normalize()*(amplitude); + } + //Gaussina force extension relation + //f=R*3/(N*b^2) + VectorDouble3 EFGauss(VectorDouble3 extensionVector){ + return extensionVector/springConstant; } //Gaussian extension force relation //R=-f/3*N*b^2 @@ -121,10 +143,13 @@ class MoveNonLinearForceEquilibrium:public MoveForceEquilibriumBase extension_force; //extension force mapping index=extension rounded to int std::vector force_extension; @@ -137,64 +162,89 @@ class MoveNonLinearForceEquilibrium:public MoveForceEquilibriumBase 0) { + int32_t number_of_neighbors(Neighbors.size()); + if (number_of_neighbors > 0) { VectorDouble3 Position(ing.getMolecules()[this->getIndex()].getVector3D()); - for (size_t i = 0; i < Neighbors.size(); i++){ - VectorDouble3 vec(Position-ing.getMolecules()[Neighbors[i].ID].getVector3D()-Neighbors[i].jump); - force+=EF(vec);//Neighbors[i].segDistance + for (size_t i = 0; i < number_of_neighbors; i++){ + VectorDouble3 vec(ing.getMolecules()[Neighbors[i].ID].getVector3D()-Neighbors[i].jump-Position); + // std::cout <<"MoveNLFE " << ing.getMolecules()[Neighbors[i].ID].getVector3D() << " " << Neighbors[i].jump<< " " << Position <getIndex()<<" "<< Neighbors[i].ID<< " " <(Neighbors.size()) )); + shift=FE(force/(static_cast(number_of_neighbors) )); } return shift; }; + + //! check is file exists: + inline bool fileExists (const std::string& name) { + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); + } }; ///////////////////////////////////////////////////////////////////////////// /////////// implementation of the members /////////////////////////////////// void MoveNonLinearForceEquilibrium::createTable(){ - std::ifstream in(filename); - uint32_t counter(0); - while(in.good() && in.peek()!=EOF){ - std::string line; - getline(in, line); - //ignore comments and blank lines - while (line.at(0) == '#' || line.empty() ) //go to next line - continue; - //read data - double force, extension; - std::stringstream ss ; - ss<< line; - ss>>force >> extension; - if(min_force > force ) min_force=force; - if(max_force < force ) max_force=force; - if(min_extension > extension ) min_extension=extension; - if(max_extension < extension ) max_extension=extension; - extension_force.insert(extension_force.end(),std::pair(force, extension)); - if(counter==1){ - force_steps=max_force-min_force; - }else if(counter>1) - counter++; - } - in.close(); - //make lookup for the extension force relation - //make a entry from 0 to max_extension in steps of 1 - for ( auto r=0;r(max_extension); r++ ){ - if(r==0) - force_extension.push_back(0.); - else{ - auto it_last=extension_force.begin(); - for (auto it=extension_force.begin(); it !=extension_force.end();it++){ - if(it->second > r){ - //at the force linear interpolated in between the two current forces - auto deltaForce(it->first-it_last->first); - auto deltaExtension(it->second-it_last->second); - auto factor( (static_cast(r)-it_last->second)/deltaExtension ); - //at interpolated force - force_extension.push_back(it_last->first+deltaForce*factor); - break; + if (fileExists(filename )){ + std::ifstream in(filename); + uint32_t counter(0); + extension_force.insert(extension_force.end(),std::pair(0., 0.)); + min_force=0.; + min_extension=0.; + while(in.good() && in.peek()!=EOF){ + std::string line; + getline(in, line); + //ignore comments and blank lines + if (line.at(0) == '#' || line.empty() ) //go to next line + continue; + //read data + double force, extension; + std::stringstream ss ; + ss<< line; + ss>>force >> extension; + if(min_force > force ) min_force=force; + if(max_force < force ) max_force=force; + if(min_extension > extension ) min_extension=extension; + if(max_extension < extension ) max_extension=extension; + extension_force.insert(extension_force.end(),std::pair(force, extension)); + if(counter==1){ + force_steps=max_force-min_force; + }else if(counter>1) + counter++; + } + in.close(); + //make lookup for the extension force relation + //make a entry from 0 to max_extension in steps of 1 + auto it_last=extension_force.begin(); + for ( auto r=0;r(max_extension/accuracy); r++ ){ + if(r==0) + force_extension.push_back(0.); + else{ + for (auto it=it_last; it !=extension_force.end();it++){ + if(it->second > static_cast(r*accuracy)){ + //at the force linear interpolated in between the two current forces + auto deltaForce(it->first-it_last->first); + auto deltaExtension(it->second-it_last->second); + auto factor( (static_cast(r*accuracy)-it_last->second)/deltaExtension ); + //at interpolated force + force_extension.push_back(it_last->first+deltaForce*factor); + break; + } + it_last=it; } - it_last=it; } } + std::cout << "MoveNonLinearForceEquilibrium::createTable() force extension" <. + +--------------------------------------------------------------------------------*/ + +/****************************************************************************** + * based on LeMonADE: https://github.com/LeMonADE-project/LeMonADE/ + * author: Toni Müller + * email: mueller-toni@ipfdd.de + * project: LeMonADE-Phantom Modulus + *****************************************************************************/ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + + +int main(int argc, char* argv[]){ + try{ + /////////////////////////////////////////////////////////////////////////////// + ///parse options/// + std::string inputBFM("init.bfm"); + std::string suffix("MolecularWeight.dat"); + bool showHelp = false; + auto parser + = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file").required() + | clara::detail::Opt( suffix, "suffix (=MolecularWeight.dat)" ) ["-o"]["--suffix" ] ("suffix name for data." ).required() + | clara::Help( showHelp ); + + auto result = parser.parse( clara::Args( argc, argv ) ); + + if( !result ) { + std::cerr << "Error in command line: " << result.errorMessage() << std::endl; + exit(1); + }else if(showHelp == true){ + parser.writeToStream(std::cout); + exit(0); + }else{ + std::cout << "inputBFM : " << inputBFM << std::endl; + std::cout << "suffix : " << suffix << std::endl; + + } + /////////////////////////////////////////////////////////////////////////////// + ///end options parsing + /////////////////////////////////////////////////////////////////////////////// + typedef LOKI_TYPELIST_2(FeatureMoleculesIOUnsaveCheck, FeatureReactiveBonds) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing ingredients; + + TaskManager taskmanager; + + taskmanager.addUpdater( new UpdaterReadBfmFile(inputBFM,ingredients, UpdaterReadBfmFile::READ_LAST_CONFIG_SAVE),0); + taskmanager.addAnalyzer( new AnalyzerMolecularWeight(ingredients, suffix) ); + //initialize and run + taskmanager.initialize(); + taskmanager.run(); + taskmanager.cleanup(); + std::cout << "Read in conformation and go on to bring it into equilibrium forces..." <. #include #include #include - +#include int main(int argc, char* argv[]){ try{ @@ -61,48 +61,66 @@ int main(int argc, char* argv[]){ std::string inputBFM("init.bfm"); std::string outputDataPos("CrosslinkPosition.dat"); std::string outputDataDist("ChainExtensionDistribution.dat"); - std::string inputConnection("BondCreationBreaking.dat"); - std::string feCurve(""); + // std::string inputConnection("BondCreationBreaking.dat"); + std::string feCurve; double relaxationParameter(10.); double threshold(0.5); double stepwidth(1.0); double minConversion(50.0); - bool custom(false); + bool custom(true); + double stretching_factor(1.0); + double dampingfactor(1.0); + double prestrainFactorX(1.0); + double prestrainFactorY(1.0); + double prestrainFactorZ(1.0); bool showHelp = false; auto parser - = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() - | clara::detail::Opt( inputConnection, "inputConnection (=BondCreationBreaking.dat)" ) ["-d"]["--inputConnection"] ("used for the time development of the topology. " ).required() - | clara::detail::Opt( outputDataPos, "outputDataPos (=CrosslinkPosition.dat)" ) ["-o"]["--outputPos" ] ("(optional) Output filename of the crosslink ID and the equilibrium Position.").optional() - | clara::detail::Opt( outputDataDist, "outputDataDist (=ChainExtensionDistribution.dat)") ["-c"]["--outputDist" ] ("(optional) Output filename of the chain extension distribution." ).optional() - | clara::detail::Opt( stepwidth, "stepwidth" ) ["-s"]["--stepwidth" ] ("(optional) Width for the increase in percentage. Default: 1%." ).optional() - | clara::detail::Opt( minConversion, "minConversion" ) ["-u"]["--minConversion" ] ("(optional) Minimum conversion to be read in. Default: 50%." ).optional() - | clara::detail::Opt( threshold, "threshold" ) ["-t"]["--threshold" ] ("(optional) Threshold of the average shift. Default 0.5 ." ).optional() - | clara::detail::Opt( feCurve, "feCurve (="")" ) ["-f"]["--feCurve" ] ("(optional) Force-Extension curve. Default \"\"." ).optional() - | clara::detail::Opt( relaxationParameter, "relaxationParameter (=10)" ) ["-r"]["--relax" ] ("(optional) Relaxation parameter. Default 10.0 ." ).optional() + = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() + // | clara::detail::Opt( inputConnection, "inputConnection (=BondCreationBreaking.dat)" ) ["-d"]["--inputConnection"] ("used for the time development of the topology. " ).required() + | clara::detail::Opt( outputDataPos, "outputDataPos (=CrosslinkPosition.dat)" ) ["-o"]["--outputPos" ] ("(optional) Output filename of the crosslink ID and the equilibrium Position.").optional() + | clara::detail::Opt( outputDataDist, "outputDataDist (=ChainExtensionDistribution.dat)") ["-c"]["--outputDist" ] ("(optional) Output filename of the chain extension distribution." ).optional() + | clara::detail::Opt( stepwidth, "stepwidth" ) ["-s"]["--stepwidth" ] ("(optional) Width for the increase in percentage. Default: 1%." ).optional() + | clara::detail::Opt( minConversion, "minConversion" ) ["-u"]["--minConversion" ] ("(optional) Minimum conversion to be read in. Default: 50%." ).optional() + | clara::detail::Opt( threshold, "threshold" ) ["-t"]["--threshold" ] ("(optional) Threshold of the average shift. Default 0.5 ." ).optional() + | clara::detail::Opt( stretching_factor, "stretching_factor (=1)" ) ["-l"]["--stretching_factor"] ("(optional) Stretching factor for uniaxial deformation. Default 1.0 ." ).optional() + | clara::detail::Opt( feCurve, "feCurve (="")" ) ["-f"]["--feCurve" ] ("(optional) Force-Extension curve. Default \"\"." ).optional() + | clara::detail::Opt( relaxationParameter, "relaxationParameter (=10)" ) ["-r"]["--relax" ] ("(optional) Relaxation parameter. Default 10.0 ." ).optional() + | clara::detail::Opt( dampingfactor, "damping (=1)" ) ["-d"]["--damping" ] ("(optional) Damping factor after 1E3MCS. Default 1.0." ).optional() + | clara::detail::Opt( prestrainFactorX, "prestrainFactorX (=1)" ) ["-x"]["--prestrainFactorX" ] ("(optional) Prestrain factor in X. Default 1.0." ).optional() + | clara::detail::Opt( prestrainFactorY, "prestrainFactorY (=1)" ) ["-y"]["--prestrainFactorY" ] ("(optional) Prestrain factor in Y. Default 1.0." ).optional() + | clara::detail::Opt( prestrainFactorZ, "prestrainFactorZ (=1)" ) ["-z"]["--prestrainFactorZ" ] ("(optional) Prestrain factor in Z. Default 1.0." ).optional() | clara::Help( showHelp ); auto result = parser.parse( clara::Args( argc, argv ) ); - + if ( feCurve.empty() ) custom=false; if( !result ) { std::cerr << "Error in command line: " << result.errorMessage() << std::endl; exit(1); }else if(showHelp == true){ - std::cout << "Simulator to connect linear chains with single monomers of certain functionality"<< std::endl; + std::cout << "Standard force equilibration for a end-linked network."<< std::endl; parser.writeToStream(std::cout); exit(0); }else{ std::cout << "outputData : " << outputDataPos << std::endl; std::cout << "outputDataDist : " << outputDataDist << std::endl; std::cout << "inputBFM : " << inputBFM << std::endl; - std::cout << "inputConnection : " << inputConnection << std::endl; + std::cout << "custom FE : " << custom << std::endl; std::cout << "stepwidth : " << stepwidth << std::endl; + std::cout << "dampingfactor : " << dampingfactor << std::endl; std::cout << "minConversion : " << minConversion << std::endl; std::cout << "threshold : " << threshold << std::endl; std::cout << "feCurve : " << feCurve << std::endl; + std::cout << "stretching_factor : " << stretching_factor << std::endl; + std::cout << "prestrainFactorX : " << prestrainFactorX << std::endl; + std::cout << "prestrainFactorY : " << prestrainFactorY << std::endl; + std::cout << "prestrainFactorZ : " << prestrainFactorZ << std::endl; } - if (! feCurve.empty()) custom=true; + + + RandomNumberGenerators rng; + rng.seedAll(); /////////////////////////////////////////////////////////////////////////////// ///end options parsing /////////////////////////////////////////////////////////////////////////////// @@ -152,24 +170,31 @@ int main(int argc, char* argv[]){ } } myIngredients2.synchronize(); - - TaskManager taskmanager2; - //read bonds and positions stepwise - taskmanager2.addUpdater( new UpdaterReadCrosslinkConnections(myIngredients2, inputConnection, stepwidth, minConversion) ); - if (custom) - taskmanager2.addUpdater( new UpdaterForceBalancedPosition(myIngredients2, threshold) ); - else { - auto updater = new UpdaterForceBalancedPosition(myIngredients2, threshold) ; - updater->setFilename(feCurve); - updater->setRelaxationParameter(relaxationParameter); - taskmanager2.addUpdater( updater ); - } - taskmanager2.addAnalyzer(new AnalyzerEquilbratedPosition(myIngredients2,outputDataPos, outputDataDist)); - //initialize and run + + auto forceUpdater = new UpdaterForceBalancedPosition(myIngredients2, threshold,dampingfactor); + if(custom) + forceUpdater->setFilename(feCurve); + forceUpdater->setRelaxationParameter(relaxationParameter); + auto forceUpdater2 = new UpdaterForceBalancedPosition(myIngredients2, threshold,dampingfactor); + auto uniaxialDeformation = new UpdaterAffineDeformation(myIngredients2, stretching_factor,prestrainFactorX,prestrainFactorY,prestrainFactorZ); + + auto analyzer = new AnalyzerEquilbratedPosition(myIngredients2,outputDataPos,outputDataDist); + + TaskManager taskmanager2; + taskmanager2.addUpdater( uniaxialDeformation,0 ); + //read bonds and positions stepwise + if(custom){ + std::cout << "Use custom force-extension curve\n"; + taskmanager2.addUpdater( forceUpdater ); + }else{ + std::cout << "Use gaussian force-extension relation\n"; + taskmanager2.addUpdater( forceUpdater2 ); + } + taskmanager2.addAnalyzer(analyzer); + //initialize and run taskmanager2.initialize(); - taskmanager2.run(); + taskmanager2.run(1); taskmanager2.cleanup(); - } catch(std::exception& e){ std::cerr<<"Error:\n" diff --git a/projects/IdealReference2ForceEquilibrium.cpp b/projects/IdealReference2ForceEquilibrium.cpp new file mode 100644 index 0000000..40a1f38 --- /dev/null +++ b/projects/IdealReference2ForceEquilibrium.cpp @@ -0,0 +1,201 @@ +/*-------------------------------------------------------------------------------- + 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 . + +--------------------------------------------------------------------------------*/ + +/****************************************************************************** + * based on LeMonADE: https://github.com/LeMonADE-project/LeMonADE/ + * author: Toni Müller + * email: mueller-toni@ipfdd.de + * project: LeMonADE-Phantom Modulus + *****************************************************************************/ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +int main(int argc, char* argv[]){ + try{ + /////////////////////////////////////////////////////////////////////////////// + ///parse options/// + // std::string inputBFM("init.bfm"); + std::string outputDataPos("CrosslinkPosition.dat"); + std::string outputDataDist("ChainExtensionDistribution.dat"); + std::string feCurve(""); + double relaxationParameter(10.); + double threshold(0.5); + double factor(0.995); + double stretching_factor(1.0); + uint32_t gauss(0); + + uint32_t nSegments(16); + uint32_t functionality(4); + uint32_t nRings(0); + + bool showHelp = false; + auto parser + // = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() + = clara::detail::Opt( outputDataPos, "outputDataPos (=CrosslinkPosition.dat)" ) ["-o"]["--outputPos" ] ("(optional) Output filename of the crosslink ID and the equilibrium Position.").optional() + | clara::detail::Opt( outputDataDist, "outputDataDist (=ChainExtensionDistribution.dat)") ["-c"]["--outputDist" ] ("(optional) Output filename of the chain extension distribution." ).optional() + | clara::detail::Opt( threshold, "threshold" ) ["-t"]["--threshold" ] ("(optional) Threshold of the average shift. Default 0.5 ." ).optional() + | clara::detail::Opt( stretching_factor, "stretching_factor (=1)" ) ["-l"]["--stretching_factor"] ("(optional) Stretching factor for uniaxial deformation. Default 1.0 ." ).optional() + | clara::detail::Opt( feCurve, "feCurve (="")" ) ["-f"]["--feCurve" ] ("(optional) Force-Extension curve. Default \"\"." ).required() + | clara::detail::Opt( relaxationParameter, "relaxationParameter (=10)" ) ["-r"]["--relax" ] ("(optional) Relaxation parameter. Default 10.0 ." ).optional() + | clara::detail::Opt( gauss, "gauss" ) ["-g"]["--gauss" ] ("(optional) Deforma with a Gaussian deformation behaviour. Default 1.0 ." ).optional() + | clara::detail::Opt( nSegments, "nSegments" ) ["-n"]["--nSegments" ] ("(optional) Number of segments for the strand." ).optional() + | clara::detail::Opt( functionality, "nStrands" ) ["-s"]["--nStrands" ] ("(optional) Functionality." ).optional() + | clara::detail::Opt( nRings, "nRings" ) ["-m"]["--nRings" ] ("(optional) number of rings." ).optional() + | clara::Help( showHelp ); + + auto result = parser.parse( clara::Args( argc, argv ) ); + + if( !result ) { + std::cerr << "Error in command line: " << result.errorMessage() << std::endl; + exit(1); + }else if(showHelp == true){ + std::cout << "Force equilibration of two a double star, where each chain is replaced effectively by a tendomer."<< std::endl; + parser.writeToStream(std::cout); + exit(0); + }else{ + std::cout << "outputData : " << outputDataPos << std::endl; + std::cout << "outputDataDist : " << outputDataDist << std::endl; + std::cout << "threshold : " << threshold << std::endl; + std::cout << "feCurve : " << feCurve << std::endl; + std::cout << "gauss : " << gauss << std::endl; + std::cout << "stretching_factor : " << stretching_factor << std::endl; + std::cout << "nSegments : " << nSegments << std::endl; + std::cout << "functionality : " << functionality << std::endl; + std::cout << "nRings : " << nRings << std::endl; + } + RandomNumberGenerators rng; + // rng.seedDefaultValuesAll(); + rng.seedAll(); + /////////////////////////////////////////////////////////////////////////////// + ///end options parsing + /////////////////////////////////////////////////////////////////////////////// + //Read in th last Config + typedef LOKI_TYPELIST_2(FeatureFixedMonomers,FeatureMoleculesIO) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + myIngredients.setBoxX(256); + myIngredients.setBoxY(256); + myIngredients.setBoxZ(256); + myIngredients.setPeriodicX(1); + myIngredients.setPeriodicY(1); + myIngredients.setPeriodicZ(1); + myIngredients.modifyBondset().addBFMclassicBondset(); + myIngredients.synchronize(); + TaskManager taskmanager; + taskmanager.addUpdater( new UpdaterAddTMDoubleStars(myIngredients,1, nSegments, nRings, functionality ),0); + taskmanager.addAnalyzer(new AnalyzerWriteBfmFile("config.bfm", myIngredients, AnalyzerWriteBfmFile::APPEND) ); + taskmanager.initialize(); + taskmanager.run(1); + taskmanager.cleanup(); + + std::cout << "Read in conformation and go on to bring it into equilibrium forces..." < Config2; + typedef Ingredients Ing2; + Ing2 myIngredients2; + + myIngredients2.setBoxX(myIngredients.getBoxX()); + myIngredients2.setBoxY(myIngredients.getBoxY()); + myIngredients2.setBoxZ(myIngredients.getBoxZ()); + myIngredients2.setPeriodicX(myIngredients.isPeriodicX()); + myIngredients2.setPeriodicY(myIngredients.isPeriodicY()); + myIngredients2.setPeriodicZ(myIngredients.isPeriodicZ()); + myIngredients2.modifyMolecules().resize(myIngredients.getMolecules().size()); + myIngredients2.modifyMolecules().setAge(myIngredients.getMolecules().getAge()); + + for(size_t i = 0; i< myIngredients.getMolecules().size();i++){ + myIngredients2.modifyMolecules()[i].setMovableTag(myIngredients.modifyMolecules()[i].getMovableTag()); + myIngredients2.modifyMolecules()[i].modifyVector3D()=myIngredients.getMolecules()[i].getVector3D(); + for (size_t j = 0 ; j < myIngredients.getMolecules().getNumLinks(i);j++){ + uint32_t neighbor(myIngredients.getMolecules().getNeighborIdx(i,j)); + if( ! myIngredients2.getMolecules().areConnected(i,neighbor) ) + myIngredients2.modifyMolecules().connect(i,neighbor); + } + } + + myIngredients2.synchronize(); + + TaskManager taskmanager2; + taskmanager2.addUpdater( new UpdaterAffineDeformation(myIngredients2, stretching_factor),0 ); + //read bonds and positions stepwise + auto updater = new UpdaterForceBalancedPosition(myIngredients2, threshold,0.95) ; + updater->setFilename(feCurve); + updater->setRelaxationParameter(relaxationParameter); + auto updater2 = new UpdaterForceBalancedPosition(myIngredients2, threshold) ; + if ( gauss == 0 ){ + std::cout << "IdealReferenceForceEquilibrium: add UpdaterForceBalancedPosition(myIngredients2, threshold) \n"; + taskmanager2.addUpdater( updater ); + }else if( gauss == 1 ){ + std::cout << "IdealReferenceForceEquilibrium: add UpdaterForceBalancedPosition(myIngredients2, threshold) \n"; + taskmanager2.addUpdater( updater2 ); + } + + taskmanager2.addAnalyzer(new AnalyzerEquilbratedPosition(myIngredients2,outputDataPos, outputDataDist)); + //initialize and run + taskmanager2.initialize(); + taskmanager2.run(); + if(gauss > 1) + taskmanager2.run(1); + taskmanager2.cleanup(); + + } + catch(std::exception& e){ + std::cerr<<"Error:\n" + <. + +--------------------------------------------------------------------------------*/ + +/****************************************************************************** + * based on LeMonADE: https://github.com/LeMonADE-project/LeMonADE/ + * author: Toni Müller + * email: mueller-toni@ipfdd.de + * project: LeMonADE-Phantom Modulus + *****************************************************************************/ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + + +double prob_q(double n, double N, double m) { + return (m+1.)/(N-n-1.); + } +int main(int argc, char* argv[]){ + try{ + /////////////////////////////////////////////////////////////////////////////// + ///parse options/// + // std::string inputBFM("init.bfm"); + std::string outputDataPos("CrosslinkPosition.dat"); + std::string outputDataDist("ChainExtensionDistribution.dat"); + std::string feCurve(""); + double relaxationParameter(10.); + double threshold(0.5); + double factor(0.995); + double stretching_factor(1.0); + uint32_t gauss(0); + + uint32_t nSegments(16); + uint32_t functionality(4); + uint32_t nRings(0); + + bool showHelp = false; + auto parser + // = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() + = clara::detail::Opt( outputDataPos, "outputDataPos (=CrosslinkPosition.dat)" ) ["-o"]["--outputPos" ] ("(optional) Output filename of the crosslink ID and the equilibrium Position.").optional() + | clara::detail::Opt( outputDataDist, "outputDataDist (=ChainExtensionDistribution.dat)") ["-c"]["--outputDist" ] ("(optional) Output filename of the chain extension distribution." ).optional() + | clara::detail::Opt( threshold, "threshold" ) ["-t"]["--threshold" ] ("(optional) Threshold of the average shift. Default 0.5 ." ).optional() + | clara::detail::Opt( stretching_factor, "stretching_factor (=1)" ) ["-l"]["--stretching_factor"] ("(optional) Stretching factor for uniaxial deformation. Default 1.0 ." ).optional() + | clara::detail::Opt( feCurve, "feCurve (="")" ) ["-f"]["--feCurve" ] ("(optional) Force-Extension curve. Default \"\"." ).required() + | clara::detail::Opt( relaxationParameter, "relaxationParameter (=10)" ) ["-r"]["--relax" ] ("(optional) Relaxation parameter. Default 10.0 ." ).optional() + | clara::detail::Opt( gauss, "gauss" ) ["-g"]["--gauss" ] ("(optional) Deforma with a Gaussian deformation behaviour. Default 1.0 ." ).optional() + | clara::detail::Opt( nSegments, "nSegments" ) ["-n"]["--nSegments" ] ("(optional) Number of segments for the strand." ).optional() + | clara::detail::Opt( functionality, "nStrands" ) ["-s"]["--nStrands" ] ("(optional) Functionality." ).optional() + | clara::detail::Opt( nRings, "nRings" ) ["-m"]["--nRings" ] ("(optional) number of rings." ).optional() + | clara::Help( showHelp ); + + auto result = parser.parse( clara::Args( argc, argv ) ); + + if( !result ) { + std::cerr << "Error in command line: " << result.errorMessage() << std::endl; + exit(1); + }else if(showHelp == true){ + std::cout << "Force equilibration of a star, where each chain is replaced effectively by a tendomer."<< std::endl; + parser.writeToStream(std::cout); + exit(0); + }else{ + std::cout << "outputData : " << outputDataPos << std::endl; + std::cout << "outputDataDist : " << outputDataDist << std::endl; + std::cout << "threshold : " << threshold << std::endl; + std::cout << "feCurve : " << feCurve << std::endl; + std::cout << "gauss : " << gauss << std::endl; + std::cout << "stretching_factor : " << stretching_factor << std::endl; + std::cout << "nSegments : " << nSegments << std::endl; + std::cout << "functionality : " << functionality << std::endl; + std::cout << "nRings : " << nRings << std::endl; + } + RandomNumberGenerators rng; + // rng.seedDefaultValuesAll(); + rng.seedAll(); + /////////////////////////////////////////////////////////////////////////////// + ///end options parsing + /////////////////////////////////////////////////////////////////////////////// + //Read in th last Config + typedef LOKI_TYPELIST_1(FeatureMoleculesIO) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + myIngredients.setBoxX(256); + myIngredients.setBoxY(256); + myIngredients.setBoxZ(256); + myIngredients.setPeriodicX(1); + myIngredients.setPeriodicY(1); + myIngredients.setPeriodicZ(1); + myIngredients.modifyBondset().addBFMclassicBondset(); + myIngredients.synchronize(); + TaskManager taskmanager; + taskmanager.addUpdater( new UpdaterAddStars(myIngredients,1, 2*nSegments+1 , functionality ),0); + taskmanager.addAnalyzer(new AnalyzerWriteBfmFile("config.bfm", myIngredients, AnalyzerWriteBfmFile::APPEND) ); + taskmanager.initialize(); + taskmanager.run(1); + taskmanager.cleanup(); + + std::cout << "Read in conformation and go on to bring it into equilibrium forces..." < Config2; + typedef Ingredients Ing2; + Ing2 myIngredients2; + + myIngredients2.setBoxX(myIngredients.getBoxX()); + myIngredients2.setBoxY(myIngredients.getBoxY()); + myIngredients2.setBoxZ(myIngredients.getBoxZ()); + myIngredients2.setPeriodicX(myIngredients.isPeriodicX()); + myIngredients2.setPeriodicY(myIngredients.isPeriodicY()); + myIngredients2.setPeriodicZ(myIngredients.isPeriodicZ()); + myIngredients2.modifyMolecules().resize(myIngredients.getMolecules().size()); + myIngredients2.modifyMolecules().setAge(myIngredients.getMolecules().getAge()); + + myIngredients2.setNumOfChains (functionality*1); + myIngredients2.setNumOfCrosslinks (functionality+1); + myIngredients2.setNumOfMonomersPerChain (2*nSegments); + myIngredients2.setNumOfMonomersPerCrosslink(1); + myIngredients2.setFunctionality (functionality); + + + for(size_t i = 0; i< myIngredients.getMolecules().size();i++){ + myIngredients2.modifyMolecules()[i].modifyVector3D()=myIngredients.getMolecules()[i].getVector3D(); + for (size_t j = 0 ; j < myIngredients.getMolecules().getNumLinks(i);j++){ + uint32_t neighbor(myIngredients.getMolecules().getNeighborIdx(i,j)); + if( ! myIngredients2.getMolecules().areConnected(i,neighbor) ) + myIngredients2.modifyMolecules().connect(i,neighbor); + } + } + + //probability distribution function + std::cout << "Calculate PDF" << std::endl; + std::vector PF((nSegments-nRings),0); + double sum(0.); + for (auto i = 1; i < (nSegments- nRings); i++ ){ + PF[i]=prob_q(i,nSegments,nRings)*(1.-sum); + // std::cout << "probabilityDisti "< convPF((nSegments- nRings)*2,0); + for (auto i=0; i < 2*(nSegments-nRings);i++){ + for(auto j=0; j < i ; j++){ + convPF[i]+=PF[j]*PF[i-j]; + } + // std::cout << "ConvprobabilityDisti "< CPF((nSegments- nRings)*2,0); + for (auto i=1; i < 2*(nSegments-nRings);i++){ + CPF[i]+=CPF[i-1]+convPF[i]; + } + //calculate the inverse cummulative convoluted probability distribution + std::cout << "Calculate inverse cummulative convoluted PDF" << std::endl; + uint32_t steps(50000); + std::vector invCPF(steps-1,0); + for (auto i =1 ; i < steps-1; i++) { + uint32_t n=0 ; + auto prob=static_cast (i) / static_cast(steps ); + while ( prob > CPF [ n ] ){ n++;} + invCPF[i]=(n-1) + static_cast(round( (n-(n-1)) *( prob- CPF[n-1] )/(CPF[n]-CPF[n-1]) )); + // std::cout << "invCPF: " << prob << " " << invCPF[i] <<" "<< n-1 <<" "<(myIngredients2, threshold) \n"; + taskmanager2.addUpdater( updater ); + }else if( gauss == 1 ){ + std::cout << "IdealReferenceForceEquilibrium: add UpdaterForceBalancedPosition(myIngredients2, threshold) \n"; + taskmanager2.addUpdater( updater2 ); + } + + taskmanager2.addAnalyzer(new AnalyzerEquilbratedPosition(myIngredients2,outputDataPos, outputDataDist)); + //initialize and run + taskmanager2.initialize(); + taskmanager2.run(); + if(gauss > 1) + taskmanager2.run(1); + taskmanager2.cleanup(); + + } + catch(std::exception& e){ + std::cerr<<"Error:\n" + <. + +--------------------------------------------------------------------------------*/ + +/****************************************************************************** + * based on LeMonADE: https://github.com/LeMonADE-project/LeMonADE/ + * author: Toni Müller + * email: mueller-toni@ipfdd.de + * project: LeMonADE-Phantom Modulus + *****************************************************************************/ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + + +int main(int argc, char* argv[]){ + try{ + /////////////////////////////////////////////////////////////////////////////// + ///parse options/// + std::string inputBFM("init.bfm"); + std::string outputBFM("Config.bfm.dat"); + std::string inputConnection("BondCreationBreaking.dat"); + bool showHelp = false; + auto parser + = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() + | clara::detail::Opt( inputConnection, "inputConnection (=BondCreationBreaking.dat)" ) ["-d"]["--inputConnection" ] ("used for the time development of the topology. " ).required() + | clara::detail::Opt( outputBFM, "outputBFM (=Config.bfm)" ) ["-o"]["--outputBFM" ] ("Output filename for the active material of the tendomer network.") + | clara::Help( showHelp ); + + auto result = parser.parse( clara::Args( argc, argv ) ); + + if( !result ) { + std::cerr << "Error in command line: " << result.errorMessage() << std::endl; + exit(1); + }else if(showHelp == true){ + std::cout << "Analyzer taking a bond table and a table specifying the activity of an object and writes the active part of the tendoemr network."<< std::endl; + std::cout << "Important: input file for the connections and the active material must fit to each other!!!"<< std::endl; + parser.writeToStream(std::cout); + exit(0); + }else{ + std::cout << "outputBFM : " << outputBFM << std::endl; + std::cout << "inputBFM : " << inputBFM << std::endl; + std::cout << "inputConnection : " << inputConnection << std::endl; + } + /////////////////////////////////////////////////////////////////////////////// + ///end options parsing + /////////////////////////////////////////////////////////////////////////////// + typedef LOKI_TYPELIST_4(FeatureMoleculesIOUnsaveCheck, FeatureReactiveBonds, FeatureAttributes<>,FeatureSystemInformationLinearMeltWithCrosslinker) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + + TaskManager taskmanager; + + taskmanager.addUpdater( new UpdaterReadBfmFile(inputBFM,myIngredients, UpdaterReadBfmFile::READ_LAST_CONFIG_SAVE),0); + taskmanager.addUpdater( new UpdaterReadCrosslinkConnections(myIngredients, inputConnection, 0., 1.0) ); + taskmanager.addAnalyzer( new AnalyzerWriteBfmFile(outputBFM, myIngredients, 1) ); + //initialize and run + taskmanager.initialize(); + taskmanager.run(1); + taskmanager.cleanup(); + std::cout << "Read in conformation and go on to bring it into equilibrium forces..." <. + +--------------------------------------------------------------------------------*/ + +/****************************************************************************** + * based on LeMonADE: https://github.com/LeMonADE-project/LeMonADE/ + * author: Toni Müller + * email: mueller-toni@ipfdd.de + * project: LeMonADE-Phantom Modulus + *****************************************************************************/ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + + +int main(int argc, char* argv[]){ + try{ + /////////////////////////////////////////////////////////////////////////////// + ///parse options/// + std::string inputBFM("init.bfm"); + std::string outputBFM("Config.bfm.dat"); + std::string inputConnection("BondCreationBreaking.dat"); + std::string activeComponent("active_00001.dat"); + bool showHelp = false; + auto parser + = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() + | clara::detail::Opt( inputConnection, "inputConnection (=BondCreationBreaking.dat)" ) ["-d"]["--inputConnection" ] ("used for the time development of the topology. " ).required() + | clara::detail::Opt( activeComponent, "activeComponents (=active_00001.dat)" ) ["-a"]["--activeComponents"] ("sets the active components . " ).required() + | clara::detail::Opt( outputBFM, "outputBFM (=Config.bfm)" ) ["-o"]["--outputBFM" ] ("Output filename for the active material of the tendomer network.") + | clara::Help( showHelp ); + + auto result = parser.parse( clara::Args( argc, argv ) ); + + if( !result ) { + std::cerr << "Error in command line: " << result.errorMessage() << std::endl; + exit(1); + }else if(showHelp == true){ + std::cout << "Analyzer taking a bond table and a table specifying the activity of an object and writes the active part of the tendoemr network."<< std::endl; + std::cout << "Important: input file for the connections and the active material must fit to each other!!!"<< std::endl; + parser.writeToStream(std::cout); + exit(0); + }else{ + std::cout << "outputBFM : " << outputBFM << std::endl; + std::cout << "inputBFM : " << inputBFM << std::endl; + std::cout << "inputConnection : " << inputConnection << std::endl; + std::cout << "activeComponent : " << activeComponent << std::endl; + } + RandomNumberGenerators rng; + rng.seedAll(); + /////////////////////////////////////////////////////////////////////////////// + ///end options parsing + /////////////////////////////////////////////////////////////////////////////// + //Read in th last Config + typedef LOKI_TYPELIST_4(FeatureMoleculesIOUnsaveCheck, FeatureLabel, FeatureReactiveBonds, FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + + TaskManager taskmanager; + + taskmanager.addUpdater( new UpdaterReadBfmFile(inputBFM,myIngredients, UpdaterReadBfmFile::READ_LAST_CONFIG_SAVE),0); + taskmanager.addUpdater( new UpdaterReadCrosslinkConnectionsTendomer(myIngredients, inputConnection, 0., 1.0) ); + + //initialize and run + taskmanager.initialize(); + taskmanager.run(1); + taskmanager.cleanup(); + std::cout << "Read in conformation and go on to bring it into equilibrium forces..." <>activeTag; + auto tmpAttribute(0); + if(objectID < nCrossLinks ){ + auto IdX(nChainMonomers+objectID); + #ifdef DEBUG + std::cout << objectID << " " << activeTag << " " << IdX << " " << IdX << std::endl; + #endif + if (activeTag == 2 ){ + tmpAttribute=1; + nActiveCrossLinks++; + } + myIngredients.modifyMolecules()[nChainMonomers+objectID].setAttributeTag(tmpAttribute); + }else{ + auto IdCStart((objectID-nCrossLinks)*nSegments); + auto IdCEnd( (objectID-nCrossLinks+1)*nSegments-1); + #ifdef DEBUG + std::cout << objectID << " " << activeTag <<" " << IdCStart << " " << IdCEnd << std::endl; + #endif + if (activeTag == 2 ) { + nActiveTendomers++; + tmpAttribute=1; + } + for (auto i =IdCStart ; i <= IdCEnd; i++ ) + myIngredients.modifyMolecules()[i].setAttributeTag(tmpAttribute); + } + objectID++; + } + } + std::cout <<"Fileend:\n" ; + // for() + myIngredients.setNumTendomers (nActiveTendomers); + myIngredients.setNumCrossLinkers (nActiveCrossLinks); + myIngredients.setNumMonomersPerChain (myIngredients.getNumMonomersPerChain()); + myIngredients.setNumLabelsPerTendomerArm(myIngredients.getNumLabelsPerTendomerArm()); + myIngredients.synchronize(); + + + MonomerGroup subgroup(myIngredients.getMolecules()); + subgroup.clear(); + hasThisType<1> predicate; + for(size_t n=0;n(outputBFM, myIngredients, 1) ); + //initialize and run + if ( nActiveTendomers > 0 ){ + taskmanager2.initialize(); + taskmanager2.run(1); + taskmanager2.cleanup(); + } + // AnalyzerWriteBfmFileSubGroup > writer(outputBFM, myIngredients, 1, hasThisType<1>()); + // writer.initialize(); + + // TaskManager taskmanager2; + // taskmanager2.addAnalyzer( new AnalyzerWriteBfmFileSubGroup >(outputBFM, myIngredients, 1, hasThisType<1>()) ); + + // //initialize and run + // taskmanager2.initialize(); + // taskmanager2.run(1); + // // taskmanager2.cleanup(); + + } + catch(std::exception& e){ + std::cerr<<"Error:\n" + <. + +--------------------------------------------------------------------------------*/ + +/****************************************************************************** + * based on LeMonADE: https://github.com/LeMonADE-project/LeMonADE/ + * author: Toni Müller + * email: mueller-toni@ipfdd.de + * project: LeMonADE-Phantom Modulus + *****************************************************************************/ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + + +int main(int argc, char* argv[]){ + try{ + /////////////////////////////////////////////////////////////////////////////// + ///parse options/// + std::string inputBFM("init.bfm"); + std::string outputBFM("Config.bfm.dat"); + std::string inputConnection("BondCreationBreaking.dat"); + std::string activeComponent("active_00001.dat"); + bool showHelp = false; + auto parser + = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() + | clara::detail::Opt( inputConnection, "inputConnection (=BondCreationBreaking.dat)" ) ["-d"]["--inputConnection" ] ("used for the time development of the topology. " ).required() + | clara::detail::Opt( activeComponent, "activeComponents (=active_00001.dat)" ) ["-a"]["--activeComponents"] ("sets the active components . " ).required() + | clara::detail::Opt( outputBFM, "outputBFM (=Config.bfm)" ) ["-o"]["--outputBFM" ] ("Output filename for the active material of the tendomer network.") + | clara::Help( showHelp ); + + auto result = parser.parse( clara::Args( argc, argv ) ); + + if( !result ) { + std::cerr << "Error in command line: " << result.errorMessage() << std::endl; + exit(1); + }else if(showHelp == true){ + std::cout << "Analyzer taking a bond table and a table specifying the activity of an object and writes the active part of the tendoemr network."<< std::endl; + std::cout << "Important: input file for the connections and the active material must fit to each other!!!"<< std::endl; + parser.writeToStream(std::cout); + exit(0); + }else{ + std::cout << "outputBFM : " << outputBFM << std::endl; + std::cout << "inputBFM : " << inputBFM << std::endl; + std::cout << "inputConnection : " << inputConnection << std::endl; + std::cout << "activeComponent : " << activeComponent << std::endl; + } + RandomNumberGenerators rng; + rng.seedAll(); + /////////////////////////////////////////////////////////////////////////////// + ///end options parsing + /////////////////////////////////////////////////////////////////////////////// + //Read in th last Config + typedef LOKI_TYPELIST_4(FeatureMoleculesIOUnsaveCheck, FeatureLabel, FeatureReactiveBonds, FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + + TaskManager taskmanager; + + taskmanager.addUpdater( new UpdaterReadBfmFile(inputBFM,myIngredients, UpdaterReadBfmFile::READ_LAST_CONFIG_SAVE),0); + taskmanager.addUpdater( new UpdaterReadCrosslinkConnectionsTendomer(myIngredients, inputConnection, 0., 1.0) ); + + //initialize and run + taskmanager.initialize(); + taskmanager.run(1); + taskmanager.cleanup(); + std::cout << "Read in conformation and go on to bring it into equilibrium forces..." <>activeTag; + auto tmpAttribute(0); + if(objectID < nCrossLinks ){ + auto IdX(nChainMonomers+objectID); + #ifdef DEBUG + std::cout << objectID << " " << activeTag << " " << IdX << " " << IdX << std::endl; + #endif + if (activeTag == 2 || activeTag == 1 ){ + tmpAttribute=1; + nActiveCrossLinks++; + } + myIngredients.modifyMolecules()[nChainMonomers+objectID].setAttributeTag(tmpAttribute); + }else{ + auto IdCStart((objectID-nCrossLinks)*nSegments); + auto IdCEnd( (objectID-nCrossLinks+1)*nSegments-1); + #ifdef DEBUG + std::cout << objectID << " " << activeTag <<" " << IdCStart << " " << IdCEnd << std::endl; + #endif + if (activeTag == 2 || activeTag == 1 ) { + nActiveTendomers++; + tmpAttribute=1; + } + for (auto i =IdCStart ; i <= IdCEnd; i++ ) + myIngredients.modifyMolecules()[i].setAttributeTag(tmpAttribute); + } + objectID++; + } + } + std::cout <<"Fileend:\n" ; + // for() + myIngredients.setNumTendomers (nActiveTendomers); + myIngredients.setNumCrossLinkers (nActiveCrossLinks); + myIngredients.setNumMonomersPerChain (myIngredients.getNumMonomersPerChain()); + myIngredients.setNumLabelsPerTendomerArm(myIngredients.getNumLabelsPerTendomerArm()); + myIngredients.synchronize(); + + + MonomerGroup subgroup(myIngredients.getMolecules()); + subgroup.clear(); + hasThisType<1> predicate; + for(size_t n=0;n(outputBFM, myIngredients, 1) ); + //initialize and run + if ( nActiveTendomers > 0 ){ + taskmanager2.initialize(); + taskmanager2.run(1); + taskmanager2.cleanup(); + } + // AnalyzerWriteBfmFileSubGroup > writer(outputBFM, myIngredients, 1, hasThisType<1>()); + // writer.initialize(); + + // TaskManager taskmanager2; + // taskmanager2.addAnalyzer( new AnalyzerWriteBfmFileSubGroup >(outputBFM, myIngredients, 1, hasThisType<1>()) ); + + // //initialize and run + // taskmanager2.initialize(); + // taskmanager2.run(1); + // // taskmanager2.cleanup(); + + } + catch(std::exception& e){ + std::cerr<<"Error:\n" + <. + +--------------------------------------------------------------------------------*/ + +/****************************************************************************** + * based on LeMonADE: https://github.com/LeMonADE-project/LeMonADE/ + * author: Toni Müller + * email: mueller-toni@ipfdd.de + * project: LeMonADE-Phantom Modulus + *****************************************************************************/ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + + +int main(int argc, char* argv[]){ + try{ + /////////////////////////////////////////////////////////////////////////////// + ///parse options/// + std::string inputBFM("init.bfm"); + std::string outputDataPos("CrosslinkPosition.dat"); + std::string outputDataDist("ChainExtensionDistribution.dat"); + std::string feCurve(""); + double relaxationParameter(10.); + double threshold(0.5); + double factor(0.995); + double stretching_factor(1.0); + uint32_t gauss(false); + double dampingfactor(1.0); + double prestrainFactorX(1.0); + double prestrainFactorY(1.0); + double prestrainFactorZ(1.0); + + bool showHelp = false; + auto parser + = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() + | clara::detail::Opt( outputDataPos, "outputDataPos (=CrosslinkPosition.dat)" ) ["-o"]["--outputPos" ] ("(optional) Output filename of the crosslink ID and the equilibrium Position.").optional() + | clara::detail::Opt( outputDataDist, "outputDataDist (=ChainExtensionDistribution.dat)") ["-c"]["--outputDist" ] ("(optional) Output filename of the chain extension distribution." ).optional() + | clara::detail::Opt( threshold, "threshold" ) ["-t"]["--threshold" ] ("(optional) Threshold of the average shift. Default 0.5 ." ).optional() + | clara::detail::Opt( stretching_factor, "stretching_factor (=1)" ) ["-l"]["--stretching_factor"] ("(optional) Stretching factor for uniaxial deformation. Default 1.0 ." ).optional() + | clara::detail::Opt( feCurve, "feCurve (="")" ) ["-f"]["--feCurve" ] ("(optional) Force-Extension curve. Default \"\"." ).required() + | clara::detail::Opt( relaxationParameter, "relaxationParameter (=10)" ) ["-r"]["--relax" ] ("(optional) Relaxation parameter. Default 10.0 ." ).optional() + | clara::detail::Opt( gauss, "gauss" ) ["-g"]["--gauss" ] ("(optional) Deforma with a Gaussian deformation behaviour. Default 1.0 ." ).optional() + | clara::detail::Opt( dampingfactor, "damping (=1)" ) ["-d"]["--damping" ] ("(optional) Damping factor after 1E3MCS. Default 1.0." ).optional() + | clara::detail::Opt( prestrainFactorX, "prestrainFactorX (=1)" ) ["-x"]["--prestrainFactorX" ] ("(optional) Prestrain factor in X. Default 1.0." ).optional() + | clara::detail::Opt( prestrainFactorY, "prestrainFactorY (=1)" ) ["-y"]["--prestrainFactorY" ] ("(optional) Prestrain factor in Y. Default 1.0." ).optional() + | clara::detail::Opt( prestrainFactorZ, "prestrainFactorZ (=1)" ) ["-z"]["--prestrainFactorZ" ] ("(optional) Prestrain factor in Z. Default 1.0." ).optional() + | clara::Help( showHelp ); + + auto result = parser.parse( clara::Args( argc, argv ) ); + + if( !result ) { + std::cerr << "Error in command line: " << result.errorMessage() << std::endl; + exit(1); + }else if(showHelp == true){ + std::cout << "Simulator to connect linear chains with single monomers of certain functionality"<< std::endl; + parser.writeToStream(std::cout); + exit(0); + }else{ + std::cout << "outputData : " << outputDataPos << std::endl; + std::cout << "outputDataDist : " << outputDataDist << std::endl; + std::cout << "inputBFM : " << inputBFM << std::endl; + std::cout << "threshold : " << threshold << std::endl; + std::cout << "dampingfactor : " << dampingfactor << std::endl; + std::cout << "feCurve : " << feCurve << std::endl; + std::cout << "gauss : " << gauss << std::endl; + std::cout << "stretching_factor : " << stretching_factor << std::endl; + std::cout << "prestrainFactorX : " << prestrainFactorX << std::endl; + std::cout << "prestrainFactorY : " << prestrainFactorY << std::endl; + std::cout << "prestrainFactorZ : " << prestrainFactorZ << std::endl; + } + RandomNumberGenerators rng; + // rng.seedDefaultValuesAll(); + rng.seedAll(); + /////////////////////////////////////////////////////////////////////////////// + ///end options parsing + /////////////////////////////////////////////////////////////////////////////// + //Read in th last Config + typedef LOKI_TYPELIST_3(FeatureMoleculesIOUnsaveCheck, FeatureLabel, FeatureReactiveBonds) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + + TaskManager taskmanager; + + taskmanager.addUpdater( new UpdaterReadBfmFile(inputBFM,myIngredients, UpdaterReadBfmFile::READ_LAST_CONFIG_SAVE),0); + + //initialize and run + taskmanager.initialize(); + taskmanager.run(1); + taskmanager.cleanup(); + std::cout << "Read in conformation and go on to bring it into equilibrium forces..." < Config2; + typedef Ingredients Ing2; + Ing2 myIngredients2; + + myIngredients2.setBoxX(myIngredients.getBoxX()); + myIngredients2.setBoxY(myIngredients.getBoxY()); + myIngredients2.setBoxZ(myIngredients.getBoxZ()); + myIngredients2.setPeriodicX(myIngredients.isPeriodicX()); + myIngredients2.setPeriodicY(myIngredients.isPeriodicY()); + myIngredients2.setPeriodicZ(myIngredients.isPeriodicZ()); + myIngredients2.modifyMolecules().resize(myIngredients.getMolecules().size()); + myIngredients2.modifyMolecules().setAge(myIngredients.getMolecules().getAge()); + + myIngredients2.setNumTendomers (myIngredients.getNumTendomers()); + myIngredients2.setNumCrossLinkers (myIngredients.getNumCrossLinkers()); + myIngredients2.setNumMonomersPerChain (myIngredients.getNumMonomersPerChain()); + myIngredients2.setNumLabelsPerTendomerArm(myIngredients.getNumLabelsPerTendomerArm()); + + for(size_t i = 0; i< myIngredients.getMolecules().size();i++){ + myIngredients2.modifyMolecules()[i].modifyVector3D()=myIngredients.getMolecules()[i].getVector3D(); + myIngredients2.modifyMolecules()[i].setReactive(myIngredients.getMolecules()[i].isReactive()); + myIngredients2.modifyMolecules()[i].setNumMaxLinks(myIngredients.getMolecules()[i].getNumMaxLinks()); + for (size_t j = 0 ; j < myIngredients.getMolecules().getNumLinks(i);j++){ + uint32_t neighbor(myIngredients.getMolecules().getNeighborIdx(i,j)); + if( ! myIngredients2.getMolecules().areConnected(i,neighbor) ) + myIngredients2.modifyMolecules().connect(i,neighbor); + } + } + myIngredients2.synchronize(); + + TaskManager taskmanager2; + taskmanager2.addUpdater( new UpdaterAffineDeformation(myIngredients2, stretching_factor,prestrainFactorX,prestrainFactorY,prestrainFactorZ),0 ); + //read bonds and positions stepwise + auto updater = new UpdaterForceBalancedPosition(myIngredients2, threshold, dampingfactor) ; + auto updater2 = new UpdaterForceBalancedPosition(myIngredients2, threshold,dampingfactor) ; + if ( gauss == 0 ){ + updater->setFilename(feCurve); + updater->setRelaxationParameter(relaxationParameter); + std::cout << "TendomerNetworkForceEquilibrium: add UpdaterForceBalancedPosition(myIngredients2, threshold) \n"; + taskmanager2.addUpdater( updater ); + }else if (gauss == 1 ){ + std::cout << "TendomerNetworkForceEquilibrium: add UpdaterForceBalancedPosition(myIngredients2, threshold) \n"; + taskmanager2.addUpdater( updater2 ); + } + taskmanager2.addAnalyzer(new AnalyzerEquilbratedPosition(myIngredients2,outputDataPos, outputDataDist)); + //initialize and run + taskmanager2.initialize(); + taskmanager2.run(); + taskmanager2.cleanup(); + + } + catch(std::exception& e){ + std::cerr<<"Error:\n" + <. + +--------------------------------------------------------------------------------*/ + +/****************************************************************************** + * based on LeMonADE: https://github.com/LeMonADE-project/LeMonADE/ + * author: Toni Müller + * email: mueller-toni@ipfdd.de + * project: LeMonADE-Phantom Modulus + *****************************************************************************/ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +#include +#include +#include +#include +#include + + +int main(int argc, char* argv[]){ + try{ + /////////////////////////////////////////////////////////////////////////////// + ///parse options/// + std::string inputBFM("init.bfm"); + std::string outputBFM("Config.bfm.dat"); + std::string inputConnection("BondCreationBreaking.dat"); + bool showHelp = false; + auto parser + = clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required() + | clara::detail::Opt( inputConnection, "inputConnection (=BondCreationBreaking.dat)" ) ["-d"]["--inputConnection" ] ("used for the time development of the topology. " ).required() + | clara::detail::Opt( outputBFM, "outputBFM (=Config.bfm)" ) ["-o"]["--outputBFM" ] ("Output filename for the active material of the tendomer network.") + | clara::Help( showHelp ); + + auto result = parser.parse( clara::Args( argc, argv ) ); + + if( !result ) { + std::cerr << "Error in command line: " << result.errorMessage() << std::endl; + exit(1); + }else if(showHelp == true){ + std::cout << "Analyzer taking a bond table and a table specifying the activity of an object and writes the active part of the tendoemr network."<< std::endl; + std::cout << "Important: input file for the connections and the active material must fit to each other!!!"<< std::endl; + parser.writeToStream(std::cout); + exit(0); + }else{ + std::cout << "outputBFM : " << outputBFM << std::endl; + std::cout << "inputBFM : " << inputBFM << std::endl; + std::cout << "inputConnection : " << inputConnection << std::endl; + } + RandomNumberGenerators rng; + rng.seedAll(); + /////////////////////////////////////////////////////////////////////////////// + ///end options parsing + /////////////////////////////////////////////////////////////////////////////// + //Read in th last Config + typedef LOKI_TYPELIST_4(FeatureMoleculesIOUnsaveCheck, FeatureLabel, FeatureReactiveBonds, FeatureAttributes<>) Features; + typedef ConfigureSystem Config; + typedef Ingredients Ing; + Ing myIngredients; + + TaskManager taskmanager; + + taskmanager.addUpdater( new UpdaterReadBfmFile(inputBFM,myIngredients, UpdaterReadBfmFile::READ_LAST_CONFIG_SAVE),0); + taskmanager.addUpdater( new UpdaterReadCrosslinkConnectionsTendomer(myIngredients, inputConnection, 0., 1.0) ); + taskmanager.addAnalyzer( new AnalyzerWriteBfmFile(outputBFM, myIngredients, 1) ); + //initialize and run + taskmanager.initialize(); + taskmanager.run(1); + taskmanager.cleanup(); + std::cout << "Read in conformation and go on to bring it into equilibrium forces..." < updater(ingredients, 0.0001); - updater.execute(); - auto vec2=LemonadeDistCalcs::MinImageVector(VectorDouble3(0.,0.,0.),ingredients.getMolecules()[0].getVector3D(),ingredients ); - REQUIRE(vec2.getX() == Approx(5.3125)); - REQUIRE(vec2.getY() == Approx(6.375)); - REQUIRE(vec2.getZ() == Approx(5.8125)); + // UpdaterForceBalancedPosition updater(ingredients, 0.0001); + // updater.execute(); + // auto vec2=LemonadeDistCalcs::MinImageVector(VectorDouble3(0.,0.,0.),ingredients.getMolecules()[0].getVector3D(),ingredients ); + // REQUIRE(vec2.getX() == Approx(5.3125)); + // REQUIRE(vec2.getY() == Approx(6.375)); + // REQUIRE(vec2.getZ() == Approx(5.8125)); } //restore cout std::cout.rdbuf(originalBuffer);