SHARE
TWEET

Untitled

a guest Nov 23rd, 2019 122 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. c ------------------------------------------------------------
  2.       subroutine nbr (ia, na, ad, ncon, nrr, nr, rc, rk, n, rmax, mc)
  3. c  routine to supply nearest neighbour data for atoms in
  4. c  a crystal structure, given
  5. c  rc(i,j)  the i'th coordinate of the j'th axis of the unit cell
  6. c  rk(i,j)  the i'th coordinate of the j'th atom in the unit cell
  7. c  nrr(ir)  the number of type-ir atoms in the unit cell
  8. c
  9. c  The information returned, for a type-ir atom, is
  10. c  ncon(ir) the number of nearest neighbour shells of a type-ir
  11. c           atom included, out to a distance of rmax, but <= mc
  12. c  ia(j,ir) the type of atoms in the j'th neighbouring shell
  13. c  na(j,ir) the number of atoms in the j'th shell
  14. c  ad(j,ir) the radius of the j'th shell initialisation
  15.  
  16.      dimension ia(mc,nr), na(mc,nr), ad(mc,nr), ncon(nr), nrr(nr)
  17.      dimension rc(3,3), rk(3,n), rj(3)
  18.      
  19.      double precision as, ax, ay, az, r, dr
  20.  
  21.      rad(a1,a2,a3) = sqrt(a1*a1 + a2*a2 + a3*a3)
  22.      rcmin = 1.0e6
  23.      do i = 1, 3
  24.        rcmin = amin1(rcmin, rad(rc(1,i), rc(2,i), rc(3,i)))
  25.      end do
  26.      ia = 0
  27.      na = 0
  28.      ad = 1.0e6
  29. c  search over adjacent unit cells to include mc nearest neighbours
  30.      itc = ifix(rmax / rcmin) + 1
  31.      limc = itc + itc + 1
  32.      as = -float(itc + 1)
  33.      ax = as
  34.      do 100 jx = 1, limc
  35.        ax = ax + 1.0d0
  36.        ay = as
  37.        do 100 jy = 1, limc
  38.        ay = ay + 1.0d0
  39.        az = as
  40.        do 100 jz = 1, limc
  41.        az = az + 1.0d0
  42.        do j = 1, 3
  43.          rj(j) = ax*rc(j,1) + ay*rc(j,2) + az*rc(j,3)
  44.        end do
  45. c  rj is current unit cell origin. For each atom in this unit cell
  46. c  find displacement r from kr-type atom in basic unit cell
  47.        j = 0
  48.        do 100 jr = 1, nr
  49.          jnr = nrr(jr)
  50.          do 100 jjr = 1, jnr
  51.            j = j + 1
  52.            k = 1
  53.            do 90 kr = 1, nr
  54.              r = rad(rj(1) + rk(1,j) - rk(1,k),
  55.     &                rj(2) + rk(2,j) - rk(2,k),
  56.     &                rj(3) + rk(3,j) - rk(3,k))
  57.              if (r > rmax) goto 90
  58. c  compare r with nearest neighbour distances already found
  59.              ic = 0
  60.   40         ic = ic + 1
  61.              if (ic > mc) goto 90
  62.              dr = r - ad(ic,kr)
  63.              if (abs(dr) < 1.0e-4) dr = 0.0d0
  64.              if (dr) 60, 50, 40
  65.   50         if (ia(ic,kr) /= jr) goto 40
  66.              na(ic,kr) = na(ic,kr) + 1
  67.              goto 90
  68.   60         if (ic == mc) goto 80
  69.              iic = ic + 1
  70.              do jjc = iic, mc
  71.                jc = mc + iic - jjc
  72.                ia(jc,kr) = ia(jc - 1,kr)
  73.                na(jc,kr) = na(jc - 1,kr)
  74.                ad(jc,kr) = ad(jc - 1,kr)
  75.              end do
  76.   80         ia(ic,kr) = jr
  77.            na(ic,kr) = 1
  78.            ad(ic,kr) = r
  79.   90       k = k + nrr(kr)
  80.  100 continue
  81.      do 120 ir = 1, nr
  82.        ncon(ir) = 0
  83.        do ic = 1, mc
  84.          if (na(ic,ir) == 0) goto 120
  85.          ncon(ir) = ncon(ir) + 1
  86.        end do
  87.  120 continue
  88.      return
  89.      end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top