1      subroutine prp_step(mdstep,stime,eww,esw,ess,fss,esk,
2     + epme,uqmd,uqmmm)
3c
4c $Id$
5c
6      implicit none
7c
8#include "prp_common.fh"
9#include "msgids.fh"
10#include "mafdecls.fh"
11#include "global.fh"
12c
13      logical frequency
14      external frequency
15c
16      integer mdstep
17      real*8 eww(mpe,2),esw(msf,mpe,2),ess(msf,msf,mpe,2),epme(3)
18      real*8 fss(msf,msf,3,2)
19      real*8 stime,esk(msf),uqmd,uqmmm
20c
21      if(me.eq.0) call prp_stat(mdstep,stime,eww,esw,ess,esk,epme(iset))
22c
23      if(frequency(mdstep,npener)) then
24      call cf_print_energy(lfnout)
25      endif
26c
27      if(frequency(mdstep,nfprop)) call prp_record()
28c
29      return
30      end
31      subroutine prp_proper(mdstep,stime,eww,esw,ess,fss,esk,epme,
32     + uqmd,uqmmm,epot,epotw,epotsw,epots,vol,dwr,dsr,ekin,etot,
33     + npolit,gsm,esa,box,xsm)
34c
35      implicit none
36c
37#include "prp_common.fh"
38#include "msgids.fh"
39#include "mafdecls.fh"
40#include "global.fh"
41c
42      external timer_wall,timer_wall_total
43      real*8 timer_wall,timer_wall_total
44c
45      integer mdstep
46      real*8 eww(mpe,2),esw(msf,mpe,2),ess(msf,msf,mpe,2),epme(3)
47      real*8 fss(msf,msf,3,2)
48      real*8 epot,epotw,epotsw,epots,vol,ekin,etot,epmec,ubias
49      real*8 eq,el,ep,ek,et,ewp,uqmd,uqmmm,box(3)
50      real*8 tempw,temps,stime,esk(msf),ewk,ep2,ep3,edrs,epmf
51      integer nwwl,nwws,nswl,nsws,nssl,nsss,nshitw,nshits,npolit
52      real*8 dwr,dsr(msm),gsm(msm,4,2),esa(nsa,2)
53      real*8 prest(3,3),virt(3,3),ekct(3,3),ep2m,ep3m,epmec2,epmec3
54      real*8 xsm(msm,3),dx(3),upmf(100)
55c
56      integer i,j,k,it,lenp
57c
58      call timer_start(54)
59c
60      if(ltwin) then
61      do 1 i=1,mpe
62      eww(i,1)=eww(i,1)+eww(i,2)
63    1 continue
64      do 2 i=1,msf
65      do 3 j=1,8
66      esw(i,j,1)=esw(i,j,1)+esw(i,j,2)
67    3 continue
68    2 continue
69      do 4 i=1,msf
70      do 5 j=1,msf
71      do 6 k=1,8
72      ess(i,j,k,1)=ess(i,j,k,1)+ess(i,j,k,2)
73    6 continue
74      fss(i,j,1,1)=fss(i,j,1,1)+fss(i,j,1,2)
75      fss(i,j,2,1)=fss(i,j,2,1)+fss(i,j,2,2)
76      fss(i,j,3,1)=fss(i,j,3,1)+fss(i,j,3,2)
77    5 continue
78    4 continue
79      endif
80c
81      el=eww(5,1)+eww(7,1)
82      eq=eww(6,1)+eww(8,1)
83      ep=zero
84      ewp=zero
85      do 9 i=1,mpe
86      ep=ep+eww(i,1)
87      ewp=ewp+eww(i,1)
88      do 10 j=1,msf
89      ep=ep+esw(j,i,1)
90      do 11 k=1,msf
91      ep=ep+ess(k,j,i,1)
92   11 continue
93   10 continue
94    9 continue
95      epmec=eww(9,1)
96      epmec2=eww(9,1)
97      epmec3=eww(9,1)
98      do 12 j=1,msf
99      el=el+esw(j,5,1)
100      eq=eq+esw(j,6,1)
101      do 13 k=1,msf
102      el=el+ess(k,j,5,1)+ess(k,j,7,1)
103      eq=eq+ess(k,j,6,1)+ess(k,j,8,1)
104      epmec=epmec+ess(k,j,9,1)
105      epmec2=epmec2+ess(k,j,10,1)
106      epmec3=epmec3+ess(k,j,11,1)
107   13 continue
108   12 continue
109c
110      if(me.eq.0) ep=ep+uqmd+uqmmm
111c
112c     if using cafe get properties from it
113c
114      if(lcafe) then
115      call cf_proper(volume,temp,tempw,temps,pres,tmpscl,tmsscl,prsscl,
116     + ewk,nwwl,nwws,nswl,nsws,nssl,nsss,nshitw,nshits,
117     + ep2,ep3,ep2m,ep3m,edrs,epmf,virial,prest,virt,ekct,ubias,upmf)
118      endif
119c
120c     if using space get properties from it
121c
122      if(lspac) then
123      endif
124c
125      if(me.eq.0) ep=ep+edrs
126      ep=ep+epmf
127c
128      if(lpme) then
129      ep=ep+epme(iset)
130      eq=eq+epme(iset)
131      ep2=ep2+epme(2)-epme(1)+epmec2-epmec
132      ep3=ep3+epme(3)-epme(1)+epmec3-epmec
133      endif
134c
135      ek=ewk
136      do 111 i=1,msf
137      ek=ek+esk(i)
138  111 continue
139      et=ep+ek
140c
141c     fill the property vector
142c
143      do 7 i=1,maxpro
144      p(i)=zero
145    7 continue
146c
147      if(me.eq.0) then
148c
149      if(ntype.eq.0) call sp_gettp(temp,pres)
150c
151      p(1)=dble(mdstep)
152      p(2)=stime
153      endif
154      p(3)=dble(nwws)
155      p(4)=dble(nwwl)
156      p(5)=dble(nsws)
157      p(6)=dble(nswl)
158      p(7)=dble(nsss)
159      p(8)=dble(nssl)
160      p(9)=dble(nshitw)
161      p(10)=dble(nshits)
162      p(11)=dble(npolit)
163      p(33)=zero
164      p(34)=zero
165      if(me.eq.0) then
166      p(12)=volume
167      p(13)=1.6605655d0*wbox/volume
168      p(14)=pres
169      p(15)=prsscl
170      p(16)=temp
171      p(17)=tempw
172      p(18)=temps
173      p(21)=tmpscl
174      p(22)=tmsscl
175      p(27)=ek
176      p(32)=ek
177      p(33)=virial
178      p(34)=pres*volume
179      p(64)=ewk
180      p(67)=ewk
181      endif
182      p(24)=eq
183      p(25)=el
184      p(26)=ep
185      p(32)=p(32)+ep
186      p(36)=ep2
187      p(37)=ep3
188      p(38)=ep2+ep2m
189      p(39)=ep3+ep3m
190      if(me.eq.0) p(40)=dfree
191      p(50)=ubias
192      p(52)=eww(6,1)
193      p(53)=eww(5,1)
194      p(54)=eww(8,1)
195      p(55)=eww(7,1)
196      p(56)=eww(1,1)
197      p(57)=eww(2,1)
198      p(58)=eww(13,1)
199      p(59)=eww(3,1)
200      p(60)=eww(4,1)
201      p(67)=ewp
202      p(68)=p(68)+ewp
203      p(70)=dwr
204      if(lnoe) p(75)=edrs
205      if(lpmf) p(76)=epmf
206      if(ntype.ne.3) then
207      p(77)=epme(iset)
208      p(80)=epmec
209      else
210      p(77)=epme(1)
211      p(78)=epme(2)-p(77)
212      p(79)=epme(3)-p(77)
213      p(80)=epme(1)+epmec
214      p(81)=epme(2)+epmec2-p(80)
215      p(82)=epme(3)+epmec3-p(80)
216      endif
217      p(84)=timer_wall(202)
218      p(85)=timer_wall(203)
219      if(me.eq.0) p(86)=timer_wall(203)
220      if(me.eq.0) then
221      p(87)=virt(1,1)
222      p(88)=virt(1,2)
223      p(89)=virt(1,3)
224      p(90)=virt(2,1)
225      p(91)=virt(2,2)
226      p(92)=virt(2,3)
227      p(93)=virt(3,1)
228      p(94)=virt(3,2)
229      p(95)=virt(3,3)
230      p(96)=prest(1,1)
231      p(97)=prest(1,2)
232      p(98)=prest(1,3)
233      p(99)=prest(2,1)
234      p(100)=prest(2,2)
235      p(101)=prest(2,3)
236      p(102)=prest(3,1)
237      p(103)=prest(3,2)
238      p(104)=prest(3,3)
239      p(105)=ekct(1,1)
240      p(106)=ekct(1,2)
241      p(107)=ekct(1,3)
242      p(108)=ekct(2,1)
243      p(109)=ekct(2,2)
244      p(110)=ekct(2,3)
245      p(111)=ekct(3,1)
246      p(112)=ekct(3,2)
247      p(113)=ekct(3,3)
248      p(114)=box(1)
249      p(115)=box(2)
250      p(116)=box(3)
251      endif
252c
253      maxp=isprop
254      it=isprop
255c
256      if(nsf.gt.0) then
257      do 33 i=1,nsf
258      it=isprop+(i-1)*30
259      p(it+2)=ess(i,i,6,1)
260      p(it+3)=ess(i,i,5,1)
261      p(it+5)=esw(i,6,1)
262      p(it+6)=esw(i,5,1)
263      p(it+7)=ess(i,i,1,1)
264      p(it+8)=ess(i,i,2,1)
265      p(it+9)=ess(i,i,13,1)
266      p(it+10)=ess(i,i,3,1)
267      p(it+11)=ess(i,i,4,1)
268      p(it+14)=ess(i,i,7,1)
269      p(it+15)=ess(i,i,8,1)
270      p(it+16)=ess(i,i,5,1)
271      p(it+17)=ess(i,i,6,1)
272      if(me.eq.0) then
273      p(it+19)=esk(i)
274      endif
275      p(it+20)=ess(i,i,5,1)+ess(i,i,6,1)+ess(i,i,7,1)+ess(i,i,8,1)+
276     + ess(i,i,1,1)+ess(i,i,2,1)+
277     + ess(i,i,3,1)+ess(i,i,4,1)+half*(esw(i,6,1)+esw(i,5,1))
278      do 34 j=i+1,nsf
279      p(it+20)=p(it+18)+half*(ess(i,j,6,1)+ess(j,i,6,1)+
280     + ess(i,j,5,1)+ess(j,i,5,1))
281   34 continue
282   33 continue
283      it=it+30
284      if(nsf.gt.1) then
285      do 35 i=1,nsf-1
286      do 36 j=i+1,nsf
287      it=it+1
288      it=it+1
289      p(it)=ess(i,j,6,1)+ess(j,i,6,1)
290      it=it+1
291      p(it)=ess(i,j,5,1)+ess(j,i,5,1)
292      it=it+1
293      it=it+1
294      p(it)=fss(i,j,1,1)-fss(j,i,1,1)
295      it=it+1
296      p(it)=fss(i,j,2,1)-fss(j,i,2,1)
297      it=it+1
298      p(it)=fss(i,j,3,1)-fss(j,i,3,1)
299      it=it+1
300      p(it)=sqrt((fss(i,j,1,1)-fss(j,i,1,1))**2+
301     + (fss(i,j,2,1)-fss(j,i,2,1))**2+(fss(i,j,3,1)-fss(j,i,3,1))**2)
302   36 continue
303   35 continue
304      endif
305      maxp=it
306c
307      if(nsm.gt.0.and.npstat.gt.0) then
308      it=maxp
309      do 37 i=1,nsm
310      it=it+1
311      if(me.eq.0) p(it)=gsm(i,iset,1)
312      if(p(it).lt.tiny) lp(it)=.false.
313   37 continue
314      do 55 i=1,nsm
315      it=it+1
316      if(me.eq.0) p(it)=gsm(i,4,1)
317      if(p(it).lt.tiny) lp(it)=.false.
318   55 continue
319      do 155 i=1,nsm-1
320      do 156 j=i+1,nsm
321      dx(1)=xsm(i,1)-xsm(j,1)
322      dx(2)=xsm(i,2)-xsm(j,2)
323      dx(3)=xsm(i,3)-xsm(j,3)
324      if(lpbc) call cf_pbc(1,dx,1,dx,1,0,1,1)
325      it=it+1
326      if(me.eq.0) p(it)=sqrt(dx(1)**2+dx(2)**2+dx(3)**2)
327  156 continue
328  155 continue
329      endif
330c
331      endif
332c
333      if(npmfi.gt.0) then
334      do 356 i=1,npmfi
335      it=it+1
336      p(it)=upmf(i)
337  356 continue
338      endif
339c
340      call timer_stop(54)
341c
342      if(iprof.eq.1) then
343      do 56 i=1,55
344      it=it+1
345      p(it)=timer_wall_total(i)
346   56 continue
347      endif
348c
349      maxp=it
350c
351      if(npener.gt.0) call cf_add_esa(esa)
352      if(np.gt.0) then
353      call ga_dgop(mrg_d44,p(9),3,'max')
354      if(me.ne.0) then
355      p(9)=zero
356      p(10)=zero
357      p(11)=zero
358      endif
359      if(maxp+mpe*(1+msf*(1+msf)).gt.maxpro) then
360      call ga_dgop(mrg_d45,p,maxp,'+')
361      call ga_dgop(mrg_d40,eww,mpe,'+')
362      call ga_dgop(mrg_d41,esw,msf*mpe,'+')
363      call ga_dgop(mrg_d42,ess,msf*msf*mpe,'+')
364      else
365      lenp=maxp
366      do 501 i=1,mpe
367      p(lenp+i)=eww(i,1)
368  501 continue
369      lenp=lenp+mpe
370      do 502 i=1,mpe
371      do 503 j=1,msf
372      lenp=lenp+1
373      p(lenp)=esw(j,i,1)
374  503 continue
375  502 continue
376      do 504 i=1,mpe
377      do 505 j=1,msf
378      do 506 k=1,msf
379      lenp=lenp+1
380      p(lenp)=ess(k,j,i,1)
381  506 continue
382  505 continue
383  504 continue
384      call ga_dgop(mrg_d45,p,lenp,'+')
385      lenp=maxp
386      do 511 i=1,mpe
387      eww(i,1)=p(lenp+i)
388  511 continue
389      lenp=lenp+mpe
390      do 512 i=1,mpe
391      do 513 j=1,msf
392      lenp=lenp+1
393      esw(j,i,1)=p(lenp)
394  513 continue
395  512 continue
396      do 514 i=1,mpe
397      do 515 j=1,msf
398      do 516 k=1,msf
399      lenp=lenp+1
400      ess(k,j,i,1)=p(lenp)
401  516 continue
402  515 continue
403  514 continue
404      endif
405cxxxxx      call ga_dgop(mrg_d43,epme,3,'+')
406      if(npener.gt.0) call ga_dgop(mrg_d47,esa,2*nsa,'+')
407      endif
408c
409      if(p(85).ne.0.0d0) then
410      p(85)=(p(85)-p(84))/p(85)
411      endif
412c
413      if(temp.eq.zero.or.me.ne.0) then
414      p(117)=zero
415      p(118)=zero
416      p(119)=zero
417      p(120)=zero
418      p(51)=zero
419      else
420      p(117)=exp(-p(36)/(rgas*temp))
421      p(118)=exp(-p(37)/(rgas*temp))
422      p(119)=exp(-p(38)/(rgas*temp))
423      p(120)=exp(-p(39)/(rgas*temp))
424      p(51)=exp(-p(50)/(rgas*temp))
425      if(nbias.gt.0) p(41)=p(40)*p(51)
426      endif
427c
428      if(nwm.gt.0) then
429      dwr=p(70)/dble(nwm)
430      p(70)=dwr
431      if(stime.ne.0) then
432      p(71)=1.0d-6*dwr/(6.0d0*stime)
433      else
434      p(71)=zero
435      endif
436      endif
437c
438      epot=p(26)
439      epots=zero
440      epotsw=zero
441      epotw=eww(5,1)+eww(6,1)+eww(8,1)
442      do 40 i=1,msf
443      epotsw=epotsw+esw(i,5,1)+esw(i,6,1)+esw(i,8,1)
444      do 41 j=1,msf
445      epots=epots+ess(i,j,1,1)+ess(i,j,2,1)+ess(i,j,3,1)+ess(i,j,4,1)
446      epots=epots+ess(i,j,5,1)+ess(i,j,6,1)+ess(i,j,8,1)
447   41 continue
448   40 continue
449c
450      epots=zero
451      epotsw=zero
452      epotw=zero
453      do 42 i=1,mpe
454      epotw=epotw+eww(i,1)
455      do 38 j=1,msf
456      epotsw=epotsw+esw(j,i,1)
457      do 39 k=1,msf
458      epots=epots+ess(k,j,i,1)
459   39 continue
460   38 continue
461   42 continue
462c
463      vol=volume
464      ekin=p(27)
465      etot=ekin+epot
466c
467      return
468      end
469      subroutine prp_stat(mdstep,stime,eww,esw,ess,esk,epme)
470c
471      implicit none
472c
473#include "prp_common.fh"
474#include "msgids.fh"
475#include "mafdecls.fh"
476#include "global.fh"
477c
478      logical frequency
479      external frequency
480c
481      integer mdstep
482      real*8 eww(mpe,2),esw(msf,mpe,2),ess(msf,msf,mpe,2),epme
483      real*8 stime,esk(msf)
484c
485      integer i,j
486      character*10 pdate,ptime
487      real*8 facs,fact,tfacs,tfact,rt
488c
489      nsum=nsum+1
490      nsumt=nsumt+1
491      nsump=nsump+1
492      do 8 i=1,maxpro
493      if(abs(p(i)).lt.tiny) p(i)=zero
494      psum(i)=psum(i)+p(i)
495      p2sum(i)=p2sum(i)+p(i)*p(i)
496      pslop(i)=pslop(i)+stime*p(i)
497      psumt(i)=psumt(i)+p(i)
498      p2sumt(i)=p2sumt(i)+p(i)*p(i)
499      pslopt(i)=pslopt(i)+stime*p(i)
500      psump(i)=psump(i)+p(i)
501    8 continue
502      tsum=tsum+stime
503      t2sum=t2sum+stime*stime
504      tsumt=tsumt+stime
505      t2sumt=t2sumt+stime*stime
506c
507      if(frequency(mdstep,nfoutp)) then
508      if(.not.lhdr) then
509      call swatch(pdate,ptime)
510      write(lfnout,1000) pdate,ptime
511 1000 format(/,' MOLECULAR DYNAMICS TIME STEP INFORMATION',T110,2A10,//,
512     + '    Time      Temp    Pres      Volume Tscalw Tscals  Pscal ',
513     + '   U(ele)    U(vdW)      U(pot)     U(kin)     U(tot) ',/,
514     + '     ps         K      Pa        nm**3                      ',
515     + '   kJ/mol    kJ/mol      kJ/mol     kJ/mol     kJ/mol ',/)
516      if(lpstep) write(lfnout,1010)
517 1010 format(14X,
518     + '   U(bnd)     U(ang)     U(dih)     U(imp)  ',
519     + '   Ui(3rd)    Ui(non) ',
520     + '   Uw(ele)    Uw(vdW)    U(kin)     U(pot)     U(tot)',/,14X,
521     + '   kJ/mol     kJ/mol     kJ/mol     kJ/mol  ',
522     + '   kJ/mol     kJ/mol  ',
523     + '   kJ/mol     kJ/mol     kJ/mol     kJ/mol     kJ/mol  ',/)
524      lhdr=.true.
525      endif
526      write(lfnout,1001) stime,temp,pres,volume,p(21),p(22),p(15),
527     + p(24),p(25),p(26),p(27),p(32)
528 1001 format(1x,f10.5,0pf8.2,1pe9.2,0pf10.3,3f7.4,5(1pe11.4),i5,i7)
529      if(lpstep) then
530      if(nwm.gt.0) then
531      rt=one/nwm
532      write(lfnout,1011) rt*p(56),rt*p(57),rt*p(58),rt*p(59),rt*p(52),
533     + rt*p(53),rt*p(54),rt*p(55),rt*p(64),rt*p(66),rt*p(67)
534 1011 format('   solvent  ',11f11.2)
535      endif
536      do 1 i=1,nsf
537      j=isprop+(i-1)*30
538      rt=p(j+7)+p(j+8)+p(j+9)+p(j+10)+p(j+2)+p(j+3)+p(j+5)+p(j+6)
539      write(lfnout,1012) i,p(j+7),p(j+8),p(j+9),p(j+10),
540     + p(j+14)+p(j+15),p(j+2)+p(j+3),p(j+5),p(j+6),p(j+19),rt,rt+p(j+19)
541 1012 format('   solute',i3,11f11.2)
542    1 continue
543      endif
544      endif
545c
546      if(frequency(mdstep,nfstat)) then
547      call swatch(pdate,ptime)
548      write(lfnout,2000) pdate,ptime,nsum,nsumt
549 2000 format(/,' MOLECULAR DYNAMICS STATISTICAL INFORMATION',t110,2a10,
550     + //,t41,2(3X,'Statistics over last ',I8,' steps',2X),/,
551     + t41,2(3X,'Average',5X,'RMS fluct',5X,'Drift/ps',3X),/)
552      facs=one/dble(nsum)
553      fact=one/dble(nsumt)
554      tfacs=one/(t2sum-facs*tsum*tsum)
555      tfact=one/(t2sumt-fact*tsumt*tsumt)
556      do 15 i=1,nprop
557      j=ixp(i)
558      if(lp(j)) write(lfnout,2001) pronam(j)(1:39),
559     + psum(j)*facs,sqrt(abs((p2sum(j)-psum(j)*psum(j)*facs)*facs)),
560     + (pslop(j)-facs*psum(j)*tsum)*tfacs,
561     + psumt(j)*fact,sqrt(abs((p2sumt(j)-psumt(j)*psumt(j)*fact)*fact)),
562     + (pslopt(j)-fact*psumt(j)*tsumt)*tfact,
563     + pronam(j)(40:50)
564 2001 format(1x,a39,t41,2(3(1pe12.5,1x),1x),a11)
565      psum(j)=zero
566      p2sum(j)=zero
567      pslop(j)=zero
568   15 continue
569      tsum=zero
570      t2sum=zero
571      nsum=0
572      lhdr=.false.
573      endif
574c
575      return
576      end
577      subroutine prp_print()
578c
579      implicit none
580c
581#include "prp_common.fh"
582c
583      character*10 pdate,ptime
584      integer i,j
585c
586      if(me.ne.0) return
587c
588      call swatch(pdate,ptime)
589      write(lfnout,1000) pdate,ptime
590 1000 format(/,' SINGLE POINT PROPERTIES',t110,2a10,/)
591      do 1  i=1,nprop
592      j=ixp(i)
593      if(lp(j)) write(lfnout,1001) pronam(j)(1:39),p(j),pronam(j)(40:50)
594 1001 format(1x,a39,t41,1pe18.9,1x,a11)
595    1 continue
596c
597      return
598      end
599      subroutine prp_record()
600c
601      implicit none
602c
603#include "prp_common.fh"
604c
605      integer i
606      character*10 pdate,ptime
607c
608      if(me.ne.0) return
609c
610      if(.not.lfhdr) then
611      call swatch(pdate,ptime)
612      write(lfnprp,3000) nprop,pdate,ptime,np,npfft
613 3000 format(i7,1x,2a10,2i5)
614      write(lfnprp,3001) (pronam(ixp(i)),i=1,nprop)
615 3001 format(a50)
616      lfhdr=.true.
617      endif
618      write(lfnprp,3002)
619 3002 format('frame')
620      if(iprop.eq.0) then
621      write(lfnprp,3003) (p(ixp(i)),i=1,nprop)
622      else
623      write(lfnprp,3003) (psump(ixp(i))/dble(nsump),i=1,nprop)
624 3003 format(4(1pe12.5))
625      endif
626      nsump=0
627      do 1 i=1,nprop
628      psump(ixp(i))=zero
629    1 continue
630c
631      return
632      end
633      subroutine prp_header()
634c
635      implicit none
636c
637#include "prp_common.fh"
638c
639      lfhdr=.false.
640c
641      return
642      end
643      logical function prp_mcti_step(ida,lda)
644c
645      implicit none
646c
647#include "prp_common.fh"
648#include "mafdecls.fh"
649#include "msgids.fh"
650#include "global.fh"
651c
652      logical prp_mcti_acc
653      external prp_mcti_acc
654c
655      integer ida,lda
656c
657      real*8 fdata(28)
658      real*8 aver,drift,stderr,corerr,ratio
659      logical done
660c
661      lerror=ida.gt.lda
662c
663      call cf_mcti(fdata)
664c
665      call ga_dgop(mrg_d44,fdata,28,'+')
666c
667      if(me.eq.0) then
668      done=prp_mcti_acc(ida,dbl_mb(i_dfr),dbl_mb(i_dfrm),fdata,
669     + aver,drift,stderr,corerr,ratio)
670      endif
671c
672      if(np.gt.1) then
673      call ga_brdcst(mrg_d46,done,ma_sizeof(mt_log,1,mt_byte),0)
674      endif
675c
676      prp_mcti_step=done
677c
678      return
679      end
680      logical function prp_mcti_acc(ida,dfr,dfrm,fdata,
681     + aver,drift,stderr,corerr,ratio)
682c
683      implicit none
684c
685#include "prp_common.fh"
686c
687      integer ida
688      real*8 dfr(mda),dfrm(mda),fdata(28)
689      real*8 aver,drift,stderr,corerr,ratio,cerror
690c
691      integer i
692      logical done
693      real*8 dfrnom
694c
695      dfree=zero
696      dfrnom=zero
697      do 1 i=1,24
698      dfree=dfree+fdata(i)
699      if(i.ne.1.and.i.ne.13) dfrnom=dfrnom+fdata(i)
700      deriv(i)=deriv(i)+fdata(i)
701    1 continue
702      nderiv=nderiv+1
703c
704      dfr(ida)=dfree
705      dfrm(ida)=dfrnom
706      if(dfree.ne.dfrnom) lfreem=.true.
707      nda=ida
708c
709    2 continue
710c
711      if(lerror) then
712      call error(lauto,lappr,1000,dfr,ida,
713     + aver,drift,stderr,corerr,ratio)
714      cerror=corerr
715      if(.not.lauto) cerror=samrat*corerr
716      if(.not.lauto.and.cerror.lt.edacq) then
717      lauto=.true.
718      lappr=.true.
719      goto 2
720      else
721      if(lauto) samrat=ratio
722      lauto=.false.
723      lappr=.false.
724      endif
725      done=cerror.lt.edacq.and.drift.lt.ddacq
726      else
727      done=.false.
728      endif
729c
730      prp_mcti_acc=done
731      return
732      end
733      subroutine prp_mcti_run(rlambd,dlambd,ndec)
734c
735      implicit none
736c
737#include "prp_common.fh"
738#include "mafdecls.fh"
739c
740      real*8 rlambd,dlambd
741      integer ndec
742c
743      call prp_mcti_r(rlambd,dlambd,ndec,dbl_mb(i_dfr),dbl_mb(i_dfrm),
744     + psumt(16),psumt(117),psumt(118),psumt(119),psumt(120),
745     + psumt(51),psumt(41))
746c
747      return
748      end
749      subroutine prp_mcti_r(rlambd,dlambd,ndec,dfr,dfrm,taver,
750     + ep2ave,ep3ave,ep2avm,ep3avm,ebias,dfbias)
751c
752      implicit none
753c
754#include "prp_common.fh"
755c
756      real*8 rlambd,dlambd,dfr(mda),dfrm(mda)
757      real*8 taver,ep2ave,ep3ave,ebias,dfbias,ep2avm,ep3avm
758      integer i,ndec
759c
760      if(me.eq.0) then
761      write(lfngib,1000) nderiv,nda,rlambd,dlambd,ndec,nsa,
762     + nbias,ebias/dble(nsumt),lfreem
763 1000 format(2i7,2f12.6,2i8,i4,e20.12,4x,l1)
764      write(lfngib,1001) deriv
765 1001 format(4e20.12)
766      write(lfngib,1002) (dfr(i),i=1,nda)
767      write(lfngib,1002) (dfrm(i),i=1,nda)
768 1002 format(4e20.12)
769      write(lfngib,1003) nsumt,taver/dble(nsumt),
770     + ep2ave/dble(nsumt),ep3ave/dble(nsumt),dfbias/ebias,
771     + ep2avm/dble(nsumt),ep3avm/dble(nsumt)
772 1003 format(i10,/,4e20.12,/,2e20.12)
773      endif
774      if(ndec.gt.0) call cf_wrtgib(lfngib)
775c
776      return
777      end
778      subroutine prp_mcti(n,filnam)
779c
780      implicit none
781c
782#include "prp_common.fh"
783#include "mafdecls.fh"
784#include "global.fh"
785c
786      integer n,i_dec,l_dec,ibl,i,igp_handle
787      character*80 filnam,string
788      character*5 gid
789c
790      if(npg.gt.1) then
791      if(me.eq.0) then
792      close(unit=lfngib,status='keep')
793      endif
794c
795c     now go global
796c
797      igp_handle=ga_pgroup_get_default()
798      call ga_pgroup_set_default(ga_pgroup_get_world())
799      call ga_sync()
800      if(ga_nodeid().eq.0) then
801c
802      ibl=index(filnam,' ')-1
803      string=filnam(1:ibl)//'.gib'
804      open(unit=lfngib,file=string,status='unknown',form='formatted')
805      rewind(lfngib)
806      do 1 i=1,npg
807      write(gid,3300) i-1
808 3300 format(i5.5)
809      string=filnam(1:ibl)//gid//'.gib'
810      open(unit=75,file=string,status='old',form='formatted')
811      rewind(unit=75)
812    2 continue
813      read(75,3301,end=99,err=99) string
814      write(lfngib,3301) string
815 3301 format(a)
816      goto 2
817   99 continue
818      close(unit=75)
819    1 continue
820      close(unit=lfngib,status='keep')
821c
822      endif
823c
824c     go back to groups
825c
826      call ga_sync()
827      call ga_pgroup_set_default(igp_handle)
828c
829      else if(ipg.eq.npg-1.or.npg.eq.1) then
830      if(me.eq.0) then
831      close(unit=lfngib,status='keep')
832      endif
833      endif
834c
835      if(ipg.eq.npg-1.or.npg.eq.1) then
836c
837      if(me.eq.0) then
838      ibl=index(filnam,' ')-1
839      string=filnam(1:ibl)//'.gib'
840      open(unit=lfngib,file=string,status='old',form='formatted')
841      rewind(lfngib)
842      endif
843c
844      npgdec=n
845c
846      if(npgdec.gt.0) then
847      if(.not.ma_push_get(mt_dbl,6*nsa,'dec',l_dec,i_dec))
848     + call md_abort('Failed to allocate dec',0)
849      else
850      if(.not.ma_push_get(mt_dbl,1,'dec',l_dec,i_dec))
851     + call md_abort('Failed to allocate dec',0)
852      endif
853c
854      call prp_mcti_s(dbl_mb(i_dfr),dbl_mb(i_dfrm),dbl_mb(i_dec))
855c
856      if(.not.ma_pop_stack(l_dec))
857     + call md_abort('Failed to deallocate dec',0)
858c
859      endif
860c      endif
861c
862      return
863      end
864      subroutine prp_mcti_s(dfr,dfrm,dec)
865c
866      implicit none
867c
868#include "prp_common.fh"
869#include "mafdecls.fh"
870c
871      real*8 dfr(mda),dfrm(mda),dec(6,nsa)
872c
873      character*10 pdate,ptime
874      integer i,j,k,number,ndec
875      real*8 rlambd,dlambd,ddrft,dsterr,dcerr,ratio
876      real*8 freeti,errti,drftti,taver,ep2ave,ep3ave,ebias,dfbias
877      real*8 ep2avm,ep3avm,dtmp
878      real*8 epr,epf,etp,slambd,rnum,fterm(24),gbias,freeb,freem
879      logical lfrm
880c
881      gbias=zero
882c
883      if(me.eq.0) then
884c
885      do 1 k=1,8
886c
887      if(k.le.4.and.npgdec.eq.0) goto 1
888c
889      call swatch(pdate,ptime)
890      if(k.eq.1) then
891      write(lfnout,1001)
892 1001 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION ',
893     + 'DECOMPOSITION',//,' Solvent derivatives in kJ/mol',//,
894     + ' Run Lambda',
895     + '      Mass Solvnt LJ Solute LJ Solvnt el Solute el     Bonds',
896     + '  Constrts    Angles Dihedrals Impropers Slvnt pol Solut pol',/)
897      elseif(k.eq.2) then
898      write(lfnout,1002)
899 1002 format(//,' Solute derivatives in kJ/mol',//,
900     + ' Run Lambda',
901     + '      Mass Solvnt LJ Solute LJ Solvnt el Solute el     Bonds',
902     + '  Constrts    Angles Dihedrals Impropers Slvnt pol Solut pol',/)
903      elseif(k.eq.3) then
904      write(lfnout,1003)
905 1003 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION ',
906     + 'DECOMPOSITION',//,' Solvent contributions in kJ/mol',//,
907     + ' Run Lambda',
908     + '      Mass Solvnt LJ Solute LJ Solvnt el Solute el     Bonds',
909     + '  Constrts    Angles Dihedrals Impropers Slvnt pol Solut pol',/)
910      elseif(k.eq.4) then
911      write(lfnout,1004)
912 1004 format(//,' Solute contributions in kJ/mol',//,
913     + ' Run Lambda',
914     + '      Mass Solvnt LJ Solute LJ Solvnt el Solute el     Bonds',
915     + '  Constrts    Angles Dihedrals Impropers Slvnt pol Solut pol',/)
916      elseif(k.eq.5) then
917      if(lfreem) then
918      write(lfnout,1005) pdate,ptime
919 1005 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION ',
920     + 'EXCLUDING MASS CONTRIBUTIONS',t110,2a10,//,
921     + ' Run Lambda   Size',
922     + '  Derivative  Derivative  Derivative Lambda',
923     + ' Free Energy Free Energy Free Energy Free Energy  Sampling',/,
924     + 18x,
925     + '     Average       Error       Drift       ',
926     + ' Accumulated       Error       Drift   Corrected     Ratio',/,
927     + 18x,
928     + '      kJ/mol      kJ/mol   kJ/mol ps       ',
929     + '      kJ/mol      kJ/mol   kJ/mol ps      kJ/mol',/)
930      else
931      write(lfnout,1006) pdate,ptime
932 1006 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION',
933     + t110,2a10,//,
934     + ' Run Lambda   Size',
935     + '  Derivative  Derivative  Derivative Lambda',
936     + ' Free Energy Free Energy Free Energy Free Energy  Sampling',/,
937     + 18x,
938     + '     Average       Error       Drift       ',
939     + ' Accumulated       Error       Drift   Corrected     Ratio',/,
940     + 18x,
941     + '      kJ/mol      kJ/mol   kJ/mol ps       ',
942     + '      kJ/mol      kJ/mol   kJ/mol ps      kJ/mol',/)
943      endif
944      elseif(k.eq.6) then
945      if(lfreem) then
946      write(lfnout,1007) pdate,ptime
947 1007 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION ',
948     + 'INCLUDING MASS CONTRIBUTIONS',t110,2a10,//,
949     + ' Run Lambda   Size',
950     + '  Derivative  Derivative  Derivative Lambda',
951     + ' Free Energy Free Energy Free Energy Free Energy  Sampling',/,
952     + 18x,
953     + '     Average       Error       Drift       ',
954     + ' Accumulated       Error       Drift   Corrected     Ratio',/,
955     + 18x,
956     + '      kJ/mol      kJ/mol   kJ/mol ps       ',
957     + '      kJ/mol      kJ/mol   kJ/mol ps      kJ/mol',/)
958      endif
959      elseif(k.eq.7) then
960      if(lfreem) then
961      write(lfnout,1008) pdate,ptime
962 1008 format(//' MULTISTEP THERMODYNAMIC PERTURBATION ',
963     + 'EXCLUDING MASS CONTRIBUTIONS',t110,2a10,//,
964     + ' Run Lambda   Size',
965     + ' Temperature     Reverse     Forward    Ensemble Lambda',
966     + ' Accumulated        Bias   Corrected',/,18x,
967     + '           K      kJ/mol      kJ/mol      kJ/mol       ',
968     + '      kJ/mol      kJ/mol      kJ/mol',/)
969      else
970      write(lfnout,1009) pdate,ptime
971 1009 format(//' MULTISTEP THERMODYNAMIC PERTURBATION',
972     + t110,2a10,//,
973     + ' Run Lambda   Size',
974     + ' Temperature     Reverse     Forward    Ensemble Lambda',
975     + ' Accumulated        Bias   Corrected',/,18x,
976     + '           K      kJ/mol      kJ/mol      kJ/mol       ',
977     + '      kJ/mol      kJ/mol      kJ/mol',/)
978      endif
979      else
980      if(lfreem) then
981      write(lfnout,1010) pdate,ptime
982 1010 format(//' MULTISTEP THERMODYNAMIC PERTURBATION ',
983     + 'INCLUDING MASS CONTRIBUTIONS',t110,2a10,//,
984     + ' Run Lambda   Size',
985     + ' Temperature     Reverse     Forward    Ensemble Lambda',
986     + ' Accumulated        Bias   Corrected',/,18x,
987     + '           K      kJ/mol      kJ/mol      kJ/mol       ',
988     + '      kJ/mol      kJ/mol      kJ/mol',/)
989      endif
990      endif
991c
992      rewind(lfngib)
993c
994      freeti=zero
995      freeb=zero
996      freem=zero
997      errti=zero
998      drftti=zero
999      etp=zero
1000      slambd=zero
1001      gbias=zero
1002      do 2 i=1,24
1003      fterm(i)=zero
1004    2 continue
1005      if(npgdec.gt.0) then
1006      do 22 i=1,nsa
1007      dec(1,i)=zero
1008      dec(2,i)=zero
1009      dec(3,i)=zero
1010      dec(4,i)=zero
1011      dec(5,i)=zero
1012      dec(6,i)=zero
1013   22 continue
1014      endif
1015      do 3 i=1,mrun
1016      read(lfngib,2000) nderiv,nda,rlambd,dlambd,ndec,nsa,nbias,ebias,
1017     + lfrm
1018      if(lfrm) lfreem=.true.
1019 2000 format(2i7,2f12.6,2i8,i4,e20.12,4x,l1)
1020      read(lfngib,2001) deriv
1021 2001 format(4e20.12)
1022      if(mda.lt.nda) then
1023      read(lfngib,2002) (dtmp,j=1,nda)
1024      read(lfngib,2002) (dtmp,j=1,nda)
1025      else
1026      read(lfngib,2002) (dfr(j),j=1,nda)
1027      read(lfngib,2002) (dfrm(j),j=1,nda)
1028      endif
1029 2002 format(4e20.12)
1030      read(lfngib,2003) number,taver,ep2ave,ep3ave,dfbias,
1031     + ep2avm,ep3avm
1032 2003 format(i10,/,4e20.12,/,2e20.12)
1033      if(nbias.gt.0) gbias=-rgas*taver*log(ebias)
1034      if(ndec.gt.0) call cf_rdgib(lfngib,dec,dlambd/dble(ndec))
1035c
1036      slambd=slambd+dlambd
1037      rnum=dlambd/dble(nderiv)
1038c
1039      do 4 j=1,24
1040      fterm(j)=fterm(j)+rnum*deriv(j)
1041    4 continue
1042c
1043      freem=freem+rnum*(deriv(1)+deriv(13))
1044c
1045      if(k.eq.1) then
1046      rnum=one/dble(nderiv)
1047      write(lfnout,1101) i,slambd,(rnum*deriv(j),j=1,12)
1048 1101 format(i4,f7.3,12f10.3)
1049      elseif(k.eq.2) then
1050      rnum=one/dble(nderiv)
1051      write(lfnout,1102) i,slambd,(rnum*deriv(j),j=13,24)
1052 1102 format(i4,f7.3,12f10.3)
1053      elseif(k.eq.3) then
1054      write(lfnout,1103) i,slambd,(fterm(j),j=1,12)
1055 1103 format(i4,f7.3,12f10.5)
1056      elseif(k.eq.4) then
1057      write(lfnout,1104) i,slambd,(fterm(j),j=13,24)
1058 1104 format(i4,f7.3,12f10.5)
1059      elseif(k.eq.5) then
1060      call error(.true.,.true.,1000,dfrm,nda,
1061     + dfree,ddrft,dsterr,dcerr,ratio)
1062c
1063      freeti=freeti+dfree*dlambd
1064      freeb=freeb+dfbias*dlambd
1065      errti=errti+(dcerr*dlambd)**2
1066      drftti=drftti+ddrft*dlambd
1067c
1068      write(lfnout,1105) i,rlambd,nderiv,dfree,dcerr,ddrft/tstep,
1069     + slambd,freeti,sqrt(errti),drftti/tstep,freeb,ratio
1070 1105 format(i4,f7.3,i7,f12.3,2f12.6,f7.3,f12.3,2f12.6,f12.3,f10.5)
1071c
1072      elseif(k.eq.6) then
1073      if(lfreem) then
1074      call error(.true.,.true.,1000,dfr,nda,
1075     + dfree,ddrft,dsterr,dcerr,ratio)
1076c
1077      freeti=freeti+dfree*dlambd
1078      freeb=freeb+dfbias*dlambd
1079      errti=errti+(dcerr*dlambd)**2
1080      drftti=drftti+ddrft*dlambd
1081c
1082      write(lfnout,1106) i,rlambd,nderiv,dfree,dcerr,ddrft/tstep,
1083     + slambd,freeti,sqrt(errti),drftti/tstep,freeb,ratio
1084 1106 format(i4,f7.3,i7,f12.3,2f12.6,f7.3,f12.3,2f12.6,f12.3,f10.5)
1085      endif
1086c
1087      elseif(k.eq.7) then
1088      epr=-rgas*taver*log(ep2ave)
1089      epf=-rgas*taver*log(ep3ave)
1090      etp=etp-epr+epf
1091      write(lfnout,1107) i,rlambd,number,taver,epr,epf,epf-epr,
1092     + slambd,etp,-gbias,etp-gbias
1093 1107 format(i4,f7.3,i7,4f12.3,f7.3,3f12.3)
1094      else
1095      if(lfreem) then
1096      epr=-rgas*taver*log(ep2avm)
1097      epf=-rgas*taver*log(ep3avm)
1098      etp=etp-epr+epf
1099      write(lfnout,1108) i,rlambd,number,taver,epr,epf,epf-epr,
1100     + slambd,etp,-gbias,etp-gbias
1101 1108 format(i4,f7.3,i7,4f12.3,f7.3,3f12.3)
1102      endif
1103      endif
1104c
1105    3 continue
1106c
1107      if(ndec.gt.0.and.k.eq.4) call cf_print_deco(lfnout,dec)
1108c
1109    1 continue
1110c
1111      close(unit=lfngib)
1112c
1113      endif
1114c
1115      return
1116      end
1117      subroutine prp_wrtmro(lfnmro,ndec)
1118c
1119      implicit none
1120c
1121#include "prp_common.fh"
1122#include "mafdecls.fh"
1123c
1124      integer lfnmro,ndec
1125c
1126      call prp_wtmro(lfnmro,ndec,dbl_mb(i_dfr),dbl_mb(i_dfrm))
1127c
1128      return
1129      end
1130      logical function prp_rdmri(lfnmri,ndec,mropt)
1131c
1132      implicit none
1133c
1134#include "prp_common.fh"
1135#include "mafdecls.fh"
1136c
1137      logical prp_rmri
1138      external prp_rmri
1139c
1140      integer lfnmri,ndec,mropt
1141c
1142      prp_rdmri=prp_rmri(lfnmri,ndec,mropt,dbl_mb(i_dfr),dbl_mb(i_dfrm))
1143c
1144      return
1145      end
1146      subroutine prp_wtmro(lfnmro,ndec,dfr,dfrm)
1147c
1148      implicit none
1149c
1150#include "prp_common.fh"
1151c
1152      integer lfnmro,ndec
1153      real*8 dfr(mda),dfrm(mda)
1154c
1155      integer i
1156c
1157      if(me.eq.0) then
1158      write(lfnmro) nderiv,nda,nprop,nsum,nsumt,ndec,maxp
1159      write(lfnmro) deriv
1160      write(lfnmro) (dfr(i),i=1,nda)
1161      write(lfnmro) (dfrm(i),i=1,nda)
1162      write(lfnmro) tsum,t2sum,tsumt,t2sumt
1163      write(lfnmro) (psum(i),i=1,maxp)
1164      write(lfnmro) (p2sum(i),i=1,maxp)
1165      write(lfnmro) (pslop(i),i=1,maxp)
1166      write(lfnmro) (psumt(i),i=1,maxp)
1167      write(lfnmro) (p2sumt(i),i=1,maxp)
1168      write(lfnmro) (pslopt(i),i=1,maxp)
1169      endif
1170      if(ndec.gt.0) call cf_wrtmro(lfnmro)
1171c
1172      return
1173      end
1174      logical function prp_rmri(lfnmri,ndec,mropt,dfr,dfrm)
1175c
1176      implicit none
1177c
1178#include "prp_common.fh"
1179c
1180      logical cf_rdmri
1181      external cf_rdmri
1182c
1183      integer lfnmri,ndec,mropt
1184      real*8 dfr(mda),dfrm(mda),dtmp
1185c
1186      integer i,nprp,nmxp
1187c
1188      read(lfnmri,err=9,end=9) nderiv,nda,nprp,nsum,nsumt,ndec,nmxp
1189      if(mropt.ne.2) then
1190      if(nprop.ne.nprp) then
1191      call md_abort('Number of properties nprop changed',0)
1192      endif
1193      if(maxp.ne.nmxp) then
1194      call md_abort('Number of properties maxp changed',0)
1195      endif
1196      endif
1197      read(lfnmri,err=9,end=9) deriv
1198      if(mda.lt.nda) then
1199      read(lfnmri,err=9,end=9) (dtmp,i=1,nda)
1200      read(lfnmri,err=9,end=9) (dtmp,i=1,nda)
1201      else
1202      read(lfnmri,err=9,end=9) (dfr(i),i=1,nda)
1203      read(lfnmri,err=9,end=9) (dfrm(i),i=1,nda)
1204      endif
1205      read(lfnmri,err=9,end=9) tsum,t2sum,tsumt,t2sumt
1206      if(mropt.ne.2) then
1207      read(lfnmri,err=9,end=9) (psum(i),i=1,maxp)
1208      read(lfnmri,err=9,end=9) (p2sum(i),i=1,maxp)
1209      read(lfnmri,err=9,end=9) (pslop(i),i=1,maxp)
1210      read(lfnmri,err=9,end=9) (psumt(i),i=1,maxp)
1211      read(lfnmri,err=9,end=9) (p2sumt(i),i=1,maxp)
1212      read(lfnmri,err=9,end=9) (pslopt(i),i=1,maxp)
1213      else
1214      do 1 i=1,6
1215      read(lfnmri,err=9,end=9) dtmp
1216    1 continue
1217      endif
1218      if(ndec.gt.0) then
1219      prp_rmri=cf_rdmri(lfnmri)
1220      else
1221      prp_rmri=.true.
1222      endif
1223      return
1224c
1225    9 continue
1226      prp_rmri=.false.
1227      return
1228      end
1229      integer function prp_dfr(ncopy)
1230c
1231      implicit none
1232c
1233#include "prp_common.fh"
1234#include "mafdecls.fh"
1235c
1236      integer prp_dfr_copy
1237      external prp_dfr_copy
1238c
1239      integer ncopy
1240c
1241      prp_dfr=prp_dfr_copy(ncopy,dbl_mb(i_dfr),dbl_mb(i_dfrm))
1242c
1243      return
1244      end
1245      integer function prp_dfr_copy(ncopy,dfr,dfrm)
1246c
1247      implicit none
1248c
1249#include "prp_common.fh"
1250c
1251      integer ncopy
1252      real*8 dfr(mda),dfrm(mda)
1253c
1254      integer i,ioff
1255c
1256      if(ncopy.lt.nda) then
1257      ioff=nda-ncopy
1258      do 1 i=1,ncopy
1259      dfr(i)=dfr(i+ioff)
1260      dfrm(i)=dfrm(i+ioff)
1261    1 continue
1262      nda=ncopy
1263      endif
1264c
1265      prp_dfr_copy=nda
1266      return
1267      end
1268