| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159 |
- 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 cudafor
- 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)
- !x=0;y=0;z=0
- !g=0
- 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);
- call nvtxStartRange("Pair Calculation")
- do iconf=1,nframes
- if (mod(iconf,1).eq.0) print*,iconf
- do concurrent(i=1 : natoms, 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
- !$acc atomic
- g(ind)=g(ind)+1.0d0
- !g(ind) = atomicadd(g(ind),1.0d0)
- endif
- enddo
- enddo
- call nvtxEndRange
-
- 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
- !print*,i+del,i,rupper,rlower,r
- !print*,r
- 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
|
