mirror of https://gitlab.com/QEF/q-e.git
335 lines
11 KiB
Fortran
335 lines
11 KiB
Fortran
!
|
|
! Copyright (C) 2003 PWSCF group
|
|
! This file is distributed under the terms of the
|
|
! GNU General Public License. See the file `License'
|
|
! in the root directory of the present distribution,
|
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
|
!
|
|
!----------------------------------------------------------------------
|
|
program dist
|
|
!----------------------------------------------------------------------
|
|
!
|
|
! find distances, nearest neighbors, angles
|
|
! taking into account periodicity
|
|
!
|
|
! Input file:
|
|
! first line contains nat, ibrav, celldm(1-6)
|
|
! if (ibrav.eq.0) :
|
|
! card CELL_PARAMETERS
|
|
! a(1,1) a(2,1) a(3,1)
|
|
! a(1,2) a(2,2) a(3,2)
|
|
! a(1,3) a(2,3) a(3,3)
|
|
! [ alat units if celldm(1) was specified, a.u. if celldm(1)=0 ]
|
|
! end if
|
|
! card ATOMIC_POSITIONS
|
|
! cards with atomic positions
|
|
! See documentation of pw.x input for more details
|
|
!
|
|
implicit none
|
|
integer nspx, ndistx, nax, nnx
|
|
parameter (nspx=6, ndistx=1000, nax=200, nnx=4)
|
|
integer nat, ibrav, nsp, ityp(nax)
|
|
integer nsp1, nsp2, nat0, na, nb, n, nd, nn, nn1, nn2, nn3
|
|
integer iout, i
|
|
integer atom1(ndistx), atom2(ndistx), index(ndistx), ndist
|
|
character*3 atm(nspx), atm1, atm2
|
|
character*80 filename, line
|
|
character*1 capital, other_cell(ndistx)
|
|
real*8 at(3,3), bg(3,3), celldm(6), omega, d(ndistx)
|
|
real*8 tau(3,nax), dr(3), dd, dn1, dn2, dn3, dmin, dmax, scale
|
|
real*8 angolo(nnx*(nnx-1)/2), drv(3), drn(3,nnx), temp, rtemp(3)
|
|
real*8 fact, pi
|
|
parameter (fact=0.529177d0, pi=3.141592653589793d0)
|
|
logical crys, matches
|
|
external capital, matches
|
|
!
|
|
!
|
|
call get_file ( filename )
|
|
open(unit=1,file=filename,form='formatted',status='old')
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Begin data reading
|
|
read(1,*) nat, ibrav, celldm
|
|
if (nat.le.0.or.nat.gt.nax) stop ' nat!'
|
|
if (ibrav.lt.0.or.ibrav.gt.14) stop ' ibrav!'
|
|
if (ibrav.eq.0) then
|
|
! read cell unit vectors from input
|
|
read(1,'(a)') line
|
|
do i=1,len(line)
|
|
line(i:i) = capital(line(i:i))
|
|
end do
|
|
if (matches('CELL_PARAMETERS',line)) then
|
|
READ (1,*) ( ( at(i,n), i=1,3 ), n=1,3 )
|
|
else
|
|
stop ' expecting card CELL_PARAMETERS'
|
|
end if
|
|
!
|
|
if (celldm (1).eq.0.d0) then
|
|
! input at are in atomic units: define alat
|
|
celldm (1) = sqrt(at(1,1)**2+at(1,2)**2+at(1,3)**2)
|
|
do i=1,3
|
|
at(i,1) = at(i,1)/celldm(1)
|
|
at(i,2) = at(i,2)/celldm(1)
|
|
at(i,3) = at(i,3)/celldm(1)
|
|
end do
|
|
end if
|
|
else
|
|
! generate cell unit vectors
|
|
call latgen(ibrav,celldm,at(1,1),at(1,2),at(1,3),omega)
|
|
do i=1,3
|
|
at(i,1) = at(i,1)/celldm(1)
|
|
at(i,2) = at(i,2)/celldm(1)
|
|
at(i,3) = at(i,3)/celldm(1)
|
|
end do
|
|
end if
|
|
call recips(at(1,1),at(1,2),at(1,3),bg(1,1),bg(1,2),bg(1,3))
|
|
!
|
|
read(1,'(a)') line
|
|
do i=1,len(line)
|
|
line(i:i) = capital(line(i:i))
|
|
end do
|
|
if (matches('ATOMIC_POSITIONS',line)) then
|
|
if ( matches('ALAT', line) ) then
|
|
scale = 1.d0
|
|
crys = .false.
|
|
else if ( matches('BOHR', line) ) then
|
|
scale = celldm(1)
|
|
crys = .false.
|
|
else if ( matches('CRYSTAL', line) ) then
|
|
scale = 1.d0
|
|
crys = .true.
|
|
else if ( matches('ANGSTROM', line) ) then
|
|
scale = celldm(1) * fact
|
|
crys = .false.
|
|
else
|
|
scale = 1.d0
|
|
crys = .false.
|
|
end if
|
|
else
|
|
read(line,*) scale
|
|
end if
|
|
if (scale.le.0.d0 .or. scale.gt.1000.d0) stop ' scale ! '
|
|
|
|
nsp = 0
|
|
do na=1,nat
|
|
read(1,*) atm1, (tau(i,na),i=1,3)
|
|
do n = 1, nsp
|
|
if ( atm1 .eq. atm(n) ) then
|
|
ityp(na) = n
|
|
go to 10
|
|
end if
|
|
end do
|
|
nsp = nsp + 1
|
|
if (nsp .gt. nspx) stop ' too many atom types!'
|
|
atm(nsp) = atm1
|
|
ityp(na) = nsp
|
|
10 continue
|
|
enddo
|
|
if (crys) then
|
|
call cryst_to_cart ( nat , tau, at, 1)
|
|
else
|
|
do na=1,nat
|
|
do i=1,3
|
|
tau(i,na)=tau(i,na)/scale
|
|
end do
|
|
end do
|
|
end if
|
|
close(unit=1)
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! End data reading
|
|
print '(''output file (CR=output) > '',$)'
|
|
read(5,'(a)',end=20,err=20) filename
|
|
if (filename.eq.' ') then
|
|
iout=6
|
|
else
|
|
iout=1
|
|
open(unit=1,file=filename,form='formatted')
|
|
end if
|
|
!
|
|
scale=fact*celldm(1)
|
|
30 continue
|
|
if (nsp.gt.1) then
|
|
do n = 1, nsp
|
|
print '(''species # '',i1,'' : '', a3)', n, atm(n)
|
|
end do
|
|
print '(''indices of species 1 and 2 > '',$)'
|
|
read(5,*,end=20,err=20) nsp1, nsp2
|
|
if (nsp1 .le. 0 .or. nsp2 .le. 0 .or.
|
|
& nsp1 .gt. nspx .or. nsp2 .gt. nspx) then
|
|
print '('' wrong indices (type control-D to exit)'')'
|
|
go to 30
|
|
end if
|
|
else
|
|
nsp1 = nsp
|
|
nsp2 = nsp
|
|
end if
|
|
atm1=atm(nsp1)
|
|
atm2=atm(nsp2)
|
|
print '(''min and max distance > '',$)'
|
|
read(5,*,end=20,err=20) dmin, dmax
|
|
if (dmin.ge.dmax.or.dmax.le.0.0) go to 20
|
|
!
|
|
!!! dmin=max(dmin,0.01)
|
|
ndist=0
|
|
do na=1,nat
|
|
if (atm(ityp(na)) .eq. atm1) then
|
|
|
|
if (atm1.eq.atm2) then
|
|
nat0=na+1
|
|
else
|
|
nat0=1
|
|
end if
|
|
do nb=nat0,nat
|
|
if (atm(ityp(nb)) .eq. atm2) then
|
|
|
|
do i=1,3
|
|
dr(i) = (tau(1,na)-tau(1,nb))*bg(1,i) +
|
|
& (tau(2,na)-tau(2,nb))*bg(2,i) +
|
|
& (tau(3,na)-tau(3,nb))*bg(3,i)
|
|
enddo
|
|
do nn1=-2,2
|
|
dn1=dr(1)-nn1
|
|
do nn2=-2,2
|
|
dn2=dr(2)-nn2
|
|
do nn3=-2,2
|
|
dn3=dr(3)-nn3
|
|
dd = scale* sqrt(
|
|
& ( dn1*at(1,1)+dn2*at(1,2)+dn3*at(1,3) )**2+
|
|
& ( dn1*at(2,1)+dn2*at(2,2)+dn3*at(2,3) )**2+
|
|
& ( dn1*at(3,1)+dn2*at(3,2)+dn3*at(3,3) )**2)
|
|
if(dd.ge.dmin.and.dd.le.dmax) then
|
|
ndist=ndist+1
|
|
if (ndist.gt.ndistx) stop ' ndist !'
|
|
atom1(ndist)=na
|
|
atom2(ndist)=nb
|
|
d(ndist)= dd
|
|
if (nn1.eq.0.and.nn2.eq.0.and.nn3.eq.0)
|
|
& then
|
|
other_cell(ndist)=' '
|
|
else
|
|
other_cell(ndist)='*'
|
|
end if
|
|
end if
|
|
end do
|
|
end do
|
|
end do
|
|
end if
|
|
end do
|
|
end if
|
|
end do
|
|
!
|
|
index(1)=0.0
|
|
if (ndist.gt.0) call hpsort(ndist,d,index)
|
|
!
|
|
if (iout.eq.1)
|
|
& write(iout,100) atm1, atm2, dmin,dmax
|
|
do nd=1,ndist
|
|
write(iout,200) atom1(index(nd)), atom2(index(nd)),
|
|
& other_cell(index(nd)), d(nd)
|
|
end do
|
|
!
|
|
go to 30
|
|
!
|
|
20 nn = nnx
|
|
print '(/''number of neighbors (max '',i1,'') > '', $)', nnx
|
|
read(5,*,end=21,err=21) n
|
|
nn = n
|
|
if (nn.gt.nnx) stop ' too many neighbors'
|
|
21 continue
|
|
!
|
|
! look for nearest neighbors
|
|
!
|
|
do na=1,nat
|
|
!!! if (atm(ityp(na)) .eq. atm1) then
|
|
!
|
|
! ndist (.le.nnx) keep tracks of how many neighbors have been found
|
|
!
|
|
ndist=0
|
|
do nd=1,nn
|
|
d(nd)=100000.0
|
|
do i=1,3
|
|
drn(i,nd)=0.0
|
|
end do
|
|
end do
|
|
do nb=1,nat
|
|
do i=1,3
|
|
dr(i)=(tau(1,na)-tau(1,nb))*bg(1,i) +
|
|
& (tau(2,na)-tau(2,nb))*bg(2,i) +
|
|
& (tau(3,na)-tau(3,nb))*bg(3,i)
|
|
end do
|
|
do nn1=-1,1
|
|
dn1=dr(1)-nn1
|
|
do nn2=-1,1
|
|
dn2=dr(2)-nn2
|
|
do nn3=-1,1
|
|
dn3=dr(3)-nn3
|
|
dd = scale* sqrt(
|
|
& ( dn1*at(1,1)+dn2*at(1,2)+dn3*at(1,3) )**2 +
|
|
& ( dn1*at(2,1)+dn2*at(2,2)+dn3*at(2,3) )**2 +
|
|
& ( dn1*at(3,1)+dn2*at(3,2)+dn3*at(3,3) )**2 )
|
|
do i=1,3
|
|
drv(i) = tau(i,na)-tau(i,nb) -
|
|
& (nn1*at(i,1)+nn2*at(i,2)+nn3*at(i,3))
|
|
end do
|
|
!
|
|
! the "first" neighbor is the atom itself
|
|
!
|
|
if (dd.gt.0.01) then
|
|
! straight insertion: look for first nn neighbors
|
|
do nd=1,nn
|
|
if (dd.lt.d(nd)) then
|
|
! swap d(nd) with dd
|
|
temp = d(nd)
|
|
d(nd)= dd
|
|
dd = temp
|
|
! do the same for delta r
|
|
do i=1,3
|
|
rtemp(i) = drn(i,nd)
|
|
end do
|
|
do i=1,3
|
|
drn(i,nd) = drv(i)
|
|
end do
|
|
do i=1,3
|
|
drv(i) = rtemp(i)
|
|
end do
|
|
!
|
|
ndist=min(ndist+1,nn)
|
|
end if
|
|
end do
|
|
end if
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
!
|
|
if (ndist.ne.nn) stop ' what?'
|
|
!
|
|
! calculate angles with nearest neighbors
|
|
!
|
|
nd=0
|
|
do nn1=1,nn
|
|
do nn2=nn1+1,nn
|
|
nd=nd+1
|
|
angolo(nd) = 360/(2*pi) * acos (scale**2 *
|
|
& ( drn(1,nn1)*drn(1,nn2) +
|
|
& drn(2,nn1)*drn(2,nn2) +
|
|
& drn(3,nn1)*drn(3,nn2) ) / d(nn1) / d(nn2) )
|
|
end do
|
|
end do
|
|
if (nd.ne.nn*(nn-1)/2) stop ' what??'
|
|
!
|
|
! dd is the distance from the origin
|
|
!
|
|
dd = sqrt(tau(1,na)**2 + tau(2,na)**2 + tau(3,na)**2)*scale
|
|
write(iout,250) atm(ityp(na)), na, (d(nn1),nn1=1,nn)
|
|
write(iout,300) dd, (angolo(nn1),nn1=1,nn*(nn-1)/2)
|
|
!!! end if
|
|
end do
|
|
!
|
|
stop
|
|
!
|
|
100 format(' species: ',a3,' - ',a3,3x,f6.2,' < D <',f6.2)
|
|
200 format(' atoms:', 2i6,a1, ' distance =',f10.5,' A')
|
|
250 format(a3,i3,': neighbors at ',4f8.3,' A')
|
|
300 format(9x,'d(center):',f6.3,' A angles :',6f8.1)
|
|
!
|
|
!
|
|
end
|