program WLC4;  { Computes and saves dsDNA extensions for 200 log-spaced
forces using power-series WLC from Brochard & Croquette, Biophysical J, 1998.
Also computes extensions for ssDNA using extensible-FJC model and values from
A. Bosco et al. Nucl.Acids Res. 2014.  Assumes [Na]= 1 Molar.  The FJC energy is
integral of F*dx to stretch 2 bases from force=zero to force=F}
{$N+}
uses crt,graph;
const
  kT = 4.114e-21; {thermal energy in joules at 25 deg.C}
  stretchMod1= 1.2e-9;{stretch modulus of dsDNA in Newtons}
  stretchMod2= 0.8e-9;{stretch modulus of ssDNA in Newtons}
  P1 = 5e-8; {persistence length of dsDNA in meters}
  b = 1.3e-9; {Kuhn length of ssDNA in meters}
  LCds = 0.34e-9;   {rise in meters per basepair of dsDNA}
  LCss = 1.14e-9;   {contour length of 2 bases in series}

var snert,tension,f,f2,x,lastX,dX,L,y,P,stretchEnergy:real;
    ch:char; i,j,n, xMax,xMid,yMax,yMid:integer;
    s:array[1..7] of real;
    outFile:text;
    message:string;
begin
clrscr;
writeln('what to call the output file? *.txt');
readln(message);
assign(outFile,message);
rewrite(outFile);
writeln(outFile,'generated by WLC4.PAS    T=25C    salt=1M NaCl');
writeln(outFile,'P(ds)=',P1*1e9:1:0,'nm   Lc(ds)=',LCds*1e9:2:2,'   S(ds)=',stretchMod1*1e12:0:0,'pN');
writeln(outFile,'b(ss)=',b*1e9:2:2,'   Lc(ss)=',LCss*1e9:2:2,'   S(ss)=',stretchMod2*1e12:0:0,'pN');
writeln(outFile,'dsDNA':10, 'ssDNA':30,'force(N)':33,'FJC_energy':33);
i:=detect;
initgraph(i,j,'C:\PASCAL\BGI');
xMax:=getMaxX; xMid :=xMax div 2; yMax:=getMaxY; yMid :=yMax div 2;
line(xMid,0,xMid,yMax-30);
rectangle(0,0,xMax,yMax);
for i:=1 to 7 do
   begin   j:= yMax-i*yMax div 7;
   line(xMid-5,j,xMid+5,j);
   str(16-i,message);
   outTextXY(xMid+7,j-3,message);
   end;
setTextStyle(1,0,1);
outTextXY(yMid-24,yMax-20,'extension / dsL');
setTextStyle(1,1,1);
outTextXY(xMid+20,yMid,'-log force (N)');
tension := 1e-15; {Newtons}
j:=1;  setColor(red);

repeat
P:=P1; snert:=kT/(4*P);{for dsDNA}
L:=(1+tension/stretchMod1);
x:=0.5;  n:=1;
repeat
      s[1]:=x/L;
      for i:=2 to 7 do s[i]:=s[1]*s[i-1];
      f := snert*(1/sqr(1-(x/L))+4*(x/L)-1-2.0656912*s[2]-10.949672*s[3]
      +64.29988*s[4]-155.50428*s[5]+157.99776*s[6]-56.70872*s[7]);
      f2 :=(2*snert/L)*(1/((1-(x/L))*(1-(x/L))*(1-(x/L)))+2);
      if (tension-f)/(n*f2) < (L-x) then x:=x+(tension-f)/(n*f2) else inc(n);
until abs(f-tension)<0.01*tension;
y:= ln(tension)/ln(10)+16;
putPixel(round(xMid*x), yMax-round(y*(yMax div 7)),green);
write(outFile, x*LCds:30);     {fractional extension times contour length}

snert:=tension*b/kT; {for ssDNA}
L:=(1+tension/stretchMod2);
if tension<1e-10
   then x:=((exp(snert)+exp(-snert))/(exp(snert)-exp(-snert))-(1/snert))*L
   else x:= (1 - 1/snert)*L;
y:= ln(tension)/ln(10)+16;
putPixel(round(xMid*x), yMax-round(y*(yMax div 7)),cyan);
write(outFile, x*LCss:30);  {extension for 2 bases in meters}
write(outFile,tension:30);  {in Newtons}
dX := LCss * (x-lastX);
stretchEnergy:=stretchEnergy + (tension * dX);
lastX:=x;
writeln(outFile,stretchEnergy:30); {in Joules}
tension:=tension*1.059253725;    {makes 200 entries}
until (tension>1e-10);
close(outFile);
repeat until keypressed;
end.