code for Rouse dynamics

Living_Polymers.DNA.digestion.nohistory.v1.c++

+/*
+v1. I try to implement a different algorithm in which I do not follow the evolution of the polymer but pick the ageing block T, perform N(T) cuts and follow relaxation with Rouse dynamics.
+NOTE: even though the recording is 10 minutes, the MSD is ~30 sec so only 15 times the relaxation dynamics.
+NOTE: I have defined the relaxation time as L^2/D0 but int reality it is L^2/(\pi^2 D0) so the survival probability plotted in units of Trep should decay as exp(-x/(1/pi^2)) and it does.
+*/
+
+#include<iostream>
+#include<stdlib.h>
+#include<cmath>
+#include<sstream>
+#include<fstream>
+#include<string>
+#include<iomanip>
+#include<vector>
+#include <algorithm>    // std::shuffle/swap/sort
+#include <array>        // std::array
+
+using namespace std;
+
+/* screen output */
+int speak=1;
+
+/* CONSTANTS */
+const int REPMAX=100;
+const int NPARTICLES=1000;
+const double dL=0.01; //
+const double dt=0.001; //integration time
+const double D0=0.01; //monomer diffusion. In units of L0^2/tau
+const double L0=1; //L0 sets the units of length = lambda-dna = 48 kb
+
+const int Trep=int(pow(L0,2.0)/D0); //Rouse time units [tau]
+int qdump=int(1.0/dt);
+
+int Tblock=int(300*Trep); //time per survival block (300x2sec=10min) units [tau]
+int Tmsd=int(15*Trep); //time of MSD is actually only 15x relaxation time
+// #of dt steps in each Tblock (= 300Trep/dt=300/dt/D0 or larger)
+const int Tblockmax=10000000;
+
+const int Nblocksmax=10;
+
+
+/* VARIABLES */
+int t;
+double L;
+
+/*RESTRICTION SITES */
+/* cumulative positions */
+double r0=L0*0.1135;
+double r1=L0*0.4607;
+double r2=L0*0.5766;
+double r3=L0*0.7113;
+double r4=L0*0.8604;
+double REsites[5]={r0,r1,r2,r3,r4};
+double REsite_iscut[5]={0,0,0,0,0}; //is_cut
+double fragment[6][2]; // up to 6 fragments -- contains 0=b1,1=b2
+
+int nREsitesmax=5;
+const int ncutmax=6; // fragments = sites + 1
+double cut[ncutmax];
+
+//start from (see below)
+//fragment[0][0]=0;
+//fragment[0][1]=L;
+int nfrags=1;
+
+/* ARRAYS */
+double pos[Tblockmax]; //this is updated every dt
+int Survival[Tblockmax]; //[Tblock] => 10 blocks of 30Trep
+//double g3ave[Nblocksmax][Tblockmax]; //[Tblock] => blocks of Tblock
+double b1; //boundary of particle
+double b2; //boundary of particle
+int alive;
+double Ncuts;
+double MeanLength;
+
+
+/* FUNCTIONS */
+double Diff(double L);
+
+/* MAIN */
+int main(int argc, char* argv[]){
+cout << "Type in 1. xi=tau_b/tau_r (around 360-720) 2. Ageing block [0:20]" <<endl;
+srand(time(NULL));
+
+double xi=atof(argv[1]);
+double c1=(1.0/Trep)*(1/xi); // [1/tau]
+double Tbreak=1/c1; // [tau]
+
+int ageblock=atoi(argv[2]);
+//each block is 10 minutes = 300 relaxation times and about 1 breakage time
+//but I am following the particle only for the first 15 relaxation times (30 seconds)
+
+/*Start loop over attempts*/
+for(int i=0; i<REPMAX; i++){
+cout << "============= " <<endl;
+cout << "== NEW REP == " <<endl;
+cout << "============= " <<endl;
+
+for(int i=0;i<5;i++)REsite_iscut[i]=0;
+
+/*Initialise arrays*/
+for(int t=0;t<Tblockmax;t++) pos[t]=0;
+for(int t=0;t<Tblockmax;t++) Survival[t]=0;
+Ncuts=0;
+
+/* CREATE&OPEN FILES */
+//stringstream namepos;
+//namepos << "Position.Xi"<<xi<<".Block"<<ageblock<<".R"<<i<<".dat";
+//ofstream writepos;
+//writepos.open(namepos.str().c_str());
+//
+stringstream namefi;
+namefi << "SURVIVAL/Survival.Xi"<<xi<<".Block"<<ageblock<<".R"<<i<<".dat";
+ofstream writeSu;
+writeSu.open(namefi.str().c_str());
+//
+stringstream namefieta;
+namefieta << "ETA/eta.Xi"<<xi<<".Block"<<ageblock<<".R"<<i<<".dat";
+ofstream writeETA;
+writeETA.open(namefieta.str().c_str());
+//
+
+int nt=ageblock;
+/* ================================================== */
+/* need to compute the fragments present at this time */
+/* ================================================== */
+    
+//probability of fragmentation is linearly propto time passed
+//divided by the breakage time -- this is true only on average
+double pfrag = (nt*Tblock)/Tbreak;
+
+//how many breaks done? this is the integer part the time passed
+// this is very "quantized" -> solve to include prob that is n+1 (see below)
+int ncuts_old,ncuts;
+if(nt==0)ncuts_old=0;
+else ncuts_old=ncuts;
+ncuts=int(pfrag); //number of cuts
+
+//a more correct way would be to draw from a poisson distribution with
+//mean that is time dependent -- let's see if still ok.
+
+double intpart;
+//this is the fraction of polymers with the new fragment
+// nad is !=0 only if nt*Tblock is not a multiple of Tbreak
+double fraction = modf(pfrag, &intpart);
+
+if(speak)cout <<"1. pfrag=" <<  pfrag << " fraction=" << fraction <<endl;
+
+//now I have ncuts done and "fraction" (The probability of finding a
+// polymer) with ncuts+1 done. Otherwise ncuts. this is to avoid harsh quantisation
+// of cuts versus time
+double pfi=rand()*1.0/RAND_MAX;
+if(pfi < fraction)ncuts=ncuts+1;
+else ncuts=ncuts;
+//bound above by the max number of restriction sites
+ncuts=min(nREsitesmax,ncuts);
+
+    
+if(speak)cout <<"2. ncuts=" << ncuts << endl;
+Ncuts+=ncuts;
+    
+//now need to pick ncuts sites
+// (making sure there aren't two cuts on the same site)
+int nc1=0;
+if(ncuts>0){
+while(1==1){
+int which_cut=rand()%5; //[0:4]
+if(speak)cout << "2.2 Surveying site " << which_cut << " " << REsite_iscut[which_cut] <<  endl;
+    if(REsite_iscut[which_cut]==0){
+    cut[nc1]=REsites[which_cut];
+    REsite_iscut[which_cut]=1;
+    if(speak)cout <<"2.5 chosen " << cut[nc1] << endl;
+    nc1++;
+    }
+if(nc1==ncuts)break;
+}
+}
+else{
+b1=0;
+b2=L0;
+}
+//
+//cin.get();
+//
+for(int n=0;n<ncuts;n++){
+    //cut[n]=REsites[nREs[n]];
+    if(speak)cout << "3. n=" << n <<" cut[n]=" << cut[n] <<endl;
+}
+    
+    
+//now the array cut contains the only cut sites
+//order them to find the pieces
+sort(cut,cut+ncuts); //yes, the arguments are (array, array+length), weird!
+double sortedcuts[ncuts];
+int nc=0;
+for(size_t i = 0; i != ncuts; ++i){
+    sortedcuts[nc]=cut[i];
+    if(speak)cout << "4. nc=" <<nc << " sortedcut[nc]=" << sortedcuts[nc] <<endl;
+    nc++;
+}
+        
+//now the cuts are stored in the array sortedcuts
+//define the fragments lengths
+double lfrag[ncuts+1][2]; //start-end of fragments
+if(ncuts==0){lfrag[0][0]=0;lfrag[0][1]=L0;}
+else{
+for(int i=0;i<ncuts+1;i++){
+    if(i==0){lfrag[i][0]=0;lfrag[i][1]=sortedcuts[i];}
+    else if(i==ncuts){lfrag[i][0]=sortedcuts[i-1];lfrag[i][1]=L0;}
+    else{lfrag[i][0]=sortedcuts[i-1];lfrag[i][1]=sortedcuts[i];}
+    if(speak)cout << "5. i=" << i <<" lfrag[i][0]=" << lfrag[i][0] << " lfrag[i][1]=" << lfrag[i][1]<<endl;
+}
+}
+
+//cin.get();
+ 
+double ml=0;
+for(int i=0;i<ncuts+1;i++){
+ml+=(lfrag[i][1]-lfrag[i][0])/(ncuts+1);
+}
+MeanLength+=ml;
+
+/* ======================== */
+/*   INITIALISE PARTICLES   */
+/* ======================== */
+
+for(int np=0;np<NPARTICLES;np++){
+
+//now that I have the fragments. I can choose them at random for the different particles
+int randfrag=int(rand()*1.0/RAND_MAX*(ncuts+1)); //[0:ncuts]
+if(speak)cout << "5.5 np=" << np << " randfrag=" << randfrag <<endl;
+
+//set my boundaries
+b1=lfrag[randfrag][0];
+b2=lfrag[randfrag][1];
+if(speak)cout << "6. np="<<np << " b1=" << b1 << " b2=" << b2 <<endl;
+
+//place particle
+pos[0]=b1+(rand()*1.0/RAND_MAX)*(b2-b1);
+cout << "=> rep="  << np << " block=" << nt << " time=" << 0 << " b1=" << b1 << " pos=" << pos[0] << " b2=" << b2 << " "<< endl;
+//writepos << np << " " << nt << " " << 0 << " " << b1[np] << " " << pos[np][0] << " " << b2[np] << endl;
+alive=1; //particle returns alive every block
+Survival[0]++;
+    
+//cin.get();
+
+/* ==================================== */
+/*          START BLOCK OF TIME         */
+/* ==================================== */
+for(int tt=1;tt<Tmsd/dt; tt++){
+
+    
+/*DIFFUSION */
+double Lnow=b2-b1;
+int dir=(rand()%2)*2-1; //+1 or -1
+//pos[np][tt]=pos[np][tt-1]+dir*sqrt(2*D0/Lnow*dt); //reptation
+pos[tt]=pos[tt-1]+dir*sqrt(2*D0*dt); //rouse
+    
+if(tt%qdump==0){
+cout << tt*dt << " " << b1 << " " << pos[tt]  << " " << b2 << " "<< i << " " << nt << " " << alive <<  endl;
+//writepos << tt*dt/Trep << " " << b1[np] << " " << pos[np][tt] << " " << b2[np] << " "<< i << " " << nt << endl;
+}
+
+//cout << "=> rep="  << np << " block=" << nt << " time=" << tt << " b1=" << b1[np] << " pos=" << pos[np][tt] << " b2=" << b2[np] << " "<< endl;
+    
+/* Survival function */
+if(pos[tt]<b2 && pos[tt]>b1 && alive==1){
+//cout << np << " " << tt << " sruvived!" <<endl;
+Survival[tt]++;
+}
+
+/* if it hits a boundary the particle dies */
+if(pos[tt]>=b2 || pos[tt]<=b1){alive=0; break;}
+    
+} //close loop over time = Tmsd/dt
+
+}//close loop over particles
+
+int b=nt;
+
+/* ==================== */
+/* PRINT SURVIVAL */
+/* ==================== */
+writeSu <<"#Block " << nt << endl;
+writeSu << "#t/Trep (Trep=L0^2/(pi^2 D0) P_survived(t) " << endl;
+for(int t=0;t<Tmsd/dt;t=t+qdump){
+writeSu << b << " " << t*dt/(Trep/(3.14*3.14)) << " " << Survival[t]*1.0/NPARTICLES << endl; //REPMAX
+}
+writeSu << endl;
+
+/* ==================== */
+/* PRINT VISCOSITY*/
+/* integral of the survival function */
+/* ==================== */
+writeETA <<"#Block " << nt << endl;
+writeETA<< "#xi AgeBlock[300Trep=10min] (AgeBlock+0.5)[30Trep=10min] eta" <<endl;
+double eta0=0;
+double eta=0;
+for(int t=0;t<Tmsd/dt;t++) eta+=Survival[t]*1.0/NPARTICLES*dt; //REPMAX
+cout << b << " "<<eta <<endl;
+writeETA << xi << " " << nt << " " << eta << " " << eta0 << endl;
+
+// xi=Tbreak/Trep
+
+/* ==================== */
+/* PRINT topo */
+/* ==================== */
+/*
+writetopo <<"#t/Trep N_cuts" << endl;
+for(int t=0;t<Ttot/dt;t=t+qdump){
+writetopo << t*dt*1.0/Trep << " " << Ncuts[t]*1.0/REPMAX << endl;
+}
+writetopo<<endl;
+writetopo<<endl;
+for(int rs=0;rs<5;rs++) writetopo<< rs << " " << REsite_iscut[rs] << endl;
+*/
+//writetopo <<"#block N_cuts" << endl;
+//for(int b=0;b<Nblocks;b++){
+//writetopo << b << " " << Ncuts*1.0/REPMAX << endl;
+//}
+
+/* ==================== */
+/* PRINT mean length */
+/* ==================== */
+/*
+writelen <<"#t/Trep Mean_length" << endl;
+for(int t=0;t<Ttot/dt;t=t+qdump){
+writelen << t*dt*1.0/Trep << " " << MeanLength[t]*1.0/REPMAX << endl;
+}
+*/
+//writelen <<"#block Mean_length" << endl;
+//for(int b=0;b<Nblocks;b++){
+//writelen << nt << " " << MeanLength*1.0/REPMAX << endl;
+//}
+
+}//close loop over replicas
+
+
+return 0;
+}
+