1      subroutine qhop_setup(lfnhdb,lfntop,
2     + mparms,nparms,matm,natm,mseq,nseq,
3     + catm,latm,cseq,lseq,mbnd,nbnd,lbnd,rbnd,mang,nang,lang,rang)
4c
5c $Id$
6c
7      implicit none
8c
9#include "qhop_common.fh"
10#include "mafdecls.fh"
11c
12      integer lfnhdb,lfntop,matm,natm,mseq,nseq,nparms,mparms
13      character*10 cseq(mseq)
14      character*6 catm(mparms,matm)
15      integer latm(11,matm),lseq(6,mseq)
16      integer mbnd,nbnd
17      integer lbnd(4,mbnd)
18      real*8 rbnd(nparms,2,mbnd)
19      integer mang,nang
20      integer lang(5,mang)
21      real*8 rang(nparms,2,mang)
22c
23      call qhop_setup2(lfnhdb,lfntop,
24     + mparms,nparms,matm,natm,mseq,nseq,
25     + catm,latm,cseq,lseq,mbnd,nbnd,lbnd,rbnd,mang,nang,lang,rang,
26     + int_mb(i_ptseq),int_mb(i_ptarat),int_mb(i_iarat),
27     + dbl_mb(i_racs),dbl_mb(i_racs+mxar),int_mb(i_ptpar),
28     + dbl_mb(i_par),dbl_mb(i_par+7*maxpar),dbl_mb(i_par+14*maxpar),
29     + dbl_mb(i_par+17*maxpar),dbl_mb(i_par+23*maxpar),
30     + dbl_mb(i_par+28*maxpar))
31c
32      return
33      end
34      subroutine qhop_setup2(lfnhdb,lfntop,
35     + mparms,nparms,matm,natm,mseq,nseq,
36     + catm,latm,cseq,lseq,mbnd,nbnd,lbnd,rbnd,mang,nang,lang,rang,
37     + ptseq,ptarat,arat,deq,aneq,ptpar,tunnel,defe12,zpef,
38     + tdsgl,tstval,e12fxy)
39c
40c $Id$
41c
42      implicit none
43c
44#include "qhop_common.fh"
45c
46c     input variables
47c     ---------------
48c
49c     matm            : dimension of atom arrays
50c     natm            : number of atoms
51c     mseq            : dimension of residue arrays
52c     nseq            : number of residues
53c     nparms          : number of parameter sets
54c     mparms          : number of parameter sets plus 1
55c     cseq            : residue names
56c     catm(1,1:natm)  : atom names
57c     catm(k,1:natm)  : atom types in parameter set k-1
58c     latm(5,1:natm)  : residue number
59c     latm(10,1:natm) : + heavy atom attached to protonatable hydrogen
60c                       - protonatable hydrogen, relative index to heavy atom
61c     lseq(5,1:nseq)  : number of protonation states of residue
62c     lseq(6,1:nseq)  : current protonation state of residue
63c
64c     mbnd            : dimension of bond index array
65c     nbnd            : number of bonds
66c     lbnd(1,1:nbnd)  : atom index i
67c     lbnd(2,1:nbnd)  : atom index j
68c
69      integer lfnhdb,lfntop,matm,natm,mseq,nseq,nparms,mparms
70      character*10 cseq(mseq)
71      character*6 catm(mparms,matm)
72      integer latm(11,matm),lseq(6,mseq)
73      integer mbnd,nbnd
74      integer lbnd(4,mbnd)
75      real*8 rbnd(nparms,2,mbnd)
76      integer mang,nang
77      integer lang(5,mang)
78      real*8 rang(nparms,2,mang)
79c
80      real*8 tunnel(maxpar,7),defe12(maxpar,7),zpef(maxpar,3),
81     + tdsgl(maxpar,6),tstval(maxpar,5),e12fxy(maxpar,3)
82      real*8 deq(mxar),aneq(mxar)
83      integer ptseq(maxseq,2),ptpar(maxpar),arat(mxar,4),ptarat(mxseq)
84c
85      integer i,j,k,l,m,n
86c      integer mxseq,maxhv,maxpar,mxar
87c      parameter(mxseq=30)
88c      parameter(maxhv=4)
89c      parameter(maxpar=10000)
90c      parameter(mxar=mxseq*30)
91c
92c      integer ptseq(mxseq,2),ptpar(maxpar),arat(mxar,4),seqat(matm)
93c      common/qh_ptrs/ptseq,ptpar,arat,seqat
94c
95c      integer shbit(5)
96c      common/sbit/shbit(5)
97      integer numseq,numpar,mnprot,mnhv,nuseq,larat,nbnat
98      integer prtin(2),flseq,flhv,nseql,eof,ptat(3)
99      integer numhv(mxseq),ihv,numprot(mxseq)
100      integer nbseq,nbhv,nbpt
101c      real*8 deq(mxar),aneq(mxar)
102      character*10 cuseq(mxseq),cseqin(2)
103      character*6 chv(mxseq,maxhv),chvin(2)
104      character*80 fileinp,fileout
105
106c      real*8 e12fxy(maxpar,3),defe12(maxpar,7),zpef(maxpar,6),
107c     + tdsgl(maxpar,6),tunnel(maxpar,7),tstval(maxpar,5)
108c
109      print*,'qhop_setup'
110c
111      fileinp='qhop.par'
112      open(unit=lfnhdb,file=fileinp(1:index(fileinp,' ')-1),err=112)
113c
114c     Initializing variables
115c
116      mnprot=0
117      mnhv=0
118      numpar=0
119      numseq=0
120      nuseq=0
121      flhv=0
122      flseq=0
123      ptat(1)=0
124      larat=0
125      do 31 i=1,mxar
126      do j=1,4
127      arat(i,j)=0
128      enddo
129      deq(i)=0
130      aneq(i)=0
131   31 continue
132      do 12 i=1,mxseq
133      ptseq(i,1)=0
134      ptseq(i,2)=0
135      ptarat(i)=0
136   12 continue
137      nbpt=0
138      nbhv=0
139      nbseq=0
140c
141c     Setting pointers to static array cuseq(numseq) and calculating max
142c     heavy atom # and prot states
143c
144c     print*,'prima loop'
145c      do 36 i=1,nbnd
146c      write(*,101)i,(lbnd(j,i),j=1,2),rbnd(1,1,i)
147c  101 format(3i5,f12.6)
148c   36 continue
149c      print*,'setup 1'
150c
151      do i=1,nang
152c      write(*,102)i,(lang(j,i),j=1,3),rang(1,1,i)
153c  102 format(4i5,f12.6)
154      enddo
155c
156      do 1 i=1,natm
157c
158c      write(*,1000) i,catm(1,i),latm(5,i),cseq(latm(5,i)),latm(10,i),
159c     + lseq(5,latm(5,i)),lseq(6,latm(5,i))
160c
161 1000 format(i5,1x,a6,1x,i5,1x,a10,3i5)
162c
163      if((i.gt.1.and.latm(5,i).ne.latm(5,i-1)).or.i.eq.1)then
164c      print*,'flhv=',flhv
165      if(flhv.eq.1)then
166      ptat(2)=i-1
167      flhv=0
168      endif
169      ptat(3)=i
170      endif
171c
172
173      if(latm(10,i).gt.0)then
174      flseq=1
175      j=1
176      do while (j.le.numseq.and.flseq.ne.0)
177      if(cuseq(j).eq.cseq(latm(5,i)))then
178      flseq=0
179      ptseq(latm(5,i),1)=j
180      ptseq(latm(5,i),2)=ptat(3)
181      endif
182      j=j+1
183      enddo
184c
185      if(flseq.ne.0)then
186c
187      flhv=1
188      if(lseq(5,latm(5,i)).gt.mnprot)mnprot=lseq(5,latm(5,i))
189      numseq=numseq+1
190      nseql=latm(5,i)
191      cuseq(numseq)=cseq(latm(5,i))
192      ptseq(latm(5,i),1)=numseq
193      ptseq(latm(5,i),2)=ptat(3)
194      numhv(numseq)=1
195      if(mnhv.eq.0)mnhv=1
196      chv(numseq,1)=catm(1,i)
197c
198      if(ptat(1).gt.0)then
199      ptarat(numseq-1)=larat+1
200c      print 1212,'ptat',numseq,ptat(1),ptat(2)
201c 1212 format(a5,3i5)
202      ihv=0
203      do 13 k=ptat(1),ptat(2)
204      if(latm(10,k).gt.0)ihv=ihv+1
205      nbnat=0
206      larat=larat+1
207      do 14 j=1,nbnd
208      if(latm(10,lbnd(1,j)).gt.0.or.latm(10,lbnd(2,j)).gt.0) then
209c      write(*,'(5i5)')
210c     + j,lbnd(1,j),lbnd(2,j),latm(10,lbnd(1,j)),latm(10,lbnd(2,j))
211      do 28 m=1,2
212      l=2**(2-m)
213      if(lbnd(m,j).eq.k.and.(latm(10,k).lt.0.or.latm(10,lbnd(l,j))
214     + .ge.0))then
215      nbnat=nbnat+1
216      arat(larat,nbnat)=lbnd(l,j)-k
217c      print 1029,'check ',numseq-1,k,lbnd(l,j),arat(larat,nbnat),
218c     + larat,nbnat
219c 1029 format(a6,6i6)
220      if(latm(10,k).lt.0)then
221      deq(larat)=rbnd(1,1,j)
222c      print 121,'larat ',numseq-1,k,arat(larat,1),
223c     + latm(10,k+arat(larat,1))
224c  121 format(a6,4i5)
225      arat(larat,nbnat+1)=latm(10,k+arat(larat,1))
226      arat(larat,nbnat+2)=ihv
227c
228      if(latm(10,k+1).ge.0)then
229      do 25 n=1,-arat(larat,nbnat)
230      arat(larat-n+1,nbnat+3)=-arat(larat,nbnat)
231   25 continue
232      endif
233c
234      if(arat(larat,nbnat+1).eq.1)then
235      do 18 n=1,nang
236c      if(k.eq.50)
237c     + print 222,'angle ',k,lbnd(l,j),k+arat(larat,1)
238c     + +arat(larat+arat(larat,1),1)
239  222 format(a6,1x,3i5)
240      if((lang(1,n).eq.k+arat(larat,1)+arat(larat+arat(larat,1),1)
241     + .and.lang(2,n).eq.lbnd(l,j).and.lang(3,n).eq.k).or.
242     + (lang(3,n).eq.k+arat(larat,1)+arat(larat+arat(larat,1),1)
243     + .and.lang(2,n).eq.lbnd(l,j).and.lang(1,n).eq.k))
244     + aneq(larat)=rang(1,1,n)
245   18 continue
246
247      endif
248      endif
249      endif
250   28 continue
251      endif
252   14 continue
253   13 continue
254      endif
255      ptat(1)=ptat(3)
256c
257      elseif(flhv.eq.1.and.latm(5,i).eq.nseql)then
258      numhv(numseq)=numhv(numseq)+1
259      if(numhv(numseq).gt.mnhv)mnhv=numhv(numseq)
260      chv(numseq,numhv(numseq))=catm(1,i)
261      else
262      flhv=0
263      endif
264      endif
265    1 continue
266c
267c     Writing array arat for last unique res
268c
269c      print*,'setup 2 ',ptat(1),ptat(2)
270      if(ptat(1).ge.0)then
271      ihv=0
272      ptarat(numseq)=larat+1
273
274      do 15 k=ptat(1),ptat(2)
275      if(latm(10,k).gt.0)ihv=ihv+1
276      nbnat=0
277      larat=larat+1
278      do 16 j=1,nbnd
279      if(latm(10,lbnd(1,j)).gt.0.or.latm(10,lbnd(2,j)).gt.0) then
280c      write(*,'(5i5)')
281c     + j,lbnd(1,j),lbnd(2,j),latm(10,lbnd(1,j)),latm(10,lbnd(2,j))
282      do 29 m=1,2
283      l=2**(2-m)
284      if(lbnd(m,j).eq.k.and.(latm(10,k).lt.0.or.latm(10,lbnd(l,j))
285     + .ge.0))then
286c      print 2012,k,latm(10,k),l,m,lbnd(l,j),lbnd(m,j)
287 2012 format('kappa ',6i5)
288      nbnat=nbnat+1
289      arat(larat,nbnat)=lbnd(l,j)-k
290      if(latm(10,k).lt.0)then
291      deq(larat)=rbnd(1,1,j)
292      arat(larat,nbnat+1)=latm(10,k+arat(larat,nbnat))
293      arat(larat,nbnat+2)=ihv
294
295      if(latm(10,k+1).ge.0)then
296      do 23 n=1,-arat(larat,nbnat)
297      arat(larat-n+1,nbnat+3)=-arat(larat,nbnat)
298   23 continue
299      endif
300c
301      if(arat(larat,nbnat+1).eq.1)then
302      do 19 n=1,nang
303      if((lang(1,n).eq.k+arat(larat,1)+arat(larat+arat(larat,1),1)
304     + .and.lang(2,n).eq.lbnd(l,j).and.lang(3,n).eq.k).or.
305     + (lang(3,n).eq.k+arat(larat,1)+arat(larat+arat(larat,1),1)
306     + .and.lang(2,n).eq.lbnd(l,j).and.lang(1,n).eq.k))
307     + aneq(larat)=rang(1,1,n)
308   19 continue
309      else if(arat(larat,nbnat+1).eq.3)then
310      aneq(larat)=(109.5/180)*3.14159
311      endif
312      endif
313      endif
314   29 continue
315      endif
316   16 continue
317   15 continue
318c      print*,'fine loop k'
319      endif
320
321c
322c     Calculating number of bits for ptrs to parameter sets
323c
324c      print*,'setup 2bis'
325      do while(2**nbpt.le.mnprot)
326      nbpt=nbpt+1
327      enddo
328      do while(2**(nbhv).le.mnhv)
329      nbhv=nbhv+1
330      enddo
331      do while(2**(nbseq).le.numseq)
332      nbseq=nbseq+1
333      enddo
334c
335      shbit(1)=nbpt
336      shbit(2)=shbit(1)+nbhv
337      shbit(3)=shbit(2)+nbseq
338      shbit(4)=shbit(3)+nbpt
339      shbit(5)=shbit(4)+nbhv
340c      print 1011,(shbit(i),i=1,5),mnprot,mnhv,numseq
341c 1011 format('shbit ',8i5)
342c
343c     reading in the relevant hopping parameters
344c
345      eof=1
346      do while (eof.ge.0)
347c
348c      print*,'setup 3'
349      read(lfnhdb,1001,end=8)cseqin(1),chvin(1),prtin(1),
350     + cseqin(2),chvin(2),prtin(2)
351 1001 format(a10,a6,i3,10x,a10,a6,i3)
352      if(cseqin(1)(1:1).eq.' ')goto 8
353      read(lfnhdb,1003) (e12fxy(numpar+1,i),i=1,3)
354 1003 format(3f12.6)
355      read(lfnhdb,1004) (defe12(numpar+1,i),i=1,7)
356 1004 format(3f12.6,f8.4,2f12.6,f12.8)
357      read(lfnhdb,1005) (zpef(numpar+1,i),i=1,6)
358 1005 format(6f12.6)
359      read(lfnhdb,1006) (tdsgl(numpar+1,i),i=1,6)
360 1006 format(6f12.6)
361      read(lfnhdb,1007) (tunnel(numpar+1,i),i=1,7)
362 1007 format(f8.2,f8.3,f10.5,e12.3,f10.5,2e12.3)
363      read(lfnhdb,1008) (tstval(numpar+1,i),i=1,5)
364 1008 format(f10.2,2f10.4,f10.6,f10.4)
365      do 4 i=1,numseq
366      if(cuseq(i).eq.cseqin(1))then
367      do 5 j=1,numseq
368      if(cuseq(j).eq.cseqin(2))then
369      do 6 k=1,numhv(i)
370      if(chvin(1).eq.chv(i,k))then
371      do 7 l=1,numhv(j)
372      if(chvin(2).eq.chv(j,l))then
373c     MANCA UN CONTROLLO SUI NUMERI CHE LEGGE SUI PROT STATES
374      numpar=numpar+1
375      ptpar(numpar)=ishft(i,shbit(5))+ishft(k,shbit(4))+
376     + ishft(prtin(1),shbit(3))+ishft(j,shbit(2))+
377     + ishft(l,shbit(1))+prtin(2)
378c      write(lfntop,1030) cseqin(1),chvin(1),prtin(1),cseqin(2),chvin(2)
379c     + ,prtin(2),ptpar(numpar)
380c 1030 format(a10,a6,i3,10x,a10,a6,i3,i10)
381c      write(lfntop,1009) i,k,prtin(1),j,l,prtin(2)
382c      write(lfntop,1003)(e12fxy(numpar,m),m=1,3)
383c      write(lfntop,1004)(defe12(numpar,m),m=1,7)
384c      write(lfntop,1005)(zpef(numpar,m),m=1,6)
385c      write(lfntop,1006)(tdsgl(numpar,m),m=1,6)
386c      write(lfntop,1007)(tunnel(numpar,m),m=1,7)
387c      write(lfntop,1008)(tstval(numpar,m),m=1,5)
388 1009 format(3i7,5x,3i7)
389      else
390c     ERR HEAVY ATOM2 NOT FIND
391      endif
392    7 continue
393      else
394c     ERR HEAVY ATOM1 NOT FIND
395      endif
396    6 continue
397      endif
398    5 continue
399      endif
400    4 continue
401      enddo
402    8 continue
403c
404      do 9 i=1,nseq
405      if(ptseq(i,1).ne.0)nuseq=nuseq+1
406    9 continue
407c
408      write(lfntop,1013)
409 1013 format('qhop')
410      write(lfntop,1014) numseq,larat,numpar,nseq
411 1014 format(2i5,i10,i5)
412      write(lfntop,1025) (shbit(i),i=1,5)
413 1025 format(5i5)
414c
415      do 11 i=1,nseq
416      write(lfntop,1010) i,cseq(i),(ptseq(i,j),j=1,2)
417 1010 format(i5,1x,a10,2i5)
418   11 continue
419c
420      do 21 i=1,numseq
421      write(lfntop,1020) i,cuseq(i),ptarat(i)
422 1020 format(i5,1x,a10,1x,i5)
423   21 continue
424c
425      do 22 i=1,larat
426      write(lfntop,1021) i,(arat(i,j),j=1,4), deq(i),aneq(i)
427 1021 format(i5,4i4,2f8.4)
428   22 continue
429c
430      do 10 i=1,numpar
431      write(lfntop,1015)ptpar(i)
432      write(lfntop,1003)(e12fxy(i,m),m=1,3)
433      write(lfntop,1004)(defe12(i,m),m=1,7)
434      write(lfntop,1005)(zpef(i,m),m=1,6)
435      write(lfntop,1006)(tdsgl(i,m),m=1,6)
436      write(lfntop,1007)(tunnel(i,m),m=1,7)
437      write(lfntop,1008)(tstval(i,m),m=1,5)
438   10 continue
439 1015 format(i10)
440c
441      close(lfnhdb)
442c      close(lfntop)
443      return
444c
445  112 write(*,1030) 'err opening input file'
446 1030 format(a40)
447      return
448      end
449