{Base-stacking thermo values from SantaLucia, PNAS 1998}

procedure writeZipComments;
begin
writeln(FXfile,'ZIPPER3F.PAS  Thermo values = SantaLucia 1998 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; {SantaLucia 1998 values in 1M NaCl at 25 deg C}
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.00 - 0.0222 *(T-37);
           'c': factor:=1.44 - 0.0224 *(T-37);
           'g': factor:=1.28 - 0.0210 *(T-37);
           't': factor:=0.88 - 0.0204 *(T-37);
           end;{case}
        'c': case Base of
           'a': factor:=1.45 - 0.0227 *(T-37);
           'c': factor:=1.84 - 0.0199 *(T-37);
           'g': factor:=2.17 - 0.0272 *(T-37);
           't': factor:=1.28 - 0.0210 *(T-37);
           end;{case}
        'g': case Base of
           'a': factor:=1.30 - 0.0222 *(T-37);
           'c': factor:=2.24 - 0.0244 *(T-37);
           'g': factor:=1.84 - 0.0199 *(T-37);
           't': factor:=1.44 - 0.0224 *(T-37);
           end;{case}
        't': case Base of
           'a': factor:=0.58 - 0.0213 *(T-37);
           'c': factor:=1.30 - 0.0222 *(T-37);
           'g': factor:=1.45 - 0.0227 *(T-37);
           't': factor:=1.00 - 0.0222 *(T-37);
           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 ~25% of the
simple-FJC energy at 20pN, which was 6pNnm. 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); {open or close basepairs at fork}
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;


