1      logical function pre_mkfrg(irtdb,source,iunit,lfnpar,
2     + lfnout,iconst,
3     + lseq,cseq,mseq,nseq,lsgm,csgm,msgm,nsgm,latm,catm,xatm,qatm,matm,
4     + natm,lato,cato,xato,qato,mato,nato,lbnd,mbnd,lbndt,mbndt,
5     + lang,mang,ldih,mdih,limp,mimp,
6     + maxscf,qscale,llnk,clnk,mlnk,nlnk)
7c
8c $Id$
9c
10c     function to prepare missing fragment definitions
11c
12c     in  : iunit     = dbase logical file number
13c           dir_t     = dbase directory name
14c           lfnout    = output file logical file number
15c           mseq      = dimension of the sequence list
16c
17c     out : lseq(1,*) = segment numbers
18c           lseq(2,*) = number of atoms
19c           lseq(3,*) = index to unique segment
20c           lsgm(1,i) = number of segments of type i
21c           lsgm(2,i) = source: 0=not found; 1=s; 2=x; 3=u; 4=t;
22c           cseq      = segment names
23c
24c
25      implicit none
26c
27#include "mafdecls.fh"
28#include "util.fh"
29#include "pre_common.fh"
30c
31      logical pre_hnames,pre_hybrid
32      external pre_hnames,pre_hybrid
33c
34      integer irtdb
35      character*80 source
36      integer mring,maxscf,iconst
37      parameter (mring=2000)
38      integer lring(6,mring)
39      logical aring(mring)
40c
41      integer pre_atnum
42      logical pre_short,pre_bnd,pre_ang,pre_dih,pre_impctr,pre_bonds
43      external pre_atnum
44      external pre_short,pre_bnd,pre_ang,pre_dih,pre_impctr,pre_bonds
45      logical pre_atype,pre_charge
46      external pre_atype,pre_charge
47c
48      integer iunit,lfnpar,lfnout
49      integer mseq,msgm,matm,mato
50      integer nseq,nsgm,natm,nato
51      integer mbnd,mang,mdih,mimp
52      integer mbndt,nbndt
53      integer lseq(6,mseq),lsgm(3,msgm),latm(5,matm),lato(5,mato)
54      character*5 cseq(2,mseq),csgm(msgm)
55      character*6 catm(3,matm),cato(3,mato)
56      real*8 xatm(3,matm),qatm(matm),xato(3,mato),qato(mato),qscale,qsum
57      integer lbnd(2,mbnd),lang(3,mang),ldih(4,mdih),limp(4,mimp)
58      integer lbndt(2,mbndt)
59      character*255 filnam
60      integer mlnk,nlnk
61      integer llnk(4,mlnk)
62      character*4 clnk(2,mlnk)
63c
64      integer lentmp,length,len_f
65c
66      integer mtyp,i_ltyp,l_ltyp
67      integer nring3,nring4,nring5,nring6
68c
69      integer i,j,natoms,ilo,ihi,iseq,jlo,jhi,ilist
70      integer ibnd,nbnd,nang,ndih,nconn,linkid
71      logical found
72      character*10 date,time
73c
74      integer i_wrk,l_wrk
75c
76      pre_mkfrg=.true.
77c
78      if(util_print('where',print_debug)) then
79      write(lfnout,2000)
80 2000 format(/,'pre_mkfrg ')
81      endif
82c
83      lentmp=index(dirpar(mdirpar),' ')-1
84      ilist=0
85c
86      nbndt=0
87      if(.not.pre_bonds(xatm,latm,matm,1,natm,lbndt,mbndt,nbndt,
88     + llnk,clnk,mlnk,nlnk)) call md_abort('Error in pre_bonds',0)
89c
90      if(iconst.gt.0) then
91      if(.not.pre_hybrid(xatm,latm,matm,natm,lbndt,mbndt,nbndt))
92     + call md_abort('Error in pre_hybrid',0)
93      endif
94c
95      do 1 i=1,nsgm
96      if(lsgm(2,i).eq.0) then
97      ilist=0
98c
99c      if(util_print('where',print_debug)) then
100      write(lfnout,2005) csgm(i)
101 2005 format(/,' Creating fragment for residue ',a5,/)
102c      endif
103c
104c     determine index first segment of this type in sequence list
105c     -----------------------------------------------------------
106c
107      iseq=0
108      do 2 j=1,mseq
109      if(lseq(2,j).eq.i) then
110      iseq=j
111      goto 3
112      endif
113    2 continue
114      pre_mkfrg=.false.
115      return
116    3 continue
117c
118c     find indices ilo : first atom of this segment
119c     ------------ ihi : last atom for this segment
120c                  jlo : first atom
121c                  jhi : last atom
122c
123      ilo=lseq(3,iseq)
124      ihi=lseq(3,iseq+1)-1
125c
126      linkid=3
127      do 55 j=ilo,ihi
128      if(latm(5,j).ge.3) then
129      latm(5,j)=linkid
130      linkid=linkid+1
131      endif
132   55 continue
133c
134      if(ilist.eq.0) then
135c
136      jlo=1
137      jhi=natm
138c      if(csgm(i)(5:5).eq.'N'.or.csgm(i)(5:5).eq.'M') jlo=ilo
139c      if(csgm(i)(5:5).eq.'C'.or.csgm(i)(5:5).eq.'M') jhi=ihi
140c
141c     determine the number of 'real' atoms
142c     ------------------------------------
143c
144      natoms=0
145      do 5 j=jlo,jhi
146      if(latm(2,j).le.0) latm(2,j)=pre_atnum(catm(2,j)(1:2))
147      if(latm(2,j).gt.0) natoms=natoms+1
148    5 continue
149c
150c     complete the list of bonds for the fragment
151c     -------------------------------------------
152c
153      nbnd=0
154      if(.not.ma_push_get(mt_int,12*mato,'wrk',l_wrk,i_wrk))
155     + call md_abort('Unable to allocate wrk array',0)
156      if(.not.pre_bnd(xatm,latm,catm,matm,natm,ilo,ihi,
157     + xato,lato,cato,mato,nato,lbnd,mbnd,nbnd,maxscf,
158     + llnk,clnk,mlnk,nlnk,iconst,mang,lang,mdih,ldih,
159     + int_mb(i_wrk)))
160     + call md_abort('pre_bnd failed',9999)
161      if(.not.ma_pop_stack(l_wrk))
162     + call md_abort('Unable to deallocate wrk array',0)
163c
164c     redetermine hydrogen names
165c     --------------------------
166      if(source(1:3).ne.'pdb') then
167      if(.not.pre_hnames(lato,cato,mato,nato,lbnd,mbnd,nbnd))
168     + call md_abort('pre_hnames failed',9999)
169      endif
170c
171c     complete the list of angles for the fragment
172c     --------------------------------------------
173c
174      nang=0
175      if(.not.pre_ang(lbnd,mbnd,nbnd,lang,mang,nang))
176     + call md_abort('pre_ang failed',9999)
177c
178c     complete the list of torsions for the fragment
179c     --------------------------------------------
180c
181      ndih=0
182      if(.not.pre_dih(lang,mang,nang,ldih,mdih,ndih))
183     + call md_abort('pre_dih failed in pre_mkfrg',9999)
184c
185c     count number of bonds for each atom
186c     -----------------------------------
187c
188      do 13 j=1,nato
189      lato(1,j)=lato(3,j)
190      lato(3,j)=0
191      if(lato(5,j).gt.0) lato(3,j)=1
192   13 continue
193      do 14 ibnd=1,nbnd
194      lato(3,lbnd(1,ibnd))=lato(3,lbnd(1,ibnd))+1
195      lato(3,lbnd(2,ibnd))=lato(3,lbnd(2,ibnd))+1
196   14 continue
197c
198      ilist=1
199      endif
200c
201c     determine center types
202c     ----------------------
203c
204      if(.not.pre_impctr(lfnout,lato,mato,1,1,nato,nato,
205     + lbnd,mbnd,nbnd,lang,mang,nang,ldih,mdih,ndih,
206     + lring,aring,mring,nring3,nring4,nring5,nring6))
207     + call md_abort('pre_impctr failed',9999)
208c
209c     memory allocation for work array ltyp
210c     -------------------------------------
211c
212      mtyp=nato+50
213      if(.not.ma_push_get(mt_int,15*mtyp,'ltyp',l_ltyp,i_ltyp))
214     + call md_abort('Memory allocation failed for ltyp',9999)
215c
216c     determine atom types
217c     --------------------
218c
219      if(.not.pre_atype(lfnout,lfnpar,
220     + lato,cato,mato,lbnd,mbnd,nbnd,
221     + 1,1,nato,nato,int_mb(i_ltyp),mtyp,lring,aring,mring,
222     + nring3,nring4,nring5,nring6,
223     + latm,matm,natm,lbndt,mbndt,nbndt)) then
224      call md_abort('pre_atype failed',9999)
225      endif
226c
227c     memory deallocation
228c     -------------------
229c
230      if(.not.ma_pop_stack(l_ltyp))
231     + call md_abort('Memory deallocation failed for ltyp',9999)
232c
233c     check if all atom types could be determined
234c     -------------------------------------------
235c
236      found=.true.
237      do 16 j=1,nato
238      if(lato(2,j).gt.0) then
239      if(cato(3,j)(1:1).eq.' ') found=.false.
240      endif
241   16 continue
242c
243c     guestimate partial charges
244c     --------------------------
245c
246      if(found) then
247      if(.not.pre_charge(irtdb,lfnout,lfnpar,
248     + source,jlo,ilo,ihi,jhi,
249     + lato,cato,xato,qato,mato,nato,lbnd,mbnd,nbnd,maxscf,qscale))
250     + call md_abort('pre_charge failed',9999)
251      else
252      if(.not.pre_charge(irtdb,lfnout,lfnpar,
253     + source,jlo,ilo,ihi,jhi,
254     + lato,cato,xato,qato,mato,nato,lbnd,mbnd,nbnd,0,qscale))
255     + call md_abort('pre_charge failed',9999)
256      endif
257c
258c     write the fragment file
259c     -----------------------
260c
261      length=index(csgm(i),' ')-1
262      if(length.le.0) length=5
263      filnam=dirpar(mdirpar)(1:lentmp)//csgm(i)(1:length)//'.frg '
264      len_f=index(filnam,' ')-1
265      if(.not.found) filnam=filnam(1:len_f)//'_TMP'
266      len_f=index(filnam,' ')-1
267c
268      if(util_print('where',print_debug)) then
269      write(lfnout,'(a,a)') 'Writing new fragment file ',
270     + filnam(1:len_f)
271      endif
272      open(unit=iunit,file=filnam(1:len_f),form='formatted',
273     + status='unknown',err=9999)
274c
275      call swatch(date,time)
276c
277      write(iunit,2001) date,time
278 2001 format('# This is an automatically generated fragment file',/,
279     + '# Atom types and connectivity were derived from coordinates',/,
280     + '# Atomic partial charges are crude estimates',/,
281     + '# ',2a10,/,'#')
282      length=index(csgm(i),' ')-1
283      if(length.le.0) length=5
284c
285c     write: csgm(i) : fragment name preceded by $
286c            nato    : number of atoms
287c            1       : number of parameter sets / protonation states
288c            1       : default protonation state
289c            0       : number of Z-matrix definitions
290c            csgm(i) : file name
291c
292      write(iunit,2002) csgm(i)(1:length),nato,1,1,0,csgm(i)(1:length)
293 2002 format('$',a,/,4i5,/,a)
294      if(util_print('sequence',print_medium)) then
295      write(lfnout,3002) csgm(i)(1:length)
296 3002 format(/,' Fragment ',t40,a,//,
297     + '  num name  type   link cntr  grp pgrp    charge     polarizab',
298     + /)
299      endif
300      found=.true.
301      qsum=0.0d0
302      do 6 j=1,nato
303      if(lato(2,j).gt.0) then
304c
305      if(util_print('sequence',print_medium)) then
306      write(lfnout,3003) j,cato(2,j),cato(3,j),
307     + lato(5,j),lato(4,j),0,1,1,qato(j),0.0
308 3003 format(i5,1x,a4,2x,a6,5i5,2f12.6)
309      qsum=qsum+qato(j)
310      endif
311      write(iunit,2003) j,cato(2,j),cato(3,j),
312     + lato(5,j),lato(4,j),0,1,1,qato(j),0.0
313 2003 format(i5,a4,2x,a6,5i5,2f12.6)
314      if(cato(3,j)(1:1).eq.' ') found=.false.
315      endif
316    6 continue
317      if(util_print('sequence',print_medium)) then
318      write(lfnout,3009) qsum
319 3009 format(43x,'------------',/,25x,'total charge ',f12.6)
320      endif
321c
322      nconn=0
323      do 15 ibnd=1,nbnd
324      if(lbnd(1,ibnd).ge.1.and.lbnd(1,ibnd).le.nato.and.
325     + lbnd(2,ibnd).ge.1.and.lbnd(2,ibnd).le.nato) then
326      write(iunit,2004) lbnd(1,ibnd),lbnd(2,ibnd)
327 2004 format(2i5)
328      nconn=nconn+1
329      endif
330   15 continue
331c
332      if(nconn.gt.0.and.util_print('sequence',print_medium)) then
333      write(lfnout,3004)
334      do 18 ibnd=1,nbnd
335      if(lbnd(1,ibnd).ge.1.and.lbnd(1,ibnd).le.nato.and.
336     + lbnd(2,ibnd).ge.1.and.lbnd(2,ibnd).le.nato)
337     + write(lfnout,3005) lbnd(1,ibnd),lbnd(2,ibnd)
338   18 continue
339 3004 format(/,' Connectivity',/)
340 3005 format(5x,i3,'-',i3)
341      endif
342c
343      close(unit=iunit)
344      if(util_print('sequence',print_medium)) then
345      write(lfnout,3006)
346 3006 format(' ')
347      endif
348c
349      if(util_print('sequence',print_medium)) then
350      write(lfnout,3007) filnam(1:len_f)
351 3007 format(' Created fragment',t40,a,/)
352      endif
353c
354      if(.not.found.and.util_print('sequence',print_none)) then
355      write(lfnout,3008) csgm(i)(1:length)
356 3008 format(' Unresolved atom types in fragment ',a,/)
357      pre_mkfrg=.false.
358c      call md_abort('Unresolved atom types',0)
359      endif
360c
361      lsgm(2,i)=-6
362c
363      endif
364    1 continue
365c
366c      pre_mkfrg=.true.
367      return
368c
369 9999 continue
370      pre_mkfrg=.false.
371      return
372      end
373