-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathcloningR.cpp
More file actions
147 lines (121 loc) · 4.47 KB
/
cloningR.cpp
File metadata and controls
147 lines (121 loc) · 4.47 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#include <iostream>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "cloningserial.hpp"
#include "env.hpp"
#include "particle.hpp"
#include "readwrite.hpp"
int main() {
// cloning parameters
double tmax = getEnvDouble("TMAX", 1); // dimensionless time to simulate
int nc = getEnvInt("NC", 10); // number of clones
double sValue = getEnvDouble("SVALUE", 0); // biasing parameter
int seed = getEnvInt("SEED", 0); // master random seed
int nRuns = getEnvInt("NRUNS", 1); // number of different runs
int cloneMethod = 2; // should keep this set to 2 (!!)
int initSim = getEnvInt("INITSIM", 1); // number of initial elementary number of iterations to "randomise" the systems
// openMP parameters
#ifdef _OPENMP
int threads = getEnvInt("THREADS", -1); // number of threads
printf("# compiled with openMP\n");
if ( threads > 0 ) {
printf("# setting threads %d\n",threads);
omp_set_num_threads(threads);
}
printf("# running on %d threads\n",omp_get_max_threads());
#endif
// physical parameters
int N = getEnvInt("N", 100); // number of rotors in the system
double Dr = getEnvDouble("DR", 1.0/2.0); // rotational diffusivity
// simulation parameters
int tau = getEnvInt("TAU", 100); // elementary number of steps
double dt = getEnvDouble("DT", 0.001); // time step
// output to file
std::string filename = getEnvString("FILE", ""); // output file name
Write output(filename); // output class
output.write<double>(tmax);
output.write<int>(nc);
output.write<double>(sValue);
output.write<int>(seed);
output.write<int>(nRuns);
output.write<int>(cloneMethod);
output.write<int>(initSim);
output.write<int>(N);
output.write<double>(Dr);
output.write<int>(tau);
output.write<double>(dt);
output.write<int>(BIAS);
// dummy system
Rotors dummy(N, Dr, dt);
printf("## CloningSerial Code: tmax %.3e numClones %d runs %d s %.3e tau %d Delta.t %.3e\n",tmax,nc,nRuns,sValue,tau,dt);
// cloning object (default cloning method is [eq])
CloningSerial<Rotors> clones(nc, tau, sValue, cloneMethod);
// set up the clones etc, using dummySystem to get system sizes, hop rates, etc
clones.init(&dummy, seed);
std::cout << "## master seed " << seed << std::endl;
#if BIAS == 0
double sFactor = N*tau*dt;
double sOffset = 0;
#endif
#if BIAS == 1
#ifdef CONTROLLED_DYNAMICS
double sFactor = sValue/Dr;
double sOffset = tau*dt;
#else
double sFactor = N*tau*dt;
double sOffset = 0;
#endif
#endif
for (int run = 0; run<nRuns;run++) {
// go! (this includes generating "random" [different] initial conditions for the clones)
clones.doCloning(tmax, initSim,
// ITERATION FUNCTION
[](Rotors* rotors, int Niter) { iterate_rotors(rotors, Niter); },
// GET WEIGHT FUNCTION
[&sFactor, &sOffset](Rotors* rotors) {
double sWeight = 0;
// biasing with order parameter
#if BIAS == 0
sWeight = rotors->getBiasingParameter() // sw = s
*rotors->getOrder() // *nu
*sFactor; // *N*tau
#endif
// biasing with squared order parameter
#if BIAS == 1
#ifdef CONTROLLED_DYNAMICS
sWeight = rotors->getBiasingParameter() // sw = s
*(sOffset - sFactor*rotors->getBiasIntegral()); // *(tau - s/Dr)*int nu^2 sin(theta-phi)^2)
#else
sWeight = rotors->getBiasingParameter() // sw = s
*rotors->getOrderSq() // *nu^2
*sFactor; // *N*tau
#endif
#endif
return sWeight;
},
// CONTROL FUNCTION
[](std::vector<Rotors*>& rotors, int pullOffset, int pushOffset) {;}
);
clones.outputOP.assign(2, 0.0);
for (int i=0; i < nc; i++) {
clones.outputOP[0] += (clones.finalSystem(i))->getTotalOrder()[0]
/(clones.finalSystem(i))->getDump()[0]; // order parameter
clones.outputOP[1] += (clones.finalSystem(i))->getTotalOrderSq()[0]
/(clones.finalSystem(i))->getDump()[0]; // squared order parameter
}
for (unsigned int j=0;j<2;j++) { clones.outputOP[j] /= nc; }
std::cout << std::endl;
std::cout << "##s " << sValue << std::endl
<< "##bias " << BIAS << std::endl
<< "#SCGF " << clones.outputPsi/N << std::endl
<< "#nu " << clones.outputOP[0] << std::endl
<< "#nu^2 " << clones.outputOP[1] << std::endl << std::endl
<< "##time " << clones.outputWalltime << std::endl;
// output to file
output.write<double>(clones.outputPsi);
output.write<double>(clones.outputOP[0]);
output.write<double>(clones.outputOP[1]);
output.write<double>(clones.outputWalltime);
}
}