C This program simulate the molecule rupture process 
C in single molecular pulling experiments. 
C Model
C   Interaction between ligand-receptor: LJ-12 potential(u0,sigma)
C Dimensionless terms 
C   time t=DT/sigma2
C   distance z=Z/sigma
C   cc=sigma*loadingrate/D
C Reduced Units 
C   energy u0 in 1KT=4E-21Joules.
C   distance sigma in 1nm
C   force in 4pN
C Parameters
C   beta 1/(KT)
C   d diffusion coefficient <0.06  
C   kc spring constant
C   idum RNG seed
C   lrate loading rate
C   trun- total number of steps for simulation
C   t(dt) dimensionless time(timestep)
C   z(dz) dimensionless position(displacement in each step)
C   dzr dimensionless random brownian displacement
C   zmin dimensionless position where PLJ is the min
C   force(flj,fs) force(LJ force, spring provided force) on the particle
C   pot(plj,ps) potential(LJ potential, spring provided potential) on the particle


      program brown
      implicit none
      integer size
      parameter (size=20000004)
      double precision fz(size),fljz(size),fsz(size)
      double precision pz(size),pljz(size),psz(size)
      double precision z(size),t(size) 
      double precision dz(size),dt(size) 
      double precision dt0,dz0,z0,t0
      integer idum,l,i,j,trun 
c     double precision lrate,d
      double precision beta,cc,kc
      double precision u0, sigma,zmin
      double precision ran1,gasdev,dzr
      double precision force,absf,flj,fs,pot,plj,ps

      open(unit=1,file='table1')
      open(unit=2,file='table2')
      open(unit=3,file='steps')
      
c      write(1,*) 'time','disp','F','Flj','Fs'
c      write(2,*) 'disp','P','Plj','Ps'



      trun=20000000
      beta=1.0D0
      idum=-1
      kc=1.0D0
      u0=5.0D0
      sigma=1.0D-1 
      zmin=1.122
c     lrate=1.0D-2
c     d=1.0D-6
Cc    cc=sigma*lrate/d 
      cc=1.0D-2

Cc Initialize
      do j=1,trun
       dz(j)=0.0D0
       z(j)=zmin
       dt(j)=0.0D0
       t(j)=0.0D0
       
       fz(j)=0.0D0
       fljz(j)=0.0D0
       fsz(j)=0.0D0
       pz(j)=0.0D0
       pljz(j)=0.0D0
       psz(j)=0.0D0
      enddo

     
Cc Move the particle
   
      do i=1,trun

         z0=z(i)
         t0=t(i)

       call potent(u0,sigma,kc,cc,t0,z0,force,flj,fs,pot,plj,ps) 
	  
         pz(i)=pot
         psz(i)=ps
         pljz(i)=plj
         fz(i)=force
         fsz(i)=fs
         fljz(i)=flj
       


Cc sig= squareroot(2Ddt)
c      sig=sqrt(d*dt0*2.0)   
c      adj=sig*gasdev(idum)
         absf=abs(force)
         dt0=0.01*0.122/beta/absf/sigma
         
         dt(i+1)=dt0
c      t(2)=t(1)+dt(2)
         t(i+1)=t(i)+dt0
         dzr=sqrt(2.0*dt0)*gasdev(idum)
         dz0=(beta*sigma*force*dt0)+dzr 
         dz(i+1)=dz0
c      z(2)=t(1)+dz(2)
         z(i+1)=z0+dz0
      enddo

Cc Regularize units in output  
    
      do l=1,trun,1000     
C       write(1,20) z(l),fz(l),fljz(l),fsz(l)
C       write(2,20) z(l),pz(l),pljz(l),psz(l)
       write(3,20) dt(l),t(l),z(l),fsz(l)
      enddo

20    format (f12.5,3(8x,f12.5))
      end 



      subroutine potent(u0,sigma,kc,cc,tt,zz,force,flj,fs,pot,plj,ps)
      double precision zz,zmin,zz0,zz2,zz6,u0,sigma,sigma2,kc,tt,cc,dl
      double precision force,flj,fs,pot,plj,ps,beta
      beta=1.
      zmin=1.122
      sigma2=sigma*sigma 
c    cc=sigma*v/d  
      dl=cc*tt-(zz-zmin)
      zz0=1/zz
      zz2=zz0*zz0
      zz6=zz2*zz2*zz2

      plj=4*u0*zz6*(zz6-1)  
      ps=0.5*kc*sigma2*dl**2 
      pot=plj+ps

      flj= 24*u0*zz6*zz0*(2*zz6-1)/sigma 
      fs= kc*sigma*dl
      force=flj+fs
      return 
      end
      


      function gasdev(idum)
      integer idum
      real gasdev
      integer iset
      real fac,gset,rsq,a1,a2,ran1
      save iset,gset
      data iset/0/
      if (idum.lt.0) iset=0
      if (iset.eq.0) then
1       a1=2.*ran1(idum)-1
        a2=2.*ran1(idum)-1
        rsq=a1**2+a2**2
        if (rsq.ge.1..or.rsq.eq.0.)go to 1
        fac=sqrt(-2.*log(rsq)/rsq)
        gset=a1*fac
        gasdev=a2*fac
        iset=1
      else   
        gasdev=gset
        iset=0
      endif  
      return
      end



      function ran1(idum)
C  Random generator from NR, pp.270-271
C  Valid for at least 10^8 calls
C  First call with idum a negative integer
      integer idum,IA,IM,IQ,IR,NTAB,NDIV
      real ran1,AM,EPS,RNMX
      parameter (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     &    NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      integer j,k,iv(NTAB),iy
      save iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
          idum=max(-idum,1)
          do j=NTAB+8,1,-1
             k=idum/IQ
             idum=IA*(idum-k*IQ)-IR*k
             if (idum.lt.0) idum=idum+IM
             if (j.le.NTAB) iv(j)=idum
          enddo
          iy=iv(1)
      endif 
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      end
