1      subroutine argos_cafe_start(irtdbi,lfnout,lfntop,filtop,
2     + ndistr,
3     + npmf,nipmf,nwmi,mwmi,nwai,mwai,nsfi,msfi,nsmi,msmi,nsai,msai,
4     + mdalg,npbt,nbxt,rcs,rcl,bx,
5     + jpme,jorder,jgx,jgy,jgz,nodpme,spmet,step,tols,mshw,mshs,nosh,
6     + iipolt,iitmps,iiprss,iipopt,rtmpx,rprsx,rtmpw,rtmps,rpres,
7     + sclq,fpmf,iislow,tempi,tempwi,tempsi,compr,ntyp,idset,ipset1,
8     + ipset2,issscl,delta,nfanal,lpbc,npgi,fldi,fvec,ffrq,npenrg,ictrl,
9     + nbiasi,mropti,incl,ltwn,nseqi,i_lseqi,nfhopi,rhopi,thopi,ndumsi,
10     + ipbtpi,lfnhopi,iradgi,nbgeti,npreci)
11c $Id$
12      implicit none
13c
14#include "argos_cafe_common.fh"
15#include "global.fh"
16#include "mafdecls.fh"
17#include "msgids.fh"
18c
19      integer lfnout,lfntop,irtdbi,lfnhopi
20      character*(*) filtop
21      integer ndistr,nwmi,mwmi,nwai,mwai,nsfi,msfi,nsmi,msmi,nsai,msai
22      integer jpme,jorder,jgx,jgy,jgz,mshw,mshs,nosh,npmf,nipmf,npgi
23      real*8 spmet,step,tols,bx(3),fldi,fvec(3),ffrq,sclq,rhopi,thopi
24      integer mdalg,npbt,nbxt,iipolt,nfanal,ictrl,incl
25      integer nfhopi,ipbtpi,iradgi,nbgeti,npreci
26      real*8 rcs,rcl,compr,delta,fpmf
27      integer iitmps,iiprss,iislow,iipopt,nseqi,i_lseqi
28      integer nodpme,ntyp,idset,ipset1,ipset2,ndumsi
29      integer issscl,npenrg,nbiasi,mropti
30      real*8 rtmpx,rprsx,rtmpw,rtmps,rpres,tempi,tempwi,tempsi
31      logical lpbc,ltwn
32c
33      integer i,itemp(2)
34      character*3 string
35c
36      call ga_sync()
37c
38      if(.not.ma_verify_allocator_stuff()) then
39      call md_abort('ERROR IN MEMORY USE AT START ARGOS_CAFE_START',0)
40      endif
41c
42      me=ga_nodeid()
43      np=ga_nnodes()
44c
45c     dimensions initially set to zero indicating non-allocated
46c
47      irtdb=irtdbi
48      lpress=lpbc
49c
50      itscal=iitmps
51      ipscal=iiprss
52      ipopt=iipopt
53      tmpext=rtmpx
54      prsext=rprsx
55      tmwrlx=rtmpw
56      tmsrlx=rtmps
57      prsrlx=rpres
58      scaleq=sclq
59      facpmf=fpmf
60      nbget=nbgeti
61      nprec=npreci
62c
63      ntype=ntyp
64      mdalgo=mdalg
65      includ=incl
66      ipbtyp=ipbtpi
67c
68      lfnhop=lfnhopi
69      nhops=0
70c
71      lfree=ntype.eq.3
72c
73      nbs=0
74      mscr=0
75      lscr=.false.
76      llst=.false.
77      lpair=.true.
78      llist=.false.
79      lpmf=npmf.gt.0
80      ndxp=0
81      nlda=0
82      maxl=0
83      lanal=nfanal.gt.0
84      npener=npenrg
85      icntrl=ictrl
86      iradgy=iradgi
87c
88      npgdec=npgi
89c
90      mwm=mwmi
91      mwa=mwai
92      msf=msfi
93      msm=msmi
94      msa=msai
95      nwm=nwmi
96      nwa=nwai
97      nsf=nsfi
98      nsm=nsmi
99      nsa=nsai
100      mscr=max(mwm+1,msa+1)
101      nwmtot=nwm
102      nsatot=nsa
103c
104      rshrt=rcs
105      rlong=rcl
106      rshrt2=rshrt*rshrt
107      rlong2=rlong*rlong
108      ltwin=ltwn
109      lssscl=issscl.ne.0
110c
111      nbxtyp=nbxt
112      npbtyp=npbt
113c
114      do 1 i=1,3
115      box(i)=bx(i)
116      boxh(i)=half*bx(i)
117    1 continue
118c
119      ipolt=iipolt
120      islow=iislow
121c
122      lstype=1
123c
124      ngc=1
125      ngl=1
126      nfrdf=99999
127      ifstep=1
128      ngrww=0
129      ngrsw=0
130      ngrss=0
131      ireact=0
132      iset=idset
133      issscl=0
134      nrwrec=0
135      isolvo=0
136c
137      ndums=ndumsi
138c
139      lpww=1
140      lpsw=1
141      lpss=1
142c
143      npww=1
144      npsw=1
145      npss=1
146      if(ltwin) then
147      npww=2
148      npsw=2
149      npss=2
150      endif
151c
152      mgc=1
153      mgl=1
154      mgr=1
155c
156      rffww=0.0d0
157      rffsw=0.0d0
158      rffss=0.0d0
159c
160      tstep=step
161      tstepi=one/tstep
162      tolsha=tols
163      mshitw=mshw
164      mshits=mshs
165      noshak=nosh
166c
167      temp=tempi
168      tempw=tempwi
169      temps=tempsi
170c
171      q14fac=0.833333d0
172      facpsc=compr*tstep/prsrlx
173c
174      field=fldi
175      fvect(1)=fvec(1)
176      fvect(2)=fvec(2)
177      fvect(3)=fvec(3)
178      ffreq=ffrq
179c
180      pi=four*atan(one)
181      twopi=two*pi
182c
183      wbox=zero
184c
185      shift0(1)=zero
186      shift0(2)=zero
187      shift0(3)=delta
188      shift0(4)=delta
189      shift0(5)=zero
190      shift0(6)=delta
191      shift1(1)=delta
192      shift1(2)=delta
193      shift1(3)=zero
194      shift1(4)=-delta
195      shift1(5)=delta
196      shift1(6)=zero
197c
198      numpmf=0
199c
200      rhop=rhopi
201      rhop2=rhop*rhop
202      thop=thopi
203      nfhop=nfhopi
204      nseq=nseqi
205      mseq=nseqi
206      i_lseq=i_lseqi
207c
208      lqhop=nfhop.ne.0
209c
210      if(lqhop.and.lfree)
211     + call md_abort('Proton hopping thermodynamics not allowed',0)
212c
213      ithint=ntype.eq.3
214      ipert2=ithint.or.(iset.eq.1.and.(ipset1.eq.2.or.ipset2.eq.2))
215      ipert3=ithint.or.(iset.eq.1.and.(ipset1.eq.2.or.ipset2.eq.2))
216      do 2 i=1,24
217      ith(i)=.false.
218      ip2(i)=.false.
219      ip3(i)=.false.
220    2 continue
221c
222      if(.not.ma_push_get(mt_byte,16*nsatot,'snam',l_snam,i_snam))
223     + call md_abort('Failed to allocate snam',me)
224c
225      if(lqhop) then
226      if(.not.ma_push_get(mt_int,mseq,'mprot',l_mprot,i_mprot))
227     + call md_abort('Failed to allocate mprot',me)
228      endif
229c
230      if(npgdec.gt.1) then
231      if(.not.ma_push_get(mt_dbl,6*nsatot,'sti',l_sti,i_sti))
232     + call md_abort('Failed to allocate sti',me)
233      else
234      if(.not.ma_push_get(mt_dbl,1,'sti',l_sti,i_sti))
235     + call md_abort('Failed to allocate sti',me)
236      endif
237c
238      call argos_cafe_rdtop(lfntop,filtop,byte_mb(i_snam))
239      call argos_cafe_topol_init(lfnout)
240c
241c     distance restraints
242c     -------------------
243c
244      if(ndistr.gt.0) then
245      if(me.eq.0) then
246      open(unit=lfntop,file=filtop(1:index(filtop,' ')-1),
247     + form='formatted',status='old',err=9999)
248      rewind(lfntop)
249    5 continue
250      read(lfntop,100,end=9999,err=9999) string
251  100 format(a3)
252      if(string.ne.'noe') goto 5
253      read(lfntop,1000) ndrs
254 1000 format(i5)
255      endif
256      call ga_brdcst(mcf_55,ndrs,ma_sizeof(mt_int,1,mt_byte),0)
257      if(ndrs.gt.0) then
258      if(.not.ma_push_get(mt_int,2*ndrs,'idrs',l_idrs,i_idrs))
259     + call md_abort('Failed to allocate idrs',0)
260      if(.not.ma_push_get(mt_dbl,6*ndrs,'rdrs',l_rdrs,i_rdrs))
261     + call md_abort('Failed to allocate rdrs',0)
262      if(.not.ma_push_get(mt_dbl,6*ndrs,'xdrs',l_xdrs,i_xdrs))
263     + call md_abort('Failed to allocate xdrs',0)
264      endif
265      call argos_cafe_rddrs(lfntop,int_mb(i_idrs),dbl_mb(i_rdrs))
266      if(me.eq.0) then
267      close(unit=lfntop)
268      endif
269      endif
270c
271c     proton hopping: donor-acceptor pair list allocation
272c     ---------------------------------------------------
273c
274      if(nhop.gt.0) then
275      if(.not.ma_push_get(mt_int,16*nhop*3,'lda',l_lda,i_lda))
276     + call md_abort('Failed to allocate lda',0)
277      if(.not.ma_push_get(mt_dbl,11*nhop*3,'rda',l_rda,i_rda))
278     + call md_abort('Failed to allocate rda',0)
279      if(.not.ma_push_get(mt_dbl,4*nhop*3,'uda',l_uda,i_uda))
280     + call md_abort('Failed to allocate uda',0)
281      if(.not.ma_push_get(mt_dbl,nhop*3,'pda',l_pda,i_pda))
282     + call md_abort('Failed to allocate pda',0)
283      if(.not.ma_push_get(mt_int,nhop*30,'lsthop',l_lsthop,i_lsthop))
284     + call md_abort('Failed to allocate lsthop',0)
285      if(.not.ma_push_get(mt_dbl,nhop*15,'timhop',l_timhop,i_timhop))
286     + call md_abort('Failed to allocate timhop',0)
287      endif
288c
289c     potential of mean force
290c     -----------------------
291c
292      if(lpmf) then
293      if(me.eq.0) then
294      open(unit=lfntop,file=filtop(1:index(filtop,' ')-1),
295     + form='formatted',status='old',err=9998)
296      rewind(lfntop)
297    6 continue
298      read(lfntop,100,end=9998,err=9998) string
299      if(string.ne.'pmf') goto 6
300      read(lfntop,2000) itemp
301 2000 format(2i5)
302      endif
303      call ga_brdcst(mcf_59,itemp,ma_sizeof(mt_int,2,mt_byte),0)
304      numpmf=itemp(1)
305      npmfa=itemp(2)
306      if(numpmf.gt.0) then
307      if(.not.ma_push_get(mt_int,8*numpmf,'ipmf',l_ipmf,i_ipmf))
308     + call md_abort('Failed to allocate ipmf',0)
309      if(.not.ma_push_get(mt_int,4*numpmf*npmfa,'jpmf',l_jpmf,i_jpmf))
310     + call md_abort('Failed to allocate jpmf',npmfa)
311      if(.not.ma_push_get(mt_dbl,18*numpmf,'rpmf',l_rpmf,i_rpmf))
312     + call md_abort('Failed to allocate rpmf',0)
313      if(.not.ma_push_get(mt_dbl,16*numpmf,'xpmf',l_xpmf,i_xpmf))
314     + call md_abort('Failed to allocate xpmf',0)
315      if(.not.ma_push_get(mt_dbl,12*numpmf,'ypmf',l_ypmf,i_ypmf))
316     + call md_abort('Failed to allocate ypmf',0)
317      if(.not.ma_push_get(mt_dbl,4*numpmf,'wpmf',l_wpmf,i_wpmf))
318     + call md_abort('Failed to allocate wpmf',0)
319      if(.not.ma_push_get(mt_dbl,numpmf,'upmf',l_upmf,i_upmf))
320     + call md_abort('Failed to allocate upmf',0)
321      endif
322      call argos_cafe_rdpmf(lfnout,lfntop,int_mb(i_ipmf),int_mb(i_jpmf),
323     + dbl_mb(i_rpmf))
324      call ga_brdcst(mcf_75,nbias,ma_sizeof(mt_int,1,mt_byte),0)
325      if(me.eq.0) then
326      close(unit=lfntop)
327      endif
328      endif
329c
330      if(ithint) then
331      do 3 i=1,24
332      ip2(i)=ith(i)
333      ip3(i)=ith(i)
334    3 continue
335      endif
336c
337c     particle-mesh Ewald initialization
338c     ----------------------------------
339c
340      ipme=jpme
341      morder=jorder
342      ngx=jgx
343      ngy=jgy
344      ngz=jgz
345      ngmax=max(ngx,ngy,ngz)
346      ngrx=ngx+morder
347      ngry=ngy+morder
348      ngrz=ngz
349      pmetol=spmet
350      if(ipme.gt.0) then
351      if(morder.gt.25) call md_abort('morder too large',0)
352      call argos_cafe_alpha
353      call argos_pme_start(alpha,morder,1,nodpme,
354     + ngx,ngy,ngz,mwm,mwa,msa,icntrl,nbget)
355      endif
356c
357      call argos_cafe_pardif(dbl_mb(i_mas),dbl_mb(i_vdw),dbl_mb(i_chg),
358     + int_mb(i_iwa),int_mb(i_iwq),
359     + mbt(1),numb(1),mbp(1),dbl_mb(i_bnd(1)),
360     + mht(1),numh(1),mhp(1),dbl_mb(i_ang(1)),
361     + mdt(1),numd(1),mdp(1),dbl_mb(i_dih(1)),
362     + mit(1),numi(1),mip(1),dbl_mb(i_imp(1)),
363     + mbt(2),mbp(2),dbl_mb(i_bnd(2)),mht(2),mhp(2),dbl_mb(i_ang(2)),
364     + mdt(2),mdp(2),dbl_mb(i_dih(2)),mit(2),mip(2),dbl_mb(i_imp(2)))
365c
366      if(me.eq.0.and.nsf.ne.nsfi) then
367      write(*,'(a,a)') ' Number of fractions differs on topology ',
368     + ' and restart files'
369      endif
370c
371c     initialize QHOP parameters
372c
373      if(nfhop.gt.0) call qhop_init(lfntop,filtop,lfnhop,me)
374c
375      nsfi=nsf
376      nbiasi=nbias
377      mropt=mropti
378      nipmf=npmfi
379c
380      return
381
382 9998 continue
383      call md_abort('Potentials of mean force input not found',0)
384 9999 continue
385      call md_abort('Distance restraints file not found',0)
386c
387      return
388      end
389