{Base-stacking thermo values from Huguet et al PNAS 2010 }

procedure writeZipComments;
begin
writeln(FXfile,'ZIPPER3G.PAS  Thermo values = Huguet 2010 for 1M NaCl');
end;

procedure  getBaseType;
begin
if endOfDataLine then
   begin readln(seqData,dataLine); lineIndex :=0; endOfDataLine:=false; end;
inc(lineIndex);
baseStr:=copy(dataLine,lineIndex,1);
baseFound:=false;   {skips line numbers and spaces}
if baseStr ='a'then begin baseFound:=true;base:='a';end;
if baseStr ='c'then begin baseFound:=true;base:='c';end;
if baseStr ='g'then begin baseFound:=true;base:='g';end;
if baseStr ='t'then begin baseFound:=true;base:='t';end;
if lineIndex >= 80 then endOfDataLine:=true;
end;

procedure loadBaseEnergies;
{values in 1M NaCl at 25 degees from J.M.Huguet, PNAS 2010}
var factor,energy,eqForce:real;
begin
baseIndex:=0; endOfDataLine:=true; previousBase:='a';
reset(seqData);
repeat
  getBaseType;
  if baseFound then
     begin
     inc(baseIndex);
     if baseIndex>startIndex then
       begin
       case previousBase of
        'a': case Base of
           'a': factor:=1.23 - 0.02028 *(T-25);
           'c': factor:=1.49 - 0.01446 *(T-25);
           'g': factor:=1.36 - 0.01289 *(T-25);
           't': factor:=1.17 - 0.01162 *(T-25);
           end;{case}
        'c': case Base of
           'a': factor:=1.66 - 0.02448 *(T-25);
           'c': factor:=1.93 - 0.02230 *(T-25);
           'g': factor:=2.37 - 0.02443 *(T-25);
           't': factor:=1.36 - 0.01289 *(T-25);
           end;{case}
        'g': case Base of
           'a': factor:=1.47 - 0.02246 *(T-25);
           'c': factor:=2.36 - 0.02596 *(T-25);
           'g': factor:=1.93 - 0.02230 *(T-25);
           't': factor:=1.49 - 0.01446 *(T-25);
           end;{case}
        't': case Base of
           'a': factor:=0.84 - 0.02506 *(T-25);
           'c': factor:=1.47 - 0.02246 *(T-25);
           'g': factor:=1.66 - 0.02448 *(T-25);
           't': factor:=1.23 - 0.02028 *(T-25);
           end;{previousBase case}
        end;{outer case}
       energy:=factor*6.95e-21; {convert kcal/mole to joules}
       deltaG[baseIndex-startIndex]:=energy;
       previousBase:=base;
       end;{baseIndex>startIndex}
     end; {baseFound}
until (baseIndex-startIndex)>=pinLength;
end;

function hump(temperature:real):real;
{Simulates the condensation plateau on ssDNA elasticity curves.
According to Fig7B (ver.18 main text) the hump represents about 25% of the
simple-FJC energy at 20pN, which was 6pN*nm. So hump = 0.25 * 6 = 1.5 pNnm.
I estimate the hump grows linearly from zero @ 35 deg to 1.5 pNnm @ 5 deg.}
var proxy:real;
begin
proxy:=1.5e-21 * (35-T)/30;
if proxy<0 then proxy:=0;
hump:=proxy;
end;

procedure zipHairpin(tension:real); {adjust hairpin length}
var dG,B:real;
begin
if forkPos<pinlength then  {try unzipping farther}
   begin
   dG:= deltaG[forkPos+1] +stretchE {+hump(T)} -tension*deltaL;
   B:=exp(-dG/2/(kB*absT));
   if B/(1+B)>random then inc(forkPos);
   end;
if forkPos>0 then   {try rezipping toward base}
   begin
   dG:= -deltaG[forkPos] -stretchE {-hump(T)} +tension*deltaL;
   B:=exp(-dG/2/(kB*absT));
   if B/(1+B)>random then dec(forkPos);
   end;
ssNum:=forkPos*2;
end;


