我想找出温度和结合能对 DNA split 的影响。 我使用了一个简单的一维漫游来描述 DNA 聚合物的位置。我已经设置了初始条件，使它们重叠。初始位置是 self 回避和单向随机游走。这些聚合物有 128 个单体，每个单体的位置由数组 a[0-127] 和 b[0-127] 给出。存在与处于同一位置的两个单体相关联的能量 E(我已将其取为 -1)。所有其他间隔距离都没有能量。 现在我已经使用 Metropolis 算法使聚合物达到平衡。 我随机选择了一个单体(从 256 个中)并翻转了它。翻转被定义为
a[i] = a[i+1] + a[i-1] + a[1]
如果是第二种聚合物，它将是“b”而不是“a”。 当然，如果选择了最终的聚合物，翻转将定义为
a[127] = 2*a[126] + a[127]
应该注意，由于翻转，位置将变为 2,0 或 -2。
现在 metropolis 算法指出，如果由于翻转而没有能量变化，则始终允许翻转(例如，如果一个已经分离的单体进一步移动，或者如果一个分离的单体靠近但没有完全结合在一起)。 当能量变化为负时，即当翻转后两个单体聚集在一起时，也总是允许翻转。 当存在正能量变化时，即当最初两个单体在一起但在翻转后它们分开时，翻转被接受的概率为
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;
}
最佳答案
我已经意识到我犯了什么错误。 首先有一个多余的b[i]=a[i]。 可能不是错误，而是我确信一个糟糕的编码习惯。 第二。我的 friend 指出我的算法将 r=0 的值称为 as 并且我的代码使用 a[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;
}
还有 where r = 127 但你明白了。 实际上，我认为单体只能从与其他聚合物相邻的位置(即 a[r] = b[r])翻转到新位置 x。但在多次迭代之后，也会出现 a[r]+4=b[r] 的情况，它会到达一个新位置 x，其中 x==b[r]-2。 而且这种翻转总是会被接受，而且概率不会小于 1。 所以新的代码行应该是
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 个案例也类似。这会照顾到这两种可能性。该代码现在工作正常。 感谢所有为解决这么长的问题而付出痛苦的人。
关于c - DNA split 算法，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/24382425/