made more efficient for long sim

NonEq_Living_Polymers_fusion.v2.c++

@@ -7,10 +7,8 @@ SET c1 instead of xi -> easier for this purpose!
 here c1=phi is the number of fusions per reptation time
 ######
 REMOVE ALL ARRAYS BC RELAXATION TIME DIVERGES SO VERY LONG SIMULATIONS ARE NEEDED
+Since there is no breakage when L>100 the relaxation is longer than the total simulation time so the sim is stopped and the longest relax time recorded.
 #######
-1. The reptation of a polymer within a tube
-2. breakage
-=> How does the viscosity depend on the mean length?
 */
 
 #include<iostream>
@@ -24,16 +22,16 @@ REMOVE ALL ARRAYS BC RELAXATION TIME DIVERGES SO VERY LONG SIMULATIONS ARE NEEDE
 using namespace std;
 
 /* CONSTANTS */
-const int REPMAX=100;
+const int REPMAX=300;
 const double dL=0.01; //
 const double dt=0.01; //integration time
 const double tau=1.0; //unit of time (100dt)
-const double L0=0.5; //L0 sets the units of length (100dL)
-const double D0=0.01; //monomer diffusion. In units of L0^2/dt
+const double L0=1.0; //L0 sets the units of length (100dL)
+const double D0=0.1; //monomer diffusion. In units of L0^2/dt
                         //if L0=1; set D0=(1/dt) to get L0^2 monomer diffusion every
                         //tau=(1/dt) steps. [D0]=0.01[L0^2/dt]=1[L0^2/tau]
 const int bin=10;//for histogram. Time units
-const int Tmaxmax=10000000; // this is bc it wants something not variable for arrays
+const int Tmaxmax=100000000; // this is bc it wants something not variable for arrays
 const int Tmaxblock=100000; // this is bc it wants something not variable for arrays
 
 /* ARRAYS */
@@ -56,10 +54,9 @@ rand();
 //cout << "Trep=" << Trep <<endl;
 
 double phi=atof(argv[1]); //phi=Trep/tau_f number of fusions in one reptation time
-double Trep=pow(L0,3.0)/D0; //reptation time (units dt)
+double Trep=pow(L0,3.0)/(D0*3.14*3.14); //reptation time (units dt)
 double Tbreak=Trep*1.0/phi;
 int Tmax=int(100000*Trep); //total experiment time (units dt)
-int Tobs=int(1*Trep); //observation time to make  (units dt)
 double c1=phi*1.0/Trep; //c1 in this case does not depend on length L0 and is 1/tau_f
 
 //ageing block (in units of Tobs=Trep)
@@ -70,10 +67,11 @@ for(int i=0;i<REPMAX;i++)RelaxTime[i]=0;
 //this is the time at which tube has relaxed; it is one number per replica
 
 /* CREATE/OPEN FILES */
-stringstream namepos;
+/*stringstream namepos;
 namepos << "Position_phi"<<phi<<"_L"<<L0<<"_Block"<<block_number <<".dat";
 ofstream writepos;
 writepos.open(namepos.str().c_str());
+*/
 //
 stringstream namefi;
 namefi << "Survival_phi"<<phi<<"_L"<<L0<<"_Block"<<block_number <<".dat";
@@ -115,7 +113,7 @@ pos_right0=-Lmean/2.*log(rand()*1.0/RAND_MAX); //sampled from expoential ditro
 
 /* Tube segment starts at 0 */
 cout << 0 << " " << -pos_left0 <<" " << pos_right0<< " "<< i << endl;
-writepos << 0 << " " << -pos_left0 <<" " << pos_right0<< " "<< i << endl;
+//writepos << 0 << " " << -pos_left0 <<" " << pos_right0<< " "<< i << endl;
 
 //cin.get();
 
@@ -126,7 +124,7 @@ double b1,b2,L;
 /* =============== */
 /* Start time. The particle is at 0.*/
 /* Pos_left is actually to negative value*/
-for(int t=1;t<Tmax/dt; t++){
+for(int t=1;t<=Tmax/dt; t++){
 
 if(t==1){
 L=pos_right0+pos_left0;
@@ -204,7 +202,7 @@ pos_left=b1-dir*sqrt(2*Diff(L)*dt);
 pos_right=b2+dir*sqrt(2*Diff(L)*dt);
 
 //cout << t*dt << " " << -pos_left << " " << pos_right << " "<< i << endl;
-writepos << t*dt << " " << -pos_left << " " << pos_right << " "<< i << endl;
+//writepos << t*dt << " " << -pos_left << " " << pos_right << " "<< i << endl;
 //if(t%10000==0)cin.get();
 
 /* Survival function */
@@ -213,12 +211,16 @@ if(pos_left>0 && pos_right>0){
 }
 
 /* if one of them hit 0 break */
-if(pos_left<0 || pos_right<0){
+if(pos_left<0 || pos_right<0 || t==Tmax/dt){
 RelaxTime[i]=t;
 cout << "relaxed @ " << t <<endl;
 break;
 }
 
+if(L>100){
+RelaxTime[i]=Tmax/dt;
+break;
+}
 //writel << i << "  " << t*dt << " " << L <<endl;
 //InstLength[t]+=L;
 
@@ -227,8 +229,8 @@ break;
 //writel <<endl;
 //writel <<endl;
 
-writepos << endl;
-writepos << endl;
+//writepos << endl;
+//writepos << endl;
 
 }//close loop over attempts
 
@@ -248,7 +250,7 @@ double mu_t=0;
     for(int i=0;i<REPMAX;i++){
     if(t<RelaxTime[i])mu_t++;
     }
-    if(mu_t>0)writeSu << t*dt*1.0/Trep << " " << mu_t*1.0/REPMAX<< endl;
+    if(mu_t>0 && t%10==0)writeSu << t*dt*1.0/Trep << " " << mu_t*1.0/REPMAX<< endl;
     eta+=dt*mu_t*1.0/REPMAX;
 }