1c X. W. Zhou, xzhou@sandia.gov
2      include 'createAtoms.h'
3      call latgen()
4      call delete()
5      call hitvalue()
6      call disturb()
7      call systemshift()
8      call colect()
9      if (ilatseed .ne. 0) call randomize()
10      call velgen()
11      call outputfiles()
12      if (float(nntype(1)) .ne. 0.0)
13     *   write(6,*)(float(nntype(i))/float(nntype(1)),i=1,ntypes)
14      write(6,*)'atoms of different species ',(nntype(i),i=1,ntypes)
15      stop
16      end
17cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
18      subroutine latgen()
19      include 'createAtoms.h'
20      character*8 lattype
21      common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3),
22     *   periodicity(3),xbound(2),ybound(2),zbound(2),strain(3),
23     *   delx(3),rcell(3),ccell(nelmax)
24      common /iran/ iseed
25      namelist /maincard/ ntypes,amass,ielement,ilatseed,
26     *   perlb,perub,iseed,xy,xz,yz
27      namelist /latcard/ lattype,alat,xrot,yrot,zrot,
28     *   periodicity,xbound,ybound,zbound,
29     *   strain,delx
30      iseed=3245566
31      natoms=0
32      ilatseed=0
33      xy=1.0e12
34      xz=1.0e12
35      yz=1.0e12
36      read(5,maincard)
37      do i = 1,3
38         perlen(i)=perub(i)-perlb(i)
39      enddo
40      do i=1,ntypes
41         nntype(i)=0
42      enddo
431     continue
44      strain(1)=0.0
45      strain(2)=0.0
46      strain(3)=0.0
47      xbound(1)=0.0
48      xbound(2)=0.0
49      ybound(1)=0.0
50      ybound(2)=0.0
51      zbound(1)=0.0
52      zbound(2)=0.0
53      delx(1)=0.0
54      delx(2)=0.0
55      delx(3)=0.0
56      lattype='none'
57      read(5,latcard)
58      if (lattype .eq. 'none') goto 2
59      call latgen1()
60      call defvalue()
61      goto 1
622     continue
63      return
64      end
65cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
66      subroutine latgen1()
67      include 'createAtoms.h'
68      character*8 lattype
69      common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3),
70     *   periodicity(3),xbound(2),ybound(2),zbound(2),strain(3),
71     *   delx(3),rcell(3),ccell(nelmax)
72      common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3),
73     *   avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100)
74      data xold/1.0,0.0,0.0/,yold/0.0,1.0,0.0/,zold/0.0,0.0,1.0/
75      namelist /subcard/ rcell,ccell
76
77      call storelat(lattype,avec,bvec,cvec)
78      if (xbound(1) .eq. xbound(2)) then
79         xbound(1)=perlb(1)
80         xbound(2)=perub(1)
81      endif
82      if (ybound(1) .eq. ybound(2)) then
83         ybound(1)=perlb(2)
84         ybound(2)=perub(2)
85      endif
86      if (zbound(1) .eq. zbound(2)) then
87         zbound(1)=perlb(3)
88         zbound(2)=perub(3)
89      endif
90      vol=alat(1)*alat(2)*alat(3)
91      xnorm=sqrt((xrot(1)*alat(2)*alat(3))**2+
92     *           (xrot(2)*alat(1)*alat(3))**2+
93     *           (xrot(3)*alat(1)*alat(2))**2)
94      ynorm=sqrt((yrot(1)*alat(2)*alat(3))**2+
95     *           (yrot(2)*alat(1)*alat(3))**2+
96     *           (yrot(3)*alat(1)*alat(2))**2)
97      znorm=sqrt((zrot(1)*alat(2)*alat(3))**2+
98     *           (zrot(2)*alat(1)*alat(3))**2+
99     *           (zrot(3)*alat(1)*alat(2))**2)
100      periodicity(1)=vol*periodicity(1)/xnorm
101      periodicity(2)=vol*periodicity(2)/ynorm
102      periodicity(3)=vol*periodicity(3)/znorm
103      do 5 i=1,3
104         xrot(i)=vol*xrot(i)/(alat(i)*xnorm)
105         yrot(i)=vol*yrot(i)/(alat(i)*ynorm)
106         zrot(i)=vol*zrot(i)/(alat(i)*znorm)
1075     continue
108      do 10 i=1,3
109      do 10 j=1,3
110         roter(i,j)=xold(i)*xrot(j)+yold(i)*yrot(j)+zold(i)*zrot(j)
11110    continue
112      do 20 i=1,3
113         avecp(i)=0.0
114         bvecp(i)=0.0
115         cvecp(i)=0.0
116      do 20 j=1,3
117         avecp(i)=avecp(i)+roter(i,j)*avec(j)*0.5*alat(j)
118         bvecp(i)=bvecp(i)+roter(i,j)*bvec(j)*0.5*alat(j)
119         cvecp(i)=cvecp(i)+roter(i,j)*cvec(j)*0.5*alat(j)
12020    continue
121      do 30 i=1,3
122         avec(i)=avecp(i)
123         bvec(i)=bvecp(i)
124         cvec(i)=cvecp(i)
12530    continue
126      call applystrain()
127      nrange=20
128      nsmall=0
129      do 300 ii=-nrange,nrange
130      do 300 jj=-nrange,nrange
131      do 300 kk=-nrange,nrange
132         xt=ii*avec(1)+jj*bvec(1)+kk*cvec(1)
133         if (xt .ge. periodicity(1)-0.1 .or. xt .lt. -0.1) goto 300
134         yt=ii*avec(2)+jj*bvec(2)+kk*cvec(2)
135         if (yt .ge. periodicity(2)-0.1 .or. yt .lt. -0.1) goto 300
136         zt=ii*avec(3)+jj*bvec(3)+kk*cvec(3)
137         if (zt .ge. periodicity(3)-0.1 .or. zt .lt. -0.1) goto 300
138         nsmall=nsmall+1
139         rvs(1,nsmall)=xt
140         rvs(2,nsmall)=yt
141         rvs(3,nsmall)=zt
142300   continue
143      nxa=xbound(1)/periodicity(1)-2
144      nxb=xbound(2)/periodicity(1)+2
145      nya=ybound(1)/periodicity(2)-2
146      nyb=ybound(2)/periodicity(2)+2
147      nza=zbound(1)/periodicity(3)-2
148      nzb=zbound(2)/periodicity(3)+2
1491     rcell(1)=-100.0
150      read(5,subcard)
151      if (rcell(1) .eq. -100.0) then
152         return
153      endif
154      do 190 m = 2,ntypes
155         ccell(m)=min(1.0,ccell(m-1)+ccell(m))
156190   continue
157      do 800 i=1,3
158         rcellp(i)=0.0
159      do 800 j=1,3
160         rcellp(i)=rcellp(i)+roter(i,j)*rcell(j)*alat(j)
161800   continue
162      do 801 i=1,3
163         rcell(i)=rcellp(i)*(1.0+strain(i))
164801   continue
165      do 250 nn=1,nsmall
166      do 251 ii=nxa,nxb
167         xt=ii*periodicity(1)+rvs(1,nn)+rcell(1)
168         if (xt .ge. xbound(2) .or. xt .lt. xbound(1)) goto 251
169      do 252 jj=nya,nyb
170         yt=jj*periodicity(2)+rvs(2,nn)+rcell(2)
171         if (yt .ge. ybound(2) .or. yt .lt. ybound(1)) goto 252
172      do 253 kk=nza,nzb
173         zt=kk*periodicity(3)+rvs(3,nn)+rcell(3)
174         if (zt .ge. zbound(2) .or. zt .lt. zbound(1)) goto 253
175         tst=ranl()
176         it=0
177      do 201 m=ntypes,1,-1
178         if (tst .le. ccell(m)) it=m
179201   continue
180         natoms=natoms+1
181         if (natoms .gt. natmax) then
182            write(6,9601)natmax
1839601        format('number of atoms exceeds maximum dimension: ',i10)
184         endif
185         itype(natoms)=it
186         nhittag(natoms)=0
187         nntype(it)=nntype(it)+1
188         rv(1,natoms)=xt+delx(1)
189         rv(2,natoms)=yt+delx(2)
190         rv(3,natoms)=zt+delx(3)
191         if (rv(1,natoms) .lt. perlb(1))
192     *      rv(1,natoms)=rv(1,natoms)+perlen(1)
193         if (rv(1,natoms) .ge. perub(1))
194     *      rv(1,natoms)=rv(1,natoms)-perlen(1)
195         if (rv(2,natoms) .lt. perlb(2))
196     *      rv(2,natoms)=rv(2,natoms)+perlen(2)
197         if (rv(2,natoms) .ge. perub(2))
198     *      rv(2,natoms)=rv(2,natoms)-perlen(2)
199         if (rv(3,natoms) .lt. perlb(3))
200     *      rv(3,natoms)=rv(3,natoms)+perlen(3)
201         if (rv(3,natoms) .ge. perub(3))
202     *      rv(3,natoms)=rv(3,natoms)-perlen(3)
203253   continue
204252   continue
205251   continue
206250   continue
207      goto 1
208      end
209cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
210      subroutine defvalue()
211      include 'createAtoms.h'
212      integer oldtype
213      character*8 lattype
214      common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3),
215     *   periodicity(3),xbound(2),ybound(2),zbound(2),strain(3),
216     *   delx(3),rcell(3),ccell(nelmax)
217      common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3),
218     *   avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100)
219      dimension cent(3),plane(3),planec(3),surface(3)
220      namelist /defcard/ xmin,xmax,ymin,ymax,zmin,zmax,
221     *   oldtype,newtype,prob,cent,plane,dist,surface,
222     *   nramp,ndirec,dismax,theta,phi,nscrew,anglemin,anglemax,
223     *   xbottom,xheight,zbottom,zheight,nrepeatsx,nrepeatsz
224      data pi/3.1415927/
225
2261     continue
227      xmin=-1.0e-20
228      xmax=1.0e20
229      ymin=-1.0e20
230      ymax=1.0e20
231      zmin=-1.0e20
232      zmax=1.0e20
233      oldtype=0
234      dist=-1.0
235      newtype=-11
236      cent(1)=0.5*(perlb(1)+perub(1))
237      cent(2)=0.5*(perlb(2)+perub(2))
238      cent(3)=0.5*(perlb(3)+perub(3))
239      surface(1)=0.0
240      surface(2)=0.0
241      surface(3)=0.0
242      nramp=1
243      nscrew=0
244      anglemin=0.0
245      anglemax=360.0
246      ndirec=3
247      dismax=0.0
248      theta=5000.0
249      phi=5000.0
250      nrepeatsx=0
251      nrepeatsz=0
252      read(5,defcard)
253      if (newtype .lt. -10 .and. dismax .eq. 0.0) return
254      if (newtype .gt. ntypes) then
255         ntypes=ntypes+1
256         if (oldtype .ne. 0) then
257            amass(ntypes)=amass(oldtype)
258         else
259            amass(ntypes)=amass(ntypes-1)
260         endif
261      endif
262      if (surface(1) .eq. 1.0) then
263         perlb(1)=perlb(1)-10.0
264         perub(1)=perub(1)+10.0
265         perlen(1)=perub(1)-perlb(1)
266      endif
267      if (surface(2) .eq. 1.0) then
268         perlb(2)=perlb(2)-10.0
269         perub(2)=perub(2)+10.0
270         perlen(2)=perub(2)-perlb(2)
271      endif
272      if (surface(3) .eq. 1.0) then
273         perlb(3)=perlb(3)-10.0
274         perub(3)=perub(3)+10.0
275         perlen(3)=perub(3)-perlb(3)
276      endif
277      vol=alat(1)*alat(2)*alat(3)
278      if (dismax .ne. 0.0) then
279         astart=anglemin
280         afinish=anglemax
281         if (anglemax .lt. anglemin) then
282            anglemin=afinish
283            anglemax=astart
284         endif
285         if (nramp .eq. 1) then
286            pstart=xmin
287            pfinish=xmax
288            if (xmax .lt. xmin) then
289               xmin=pfinish
290               xmax=pstart
291            endif
292         endif
293         if (nramp .eq. 2) then
294            pstart=ymin
295            pfinish=ymax
296            if (ymax .lt. ymin) then
297               ymin=pfinish
298               ymax=pstart
299            endif
300         endif
301         if (nramp .eq. 3) then
302            pstart=zmin
303            pfinish=zmax
304            if (zmax .lt. zmin) then
305               zmin=pfinish
306               zmax=pstart
307            endif
308         endif
309         do 33 i=1,natoms
310            if (rv(1,i) .lt. xmin) goto 33
311            if (rv(1,i) .gt. xmax) goto 33
312            if (rv(2,i) .lt. ymin) goto 33
313            if (rv(2,i) .gt. ymax) goto 33
314            if (rv(3,i) .lt. zmin) goto 33
315            if (rv(3,i) .gt. zmax) goto 33
316            if (nscrew .eq. 0) then
317               if (nramp .eq. 0) then
318                  rv(ndirec,i)=rv(ndirec,i)+dismax
319               else
320                  rv(ndirec,i)=rv(ndirec,i)+
321     *               (rv(nramp,i)-pstart)/(pfinish-pstart)*dismax
322               endif
323            else
324               if (ndirec .eq. 1) then
325                  dx=rv(2,i)-cent(2)
326                  dy=rv(3,i)-cent(3)
327                  dxmax=abs(ymax-cent(2))
328                  dxmin=abs(ymin-cent(2))
329                  rmax=min(dxmin,dxmax)
330               endif
331               if (ndirec .eq. 2) then
332                  dx=rv(3,i)-cent(3)
333                  dy=rv(1,i)-cent(1)
334                  dxmax=abs(zmax-cent(3))
335                  dxmin=abs(zmin-cent(3))
336                  rmax=min(dxmin,dxmax)
337               endif
338               if (ndirec .eq. 3) then
339                  dx=rv(1,i)-cent(1)
340                  dy=rv(2,i)-cent(2)
341                  dxmax=abs(xmax-cent(1))
342                  dxmin=abs(xmin-cent(1))
343                  rmax=min(dxmin,dxmax)
344               endif
345               dr=sqrt(dx**2+dy**2)
346               if (dr .eq. 0.0) then
347                  angle=0.0
348               else
349                  angle=acos(dx/dr)*180.0/pi
350                  if (dy .lt. 0.0) angle=360.0-angle
351               endif
352               if (angle .lt. anglemin) goto 33
353               if (angle .gt. anglemax) goto 33
354               if (nscrew .eq. 1) then
355                  astart1=astart
356                  afinish1=afinish
357               else if (nscrew .eq. 2) then
358                  if (astart .eq. 0.0 .and.
359     *               afinish .eq. 180.0) then
360                     if (dr .le. dxmax) then
361                        astart1=astart
362                     else
363                        astart1=acos(dxmax/dr)*180.0/pi
364                     endif
365                     if (dr .le. dxmin) then
366                        afinish1=afinish
367                     else
368                        afinish1=180.0-acos(dxmin/dr)*180.0/pi
369                     endif
370                  else if (astart .eq. 180.0 .and.
371     *               afinish .eq. 0.0) then
372                     if (dr .le. dxmin) then
373                        astart1=astart
374                     else
375                        astart1=180.0-acos(dxmin/dr)*180.0/pi
376                     endif
377                     if (dr .le. dxmax) then
378                        afinish1=afinish
379                     else
380                        afinish1=acos(dxmax/dr)*180.0/pi
381                     endif
382                  else if (astart .eq. 180.0 .and.
383     *               afinish .eq. 360.0) then
384                     if (dr .le. dxmin) then
385                        astart1=astart
386                     else
387                        astart1=180.0+acos(dxmin/dr)*180.0/pi
388                     endif
389                     if (dr .le. dxmax) then
390                        afinish1=afinish
391                     else
392                        afinish1=360.0-acos(dxmax/dr)*180.0/pi
393                     endif
394                  else if (astart .eq. 360.0 .and.
395     *               afinish .eq. 180.0) then
396                     if (dr .le. dxmax) then
397                        astart1=astart
398                     else
399                        astart1=360.0-acos(dxmax/dr)*180.0/pi
400                     endif
401                     if (dr .le. dxmin) then
402                        afinish1=afinish
403                     else
404                        afinish1=180.0+acos(dxmin/dr)*180.0/pi
405                     endif
406                  else
407                     write(6,*)'problem in nscrew=2'
408                     stop
409                  endif
410               endif
411               if (dr .gt. rmax) dr=rmax
412               rv(ndirec,i)=rv(ndirec,i)+dr/rmax*
413     *            (angle-astart1)/(afinish1-astart1)*dismax
414            endif
41533       continue
416      else if (theta .eq. 5000.0 .and. phi .eq. 5000.0 .and.
417     *   nrepeatsx*nrepeatsz .eq. 0.0) then
418         if (dist .lt. 0.0) then
419            do 2 i=1,natoms
420               if (rv(1,i) .lt. xmin) goto 2
421               if (rv(1,i) .gt. xmax) goto 2
422               if (rv(2,i) .lt. ymin) goto 2
423               if (rv(2,i) .gt. ymax) goto 2
424               if (rv(3,i) .lt. zmin) goto 2
425               if (rv(3,i) .gt. zmax) goto 2
426               if (oldtype .eq. itype(i) .or. oldtype .eq. 0) then
427                  if (ranl() .lt. prob) itype(i)=newtype
428               endif
4292           continue
430         else
431            if (theta .eq. 5000.0) then
432               pnorm=sqrt((plane(1)*alat(2)*alat(3))**2+
433     *                    (plane(2)*alat(1)*alat(3))**2+
434     *                    (plane(3)*alat(1)*alat(2))**2)
435               do 30 n=1,3
436                  plane(n)=vol*plane(n)/(alat(n)*pnorm)
43730             continue
438               do 20 n=1,3
439                  planec(n)=0.0
440               do 20 m=1,3
441                  planec(n)=planec(n)+roter(n,m)*plane(m)
44220             continue
443            else
444               planec(1)=sin(pi/180.0*theta)*cos(pi/180.0*phi)
445               planec(3)=sin(pi/180.0*theta)*sin(pi/180.0*phi)
446               planec(2)=cos(pi/180.0*theta)
447            endif
448            do 22 i=1,natoms
449               dis1=rv(1,i)-cent(1)
450               dis2=rv(2,i)-cent(2)
451               dis3=rv(3,i)-cent(3)
452               if (dis1*planec(1)+dis2*planec(2)+dis3*planec(3)
453     *         .gt. dist) then
454                  if (oldtype .eq. itype(i) .or. oldtype .eq. 0)
455     *               itype(i)=newtype
456               endif
45722          continue
458         endif
459      else if (nrepeatsx*nrepeatsz .ne. 0) then
460         do 4 i=1,natoms
461            if (rv(1,i) .le. xmin .or. rv(1,i) .ge. xmax) then
462               ylocalx = xbottom
463            else
464               ylocalx = xbottom + 0.5*xheight*(1.0+
465     *            cos(2.0*pi*nrepeatsx/(xmax-xmin)*
466     *            (rv(1,i)-(xmin+xmax)/(2.0*nrepeatsx))))
467            endif
468            if (rv(3,i) .le. zmin .or. rv(3,i) .ge. zmax) then
469               ylocalz = zbottom
470            else
471               ylocalz = zbottom + 0.5*zheight*(1.0+
472     *            cos(2.0*pi*nrepeatsz/(zmax-zmin)*
473     *            (rv(3,i)-(zmin+zmax)/(2.0*nrepeatsz))))
474            endif
475            ylocal = ylocalx
476            if (ylocal .gt. ylocalz) ylocal = ylocalz
477            if (rv(2,i) .gt. ylocal) then
478               itype(i)=newtype
479            endif
4804        continue
481      endif
482
483      goto 1
484      end
485cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
486      subroutine storelat(latty,avec,bvec,cvec)
487      character*8 lstored(10),latty
488      dimension avecs(3,10),bvecs(3,10),cvecs(3,10)
489      dimension avec(3),bvec(3),cvec(3)
490      data lstored/'fcc     ','bcc     ','sc      ',7*'        '/
491      data avecs/1.0,1.0,0.0, 1.0,1.0,-1.0, 2.0,0.0,0.0, 21*0.0/
492      data bvecs/0.0,1.0,1.0, -1.0,1.0,1.0, 0.0,2.0,0.0, 21*0.0/
493      data cvecs/1.0,0.0,1.0, 1.0,-1.0,1.0, 0.0,0.0,2.0, 21*0.0/
494      imatch=0
495      do 10 i=1,10
496         if (latty .ne. lstored(i)) goto 10
497         imatch = 1
498         do 5 j=1,3
499            avec(j)=avecs(j,i)
500            bvec(j)=bvecs(j,i)
501            cvec(j)=cvecs(j,i)
502 5       continue
50310    continue
504      if (imatch .eq. 1) return
505      write(6,15) latty
50615    format('could not find lattice type ')
507      stop
508      end
509cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
510      subroutine applystrain()
511      include 'createAtoms.h'
512      character*8 lattype
513      common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3),
514     *   periodicity(3),xbound(2),ybound(2),zbound(2),strain(3),
515     *   delx(3),rcell(3),ccell(nelmax)
516      common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3),
517     *   avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100)
518      do i=1,3
519         periodicity(i)=(1.0+strain(i))*periodicity(i)
520         avec(i)=(1.0+strain(i))*avec(i)
521         bvec(i)=(1.0+strain(i))*bvec(i)
522         cvec(i)=(1.0+strain(i))*cvec(i)
523      enddo
524      return
525      end
526cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
527      subroutine hitvalue()
528      include 'createAtoms.h'
529      namelist /hitcard/ xmin,xmax,ymin,ymax,zmin,zmax,nhit
530      nhitcards=0
5311001  nhit=0
532      xmin=-10000.0
533      xmax=10000.0
534      ymin=-10000.0
535      ymax=10000.0
536      zmin=-10000.0
537      zmax=10000.0
538      read(5,hitcard)
539      if (nhit .eq. 0) then
540         return
541      endif
542      nchange=0
543      do 1002 i=1,natoms
544         if (rv(1,i) .lt. xmin) goto 1002
545         if (rv(1,i) .gt. xmax) goto 1002
546         if (rv(2,i) .lt. ymin) goto 1002
547         if (rv(2,i) .gt. ymax) goto 1002
548         if (rv(3,i) .lt. zmin) goto 1002
549         if (rv(3,i) .gt. zmax) goto 1002
550         nchange=nchange+1
551         nhittag(i)=nhit
5521002  continue
553      if (nchange .ne. 0) nhitcards=nhitcards+1
554      goto 1001
555      end
556cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
557      subroutine systemshift()
558      include 'createAtoms.h'
559      namelist /shiftcard/ mode
560      mode=0
561      read(5,shiftcard)
562      if (mode .eq. 0) return
563      if (mode .eq. 1) then
564         top=perub(2)
565         bottom=perlb(2)
566         sdely=0.0
567         do i=1,natoms
568            if (nhittag(i) .eq. 3 .and. top .gt. rv(2,i)) top=rv(2,i)
569            if (nhittag(i) .eq. 4 .and. bottom .lt. rv(2,i))
570     *         bottom=rv(2,i)
571         enddo
572         sdely=-0.5*(top+bottom)
573         do i=1,natoms
574            rv(2,i)=rv(2,i)+sdely
575         enddo
576      endif
577      if (mode .eq. 2) then
578c     This shift only assumes a maximum of two plane spacings.
579         xmin=10.0
580         do i=1,natoms
581            if (rv(1,i) .lt. xmin) xmin=rv(1,i)
582         enddo
583         shiftneg=-10.0
584         shiftpos=10.0
585         do i=1,natoms
586            shift=rv(1,i)-xmin
587            shift=shift-perlen(1)*nint(shift/perlen(1))
588            if (shift .gt. 0.01 .and. shift .lt. shiftpos)
589     *         shiftpos=shift
590            if (shift .lt. -0.01 .and. shift .gt. shiftneg)
591     *         shiftneg=shift
592         enddo
593         if (-shiftneg .gt. shiftpos) then
594            shift=0.5*shiftneg
595         else
596            shift=0.5*shiftpos
597         endif
598         perlb(1)=xmin+shift
599         perub(1)=perlb(1)+perlen(1)
600      endif
601      end
602cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
603      subroutine colect()
604      include 'createAtoms.h'
605      do i=1,natoms
606         if (rv(1,i) .lt. perlb(1)) rv(1,i)=rv(1,i)+perlen(1)
607         if (rv(1,i) .ge. perub(1)) rv(1,i)=rv(1,i)-perlen(1)
608         if (rv(2,i) .lt. perlb(2)) rv(2,i)=rv(2,i)+perlen(2)
609         if (rv(2,i) .ge. perub(2)) rv(2,i)=rv(2,i)-perlen(2)
610         if (rv(3,i) .lt. perlb(3)) rv(3,i)=rv(3,i)+perlen(3)
611         if (rv(3,i) .ge. perub(3)) rv(3,i)=rv(3,i)-perlen(3)
612      enddo
613      return
614      end
615cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
616      real function ranl()
617      common /iran/ iseed
618      iseed=iseed*125
619      iseed=iseed-2796203*(iseed/2796203)
620      ranl=abs(iseed/2796203.0)
621      return
622      end
623cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
624      subroutine delete()
625      include 'createAtoms.h'
626      nw_del=0
627      do i=1,ntypes
628         nntype(i)=0
629      enddo
630      natoms0=natoms
631      natoms=0
632      do 1 i=1,natoms0
633         if (itype(i) .eq. -1) nw_del=nw_del+1
634         if (itype(i) .le. 0) goto 1
635         natoms=natoms+1
636         ntag(natoms)=natoms
637         itype(natoms)=itype(i)
638         nhittag(natoms)=nhittag(i)
639         nntype(itype(natoms))=nntype(itype(natoms))+1
640         rv(1,natoms)=rv(1,i)
641         rv(2,natoms)=rv(2,i)
642         rv(3,natoms)=rv(3,i)
6431     continue
644      return
645      end
646cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
647      subroutine randomize()
648      include 'createAtoms.h'
649      dimension rv0(6,natmax),itype0(natmax),nhittag0(natmax)
650      do 1 i=1,natoms
651         itype0(i)=itype(i)
652         if (nhitcards .gt. 0) nhittag0(i)=nhittag(i)
653         rv0(1,i)=rv(1,i)
654         rv0(2,i)=rv(2,i)
655         rv0(3,i)=rv(3,i)
656         rv0(4,i)=rv(4,i)
657         rv0(5,i)=rv(5,i)
658         rv0(6,i)=rv(6,i)
6591     continue
660      im=1
66110    im=im*2
662      if (im .lt. natoms) goto 10
663      ia=365
664      ic=8161
665      ivalue=mod(ilatseed,natoms)
666      do 2 i=1,natoms
66720       ivalue=mod(ia*ivalue+ic,im)
668         if (ivalue .ge. natoms) goto 20
669         ntmp=ivalue+1
670         if (ntmp .ge. 1 .and. ntmp .le. natoms) then
671            ntag(ntmp)=i
672            itype(ntmp)=itype0(i)
673            if (nhitcards .gt. 0) nhittag(ntmp)=nhittag0(i)
674            rv(1,ntmp)=rv0(1,i)
675            rv(2,ntmp)=rv0(2,i)
676            rv(3,ntmp)=rv0(3,i)
677            rv(4,ntmp)=rv0(4,i)
678            rv(5,ntmp)=rv0(5,i)
679            rv(6,ntmp)=rv0(6,i)
680         endif
6812     continue
682      return
683      end
684cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
685      subroutine write0(dumpfile)
686      include 'createAtoms.h'
687      character*32 dumpfile, typ
688      open (unit=2,file=dumpfile)
689      write(2,*) natoms
690      write(2,100) 'GaN lattice'
691      do i=1,natoms
692         if (itype(i).eq.1) then
693            typ = 'Ga'
694         else
695            typ = 'N'
696         endif
697         write(2,150) typ(1:length(typ)),rv(1,i),rv(2,i),rv(3,i)
698      enddo
699 50   format(i8)
700 100  format(a)
701 150  format(a3,3f13.5)
702      close(2)
703      return
704      end
705cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
706      subroutine write1(dumpfile)
707      include 'createAtoms.h'
708      character*32 dumpfile
709      zero=0.0
710      open (unit=2,file=dumpfile)
711      write(2,*) 'ITEM: NUMBER OF ATOMS'
712      write(2,*) natoms
713      write(2,*) 'ITEM: BOX BOUNDS'
714      write(2,*) perlb(1),perub(1)
715      write(2,*) perlb(2),perub(2)
716      write(2,*) perlb(3),perub(3)
717      write(2,*) 'ITEM: TIME'
718      write(2,*) 0.0
719      write(2,*) 'ITEM: ATOMS'
720      if (nhitcards .eq. 0) then
721         do i=1,natoms
722            write(2,*) ntag(i),itype(i),rv(1,i),rv(2,i),rv(3,i)
723         enddo
724      else
725         do i=1,natoms
726            write(2,*) ntag(i),itype(i),nhittag(i),
727     *         rv(1,i),rv(2,i),rv(3,i)
728         enddo
729      endif
730      close(2)
731      open(unit=2,file=dumpfile(1:index(dumpfile,' ')-1)//'.vel')
732      write(2,*) 'ITEM: NUMBER OF ATOMS'
733      write(2,*) natoms
734      write(2,*) 'ITEM: BOUNDARY VELOCITY'
735      write(2,*) zero,zero,zero
736      write(2,*) 'ITEM: TIME'
737      write(2,*) zero
738      write(2,*) 'ITEM: VELOCITIES'
739      do i=1,natoms
740         write(2,*) ntag(i),rv(4,i),rv(5,i),rv(6,i)
741      enddo
742      close(2)
743      return
744      end
745cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
746      subroutine write2(dumpfile)
747      include 'createAtoms.h'
748      data conmas/1.0365e-4/
749      character*32 dumpfile
750      open (unit=2,file=dumpfile)
751      zero=0.0
752      write(2,*) 'using dynamo'
753      write(2,9502) natoms,ntypes,zero
7549502  format(2i10,e15.8)
755      write(2,9503) (perub(i),i=1,3),(perlb(i),i=1,3)
7569503  format(3e25.16)
757      write(2,9504) (amass(i)*conmas,ielement(i),i=1,ntypes)
7589504  format(e25.16,i10)
759      write(2,9505) ((rv(i,j),i=1,6),itype(j),j=1,natoms)
7609505  format(3e25.16/3e25.16/i10)
761      write(2,9503) zero,zero,zero
762      close(2)
763      return
764      end
765cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
766      subroutine write3(dumpfile)
767      include 'createAtoms.h'
768      character*32 dumpfile
769      open (unit=2,file=dumpfile)
770      area=(1.0-float(nw_del)/float(natoms0))*perlen(2)*perlen(3)
771      write(2,*) 'Start File for LAMMPS ',area
772      write(2,*)
773      write(2,*) natoms,' atoms'
774      write(2,*)
775      write(2,*) ntypes,' atom types'
776      write(2,*)
777      write(2,*) perlb(1),perub(1),' xlo xhi'
778      write(2,*) perlb(2),perub(2),' ylo yhi'
779      write(2,*) perlb(3),perub(3),' zlo zhi'
780      if (xy .ge. 1.0e8 .or. xz .ge. 1.0e8 .or. yz .ge. 1.0e8) then
781         write(2,*)
782      else
783         write(2,*)
784         write(2,*)xy,xz,yz,' xy xz yz'
785         write(2,*)
786      endif
787      write(2,*) 'Masses'
788      write(2,*)
789      do i=1,ntypes
790         write(2,*)i,amass(i)
791      enddo
792      write(2,*)
793      write(2,*) 'Atoms'
794      write(2,*)
795      do i=1,natoms
796         write(2,*) i,itype(i),rv(1,i),rv(2,i),rv(3,i)
797      enddo
798      write(2,*)
799      write(2,*) 'Velocities'
800      write(2,*)
801      zero=0.0
802      do i=1,natoms
803         write(2,*) i,rv(4,i),rv(5,i),rv(6,i)
804      enddo
805      close(2)
806      return
807      end
808cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
809      subroutine write4(dumpfile)
810      include 'createAtoms.h'
811      character*32 dumpfile,full
812      dumpfile='ens'
813      full=dumpfile(1:index(dumpfile,' ')-1)//'.case'
814      open(unit=2,file=full)
815      write(2,101)
816101   format('FORMAT')
817      write(2,102)
818102   format('type: ensight')
819      write(2,103)
820103   format('GEOMETRY')
821      write(2,104)'ens.case.****.geo'
822104   format('model: 1                     ',a30)
823      write(2,105)
824105   format('VARIABLE')
825      write(2,111)
826111   format('TIME')
827      write(2,112)
828112   format('time set: 1')
829      write(2,113)
830113   format('number of steps: 1')
831      write(2,114)
832114   format('filename start number: 1')
833      write(2,115)
834115   format('filename increment: 1')
835      write(2,116)
836116   format('time values:')
837      write(2,117)
838117   format('1')
839      close(2)
840      full=dumpfile(1:index(dumpfile,' ')-1)//'.case.0001.geo'
841      open(unit=2,file=full)
842      write(2,201)
843      write(2,202)
844      write(2,203)
845      write(2,204)
846      write(2,205)
847      write(2,206)natoms
848      do j=1,natoms
849         write(2,207)j,rv(1,j),rv(2,j),rv(3,j)
850      enddo
851      do j=1,ntypes
852         write(2,208)j
853         write(2,209)j
854         write(2,210)
855         write(2,206)nntype(j)
856         do k=1,natoms
857            if (itype(k) .eq. j) write(2,211)k,k
858         enddo
859      enddo
860201   format('spparks input')
861202   format('geometry for element group 1')
862203   format('node id given')
863204   format('element id given')
864205   format('coordinates')
865206   format(i8)
866207   format(i8,3e12.5)
867208   format('part',i2)
868209   format('atombox',i2)
869210   format('point')
870211   format(2i8)
871      close(2)
872      return
873      end
874cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
875      subroutine write5(dumpfile)
876      include 'createAtoms.h'
877      character*32 dumpfile,full
878      character*200 sentence
879      full=dumpfile(1:index(dumpfile,' ')-1)//'.lat'
880      open (unit=2,file=full)
881         call neigen()
882         nmax=numneigh(1)
883         do i=1,natoms
884            if (nmax .lt. numneigh(i)) nmax=numneigh(i)
885         enddo
886         write(2,*) 'spparks input lattice'
887         write(2,*)
888         write(2,*) '3 dimension'
889         write(2,*) natoms,' vertices'
890         write(2,*) nmax,' max connectivity'
891         write(2,*) perlb(1),perub(1),' xlo xhi'
892         write(2,*) perlb(2),perub(2),' ylo yhi'
893         write(2,*) perlb(3),perub(3),' zlo zhi'
894         write(2,*)
895         write(2,*) 'Vertices'
896         write(2,*)
897         do i=1,natoms
898            write(2,*) i,rv(1,i),rv(2,i),rv(3,i)
899         enddo
900         write(2,*)
901         write(2,*) 'Edges'
902         write(2,*)
903         do i=1,natoms
904            write(sentence,*) i,(neigh(k,i),k=1,numneigh(i))
905            write(2,*)sentence(1:len_trim(sentence))
906         enddo
907      close(2)
908      full=dumpfile(1:index(dumpfile,' ')-1)//'.conf'
909      open (unit=2,file=full)
910      nvalues=1
911      write(2,*) 'spparks input configuration'
912      write(2,*) natoms,nvalues
913      write(2,*)
914      do i=1,natoms
915         write(2,*) i,itype(i)
916      enddo
917      close(2)
918      return
919      end
920cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
921      subroutine disturb()
922      include 'createAtoms.h'
923      dimension strain(3)
924      namelist /disturbcard/ dismax,strain
925      dismax=0.0
926      strain(1)=0.0
927      strain(2)=0.0
928      strain(3)=0.0
929      read(5,disturbcard)
930      if (dismax .ne. 0.0) then
931         do 1 i=1,natoms
932            rv(1,i)=rv(1,i)+dismax*(ranl()-0.5)
933            rv(2,i)=rv(2,i)+dismax*(ranl()-0.5)
934            rv(3,i)=rv(3,i)+dismax*(ranl()-0.5)
9351        continue
936      endif
937      if (strain(1)**2+strain(2)**2+strain(3)**2 .ne. 0.0) then
938         do 2 i=1,3
939            perlen(i)=perlen(i)*(1.0+strain(i))
940            perub(i)=perlb(i)+perlen(i)
9412        continue
942         do 3 i=1,natoms
943         do 3 j=1,3
944            rv(j,i)=perlb(j)+(rv(j,i)-perlb(j))*(1.0+strain(j))
9453        continue
946      endif
947      return
948      end
949cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
950      subroutine outputfiles()
951      character*32 dynamo,paradyn,lammps,xyz,ensight,spparks
952      namelist /filecard/ dynamo,paradyn,lammps,xyz,ensight,spparks
953      xyz="none"
954      paradyn="none"
955      dynamo="none"
956      lammps="none"
957      ensight="none"
958      spparks="none"
959      read(5,filecard)
960      if (xyz .ne. "none") call write0(xyz)
961      if (paradyn .ne. "none") call write1(paradyn)
962      if (dynamo .ne. "none") call write2(dynamo)
963      if (lammps .ne. "none") call write3(lammps)
964      if (ensight .ne. "none") call write4(ensight)
965      if (spparks .ne. "none") call write5(spparks)
966      return
967      end
968cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
969      integer function length(str)
970      character*(*) str
971      n = len(str)
972 10   if (n.gt.0) then
973        if (str(n:n).eq.' ') then
974          n = n - 1
975          goto 10
976        endif
977      endif
978      length = n
979      return
980      end
981cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
982      subroutine velgen()
983      include 'createAtoms.h'
984      data boltz/8.617e-5/,conmas/1.0365e-4/
985      namelist /velcard/ Tbnd,Tmid,naxis,nmesh,equal,ensureTave,
986     *   steeper
987      dimension Tmesh(200),pmesh(200),vcm(3,200),tmass(200),ekin(200),
988     *   scale(200),ncount(200)
989      Tbnd=0.0
990      Tmid=0.0
991      naxis=1
992      nmesh=100
993      equal=1.0
994      ensureTave=0.0
995      steeper=0.0
996      read(5,velcard)
997      Tbnd=equal*Tbnd
998      Tmid=equal*Tmid
999      Tave=0.5*(Tbnd+Tmid)
1000      Tbnd=Tave+(Tbnd-Tave)*(1.0+steeper)
1001      Tmid=Tave+(Tmid-Tave)*(1.0+steeper)
1002
1003      if (Tbnd+Tmid .le. 0.0) then
1004         do i=1,natoms
1005            rv(4,i)=0.0
1006            rv(5,i)=0.0
1007            rv(6,i)=0.0
1008         enddo
1009      else
1010         wmesh=perlen(naxis)/nmesh
1011         do i=1,nmesh
1012            pmesh(i)=(i-0.5)*wmesh
1013            Tmesh(i)=Tmid+(Tbnd-Tmid)*abs(pmesh(i)-0.5*perlen(naxis))/
1014     *         (0.5*perlen(naxis))
1015            vcm(1,i)=0.0
1016            vcm(2,i)=0.0
1017            vcm(3,i)=0.0
1018            tmass(i)=0.0
1019            ekin(i)=0.0
1020            ncount(i)=0
1021         enddo
1022         do i=1,natoms
1023            index=(rv(naxis,i)-perlb(naxis))/wmesh+1
1024            it=itype(i)
1025            vnorm=sqrt(Tmesh(index)*boltz/(amass(it)*conmas))
1026            do j=1,3
1027               rv(j+3,i)=vnorm*rgauss()
1028               vcm(j,index)=vcm(j,index)+amass(it)*conmas*rv(j+3,i)
1029            enddo
1030            tmass(index)=tmass(index)+amass(it)*conmas
1031            ncount(index)=ncount(index)+1
1032         enddo
1033         do i=1,nmesh
1034            do j=1,3
1035               vcm(j,i)=vcm(j,i)/tmass(i)
1036            enddo
1037         enddo
1038         do i=1,natoms
1039            index=(rv(naxis,i)-perlb(naxis))/wmesh+1
1040            do j=1,3
1041               rv(j+3,i)=rv(j+3,i)-vcm(j,index)
1042            enddo
1043         enddo
1044         do i=1,natoms
1045            index=(rv(naxis,i)-perlb(naxis))/wmesh+1
1046            it=itype(i)
1047            ekin(index)=ekin(index)+0.5*amass(it)*conmas*
1048     *         (rv(4,i)**2+rv(5,i)**2+rv(6,i)**2)
1049         enddo
1050         do i=1,nmesh
1051            treal=2.0/(3.0*boltz)*ekin(i)/ncount(i)
1052            if (treal .eq. 0.0) then
1053               scale(i)=0.0
1054            else
1055               scale(i)=sqrt(Tmesh(i)/treal)
1056            endif
1057         enddo
1058         do i=1,natoms
1059            index=(rv(naxis,i)-perlb(naxis))/wmesh+1
1060            do j=4,6
1061               rv(j,i)=rv(j,i)*scale(index)
1062            enddo
1063         enddo
1064         if (ensureTave .eq. 1.0) then
1065            ekint=0.0
1066            do i=1,natoms
1067               it=itype(i)
1068               ekint=ekint+0.5*amass(it)*conmas*
1069     *            (rv(4,i)**2+rv(5,i)**2+rv(6,i)**2)
1070            enddo
1071            trealt=2.0/(3.0*boltz)*ekint/natoms
1072            if (trealt .eq. 0.0) then
1073               scalet=0.0
1074            else
1075               scalet=sqrt(Tave/trealt)
1076            endif
1077            do i=1,natoms
1078               do j=4,6
1079                  rv(j,i)=rv(j,i)*scalet
1080               enddo
1081            enddo
1082         endif
1083      endif
1084      return
1085      end
1086cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1087      function rgauss()
1088      data twopi/6.283185308/
10891     continue
1090      a1 = ranl()
1091      if (a1 .eq. 0.0) goto 1
1092      a2 = ranl()
1093      rgauss = sqrt(-2.0*log(a1))*cos(twopi*a2)
1094      return
1095      end
1096cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1097      subroutine neigen()
1098      include 'createAtoms.h'
1099      dimension idbox(100,100,100,100),nbox(100,100,100)
1100      cutoff=100.0
1101      xmin=rv(1,1)
1102      xmax=rv(1,1)
1103      ymin=rv(2,1)
1104      ymax=rv(2,1)
1105      zmin=rv(3,1)
1106      zmax=rv(3,1)
1107      do i=2,natoms
1108         if (xmin .gt. rv(1,i)) xmin=rv(1,i)
1109         if (xmax .lt. rv(1,i)) xmax=rv(1,i)
1110         if (ymin .gt. rv(2,i)) ymin=rv(2,i)
1111         if (ymax .lt. rv(2,i)) ymax=rv(2,i)
1112         if (zmin .gt. rv(3,i)) zmin=rv(3,i)
1113         if (zmax .lt. rv(3,i)) zmax=rv(3,i)
1114         dx=rv(1,i)-rv(1,1)
1115         dy=rv(2,i)-rv(2,1)
1116         dz=rv(3,i)-rv(3,1)
1117         dx=dx-perlen(1)*nint(dx/perlen(1))
1118         dy=dy-perlen(2)*nint(dy/perlen(2))
1119         dz=dz-perlen(3)*nint(dz/perlen(3))
1120         tmp=dx**2+dy**2+dz**2
1121         if (tmp .lt. cutoff) cutoff=tmp
1122      enddo
1123      cutoff=sqrt(cutoff)+0.1
1124      nx=(xmax-xmin)/cutoff+1
1125      ny=(ymax-ymin)/cutoff+1
1126      nz=(zmax-zmin)/cutoff+1
1127      if (nx .gt. 100) nx=100
1128      if (ny .gt. 100) ny=100
1129      if (nz .gt. 100) nz=100
1130      delx=(xmax-xmin)/nx
1131      dely=(ymax-ymin)/ny
1132      delz=(zmax-zmin)/nz
1133      do n1=1,nx
1134      do n2=1,ny
1135      do n3=1,nz
1136         nbox(n1,n2,n3)=0
1137      enddo
1138      enddo
1139      enddo
1140      do j=1,natoms
1141         numneigh(j)=0
1142         n1=(rv(1,j)-xmin)/delx+1
1143         if (n1 .lt. 1) n1=1
1144         if (n1 .gt. nx) n1=nx
1145         n2=(rv(2,j)-ymin)/dely+1
1146         if (n2 .lt. 1) n2=1
1147         if (n2 .gt. ny) n2=ny
1148         n3=(rv(3,j)-zmin)/delz+1
1149         if (n3 .lt. 1) n3=1
1150         if (n3 .gt. nz) n3=nz
1151         nbox(n1,n2,n3)=nbox(n1,n2,n3)+1
1152         idbox(n1,n2,n3,nbox(n1,n2,n3))=j
1153      enddo
1154      do j=1,natoms
1155         n1=(rv(1,j)-xmin)/delx+1
1156         if (n1 .lt. 1) n1=1
1157         if (n1 .gt. nx) n1=nx
1158         n2=(rv(2,j)-ymin)/dely+1
1159         if (n2 .lt. 1) n2=1
1160         if (n2 .gt. ny) n2=ny
1161         n3=(rv(3,j)-zmin)/delz+1
1162         if (n3 .lt. 1) n3=1
1163         if (n3 .gt. nz) n3=nz
1164      do na=n1-1,n1+1
1165         if (na .lt. 1) then
1166            naa=nx
1167         else if (na .gt. nx) then
1168            naa=1
1169         else
1170            naa=na
1171         endif
1172      do nb=n2-1,n2+1
1173         if (nb .lt. 1) then
1174            nbb=ny
1175         else if (nb .gt. ny) then
1176            nbb=1
1177         else
1178            nbb=nb
1179         endif
1180      do nc=n3-1,n3+1
1181         if (nc .lt. 1) then
1182            ncc=nz
1183         else if (nc .gt. nz) then
1184            ncc=1
1185         else
1186            ncc=nc
1187         endif
1188      do 100 ndd=1,nbox(naa,nbb,ncc)
1189         k=idbox(naa,nbb,ncc,ndd)
1190         if (k .ge. j) goto 100
1191         dx=rv(1,k)-rv(1,j)
1192         dy=rv(2,k)-rv(2,j)
1193         dz=rv(3,k)-rv(3,j)
1194         dx=dx-perlen(1)*nint(dx/perlen(1))
1195         dy=dy-perlen(2)*nint(dy/perlen(2))
1196         dz=dz-perlen(3)*nint(dz/perlen(3))
1197         tmp=sqrt(dx**2+dy**2+dz**2)
1198         if (tmp .le. cutoff) then
1199            numneigh(j)=numneigh(j)+1
1200            neigh(numneigh(j),j)=k
1201            numneigh(k)=numneigh(k)+1
1202            neigh(numneigh(k),k)=j
1203         endif
1204100   continue
1205      enddo
1206      enddo
1207      enddo
1208      enddo
1209      return
1210      end
1211