Skip to content

Commit

Permalink
Merge pull request #9 from tonimueller/develop
Browse files Browse the repository at this point in the history
Develop
  • Loading branch information
tonimueller authored Jan 31, 2022
2 parents 62a24ad + 10fbe21 commit c739341
Show file tree
Hide file tree
Showing 28 changed files with 3,751 additions and 169 deletions.
20 changes: 17 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
42 changes: 24 additions & 18 deletions include/LeMonADE_PM/analyzer/AnalyzerEquilbratedPosition.h
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,13 @@ template < class IngredientsType > class AnalyzerEquilbratedPosition : public Ab
std::vector< std::vector<double> > CalculateDistance();

//! just collects the id and the position for the cross links
std::vector<std::vector<int> > CollectAveragePositions();
std::vector<std::vector<double> > CollectAveragePositions();

//! setter for the filenames
void setFilenames(std::string outAvPosBasename_, std::string outDistBasename_ ){
outAvPosBasename= outAvPosBasename_;
outDistBasename=outDistBasename_;
}
};

/*************************************************************************
Expand All @@ -102,18 +108,18 @@ std::vector< std::vector<double> > AnalyzerEquilbratedPosition<IngredientsType>
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=("<<ingredients.getMolecules()[IDx].getVector3D() <<") \n"
<< "position2=("<<ingredients.getMolecules()[neighbors[j].ID].getVector3D() <<") \n"
<< std::endl;
throw std::runtime_error(error_message.str());
}
// 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=("<<ingredients.getMolecules()[IDx].getVector3D() <<") \n"
// << "position2=("<<ingredients.getMolecules()[neighbors[j].ID].getVector3D() <<") \n"
// << std::endl;
// throw std::runtime_error(error_message.str());
// }
dist[0].push_back(IDx);
dist[1].push_back(neighbors[j].ID);
dist[2].push_back(vec.getX());
Expand All @@ -128,8 +134,8 @@ std::vector< std::vector<double> > AnalyzerEquilbratedPosition<IngredientsType>
}
////////////////////////////////////////////////////////////////////////////////
template< class IngredientsType >
std::vector< std::vector<int> > AnalyzerEquilbratedPosition<IngredientsType>::CollectAveragePositions(){
std::vector<std::vector<int> > AveragePosition(4, std::vector<int>());
std::vector< std::vector<double> > AnalyzerEquilbratedPosition<IngredientsType>::CollectAveragePositions(){
std::vector<std::vector<double> > AveragePosition(4, std::vector<double>());
auto crosslinkID(ingredients.getCrosslinkIDs());
for (size_t i = 0 ; i < crosslinkID.size(); i++){
auto IDx(crosslinkID[i]);
Expand Down Expand Up @@ -189,14 +195,14 @@ void AnalyzerEquilbratedPosition<IngredientsType>::dumpData()
std::cout << "////////////////////////////////////"<<std::endl;

//output for the equilibrated positions
std::vector< std::vector<int> > CrossLinkPositions=CollectAveragePositions() ;
std::vector< std::vector<double> > CrossLinkPositions=CollectAveragePositions() ;
std::stringstream commentAveragePosition;
commentAveragePosition<<"Created by AnalyzerEquilbratedPosition\n";
commentAveragePosition<<"ID's start at 0 \n";
commentAveragePosition<<"conversion="<<conversion<<"\n";
commentAveragePosition<<"ID equilibrated position\n";
std::stringstream outAvPos;
outAvPos<< std::setprecision(2)<< "C" << conversion;
outAvPos<< std::setprecision(3)<< "C" << conversion;
outAvPos << "_" << outAvPosBasename;


Expand All @@ -217,7 +223,7 @@ void AnalyzerEquilbratedPosition<IngredientsType>::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(
Expand Down
129 changes: 129 additions & 0 deletions include/LeMonADE_PM/analyzer/AnalyzerMolecularWeight.h
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
--------------------------------------------------------------------------------*/

#ifndef LEMONADE_PM_ANALYZER_MOLECULAR_WEIGTH_H
#define LEMONADE_PM_ANALYZER_MOLECULAR_WEIGTH_H

#include <string>

#include <LeMonADE/utility/Vector3D.h>
#include <LeMonADE/analyzer/AbstractAnalyzer.h>
#include <LeMonADE/utility/ResultFormattingTools.h>
#include <LeMonADE/utility/MonomerGroup.h>
#include <LeMonADE/utility/DepthIterator.h>
#include <LeMonADE/analyzer/AnalyzerAbstractDump.h>


/*************************************************************************
* 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<IngredientsType,double>
{
typedef AnalyzerAbstractDump<IngredientsType,double> 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<class IngredientsType>
AnalyzerMolecularWeight<IngredientsType>::AnalyzerMolecularWeight(
const IngredientsType& ing,std::string fileSuffix_)
:BaseClass(ing,fileSuffix_) {}

template <class IngredientsType>
void AnalyzerMolecularWeight<IngredientsType>::initialize()
{
basename=BaseClass::getOutputFilename() ;
BaseClass::setBufferSize(1);
BaseClass::setNumberOfColumns(1);
execute();
}
template< class IngredientsType >
bool AnalyzerMolecularWeight<IngredientsType>::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<MonomerGroup<molecules_type> > groups;
fill_connected_groups(ingredients.getMolecules(), groups, MonomerGroup<molecules_type>(ingredients.getMolecules()), alwaysTrue());
//key is the size of the molecule and the value is the number of occurence
std::map<uint32_t,uint32_t> 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
12 changes: 11 additions & 1 deletion include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUp.h
Original file line number Diff line number Diff line change
Expand Up @@ -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<neighborX> 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);
};

Expand Down
Loading

0 comments on commit c739341

Please sign in to comment.