implicit real*8(a-h,o-z) parameter (nat=2000) parameter (ntypd = 4) parameter (nshelld=100) parameter (natmd=10) dimension num(5) dimension a(3,3),natom(3),neigh(ntypd,natmd,ntypd,natmd) dimension r(ntypd,natmd,ntypd,natmd,nat),rr(nshelld) dimension no(ntypd,natmd,ntypd,natmd,nat),mo(nshelld) dimension itype(ntypd,natmd,ntypd,natmd,nat) dimension t(ntypd,natmd,3) character*4 fname(ntypd) open(22,file='neig.dat') write(6,*) ' enter the input file POSCAR/1 CONTCAR/2 ' read(5,*) ifile if (ifile .eq. 1) open(8,file='POSCAR') if (ifile .eq. 2) open(8,file='CONTCAR') write(6,*) 'enter the type of atom' read(5,*) ntype if (ntype .gt. ntypd) stop ' ntype too large ' write(6,*) ' enter the name of each atom (10a4) ' read(5,1) (fname(j),j=1,ntype) read(8,600) temp 600 format(20a4) read (8,*) aa c *** read lattice constant from POSCAR** do i=1,3 read (8,*) (a(i,j),j=1,3) 500 format (3f12.8) enddo do i=1,3 do j=1,3 a(i,j)=aa*a(i,j) enddo enddo read(8,*) (natom(i),i=1,ntype) c read(8,1) (fname(j),j=1,ntype) write(6,1) (fname(j),j=1,ntype) 1 format(10a4) do i = 1,ntype if (natom(i) .gt. natmd) stop ' natom too large ' enddo c write(6,*) ' natom(i) =',(natom(i),i=1,ntype) read(8,600) temp do i=1,ntype do j=1,natom(i) read(8,*) t(i,j,1),t(i,j,2),t(i,j,3) enddo enddo c *** read lattice vector from POSCAR*** do 200 i = 1,ntype do 200 j = 1,natom(i) x = t(i,j,1)*a(1,1)+t(i,j,2)*a(2,1)+t(i,j,3)*a(3,1) y = t(i,j,1)*a(1,2)+t(i,j,2)*a(2,2)+t(i,j,3)*a(3,2) z = t(i,j,1)*a(1,3)+t(i,j,2)*a(2,3)+t(i,j,3)*a(3,3) do 100 ip = 1,ntype do 100 jp = 1,natom(ip) nn = 0 do 50 n1=-5,5 do 50 n2=-5,5 do 50 n3=-5,5 dtij1 = t(i,j,1) - t(ip,jp,1) dtij2 = t(i,j,2) - t(ip,jp,2) dtij3 = t(i,j,3) - t(ip,jp,3) c write(6,*) ' test1 ',dtij1,dtij2,dtij3 xp = (n1+dtij1)*a(1,1) + (n2+dtij2)*a(2,1) + & (n3+dtij3)*a(3,1) yp = (n1+dtij1)*a(1,2) + (n2+dtij2)*a(2,2) + & (n3+dtij3)*a(3,2) zp = (n1+dtij1)*a(1,3) + (n2+dtij2)*a(2,3) + & (n3+dtij3)*a(3,3) c write(6,*) ' test2 ',xp,yp,zp dr = sqrt(xp**2+yp**2+zp**2) c write(6,1234) i,j,ip,jp,dr,n1,n2,n3,xp,yp,zp 1234 format(4i3,f10.3,3i3,3f10.3) c write(6,'(4f10.3)') xp,yp,zp,dr if ((dr .gt. 15.0) .or. (dr .eq. 0.0)) go to 50 nn = nn+ 1 if (nn .gt. nat) stop ' nn too large ' r(i,j,ip,jp,nn) = dr itype(i,j,ip,jp,nn) = ip 50 continue neigh(i,j,ip,jp) = nn c do n = 1,neigh(i,j,ip,jp) c write(6,*) i,j,ip,jp,r(i,j,ip,jp,n) c enddo c write(6,*) nn c write(17,*) i,j,ip,jp,neigh(i,j,ip,jp) 100 continue 200 continue c close(17) do 1000 i = 1,ntype do 1000 j = 1,natom(i) do 1000 ip = 1,ntype do 1000 jp = 1,natom(ip) c do n = 1,neigh(i,j,ip,jp) c write(6,*) ' test ',i,j,r(i,j,ip,jp,n),itype(i,j,ip,jp,n) c enddo c write(6,*) ' test4 ',i,j,ip,jp,neigh(i,j,ip,jp) do 60 m = 1,neigh(i,j,ip,jp)-1 do 60 n = m+1,neigh(i,j,ip,jp) if (r(i,j,ip,jp,n) .lt. r(i,j,ip,jp,m)) then dump = r(i,j,ip,jp,n) r(i,j,ip,jp,n) = r(i,j,ip,jp,m) r(i,j,ip,jp,m) = dump idump = itype(i,j,ip,jp,n) itype(i,j,ip,jp,n) = itype(i,j,ip,jp,m) itype(i,j,ip,jp,m) = idump endif 60 continue c do n = 1,neigh(i,j,ip,jp) c write(6,*) i,j,ip,jp,r(i,j,ip,jp,n) c enddo c write(6,*) ' test5 ' c write(6,*) ' test6 ' 1000 continue c do n=1,neigh(1,1,1,1) c write(6,*) r(1,1,1,1,n) c enddo do 2200 i = 1,ntype do 2200 j = 1,natom(i) write(22,111) fname(i),j do 2000 ip = 1,ntype do 2000 jp = 1,natom(ip) c write(6,*) ' i,j,ip,jp ',i,j,ip,jp,neigh(i,j,ip,jp) c write(13,*) ' i,j,ip,jp ',i,j,ip,jp c write(21,*) ' i,j,ip,jp ',i,j,ip,jp nshell = 1 do n= 2,neigh(i,j,ip,jp) if (abs(r(i,j,ip,jp,n)-r(i,j,ip,jp,n-1)) .gt. 0.0001) then nshell = nshell + 1 if(nshell .gt. nshelld) stop ' nshell too large ' rr(nshell) = r(i,j,ip,jp,n) no(i,j,ip,jp,n) = n mo(nshell) = no(i,j,ip,jp,n) c write(6,*) rr(nshell),mo(nshell) endif enddo rr(1) = r(i,j,ip,jp,(mo(2)-1)) mo(1) = mo(2)-1 mo(nshell+1) = neigh(i,j,ip,jp) c write(6,*) rr(1),mo(1),nshell,mo(nshell+1) c do ii = 1,nshell c write(21,*) ii,i,j,ip,jp,mo(ii),rr(ii) c write(6,*) ii,rr(ii) c enddo num(1) = mo(1) do m = 2,nshell-1 num(m) = mo(m+1)-mo(m) num(nshell) = mo(nshell+1)-mo(nshell)+1 c write(6,*) num(1),rr(1),num(m),rr(m),num(nshell),rr(nshell) 111 format(/,a4,'-',i1 & ,' ********************************************************* ') 112 format(2x,a4,'-',i1,7(f7.3,i3)) enddo write(22,112) fname(ip),jp,(rr(ii),num(ii),ii=1,7) 2000 continue 2200 continue stop end