• API
• FAQ
• Tools
• Archive
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
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
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)
75.              end do
76.   80         ia(ic,kr) = jr
77.            na(ic,kr) = 1
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.
Top