Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- subroutine gr(switch)
- if (switch.eq.0) then %switch = 0 initialiazation
- ngr=0 %switch = 1 sample
- delg =box/(2*nihs) %switch =2 results
- do i= 0,nihs
- g(i)=0
- end do
- else if (switch.eq.1)then %sample
- ngr=ngr+1
- do i=1,npart-1
- do j=i+1,npart-1 %loop all over pairs
- xr=x(i)-x(j)
- xr= xr-box*niht(xr/box) %periodic boundary conditions
- r =sqrt(xr**2)
- if(r.lt.box/2)then % only within half of the box
- ig=int(r/delg)
- g(ig)=g(ig)+2 %contribution for particle i and j
- endif
- enddo
- enddo
- else if (switch.eq.2)then % determine g(r)
- do i=1,nihs
- r=delg*(i+0.5) % distance r
- vb=((i+1)**3-i**3)delg**3 %volume betweeb bin i+1 and i
- nid = (4/3)*pi*vb*rho
- g(i)=g(i)/(ngr*npart*nid) % normalize g(r)
- enddo
- endif
- return
- end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement