src/core/chain_model.cpp
00001 #include <cmath>
00002 #include <stdlib.h>
00003 #include "chain_model.h"
00004 #include "parsing.h"
00005 
00006 namespace BioLCCC
00007 {
00008 
00009 std::vector<double> calculateBoltzmannFactorProfile(
00010     const std::vector<double> &effectiveEnergyProfile)
00011 {
00012     std::vector<double> boltzmannFactorProfile; 
00013 
00014     for (std::vector<double>::const_iterator energy =
00015                 effectiveEnergyProfile.begin();
00016             energy != effectiveEnergyProfile.end();
00017             ++energy)
00018     {
00019         boltzmannFactorProfile.push_back(exp(*energy));
00020     }
00021 
00022     return boltzmannFactorProfile;
00023 }
00024 
00025 double calculateKdChain(
00026                                                 const std::vector<ChemicalGroup> &parsedSequence,
00027                                                 const double secondSolventConcentration,
00028                                                 const ChemicalBasis &chemBasis,
00029                                                 const double columnPoreSize,
00030                                                 const double columnRelativeStrength,
00031                                                 const double temperature
00032                                                 ) throw (BioLCCCException)
00033 {
00034     // At first, we need to convert the energy profile to a profile of 
00035     // distribution probabilities. Probability = exp(E_effective),
00036     // where E_effective = E_of_residue - Eab,
00037     // and Eab is an energy of binding for a binary solvent
00038     // and Eab = Ea + ln ( 1 + Nb + Nb*exp (Ea - Eb) )
00039     // also corrections of energies due to temperature (energies in exponents
00040     // are scaled to the new temperature) and the relative strength 
00041     // are introduced.
00042     std::vector<std::vector<double> > boltzmannFactorProfiles;
00043     for (std::vector<double>::const_iterator adsorptionLayerFactor = 
00044             chemBasis.adsorptionLayerFactors().begin();
00045          adsorptionLayerFactor != chemBasis.adsorptionLayerFactors().end();
00046          ++adsorptionLayerFactor) 
00047     {
00048         boltzmannFactorProfiles.push_back(
00049             calculateBoltzmannFactorProfile(
00050                 calculateSegmentEnergyProfile(
00051                     calculateMonomerEnergyProfile(
00052                         parsedSequence,
00053                         chemBasis,
00054                         secondSolventConcentration,
00055                         columnRelativeStrength * (*adsorptionLayerFactor),
00056                         temperature),
00057                     chemBasis.monomerLength(),
00058                     chemBasis.kuhnLength())));
00059     }
00060 
00061     // The size of the lattice must be greater than 
00062     // (number of adsorbing layers) * 2.
00063     // double round (double x) {return floor(x+0.5);}
00064     unsigned int latticeSize = 
00065         floor(columnPoreSize / chemBasis.kuhnLength() + 0.5);
00066 
00067     // If we want to neglect the partially desorbed states, we need to insert
00068     // two impenetrable layers right after the near-wall ones.
00069     if (chemBasis.neglectPartiallyDesorbedStates())
00070     {
00071         boltzmannFactorProfiles.push_back(
00072             std::vector<double>(
00073                 boltzmannFactorProfiles.back().size(), 0.0));
00074         latticeSize += 2;
00075     }
00076 
00077     if (latticeSize < boltzmannFactorProfiles.size() * 2)
00078     {
00079         throw BioLCCCException(
00080             "The pore size is too small for the given number of adsorbing "
00081             "layers.");
00082     }
00083 
00084     // The density vector correspond to a probability of n-th residue to be in
00085     // a certain layer between pore walls.
00086     double *density;
00087 
00088     // The transition matrix used to calculate a density vector of n-th
00089     // particle from a density vector of (n-1)-th particle.
00090     double *transitionMatrix;
00091 
00092     // The density buffer vector is used during matrix calculations.
00093     double *densityBuffer;
00094 
00095     // Memory managment.
00096     try
00097     {
00098         density = new double[latticeSize];
00099         densityBuffer = new double[latticeSize];
00100         transitionMatrix = new double[latticeSize*latticeSize];
00101     }
00102     catch (...)
00103     {
00104         throw BioLCCCException("Cannot allocate memory for calculations");
00105     }
00106 
00107     // Constructing a density vector for the first amino acid residue.
00108     // A density is distributed uniformly over all non-adsorbing layers of 
00109     // the lattice.
00110     // The density in adsorbing layers is multiplied by a Boltzmann factor of
00111     // the first segment.
00112 
00113     for (unsigned int i = 0; i < latticeSize; i++)
00114     {
00115         density[i] = 1.0;
00116     }
00117 
00118     for (unsigned int i = 0; i < boltzmannFactorProfiles.size(); ++i) 
00119     {
00120         density[i] = boltzmannFactorProfiles[i][0];
00121         density[latticeSize - i - 1] = boltzmannFactorProfiles[i][0];
00122     }
00123 
00124     // Debugging facilities.
00125     //for (unsigned int i = 0; i < latticeSize; i++)
00126     //{
00127     //    std::cout << density[i] << " ";
00128     //}
00129     //std::cout << std::endl;
00130     //std::cout << std::endl;
00131 
00132     // Than we construct a basis for the transition matrix. The basis is
00133     // a diagonal matrix with 4.0/6.0 on the main diagonal and 1.0/6.0 on
00134     // the side diagonals.
00135 
00136     // Filling the matrix.
00137     for (unsigned int i = 0; i < latticeSize; i++)
00138     {
00139         for (unsigned int j = 0; j < latticeSize; j++)
00140         {
00141             switch ((int)abs((int)(j - i)))
00142             {
00143             case 0:
00144             {
00145                 transitionMatrix[j + latticeSize * i] = 4.0/6.0;
00146                 break;
00147             }
00148             case 1:
00149             {
00150                 transitionMatrix[j + latticeSize * i] = 1.0/6.0;
00151                 break;
00152             }
00153             default:
00154                 transitionMatrix[j + latticeSize * i] = 0.0;
00155             }
00156         }
00157     }
00158 
00159     // On each step we obtain the density vector for the n-th segment
00160     // multiplying the transition matrix and the density vector of the 
00161     // (n-1)th residue.
00162     // The multiplication starts from the second segment.
00163     for (unsigned int segmentIndex = 1; 
00164          segmentIndex < boltzmannFactorProfiles[0].size();
00165          ++segmentIndex) 
00166     {
00167         // Filling the matrix elements that describe the adsorption.
00168         for (unsigned int layerIndex = 0;
00169              layerIndex < boltzmannFactorProfiles.size();
00170              ++layerIndex) 
00171         {
00172             int indexShift = layerIndex * ( latticeSize + 1 );
00173             double boltzmannFactor = 
00174                 boltzmannFactorProfiles[layerIndex][segmentIndex];
00175 
00176             transitionMatrix[indexShift + 0] = 4.0/6.0 * boltzmannFactor;
00177             transitionMatrix[indexShift + 1] = 1.0/6.0 * boltzmannFactor;
00178             transitionMatrix[latticeSize*latticeSize - 1 - indexShift - 0] =
00179                     4.0 / 6.0 * boltzmannFactor;
00180             transitionMatrix[latticeSize*latticeSize - 1 - indexShift - 1] =
00181                     1.0 / 6.0 * boltzmannFactor;
00182 
00183             // A segment may enter the second and further adsorption layers from
00184             // the inner layer (i.e. the layer lying closer to the walls).
00185             if (layerIndex > 0) 
00186             {
00187                 transitionMatrix[indexShift - 1] = 1.0/6.0 * boltzmannFactor;
00188                 transitionMatrix[latticeSize*latticeSize - 1 - indexShift + 1] =
00189                         1.0 / 6.0 * boltzmannFactor;
00190             }
00191         }
00192 
00193         // Zeroing the calculation buffer.
00194         for (unsigned int i = 0; i < latticeSize; i++)
00195         {
00196             densityBuffer[i] = 0.0;
00197         }
00198 
00199         // Multiplying the transition matrix by the density vector. The result
00200         // is stored in the buffer vector.
00201         for (unsigned int i = 0; i < latticeSize; i++)
00202         {
00203             for (unsigned int j = 0; j < latticeSize; j++)
00204             {
00205                 densityBuffer[i] = densityBuffer[i] + density[j] *
00206                                    transitionMatrix[j + i * latticeSize];
00207             }
00208         }
00209 
00210         // Transferring the results from the density vector.
00211         for (unsigned int i = 0; i < latticeSize; i++)
00212         {
00213             density[i] = densityBuffer[i];
00214         }
00215     }
00216 
00217     // Kd is a sum of elements of the density vector, normalized to the size 
00218     // of the lattice.
00219     double Kd=0;
00220     for (unsigned int i=0; i < latticeSize; i++)
00221     {
00222         Kd += density[i];
00223     }
00224 
00225     if (chemBasis.neglectPartiallyDesorbedStates())
00226     {
00227         Kd = Kd / (double)(latticeSize - 2);
00228     }
00229     else
00230     {
00231         Kd = Kd / (double)(latticeSize);
00232     }
00233 
00234     // Cleaning memory.
00235     try
00236     {
00237         delete[] density;
00238         delete[] densityBuffer;
00239         delete[] transitionMatrix;
00240     }
00241     catch (...)
00242     {
00243         throw BioLCCCException("Cannot allocate memory for calculations");
00244     }
00245 
00246     return Kd;
00247 }
00248 
00249 }
 All Classes Namespaces Functions Variables Enumerations Enumerator