| 
					
				 | 
			
			
				@@ -0,0 +1,166 @@ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+!///////////////////////////////////////////////////////////////////////////////////////// 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+!// Author: Manish Agarwal and Gourav Shrivastava  , IIT Delhi 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+!///////////////////////////////////////////////////////////////////////////////////////// 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+! Copyright (c) 2021 NVIDIA Corporation.  All rights reserved. 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+module readdata 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      contains 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      subroutine readdcd(maxframes,maxatoms,x,y,z,xbox,ybox,zbox,natoms,nframes) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      integer i,j 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      integer maxframes,maxatoms 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      double precision d(6),xbox,ybox,zbox 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      real*4, allocatable   :: x(:,:) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      real*4, allocatable   :: y(:,:) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      real*4, allocatable   :: z(:,:) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      real*4 dummyr 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      integer*4 nset, natoms, dummyi,nframes,tframes 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      character*4 dummyc 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+       
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      open(10,file='../input/alk.traj.dcd',status='old',form='unformatted') 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      read(10) dummyc, tframes,(dummyi,i=1,8),dummyr, (dummyi,i=1,9) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      read(10) dummyi, dummyr,dummyr 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      read(10) natoms 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      print*,"Total number of frames and atoms are",tframes,natoms 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      allocate ( x(maxframes,natoms) ) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      allocate ( y(maxframes,natoms) ) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      allocate ( z(maxframes,natoms) ) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      do i = 1,nframes 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+           read(10) (d(j),j=1, 6) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+           read(10) (x(i,j),j=1,natoms) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+           read(10) (y(i,j),j=1,natoms) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+           read(10) (z(i,j),j=1,natoms) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      end do 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+       
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      xbox=d(1) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      ybox=d(3) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      zbox=d(6) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+       
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      print*,"File reading is done: xbox,ybox,zbox",xbox,ybox,zbox 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      return 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      end subroutine readdcd 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ end module readdata 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+program rdf 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      use readdata 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      use nvtx 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      implicit none 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      integer n,i,j,iconf,ind 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      integer natoms,nframes,nbin 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      integer maxframes,maxatoms 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      parameter (maxframes=10,maxatoms=60000,nbin=2000) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      real*4, allocatable   :: x(:,:) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      real*4, allocatable   :: y(:,:) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      real*4, allocatable   :: z(:,:) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      double precision dx,dy,dz 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      double precision xbox,ybox,zbox,cut 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      double precision vol,r,del,s2,s2bond 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      double precision, allocatable   ::  g(:) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      double precision rho,gr,lngr,lngrbond,pi,const,nideal,rf 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      double precision rlower,rupper 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      character  atmnm*4 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      real*4 start,finish 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+         
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      open(23,file='RDF.dat',status='unknown') 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      open(24,file='Pair_entropy.dat',status='unknown') 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      nframes=10 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call cpu_time(start) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      print*,"Going to read coordinates" 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call nvtxStartRange("Read File") 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call readdcd(maxframes,maxatoms,x,y,z,xbox,ybox,zbox,natoms,nframes) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call nvtxEndRange 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      allocate ( g(nbin) ) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      g = 0.0d0 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+  
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      pi=dacos(-1.0d0) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      vol=xbox*ybox*zbox 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      rho=dble(natoms)/vol 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      del=xbox/dble(2.0*nbin) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      write(*,*) "bin width is : ",del 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      cut = dble(xbox * 0.5); 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      !pair calculation 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      !$omp target data map(x(:,:), y (:,:), z (:,:), g (:)) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call nvtxStartRange("Pair Calculation") 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      do iconf=1,nframes 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+         if (mod(iconf,1).eq.0) print*,iconf 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+         !$omp target teams distribute num_teams(65535) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+         do i=1,natoms 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            !$omp parallel do private(dx,dy,dz,r,ind) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            do j=1,natoms 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               dx=x(iconf,i)-x(iconf,j) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               dy=y(iconf,i)-y(iconf,j) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               dz=z(iconf,i)-z(iconf,j) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               dx=dx-nint(dx/xbox)*xbox 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               dy=dy-nint(dy/ybox)*ybox 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               dz=dz-nint(dz/zbox)*zbox 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+    
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               r=dsqrt(dx**2+dy**2+dz**2) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               ind=int(r/del)+1 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               !if (ind.le.nbin) then 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               if(r<cut)then 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+                  !$omp atomic 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+                  g(ind)=g(ind)+1.0d0 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+               endif 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            enddo 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+         enddo 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      enddo 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call nvtxEndRange 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      !$omp end target data 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      !entropy calculation 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      s2=0.01d0 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      s2bond=0.01d0  
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      const=(4.0d0/3.0d0)*pi*rho 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call nvtxStartRange("Entropy Calculation") 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      do i=1,nbin 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          rlower=dble((i-1)*del) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          rupper=rlower+del 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          nideal=const*(rupper**3-rlower**3) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          g(i)=g(i)/(dble(nframes)*dble(natoms)*nideal) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          r=dble(i)*del 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          if (r.lt.2.0) then 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            gr=0.0 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          else 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            gr=g(i) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          endif 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          if (gr.lt.1e-5) then 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            lngr=0.0 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          else 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            lngr=dlog(gr) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          endif 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          if (g(i).lt.1e-6) then 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            lngrbond=0.01 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          else 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+            lngrbond=dlog(g(i)) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          endif 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          s2=s2-2*pi*rho*((gr*lngr)-gr+1)*del*r**2.0 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          s2bond=s2bond-2*pi*rho*((g(i)*lngrbond)-g(i)+1)*del*r*r 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+           
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          rf=dble(i-.5)*del 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+          write(23,*) rf,g(i) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      enddo 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call nvtxEndRange 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+ 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      write(24,*)"s2      : ",s2 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      write(24,*)"s2bond  : ",s2bond 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      call cpu_time(finish) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      print*,"starting at time",start,"and ending at",finish 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      stop 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+      deallocate(x,y,z,g) 
			 | 
		
	
		
			
				 | 
				 | 
			
			
				+end 
			 |