src/core/rod_model.cpp
00001 #include <cmath>
00002 #include <cfloat>
00003 #include <algorithm>
00004 #include <functional>
00005 #include <numeric>
00006 #include <vector>
00007 #include <utility>
00008 #include "rod_model.h"
00009 #include "parsing.h"
00010 
00011 #include <iostream>
00012 #include <iomanip>
00013 
00014 #define PI 3.14159265
00015 
00016 namespace BioLCCC
00017 {
00018 
00019 double rodAdsorptionEnergy(const std::vector<double> & rodEnergyProfile,
00020                            int n1, int n2) throw(BioLCCCException)
00021 {
00022     if ((n1 < 0) || (n1 > rodEnergyProfile.size())
00023         || (n2 < 0) || (n2 > rodEnergyProfile.size()))
00024     {
00025         throw BioLCCCException("Index is out of range.");
00026     }
00027 
00028     double init = 0;
00029     return std::accumulate(rodEnergyProfile.begin(),
00030                            rodEnergyProfile.begin()+n1,
00031                            init)
00032            + std::accumulate(rodEnergyProfile.end()-n2,
00033                            rodEnergyProfile.end(),
00034                            init);
00035 }
00036 
00037 double partitionFunctionRodFreeSlit(double rodLength,
00038                                     double slitWidth)
00039 {
00040     // the equation works only if the slit is wider than the rod
00041     if (rodLength <= slitWidth)
00042     {
00043         // full volume without exclusion caused by walls
00044         return (4 * PI * slitWidth * rodLength * rodLength -
00045                 // minus volume excluded by walls
00046                 2 * PI * rodLength * rodLength * rodLength);
00047     }
00048     else
00049     {
00050         // This case is considered in the paper as Zc.
00051         return (2 * PI * slitWidth * slitWidth * rodLength);
00052     }
00053 
00054 }
00055 
00056 namespace {
00057     std::pair<double, double> intersection(
00058         std::pair<double, double> first_line,
00059         std::pair<double, double> second_line)
00060     {
00061         std::pair<double, double> output;
00062         if (first_line.first != second_line.first)
00063         {
00064             double x = (second_line.second - first_line.second) 
00065                        / (first_line.first - second_line.first);
00066             x = ceil(x * 1.0e10) / 1.0e10;
00067             output = std::make_pair(x, first_line.first*x + first_line.second);
00068         }
00069         else
00070         {
00071             if (first_line.first >= 0.0)
00072             {
00073                 output = std::make_pair(-DBL_MAX, -DBL_MAX);
00074             }
00075             else
00076             {
00077                 output = std::make_pair(-DBL_MAX, DBL_MAX);
00078             }
00079         }
00080         return output;
00081     }
00082 }
00083 
00084 double partitionFunctionRodPartiallySubmergedTermGeneral(
00085     double segmentLength, double slitWidth, double layerWidth,
00086     int N, int n1, int n2)
00087 {
00088     const double rodLength = (N - 1) * segmentLength;
00089     std::vector<std::pair<double, double> > lowerLimits;
00090     std::vector<std::pair<double, double> > upperLimits;
00091 
00092     // Calculating the limits of integration over theta.
00093     // A limit is the minimal or maximal value of cos(theta) at which 
00094     // a rod conformation is still acceptable. 
00095 
00096     // The maximal angle at which the bead next to the n1-th one is still out
00097     // of the adsorbing layer.
00098     upperLimits.push_back(std::make_pair(
00099         -1.0 / n1 / segmentLength,
00100         layerWidth / n1 / segmentLength));
00101 
00102     // The minimal angle theta at which the bead next to the n2-th one is still
00103     // out of the adsorbing layer.
00104     lowerLimits.push_back(std::make_pair(
00105         -1.0 / (N - n2 - 1) / segmentLength,
00106         (slitWidth - layerWidth) / (N - n2 - 1) / segmentLength));
00107 
00108     if (n1 > 1)
00109     {
00110         // The minimal angle at which the n1-th bead is still in the adsorbing
00111         // layer.
00112         lowerLimits.push_back(std::make_pair(
00113             -1.0 / (n1 - 1) / segmentLength,
00114             layerWidth / (n1 - 1) / segmentLength));
00115     }
00116     if (n2 != 0)
00117     {
00118         // The minimal angle at which the last bead does not touch the wall.  
00119         lowerLimits.push_back(std::make_pair(
00120             -1.0 / (N - 1) / segmentLength,
00121             slitWidth / (N - 1) / segmentLength));
00122         // The maximal angle at which the n2-th bead is still in the adsorbing
00123         // layer.  
00124         upperLimits.push_back(std::make_pair(
00125             -1.0 / (N - n2) / segmentLength,
00126             (slitWidth - layerWidth) / (N - n2) / segmentLength));
00127     }
00128 
00129     std::vector<std::pair<double, double> > allLimits;
00130     allLimits.insert(allLimits.begin(), lowerLimits.begin(), lowerLimits.end());
00131     allLimits.insert(allLimits.begin(), upperLimits.begin(), upperLimits.end());
00132 
00133     // Finding the point at which the dependency of the limits of integration
00134     // over theta on z may change.
00135     std::vector<double> critPoints;
00136     for (std::vector<std::pair<double, double> >::const_iterator lineIter1 =
00137             allLimits.begin();
00138          lineIter1 < allLimits.end();
00139          ++lineIter1)
00140     {
00141         for (std::vector<std::pair<double, double> >::const_iterator lineIter2 =
00142                 lineIter1 + 1;
00143              lineIter2 < allLimits.end();
00144              ++lineIter2)
00145         {
00146             critPoints.push_back(intersection(*lineIter1, *lineIter2).first);
00147         }
00148         critPoints.push_back(
00149             intersection(*lineIter1, std::make_pair(0.0, 1.0)).first);
00150     }
00151     critPoints.push_back(layerWidth);
00152     critPoints.push_back(0.0);
00153 
00154     // Removing points lying out the range [0; layerWidth].
00155     critPoints.erase(std::remove_if(critPoints.begin(), critPoints.end(),
00156         bind2nd(std::less <double>(), 0)), critPoints.end());
00157 
00158     critPoints.erase(std::remove_if(critPoints.begin(), critPoints.end(),
00159         bind2nd(std::greater <double>(), layerWidth)), critPoints.end());
00160 
00161     // Excluding repeating points.
00162     std::sort(critPoints.begin(), critPoints.end());
00163     critPoints.erase(
00164         std::unique(critPoints.begin(), critPoints.end()), critPoints.end());
00165 
00166     // Finding the exact dependencies for each subrange of z.
00167     std::vector<std::pair<double, double> > terms;
00168     for (int i = 0; i < critPoints.size() - 1; i++) 
00169     {
00170         std::pair<double, double> lowerLimit = *(lowerLimits.begin());
00171         for (std::vector<std::pair<double, double> >::const_iterator lineIter =
00172                 lowerLimits.begin() + 1;
00173              lineIter < lowerLimits.end();
00174              ++lineIter)
00175         {
00176             if ((lineIter->first * (critPoints[i] + critPoints[i+1]) / 2.0
00177                    + lineIter->second) <
00178                 (lowerLimit.first * (critPoints[i] + critPoints[i+1]) / 2.0
00179                    + lowerLimit.second))
00180             {
00181                 lowerLimit = *(lineIter);
00182             }
00183 
00184         }
00185         if ((lowerLimit.first * (critPoints[i] + critPoints[i+1]) / 2.0
00186                + lowerLimit.second) > 1.0)
00187         {
00188             lowerLimit = std::make_pair(0.0, 1.0);
00189         }
00190 
00191         std::pair<double, double> upperLimit = *(upperLimits.begin());
00192         for (std::vector<std::pair<double, double> >::const_iterator lineIter =
00193                 upperLimits.begin() + 1;
00194              lineIter < upperLimits.end();
00195              ++lineIter)
00196         {
00197             if ((lineIter->first * (critPoints[i] + critPoints[i+1]) / 2.0
00198                    + lineIter->second) >
00199                 (upperLimit.first * (critPoints[i] + critPoints[i+1]) / 2.0
00200                    + upperLimit.second))
00201             {
00202                 upperLimit = *(lineIter);
00203             }
00204             
00205         }
00206         if ((upperLimit.first * (critPoints[i] + critPoints[i+1]) / 2.0
00207                + upperLimit.second) > 1.0)
00208         {
00209             upperLimit = std::make_pair(0.0, 1.0);
00210         }
00211 
00212         terms.push_back(
00213             std::make_pair(
00214                 lowerLimit.first - upperLimit.first,
00215                 lowerLimit.second - upperLimit.second));
00216     }
00217 
00218     // Calculating the integral.
00219     double integral = 0.0;
00220     for (int i = 0; i < critPoints.size() - 1; i++) 
00221     {
00222         double term = 
00223             (terms[i].first * (critPoints[i+1] + critPoints[i]) / 2.0
00224                 + terms[i].second)
00225             * (critPoints[i+1] - critPoints[i]);
00226         if (term > 0.0)
00227         {
00228             integral += term;
00229         }
00230     }
00231 
00232     return 2.0 * PI * rodLength * rodLength * integral;
00233 }
00234 
00235 double partitionFunctionRodPartiallySubmergedGeneral(
00236     double segmentLength,
00237     double slitWidth,
00238     double layerWidth,
00239     const std::vector<double> & rodEnergyProfile,
00240     bool reversed) throw(BioLCCCException)
00241 {
00242     double partitionFunction = 0.0;
00243     const double N = rodEnergyProfile.size();
00244     for (int n1 = 1; n1 < N; n1++)
00245     {
00246         for (int n2 = 0; n2 <= N-n1; n2++)
00247         {
00248             if (reversed)
00249             {
00250                 partitionFunction +=
00251                     ( (n2 == 0 ) ? 1.0 : 0.5 ) *
00252                     partitionFunctionRodPartiallySubmergedTermGeneral(
00253                         segmentLength, slitWidth, layerWidth, N, n1, n2)
00254                     * exp(rodAdsorptionEnergy(rodEnergyProfile, n2, n1));
00255             }
00256             else
00257             {
00258                 partitionFunction +=
00259                     ( (n2 == 0 ) ? 1.0 : 0.5 ) *
00260                     partitionFunctionRodPartiallySubmergedTermGeneral(
00261                         segmentLength, slitWidth, layerWidth, N, n1, n2)
00262                     * exp(rodAdsorptionEnergy(rodEnergyProfile, n1, n2));
00263             }
00264         }
00265     }
00266     return partitionFunction;
00267 }
00268 
00269 double partitionFunctionRodPartiallySubmergedTermSpecial(
00270     double segmentLength, double slitWidth, double layerWidth,
00271     int N, int n1)
00272 {
00273     double output;
00274     double rodLength = (N-1) * segmentLength;
00275     if (layerWidth >= segmentLength * double(n1))
00276     {
00277         output = 2.0 * PI * rodLength * rodLength *
00278                              segmentLength / 2.0;
00279     }
00280     else if ((segmentLength * (double)(n1-1) < layerWidth) &&
00281              (layerWidth < segmentLength * double(n1)))
00282     {
00283         output = 2.0 * PI * rodLength * rodLength *
00284                              (layerWidth
00285                               - segmentLength * (double)(n1-1) / 2.0
00286                               - layerWidth * layerWidth / 2.0 / 
00287                                   (double)n1 / segmentLength);
00288     }
00289     else
00290     {
00291         output = 2.0 * PI * rodLength * rodLength 
00292                  * layerWidth * layerWidth / 2.0 / double(n1)
00293                  / double(n1-1) / segmentLength;
00294     }
00295     return output;
00296 }
00297 
00298 double partitionFunctionRodPartiallySubmergedSpecial(
00299     double segmentLength,
00300     double slitWidth,
00301     double layerWidth,
00302     const std::vector<double> & rodEnergyProfile,
00303     bool reversed) throw(BioLCCCException)
00304 {
00305     double partitionFunction = 0.0;
00306     for (unsigned int n1 = 1; n1 < rodEnergyProfile.size(); ++n1)
00307     {
00308         if (reversed)
00309         {
00310             partitionFunction +=
00311                 partitionFunctionRodPartiallySubmergedTermSpecial(
00312                     segmentLength, slitWidth, layerWidth, 
00313                     rodEnergyProfile.size(), n1)
00314                 * exp(rodAdsorptionEnergy(rodEnergyProfile, 0, n1));
00315         }
00316         else
00317         {
00318             partitionFunction += 
00319                 partitionFunctionRodPartiallySubmergedTermSpecial(
00320                     segmentLength, slitWidth, layerWidth, 
00321                     rodEnergyProfile.size(), n1)
00322                 * exp(rodAdsorptionEnergy(rodEnergyProfile, n1, 0));
00323         }
00324     }
00325     return partitionFunction;
00326 }
00327 
00328 double partitionFunctionRodFreeVolume(double rodLength,
00329                                       double slitWidth)
00330 {
00331     return (4 * PI * slitWidth * rodLength * rodLength);
00332 }
00333 
00334 double calculateKdRod(
00335                                           const std::vector<ChemicalGroup> &parsedSequence,
00336                                           const double secondSolventConcentration,
00337                                           const ChemicalBasis &chemBasis,
00338                                           const double columnPoreSize,
00339                                           const double columnRelativeStrength,
00340                                           const double temperature
00341                                           ) throw(BioLCCCException)
00342 {
00343     if (parsedSequence.size() == 0)
00344     {
00345         return 0.0;
00346     }
00347 
00348     std::vector<double> segmentEnergyProfile = 
00349         calculateSegmentEnergyProfile(
00350             calculateMonomerEnergyProfile(
00351                 parsedSequence,
00352                 chemBasis,
00353                 secondSolventConcentration,
00354                 columnRelativeStrength,
00355                 temperature),
00356             chemBasis.monomerLength(),
00357             chemBasis.kuhnLength());
00358 
00359     double rodLength = chemBasis.kuhnLength() *
00360         (segmentEnergyProfile.size() - 1);
00361 
00362     double Kd =
00363             partitionFunctionRodFreeSlit(
00364                 rodLength,
00365                 columnPoreSize - 2.0 * chemBasis.adsorptionLayerWidth())
00366 
00367             + 2.0 * partitionFunctionRodFreeSlit(
00368                 rodLength,
00369                 chemBasis.adsorptionLayerWidth())
00370             * exp(rodAdsorptionEnergy(
00371                 segmentEnergyProfile,
00372                 segmentEnergyProfile.size(), 
00373                 0));
00374 
00375     if (!chemBasis.neglectPartiallyDesorbedStates())
00376     {
00377         if (chemBasis.specialRodModel())
00378         {
00379             Kd += 2.0 * partitionFunctionRodPartiallySubmergedSpecial(
00380                     chemBasis.kuhnLength(),
00381                     columnPoreSize,
00382                     chemBasis.adsorptionLayerWidth(),
00383                     segmentEnergyProfile,
00384                     false)
00385 
00386                   + 2.0 * partitionFunctionRodPartiallySubmergedSpecial(
00387                       chemBasis.kuhnLength(),
00388                       columnPoreSize,
00389                       chemBasis.adsorptionLayerWidth(),
00390                       segmentEnergyProfile,
00391                       true);
00392         }
00393         else
00394         {
00395             Kd += 2.0 * partitionFunctionRodPartiallySubmergedGeneral(
00396                     chemBasis.kuhnLength(),
00397                     columnPoreSize,
00398                     chemBasis.adsorptionLayerWidth(),
00399                     segmentEnergyProfile,
00400                     false)
00401 
00402                   + 2.0 * partitionFunctionRodPartiallySubmergedGeneral(
00403                       chemBasis.kuhnLength(),
00404                       columnPoreSize,
00405                       chemBasis.adsorptionLayerWidth(),
00406                       segmentEnergyProfile,
00407                       true);
00408         }
00409     }
00410 
00411     Kd /= partitionFunctionRodFreeVolume(
00412             rodLength,
00413             columnPoreSize);
00414 
00415     return Kd;
00416 }
00417 
00418 }
 All Classes Namespaces Functions Variables Enumerations Enumerator