!RADIAL DISTRIBUTION FUNCTION CODE (WORK FOR BOTH NVT & NPT)
       program gofr
!===============================================================================
! Gromacs trajectroy reading utility available on github
!===============================================================================
       use gmxfort_trajectory
       use gmxfort_utils        
       
       implicit none
       
       type(Trajectory) :: trj
!=================================================================================
! This section is for the trajectory file input variables
!=================================================================================
       integer :: ninput,ntot,inp,nsites,nframes ! input for gofr
       integer,parameter :: maxinp=150      ! maximum number of trajectory files as an input
       character(len=100) :: infile(maxinp)
       integer :: ibegin(maxinp),iend(maxinp),nskip,nread,nused,intv

       integer :: i1,j1,u1,n1
       integer :: i,j,k,ia,im,ja,jm
       integer :: total_atom,natms,nstep
       integer :: nmol_cta,nmol_br,nmol_water,natom_cta,natom_br,natom_water
       integer :: atom_cbw,ifile
       integer :: total_frame
       integer :: iatom_cta,iatom_br,iatom_water,imat,jmat
       integer :: imol_cta,imol_br,imol_water
       integer :: no_atname1
       integer :: ibin,nprtot,npr(1000),nbins
       integer :: iflg,keytrj,imcon,nmol_diff,new_file
       
       real :: myatom(3)
       real(kind=8),allocatable :: chge(:),weight(:)!,yyy(:),zzz(:)
       real(kind=8),allocatable :: vxx(:),vyy(:),vzz(:),fxx(:),fyy(:),fzz(:)
       real(kind=8):: xx1(1),yy1(1),zz1(1),xxx(1),yyy(1),zzz(1)
       real(kind=8),allocatable :: rw(:,:,:),rcta(:,:,:),rbr(:,:,:)
       real(kind=8) :: boxl(3,3),dr(3),rmsx,rmsy,rmsz
       real(kind=8) :: timstp
       real(kind=8),parameter :: au=0.5291770
       !================================================================
       real(kind=8),parameter :: drg=0.05            ! spacing in g(r)
       !================================================================
       real(kind=8) :: cell(3,3),rcell(9),det,dt0
       !================================================================
       real(kind=8) :: rmax,rdr,rsq,r,rr,rrs,cn,g,gavg,fact,cutoff_vol
       real(kind=8) :: pi
       real(kind=8),allocatable :: rs(:),gs(:)
       real :: start,finish
       
       character(len=10) :: nchar
       character(len=100) :: outfile1,outfile2,a
       character(len=8),allocatable :: atnm_rw(:,:),atnm_rcta(:,:),atnm_rbr(:,:),x(:,:)      
       character(len=8) :: atname1,atname2
       character(len=8),allocatable :: atname(:)     
       character(len=80) :: title
      
      ! write (*,*) "type the trjectory file name"
      call cpu_time(start)
!=======================================================================
      open(unit=2,file='br_br_gofr.inp',action='read')
       read(2,*)outfile1
       read(2,*)outfile2
       open(unit=3,form='formatted',file=outfile1,&
       action='write')
       open(unit=4,form='formatted',file=outfile2,&
       action='write')
!=======================================================================
!=======================================================================
! Begining of the code for gofr calculation
!=======================================================================

       write(*,*)'Enter the no of trajectory file you want to analyze'
       read(2,*)ninput    												! (maxinp=150)
       write(*,*)'Along with this also enter the number of frames you',&
               ' want to work with.'
       read(2,*)nframes    												! number of time-steps used for binning of atoms
       write(*,*)'Enter the number of different molecules types'
       read(2,*)nmol_diff ! (1)
       write(*,*)'Enter the no of molecules of all types'
       read(2,*)nmol_water ! (884) 
       write(*,*)'Enter the no of atom in each of the diff molecules'
       read(2,*) natom_water ! (3)
               if(ninput>maxinp)then
               write(*,*)'Too many inputs'
       stop
       endif
               write(*,*)'Pair correlation atom pair X-Y: X'
               read(2,*)atname1
               write(*,*)'Pair correlation atom pair X-Y: Y'
               read(2,*)atname2
               write(*,*)'Enter the number of X atom'
               read(2,*)no_atname1
               
!=======================================================================
            atom_cbw=0
            pi=datan(1.0d0)*4
            intv=10
            write(*,*)'pi====>',pi
            total_atom=nmol_water*natom_water
            natms=total_atom
!=======================================================================
            ntot=0                                                      ! total number of inputs in trajectory             
!=======================================================================
       write(*,*)'enter the skip value'
       read(2,*)nskip
       read(2,*)rmax                                                    ! half the box length
       !write(4,*)'nskip,rmax=======>',nskip,rmax
       if(nskip/=0)nskip=1
       !ntot=ntot/nskip
       !total_frame=ntot
       write(*,*)'Finished reading inputs'
!=======================================================================
       !GOFR CALCULATION
!=======================================================================
       rdr=1.0d0/drg
       nbins=nint(rmax*rdr)+1                                           !number of bins for calculation has been evaluated
       if(real(nbins)*drg>rmax)nbins=nbins-1
       do j=1,nbins
       npr(j)=0
       enddo
       nprtot=0
       !write(4,*)'rmax =',rmax,'nbins =',nbins
!=======================================================================
! from this line of code below trajectory file name was taken directly 
! from the file
!=======================================================================
fread:       do inp=1,ninput,1
             write(*,*) 'Enter the name of trajectory file',inp
             read (2,*) infile(inp)
             read (2,*) ibegin(inp),iend(inp)
!=======================================================================
!Number of atoms was going to be read and passed to a variable name nsites
!=======================================================================
        call trj%open(infile(inp))
        nsites = trj%natoms()
!-----------------------------------------------------------------------
			allocate(rw(nmol_water,natom_water,3))!,rcta(nmol_cta,natom_cta,3),rbr(nmol_br,natom_br,3))
            allocate(atnm_rw(nmol_water,natom_water))!,atnm_rcta(nmol_cta,natom_cta),atnm_rbr(nmol_br,natom_br))
           ! allocate(atname(nsites))
            !allocate(xx1(nsites),yy1(nsites),zz1(nsites))
            !allocate(xxx(nsites),yyy(nsites),zzz(nsites))
!-----------------------------------------------------------------------            
        n1 = trj%read_next(nframes)                                     !trj%read_next lib allow us to read frame wise data
!=======================================================================
               !write(4,*)'number of sites = ',nsites
               !write(4,*)'Total atom=',total_atom
             !  allocate(xxx(nsites),yyy(nsites),zzz(nsites))
              ! allocate(xx1(nsites),yy1(nsites),zz1(nsites))
!=======================================================================
!here the code will loop over the provided number of frames you want to 
!work with as an input trajectory
!=======================================================================
              
nframe:       do i1=1, n1                                               ! for nframes frames as mentioned
            !  write (3,*)"frame=========",i1
             ! write(*,*) i1,"frame============"
              cell= trj%box(i1)
              
              call invert (cell,rcell,det)
              
site:         do j1=1,nsites
              myatom = trj%x(i1,j1)
              xxx = trj%frameArray(i1)%xyz(1,j1)
              yyy = trj%frameArray(i1)%xyz(2,j1)
              zzz = trj%frameArray(i1)%xyz(3,j1)
                    write(3,*) j1,"atom-----------------------------"
                     write(3,*) xxx , yyy, zzz 
              !write(*,*)'h1',j1,xxx(j1)
                     xx1=xxx(j1)*rcell(1)+yyy(j1)*rcell(4)+zzz(j1)*rcell(7)      !
                     yy1=xxx(j1)*rcell(2)+yyy(j1)*rcell(5)+zzz(j1)*rcell(8)      !
                     zz1=xxx(j1)*rcell(3)+yyy(j1)*rcell(6)+zzz(j1)*rcell(9)
                     enddo site 
                     !write(3,*) j1,a
                    ! write(4,*) j1,"atom  ","  frame",i1
              !
                     !write(4,*) j1," atom ",xx1,yy1,zz1
!100  format (i5,2x,)                        !
               
wt: 		   do imol_water=1,nmol_water
				do iatom_water=1,natom_water
					atom_cbw=atom_cbw+1
					rw(imol_water,iatom_water,1)=xx1(atom_cbw)
					rw(imol_water,iatom_water,2)=yy1(atom_cbw)
					rw(imol_water,iatom_water,3)=zz1(atom_cbw)
					atnm_rw(imol_water,iatom_water)=atname(atom_cbw)
					!write(*,*) rw
					enddo
					enddo wt
					
					do imat=1,3
					do jmat=1,3
					boxl(imat,jmat)=cell(imat,jmat)
					enddo
					enddo
					
br1:   do im=1,nmol_water-1
       do ia=1,natom_water
       if(atnm_rw(im,ia)==atname1)then                                  !pair function atm1 (cta)
br2:   do jm=(im+1),nmol_water
       do ja=1,natom_water
       if(atnm_rw(jm,ja)==atname2)then                                  !pair function atm2 (cta)
       
       dr(:)=rw(im,ia,:)-rw(jm,ja,:)
       dr(:)=dr(:)-anint(dr(:))                                         !Box length L=1
       
       !CALCULATE CORRECTED REAL SPACE COORDINATES
       rmsx=dr(1)*cell(1,1)+dr(2)*cell(2,1)+dr(3)*cell(3,1)
       rmsy=dr(1)*cell(1,2)+dr(2)*cell(2,2)+dr(3)*cell(3,2)
       rmsz=dr(1)*cell(1,3)+dr(2)*cell(2,3)+dr(3)*cell(3,3)
       
       dr(1)=rmsx;dr(2)=rmsy;dr(3)=rmsz
       rsq=dot_product(dr,dr)
       r=(sqrt(rsq))
       ibin=int(r*rdr)+1
       
       if(ibin<=nbins)then
       npr(ibin)=npr(ibin)+1
       nprtot=nprtot+1
       endif
       
       endif
       enddo
       enddo br2
       endif
       enddo
       enddo br1
               
               enddo nframe
               !
               call trj%close()
               enddo fread
100 		   continue
!-----------------------------------------------------------------------
				cutoff_vol=((4.0d0/3.0d0)*pi*(rmax**3))
!-----------------------------------------------------------------------
			cn=0.0
			fact=(cutoff_vol)/(4.0d0*pi*(drg**3)*real(nprtot))
			allocate(rs(nbins),gs(nbins))
			do i=1,nbins
			r=(real(i)-0.5d0)*drg
			rr=((real(i)+0.5d0)**2+(1.d0/12.d0))
			g=fact*(real(npr(i))/(rr))
!--------------------------------------------------------
       cn=cn+real(npr(i))/real(nused*(no_atname1-1))
!--------------------------------------------------------
       rs(i)=(real(i)-0.5d0)*drg
       rrs=((real(i)+0.5d0)**2+(1.d0/12.d0))
       gs(i)=fact*(real(npr(i))/(rrs))
       write(3,2)r,g,cn
       enddo
!=========================================================
2      format(2x,f10.5,f10.5,f10.5)
3      format(2x,f10.5,f10.5)
!-----------------------------------------
       if(a=='smooth')then
       write(4,3)rs(1),gs(1)
       write(4,3)rs(2),gs(2)
       do i=3,nbins-2
       gavg=0.0d0
       do j=i-2,i+2
       gavg=gavg+gs(j)
       enddo
       gavg=(gavg/5.0)
       write(4,3)rs(i),gavg
       enddo
       write(4,3)rs(nbins-1),gs(nbins-1)
       write(4,3)rs(nbins),gs(nbins)
       endif        !For 5 point smoothing
!-----------------------------------------------------------------------

              ! deallocate(xxx,yyy,zzz,xx1,yy1,zz1)
               deallocate(rw,atnm_rw,atname)
               call cpu_time(finish)
               print'("Time= ",f6.3,"seconds")',finish-start,finish,start
               end program gofr
!---------------------------------------------------------------
      subroutine invert(box,b,d) !subroutine sname (dummy arugments)
!---------------------------------------------------------------
!     dl_poly subroutine to invert a 3*3 matrix using cofactors
!     author - w. smith   april 1992
!---------------------------------------------------------------
      real(kind=8):: box(3,3) ! data passed to subroutine that cann't be changed 
      real(kind=8)::b(9),d
      real(kind=8):: r ! r is a temporary variable, not accessible to main program
      write(4,*) "in subroutine"
       do i=1,3
              write(4,*) (box(i,j),j=1,3)
              enddo
!     calculate adjoint matrix cofactors with simultaneous doing transpose of mtx
      b(1)=box(2,2)*box(3,3)-(box(2,3)*box(3,2))
      b(2)=box(3,2)*box(1,3)-(box(1,2)*box(3,3))
      b(3)=box(1,2)*box(2,3)-(box(2,2)*box(1,3))
      b(4)=box(3,1)*box(2,3)-(box(2,1)*box(3,3))
      b(5)=box(1,1)*box(3,3)-(box(3,1)*box(1,3))
      b(6)=box(2,1)*box(1,3)-(box(1,1)*box(2,3))
      b(7)=box(2,1)*box(3,2)-(box(2,2)*box(1,3))
      b(8)=box(1,2)*box(3,1)-(box(1,1)*box(3,2))
      b(9)=box(1,1)*box(2,2)-(box(2,1)*box(1,2))
!     calculate determinant
      d=box(1,1)*b(1)+box(1,2)*b(4)+box(1,3)*b(7)
      r=0.d0
      if(abs(d)>0.d0)r=1.d0/d
!     complete inverse matrix with all elements of matrix calculation
      b(1)=r*b(1)
      b(2)=r*b(2)
      b(3)=r*b(3)
      b(4)=r*b(4)
      b(5)=r*b(5)
      b(6)=r*b(6)
      b(7)=r*b(7)
      b(8)=r*b(8)
      b(9)=r*b(9)
      do i=1,9
              write(4,*)i," element of inverse matrix ",b(i)
              enddo
              write(4,*)"determinant",d
      return     ! to return to the main program
      end subroutine invert ! subroutine program was closed
      !-g -fcheck=all -Wall -fbacktrace command to check error
!----------------------------------------------------------------------