﻿!-----------------------------------------------------------
!    programm kram5sd: Kramers chain with up to 5 beads,
!    with stochastic impulsion. Ryckaert's SHAKE method.
!    F. J. Vesely 2012
!___________________________________________________________
      program kram5sd
      implicit double precision (a-h,o-z)
      character plot,blnk
      common /par/pi,dr,xlq,sig,stol,shist,sish,scsh,nmol,nc,nsh,nhist
      common /vct/x(5),y(5),xn(5),yn(5)
      common /hist/sphi(50),plot(50)
      integer(8) ncm
!
      pi=4.*datan(1.d0);
      blnk=" "
      open(unit=7,file='kram5sd.out',form='formatted'
     1 ,status='unknown')
      open(unit=77,file='kram5sd.dat',form='formatted'
     1 ,status='unknown') ! for gnuplot
      call srand(iran);
      nhist=25;
      do 37 ii=1,nhist
        sphi(ii)=0.;
        plot(ii)=blnk ! ascii plot array
37    continue
!---------------------------------------
! Input parameters:
      nmol=3          ! no of beads (3..5)
      nmol=min(5,max(3,nmol))
      ncm=20*10**6;    ! no of random steps
      nsh=50;         ! max no of chain iterations
      dr=0.25d0        ! step size (max 0.25) 
      dr=min(0.25d0,dr)
      stol=1.d-12;    ! tolerance for sig 
100   format(80a1/3i10)  ! later: read in parameters!
      write(7,776)
776   format(1h1//5x,'Programm kram5sd (stochastic constraint'
     1 ,' dynamics)'/5x,'Input:')
      ncp=ncm/100;  ! print progress to console every ncd steps
      nav=10;       ! histogram is incremented every nav steps
      dphi=pi/float(nhist) ! angle channel width
      xl=1.d0;      ! bond length
      xlq=xl*xl;
!---------------------------------------
!
!  kram5sd simulates the 2-dimensional diffusive motion of a short polymer
!  which is thermally activated by stochastic impulses
!
!----------------------------------------
      write(7,908)
908   format(/4x,'*** Constrained dynamics a la Ryckaert ***')
      write(7,101) 
101   format(/5x,'Diffusing Kramers chain'//)
      do 11 i=1,nmol
        x(i)=i*xl; y(i)=0.;
11    continue
      shist=0.d0     ! total no of histogram increments
      sish=0.d0; scsh=0.d0; ! for loop statistics
!
!  Main loop:
      do 3 nc=1,ncm  
        call step   ! random move for each bead
        call shake  ! repair bond lengths a la Ryckaert
        if((nc/nav)*nav.eq.nc) then
          call incrhist ! increment p(phi) histogram
        endif
3     continue        ! end main loop
!
      call drawhist(dphi)
      write(*,102) sish/scsh
      write(7,102) sish/scsh
102   format(//1x,'Avrg. no of SHAKE iterations: ',f10.3//)
      stop
      end
!
!-----------------------------------------------
!
      subroutine step
! Random move for each bead 
      implicit double precision (a-h,o-z)
      common /par/pi,dr,xlq,sig,stol,shist,sish,scsh,nmol,nc,nsh,nhist
      common /vct/x(5),y(5),xn(5),yn(5)
!
      do 555 i=1,nmol  ! loop over beads: stochastic step
        xgr=rand()-0.5d0; xn(i)=x(i)+dr*xgr;
        xgr=rand()-0.5d0; yn(i)=y(i)+dr*xgr;
555   continue
      return
      end
!
!------------------------------------
!
      subroutine shake
! Repair bond lengths a la Ryckaert
      implicit double precision (a-h,o-z)
      common /par/pi,dr,xlq,sig,stol,shist,sish,scsh,nmol,nc,nsh,nhist
      common /vct/x(5),y(5),xn(5),yn(5)
      logical lbreak
!
      do 1 ish=1,nsh !   // Loop over SHAKE iterations
        lbreak=.true.;
        do 2 i=2,nmol-1
          do 3 k=i-1,i+1,2
            dxn=xn(k)-xn(i); dyn=yn(k)-yn(i);
            sig=xlq-dxn*dxn-dyn*dyn;
 !        print *,'sh: nc,ish,dxn,dyn=',nc,ish,dxn,dyn
            if(dabs(sig).gt.stol) then
              lbreak=.false.
              dx=x(k)-x(i);    dy=y(k)-y(i);
              rn=dxn*dx+dyn*dy;
              cf=-sig/4./rn; ! forward correction
              cx=cf*dx; cy=cf*dy;
              xn(i)=xn(i)+cx; yn(i)=yn(i)+cy;
              xn(k)=xn(k)-cx; yn(k)=yn(k)-cy;
            endif
3         continue
2       continue
        if(lbreak) exit
1     continue
      sish=sish+ish; scsh=scsh+1.d0;
!
      do 4 i=1,nmol
        x(i)=xn(i); y(i)=yn(i);
4     continue
!
      return
      end
!
!------------------------------------
!
      subroutine incrhist
! Increment histogram of prob. dens. p(phi)
      implicit double precision (a-h,o-z)
      character plot
      common /par/pi,dr,xlq,sig,stol,shist,sish,scsh,nmol,nc,nsh,nhist
      common /vct/x(5),y(5),xn(5),yn(5)
      common /hist/sphi(50),plot(50)
!   increment histogram:
      iv=(nmol+1)/2  ! somewhere in mid-chain
      dxb=xn(iv)-xn(iv-1); dyb=yn(iv)-yn(iv-1); 
      dnb=dsqrt(dxb*dxb+dyb*dyb);
      dxb=dxb/dnb; dyb=dyb/dnb; ! unit vector along (1-2) 
      dxf=xn(iv+1)-xn(iv); dyf=yn(iv+1)-yn(iv); 
      dnf=dsqrt(dxf*dxf+dyf*dyf);
      dxf=dxf/dnf; dyf=dyf/dnf; ! unit vector along (2-3) 
!
      arg=dxb*dxf+dyb*dyf
      phi=dacos(arg); ! bond angle 1-2-3 or similar

      ic=int(nhist*phi/pi)+1
      sphi(ic)=sphi(ic)+1.d0
      shist=shist+1.d0
!
      return
      end

!------------------------------------
!
      subroutine drawhist(dphi)
! Sketch histogram of prob. dens. p(phi)
      implicit double precision (a-h,o-z)
      character plot
      common /par/pi,dr,xlq,sig,stol,shist,sish,scsh,nmol,nc,nsh,nhist
      common /vct/x(5),y(5),xn(5),yn(5)
      common /hist/sphi(50),plot(50)
!  Draw distribution of bond angles:
      write(7,1666)
1666  format(1h1//5x,'Distribution density of angle phi:'//)
      pmax=0.d0
      do 12 i=1,nhist
        pmax=max(pmax,sphi(i))
12    continue
      phi=-dphi/2.
      do 499 i=1,nhist
        phi=phi+dphi
        cxp=sphi(i)/pmax
        cth=1.-0.25*dcos(phi)**2
        cth=dsqrt(cth)
        ipe=int(nhist*cxp)+1
        if(ipe.gt.nhist) ipe=nhist
        ipt=int(nhist*cth)+1
        if(ipt.gt.nhist) ipt=nhist
        plot(1)='|'; plot(ipe)='1'; plot(ipt)='0'; plot(nhist)='|';
        write(*,5000) i,phi,cxp,cth,cxp-cth,(plot(ii),ii=1,nhist)
        write(7,5000) i,phi,cxp,cth,cxp-cth,(plot(ii),ii=1,nhist)
        write(77,5000) i,phi,cxp,cth,cxp-cth
5000   format(1x,i5,4f8.4,50a1)
        plot(ipe)=' ';
        plot(ipt)=' ';
        plot(nhist)=' ';
499   continue
      return
      end
!------------------------------------------
