c============================================================================
c
c  // Dynamic Monte Carlo Code (DMCC) of electrons in liguid water //
c  // version: chemical; single spar //
c
c                            developed by Takeshi Kai and Tomohiro Toigawa
c
c============================================================================
      implicit real*8(a-h,o-z)
     
      include 'param.inc'
 
      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id
      common / main2 / dif1(nox), dif2(nox)
      common / main3 / dx(3,3), dy(3,3), dz(3,3), dr(3,3)

      common / analyze1 / xdis(-nrg:nrg), xps(3,10,-nrg:nrg), nxpnt
      common / analyze2 / itime, ntime, dt, t0, ntrack
      common / analyze3 / prb(3,nit)
      common / analyze4 / reacd1(3), reacp1(3), ipflg(3), iflg
      common / analyze5 / rdis(0:nrg), rps(3,10,0:nrg), irpoint
      common / analyze6 / reac, scav, escp, DNAe, DNAo, DNAc
      common / analyze7 / iflge, iflgo, iflgc

      common / randam / idum

      character*12 fname1, fname2

      real ran2


* (distribution: 1:delta, 2:gauss, 3:exponential, 4:dmcc_phys)
      id  = 2   

* (radical scavenging on/off = 1/0)
      ird = 0   

* (geometry on/off = 1/0)
      igd = 1


* (input)
      ntrack    =  500000                
      itime     =  1000000        ! 100 (ns)
      ntime     =  100      

      dt        =  100.d-15       ! 100 (fs)
      t0        =    1.d-12       !   1 (ps)
      rs        =    0.3d-9 


      idum    =         15
      fname1  =  'rd12_015.out'
      fname2  =  'rd13_015.out'


          call setup1             ! [analyze]


      do 1000 itrack = 1, ntrack 

          call setup2             ! [analyze]


        do 2000 it = 1, itime


          t = t0 + dble( it - 1 ) * dt


        if( it .eq. 1 ) then
          call reaction0( it )
              if( iflg .eq. 1 ) then
                goto 1001
              else
                goto 2000
              endif 
        endif


          if( ird .eq. 1 )then
          call scavenging( it )
              if( iflg .eq. 2 ) goto 1001
          endif


          if( igd .eq. 1 )then
          call geometry( it )
              if( iflge .eq. 1 .and. iflgo .eq. 0 ) then
                  iflg = 3
                  goto 1002
              endif
              if( iflge .eq. 0 .and. iflgo .eq. 1 ) then
                  iflg = 4
                  goto 1002
              endif
              if( iflge .eq. 1 .and. iflgo .eq. 1 ) then
                  iflg = 5 
                  goto 1001
              endif
          endif
 1002   continue

          call reaction1( it )
              if( iflg .eq. 1 ) goto 1001


          call diffusion( it )
            do j = 2, 3
              sdist = reacd1(j) + rs
              if( dr(1,j) .gt. sdist ) ipflg(j) = 0
            enddo


          call analyze( it )


 2000   continue ! (time loop)


 1001   continue
        if(     iflg .eq. 0 ) then
          escp = escp + 1.d0 / dble(ntrack)
        elseif( iflg .eq. 1 ) then
          reac = reac + 1.d0 / dble(ntrack)
        elseif( iflg .eq. 2 ) then
          scav = scav + 1.d0 / dble(ntrack)
        elseif( iflg .eq. 3 ) then
          DNAe = DNAe + 1.d0 / dble(ntrack)
        elseif( iflg .eq. 4 ) then
          DNAo = DNAo + 1.d0 / dble(ntrack)
        elseif( iflg .eq. 5 ) then
          DNAc = DNAc + 1.d0 / dble(ntrack)
        endif


 1000 continue ! (track loop)


      call output( fname1, fname2 )             ! [analyze]


      end 


c================================================
      subroutine geometry( it )
c===============================================
      implicit real*8(a-h,o-z)

      include 'param.inc'

      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id
      common / main2 / dif1(nox), dif2(nox)
      common / main3 / dx(3,3), dy(3,3), dz(3,3), dr(3,3)

      common / analyze1 / xdis(-nrg:nrg), xps(3,10,-nrg:nrg), nxpnt
      common / analyze2 / itime, ntime, dt, t0, ntrack
      common / analyze3 / prb(3,nit)
      common / analyze4 / reacd1(3), reacp1(3), ipflg(3), iflg
      common / analyze5 / rdis(0:nrg), rps(3,10,0:nrg), irpoint
      common / analyze7 / iflge, iflgo, iflgc

      common / randam / idum

      real ran2


c===================================================
c      rate = 1.9d+8   ! H + AMP (Buxton,p649)
c      rate = 2.9d+8   ! H + TMP (Buxton,p673)
c      rate = 2.1d+10  ! e + 5-Bromouracil (Buxton,p584)
c===================================================


      rate_eA = 3.8d+9   ! e + AMP (Buxton,p574)
      rate_eC = 6.8d+9   ! e + CMP (Buxton,p591) 
      rate_eG = 1.5d+9   ! e + GMP (Buxton,p606)
      rate_eT = 1.5d+9   ! e + TMP (Buxton,p631)


      ratioA = 1140.d0 / 4731.d0
      ratioC = 1284.d0 / 4731.d0
      ratioG = 1243.d0 / 4731.d0
      ratioT = 1064.d0 / 4731.d0

 
      cons1b = 3.02E-03    ! [mol/L] 
      consA = cons1b * ratioA  
      consC = cons1b * ratioC  
      consG = cons1b * ratioG  
      consT = cons1b * ratioT  


      alpha_e = consA * rate_eA + consC * rate_eC + 
     .          consG * rate_eG + consT * rate_eT 


      rdne = dsqrt( x(1)**2 + y(1)**2 ) / xnm
      dna1 =  1.7d0 * xnm 
      dna2 = -1.7d0 * xnm 


      prbs = 1.d0 - dexp( -alpha_e * dt )   ! from 0 to 1
      rans = dble(ran2(idum))

     
      if( iflge .eq. 1 ) goto 100 
      if( rdne .le. 1.d0 .and. 
     .    z(1) .le. dna1 .and. 
     .    z(1) .ge. dna2 .and. 
     .    rans .lt. prbs ) then
        do i = 1, 1
          do nt = it, itime
            prb(i,nt) = prb(i,nt) - 1.d0 / dble(ntrack)
            if(prb(i,nt) .lt. 0.d0 ) then
               prb(i,nt) = 0.d0
            endif
          enddo
        enddo
        iflge = 1
      endif
  100 continue



      rate_OHA = 4.1d+9   ! OH + AMP (Buxton,p631)
      rate_OHC = 4.7d+9   ! OH + CMP (Buxton,p710)
      rate_OHG = 4.7d+9   ! OH + GMP (Buxton,p727)
      rate_OHT = 5.2d+9   ! OH + TMP (Buxton,p754)


      alpha_o = consA * rate_OHA + consC * rate_OHC +
     .          consG * rate_OHG + consT * rate_OHT

c      time0 = 1.d-9
c      write(*,300)alpha_e*time0, alpha_o*time0
c  300 format(10(1pe15.7))
c      stop

      rdno = dsqrt( x(2)**2 + y(2)**2 ) / xnm


      prbs = 1.d0 - dexp( -alpha_o * dt )   ! from 0 to 1
      rans = dble(ran2(idum))


      if( iflgo .eq. 1 ) goto 200
      if( rdno .le. 1.d0 .and.
     .    z(2) .le. dna1 .and.
     .    z(2) .ge. dna2 .and.
     .    rans .lt. prbs ) then
        do i = 2, 2
          do nt = it, itime
            prb(i,nt) = prb(i,nt) - 1.d0 / dble(ntrack)
            if(prb(i,nt) .lt. 0.d0 ) then
               prb(i,nt) = 0.d0
            endif
          enddo
        enddo
        iflgo = 1
      endif
  200 continue


      end

c================================================
      subroutine scavenging( it )
c===============================================
      implicit real*8(a-h,o-z)

      include 'param.inc'

      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id
      common / main2 / dif1(nox), dif2(nox)
      common / main3 / dx(3,3), dy(3,3), dz(3,3), dr(3,3)

      common / analyze1 / xdis(-nrg:nrg), xps(3,10,-nrg:nrg), nxpnt
      common / analyze2 / itime, ntime, dt, t0, ntrack
      common / analyze3 / prb(3,nit)
      common / analyze4 / reacd1(3), reacp1(3), ipflg(3), iflg
      common / analyze5 / rdis(0:nrg), rps(3,10,0:nrg), irpoint

      common / randam / idum

      real ran2


      rate = 2.1d+10
      dens = 1.d-3
      prbs = dexp( -dens * rate * dt )
      rans = dble(ran2(idum))


      if( rans .gt. prbs ) then
        do i = 1, 1
          do nt = it, itime
            prb(i,nt) = prb(i,nt) - 1.d0 / dble(ntrack)
            if(prb(i,nt) .lt. 0.d0 ) then
               prb(i,nt) = 0.d0
            endif
          enddo
        enddo
        iflg = 2
      endif


      end


c================================================
      subroutine reaction0( it )
c===============================================
      implicit real*8(a-h,o-z)

      include 'param.inc'

      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id
      common / main2 / dif1(nox), dif2(nox)
      common / main3 / dx(3,3), dy(3,3), dz(3,3), dr(3,3)

      common / analyze1 / xdis(-nrg:nrg), xps(3,10,-nrg:nrg), nxpnt
      common / analyze2 / itime, ntime, dt, t0, ntrack
      common / analyze3 / prb(3,nit)
      common / analyze4 / reacd1(3), reacp1(3), ipflg(3), iflg
      common / analyze5 / rdis(0:nrg), rps(3,10,0:nrg), irpoint 

      common / randam / idum

      real ran2


      do i = 1, 3
      do j = 1, 3
        dx(i,j) =  ( x(i) - x(j) )
        dy(i,j) =  ( y(i) - y(j) )
        dz(i,j) =  ( z(i) - z(j) )
        dr(i,j) =  dsqrt( dx(i,j)**2 + dy(i,j)**2 + dz(i,j)**2 )
      enddo
      enddo


      if( dr(1,3) .le. 0.75d-9 ) then
        do i = 1, 3
          do nt = it, itime
            prb(i,nt) = prb(i,nt) - 1.d0 / dble(ntrack)
          enddo
        enddo
        iflg = 1
      endif


      end


c================================================
      subroutine reaction1( it )
c===============================================
      implicit real*8(a-h,o-z)

      include 'param.inc'

      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id
      common / main2 / dif1(nox), dif2(nox)
      common / main3 / dx(3,3), dy(3,3), dz(3,3), dr(3,3)

      common / analyze1 / xdis(-nrg:nrg), xps(3,10,-nrg:nrg), nxpnt
      common / analyze2 / itime, ntime, dt, t0, ntrack
      common / analyze3 / prb(3,nit)
      common / analyze4 / reacd1(3), reacp1(3), ipflg(3), iflg
      common / analyze5 / rdis(0:nrg), rps(3,10,0:nrg), irpoint

      common / randam / idum

      real ran2


      do i = 1, 3
      do j = 1, 3
        dx(i,j) =  ( x(i) - x(j) )
        dy(i,j) =  ( y(i) - y(j) )
        dz(i,j) =  ( z(i) - z(j) )
        dr(i,j) =  dsqrt( dx(i,j)**2 + dy(i,j)**2 + dz(i,j)**2 )
      enddo
      enddo


      do i = 1, 1
      do j = 2, 3
        rnd = dble( ran2(idum) )
        if( dr(i,j) .le. reacd1(j) ) then

          if( rnd .le. reacp1(j) .and. ipflg(j) .eq. 0 ) then
            do nt = it, itime
              prb(i,nt) = prb(i,nt) - 1.d0/dble(ntrack)
              prb(j,nt) = prb(j,nt) - 1.d0/dble(ntrack)
            enddo

            ipflg(j) = 1
            iflg = 1
            goto 1000

          else
            ipflg(j) = 1
          endif

        endif
      enddo
      enddo
 1000 continue


      end


c================================================
      subroutine diffusion( it )
c===============================================
      implicit real*8(a-h,o-z)

      include 'param.inc'

      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id
      common / main2 / dif1(nox), dif2(nox)
      common / main3 / dx(3,3), dy(3,3), dz(3,3), dr(3,3)

      common / analyze2 / itime, ntime, dt, t0, ntrack

      dimension r(3)

      t  =  t0 + dble( it - 1 ) * dt

      call eps_rela(t,er)

        x1  =  chge**2
        x2  =  4.d0 * pi * er * eps0
        x3  =  bltz * temp
        rc  =  x1 / x2 / x3     


      dna1 =  1.7d0 * xnm
      dna2 = -1.7d0 * xnm

      do i = 1, 2
        r(i) = dsqrt( x(i)**2 + y(i)**2 ) / xnm

          if( r(i) .le. 1.d0 .and.
     .        z(i) .le. dna1 .and.
     .        z(i) .ge. dna2 ) then 
    
              dif1(i) = dif1(i) / 1.d+2
              dif2(i) = dsqrt( 4.d0 * dif1(i) * dt )

          endif

      enddo


      do i = 1, 3
        dxd  =  dif2(i) * sample( dble( ran2(idum) ) )
        dyd  =  dif2(i) * sample( dble( ran2(idum) ) )
        dzd  =  dif2(i) * sample( dble( ran2(idum) ) )

        dxc = 0.d0
        dyc = 0.d0
        dzc = 0.d0

      do j = 1, 3
        if( i .ne. j ) then
 
          if( is(i) .eq. 2 ) goto 1000 
          if( is(j) .eq. 2 ) goto 1000 

          drc  =  rc  * dif1(i) * dt  / dr(i,j)**2
          dxc  =  dxc + drc * dx(i,j) / dr(i,j)
          dyc  =  dyc + drc * dy(i,j) / dr(i,j)
          dzc  =  dzc + drc * dz(i,j) / dr(i,j)

        endif
 1000 continue
      enddo

        x(i)  =  x(i)  +  dxd - dxc
        y(i)  =  y(i)  +  dyd - dyc
        z(i)  =  z(i)  +  dzd - dzc

c        x(i)  =  x(i)  +  dxd 
c        y(i)  =  y(i)  +  dyd 
c        z(i)  =  z(i)  +  dzd
 
      enddo


      dif1(1)  =  4.9d-9    ! (m2/s)
      dif1(2)  =  2.2d-9    ! (m2/s)
      dif1(3)  =  9.5d-9    ! (m2/s)


      do i = 1, 3
        dif2(i)  =  dsqrt( 4.d0 * dif1(i) * dt )
      enddo


  100 format(8(1pe15.7))
      end 


c================================================
      subroutine eps_rela(t,er)
c===============================================
      implicit real*8(a-h,o-z)
      parameter ( pi  =  3.1415926535897932385d0 )


      eps1 = 74.9d0
      eps2 = 1.67d0
      eps8 = 2.00d0

      tau2 = 248.d-15

      as   = 31.5d+24 * ( 2.d0 * pi )**2
      ws   = 5.30d+12 * ( 2.d0 * pi )
      gs   = 5.35d+12 * ( 2.d0 * pi )

      al   = 108.d+24 * ( 2.d0 * pi )**2
      wl   = 14.7d+12 * ( 2.d0 * pi )
      gl   = 8.08d+12 * ( 2.d0 * pi )

      temp = 293.d0
      tau1 = 1.08d0 * (temp/228.d0 - 1)**(-1.73d0) * 1.d-12



      x1 = 1.d0 - dexp( -t / tau1 )
      x2 = 1.d0 - dexp( -t / tau2 )
      x3 = 1.d0 - dcos( ws * t ) * dexp( -gs / 2.d0 * t)
      x4 = 1.d0 - dcos( wl * t ) * dexp( -gl / 2.d0 * t)

      s1 = eps1       * x1
      s2 = eps2       * x2
      s3 = as / ws**2 * x3
      s4 = al / wl**2 * x4
      s5 = 2.d0

      er = s1 + s2 + s3 + s4 + s5
c      er = s3 + s4 + s5
c      er = 80.d0

      end


c====================================
      function sample(randam)
c=====================================
      implicit real*8 (a-h,o-z)

       y = 2.d0 - 2.d0 * randam 
       sample = dierfc(y)

      end

c===============================================
! inverse of error function in double precision
!
      function dierfc(y)
c===============================================
      implicit real*8 (a - h, o - z)
      parameter (
     &    qa = 9.16461398268964d-01, 
     &    qb = 2.31729200323405d-01, 
     &    qc = 4.88826640273108d-01, 
     &    qd = 1.24610454613712d-01, 
     &    q0 = 4.99999303439796d-01, 
     &    q1 = 1.16065025341614d-01, 
     &    q2 = 1.50689047360223d-01, 
     &    q3 = 2.69999308670029d-01, 
     &    q4 = -7.28846765585675d-02)
      parameter (
     &    pa = 3.97886080735226000d+00, 
     &    pb = 1.20782237635245222d-01, 
     &    p0 = 2.44044510593190935d-01, 
     &    p1 = 4.34397492331430115d-01, 
     &    p2 = 6.86265948274097816d-01, 
     &    p3 = 9.56464974744799006d-01, 
     &    p4 = 1.16374581931560831d+00, 
     &    p5 = 1.21448730779995237d+00, 
     &    p6 = 1.05375024970847138d+00, 
     &    p7 = 7.13657635868730364d-01, 
     &    p8 = 3.16847638520135944d-01, 
     &    p9 = 1.47297938331485121d-02, 
     &    p10 = -1.05872177941595488d-01, 
     &    p11 = -7.43424357241784861d-02)
      parameter (
     &    p12 = 2.20995927012179067d-03, 
     &    p13 = 3.46494207789099922d-02, 
     &    p14 = 1.42961988697898018d-02, 
     &    p15 = -1.18598117047771104d-02, 
     &    p16 = -1.12749169332504870d-02, 
     &    p17 = 3.39721910367775861d-03, 
     &    p18 = 6.85649426074558612d-03, 
     &    p19 = -7.71708358954120939d-04, 
     &    p20 = -3.51287146129100025d-03, 
     &    p21 = 1.05739299623423047d-04, 
     &    p22 = 1.12648096188977922d-03)
      z = y
      if (y .gt. 1) z = 2 - y
      w = qa - log(z)
      u = sqrt(w)
      s = (qc + log(u)) / w
      t = 1 / (u + qb)
      x = u * (1 - s * (0.5d0 + s * qd)) - 
     &    ((((q4 * t + q3) * t + q2) * t + q1) * t + q0) * t
      t = pa / (pa + x)
      u = t - 0.5d0
      s = (((((((((p22 * u + p21) * u + p20) * u + 
     &    p19) * u + p18) * u + p17) * u + p16) * u + 
     &    p15) * u + p14) * u + p13) * u + p12
      s = ((((((((((((s * u + p11) * u + p10) * u + 
     &    p9) * u + p8) * u + p7) * u + p6) * u + p5) * u + 
     &    p4) * u + p3) * u + p2) * u + p1) * u + p0) * t - 
     &    z * exp(x * x - pb)
      x = x + s * (1 + x * s)
      if (y .gt. 1) x = -x
      dierfc = x
      end



C=================================================================
      FUNCTION ran2(idum)
C=================================================================
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END
