1      logical function cf_hopping(lpbc,lpbcs,stimei,
2     + isl,issgm,isgan,isq3,ishop,xs,nsaloc)
3c
4c $Id$
5c
6      implicit none
7c
8#include "cf_common.fh"
9#include "global.fh"
10#include "mafdecls.fh"
11#include "msgids.fh"
12c
13      logical cf_hop
14      external cf_hop
15c
16      logical lpbc,lpbcs
17      integer isl(msa,mis2)
18      integer issgm(msa),nsaloc,isgan(msa),ishop(msa),isq3(msa)
19      real*8 xs(msa,3)
20      real*8 stimei
21c
22      integer i_itmp,l_itmp,i_dtmp,l_dtmp
23c
24      stime=stimei
25      cf_hopping=.false.
26c
27      if(lpair.and.nhop.gt.0) then
28      lhop=.true.
29c
30      if(.not.ma_push_get(mt_int,np,'itmp',l_itmp,i_itmp))
31     + call md_abort('Failed to allocate itmp',me)
32c
33      call cf_hoplist(issgm,nsaloc,int_mb(i_itmp),int_mb(i_lda),
34     + dbl_mb(i_rda))
35c
36      if(.not.ma_pop_stack(l_itmp))
37     + call md_abort('Failed to deallocate itmp',me)
38c
39      cf_hopping=.true.
40      elseif(lhop) then
41c
42c     evaluate the E12
43c
44      call cf_uhop(lpbc,lpbcs,int_mb(i_lseq),isl,issgm,isgan,isq3,
45     + xs,nsaloc,
46     + dbl_mb(i_chg),int_mb(i_lda),dbl_mb(i_rda),dbl_mb(i_uda))
47c
48      if(.not.ma_push_get(mt_dbl,nldat,'dtmp',l_dtmp,i_dtmp))
49     + call md_abort('Failed to allocate dtmp',me)
50      if(.not.ma_push_get(mt_int,nldat,'itmp',l_itmp,i_itmp))
51     + call md_abort('Failed to allocate itmp',me)
52c
53      cf_hopping=cf_hop(int_mb(i_lseq),issgm,
54     + int_mb(i_lda),dbl_mb(i_rda),dbl_mb(i_uda),
55     + int_mb(i_itmp),dbl_mb(i_dtmp),isgan,ishop,xs,nsaloc,
56     + int_mb(i_lsthop),dbl_mb(i_timhop))
57c
58      if(.not.ma_pop_stack(l_itmp))
59     + call md_abort('Failed to deallocate itmp',me)
60      if(.not.ma_pop_stack(l_dtmp))
61     + call md_abort('Failed to deallocate dtmp',me)
62c
63      endif
64c
65      return
66      end
67      subroutine cf_hoplist(issgm,nsaloc,nhopl,lda,rda)
68c
69      implicit none
70c
71#include "cf_common.fh"
72#include "global.fh"
73#include "mafdecls.fh"
74#include "msgids.fh"
75#include "util.fh"
76c
77      integer issgm(msa),nsaloc,nhopl(np)
78      integer lda(16,*)
79      real*8 rda(11,*)
80c
81      integer i,j,ioff
82c
83      do 1 i=1,np
84      nhopl(i)=0
85    1 continue
86      nhopl(me+1)=nlda
87c      write(*,'(i5,a,i5)') me,' nlda= ',nlda
88c
89      call ga_igop(mcf_76,nhopl,np,'+')
90c
91      do 2 i=2,np
92      nhopl(i)=nhopl(i)+nhopl(i-1)
93    2 continue
94c      write(*,'(i5,a,10i5)') me,' nhopl ',(nhopl(i),i=1,np)
95c
96c
97      ioff=0
98      if(me.gt.0) then
99      ioff=nhopl(me)
100      do 3 i=1,nlda
101      do 4 j=1,16
102      lda(j,i+ioff)=lda(j,i)
103    4 continue
104      do 5 j=1,8
105      rda(j,i+ioff)=rda(j,i)
106    5 continue
107    3 continue
108      do 6 i=1,nhopl(me)
109      do 7 j=1,16
110      lda(j,i)=0
111    7 continue
112      do 8 j=1,8
113      rda(j,i)=0.0d0
114    8 continue
115    6 continue
116      endif
117      if(me.lt.np-1) then
118      do 9 i=nhopl(me+1)+1,nhopl(np)
119      do 10 j=1,16
120      lda(j,i)=0
121   10 continue
122      do 11 j=1,8
123      rda(j,i)=0.0d0
124   11 continue
125    9 continue
126      endif
127c
128      nldat=nhopl(np)
129      call ga_igop(mcf_77,lda,16*nldat,'+')
130c
131c      do 322 i=1,nldat
132c      write(*,'(18i5)') me,i,(lda(j,i),j=1,16)
133c  322 continue
134c
135      do 201 i=1,nldat
136      do 202 j=1,nsaloc
137      if(lda(3,i).eq.issgm(j)) then
138      if(lda(12,i).eq.0) lda(12,i)=j
139      lda(13,i)=j
140      endif
141      if(lda(8,i).eq.issgm(j)) then
142      if(lda(15,i).eq.0) lda(15,i)=j
143      lda(16,i)=j
144      endif
145  202 continue
146  201 continue
147c
148      if(me.gt.0) then
149      do 221 j=1,nldat
150      do 222 i=1,11
151      lda(i,j)=0
152  222 continue
153      lda(14,j)=0
154  221 continue
155      endif
156c
157      call ga_igop(mcf_77,lda,16*nldat,'+')
158      call ga_dgop(mcf_78,rda,11*nldat,'+')
159c
160c      do 2222 i=1,nlda
161c      write(*,'(18i5)') me,i,(lda(j,i),j=1,16)
162c 2222 continue
163c
164      nlda=0
165      do 12 i=1,nldat
166      do 13 j=1,nlda
167      if(lda(1,i).eq.lda(1,j).and.lda(6,i).eq.lda(6,j)) goto 12
168   13 continue
169      nlda=nlda+1
170      do 14 j=1,16
171      lda(j,nlda)=lda(j,i)
172   14 continue
173      do 15 j=1,8
174      rda(j,nlda)=rda(j,i)
175   15 continue
176   12 continue
177      nldat=nlda
178c
179c
180      if(me.eq.0) then
181      if(util_print('qhop',print_high)) then
182      write(lfnhop,110) nldat
183  110 format(/,'Number of donor-acceptor pairs is ',i5,/,
184     + '---------- Donor --------- --------- Acceptor ----',
185     + ' --- Donor ---  ---Acceptor---',/,
186     + ' glob  loc  sgm prot  off glob  loc  sgm prot  off',
187     + '  own from   to  own from   to',/)
188      do 212 i=1,nldat
189      write(lfnhop,111) (lda(j,i),j=1,16),(rda(j,i),j=1,3)
190  111 format(16i5,3f12.6)
191  212 continue
192      endif
193      endif
194c
195      return
196      end
197      logical function cf_hop(lseq,issgm,lda,rda,uda,ndx,pda,
198     + isgan,ishop,xs,nsaloc,lsthop,timhop)
199c
200      implicit none
201c
202#include "cf_common.fh"
203#include "global.fh"
204#include "mafdecls.fh"
205#include "msgids.fh"
206#include "util.fh"
207#include "bitops_decls.fh"
208#include "bitops_funcs.fh"
209c
210      integer lseq(mseq),issgm(msa)
211      integer lda(16,*),nsaloc,isgan(msa),ishop(msa)
212      real*8 rda(11,*),uda(4,*),pda(*),xs(msa,3),x(3)
213      integer ndx(*)
214      integer lsthop(2,*)
215      real*8 timhop(*)
216c
217      integer i,j,k,ndxrel
218      real*8 dprob,drand,e12,e120
219c
220      call ga_dgop(mcf_79,uda,4*nldat,'+')
221c
222      if(me.eq.0) then
223      if(util_print('qhop',print_high)) then
224      write(lfnhop,110) stime
225  110 format(/,' Time ',f12.6,/,
226     + ' QHOP  ----donor---   --acceptor-- ',
227     + '     pre-hop energies        post-hop energies  ',/,
228     + ' numb atom  res prot atom  res prot',
229     + '    donor      acceptor     donor      acceptor ',
230     + '     e12        e120        distance     angle   probability',/)
231      endif
232      do 1 i=1,nldat
233c      write(*,'(i5,8f12.6)') i,(rda(j,i),j=1,8)
234      call qhop_prob(lda(1,i),lda(3,i),lseq(lda(3,i)),
235     + lda(6,i),lda(8,i),lseq(lda(8,i)),uda(1,i),rda(4,i),pda(i),e120)
236c      write(*,'(i5,8f12.6)') i,(rda(j,i),j=1,8)
237c      pda(i)=0.0d0
238c
239      if(util_print('qhop',print_high)) then
240      e12=uda(3,i)+uda(4,i)-uda(1,i)-uda(2,i)
241      write(lfnhop,111) i,lda(1,i),lda(3,i),lseq(lda(3,i)),
242     + lda(6,i),lda(8,i),lseq(lda(8,i)),
243     + (uda(j,i),j=1,4),e12,e120*4.184,rda(4,i),rda(5,i),pda(i)
244  111 format(7i5,8f12.6,f12.9)
245      endif
246c
247    1 continue
248c
249      do 2 i=1,nldat
250      ndx(i)=i
251    2 continue
252c
253      do 3 i=1,nldat-1
254      do 4 j=i+1,nldat
255      if(pda(ndx(i)).lt.pda(ndx(j))) then
256      k=ndx(i)
257      ndx(i)=ndx(j)
258      ndx(j)=k
259      endif
260    4 continue
261    3 continue
262c
263      do 5 i=1,nldat
264      drand=util_random(0)
265      dprob=pda(ndx(i))
266c      write(*,'(a,2f12.6)') 'RAND/PROB ',drand,dprob
267      if(drand.lt.dprob) then
268      write(lfnhop,'(/,a,i5,a,i5,a,f12.6/)')
269     + ' PROTON HOP FROM ',lda(3,ndx(i)),' TO ',lda(8,ndx(i)),
270     + ' AT TIME ',stime
271      pda(ndx(i))=1.0d0
272      do 6 j=i+1,nldat
273      if(lda(3,ndx(i)).eq.lda(3,ndx(j))) pda(ndx(j))=0.0d0
274      if(lda(8,ndx(i)).eq.lda(8,ndx(j))) pda(ndx(j))=0.0d0
275      if(lda(3,ndx(i)).eq.lda(8,ndx(j))) pda(ndx(j))=0.0d0
276      if(lda(8,ndx(i)).eq.lda(3,ndx(j))) pda(ndx(j))=0.0d0
277    6 continue
278      else
279      pda(ndx(i))=0.0d0
280      endif
281    5 continue
282c
283      endif
284c
285      call ga_brdcst(mcf_80,pda,nldat*ma_sizeof(mt_dbl,1,mt_byte),0)
286c
287      cf_hop=.false.
288      do 7 i=1,nldat
289      if(pda(i).gt.0.5d0) then
290      nhops=nhops+1
291      lsthop(1,nhops)=lda(3,i)
292      lsthop(2,nhops)=lda(8,i)
293      timhop(nhops)=stime
294      lseq(lda(3,i))=lda(4,i)
295      lseq(lda(8,i))=lda(9,i)
296      do 8 j=1,nsaloc
297      if(lda(1,i).eq.isgan(j)) then
298      call qhop_dsite(isgan(j),issgm(j),ndxrel)
299c      xs(j,1)=0.9d0*xs(j+ndxrel,1)+0.1d0*xs(j,1)
300c      xs(j,2)=0.9d0*xs(j+ndxrel,2)+0.1d0*xs(j,2)
301c      xs(j,3)=0.9d0*xs(j+ndxrel,3)+0.1d0*xs(j,3)
302      ishop(j)=ior(ishop(j),1)
303      if(lda(5,i).ne.0) then
304      x(1)=xs(j,1)
305      x(2)=xs(j,2)
306      x(3)=xs(j,3)
307      xs(j,1)=xs(j+lda(5,i),1)
308      xs(j,2)=xs(j+lda(5,i),2)
309      xs(j,3)=xs(j+lda(5,i),3)
310      xs(j+lda(5,i),1)=x(1)
311      xs(j+lda(5,i),2)=x(2)
312      xs(j+lda(5,i),3)=x(3)
313      k=ishop(j)
314      ishop(j)=ishop(j+lda(5,i))
315      ishop(j+lda(5,i))=k
316      endif
317      endif
318      if(lda(6,i).eq.isgan(j)) then
319      xs(j,1)=rda(1,i)
320      xs(j,2)=rda(2,i)
321      xs(j,3)=rda(3,i)
322      ishop(j)=ieor(ishop(j),1)
323      if(lda(10,i).ne.0) then
324      x(1)=xs(j,1)
325      x(2)=xs(j,2)
326      x(3)=xs(j,3)
327      xs(j,1)=xs(j+lda(10,i),1)
328      xs(j,2)=xs(j+lda(10,i),2)
329      xs(j,3)=xs(j+lda(10,i),3)
330      xs(j+lda(10,i),1)=x(1)
331      xs(j+lda(10,i),2)=x(2)
332      xs(j+lda(10,i),3)=x(3)
333      k=ishop(j)
334      ishop(j)=ishop(j+lda(10,i))
335      ishop(j+lda(10,i))=k
336      endif
337      endif
338    8 continue
339      cf_hop=.true.
340      endif
341    7 continue
342c
343      return
344      end
345