subroutine friction_force
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.MD'
#ifndef LANG0
      include 'COMMON.LANGEVIN'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
      include 'COMMON.IOUNITS'
      double precision gamvec(MAXRES6)
      common /syfek/ gamvec
      double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
     & ginvfric(maxres2,maxres2)
      common /przechowalnia/ ginvfric
      
      logical lprn /.false./, checkmode /.false./

      do i=0,MAXRES2
        do j=1,3
          friction(j,i)=0.0d0
        enddo
      enddo
  
      do j=1,3
        d_t_work(j)=d_t(j,0)
      enddo
      ind=3
      do i=nnt,nct-1
        do j=1,3
          d_t_work(ind+j)=d_t(j,i)
        enddo
        ind=ind+3
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          do j=1,3
            d_t_work(ind+j)=d_t(j,i+nres)
          enddo
          ind=ind+3
        endif
      enddo

      call fricmat_mult(d_t_work,fric_work)
      
      if (.not.checkmode) return

      if (lprn) then
        write (iout,*) "d_t_work and fric_work"
        do i=1,3*dimen
          write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
        enddo
      endif
      do j=1,3
        friction(j,0)=fric_work(j)
      enddo
      ind=3
      do i=nnt,nct-1
        do j=1,3
          friction(j,i)=fric_work(ind+j)
        enddo
        ind=ind+3
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          do j=1,3
            friction(j,i+nres)=fric_work(ind+j)
          enddo
          ind=ind+3
        endif
      enddo
      if (lprn) then
        write(iout,*) "Friction backbone"
        do i=0,nct-1
          write(iout,'(i5,3e15.5,5x,3e15.5)') 
     &     i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
        enddo
        write(iout,*) "Friction side chain"
        do i=nnt,nct
          write(iout,'(i5,3e15.5,5x,3e15.5)') 
     &     i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
        enddo   
      endif
      if (lprn) then
        do j=1,3
          vv(j)=d_t(j,0)
        enddo
        do i=nnt,nct
          do j=1,3
            vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
            vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
            vv(j)=vv(j)+d_t(j,i)
          enddo
        enddo
        write (iout,*) "vvtot backbone and sidechain"
        do i=nnt,nct
          write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
     &     (vvtot(j,i+nres),j=1,3)
        enddo
        ind=0
        do i=nnt,nct-1
          do j=1,3
            v_work(ind+j)=vvtot(j,i)
          enddo
          ind=ind+3
        enddo
        do i=nnt,nct
          do j=1,3
            v_work(ind+j)=vvtot(j,i+nres)
          enddo
          ind=ind+3
        enddo
        write (iout,*) "v_work gamvec and site-based friction forces"
        do i=1,dimen1
          write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
     &      gamvec(i)*v_work(i) 
        enddo
c        do i=1,dimen
c          fric_work1(i)=0.0d0
c          do j=1,dimen1
c            fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
c          enddo
c        enddo  
c        write (iout,*) "fric_work and fric_work1"
c        do i=1,dimen
c          write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
c        enddo 
        do i=1,dimen
          do j=1,dimen
            ginvfric(i,j)=0.0d0
            do k=1,dimen
              ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
            enddo
          enddo
        enddo
        write (iout,*) "ginvfric"
        do i=1,dimen
          write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
        enddo
        write (iout,*) "symmetry check"
        do i=1,dimen
          do j=1,i-1
            write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
          enddo   
        enddo
      endif 
      return
      end
c-----------------------------------------------------
      subroutine stochastic_force(stochforcvec)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.MD'
      include 'COMMON.TIME1'
#ifndef LANG0
      include 'COMMON.LANGEVIN'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
      include 'COMMON.IOUNITS'
      
      double precision x,sig,lowb,highb,
     & ff(3),force(3,0:MAXRES2),zeta2,lowb2,
     & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6)
      logical lprn /.false./
      do i=0,MAXRES2
        do j=1,3
          stochforc(j,i)=0.0d0
        enddo
      enddo
      x=0.0d0 	

#ifdef MPI
      time00=MPI_Wtime()
#else
      time00=tcpu()
#endif
c Compute the stochastic forces acting on bodies. Store in force.
      do i=nnt,nct-1
        sig=stdforcp(i)
        lowb=-5*sig
        highb=5*sig
        do j=1,3
          force(j,i)=anorm_distr(x,sig,lowb,highb)
        enddo
      enddo
      do i=nnt,nct
        sig2=stdforcsc(i)
        lowb2=-5*sig2
        highb2=5*sig2
        do j=1,3
          force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
        enddo
      enddo
#ifdef MPI
      time_fsample=time_fsample+MPI_Wtime()-time00
#else
      time_fsample=time_fsample+tcpu()-time00
#endif
c Compute the stochastic forces acting on virtual-bond vectors.
      do j=1,3
        ff(j)=0.0d0
      enddo
      do i=nct-1,nnt,-1
        do j=1,3
          stochforc(j,i)=ff(j)+0.5d0*force(j,i)
        enddo
        do j=1,3
          ff(j)=ff(j)+force(j,i)
        enddo
        if (itype(i+1).ne.21) then
          do j=1,3
            stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
            ff(j)=ff(j)+force(j,i+nres+1)
          enddo
        endif
      enddo 
      do j=1,3
        stochforc(j,0)=ff(j)+force(j,nnt+nres)
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          do j=1,3
            stochforc(j,i+nres)=force(j,i+nres)
          enddo
        endif
      enddo 

      do j=1,3
        stochforcvec(j)=stochforc(j,0)
      enddo
      ind=3
      do i=nnt,nct-1
        do j=1,3 
          stochforcvec(ind+j)=stochforc(j,i)
        enddo
        ind=ind+3
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          do j=1,3
            stochforcvec(ind+j)=stochforc(j,i+nres)
          enddo
          ind=ind+3
        endif
      enddo
      if (lprn) then
        write (iout,*) "stochforcvec"
        do i=1,3*dimen
          write(iout,'(i5,e15.5)') i,stochforcvec(i)
        enddo
        write(iout,*) "Stochastic forces backbone"
        do i=0,nct-1
          write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
        enddo
        write(iout,*) "Stochastic forces side chain"
        do i=nnt,nct
          write(iout,'(i5,3e15.5)') 
     &      i,(stochforc(j,i+nres),j=1,3)
        enddo   
      endif

      if (lprn) then

      ind=0
      do i=nnt,nct-1
        write (iout,*) i,ind
        do j=1,3
          forcvec(ind+j)=force(j,i)
        enddo
        ind=ind+3
      enddo
      do i=nnt,nct
        write (iout,*) i,ind
        do j=1,3
          forcvec(j+ind)=force(j,i+nres)
        enddo
        ind=ind+3
      enddo 

      write (iout,*) "forcvec"
      ind=0
      do i=nnt,nct-1
        do j=1,3
          write (iout,'(2i3,2f10.5)') i,j,force(j,i),
     &      forcvec(ind+j)
        enddo
        ind=ind+3
      enddo
      do i=nnt,nct
        do j=1,3
          write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
     &     forcvec(ind+j)
        enddo
        ind=ind+3
      enddo

      endif

      return
      end
c------------------------------------------------------------------
      subroutine setup_fricmat
      implicit real*8 (a-h,o-z)
#ifdef MPI
      include 'mpif.h'
#endif
      include 'DIMENSIONS'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.MD'
      include 'COMMON.SETUP'
      include 'COMMON.TIME1'
c      integer licznik /0/
c      save licznik
#ifndef LANG0
      include 'COMMON.LANGEVIN'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
      include 'COMMON.IOUNITS'
      integer IERROR
      integer i,j,ind,ind1,m
      logical lprn /.false./
      double precision dtdi,gamvec(MAXRES2),
     &  ginvfric(maxres2,maxres2),Ghalf(mmaxres2),fcopy(maxres2,maxres2)
      common /syfek/ gamvec
      double precision work(8*maxres2)
      integer iwork(maxres2)
      common /przechowalnia/ ginvfric,Ghalf,fcopy
#ifdef MPI
      if (fg_rank.ne.king) goto 10
#endif
c  Zeroing out fricmat
      do i=1,dimen
        do j=1,dimen
          fricmat(i,j)=0.0d0
        enddo	 
      enddo
c  Load the friction coefficients corresponding to peptide groups
      ind1=0
      do i=nnt,nct-1
        ind1=ind1+1
        gamvec(ind1)=gamp
      enddo
c  Load the friction coefficients corresponding to side chains
      m=nct-nnt
      ind=0
      do i=nnt,nct
        ind=ind+1
        ii = ind+m
        iti=itype(i)
        gamvec(ii)=gamsc(iti)
      enddo
      if (surfarea) call sdarea(gamvec)
c      if (lprn) then
c        write (iout,*) "Matrix A and vector gamma"
c        do i=1,dimen1
c          write (iout,'(i2,$)') i
c          do j=1,dimen
c            write (iout,'(f4.1,$)') A(i,j)
c          enddo
c          write (iout,'(f8.3)') gamvec(i)
c        enddo
c      endif
      if (lprn) then
        write (iout,*) "Vector gamvec"
        do i=1,dimen1
          write (iout,'(i5,f10.5)') i, gamvec(i)
        enddo
      endif
        
c The friction matrix       
      do k=1,dimen
       do i=1,dimen
         dtdi=0.0d0
         do j=1,dimen1
           dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
         enddo
         fricmat(k,i)=dtdi
       enddo
      enddo 

      if (lprn) then
        write (iout,'(//a)') "Matrix fricmat"
        call matout2(dimen,dimen,maxres2,maxres2,fricmat)
      endif 
      if (lang.eq.2 .or. lang.eq.3) then
c Mass-scale the friction matrix if non-direct integration will be performed
      do i=1,dimen
        do j=1,dimen
          Ginvfric(i,j)=0.0d0
          do k=1,dimen
            do l=1,dimen
              Ginvfric(i,j)=Ginvfric(i,j)+
     &          Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
            enddo
          enddo
        enddo
      enddo
c Diagonalize the friction matrix
      ind=0
      do i=1,dimen
        do j=1,i
          ind=ind+1
          Ghalf(ind)=Ginvfric(i,j)
        enddo
      enddo
      call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
     &  ierr,iwork)
      if (lprn) then
        write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
     &    " mass-scaled friction matrix"
        call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
      endif
c Precompute matrices for tinker stochastic integrator
#ifndef LANG0
      do i=1,dimen
        do j=1,dimen
          mt1(i,j)=0.0d0
          mt2(i,j)=0.0d0
          do k=1,dimen
            mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
            mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)		   
          enddo
          mt3(j,i)=mt1(i,j)
        enddo
      enddo
#endif
      else if (lang.eq.4) then
c Diagonalize the friction matrix
      ind=0
      do i=1,dimen
        do j=1,i
          ind=ind+1
          Ghalf(ind)=fricmat(i,j)
        enddo
      enddo
      call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
     &  ierr,iwork)
      if (lprn) then
        write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
     &    " friction matrix"
        call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
      endif
c Determine the number of zero eigenvalues of the friction matrix
      nzero=max0(dimen-dimen1,0)
c      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
c        nzero=nzero+1
c      enddo
      write (iout,*) "Number of zero eigenvalues:",nzero
      do i=1,dimen
        do j=1,dimen
          fricmat(i,j)=0.0d0
          do k=nzero+1,dimen
            fricmat(i,j)=fricmat(i,j)
     &        +fricvec(i,k)*fricvec(j,k)/fricgam(k) 
          enddo
        enddo
      enddo
      if (lprn) then
        write (iout,'(//a)') "Generalized inverse of fricmat"
        call matout(dimen,dimen,maxres6,maxres6,fricmat)
      endif 
      endif
#ifdef MPI
  10  continue
      if (nfgtasks.gt.1) then
        if (fg_rank.eq.0) then
c The matching BROADCAST for fg processors is called in ERGASTULUM
#ifdef MPI
          time00=MPI_Wtime()
#else
          time00=tcpu()
#endif
          call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
#ifdef MPI
          time_Bcast=time_Bcast+MPI_Wtime()-time00
#else
          time_Bcast=time_Bcast+tcpu()-time00
#endif
c          print *,"Processor",myrank,
c     &       " BROADCAST iorder in SETUP_FRICMAT"
        endif
c	licznik=licznik+1
c        write (iout,*) "setup_fricmat licznik",licznik
#ifdef MPI
        time00=MPI_Wtime()
#else
        time00=tcpu()
#endif
c Scatter the friction matrix
        call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
     &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
#ifdef TIMING
#ifdef MPI
        time_scatter=time_scatter+MPI_Wtime()-time00
        time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
#else
        time_scatter=time_scatter+tcpu()-time00
        time_scatter_fmat=time_scatter_fmat+tcpu()-time00
#endif
#endif
        do i=1,dimen
          do j=1,2*my_ng_count
            fricmat(j,i)=fcopy(i,j)
          enddo
        enddo
c        write (iout,*) "My chunk of fricmat"
c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
      endif
#endif
      return
      end
c-------------------------------------------------------------------------------
      subroutine sdarea(gamvec)
c
c Scale the friction coefficients according to solvent accessible surface areas
c Code adapted from TINKER
c AL 9/3/04
c
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CONTROL'
      include 'COMMON.VAR'
      include 'COMMON.MD'
#ifndef LANG0
      include 'COMMON.LANGEVIN'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.IOUNITS'
      include 'COMMON.NAMES'
      double precision radius(maxres2),gamvec(maxres6)
      parameter (twosix=1.122462048309372981d0)
      logical lprn /.false./
c
c     determine new friction coefficients every few SD steps
c
c     set the atomic radii to estimates of sigma values
c
c      print *,"Entered sdarea"
      probe = 0.0d0
      
      do i=1,2*nres
        radius(i)=0.0d0
      enddo
c  Load peptide group radii
      do i=nnt,nct-1
        radius(i)=pstok
      enddo
c  Load side chain radii
      do i=nnt,nct
        iti=itype(i)
        radius(i+nres)=restok(iti)
      enddo
c      do i=1,2*nres
c        write (iout,*) "i",i," radius",radius(i) 
c      enddo
      do i = 1, 2*nres
         radius(i) = radius(i) / twosix
         if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
      end do
c
c     scale atomic friction coefficients by accessible area
c
      if (lprn) write (iout,*) 
     &  "Original gammas, surface areas, scaling factors, new gammas, ",
     &  "std's of stochastic forces"
      ind=0
      do i=nnt,nct-1
        if (radius(i).gt.0.0d0) then
          call surfatom (i,area,radius)
          ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
          if (lprn) write (iout,'(i5,3f10.5,$)') 
     &      i,gamvec(ind+1),area,ratio
          do j=1,3
            ind=ind+1
            gamvec(ind) = ratio * gamvec(ind)
          enddo
          stdforcp(i)=stdfp*dsqrt(gamvec(ind))
          if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
        endif
      enddo
      do i=nnt,nct
        if (radius(i+nres).gt.0.0d0) then
          call surfatom (i+nres,area,radius)
          ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
          if (lprn) write (iout,'(i5,3f10.5,$)') 
     &      i,gamvec(ind+1),area,ratio
          do j=1,3
            ind=ind+1 
            gamvec(ind) = ratio * gamvec(ind)
          enddo
          stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
          if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
        endif
      enddo

      return
      end