DNA分裂算法

Question

在第二种聚合物的情况下，它将是'b'而不是'a'。当然，如果选择了最终的聚合物，则翻转的定义如下

A[127]=2*A[126]+A[127]

需要注意的是，由于翻转，位置将改变为2、0或-2。

不管怎么说，现在我的问题是，平衡距离比应该的要高得多。理想情况下，我被告知它应该是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;
}

Answer 1

我已经意识到我犯了什么错误。首先，存在一个冗余的b[i]=a[i]。可能不是一个错误，但我肯定是一个坏的编码习惯。第二。我的一个朋友指出，我的算法把r=0的值称为，我的代码使用了一个[r-1]，它会给出随机数。事实上，我必须确保r不等于0，因为两种聚合物的第一单体都固定在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;
                }

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;
                }
            }