rod_model.cpp
1 #include <cmath>
2 #include <cfloat>
3 #include <algorithm>
4 #include <functional>
5 #include <numeric>
6 #include <vector>
7 #include <utility>
8 #include "rod_model.h"
9 #include "parsing.h"
10 
11 #include <iostream>
12 #include <iomanip>
13 
14 #define PI 3.14159265
15 
16 namespace BioLCCC
17 {
18 
19 double rodAdsorptionEnergy(const std::vector<double> & rodEnergyProfile,
20  int n1, int n2) throw(BioLCCCException)
21 {
22  if ((n1 < 0) || (n1 > rodEnergyProfile.size())
23  || (n2 < 0) || (n2 > rodEnergyProfile.size()))
24  {
25  throw BioLCCCException("Index is out of range.");
26  }
27 
28  double init = 0;
29  return std::accumulate(rodEnergyProfile.begin(),
30  rodEnergyProfile.begin()+n1,
31  init)
32  + std::accumulate(rodEnergyProfile.end()-n2,
33  rodEnergyProfile.end(),
34  init);
35 }
36 
37 double partitionFunctionRodFreeSlit(double rodLength,
38  double slitWidth)
39 {
40  // the equation works only if the slit is wider than the rod
41  if (rodLength <= slitWidth)
42  {
43  // full volume without exclusion caused by walls
44  return (4 * PI * slitWidth * rodLength * rodLength -
45  // minus volume excluded by walls
46  2 * PI * rodLength * rodLength * rodLength);
47  }
48  else
49  {
50  // This case is considered in the paper as Zc.
51  return (2 * PI * slitWidth * slitWidth * rodLength);
52  }
53 
54 }
55 
56 namespace {
57  std::pair<double, double> intersection(
58  std::pair<double, double> first_line,
59  std::pair<double, double> second_line)
60  {
61  std::pair<double, double> output;
62  if (first_line.first != second_line.first)
63  {
64  double x = (second_line.second - first_line.second)
65  / (first_line.first - second_line.first);
66  x = ceil(x * 1.0e10) / 1.0e10;
67  output = std::make_pair(x, first_line.first*x + first_line.second);
68  }
69  else
70  {
71  if (first_line.first >= 0.0)
72  {
73  output = std::make_pair(-DBL_MAX, -DBL_MAX);
74  }
75  else
76  {
77  output = std::make_pair(-DBL_MAX, DBL_MAX);
78  }
79  }
80  return output;
81  }
82 }
83 
85  double segmentLength, double slitWidth, double layerWidth,
86  int N, int n1, int n2)
87 {
88  const double rodLength = (N - 1) * segmentLength;
89  std::vector<std::pair<double, double> > lowerLimits;
90  std::vector<std::pair<double, double> > upperLimits;
91 
92  // Calculating the limits of integration over theta.
93  // A limit is the minimal or maximal value of cos(theta) at which
94  // a rod conformation is still acceptable.
95 
96  // The maximal angle at which the bead next to the n1-th one is still out
97  // of the adsorbing layer.
98  upperLimits.push_back(std::make_pair(
99  -1.0 / n1 / segmentLength,
100  layerWidth / n1 / segmentLength));
101 
102  // The minimal angle theta at which the bead next to the n2-th one is still
103  // out of the adsorbing layer.
104  lowerLimits.push_back(std::make_pair(
105  -1.0 / (N - n2 - 1) / segmentLength,
106  (slitWidth - layerWidth) / (N - n2 - 1) / segmentLength));
107 
108  if (n1 > 1)
109  {
110  // The minimal angle at which the n1-th bead is still in the adsorbing
111  // layer.
112  lowerLimits.push_back(std::make_pair(
113  -1.0 / (n1 - 1) / segmentLength,
114  layerWidth / (n1 - 1) / segmentLength));
115  }
116  if (n2 != 0)
117  {
118  // The minimal angle at which the last bead does not touch the wall.
119  lowerLimits.push_back(std::make_pair(
120  -1.0 / (N - 1) / segmentLength,
121  slitWidth / (N - 1) / segmentLength));
122  // The maximal angle at which the n2-th bead is still in the adsorbing
123  // layer.
124  upperLimits.push_back(std::make_pair(
125  -1.0 / (N - n2) / segmentLength,
126  (slitWidth - layerWidth) / (N - n2) / segmentLength));
127  }
128 
129  std::vector<std::pair<double, double> > allLimits;
130  allLimits.insert(allLimits.begin(), lowerLimits.begin(), lowerLimits.end());
131  allLimits.insert(allLimits.begin(), upperLimits.begin(), upperLimits.end());
132 
133  // Finding the point at which the dependency of the limits of integration
134  // over theta on z may change.
135  std::vector<double> critPoints;
136  for (std::vector<std::pair<double, double> >::const_iterator lineIter1 =
137  allLimits.begin();
138  lineIter1 < allLimits.end();
139  ++lineIter1)
140  {
141  for (std::vector<std::pair<double, double> >::const_iterator lineIter2 =
142  lineIter1 + 1;
143  lineIter2 < allLimits.end();
144  ++lineIter2)
145  {
146  critPoints.push_back(intersection(*lineIter1, *lineIter2).first);
147  }
148  critPoints.push_back(
149  intersection(*lineIter1, std::make_pair(0.0, 1.0)).first);
150  }
151  critPoints.push_back(layerWidth);
152  critPoints.push_back(0.0);
153 
154  // Removing points lying out the range [0; layerWidth].
155  critPoints.erase(std::remove_if(critPoints.begin(), critPoints.end(),
156  bind2nd(std::less <double>(), 0)), critPoints.end());
157 
158  critPoints.erase(std::remove_if(critPoints.begin(), critPoints.end(),
159  bind2nd(std::greater <double>(), layerWidth)), critPoints.end());
160 
161  // Excluding repeating points.
162  std::sort(critPoints.begin(), critPoints.end());
163  critPoints.erase(
164  std::unique(critPoints.begin(), critPoints.end()), critPoints.end());
165 
166  // Finding the exact dependencies for each subrange of z.
167  std::vector<std::pair<double, double> > terms;
168  for (int i = 0; i < critPoints.size() - 1; i++)
169  {
170  std::pair<double, double> lowerLimit = *(lowerLimits.begin());
171  for (std::vector<std::pair<double, double> >::const_iterator lineIter =
172  lowerLimits.begin() + 1;
173  lineIter < lowerLimits.end();
174  ++lineIter)
175  {
176  if ((lineIter->first * (critPoints[i] + critPoints[i+1]) / 2.0
177  + lineIter->second) <
178  (lowerLimit.first * (critPoints[i] + critPoints[i+1]) / 2.0
179  + lowerLimit.second))
180  {
181  lowerLimit = *(lineIter);
182  }
183 
184  }
185  if ((lowerLimit.first * (critPoints[i] + critPoints[i+1]) / 2.0
186  + lowerLimit.second) > 1.0)
187  {
188  lowerLimit = std::make_pair(0.0, 1.0);
189  }
190 
191  std::pair<double, double> upperLimit = *(upperLimits.begin());
192  for (std::vector<std::pair<double, double> >::const_iterator lineIter =
193  upperLimits.begin() + 1;
194  lineIter < upperLimits.end();
195  ++lineIter)
196  {
197  if ((lineIter->first * (critPoints[i] + critPoints[i+1]) / 2.0
198  + lineIter->second) >
199  (upperLimit.first * (critPoints[i] + critPoints[i+1]) / 2.0
200  + upperLimit.second))
201  {
202  upperLimit = *(lineIter);
203  }
204 
205  }
206  if ((upperLimit.first * (critPoints[i] + critPoints[i+1]) / 2.0
207  + upperLimit.second) > 1.0)
208  {
209  upperLimit = std::make_pair(0.0, 1.0);
210  }
211 
212  terms.push_back(
213  std::make_pair(
214  lowerLimit.first - upperLimit.first,
215  lowerLimit.second - upperLimit.second));
216  }
217 
218  // Calculating the integral.
219  double integral = 0.0;
220  for (int i = 0; i < critPoints.size() - 1; i++)
221  {
222  double term =
223  (terms[i].first * (critPoints[i+1] + critPoints[i]) / 2.0
224  + terms[i].second)
225  * (critPoints[i+1] - critPoints[i]);
226  if (term > 0.0)
227  {
228  integral += term;
229  }
230  }
231 
232  return 2.0 * PI * rodLength * rodLength * integral;
233 }
234 
236  double segmentLength,
237  double slitWidth,
238  double layerWidth,
239  const std::vector<double> & rodEnergyProfile,
240  bool reversed) throw(BioLCCCException)
241 {
242  double partitionFunction = 0.0;
243  const double N = rodEnergyProfile.size();
244  for (int n1 = 1; n1 < N; n1++)
245  {
246  for (int n2 = 0; n2 <= N-n1; n2++)
247  {
248  if (reversed)
249  {
250  partitionFunction +=
251  ( (n2 == 0 ) ? 1.0 : 0.5 ) *
253  segmentLength, slitWidth, layerWidth, N, n1, n2)
254  * exp(rodAdsorptionEnergy(rodEnergyProfile, n2, n1));
255  }
256  else
257  {
258  partitionFunction +=
259  ( (n2 == 0 ) ? 1.0 : 0.5 ) *
261  segmentLength, slitWidth, layerWidth, N, n1, n2)
262  * exp(rodAdsorptionEnergy(rodEnergyProfile, n1, n2));
263  }
264  }
265  }
266  return partitionFunction;
267 }
268 
270  double segmentLength, double slitWidth, double layerWidth,
271  int N, int n1)
272 {
273  double output;
274  double rodLength = (N-1) * segmentLength;
275  if (layerWidth >= segmentLength * double(n1))
276  {
277  output = 2.0 * PI * rodLength * rodLength *
278  segmentLength / 2.0;
279  }
280  else if ((segmentLength * (double)(n1-1) < layerWidth) &&
281  (layerWidth < segmentLength * double(n1)))
282  {
283  output = 2.0 * PI * rodLength * rodLength *
284  (layerWidth
285  - segmentLength * (double)(n1-1) / 2.0
286  - layerWidth * layerWidth / 2.0 /
287  (double)n1 / segmentLength);
288  }
289  else
290  {
291  output = 2.0 * PI * rodLength * rodLength
292  * layerWidth * layerWidth / 2.0 / double(n1)
293  / double(n1-1) / segmentLength;
294  }
295  return output;
296 }
297 
299  double segmentLength,
300  double slitWidth,
301  double layerWidth,
302  const std::vector<double> & rodEnergyProfile,
303  bool reversed) throw(BioLCCCException)
304 {
305  double partitionFunction = 0.0;
306  for (unsigned int n1 = 1; n1 < rodEnergyProfile.size(); ++n1)
307  {
308  if (reversed)
309  {
310  partitionFunction +=
312  segmentLength, slitWidth, layerWidth,
313  rodEnergyProfile.size(), n1)
314  * exp(rodAdsorptionEnergy(rodEnergyProfile, 0, n1));
315  }
316  else
317  {
318  partitionFunction +=
320  segmentLength, slitWidth, layerWidth,
321  rodEnergyProfile.size(), n1)
322  * exp(rodAdsorptionEnergy(rodEnergyProfile, n1, 0));
323  }
324  }
325  return partitionFunction;
326 }
327 
328 double partitionFunctionRodFreeVolume(double rodLength,
329  double slitWidth)
330 {
331  return (4 * PI * slitWidth * rodLength * rodLength);
332 }
333 
335  const std::vector<ChemicalGroup> &parsedSequence,
336  const double secondSolventConcentration,
337  const ChemicalBasis &chemBasis,
338  const double columnPoreSize,
339  const double columnRelativeStrength,
340  const double temperature
341  ) throw(BioLCCCException)
342 {
343  if (parsedSequence.size() == 0)
344  {
345  return 0.0;
346  }
347 
348  std::vector<double> segmentEnergyProfile =
351  parsedSequence,
352  chemBasis,
353  secondSolventConcentration,
354  columnRelativeStrength,
355  temperature),
356  chemBasis.monomerLength(),
357  chemBasis.kuhnLength());
358 
359  double rodLength = chemBasis.kuhnLength() *
360  (segmentEnergyProfile.size() - 1);
361 
362  double Kd =
364  rodLength,
365  columnPoreSize - 2.0 * chemBasis.adsorptionLayerWidth())
366 
368  rodLength,
369  chemBasis.adsorptionLayerWidth())
370  * exp(rodAdsorptionEnergy(
371  segmentEnergyProfile,
372  segmentEnergyProfile.size(),
373  0));
374 
375  if (!chemBasis.neglectPartiallyDesorbedStates())
376  {
377  if (chemBasis.specialRodModel())
378  {
380  chemBasis.kuhnLength(),
381  columnPoreSize,
382  chemBasis.adsorptionLayerWidth(),
383  segmentEnergyProfile,
384  false)
385 
387  chemBasis.kuhnLength(),
388  columnPoreSize,
389  chemBasis.adsorptionLayerWidth(),
390  segmentEnergyProfile,
391  true);
392  }
393  else
394  {
396  chemBasis.kuhnLength(),
397  columnPoreSize,
398  chemBasis.adsorptionLayerWidth(),
399  segmentEnergyProfile,
400  false)
401 
403  chemBasis.kuhnLength(),
404  columnPoreSize,
405  chemBasis.adsorptionLayerWidth(),
406  segmentEnergyProfile,
407  true);
408  }
409  }
410 
412  rodLength,
413  columnPoreSize);
414 
415  return Kd;
416 }
417 
418 }
An instance of ChemicalBasis contains a set of BioLCCC constants.
Definition: chemicalbasis.h:86
double rodAdsorptionEnergy(const std::vector< double > &rodEnergyProfile, int n1, int n2)
Calculates the adsorption energy of the first n segments of a rod.
Definition: rod_model.cpp:19
Apart from classes, BioLCCC contains calculation methods and constants.
Definition: auxiliary.h:4
double partitionFunctionRodPartiallySubmergedTermSpecial(double segmentLength, double slitWidth, double layerWidth, int N, int n1)
Calculates a term in the special expression for the partition function.
Definition: rod_model.cpp:269
std::vector< double > calculateMonomerEnergyProfile(const std::vector< ChemicalGroup > &parsedSequence, const ChemicalBasis &chemBasis, const double secondSolventConcentration, const double columnRelativeStrength, const double temperature)
Calculates the effective energy profile of monomers of the polymer chain.
Definition: parsing.cpp:11
Base class for all BioLCCC exceptions. Can be used by itself.
double partitionFunctionRodPartiallySubmergedSpecial(double segmentLength, double slitWidth, double layerWidth, const std::vector< double > &rodEnergyProfile, bool reversed=false)
Calculates Z of the rod partially submerged into an adsorbing layer.
Definition: rod_model.cpp:298
double partitionFunctionRodPartiallySubmergedTermGeneral(double segmentLength, double slitWidth, double layerWidth, int N, int n1, int n2)
Calculates a term in the general expression for the partition function.
Definition: rod_model.cpp:84
double partitionFunctionRodFreeSlit(double rodLength, double slitWidth)
Calculates the partition function of a rod in a slit with impenetrable walls.
Definition: rod_model.cpp:37
double partitionFunctionRodFreeVolume(double rodLength, double slitWidth)
Calculates the partition function of a rod in a slit without walls.
Definition: rod_model.cpp:328
std::vector< double > calculateSegmentEnergyProfile(const std::vector< double > &monomerEnergyProfile, const double monomerLength, const double kuhnLength)
Calculates the effective energy profile of segments of the polymer chain.
Definition: parsing.cpp:87
double partitionFunctionRodPartiallySubmergedGeneral(double segmentLength, double slitWidth, double layerWidth, const std::vector< double > &rodEnergyProfile, bool reversed=false)
Calculates Z of a rod partially submerged into the adsorbing layer.
Definition: rod_model.cpp:235
double calculateKdRod(const std::vector< ChemicalGroup > &parsedSequence, const double secondSolventConcentration, const ChemicalBasis &chemBasis, const double columnPoreSize, const double columnRelativeStrength, const double temperature)
Calculates the coefficient of distribution of a polymer using the rod model.
Definition: rod_model.cpp:334