Hi every one
I have this code I have co developed,now I want to add a subroutine which calculated or do binning for different range of calculated angles. The code works fine but I struggling to add a subroutine which will do binning for this code. To make this easier I would like to do do binning at line 216.If some one can show me by example how to this on this code help will be much appreciated as I would like to use it in other codes in future.
[/ program adist
c*********************************************************************
c
c routine to calculate angle distributions in silica-based melts
c
c The following angles are supported:
c O-M-O
c M-O-M
c
c*********************************************************************
implicit real*8(a-h,o-z)
parameter (mxatms=2377)
character*1 fmt
character*8 atname(3)
character*40 history
character*80 title
character*8 atmnam(mxatms)
dimension xxx(mxatms),yyy(mxatms),zzz(mxatms)
dimension d(mxatms,mxatms)
dimension nbl(mxatms,mxatms)
dimension nb(mxatms), nnlist(mxatms)
dimension shell(2)
real aaa,bbb,ccc,pi,rtod
dimension weight(mxatms),chge(mxatms)
dimension avcell(9),cell(9),rcell(9)
logical new
data new/.true./,nhist/77/
parameter (pi=3.1415926536d0)
parameter (rtod=180.d0/pi)
natms=2377
nconf=500
history='HISTORY'
c list of coordinated atoms
open(9,file='clist.dat')
write(9,*)'atom type atom type dist'
c list of bridging and nonbridging oxygens
open(11,file='bonbo.dat')
write(11,*)'bridging nonbridging ratio'
c make sure all the arrays are .0d0 (real) or 0 (integer)
do iconf=1,nconf
do i = 1,mxatms
nb(i) = 0
do j = i,mxatms
d(i,j) = 0.0d0
nbl(i,j) = 0
enddo
enddo
if(new)then
open(nhist,file=history,status='old',err=100)
read(nhist,'(a80)',err=200) title
write(*,'(a,a)')'# History file header: ',title
read(nhist,'(2i10)',end=200) ktrj,imcon
new=.false.
endif
read(nhist,'(a8,4i10,f12.6)',end=200)
x step,nstep,matms,ktrj,imcon,timstp
c if(natms.ne.matms)then
c
c write(*,'(a)')'# error - incorrect number of atoms in file'
c write(*,'(a,i6,a)')'# file contains',matms,' atoms'
c stop
c
c endif
if(imcon.gt.0) read(nhist,'(3g12.4)',end=200) cell
if (imcon.ne.1) then
write(*,'(a)')'# error - Cubic PBC not used.'
stop
endif
do i = 1,natms
read(nhist,'(a8,i10,2f12.6)',end=200)
x atmnam(i),j,weight(i),chge(i)
read(nhist,'(1p,3e12.4)',end=200) xxx(i),yyy(i),zzz(i)
if(ktrj.ge.1)then
read(nhist,'(1p,3e12.4)',end=200) aaa,bbb,ccc
endif
if(ktrj.ge.2)then
read(nhist,'(1p,3e12.4)',end=200) aaa,bbb,ccc
endif
enddo
c do distance checking using minimum image convention
c see do i = 1,natms
do j = i,natms
dx = xxx(j) - xxx(i)
if (abs(dx) .gt. cell(1) * 0.5) then
dx = dx - sign(cell(1),dx)
endif
dy = yyy(j) - yyy(i)
if (abs(dy) .gt. cell(5) * 0.5) then
dy = dy - sign(cell(5),dy)
endif
dz = zzz(j) - zzz(i)
if (abs(dz) .gt. cell(9) * 0.5) then
dz = dz - sign(cell(9),dz)
endif
d(i,j) = sqrt(dx**2 + dy**2 + dz**2)
d(j,i) = d(i,j)
enddo
enddo
c this is a test. writes out first 5 atoms of each frame
c do i = 1,5
c write(7,'(5e12.4)') (d(i,j),j=1,5)
c enddo
c write(7,*)
c count neighbours within specified cutoff in each frame
c these variables should finally be input and not hard coded
atname(1) = 'Si'
atname(2) = 'Al'
atname(3) = 'O'
shell(1) = 2.275d0
shell(2) = 2.475d0
write(9,*)' Start Frame ',iconf,atname(1),atname(3)
do i = 1,natms
if (atmnam(i) .eq. atname(1)) then
do j = 1,natms
if (atmnam(j) .eq. atname(3)) then
if (d(i,j) .le. shell(1)) then
write(9,*)i,j,d(i,j)
nbl(i,j) = 1
c nb(i)=nb(i)+1
endif
endif
enddo
endif
enddo
write(9,*)' End Frame '
write(9,*)' Start Frame ',iconf,atname(2),atname(3)
do i = 1,natms
if (atmnam(i) .eq. atname(2)) then
do j = 1,natms
if (atmnam(j) .eq. atname(3)) then
if (d(i,j) .le. shell(2)) then
write(9,*)i,j,d(i,j)
nbl(i,j) = 1
c nb(i)=nb(i)+1
endif
endif
enddo
endif
enddo
write(9,*)' End Frame '
c counters/identifiers:
c i is the chosen central atom (Si,Al)
c j will count over neighbours, but is not used again after this
c k is the index in the neighbour list
c nnn is the number of nearest neighbours
c nnlist is the list of nearest neighbours
do i = 1,natms
if ((atmnam(i) .eq. atname(1)) .or.
x (atmnam(i) .eq. atname(2))) then
nnn = 0
k = 0
do j = 1, natms
if (nbl(i,j) .eq. 1) then
nnn = nnn + 1
k = k + 1
nnlist(k) = j
endif
enddo
write(*,*) 'test'
write(*,*) 'nnn',nnn
angsum = 0.0d0
do l = 1,nnn
do m = l+1,nnn
dxa = xxx(nnlist(l)) - xxx(i)
if (abs(dxa) .gt. cell(1) * 0.5) then
dxa = dxa - sign(cell(1),dxa)
endif
dya = yyy(nnlist(l)) - yyy(i)
if (abs(dya) .gt. cell(5) * 0.5) then
dya = dya - sign(cell(5),dya)
endif
dza = zzz(nnlist(l)) - zzz(i)
if (abs(dza) .gt. cell(9) * 0.5) then
dza = dza - sign(cell(9),dza)
endif
dxb = xxx(nnlist(m)) - xxx(i)
if (abs(dxb) .gt. cell(1) * 0.5) then
dxb = dxb - sign(cell(1),dxb)
endif
dyb = yyy(nnlist(m)) - yyy(i)
if (abs(dyb) .gt. cell(5) * 0.5) then
dyb = dyb - sign(cell(5),dyb)
endif
dzb = zzz(nnlist(m)) - zzz(i)
if (abs(dzb) .gt. cell(9) * 0.5) then
dzb = dzb - sign(cell(9),dzb)
endif
abdot = dxa*dxb + dya*dyb + dza*dzb
anorm = dxa**2 + dya**2 + dza**2
bnorm = dxb**2 + dyb**2 + dzb**2
ang = acos(abdot/(sqrt(anorm*bnorm)))
c ---------------------------------------------------
c DO BINNING HERE
c ---------------------------------------------------
angsum = angsum + ang
write(*,*) ang*rtod
enddo
enddo
angave = angsum/(0.5d0*(dble(nnn*(nnn-1))))
write(*,*)'average angle',angave*rtod
endif
enddo
c ibo = 0
c inbo = 0
c do i = 1,natms
c if (atmnam(i) .eq. atname(3)) then
c nsum = 0
c do j = 1, natms
c nsum = nsum + nbl(j,i)
c enddo
c if (nsum .ge. 2) then
c ibo = ibo + 1
c else
c inbo = inbo + 1
c endif
c endif
c enddo
c write(11,*)ibo,inbo,real(ibo)/real(inbo)
c END LOOP OVER FRAMES
enddo
stop
100 continue
write(*,'(a)')'# error - History file not found'
close (nhist)
stop
200 continue
write(*,'(a)')'# warning - end of History file encountered'
close (nhist)
end
subroutine simpson(n,del,sum,aaa)
c*********************************************************************
c
c dl_poly subroutine for performing integration using
c Simpson's rule.
c***********************************************************************
implicit real*8(a-h,o-z)
dimension aaa(*)
j=n
m=2*(n/2)
if(m.eq.n)j=n-1
sum=(aaa(1)+aaa(j))/2.d0
do i=2,j,2
sum=sum+2.d0*aaa(i)+aaa(i+1)
enddo
sum=2.d0*del*sum/3.d0
if(m.eq.n)sum=sum+del*(aaa(j)+aaa
)/2.d0
return
end]
I have this code I have co developed,now I want to add a subroutine which calculated or do binning for different range of calculated angles. The code works fine but I struggling to add a subroutine which will do binning for this code. To make this easier I would like to do do binning at line 216.If some one can show me by example how to this on this code help will be much appreciated as I would like to use it in other codes in future.
[/ program adist
c*********************************************************************
c
c routine to calculate angle distributions in silica-based melts
c
c The following angles are supported:
c O-M-O
c M-O-M
c
c*********************************************************************
implicit real*8(a-h,o-z)
parameter (mxatms=2377)
character*1 fmt
character*8 atname(3)
character*40 history
character*80 title
character*8 atmnam(mxatms)
dimension xxx(mxatms),yyy(mxatms),zzz(mxatms)
dimension d(mxatms,mxatms)
dimension nbl(mxatms,mxatms)
dimension nb(mxatms), nnlist(mxatms)
dimension shell(2)
real aaa,bbb,ccc,pi,rtod
dimension weight(mxatms),chge(mxatms)
dimension avcell(9),cell(9),rcell(9)
logical new
data new/.true./,nhist/77/
parameter (pi=3.1415926536d0)
parameter (rtod=180.d0/pi)
natms=2377
nconf=500
history='HISTORY'
c list of coordinated atoms
open(9,file='clist.dat')
write(9,*)'atom type atom type dist'
c list of bridging and nonbridging oxygens
open(11,file='bonbo.dat')
write(11,*)'bridging nonbridging ratio'
c make sure all the arrays are .0d0 (real) or 0 (integer)
do iconf=1,nconf
do i = 1,mxatms
nb(i) = 0
do j = i,mxatms
d(i,j) = 0.0d0
nbl(i,j) = 0
enddo
enddo
if(new)then
open(nhist,file=history,status='old',err=100)
read(nhist,'(a80)',err=200) title
write(*,'(a,a)')'# History file header: ',title
read(nhist,'(2i10)',end=200) ktrj,imcon
new=.false.
endif
read(nhist,'(a8,4i10,f12.6)',end=200)
x step,nstep,matms,ktrj,imcon,timstp
c if(natms.ne.matms)then
c
c write(*,'(a)')'# error - incorrect number of atoms in file'
c write(*,'(a,i6,a)')'# file contains',matms,' atoms'
c stop
c
c endif
if(imcon.gt.0) read(nhist,'(3g12.4)',end=200) cell
if (imcon.ne.1) then
write(*,'(a)')'# error - Cubic PBC not used.'
stop
endif
do i = 1,natms
read(nhist,'(a8,i10,2f12.6)',end=200)
x atmnam(i),j,weight(i),chge(i)
read(nhist,'(1p,3e12.4)',end=200) xxx(i),yyy(i),zzz(i)
if(ktrj.ge.1)then
read(nhist,'(1p,3e12.4)',end=200) aaa,bbb,ccc
endif
if(ktrj.ge.2)then
read(nhist,'(1p,3e12.4)',end=200) aaa,bbb,ccc
endif
enddo
c do distance checking using minimum image convention
c see do i = 1,natms
do j = i,natms
dx = xxx(j) - xxx(i)
if (abs(dx) .gt. cell(1) * 0.5) then
dx = dx - sign(cell(1),dx)
endif
dy = yyy(j) - yyy(i)
if (abs(dy) .gt. cell(5) * 0.5) then
dy = dy - sign(cell(5),dy)
endif
dz = zzz(j) - zzz(i)
if (abs(dz) .gt. cell(9) * 0.5) then
dz = dz - sign(cell(9),dz)
endif
d(i,j) = sqrt(dx**2 + dy**2 + dz**2)
d(j,i) = d(i,j)
enddo
enddo
c this is a test. writes out first 5 atoms of each frame
c do i = 1,5
c write(7,'(5e12.4)') (d(i,j),j=1,5)
c enddo
c write(7,*)
c count neighbours within specified cutoff in each frame
c these variables should finally be input and not hard coded
atname(1) = 'Si'
atname(2) = 'Al'
atname(3) = 'O'
shell(1) = 2.275d0
shell(2) = 2.475d0
write(9,*)' Start Frame ',iconf,atname(1),atname(3)
do i = 1,natms
if (atmnam(i) .eq. atname(1)) then
do j = 1,natms
if (atmnam(j) .eq. atname(3)) then
if (d(i,j) .le. shell(1)) then
write(9,*)i,j,d(i,j)
nbl(i,j) = 1
c nb(i)=nb(i)+1
endif
endif
enddo
endif
enddo
write(9,*)' End Frame '
write(9,*)' Start Frame ',iconf,atname(2),atname(3)
do i = 1,natms
if (atmnam(i) .eq. atname(2)) then
do j = 1,natms
if (atmnam(j) .eq. atname(3)) then
if (d(i,j) .le. shell(2)) then
write(9,*)i,j,d(i,j)
nbl(i,j) = 1
c nb(i)=nb(i)+1
endif
endif
enddo
endif
enddo
write(9,*)' End Frame '
c counters/identifiers:
c i is the chosen central atom (Si,Al)
c j will count over neighbours, but is not used again after this
c k is the index in the neighbour list
c nnn is the number of nearest neighbours
c nnlist is the list of nearest neighbours
do i = 1,natms
if ((atmnam(i) .eq. atname(1)) .or.
x (atmnam(i) .eq. atname(2))) then
nnn = 0
k = 0
do j = 1, natms
if (nbl(i,j) .eq. 1) then
nnn = nnn + 1
k = k + 1
nnlist(k) = j
endif
enddo
write(*,*) 'test'
write(*,*) 'nnn',nnn
angsum = 0.0d0
do l = 1,nnn
do m = l+1,nnn
dxa = xxx(nnlist(l)) - xxx(i)
if (abs(dxa) .gt. cell(1) * 0.5) then
dxa = dxa - sign(cell(1),dxa)
endif
dya = yyy(nnlist(l)) - yyy(i)
if (abs(dya) .gt. cell(5) * 0.5) then
dya = dya - sign(cell(5),dya)
endif
dza = zzz(nnlist(l)) - zzz(i)
if (abs(dza) .gt. cell(9) * 0.5) then
dza = dza - sign(cell(9),dza)
endif
dxb = xxx(nnlist(m)) - xxx(i)
if (abs(dxb) .gt. cell(1) * 0.5) then
dxb = dxb - sign(cell(1),dxb)
endif
dyb = yyy(nnlist(m)) - yyy(i)
if (abs(dyb) .gt. cell(5) * 0.5) then
dyb = dyb - sign(cell(5),dyb)
endif
dzb = zzz(nnlist(m)) - zzz(i)
if (abs(dzb) .gt. cell(9) * 0.5) then
dzb = dzb - sign(cell(9),dzb)
endif
abdot = dxa*dxb + dya*dyb + dza*dzb
anorm = dxa**2 + dya**2 + dza**2
bnorm = dxb**2 + dyb**2 + dzb**2
ang = acos(abdot/(sqrt(anorm*bnorm)))
c ---------------------------------------------------
c DO BINNING HERE
c ---------------------------------------------------
angsum = angsum + ang
write(*,*) ang*rtod
enddo
enddo
angave = angsum/(0.5d0*(dble(nnn*(nnn-1))))
write(*,*)'average angle',angave*rtod
endif
enddo
c ibo = 0
c inbo = 0
c do i = 1,natms
c if (atmnam(i) .eq. atname(3)) then
c nsum = 0
c do j = 1, natms
c nsum = nsum + nbl(j,i)
c enddo
c if (nsum .ge. 2) then
c ibo = ibo + 1
c else
c inbo = inbo + 1
c endif
c endif
c enddo
c write(11,*)ibo,inbo,real(ibo)/real(inbo)
c END LOOP OVER FRAMES
enddo
stop
100 continue
write(*,'(a)')'# error - History file not found'
close (nhist)
stop
200 continue
write(*,'(a)')'# warning - end of History file encountered'
close (nhist)
end
subroutine simpson(n,del,sum,aaa)
c*********************************************************************
c
c dl_poly subroutine for performing integration using
c Simpson's rule.
c***********************************************************************
implicit real*8(a-h,o-z)
dimension aaa(*)
j=n
m=2*(n/2)
if(m.eq.n)j=n-1
sum=(aaa(1)+aaa(j))/2.d0
do i=2,j,2
sum=sum+2.d0*aaa(i)+aaa(i+1)
enddo
sum=2.d0*del*sum/3.d0
if(m.eq.n)sum=sum+del*(aaa(j)+aaa
)/2.d0return
end]