program Tweez13E; {by Steve Smith. Written 2013, modified August 2014
Simulated optical trap unzips DNA hairpins at various temperatures while
recording Force-Position data. Sequences loaded from GENOME folder and
reports written to OUTPUT folder. Basepair energies taken from literature.
Single-stranded stretch energy derived from extensible-FJC model}
uses crt,graph,dos;
const dt = 1e-6;     {simulated iteration time in seconds}
      kB = 1.38066e-23;   {Boltzmann constant J/degK}
      r = 1.5e-6;        {radius of trap bead in meters}
      kTrap = 1e-4;     {trap spring constant in Newtons/meter}
      forceSensitivity = 6e13; {force sensitivity in ADC counts per Newton}
      LLsensitivity = 7e8; {lightLever sensitivity in ADC-units/meter}
      filterConst = 1e-4 ; {pre-amp capacitor time const in seconds}
      cyclePeriod = 1e-3; {program cycle time (per ADC reading) in seconds}
      plotPeriod = 1e-2; {seconds per pixel plotted....also bandwidth}
      minStep = 1.7e-10; {size of single piezo step in meters}
      initDsNum = 30;  {initial dsDNA basepairs for "handles"}
      initSsNum = 0;  {initial ssDNA bases}

type wlcType = array[0..200] of real;
var
f,i,j,k,xMax,yMax,xMid,yMid,oldX,oldY,oldXP,oldYP,sumCount,pixR,pixC,pinLength,
oldPixR,oldPixC,plotColor,forkPos,oldEx,oldF,pullNum,xDir,yDir:integer;

dsNum, ssNum:longInt;

xForce,yForce,dragCoeff,meanThermForce,maxForce,minForce,tension,weight,
lightForceSum,lightLeverSum,fltLightForce,fltLightLever,oldLightForce,
initLightLever,oldLightLever,avLightForce, avLightLever,beadPosX,beadPosY,
pipettePosX,pipettePosY,time,resetTime,filterTime,plotTime,cycleTime,
pixScale,ext,extX,extY,forcePix,initForce,avFsqr,T,absT, viscosity,
avTens,avdTenSqr,pullRate,piezoStep,Lmax,piezoPosX,avXPos,dXsqr,
Tmin,Tmax,Tstep,deltaL,stretchE:real;

dsDNA,ssDNA,molForce,fjcStretch:wlcType; {holds the fractional extensions
of ds and ssDNA for 200 log-spaced forces, plus FJC-stretch energies}
deltaG:array[0..7000] of real; {basepair energies along hairpin}
oldYf,oldYL:array[0..640] of integer;{erases previous force or lever lines}
reading: array[0..101,1..2] of integer; {for force averaging}

holdingConstForce,rampingForce,movingForward,showingForce,showingBeads,
showingReadings,previouslyMovingForward,FXfileOpen: boolean;

ch:char;  dataLine,message,seqName:string;  forceFile, FXfile:text;
    fileIndex,lineIndex:integer;
    baseIndex,startIndex,genomeSize:longInt;
    baseStr:string[1];
    reply:string;
    base,previousBase:char;
    s,fileName:pathStr;
    dird:dirStr; name:nameStr; fileExt:extStr;
    seqdata:text;
    fileopen,endOfDataline,baseFound:boolean;

{$I display4.pas} {procedures to show beads, forces, etc.}
{$I zipper4G.pas} {energy at each basepair 4F=SantaLucia 4G=Huguet}
{$I filestff.pas} {file handling procedures}

procedure initializeMolecule;
begin;
dsNum:=initdsNum;   ssNum:=initssNum;
Lmax := 1.1e-9 * pinLength + 2.0e-11 / kTrap; {sets screen width for F-X graph}
beadPosX:=0; beadPosY:=0;
piezoPosX:=(dsNum+ssNum)*3.4e-10 * 0.90 + 2*r; {initial bead position}
pipettePosX:=piezoPosX; pipettePosY:=0; time:=0; lightForceSum:=0;
oldLightForce:=0;sumCount:=0; filterTime:=0; plotTime:=0;
holdingConstForce:=false; cycleTime:=0; movingForward:=true;
end;

function sumLength(f:integer):real;  {figures length of 2 chains in series}
begin
sumLength:=dsNum*dsDNA[f]+ssNum*ssDNA[f]/2;
{ssDNA length in lookup file is for 2 bases so we must divide ssNum by 2}
end;

procedure computeForces;  {forces on the bead in the trap}
var fraction:real;
begin
extX:=pipettePosX-beadPosX;  {compute extension of molecule}
extY:=pipettePosY-beadPosY;
ext:=sqrt(extX*extX+extY*extY)-2*r;  {end-to-end extension assumes both beads rotate freely}
if ext<0 then tension := 1.0e-2 *ext; {beads bump together}
if ext>=0 then     {find tension for 2 chains in series}
  begin
  while  (sumLength(f+1) < ext) and (f<199) do inc(f);
  while  (sumLength(f) > ext) and (f>0) do dec(f);
  fraction:=(ext-sumLength(f))/(sumLength(f+1)-sumLength(f));
  {the fraction interpolates between two data lines in WLC table}
  tension:=(molForce[f]+ (molForce[f+1]-molForce[f])*fraction)*(T+273.15)/298;
  {T correction because the WLC table was made assuming absT=298 deg}
  deltaL:= (1-fraction)*ssDNA[f]+fraction*ssDNA[f+1];
  stretchE := (1-fraction)*fjcStretch[f]+fraction*fjcStretch[f+1]
  end;
avTens := 0.9999*avTens + 0.0001*tension; {average tension for window display}
avdTenSqr:= 0.999* avdTenSqr+0.001 *sqr(tension-avTens);
xDir:=2*random(2)-1;  {randomize direction of Brownian impact}
yDir:=2*random(2)-1;  {produces either +1 or -1}
{add optical and thermal forces}
xForce:= tension *extX/(ext+2*r) - (kTrap * beadPosX) + xDir*meanThermForce;
yForce:= tension* extY/(ext+2*r) - (kTrap * beadPosY) + yDir*meanThermForce;
end;

procedure adjustTrapBead; {Langevin eqn. for trap bead}
begin
beadPosX:=beadPosX+(xForce/dragCoeff)*dt;
beadPosY:=beadPosY+(yForce/dragCoeff)*dt;
avXPos:=0.999*avXPos + 0.001*beadPosX;   {average position for window display}
dXsqr:= 0.999*dXsqr  + 0.001*sqr(beadPosX-avXPos);
end;

procedure trackLightForce;  {simulates optical force/distance measurements}
begin
lightForceSum:=lightForceSum + (beadPosX * kTrap); {linear trap approximation}
lightLeverSum:=lightLeverSum + pipettePosx;  {measures pipette position}
inc(sumCount);
if time>filterTime then
  begin   {simulate filter capacitor at input of A/D converter}
  fltLightForce:= lightForceSum/sumCount;
  fltLightLever:= lightLeverSum/sumCount;
  fltLightForce:= 0.63*fltLightForce + 0.37*oldLightForce;
  fltLightLever:= 0.63*fltLightLever + 0.37*oldLightLever;
  oldLightForce:=fltLightForce;
  oldLightLever:=fltLightLever;
  lightForceSum:=0; lightLeverSum:=0; sumCount:=0;
  filterTime:=time+filterConst;
  avFsqr:=0.99*avFsqr+0.01*sqr(fltLightForce);  {for window display}
  end;
end;

procedure adjustPipettePosition; {move pipette to pull/relax molecule}
begin
if holdingConstForce then
     begin
     if fltLightForce>weight then piezoPosX:=piezoPosX-minStep;
     if fltLightForce<weight then piezoPosX:=piezoPosX+minStep;
     end
else if rampingForce then
     begin
     if movingForward then piezoPosX:=piezoPosX+piezoStep
                      else piezoPosX:=piezoPosX-piezoStep;
     if fltLightForce>maxForce then movingForward:=false;
     if fltLightForce<minForce then movingForward:=true;
     if previouslyMovingForward <> movingForward then
        begin inc(pullNum); changeColor; end;
     previouslyMovingForward := movingForward;
     end;
pipettePosX:=piezoPosX; {assumes perfect piezo moves pipette instantly
                        Modify here for delay, overshoot, creep}
end;

procedure holdLightReadings;  {get ADC readings and store in array}
begin
if reading[0,1]<100 then inc(reading[0,1]);{ index of current reading }
reading[reading[0,1],1]:=round(fltLightForce * forceSensitivity); { to form a running  average}
reading[reading[0,1],2]:=round(fltLightLever * LLsensitivity); { of recent defl. readings }
end;

procedure averageReadings;
{Average over time between data lines. Further reduce bandwidth}
var i:integer;
    temp1,temp2:longint;
begin
temp1:=0; temp2:=0;
for i:=1 to reading[0,1] do
  begin
  temp1:=temp1+reading[i,1];
  temp2:=temp2+reading[i,2];
  end;
if reading[0,1]>0 then
  begin
  avLightForce:=(temp1 / reading[0,1])/forceSensitivity;
  avLightLever:=temp2 / reading[0,1]/LLsensitivity;
  end;
reading[0,1]:=0;
end;

begin   {main}
clrscr;
assign(forceFile,'xfjc.txt');
readWLC;  {read in F/X tables for ssDNA and dsDNA from WLC model}
pinLength := 7000;
promptSequence;
startGraphics;
showingBeads:=false; showingForce:=true; showingReadings:=false;
rampingForce:=false;
T:=25; {temperture degrees celsius}
initializeMolecule;
repeat
  loadBaseEnergies;  {bp-stacking energies using latest temperature}
  getViscosity;  {water viscosity depends on T}
  repeat
    time:=time+dt;
    computeForces;   {thermal + molecule tension + optical forces}
    adjustTrapBead;  {using Langevin Eqn}
    computeForces;
    trackLightForce;     {filter holds light force between ADC readings}
    zipHairpin(tension); {unfold or refold haipin}
    if time > cycleTime then  {instrument program cycle time}
      begin
      adjustPipettePosition; {for const-force or ramp-force}
      holdLightReadings; {gets and holds ADC readings}
      if time>plotTime then
         begin
         averageReadings;
         if showingForce then plotForce;
         if showingBeads then showBeads;
         if showingReadings then displayReadings;
         if FXfileOpen then  writeln(FXfile,pullNum,avLightForce*1e12:20:3,
                             (piezoPosX-(2*r))*1e6:15:3,forkPos:15);
         plotTime:=time+plotPeriod;
         end;  {plotTime}
      cycleTime:=time+cyclePeriod;
      end;  {cycleTime}
    if keypressed then consultMenu;
  until pullNum>2;  {makes 1 pull/relax cycle}
  T:=T+Tstep;
  if T <= Tmax then incrementFile;  {opens new file for each new temperature}
until (ch='q') or (T > Tmax);
closeGraph;
close(FXfile);
end.
