1      subroutine argos_space_trvl(xw,vw,xwcr,gw,iwl,iwlp,
2     + xs,vs,gs,isl,islp,
3     + boxsiz,ibownr,ipl,ndx,itmp,rtmp,lenx)
4c
5      implicit none
6c
7#include "argos_space_common.fh"
8#include "global.fh"
9#include "msgids.fh"
10c
11      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3)
12      real*8 xs(msa,3),vs(msa,3)
13      real*8 gw(mwm,3,mwa),gs(msa,3)
14      integer iwl(mwm,miw2),iwlp(mwm,npackw)
15      integer isl(msa,mis2),islp(msa,npack)
16      integer lenx
17      real*8 boxsiz(maxbox,3)
18      integer ipl(mbox,mip2),ibownr(maxbox,3)
19      integer ndx(lenx),itmp(lenx)
20      real*8 rtmp(lenx)
21      logical lrec(27)
22c
23      integer i,indexw,indexs,j,k,ibx,iby,ibz,ipx,ipy,ipz
24      integer isbox,isnod,nrbox
25      integer ilp,ihp,jlp,jhp
26      integer il,ih,jl,jh,ilw,ihw,jlw,jhw,ils,ihs,jls,jhs
27      integer iliw,ihiw,jliw,jhiw
28      integer ilis,ihis,jlis,jhis
29      integer iwm,iwstay,jwstay,lwstay,nwgo,nwgosm
30      integer nwgtsm
31      integer nwnew,nwstay,iwmloc,jwmloc,lwmloc,irw
32      integer isa,jsa,isstay,jsstay,icsgm,ifsgm,ilsgm
33      integer nsnew,nsstay,lsstay,isaloc,jsaloc,lsaloc,irs
34      integer nsgo,jnode,iwfr,iwto,isfr,isto
35      real*8 factor,xscx,xscy,xscz,boxi(3)
36      integer itemps,nfold
37      logical lend
38      character*255 string
39c
40      boxi(1)=one/box(1)
41      boxi(2)=one/box(2)
42      boxi(3)=one/box(3)
43      nfold=0
44      lpbc9=.false.
45c
46      nwstay=0
47c
48c     order the solvent molecules
49c
50      if(nwmloc.gt.0) then
51      do 1 i=1,nwmloc
52      ndx(i)=i
53    1 continue
54      endif
55      if(nwmloc.gt.1) then
56      lwmloc=nwmloc/2+1
57      irw=nwmloc
58    2 continue
59      if(lwmloc.gt.1) then
60      lwmloc=lwmloc-1
61      itemps=ndx(lwmloc)
62      else
63      itemps=ndx(irw)
64      ndx(irw)=ndx(1)
65      irw=irw-1
66      if(irw.eq.1) then
67      ndx(1)=itemps
68      goto 3
69      endif
70      endif
71      iwmloc=lwmloc
72      jwmloc=lwmloc+lwmloc
73    4 continue
74      if(jwmloc.le.irw) then
75      if(jwmloc.lt.irw) then
76      if((iwl(ndx(jwmloc),lwnod).eq.iwl(ndx(jwmloc+1),lwnod).and.
77     + iwl(ndx(jwmloc),lwbox).le.iwl(ndx(jwmloc+1),lwbox)).or.
78     + ((iwl(ndx(jwmloc),lwnod).eq.me.or.
79     + (iwl(ndx(jwmloc),lwnod).ne.me.and.
80     + iwl(ndx(jwmloc),lwnod).le.iwl(ndx(jwmloc+1),lwnod))).and.
81     + iwl(ndx(jwmloc+1),lwnod).ne.me)) jwmloc=jwmloc+1
82      endif
83      if((iwl(itemps,lwnod).eq.iwl(ndx(jwmloc),lwnod).and.
84     + iwl(itemps,lwbox).le.iwl(ndx(jwmloc),lwbox)).or.
85     + ((iwl(itemps,lwnod).eq.me.or. (iwl(itemps,lwnod).ne.me.and.
86     + iwl(itemps,lwnod).le.iwl(ndx(jwmloc),lwnod))).and.
87     + iwl(ndx(jwmloc),lwnod).ne.me)) then
88      ndx(iwmloc)=ndx(jwmloc)
89      iwmloc=jwmloc
90      jwmloc=jwmloc+jwmloc
91      else
92      jwmloc=irw+1
93      endif
94      goto 4
95      endif
96      ndx(iwmloc)=itemps
97      goto 2
98    3 continue
99c
100      do 5 k=1,3
101      do 8 i=1,nwmloc
102      rtmp(i)=xwcr(i,k)
103    8 continue
104      do 9 i=1,nwmloc
105      xwcr(i,k)=rtmp(ndx(i))
106    9 continue
107      do 10 j=1,nwa
108      do 11 i=1,nwmloc
109      rtmp(i)=xw(i,k,j)
110   11 continue
111      do 12 i=1,nwmloc
112      xw(i,k,j)=rtmp(ndx(i))
113   12 continue
114      do 13 i=1,nwmloc
115      rtmp(i)=vw(i,k,j)
116   13 continue
117      do 14 i=1,nwmloc
118      vw(i,k,j)=rtmp(ndx(i))
119   14 continue
120      if(iguide.gt.0) then
121      do 113 i=1,nwmloc
122      rtmp(i)=gw(i,k,j)
123  113 continue
124      do 114 i=1,nwmloc
125      gw(i,k,j)=rtmp(ndx(i))
126  114 continue
127      endif
128   10 continue
129    5 continue
130      do 18 k=1,miw2
131      do 19 i=1,nwmloc
132      itmp(i)=iwl(i,k)
133   19 continue
134      do 20 i=1,nwmloc
135      iwl(i,k)=itmp(ndx(i))
136   20 continue
137   18 continue
138      endif
139c
140      if(nwmloc.gt.0) then
141      do 21 iwm=1,nwmloc
142      if(iwl(iwm,lwnod).eq.me) then
143      nwstay=iwm
144      else
145c
146c     check if moving atoms go to neighboring processor
147c
148      do 222 k=1,27
149      if(iwl(iwm,lwnod).eq.neighb(k,1)) goto 223
150  222 continue
151      write(string,'(a,i4,a,i4,a,i4,9f6.2)')
152     +  'argos_space_travel: solvent molecule ',
153     + iwl(iwm,lwgmn),' moving to non-neighbor ',iwl(iwm,lwnod),
154     + ' from ',me,((xw(iwm,i,j),i=1,3),j=1,3)
155      call md_abort(string,me)
156  223 continue
157c
158      endif
159c
160c     testcode
161c
162      if(iand(idebug,8).eq.8) then
163      if(iwl(iwm,lwnod).ne.me) write(lfndbg,'(a,3i5)')
164     +  'Travel w fnd ',me,iwl(iwm,lwnod),iwl(iwm,lwgmn)
165      endif
166c
167c     end test code
168c
169   21 continue
170      endif
171c
172c     order the solute atoms
173c
174c     isl(isa,lsbox) : box
175c     isl(isa,lsnod) : node
176c     isl(isa,lssgm) : segment
177c
178      nsstay=0
179      if(nsaloc.gt.0) then
180      do 22 i=1,nsaloc
181      ndx(i)=i
182   22 continue
183      endif
184c
185      if(nsaloc.gt.1) then
186      lsaloc=nsaloc/2+1
187      irs=nsaloc
188   23 continue
189      if(lsaloc.gt.1) then
190      lsaloc=lsaloc-1
191      itemps=ndx(lsaloc)
192      else
193      itemps=ndx(irs)
194      ndx(irs)=ndx(1)
195      irs=irs-1
196      if(irs.eq.1) then
197      ndx(1)=itemps
198      goto 24
199      endif
200      endif
201      isaloc=lsaloc
202      jsaloc=lsaloc+lsaloc
203   25 continue
204      if(jsaloc.le.irs) then
205      if(jsaloc.lt.irs) then
206      if((isl(ndx(jsaloc),lsnod).eq.isl(ndx(jsaloc+1),lsnod).and.
207     + (isl(ndx(jsaloc),lsbox).lt.isl(ndx(jsaloc+1),lsbox).or.
208     + (isl(ndx(jsaloc),lsbox).eq.isl(ndx(jsaloc+1),lsbox).and.
209     + isl(ndx(jsaloc),lssgm).le.isl(ndx(jsaloc+1),lssgm)))).or.
210     + ((isl(ndx(jsaloc),lsnod).eq.me.or.
211     + (isl(ndx(jsaloc),lsnod).ne.me.and.
212     + isl(ndx(jsaloc),lsnod).le.isl(ndx(jsaloc+1),lsnod))).and.
213     + isl(ndx(jsaloc+1),lsnod).ne.me)) jsaloc=jsaloc+1
214      endif
215      if((isl(itemps,lsnod).eq.isl(ndx(jsaloc),lsnod).and.
216     + (isl(itemps,lsbox).lt.isl(ndx(jsaloc),lsbox).or.
217     + (isl(itemps,lsbox).eq.isl(ndx(jsaloc),lsbox).and.
218     + isl(itemps,lssgm).le.isl(ndx(jsaloc),lssgm)))).or.
219     + ((isl(itemps,lsnod).eq.me.or. (isl(itemps,lsnod).ne.me.and.
220     + isl(itemps,lsnod).le.isl(ndx(jsaloc),lsnod))).and.
221     + isl(ndx(jsaloc),lsnod).ne.me)) then
222      ndx(isaloc)=ndx(jsaloc)
223      isaloc=jsaloc
224      jsaloc=jsaloc+jsaloc
225      else
226      jsaloc=irs+1
227      endif
228      goto 25
229      endif
230      ndx(isaloc)=itemps
231      goto 23
232   24 continue
233c
234      do 26 k=1,3
235      do 27 i=1,nsaloc
236      rtmp(i)=xs(i,k)
237   27 continue
238      do 28 i=1,nsaloc
239      xs(i,k)=rtmp(ndx(i))
240   28 continue
241      do 29 i=1,nsaloc
242      rtmp(i)=vs(i,k)
243   29 continue
244      do 30 i=1,nsaloc
245      vs(i,k)=rtmp(ndx(i))
246   30 continue
247      if(iguide.gt.0) then
248      do 2129 i=1,nsaloc
249      rtmp(i)=gs(i,k)
250 2129 continue
251      do 2130 i=1,nsaloc
252      gs(i,k)=rtmp(ndx(i))
253 2130 continue
254      endif
255   26 continue
256      do 40 k=1,mis2
257      do 41 i=1,nsaloc
258      itmp(i)=isl(i,k)
259   41 continue
260      do 42 i=1,nsaloc
261      isl(i,k)=itmp(ndx(i))
262   42 continue
263   40 continue
264      endif
265c
266      if(nsa.gt.0) then
267      do 43 isa=1,nsaloc
268      if(isl(isa,lsnod).eq.me) then
269      nsstay=isa
270      else
271c
272c     check if moving atoms go to neighboring processor
273c
274      do 444 k=1,27
275      if(isl(isa,lsnod).eq.neighb(k,1)) goto 445
276  444 continue
277      write(string,'(a,i4,a,i4,a,i4,3f6.2)')
278     +  'argos_space_travel: solute segment ',
279     + isl(isa,lssgm),' moving to non-neighbor ',isl(isa,lsnod),
280     + ' from ',me,(xs(isa,i),i=1,3)
281      call md_abort(string,me)
282  445 continue
283c
284      endif
285   43 continue
286      endif
287c
288c     make packages ready for shipment
289c
290c     loop over all neighboring nodes
291c
292      call ga_distribution(ga_iw,me,iliw,ihiw,jliw,jhiw)
293      call ga_distribution(ga_w,me,ilw,ihw,jlw,jhw)
294      call ga_distribution(ga_is,me,ilis,ihis,jlis,jhis)
295      call ga_distribution(ga_s,me,ils,ihs,jls,jhs)
296c
297      indexw=0
298      indexs=0
299      nwgosm=0
300c
301      do 70 i=1,27
302      jnode=neighb(i,1)
303      if(jnode.ge.0.and.jnode.ne.me) then
304c
305c     for the solvent
306c
307      iwfr=0
308      iwto=0
309      do 71 iwm=nwstay+1,nwmloc
310      if(iwl(iwm,lwnod).eq.jnode) then
311      if(iwfr.eq.0) iwfr=iwm
312      iwto=iwm
313c
314c     testcode
315c
316      if(iand(idebug,8).eq.8) then
317      if(iwl(iwm,lwnod).ne.me) write(lfndbg,'(a,3i5)')
318     +  'Travel w snd ',me,iwl(iwm,lwnod),iwl(iwm,lwgmn)
319      endif
320c
321c     end test code
322c
323      endif
324   71 continue
325c
326c     if molecules need to travel copy coordinates etc into global array
327c
328      nwgo=iwto-iwfr+1
329      if(iwfr.eq.0) nwgo=0
330      ipl(1,1)=0
331      ipl(1,2)=0
332c
333      if(nwgo.gt.0) then
334      nwgosm=nwgosm+nwgo
335      il=iliw+indexw
336      ih=il+nwgo-1
337      if(npackw.eq.0) then
338      call ga_put(ga_iw,il,ih,jliw,jhiw,iwl(iwfr,1),mwm)
339      else
340      call argos_space_packw(ih-il+1,iwl(iwfr,1),iwlp(iwfr,1))
341      call ga_put(ga_iw,il,ih,jliw,jliw+npackw-1,iwlp(iwfr,1),mwm)
342      endif
343      il=ilw+indexw
344      ih=il+nwgo-1
345      call ga_put(ga_w,il,ih,jlw,jlw+3*mwa-1,xw(iwfr,1,1),mwm)
346      call ga_put(ga_w,il,ih,jlw+3*mwa,jlw+6*mwa-1,vw(iwfr,1,1),mwm)
347      call ga_put(ga_w,il,ih,jlw+6*mwa,jlw+6*mwa+2,xwcr(iwfr,1),mwm)
348      if(iguide.gt.0) then
349      call ga_put(ga_w,il,ih,jlw+6*mwa+3,jlw+9*mwa+2,gw(iwfr,1,1),mwm)
350      endif
351      ipl(1,1)=indexw+1
352      ipl(1,2)=indexw+nwgo
353      indexw=indexw+nwgo
354      endif
355c
356c     for the solute
357c
358      isfr=0
359      isto=0
360      do 72 isa=nsstay+1,nsaloc
361      if(isl(isa,lsnod).eq.jnode) then
362      if(isfr.eq.0) isfr=isa
363      isto=isa
364      endif
365   72 continue
366      nsgo=isto-isfr+1
367      if(isfr.eq.0) nsgo=0
368      ipl(1,3)=0
369      ipl(1,4)=0
370      if(nsgo.gt.0) then
371      il=ilis+indexs
372      ih=il+nsgo-1
373      if(npack.eq.0) then
374      call ga_put(ga_is,il,ih,jlis,jhis,isl(isfr,1),msa)
375      else
376      call argos_space_pack(ih-il+1,isl(isfr,1),islp(isfr,1))
377      call ga_put(ga_is,il,ih,jlis,jlis+npack-1,islp(isfr,1),msa)
378      endif
379      call ga_put(ga_s,il,ih,jls,jls+2,xs(isfr,1),msa)
380      call ga_put(ga_s,il,ih,jls+3,jls+5,vs(isfr,1),msa)
381      if(iguide.gt.0) then
382      call ga_put(ga_s,il,ih,jls+6,jls+8,gs(isfr,1),msa)
383      endif
384      ipl(1,3)=indexs+1
385      ipl(1,4)=indexs+nsgo
386      indexs=indexs+nsgo
387      endif
388c
389c     inform other node of number of molecules to get
390c
391      if(ipl(1,1).gt.0.or.ipl(1,3).gt.0) then
392      call ga_distribution(ga_ip,jnode,ilp,ihp,jlp,jhp)
393      ilp=ilp+2+i
394      call ga_put(ga_ip,ilp,ilp,jlp,jhp,ipl,mbox)
395      endif
396      endif
397      lrec(i)=.false.
398   70 continue
399c
400      call ga_sync()
401c
402c     receive molecules from other nodes
403c
404      nwgtsm=0
405c
406      call ga_distribution(ga_ip,me,ilp,ihp,jlp,jhp)
407      call ga_get(ga_ip,ilp,ilp+30,jlp,jhp,ipl,mbox)
408c
409      do 74 i=1,27
410      jnode=neighb(i,2)
411      if(jnode.ge.0.and.jnode.ne.me.and..not.lrec(i)) then
412c
413      iwfr=ipl(3+i,1)
414      iwto=ipl(3+i,2)
415      isfr=ipl(3+i,3)
416      isto=ipl(3+i,4)
417c
418      nwnew=iwto-iwfr+1
419      nsnew=isto-isfr+1
420c
421      if(iwfr.eq.0) nwnew=0
422      if(isfr.eq.0) nsnew=0
423c
424      if(nwstay+nwnew.gt.mwm) then
425      write(string,'(a,i7,a,i7)')
426     + 'Travel: mwm needs increase with ',nwnew,' to ',nwstay+nwnew
427      call md_abort(string,me)
428      endif
429      if(nsstay+nsnew.gt.msa) then
430      write(string,'(a,i7,a,i7)')
431     + 'Travel: msa needs increase with ',nsnew,' to ',nsstay+nsnew
432      call md_abort(string,me)
433      endif
434c
435      lrec(i)=.true.
436c
437      if(iwfr.gt.0) then
438      nwgtsm=nwgtsm+nwnew
439      iwto=ipl(3+i,2)
440      call ga_distribution(ga_iw,jnode,iliw,ihiw,jliw,jhiw)
441      call ga_distribution(ga_w,jnode,ilw,ihw,jlw,jhw)
442c
443c     get data for additional molecules
444c
445      il=iliw+iwfr-1
446      ih=iliw+iwto-1
447      if(npackw.eq.0) then
448      call ga_get(ga_iw,il,ih,jliw,jhiw,iwl(nwstay+1,1),mwm)
449      else
450      call ga_get(ga_iw,il,ih,jliw,jliw+npackw-1,iwlp(nwstay+1,1),mwm)
451      call argos_space_unpackw(ih-il+1,iwl(nwstay+1,1),iwlp(nwstay+1,1))
452      endif
453      call ga_get(ga_w,il,ih,jlw,jlw+3*mwa-1,xw(nwstay+1,1,1),mwm)
454      call ga_get(ga_w,il,ih,jlw+3*mwa,jlw+6*mwa-1,vw(nwstay+1,1,1),mwm)
455      call ga_get(ga_w,il,ih,jlw+6*mwa,jlw+6*mwa+2,xwcr(nwstay+1,1),mwm)
456      if(iguide.gt.0) then
457      call ga_get(ga_w,il,ih,jlw+6*mwa+3,jlw+9*mwa+2,
458     + gw(nwstay+1,1,1),mwm)
459      endif
460c
461c     testcode
462c
463      if(iand(idebug,8).eq.8) then
464      write(lfndbg,'(a,3i5)')
465     +  ('Travel w rcv ',me,jnode,iwl(nwstay+k,lwgmn),k=1,nwnew)
466      endif
467c
468c     end test code
469c
470c
471c     update number of local solvent molecules
472c
473      nwstay=nwstay+nwnew
474c
475      endif
476c
477c     for the solute
478c
479      if(isfr.gt.0) then
480      call ga_distribution(ga_is,jnode,ilis,ihis,jlis,jhis)
481      call ga_distribution(ga_s,jnode,ils,ihs,jls,jhs)
482      il=ilis+isfr-1
483      ih=ilis+isto-1
484      jl=jlis
485      jh=jhis
486      if(npack.eq.0) then
487      call ga_get(ga_is,il,ih,jlis,jhis,isl(nsstay+1,1),msa)
488      else
489      call ga_get(ga_is,il,ih,jlis,jlis+npack-1,islp(nsstay+1,1),msa)
490      call argos_space_unpack(ih-il+1,isl(nsstay+1,1),islp(nsstay+1,1))
491      endif
492      call ga_get(ga_s,il,ih,jls,jls+2,xs(nsstay+1,1),msa)
493      call ga_get(ga_s,il,ih,jls+3,jls+5,vs(nsstay+1,1),msa)
494      if(iguide.gt.0) then
495      call ga_get(ga_s,il,ih,jls+6,jls+8,gs(nsstay+1,1),msa)
496      endif
497c
498      nsstay=nsstay+nsnew
499      endif
500c
501      endif
502c
503c     reset the pointers to zero
504c
505      ipl(3+i,1)=0
506      ipl(3+i,2)=0
507      ipl(3+i,3)=0
508      ipl(3+i,4)=0
509c
510   74 continue
511c
512c     reset ipl in global array
513c
514      call ga_put(ga_ip,ilp,ilp+30,jlp,jhp,ipl,mbox)
515c
516c     order the solvent molecules according to subbox and
517c     store indices into ip
518c
519c     ip(1,1)    : number of boxes on this node
520c     ip(1,2)    : number of solvent molecules on this node
521c     ip(2,2)    : number of solute atoms on this node
522c
523c     ip(3+i,1)  : index for solvents to be moved to the i-th neighbor
524c
525c     ip(30+i,1) : number of i-th box on this node
526c     ip(30+i,2) : index to first solvent in i-th box
527c     ip(30+i,3) : index to lasst solvent in i-th box
528c
529      if(nwstay.gt.0.and.(nwgosm.gt.0.or.nwgtsm.gt.0)) then
530      do 81 i=1,nwstay
531      ndx(i)=i
532   81 continue
533      if(nwstay.gt.1) then
534      lwstay=nwstay/2+1
535      irw=nwstay
536   82 continue
537      if(lwstay.gt.1) then
538      lwstay=lwstay-1
539      itemps=ndx(lwstay)
540      else
541      itemps=ndx(irw)
542      ndx(irw)=ndx(1)
543      irw=irw-1
544      if(irw.eq.1) then
545      ndx(1)=itemps
546      goto 83
547      endif
548      endif
549      iwstay=lwstay
550      jwstay=lwstay+lwstay
551   84 continue
552      if(jwstay.le.irw) then
553      if(jwstay.lt.irw) then
554      if(iwl(ndx(jwstay),lwbox).le.iwl(ndx(jwstay+1),lwbox))
555     + jwstay=jwstay+1
556      endif
557      if(iwl(itemps,lwbox).le.iwl(ndx(jwstay),lwbox)) then
558      ndx(iwstay)=ndx(jwstay)
559      iwstay=jwstay
560      jwstay=jwstay+jwstay
561      else
562      jwstay=irw+1
563      endif
564      goto 84
565      endif
566      ndx(iwstay)=itemps
567      goto 82
568   83 continue
569c
570      do 85 k=1,3
571      do 88 i=1,nwstay
572      rtmp(i)=xwcr(i,k)
573   88 continue
574      do 89 i=1,nwstay
575      xwcr(i,k)=rtmp(ndx(i))
576   89 continue
577      do 90 j=1,nwa
578      do 91 i=1,nwstay
579      rtmp(i)=xw(i,k,j)
580   91 continue
581      do 92 i=1,nwstay
582      xw(i,k,j)=rtmp(ndx(i))
583   92 continue
584      do 93 i=1,nwstay
585      rtmp(i)=vw(i,k,j)
586   93 continue
587      do 94 i=1,nwstay
588      vw(i,k,j)=rtmp(ndx(i))
589   94 continue
590      if(iguide.gt.0) then
591      do 193 i=1,nwstay
592      rtmp(i)=gw(i,k,j)
593  193 continue
594      do 194 i=1,nwstay
595      gw(i,k,j)=rtmp(ndx(i))
596  194 continue
597      endif
598   90 continue
599   85 continue
600      do 98 k=1,miw2
601      do 99 i=1,nwstay
602      itmp(i)=iwl(i,k)
603   99 continue
604      do 100 i=1,nwstay
605      iwl(i,k)=itmp(ndx(i))
606  100 continue
607   98 continue
608c
609      endif
610      endif
611c
612c     order the solute according to segment
613c
614      if(nsstay.gt.0) then
615      do 122 i=1,nsstay
616      ndx(i)=i
617  122 continue
618      if(nsstay.gt.1) then
619      lsstay=nsstay/2+1
620      irs=nsstay
621  123 continue
622      if(lsstay.gt.1) then
623      lsstay=lsstay-1
624      itemps=ndx(lsstay)
625      else
626      itemps=ndx(irs)
627      ndx(irs)=ndx(1)
628      irs=irs-1
629      if(irs.eq.1) then
630      ndx(1)=itemps
631      goto 124
632      endif
633      endif
634      isstay=lsstay
635      jsstay=lsstay+lsstay
636  125 continue
637      if(jsstay.le.irs) then
638      if(jsstay.lt.irs) then
639      if(isl(ndx(jsstay),lssgm).le.isl(ndx(jsstay+1),lssgm))
640     + jsstay=jsstay+1
641      endif
642      if(isl(itemps,lssgm).le.isl(ndx(jsstay),lssgm)) then
643      ndx(isstay)=ndx(jsstay)
644      isstay=jsstay
645      jsstay=jsstay+jsstay
646      else
647      jsstay=irs+1
648      endif
649      goto 125
650      endif
651      ndx(isstay)=itemps
652      goto 123
653  124 continue
654      endif
655c
656c     for each segment : 1. determine box number
657c                        2. assign box number to each atom
658c                        3. when box not owned by node:
659c                           a. assign box number
660c                           b. assign correct node number
661c
662      goto 666
663      icsgm=isl(ndx(1),lssgm)
664      ifsgm=1
665      ilsgm=1
666      do 126 isa=2,nsstay+1
667c
668c     if isa is first atom of a new segment or very last atom
669c
670
671      if(isa.le.nsstay) then
672      lend=isl(ndx(isa),lssgm).ne.icsgm
673      else
674      lend=.true.
675      endif
676      if(lend) then
677      if(isa.gt.nsstay) ilsgm=nsstay
678      if(ifsgm.gt.0.and.ilsgm.ge.ifsgm) then
679      xscx=zero
680      xscy=zero
681      xscz=zero
682      do 127 jsa=ifsgm,ilsgm
683      xscx=xscx+xs(ndx(jsa),1)
684      xscy=xscy+xs(ndx(jsa),2)
685      xscz=xscz+xs(ndx(jsa),3)
686  127 continue
687      factor=one/dble(ilsgm-ifsgm+1)
688      xscx=factor*xscx
689      xscy=factor*xscy
690      xscz=factor*xscz
691      if(npbtyp.ne.0) then
692      if(abs(xscx).gt.boxh(1)) then
693      xscx=xscx-nint(xscx*boxi(1))*box(1)
694      nfold=1
695      endif
696      if(abs(xscy).gt.boxh(2)) then
697      xscy=xscy-nint(xscy*boxi(2))*box(2)
698      nfold=1
699      endif
700      if(abs(xscz).gt.boxh(3)) then
701      xscz=xscz-nint(xscz*boxi(3))*box(3)
702      nfold=1
703      endif
704      endif
705c
706c     determine the box number
707c
708      ibx=0
709      iby=0
710      ibz=0
711      do 128 i=1,nbx-1
712      if(xscx+boxh(1).gt.boxsiz(i,1)) ibx=i
713  128 continue
714      do 129 i=1,nby-1
715      if(xscy+boxh(2).gt.boxsiz(i,2)) iby=i
716  129 continue
717      do 1130 i=1,nbz-1
718      if(xscz+boxh(3).gt.boxsiz(i,3)) ibz=i
719 1130 continue
720      if(npbtyp.gt.0) then
721      if(ibx.ge.nbx) ibx=ibx-nbx
722      if(iby.ge.nby) iby=iby-nby
723      if(ibx.lt.0) ibx=ibx+nbx
724      if(iby.lt.0) iby=iby+nby
725      if(npbtyp.eq.1) then
726      if(ibz.ge.nbz) ibz=ibz-nbz
727      if(ibz.lt.0) ibz=ibz+nbz
728      else
729      if(ibz.ge.nbz) ibz=nbz-1
730      if(ibz.lt.0) ibz=0
731      endif
732      else
733      if(ibx.ge.nbx) ibx=nbx-1
734      if(iby.ge.nby) iby=nby-1
735      if(ibz.ge.nbz) ibz=nbz-1
736      if(ibx.lt.0) ibx=0
737      if(iby.lt.0) iby=0
738      if(ibz.lt.0) ibz=0
739      endif
740      ipx=ibownr(ibx+1,1)
741      ipy=ibownr(iby+1,2)
742      ipz=ibownr(ibz+1,3)
743      isbox=(ibz*nby+iby)*nbx+ibx
744      isnod=(ipz*npy+ipy)*npx+ipx
745c
746c     assign box and node numbers
747c
748      do 1131 jsa=ifsgm,ilsgm
749      isl(ndx(jsa),lsbox)=isbox
750      isl(ndx(jsa),lsnod)=isnod
751 1131 continue
752c
753      endif
754      if(isa.le.nsstay) icsgm=isl(ndx(isa),lssgm)
755      ifsgm=isa
756      else
757      ilsgm=isa
758      endif
759  126 continue
760  666 continue
761c
762c     order solute according to box, segment, charge group, atom number
763c
764      if(nsstay.gt.1) then
765      lsstay=nsstay/2+1
766      irs=nsstay
767  132 continue
768      if(lsstay.gt.1) then
769      lsstay=lsstay-1
770      itemps=ndx(lsstay)
771      else
772      itemps=ndx(irs)
773      ndx(irs)=ndx(1)
774      irs=irs-1
775      if(irs.eq.1) then
776      ndx(1)=itemps
777      goto 133
778      endif
779      endif
780      isstay=lsstay
781      jsstay=lsstay+lsstay
782  134 continue
783      if(jsstay.le.irs) then
784      if(jsstay.lt.irs) then
785      if(isl(ndx(jsstay),lsbox).lt.isl(ndx(jsstay+1),lsbox).or.
786     + (isl(ndx(jsstay),lsbox).eq.isl(ndx(jsstay+1),lsbox).and.
787     + (isl(ndx(jsstay),lssgm).lt.isl(ndx(jsstay+1),lssgm).or.
788     + (isl(ndx(jsstay),lssgm).eq.isl(ndx(jsstay+1),lssgm).and.
789     + (isl(ndx(jsstay),lsgrp).lt.isl(ndx(jsstay+1),lsgrp).or.
790     + (isl(ndx(jsstay),lsgrp).eq.isl(ndx(jsstay+1),lsgrp).and.
791     + isl(ndx(jsstay),lsgan).le.isl(ndx(jsstay+1),lsgan)))))))
792     + jsstay=jsstay+1
793      endif
794      if(isl(itemps,lsbox).lt.isl(ndx(jsstay),lsbox).or.
795     + (isl(itemps,lsbox).eq.isl(ndx(jsstay),lsbox).and.
796     + (isl(itemps,lssgm).lt.isl(ndx(jsstay),lssgm).or.
797     + (isl(itemps,lssgm).eq.isl(ndx(jsstay),lssgm).and.
798     + (isl(itemps,lsgrp).lt.isl(ndx(jsstay),lsgrp).or.
799     + (isl(itemps,lsgrp).eq.isl(ndx(jsstay),lsgrp).and.
800     + isl(itemps,lsgan).le.isl(ndx(jsstay),lsgan))))))) then
801      ndx(isstay)=ndx(jsstay)
802      isstay=jsstay
803      jsstay=jsstay+jsstay
804      else
805      jsstay=irs+1
806      endif
807      goto 134
808      endif
809      ndx(isstay)=itemps
810      goto 132
811  133 continue
812      endif
813c
814      do 135 k=1,3
815      do 136 i=1,nsstay
816      rtmp(i)=xs(i,k)
817  136 continue
818      do 137 i=1,nsstay
819      xs(i,k)=rtmp(ndx(i))
820  137 continue
821      do 138 i=1,nsstay
822      rtmp(i)=vs(i,k)
823  138 continue
824      do 139 i=1,nsstay
825      vs(i,k)=rtmp(ndx(i))
826  139 continue
827      if(iguide.gt.0) then
828      do 1138 i=1,nsstay
829      rtmp(i)=gs(i,k)
830 1138 continue
831      do 1139 i=1,nsstay
832      gs(i,k)=rtmp(ndx(i))
833 1139 continue
834      endif
835  135 continue
836      do 149 k=1,mis2
837      do 150 i=1,nsstay
838      itmp(i)=isl(i,k)
839  150 continue
840      do 151 i=1,nsstay
841      isl(i,k)=itmp(ndx(i))
842  151 continue
843  149 continue
844c
845      endif
846c
847      do 200 i=1,ipl(1,1)
848      ipl(30+i,2)=0
849      ipl(30+i,3)=0
850      ipl(30+i,4)=0
851      ipl(30+i,5)=0
852  200 continue
853c
854      do 201 i=1,ipl(1,1)
855      nrbox=ipl(30+i,1)
856      if(nwstay.gt.0) then
857      do 202 iwm=1,nwstay
858      if(iwl(iwm,lwbox).eq.nrbox) then
859      if(ipl(30+i,2).eq.0) ipl(30+i,2)=iwm
860      ipl(30+i,3)=iwm
861      endif
862  202 continue
863      endif
864      if(nsstay.gt.0) then
865      do 203 isa=1,nsstay
866      if(isl(isa,lsbox).eq.nrbox) then
867      if(ipl(30+i,4).eq.0) ipl(30+i,4)=isa
868      ipl(30+i,5)=isa
869      endif
870  203 continue
871      endif
872  201 continue
873c
874      nwmloc=nwstay
875      ipl(1,2)=nwmloc
876      nsaloc=nsstay
877      ipl(2,2)=nsaloc
878c
879      call ga_igop(msp_23,nfold,1,'+')
880      lpbc9=nfold.gt.0
881c
882      return
883      end
884c $Id$
885