2014-06-24 49 views
1

我想知道温度和结合能对DNA分裂的影响。 我用一个简单的一维步行来描述DNA聚合物的位置。我已经设定了它们重叠的初始条件。初始位置是一种自我避免和单向随机游走。该聚合物具有128个单体,每个单体的位置由阵列a [0-127]和b [0-127]给出。有两种单体在相同位置上的能量E(我已经采用-1)。所有其他分隔距离都没有能量。 现在我已经使用Metropolis算法来使聚合物达到平衡。 我已经随机选择了一个单体(256个)并翻转了它。翻转已被定义为DNA分裂算法

A [1] = A [1 + 1] + A [I-1] + A [1]

这将是 'B' 而不是“一个'在第二种聚合物的情况下。 当然,如果所选择的最终聚合物被倒装将由

来定义[127] = 2 *一个[126] + A [127]

应当注意的是,由于翻转位置将改变为2,0或-2。

现在,大都会算法规定,如果由于翻转而没有能量变化(例如,如果已经分离的单体进一步变得更远,或者如果分离的单体变得更近但不完全在一起),将总是允许翻盖。 当能量变化为负值时,也就是在翻转两个单体后,总是可以翻转。 当有一个积极的能量变化即,当最初两种单体在一起但翻转后,他们然后分离该倒装被接受与

函数powf(M_E,(E/T))的概率

此时T也取1。

该算法迭代很多次,直到达到平衡已经达到了结束间隔距离,即b [127] -a [127]。 为了生成随机数,我使用了我在代码中定义的drand函数。由于有人告诉我这可能不是一个非常好的随机数生成器,我还尝试使用函数ran2,我刚刚将代码复制到代码中,但没有如何工作。

无论如何,现在我的问题是平衡距离出来远远高于应该是。理想情况下,我被告知它应该最多为0或者2和4。比这更不可能。但是我的代码非常频繁地给出像22,30等值。 有人能告诉我什么是错的?随意要求进一步澄清。

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#define IM1 2147483563 
#define IM2 2147483399 
#define AM (1.0/IM1) 
#define IMM1 (IM1-1) 
#define IA1 40014 
#define IA2 40692 
#define IQ1 53668 
#define IQ2 52774 
#define IR1 12211 
#define IR2 3791 
#define NTAB 32 
#define NDIV (1+IMM1/NTAB) 
#define EPS 1.2e-7 
#define RNMX (1.0-EPS) 


float drand() 
{ 
    float f, r, randmax; 
    r = rand(); 
    randmax = RAND_MAX; 
    f = r/(randmax+1); 
    return(f); 
} 

double ran2(long *idum) 
{ 
    int j; 
    long k; 
    static long idum2=123456789; 
    static long iy=0; 
    static long iv[NTAB]; 
    double temp; 

    if (*idum <= 0)   /* Initialize. */ 
    { 
     if (-(*idum) < 1) *idum=1; /*Be sure to prevent idum = 0. */ 
     else *idum = -(*idum); 
     idum2=(*idum); 
     for (j=NTAB+7; j>=0; j--) /* Load the shuffle table (after 8 warm-ups).*/ 
     { 
      k=(*idum)/IQ1; 
      *idum=IA1*(*idum-k*IQ1)-k*IR1; 
      if (*idum < 0) *idum += IM1; 
      if (j < NTAB) iv[j] = *idum; 
     } 
     iy=iv[0]; 
    } 
    k=(*idum)/IQ1; /* Start here when not initializing.*/ 
    *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without 
            overflows by Schrage's method. */ 
    if (*idum < 0) *idum += IM1; 
    k=idum2/IQ2; 
    idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise. */ 
    if (idum2 < 0) idum2 += IM2; 
    j=iy/NDIV;     /* Will be in the range 0..NTAB-1. */ 
    iy=iv[j]-idum2;    /* Here idum is shuffled, idum and idum2 are 
            combined to generate output. */ 
    iv[j] = *idum; 
    if (iy < 1) iy += IMM1; 
    if ((temp=AM*iy) > RNMX) 
     return RNMX;    /* Because users don't expect endpoint values. */ 
    else return temp; 
} 


int main() 
{ 
    int a[128],b[128]; /*array defining position of polymer*/ 
    long int i, j;   /* integers defined for iteration purposes*/ 
    int r;    /* The rth random monomer of the polymer while conducting the MC algorithm*/ 
    int x;    /* The new position of the monomer if it overcomes the probability*/ 
    float E = -1;  /* Energy associated with overlapping monomers*/ 
    float T = 1;  /* Temperature*/ 
    int t;    /*separation between final monomers*/ 
    long idum = time(NULL); 
    srand (time(NULL)); /*set seed for the random number*/ 
    a[0]=0; 
    b[0]=0; 
    for (i=1; i<128; i++) /*Defining a random but overlapping initial position for the strands*/ 
    { 
     if (ran2(&idum)<0.5) 
     { 
      a[i]=a[i-1]+1; 
      b[i]=a[i]; 
     } 
     else 
     { 
      a[i]=b[i]=a[i-1]-1; 
      b[i]=a[i]; 
     } 
    } 

    /* Following is the metropolis algorithm*/ 
    for (j=1; j<1000000; j=j+1) 
    { 
     r = floor(ran2(&idum)*128); 
     if (ran2(&idum)<0.5) 
     { 
      if (r<=126) 
      { 
       x=a[r+1]+a[r-1]-a[r]; 
       if (x==b[r]) 
       { 
        a[r]=x; 
       } 
       else if (x==b[r]-2) 
       { 
        if (ran2(&idum)<powf(M_E,(E/T))) 
        { 
         a[r]=x; 
        } 
       } 
       else if (x<b[r]-2) 
       { 
        a[r]=x; 
       } 
      } 
      else if (r==127) 
      { 
       x=2*a[126]-a[127]; 
       if (x==b[127]) 
       { 
        a[127]=x; 
       } 
       else if (x==b[127]-2) 
       { 
        if (ran2(&idum)<powf(M_E,(E/T))) 
        { 
         a[127]=x; 
        } 
       } 
       else if (x<b[127]-2) 
       { 
        a[127]=x; 
       } 
      } 
     } 
     else 
     { 
      if (r<=126) 
      { 
       x=b[r+1]+b[r-1]-b[r]; 
       if (x==a[r]) 
       { 
        b[r]=x; 
       } 
       else if (x==a[r]+2) 
       { 
        if (ran2(&idum)<powf(M_E,(E/T))) 
        { 
         b[r]=x; 
        } 
       } 
       else if (x>a[r]+2) 
       { 
        b[r]=x; 
       } 
      } 
      else if (r==127) 
      { 
       x=2*b[126]-b[127]; 
       if (x==a[127]) 
       { 
        b[127]=x; 
       } 
       else if (x==a[127]+2) 
       { 
        if (ran2(&idum)<powf(M_E,(E/T))) 
        { 
         b[127]=x; 
        } 
       } 
       else if (x>a[127]+2) 
       { 
        b[127]=x; 
       } 
      } 
     } 
     t = b[127]-a[127]; 
     if (j%(25600)==0) 
     { 
      printf("%d\n", t); 
     } 
    } 
    printf("%f\n", powf(M_E,(E/T))); 
    return 0; 
} 
+0

所以你问我们分析你的算法,然后找到它的问题?我认为这个问题太多了,因为你的算法非常广泛。你的代码看起来非常好/干净;做得好!你不能放大,或削减部分处理和检查中间结果? –

+0

是的。我确实意识到代码很长,但在我发布之前尽量找到问题。 – aniruud

+0

我确实尝试隔离某些部件,以便找到问题。随机数发生器似乎没有任何问题。初始设置聚合物也是如此。问题必须在大都会算法中。当迭代次数较低时,结果更加明智。但是当我们增加迭代次数时,尽管事实上它们应该在某个点达到平衡,并且在后续迭代之后才停留在那里,但值会变得更加怪异和怪异。 – aniruud

回答

1

我已经意识到我一直在做什么错误。 首先有一个冗余b [i] = a [i]。 可能不是一个错误,但我相信一个不好的编码习惯。其次。我的朋友指出,我的算法称为r = 0的值,而我的代码使用[r-1],然后发出随机数。由于两种聚合物的第一单体保持固定在0,所以我实际上必须确保r不占0的值。

最后也是最重要的问题,是文中陈述

else if (x==a[127]+2) 
      { 
       if (ran2(&idum)<powf(M_E,(E/T))) 
       { 
        b[127]=x; 
       } 
      } 

else if (x==b[r]-2) 
      { 
       if (ran2(&idum)<powf(M_E,(E/T))) 
       { 
        a[r]=x; 
       } 

而且其中r = 127,但你明白了吧。实际上我认为单体只能从与其他聚合物相邻的位置翻转,即a [r] = b [r]到新的位置x。但经过多次迭代之后,还会出现a [r] + 4 = b [r]的情况,并且它会进入x = b [r] -2的新位置x。 而这个翻盖总是会被接受,而不是以一定的概率。 因此,代码的新的生产线应该是

else if (x==b[r]-2) 
      { 
       if (a[r]==b[r]) 
       { 
        if (drand()<powf(M_E,(E/T))) 
       { 
        a[r]=x; 
       } 
       } 
       else 
       { 
        a[r]=x; 
       } 
      } 

,类似的还有另外3例。这考虑到了两种可能性。代码现在工作正常。 感谢所有那些经历了这么长的问题而痛苦的人。