1      subroutine argos_cafe_rdtop(lfntop,filtop,snam)
2c $Id$
3      implicit none
4c
5#include "argos_cafe_common.fh"
6#include "mafdecls.fh"
7#include "global.fh"
8#include "msgids.fh"
9c
10      integer lfntop
11      character*(*) filtop
12c
13      character*80 title(3)
14      character*10 topdat,toptim
15      real*8 releps
16      character*16 stemp,snam(nsatot)
17      integer natyps,nqtyps
18      real*8 rdata(24),vdata(24,4),qdata(24,3),pdata(4,24)
19      integer idata(24)
20      character*6 cdata(24)
21      integer i,ia,iq,j,k,l,m,m2,m3,m4,itemp(24),i_i,i_j,l_i,l_j
22      integer naw,nbw,nhw,ndw,now,ntw,nnw
23      integer nhp
24      integer nas,nbt,nhs,nds,nos,nts,nxs
25      integer numbl,numbi,numbd,iqfr,iqto,nsad,maxpro
26      character*10 string
27c
28      real*8 argos_cafe_charge
29      external argos_cafe_charge
30c
31      if(.not.ma_verify_allocator_stuff()) then
32      call md_abort('ERROR IN MEMORY USE AT START RDTOP',0)
33      endif
34c
35      if(me.eq.0) then
36c
37      open(unit=lfntop,file=filtop(1:index(filtop,' ')-1),
38     + form='formatted',status='old',err=9999)
39c
40      nwc=0
41      nsc=0
42      nmult(1)=0
43      nmult(2)=0
44      nmult(3)=0
45      nmult(4)=0
46      totchg=0.0d0
47      nhop=0
48c
49 1000 format(a)
50      read(lfntop,1001,end=9997,err=9998) title
51 1001 format(a80)
52      read(lfntop,1002,end=9997,err=9998) topdat,toptim,ffield
53 1002 format(12x,3a10)
54      nhp=2
55      if(ffield(1:6).eq.'charmm') nhp=4
56      read(lfntop,1003,end=9997,err=9998) nparms
57      mset=nparms
58      if(lfree) mset=6
59      read(lfntop,1003,end=9997,err=9998) natyps
60      read(lfntop,1003,end=9997,err=9998) nqtyps
61      read(lfntop,1003,end=9997,err=9998) nseq
62 1003 format(i5)
63      read(lfntop,1004,end=9997,err=9998) q14fac,releps
64 1004 format(2f12.6)
65      if(q14fac.lt.small) q14fac=one
66      if(releps.lt.small) releps=one
67      qfac=sqrt(1.389354428d+02/releps)
68      call argos_cafe_inita(natyps,4,nqtyps,4)
69      do 1 i=1,natyps
70      read(lfntop,1005,end=9997,err=9998)
71     + (idata(k),cdata(k),rdata(k),k=1,nparms)
72 1005 format(5x,i5,1x,a6,f12.6)
73      call argos_cafe_para(i,cdata,rdata,idata)
74  101 continue
75    1 continue
76      if(.not.ma_verify_allocator_stuff()) then
77      call md_abort('ERROR IN MEMORY IN RDTOP ATOM TYPES',0)
78      endif
79      do 2 i=1,natyps
80      do 3 j=i,natyps
81      read(lfntop,1006,end=9997,err=9998)
82     + ((vdata(k,l),l=1,4),k=1,nparms)
83 1006 format(10x,4e12.5)
84      call argos_cafe_parv(i,j,vdata)
85    3 continue
86    2 continue
87      if(.not.ma_verify_allocator_stuff()) then
88      call md_abort('ERROR IN MEMORY IN RDTOP VDWAALS',0)
89      endif
90      do 4 i=1,nqtyps
91      read(lfntop,1007,end=9997,err=9998)
92     + ((qdata(k,l),l=1,3),k=1,nparms)
93 1007 format(5x,f12.6,e12.5,f12.6)
94      call argos_cafe_parq(i,qdata)
95    4 continue
96      if(.not.ma_verify_allocator_stuff()) then
97      call md_abort('ERROR IN MEMORY IN RDTOP CHARGE TYPES',0)
98      endif
99      do 4004 i=1,nseq
100      read(lfntop,4008,end=9997,err=9998) string,maxpro
101 4008 format(a10,33x,i5)
102      call argos_cafe_parseq(i,maxpro)
103 4004 continue
104      read(lfntop,1008,end=9997,err=9998) naw,nbw,nhw,ndw,now,ntw,nnw
105 1008 format(5i7,2i10)
106      call argos_cafe_initb(1,1,nbw,2,nhw,nhp,ndw,3,now,3,
107     + ntw,2,nnw,2,naw)
108      read(lfntop,1009,end=9997,err=9998) nas,nbt,nhs,nds,nos,nts,nxs
109 1009 format(5i7,3i10,2i5)
110      call argos_cafe_initb(2,nas,nbt,2,nhs,nhp,nds,3,nos,3,
111     + nts,2,nxs,2,0)
112      do 5 i=1,naw
113      read(lfntop,1010,end=9997,err=9998) wnam(i),ia,iq
114 1010 format(a16,25x,2i5)
115      call argos_cafe_parwiq(i,ia,iq)
116    5 continue
117      idata(4)=0
118      do 6 i=1,nbw
119      read(lfntop,1011,end=9997,err=9998) (idata(j),j=1,3)
120 1011 format(3i7)
121      read(lfntop,2011,end=9997,err=9998)
122     + ((pdata(k,j),k=1,2),j=1,nparms)
123 2011 format(f12.6,e12.5)
124      if(noshak.eq.1.or.noshak.eq.3) idata(3)=0
125      call argos_cafe_parbnd(1,i,idata,pdata)
126      if(idata(3).gt.0) nwc=nwc+1
127    6 continue
128      idata(5)=0
129      do 7 i=1,nhw
130      read(lfntop,1012,end=9997,err=9998) (idata(j),j=1,3)
131 1012 format(3i7)
132      do 70 j=1,nparms
133      read(lfntop,2012,end=9997,err=9998) (pdata(k,j),k=1,nhp)
134 2012 format(2(f10.6,e12.5))
135   70 continue
136      call argos_cafe_parang(1,i,idata,pdata)
137    7 continue
138      idata(6)=0
139      do 8 i=1,ndw
140      read(lfntop,1013,end=9997,err=9998) (idata(j),j=1,4)
141 1013 format(4i7)
142      read(lfntop,2013,end=9997,err=9998)
143     + (itemp(j),(pdata(k,j),k=2,3),j=1,nparms)
144 2013 format(i3,f10.6,e12.5)
145      do 2008 j=1,nparms
146      pdata(1,j)=dble(itemp(j))
147 2008 continue
148      call argos_cafe_pardih(1,i,idata,pdata)
149    8 continue
150      do 9 i=1,now
151      read(lfntop,1014,end=9997,err=9998) (idata(j),j=1,4)
152 1014 format(4i7)
153      read(lfntop,2014,end=9997,err=9998)
154     + ((pdata(k,j),k=2,3),j=1,nparms)
155 2014 format(3x,f10.6,e12.5)
156      do 2009 j=1,nparms
157      pdata(1,j)=0.0d0
158 2009 continue
159      call argos_cafe_parimp(1,i,idata,pdata)
160    9 continue
161      if(.not.ma_push_get(mt_int,max(ntw,nnw,nts,nxs),'i',l_i,i_i))
162     + call md_abort('Failed to allocate temp array',0)
163      if(.not.ma_push_get(mt_int,max(ntw,nnw,nts,nxs),'j',l_j,i_j))
164     + call md_abort('Failed to allocate temp array',0)
165      if(ntw.gt.0) then
166      read(lfntop,1015,end=9997,err=9998) (int_mb(i_i+i-1),i=1,ntw)
167 1015 format(11i7)
168      read(lfntop,1015,end=9997,err=9998) (int_mb(i_j+i-1),i=1,ntw)
169      call argos_cafe_ndxtrd(1,int_mb(i_i),int_mb(i_j),ntw)
170      endif
171      if(nnw.gt.0) then
172      read(lfntop,1016,end=9997,err=9998) (int_mb(i_i+i-1),i=1,nnw)
173 1016 format(11i7)
174      read(lfntop,1016,end=9997,err=9998) (int_mb(i_j+i-1),i=1,nnw)
175      call argos_cafe_ndxxcl(1,int_mb(i_i),int_mb(i_j),nnw)
176      endif
177      read(lfntop,8888) string
178 8888 format(a)
179      read(lfntop,1024,end=9997,err=9998) (ewc(i),i=1,nparms)
180 1024 format(f12.6)
181      if(.not.ma_verify_allocator_stuff()) then
182      call md_abort('ERROR IN MEMORY IN RDTOP SOLVENT',0)
183      endif
184      do 10 i=1,nas
185      read(lfntop,1017,end=9997,err=9998) stemp,j,k,l,m,m2,m3,m4
186 1017 format(a16,i3,2i7,19x,3i5,5x,i5)
187c 1017 format(a16,3i5,15x,3i5,5x,i5)
188      if(i.eq.1) then
189      iqfr=min(m,m2,m3)
190      iqto=max(m,m2,m3)
191      else
192      iqfr=min(iqfr,m,m2,m3)
193      iqto=max(iqto,m,m2,m3)
194      endif
195      write(snam(i),'(a5,a5,i6)') stemp(1:5),stemp(11:15),l
196c      snam(i)=stemp
197      if(j.gt.nsf) nsf=j
198c      if(j.gt.nsm) nsm=j
199      totchg=totchg+argos_cafe_charge(m)
200      if(k.gt.msm)
201     + call md_abort('Topology and Restart are incompatible',0)
202      if(m4.lt.0) nhop=nhop+1
203   10 continue
204      if(scaleq.ge.zero) call argos_cafe_scaleq(iqfr,iqto)
205      idata(4)=0
206      do 11 i=1,nbt
207      read(lfntop,1018,end=9997,err=9998) (idata(j),j=1,3)
208 1018 format(3i7)
209      read(lfntop,2018,end=9997,err=9998)
210     + ((pdata(k,j),k=1,2),j=1,nparms)
211 2018 format(f12.6,e12.5)
212      if(noshak.eq.2.or.noshak.eq.3) idata(3)=0
213      call argos_cafe_parbnd(2,i,idata,pdata)
214      if(idata(3).gt.0) nsc=nsc+1
215   11 continue
216      idata(5)=0
217      do 12 i=1,nhs
218      read(lfntop,1019,end=9997,err=9998) (idata(j),j=1,3)
219 1019 format(3i7)
220      do 120 j=1,nparms
221      read(lfntop,2019,end=9997,err=9998) (pdata(k,j),k=1,nhp)
222 2019 format(2(f10.6,e12.5))
223  120 continue
224      call argos_cafe_parang(2,i,idata,pdata)
225   12 continue
226      idata(6)=0
227      do 13 i=1,nds
228      read(lfntop,1020,end=9997,err=9998) (idata(j),j=1,4)
229 1020 format(4i7)
230      read(lfntop,2020,end=9997,err=9998)
231     + (itemp(j),(pdata(k,j),k=2,3),j=1,nparms)
232 2020 format(i3,f10.6,e12.5)
233      do 3013 j=1,nparms
234      pdata(1,j)=dble(abs(itemp(j)))
235      pdata(3,j)=abs(pdata(3,j))
236 3013 continue
237      call argos_cafe_pardih(2,i,idata,pdata)
238   13 continue
239      do 14 i=1,nos
240      read(lfntop,1021,end=9997,err=9998) (idata(j),j=1,4)
241 1021 format(4i7)
242      read(lfntop,2021,end=9997,err=9998)
243     + ((pdata(k,j),k=2,3),j=1,nparms)
244 2021 format(3x,f10.6,e12.5)
245      do 3014 j=1,nparms
246      pdata(1,j)=dble(abs(itemp(j)))
247      pdata(3,j)=abs(pdata(3,j))
248 3014 continue
249      call argos_cafe_parimp(2,i,idata,pdata)
250   14 continue
251      if(nts.gt.0) then
252      read(lfntop,1022,end=9997,err=9998) (int_mb(i_i+i-1),i=1,nts)
253 1022 format(11i7)
254      read(lfntop,1022,end=9997,err=9998) (int_mb(i_j+i-1),i=1,nts)
255      endif
256      call argos_cafe_ndxtrd(2,int_mb(i_i),int_mb(i_j),nts)
257      if(nxs.gt.0) then
258      read(lfntop,1023,end=9997,err=9998) (int_mb(i_i+i-1),i=1,nxs)
259 1023 format(11i7)
260      read(lfntop,1023,end=9997,err=9998) (int_mb(i_j+i-1),i=1,nxs)
261      endif
262      call argos_cafe_ndxxcl(2,int_mb(i_i),int_mb(i_j),nxs)
263      if(.not.ma_pop_stack(l_j)) call md_abort('Dealloc failed',0)
264      if(.not.ma_pop_stack(l_i)) call md_abort('Dealloc failed',0)
265c
266      close(unit=lfntop)
267c
268      if(.not.ma_verify_allocator_stuff()) then
269      call md_abort('ERROR IN MEMORY IN RDTOP',1)
270      endif
271c
272      endif
273c
274      if(np.gt.1) then
275c
276      numbl=ma_sizeof(mt_log,1,mt_byte)
277      numbi=ma_sizeof(mt_int,1,mt_byte)
278      numbd=ma_sizeof(mt_dbl,1,mt_byte)
279c
280c     broadcast dimensions
281c
282      idata(1)=natyps
283      idata(2)=nqtyps
284      idata(3)=nbw
285      idata(4)=nhw
286      idata(5)=ndw
287      idata(6)=now
288      idata(7)=ntw
289      idata(8)=nnw
290      idata(9)=naw
291      idata(10)=nas
292      idata(11)=nbt
293      idata(12)=nhs
294      idata(13)=nds
295      idata(14)=nos
296      idata(15)=nts
297      idata(16)=nxs
298      idata(17)=mset
299      idata(18)=nparms
300      idata(19)=nhop
301      idata(20)=nhp
302      call ga_brdcst(mcf_01,idata,20*numbi,0)
303      natyps=idata(1)
304      nqtyps=idata(2)
305      nbw=idata(3)
306      nhw=idata(4)
307      ndw=idata(5)
308      now=idata(6)
309      ntw=idata(7)
310      nnw=idata(8)
311      naw=idata(9)
312      nas=idata(10)
313      nbt=idata(11)
314      nhs=idata(12)
315      nds=idata(13)
316      nos=idata(14)
317      nts=idata(15)
318      nxs=idata(16)
319      mset=idata(17)
320      nparms=idata(18)
321      nhop=idata(19)
322      nhp=idata(20)
323c
324c     initialize on nodes other than 0
325c
326      if(me.ne.0) then
327      call argos_cafe_inita(natyps,4,nqtyps,4)
328      call argos_cafe_initb(1,1,nbw,2,nhw,nhp,ndw,3,now,3,
329     + ntw,2,nnw,2,naw)
330      call argos_cafe_initb(2,nas,nbt,2,nhs,nhp,nds,3,nos,3,
331     + nts,2,nxs,2,0)
332      endif
333c
334c     broadcast force field parameters
335c
336c      call ga_brdcst(mcf_02,byte_mb(i_nam),16*mat,0)
337c
338      call ga_brdcst(mcf_02,int_mb(i_typ),nparms*mat*numbi,0)
339      call ga_brdcst(mcf_03,dbl_mb(i_mas),mset*mat*numbd,0)
340      call ga_brdcst(mcf_04,int_mb(i_num),nparms*mat*numbi,0)
341      call ga_brdcst(mcf_05,dbl_mb(i_vdw),mset*mat*mat*map*numbd,0)
342      call ga_brdcst(mcf_06,dbl_mb(i_chg),mset*mqt*mqp*numbd,0)
343      call ga_brdcst(mcf_07,int_mb(i_iwa),mwa*numbi,0)
344      call ga_brdcst(mcf_08,int_mb(i_iwq),mwa*numbi,0)
345      call ga_brdcst(mcf_09,int_mb(i_ibnd(1)),4*mbt(1)*numbi,0)
346      call ga_brdcst(mcf_10,int_mb(i_ibnd(2)),4*mbt(2)*numbi,0)
347      call ga_brdcst(mcf_11,dbl_mb(i_bnd(1)),mset*mbp(1)*mbt(1)*numbd,0)
348      call ga_brdcst(mcf_12,dbl_mb(i_bnd(2)),mset*mbp(2)*mbt(2)*numbd,0)
349      call ga_brdcst(mcf_13,int_mb(i_iang(1)),5*mht(1)*numbi,0)
350      call ga_brdcst(mcf_14,int_mb(i_iang(2)),5*mht(2)*numbi,0)
351      call ga_brdcst(mcf_15,dbl_mb(i_ang(1)),mset*mhp(1)*mht(1)*numbd,0)
352      call ga_brdcst(mcf_16,dbl_mb(i_ang(2)),mset*mhp(2)*mht(2)*numbd,0)
353      call ga_brdcst(mcf_17,int_mb(i_idih(1)),6*mdt(1)*numbi,0)
354      call ga_brdcst(mcf_18,int_mb(i_idih(2)),6*mdt(2)*numbi,0)
355      call ga_brdcst(mcf_19,dbl_mb(i_dih(1)),mset*mdp(1)*mdt(1)*numbd,0)
356      call ga_brdcst(mcf_20,dbl_mb(i_dih(2)),mset*mdp(2)*mdt(2)*numbd,0)
357      call ga_brdcst(mcf_21,int_mb(i_iimp(1)),6*mit(1)*numbi,0)
358      call ga_brdcst(mcf_22,int_mb(i_iimp(2)),6*mit(2)*numbi,0)
359      call ga_brdcst(mcf_23,dbl_mb(i_imp(1)),mset*mip(1)*mit(1)*numbd,0)
360      call ga_brdcst(mcf_24,dbl_mb(i_imp(2)),mset*mip(2)*mit(2)*numbd,0)
361      call ga_brdcst(mcf_25,int_mb(i_itrd(1)),2*(mtt(1)+1)*numbi,0)
362      call ga_brdcst(mcf_26,int_mb(i_itrd(2)),2*(mtt(2)+1)*numbi,0)
363      call ga_brdcst(mcf_27,int_mb(i_ixcl(1)),2*(mxt(1)+1)*numbi,0)
364      call ga_brdcst(mcf_28,int_mb(i_ixcl(2)),2*(mxt(2)+1)*numbi,0)
365      call ga_brdcst(mcf_29,nmult,4*numbi,0)
366      call ga_brdcst(mcf_30,ith,24*numbl,0)
367      call ga_brdcst(mcf_31,ip2,24*numbl,0)
368      call ga_brdcst(mcf_32,ip3,24*numbl,0)
369      call ga_brdcst(mcf_56,qfac,numbd,0)
370c
371      do 15 i=1,nsatot
372      stemp=snam(i)
373      call util_char_ga_brdcst(mcf_66,stemp,0)
374      if(me.ne.0) snam(i)=stemp
375   15 continue
376c      call ga_brdcst(mcf_66,byte_mb(i_snam),16*nsatot,0)
377c
378c      if(lanal) then
379c      call ana_select(byte_mb(i_snam))
380c      call ana_initx()
381c      endif
382c
383      itemp(1)=nwc
384      itemp(2)=nsc
385      itemp(3)=nsf
386      call ga_brdcst(mcf_33,itemp,3*ma_sizeof(mt_int,1,mt_byte),0)
387      nwc=itemp(1)
388      nsc=itemp(2)
389      nsf=itemp(3)
390c
391      endif
392c
393      mmult=2*nmult(1)+3*nmult(2)+4*(nmult(3)+nmult(4))
394      mmuli=nmult(1)+nmult(2)+nmult(3)+nmult(4)
395      if(mmult.gt.0) then
396      if(.not.ma_push_get(mt_int,mmuli,'ixmul',l_ixmul,i_ixmul))
397     + call md_abort('Failed to allocate memory for ixmul',0)
398      if(.not.ma_push_get(mt_int,4*mmult,'imul',l_imul,i_imul))
399     + call md_abort('Failed to allocate memory for imul',0)
400      if(.not.ma_push_get(mt_dbl,3*mmult,'xmul',l_xmul,i_xmul))
401     + call md_abort('Failed to allocate memory for xmul',0)
402      if(.not.ma_push_get(mt_dbl,3*mmult,'fmul',l_fmul,i_fmul))
403     + call md_abort('Failed to allocate memory for fmul',0)
404      call argos_cafe_lstmul(int_mb(i_ixmul),int_mb(i_imul),
405     + mbt(2),int_mb(i_ibnd(2)),mht(2),int_mb(i_iang(2)),
406     + mdt(2),int_mb(i_idih(2)),mit(2),int_mb(i_iimp(2)))
407      endif
408c
409      factmw=zero
410      factms=zero
411      factmp=zero
412      if(noshak.eq.2.or.noshak.eq.3) then
413      nsad=3*nsa-3*ndums
414      else
415      nsad=3*nsa-2*ndums
416      endif
417      if(nwm*(3*nwa-nwc)-3*islow.gt.0)
418     + factmw=two/(rgas*dble(nwm*(3*nwa-nwc)-3*islow))
419      if(nsad-nsc-3*islow.gt.0)
420     + factms=two/(rgas*dble(nsad-nsc-3*islow))
421      if(nwm*(3*nwa-nwc)+nsad-nsc-3*islow.gt.0)
422     + factmp=two/(rgas*dble(nwm*(3*nwa-nwc)+nsad-nsc-3*islow))
423c
424      if(.not.ma_verify_allocator_stuff()) then
425      call md_abort('ERROR IN MEMORY USE AFTER RDTOP',0)
426      endif
427c
428      return
429c
430 9997 continue
431      call md_abort('EOF reading topology file',0)
432 9998 continue
433      call md_abort('Error reading topology file',0)
434 9999 continue
435      call md_abort('Failed to open topology file',0)
436      return
437      end
438