Skip to content

WIP: Improve analyze_movie #38

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
danielhollas wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

WIP: Improve analyze_movie #38

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
92 changes: 65 additions & 27 deletions UTILS/analyze_movie.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,8 @@
! Simple program for analysis of bonds, angles and dihedrals
! from xyz trajectories.
! We pull out time evolution as well as histograms.
! from XYZ trajectories.
! Writes out time evolution as well as histograms.

! Daniel Hollas 2014

! 2016, RDF funcionality added
! Author: Daniel Hollas

module mod_analyze
implicit none
Expand All @@ -14,7 +12,7 @@ module mod_analyze

CONTAINS

! The following function should respect the PBC conditions
! The following function also handles the PBC conditions
real*8 function get_distance( at1, at2, boxx, boxy, boxz)
implicit none
integer,intent(in) :: at1, at2
Expand All @@ -29,8 +27,8 @@ real*8 function get_distance( at1, at2, boxx, boxy, boxz)
dx = x(at1) - x(at2)
dy = y(at1) - y(at2)
dz = z(at1) - z(at2)

end if

if (boxx.gt.0)then
dx = dx - boxx * nint(dx/boxx)
dy = dy - boxy * nint(dy/boxy)
Expand Down Expand Up @@ -154,12 +152,6 @@ program analyze_movie
call Get_cmdline(chmovie, nbin_dist, nbin_ang, nbin_dih, distmin, distmax, shiftdih, lecho, &
imod, boxx, boxy, boxz)

! the following formats were not a good idea and are now obsolete
10 format(I3)
20 format(2I3)
30 format(3I3)
40 format(4I3)

! Now read from stdin, what to analyze
select case (imod)
case (0)
Expand Down Expand Up @@ -195,7 +187,7 @@ program analyze_movie
end if

if (lecho) write(*,*)'How many dihedrals?'
read(*, 10,IOSTAT = iost) ndih
read(*, '(I3)', IOSTAT = iost) ndih
if (iost.ne.0) call PrintInputError()

if (ndih.gt.0)then
Expand Down Expand Up @@ -361,6 +353,10 @@ program analyze_movie

r(idist) = get_distance(dists(1,idist), dists(2,idist), boxx, boxy, boxz)

! TODO: Put histograming at the end
! so we can calculate optimal bin width based on Scott's rule
! https://en.wikipedia.org/wiki/Histogram#Number_of_bins_and_width
! and can calculate max and min data
ipom = ceiling( ( (r(idist)) - distmin )/dx )
if(ipom.gt.nbin_dist.or.ipom.le.0)then
write(*,*)'Problems with bond distribution function'
Expand Down Expand Up @@ -477,7 +473,7 @@ program analyze_movie
anorm=0.0d0

do ian=1,nbin_ang
anorm=anorm+bins_ang(ian,idist)
anorm = anorm + bins_ang(ian,idist)
enddo

do ian=1,nbin_ang
Expand Down Expand Up @@ -533,24 +529,18 @@ program analyze_movie

end program








subroutine Get_cmdline(chmovie, nbin_dist, nbin_ang, nbin_dih, distmin, distmax, shiftdih, lecho, &
imod, boxx, boxy, boxz)
implicit none
real*8,intent(inout) :: distmin, distmax, boxx, boxy, boxz
real*8,intent(inout) :: shiftdih
real*8,intent(inout) :: distmin, distmax, boxx, boxy, boxz
real*8,intent(inout) :: shiftdih
integer,intent(inout) :: nbin_dist,nbin_ang,nbin_dih
integer, intent(inout) :: imod
character(len=100),intent(out) :: chmovie
character(len=100) :: arg
integer :: i
logical, intent(inout) :: lecho
character(len=100),intent(out) :: chmovie
logical, intent(inout) :: lecho
character(len=100) :: arg
integer :: i, ngeom

i=0
do while (i < command_argument_count())
Expand Down Expand Up @@ -587,6 +577,12 @@ subroutine Get_cmdline(chmovie, nbin_dist, nbin_ang, nbin_dih, distmin, distmax,
i=i+1
call get_command_argument(i, arg)
read(arg,*)distmax
case ('-optimal_nbin')
i=i+1
call get_command_argument(i, arg)
read(arg,*)ngeom
call print_bin_number_estimators(ngeom)
stop 0
case ('-box')
i=i+1
call get_command_argument(i, arg)
Expand Down Expand Up @@ -642,4 +638,46 @@ subroutine PrintInputError()
stop 1
end subroutine PrintInputError

! Various estimates to number of bins
! See https://en.wikipedia.org/wiki/Histogram#Number_of_bins_and_width
function square_root_bin_number(n) result(k)
integer, intent(in) :: n
integer :: k
k = ceiling(sqrt(REAL(n)))
end function square_root_bin_number

function sturges_bin_number(n) result(k)
integer, intent(in) :: n
integer :: k
k = ceiling(log(REAL(n))/log(2.0d0)) + 1
end function sturges_bin_number

function rice_bin_number(n) result(k)
integer, intent(in) :: n
integer :: k
k = ceiling(2*(REAL(n)**(1/3.d0)))
end function rice_bin_number

function scott_bin_width(n, sigma) result(h)
integer, intent(in) :: n
real*8, intent(in) :: sigma
real*8 :: h
h = 3.5d0 * sigma / n**(1/3.0d0)
end function scott_bin_width

! TODO Freedman-Diaconis choice


subroutine print_bin_number_estimators(ngeom)
integer, intent(in) :: ngeom
integer :: square_root_bin_number, sturges_bin_number, rice_bin_number

write(*,'(A)')'Printing various estimators for optimal number of bins'
write(*,'(A)')'Square root estimate'
write(*,'(A4,I5)')'k = ', square_root_bin_number(ngeom)
write(*,'(A)')'Sturges estimate'
write(*,'(A4,I5)')'k = ', sturges_bin_number(ngeom)
write(*,'(A)')'Rice estimate'
write(*,'(A4,I5)')'k = ', rice_bin_number(ngeom)
end subroutine print_bin_number_estimators