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