c - Algorithm for DNA splitting -

c - Algorithm for DNA splitting -

i wanted find out effect of temperature , binding energy on splitting of dna. have used simple 1 dimensional walk describe position of deoxyribonucleic acid polymers. have set initial status such overlapping. initial position self avoiding , unidirectional random walk. polymers have 128 monomers each positions given array a[0-127] , b[0-127]. there energy e (which have taken -1) associated 2 monomers beingness @ same position. no energy in other separation distance. have used metropolis algorithm bring polymers equilibrium. have randomly selected monomer (out of 256) , flipped it. flipping has been defined as

a[i] = a[i+1] + a[i-1] + a[1]

it 'b' instead of 'a' in case of sec polymer. of course of study if final polymer selected flip defined by

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

it should noted due flip position alter either 2,0 or -2.

now metropolis algorithm states flip allowed if there no energy alter due flip (for illustration if separated monomer goes farther or if separated monomer comes closer not together). flip allowed when energy alter negative ie when after flipping 2 monomers come together. when there positive energy alter ie when 2 monomers after flip separate flip accepted probability of

powf(m_e, (e/t))

t taken 1 moment.

this algorithm iterated many times until equilibrium has been reached regards end separation distance ie b[127]-a[127]. generate random numbers, have used drand function have defined in code. since told me not random number generator tried using function ran2 have copied onto code without how working.

anyway problem equilibrium distance coming out much higher should be. ideally told should come out 0 or maybe 2 , 4 @ maximum. more unlikely. code giving values 22, 30 , on. tell me wrong? sense free inquire farther clarifications.

#include <stdio.h> #include <stdlib.h> #include <math.h> #define im1 2147483563 #define im2 2147483399 #define (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 prevent idum = 0. */ else *idum = -(*idum); idum2=(*idum); (j=ntab+7; j>=0; j--) /* load 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 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; /* in range 0..ntab-1. */ iy=iv[j]-idum2; /* here idum shuffled, idum , idum2 combined generate output. */ iv[j] = *idum; if (iy < 1) iy += imm1; if ((temp=am*iy) > rnmx) homecoming rnmx; /* because users don't expect endpoint values. */ else homecoming temp; } int main() { int a[128],b[128]; /*array defining position of polymer*/ long int i, j; /* integers defined iteration purposes*/ int r; /* rth random monomer of polymer while conducting mc algorithm*/ int x; /* new position of monomer if overcomes probability*/ float e = -1; /* energy associated overlapping monomers*/ float t = 1; /* temperature*/ int t; /*separation between final monomers*/ long idum = time(null); srand (time(null)); /*set seed random number*/ a[0]=0; b[0]=0; (i=1; i<128; i++) /*defining random overlapping initial position 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]; } } /* next metropolis algorithm*/ (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))); homecoming 0; }

i have realised error had been making. first there redundant b[i]=a[i]. not error sure bad coding habit. second. friend pointed out algorithm called value of r=0 , code uses a[r-1] give out random numbers. have create sure r doesn't take value of 0 since first monomers of both polymers kept fixed @ 0.

the final , major problem, in statements

else if (x==a[127]+2) { if (ran2(&idum)<powf(m_e,(e/t))) { b[127]=x; } }

and

else if (x==b[r]-2) { if (ran2(&idum)<powf(m_e,(e/t))) { a[r]=x; }

and r = 127 point. thought monomer flipping position next other polymer ie a[r] = b[r] new position x. after many iterations there case a[r]+4=b[r] , goes new position x x==b[r]-2. , flip accepted , not less 1 probability. new lines of code should be

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

and other 3 cases. takes care of both possibilities. code working fine now. took pain in going through such long question.

c algorithm statistics