Skip to content

Commit 62a24ad

Browse files
tonimuellertonimueller
authored
Merge pull request #7 from LeMonADE-project/develop
Develop
2 parents c1ce655 + aa3bbf4 commit 62a24ad

10 files changed

Lines changed: 332 additions & 97 deletions

‎include/LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUp.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -125,7 +125,7 @@ void FeatureCrosslinkConnectionsLookUp::fillTables(IngredientsType& ingredients)
125125
VectorDouble3 jumpVector(posHead-bond-posX); // tracks if one bond jumps across periodic images
126126

127127
//direct connection of two cross links
128-
if (molecules.getNumLinks(head) > 2) {
128+
if ( molecules[head].isReactive() && molecules.getNumLinks(head) > 2) {
129129
NeighborIDs.push_back( neighborX(head, 1, jumpVector) );
130130
}else{
131131
uint32_t nSegments(1);

‎include/LeMonADE_PM/updater/UpdaterForceBalancedPosition.h‎

Lines changed: 22 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -29,24 +29,26 @@ along with LeMonADE. If not, see <http://www.gnu.org/licenses/>.
2929

3030

3131
#include <LeMonADE/updater/AbstractUpdater.h>
32-
#include <LeMonADE_PM/updater/moves/MoveForceEquilibrium.h>
3332
#include <vector>
3433
/**
3534
* @class UpdaterForceBalancedPosition
3635
* @tparam IngredientsType
3736
*/
3837

39-
template <class IngredientsType>
38+
template <class IngredientsType, class moveType >
4039
class UpdaterForceBalancedPosition:public AbstractUpdater
4140
{
4241
public:
4342
//! constructor for UpdaterForceBalancedPosition
44-
UpdaterForceBalancedPosition(IngredientsType& ing_, double threshold_):ing(ing_),threshold(threshold_){};
43+
UpdaterForceBalancedPosition(IngredientsType& ing_, double threshold_ ):
44+
ing(ing_),threshold(threshold_){};
4545

46-
virtual void initialize(){};
46+
virtual void initialize(){}// init(move);};
4747
bool execute();
4848
virtual void cleanup(){};
4949

50+
void setFilename(const std::string filename) {move.setFilename(filename); }
51+
void setRelaxationParameter( const double relax ) {move.setRelaxationParameter(relax);}
5052
private:
5153
//!copy of the main container for the system informations
5254
IngredientsType& ing;
@@ -55,13 +57,26 @@ class UpdaterForceBalancedPosition:public AbstractUpdater
5557
double threshold;
5658

5759
//! move to equilibrate the cross links by force equilibrium
58-
MoveForceEquilibrium move;
60+
moveType move;
5961

6062
//! random number generator
6163
RandomNumberGenerators rng;
64+
65+
// //! input filename for the MoveNonLinearForceEquilibrium
66+
// std::string filename;
67+
// //! relaxation parameter
68+
// double relaxationParameter;
69+
// //!initialize of the move : MoveForceEquilibrium
70+
// void init ( MoveForceEquilibrium& move ) {std::cout << "Nothing to initialize for this move type.\n";}
71+
// //!initialize of the move :MoveNonLinearForceEquilibrium
72+
// void init (MoveNonLinearForceEquilibrium& move ) {
73+
// move.setRelaxationParameter(relaxationParameter);
74+
// move.setFilename(filename);
75+
// };
76+
6277
};
63-
template <class IngredientsType>
64-
bool UpdaterForceBalancedPosition<IngredientsType>::execute(){
78+
template <class IngredientsType, class moveType>
79+
bool UpdaterForceBalancedPosition<IngredientsType,moveType>::execute(){
6580
std::cout << "UpdaterForceBalancedPosition::execute(): Start equilibration" <<std::endl;
6681
double avShift(threshold*1.1);
6782
uint32_t StartMCS(ing.getMolecules().getAge());

‎include/LeMonADE_PM/updater/moves/MoveForceEquilibrium.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -64,12 +64,12 @@ class MoveForceEquilibrium:public MoveForceEquilibriumBase<MoveForceEquilibrium>
6464
//Gaussina force extension relation
6565
//f=R*3/(N*b^2)
6666
VectorDouble3 FE(VectorDouble3 extensionVector, double nSegs){
67-
return extensionVector*3./(std::sqrt(nSegs)*bondlength*bondlength);
67+
return extensionVector*3./((nSegs)*bondlength*bondlength);
6868
}
6969
//Gaussian extension force relation
7070
//R=-f/3*N*b^2
7171
VectorDouble3 EF(VectorDouble3 force, double nSegs){
72-
return force/(3.)*std::sqrt(nSegs)*bondlength*bondlength;
72+
return force/(3.)*(nSegs)*bondlength*bondlength;
7373
}
7474

7575
//calculate the shift for the cross link

‎include/LeMonADE_PM/updater/moves/MoveNonLinearForceEquilibrium.h‎

Lines changed: 67 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ along with LeMonADE. If not, see <http://www.gnu.org/licenses/>.
2929
#define LEMONADE_PM_UPDATER_MOVES_MOVENONLINEARFORCEEQUILIBRIUM_H
3030
#include <limits>
3131
#include <fstream>
32-
#include <LeMonADE_PM/updater/moves/MoveNonLinearForceEquilibriumBase.h>
32+
#include <LeMonADE_PM/updater/moves/MoveForceEquilibriumBase.h>
3333
#include <LeMonADE/utility/DistanceCalculation.h>
3434
#include <LeMonADE_PM/utility/neighborX.h>
3535

@@ -48,8 +48,12 @@ along with LeMonADE. If not, see <http://www.gnu.org/licenses/>.
4848

4949
class MoveNonLinearForceEquilibrium:public MoveForceEquilibriumBase<MoveNonLinearForceEquilibrium>{
5050
public:
51-
MoveNonLinearForceEquilibrium(std::filename input_, double relaxationChain_):input(input_),relaxationChain(relaxationChain_),bondlength(2.68){};
52-
51+
//! constructor for the class MoveNonLinearForceEquilibrium taking the filename of the force extension relation and the realaxation parameter for the chain
52+
MoveNonLinearForceEquilibrium(std::string filename_="", double relaxationChain_=1.):
53+
filename(filename_),bondlength(2.68){
54+
if( !filename.empty() ) createTable();
55+
setRelaxationParameter(relaxationChain_);
56+
};
5357
// overload initialise function to be able to set the moves index and direction if neccessary
5458
template <class IngredientsType> void init(const IngredientsType& ing);
5559
template <class IngredientsType> void init(const IngredientsType& ing, uint32_t index);
@@ -58,50 +62,74 @@ class MoveNonLinearForceEquilibrium:public MoveForceEquilibriumBase<MoveNonLinea
5862
template <class IngredientsType> bool check(IngredientsType& ing);
5963
template< class IngredientsType> void apply(IngredientsType& ing);
6064

65+
//!read file to create a lookup table of the force extension curve
66+
void createTable();
67+
// //!read file to create a lookup table of the force extension curve
68+
// void createTable(std::string filename_){
69+
// setFilename(filename);
70+
// createTable();}
71+
//! set the filename for the force extension data
72+
void setFilename(std::string filename_){filename=filename_;createTable();}
73+
//! get the filename for the force extension data
74+
std::string const getFilename(){return filename;}
75+
76+
//! set the relaxation parameter for the cross link
77+
void setRelaxationParameter(double relaxationChain_){
78+
relaxationChain=relaxationChain_;
79+
springConstant = relaxationChain*bondlength*bondlength/(-3.);
80+
}
81+
//! get the relaxation parameter for the cross link
82+
double getRelaxationParameter(){return relaxationChain;}
83+
//uses the read force extension relation
84+
VectorDouble3 EF(VectorDouble3 extensionVector){
85+
uint32_t length( static_cast<uint32_t>(round(extensionVector.getLength())) );
86+
if ( max_extension < length ){
87+
std::stringstream errormessage;
88+
errormessage << "The length of the extension vector is greater than the maximum given in the file: \n"
89+
<< "length is " << length
90+
<< " maximum is " << max_extension << "\n";
91+
throw std::runtime_error(errormessage.str());
92+
}
93+
if (length == 0 )
94+
return VectorDouble3(0.,0.,0.);
95+
return force_extension[length]*extensionVector.normalize();
96+
97+
}
98+
//Gaussian extension force relation
99+
//R=-f/3*N*b^2
100+
VectorDouble3 FE(VectorDouble3 const force ) const {
101+
return force*springConstant;
102+
}
103+
61104
private:
62105
//average square bond length
63106
const double bondlength;
64-
//input filename for the force extension curve
65-
std::string input;
107+
//filename filename for the force extension curve
108+
std::string filename;
66109

67110
//minimum force in the file
68-
double min_force:
111+
double min_force;
69112
//maximum force in the file
70113
double max_force;
71114
//force steps
72115
double force_steps;
73-
116+
74117
//min_extension in the file
75118
double min_extension;
76119
//maximum extensions in the file
77120
double max_extension;
78121
//extension steps
79122
double extension_steps;
80-
123+
124+
//spring constant for the equivalent chain used for the relaxation of the cross links
125+
double springConstant;
126+
81127
//force extension mapping
82-
std::map<double, double> extension_force
83-
//extension force mapping
128+
std::map<double, double> extension_force;
129+
//extension force mapping index=extension rounded to int
84130
std::vector<double> force_extension;
85-
86131
//an equivalent chain which relaxes the cross link
87132
double relaxationChain;
88-
89-
90-
91-
//read file to create a lookup table of the force extension curve
92-
void createTable();
93-
94-
//Gaussina force extension relation
95-
//f=R*3/(N*b^2)
96-
VectorDouble3 FE(VectorDouble3 extensionVector, double nSegs){
97-
return extensionVector*3./(std::sqrt(nSegs)*bondlength*bondlength);
98-
}
99-
//Gaussian extension force relation
100-
//R=-f/3*N*b^2
101-
VectorDouble3 EF(VectorDouble3 force, double nSegs){
102-
return force/(-3.)*std::sqrt(nSegs)*bondlength*bondlength;
103-
}
104-
105133
//calculate the shift for the cross link
106134
template< class IngredientsType >
107135
VectorDouble3 CalculateShift(IngredientsType& ing ){
@@ -111,42 +139,35 @@ class MoveNonLinearForceEquilibrium:public MoveForceEquilibriumBase<MoveNonLinea
111139
double avNSegments(0.);
112140
if (Neighbors.size() > 0) {
113141
VectorDouble3 Position(ing.getMolecules()[this->getIndex()].getVector3D());
114-
// std::cout << "CorssLinkPos=" << Position << std::endl;
115142
for (size_t i = 0; i < Neighbors.size(); i++){
116143
VectorDouble3 vec(Position-ing.getMolecules()[Neighbors[i].ID].getVector3D()-Neighbors[i].jump);
117-
// avNSegments+=1./Neighbors[i].segDistance;
118-
// force+=FE(vec,Neighbors[i].segDistance);
119-
force+=force_extension[vec.getLength()]*vec.normalize();
144+
force+=EF(vec);//Neighbors[i].segDistance
120145
}
121-
// force/=(1.*Neighbors.size());
122-
shift=EF(force,relaxationChain);
146+
shift=FE(force/(static_cast<double>(Neighbors.size()) ));
123147
}
124148
return shift;
125149
};
126-
127-
128150
};
129151
/////////////////////////////////////////////////////////////////////////////
130152
/////////// implementation of the members ///////////////////////////////////
131153
void MoveNonLinearForceEquilibrium::createTable(){
132-
std::ifstream in(input);
154+
std::ifstream in(filename);
133155
uint32_t counter(0);
134-
while(in.good()){
156+
while(in.good() && in.peek()!=EOF){
135157
std::string line;
136-
getline(stream, line);
158+
getline(in, line);
137159
//ignore comments and blank lines
138-
if (line.empty())
139-
findStartofData = true;
140-
else if (line.at(0) == '#') //go to next line
160+
while (line.at(0) == '#' || line.empty() ) //go to next line
141161
continue;
142162
//read data
143163
double force, extension;
144-
std::stringstream ss << line;
164+
std::stringstream ss ;
165+
ss<< line;
145166
ss>>force >> extension;
146167
if(min_force > force ) min_force=force;
147168
if(max_force < force ) max_force=force;
148169
if(min_extension > extension ) min_extension=extension;
149-
if(min_extension < extension ) max_extension=extension;
170+
if(max_extension < extension ) max_extension=extension;
150171
extension_force.insert(extension_force.end(),std::pair<double,double>(force, extension));
151172
if(counter==1){
152173
force_steps=max_force-min_force;
@@ -158,9 +179,9 @@ void MoveNonLinearForceEquilibrium::createTable(){
158179
//make a entry from 0 to max_extension in steps of 1
159180
for ( auto r=0;r<static_cast<uint32_t>(max_extension); r++ ){
160181
if(r==0)
161-
extension_force.push_back(0.);
182+
force_extension.push_back(0.);
162183
else{
163-
auto it_last=extension_force.begin()
184+
auto it_last=extension_force.begin();
164185
for (auto it=extension_force.begin(); it !=extension_force.end();it++){
165186
if(it->second > r){
166187
//at the force linear interpolated in between the two current forces
@@ -173,12 +194,9 @@ void MoveNonLinearForceEquilibrium::createTable(){
173194
}
174195
it_last=it;
175196
}
176-
177197
}
178198
}
179199
}
180-
181-
182200
/*****************************************************************************/
183201
/**
184202
* @brief Initialize the move.

‎projects/ForceEquilibrium.cpp‎

Lines changed: 25 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,8 @@ along with LeMonADE. If not, see <http://www.gnu.org/licenses/>.
4848

4949
#include <LeMonADE_PM/updater/UpdaterForceBalancedPosition.h>
5050
#include <LeMonADE_PM/updater/UpdaterReadCrosslinkConnections.h>
51+
#include <LeMonADE_PM/updater/moves/MoveForceEquilibrium.h>
52+
#include <LeMonADE_PM/updater/moves/MoveNonLinearForceEquilibrium.h>
5153
#include <LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUp.h>
5254
#include <LeMonADE_PM/analyzer/AnalyzerEquilbratedPosition.h>
5355

@@ -60,19 +62,24 @@ int main(int argc, char* argv[]){
6062
std::string outputDataPos("CrosslinkPosition.dat");
6163
std::string outputDataDist("ChainExtensionDistribution.dat");
6264
std::string inputConnection("BondCreationBreaking.dat");
65+
std::string feCurve("");
66+
double relaxationParameter(10.);
6367
double threshold(0.5);
6468
double stepwidth(1.0);
6569
double minConversion(50.0);
70+
bool custom(false);
6671

6772
bool showHelp = false;
6873
auto parser
69-
= clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required()
70-
| clara::detail::Opt( inputConnection, "inputConnection (=BondCreationBreaking.dat)" ) ["-d"]["--inputConnection"] ("used for the time development of the topology. " ).required()
71-
| clara::detail::Opt( outputDataPos, "outputDataPos (=CrosslinkPosition.dat)" ) ["-o"]["--outputPos" ] ("(optional) Output filename of the crosslink ID and the equilibrium Position.").optional()
72-
| clara::detail::Opt( outputDataDist, "outputDataDist (=ChainExtensionDistribution.dat)") ["-c"]["--outputDist" ] ("(optional) Output filename of the chain extension distribution." ).optional()
73-
| clara::detail::Opt( stepwidth, "stepwidth" ) ["-s"]["--stepwidth" ] ("(optional) Width for the increase in percentage. Default: 1%." ).optional()
74-
| clara::detail::Opt( minConversion, "minConversion" ) ["-u"]["--minConversion" ] ("(optional) Minimum conversion to be read in. Default: 50%." ).optional()
75-
| clara::detail::Opt( threshold, "threshold" ) ["-t"]["--threshold" ] ("(optional) Threshold of the average shift. Default 0.5 ." ).optional()
74+
= clara::detail::Opt( inputBFM, "inputBFM (=inconfig.bfm)" ) ["-i"]["--input" ] ("(required)Input filename of the bfm file" ).required()
75+
| clara::detail::Opt( inputConnection, "inputConnection (=BondCreationBreaking.dat)" ) ["-d"]["--inputConnection"] ("used for the time development of the topology. " ).required()
76+
| clara::detail::Opt( outputDataPos, "outputDataPos (=CrosslinkPosition.dat)" ) ["-o"]["--outputPos" ] ("(optional) Output filename of the crosslink ID and the equilibrium Position.").optional()
77+
| clara::detail::Opt( outputDataDist, "outputDataDist (=ChainExtensionDistribution.dat)") ["-c"]["--outputDist" ] ("(optional) Output filename of the chain extension distribution." ).optional()
78+
| clara::detail::Opt( stepwidth, "stepwidth" ) ["-s"]["--stepwidth" ] ("(optional) Width for the increase in percentage. Default: 1%." ).optional()
79+
| clara::detail::Opt( minConversion, "minConversion" ) ["-u"]["--minConversion" ] ("(optional) Minimum conversion to be read in. Default: 50%." ).optional()
80+
| clara::detail::Opt( threshold, "threshold" ) ["-t"]["--threshold" ] ("(optional) Threshold of the average shift. Default 0.5 ." ).optional()
81+
| clara::detail::Opt( feCurve, "feCurve (="")" ) ["-f"]["--feCurve" ] ("(optional) Force-Extension curve. Default \"\"." ).optional()
82+
| clara::detail::Opt( relaxationParameter, "relaxationParameter (=10)" ) ["-r"]["--relax" ] ("(optional) Relaxation parameter. Default 10.0 ." ).optional()
7683
| clara::Help( showHelp );
7784

7885
auto result = parser.parse( clara::Args( argc, argv ) );
@@ -92,7 +99,10 @@ int main(int argc, char* argv[]){
9299
std::cout << "stepwidth : " << stepwidth << std::endl;
93100
std::cout << "minConversion : " << minConversion << std::endl;
94101
std::cout << "threshold : " << threshold << std::endl;
102+
std::cout << "feCurve : " << feCurve << std::endl;
95103
}
104+
if (! feCurve.empty()) custom=true;
105+
96106
///////////////////////////////////////////////////////////////////////////////
97107
///end options parsing
98108
///////////////////////////////////////////////////////////////////////////////
@@ -146,7 +156,14 @@ int main(int argc, char* argv[]){
146156
TaskManager taskmanager2;
147157
//read bonds and positions stepwise
148158
taskmanager2.addUpdater( new UpdaterReadCrosslinkConnections<Ing2>(myIngredients2, inputConnection, stepwidth, minConversion) );
149-
taskmanager2.addUpdater( new UpdaterForceBalancedPosition<Ing2>(myIngredients2, threshold) );
159+
if (custom)
160+
taskmanager2.addUpdater( new UpdaterForceBalancedPosition<Ing2,MoveForceEquilibrium>(myIngredients2, threshold) );
161+
else {
162+
auto updater = new UpdaterForceBalancedPosition<Ing2,MoveNonLinearForceEquilibrium>(myIngredients2, threshold) ;
163+
updater->setFilename(feCurve);
164+
updater->setRelaxationParameter(relaxationParameter);
165+
taskmanager2.addUpdater( updater );
166+
}
150167
taskmanager2.addAnalyzer(new AnalyzerEquilbratedPosition<Ing2>(myIngredients2,outputDataPos, outputDataDist));
151168
//initialize and run
152169
taskmanager2.initialize();

‎tests/feature/TestFeatureCrosslinkConnectionsLookUp.cpp‎

Lines changed: 16 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -33,14 +33,14 @@ along with LeMonADE. If not, see <http://www.gnu.org/licenses/>.
3333
#include <LeMonADE/core/Ingredients.h>
3434
#include <LeMonADE/feature/FeatureBox.h>
3535
#include <LeMonADE/utility/Vector3D.h>
36-
36+
#include <LeMonADE/feature/FeatureSystemInformationLinearMeltWithCrosslinker.h>
3737
#include <extern/catch.hpp>
3838

3939
#include <LeMonADE_PM/feature/FeatureCrosslinkConnectionsLookUp.h>
4040

4141
TEST_CASE( "Test class FeatureCrosslinkConnectionsLookUp" )
4242
{
43-
typedef LOKI_TYPELIST_2(FeatureBox, FeatureCrosslinkConnectionsLookUp ) Features;
43+
typedef LOKI_TYPELIST_3(FeatureBox, FeatureCrosslinkConnectionsLookUp,FeatureSystemInformationLinearMeltWithCrosslinker ) Features;
4444
typedef ConfigureSystem<VectorDouble3,Features,4> Config;
4545
typedef Ingredients<Config> IngredientsType;
4646

@@ -94,6 +94,20 @@ TEST_CASE( "Test class FeatureCrosslinkConnectionsLookUp" )
9494
ingredients.modifyMolecules().addMonomer(8.,6.,6.);
9595
ingredients.modifyMolecules().connect(4,11);
9696
ingredients.modifyMolecules().connect(4,12);
97+
ingredients.setNumOfChains(0);
98+
ingredients.setNumOfMonomersPerChain(0);
99+
ingredients.modifyMolecules()[0].setReactive(true);
100+
ingredients.modifyMolecules()[0].setNumMaxLinks(4);
101+
ingredients.modifyMolecules()[1].setReactive(true);
102+
ingredients.modifyMolecules()[1].setNumMaxLinks(4);
103+
ingredients.modifyMolecules()[2].setReactive(true);
104+
ingredients.modifyMolecules()[2].setNumMaxLinks(4);
105+
ingredients.modifyMolecules()[3].setReactive(true);
106+
ingredients.modifyMolecules()[3].setNumMaxLinks(4);
107+
ingredients.modifyMolecules()[4].setReactive(true);
108+
ingredients.modifyMolecules()[4].setNumMaxLinks(4);
109+
110+
97111

98112
REQUIRE(ingredients.getMolecules().size()==13 );
99113
REQUIRE_NOTHROW(ingredients.synchronize(ingredients));

0 commit comments

Comments
 (0)