1 subroutine read_constr_homology
4 include 'DIMENSIONS.ZSCOPT'
5 include 'DIMENSIONS.FREE'
10 include 'COMMON.CONTROL'
11 include 'COMMON.CHAIN'
12 include 'COMMON.IOUNITS'
14 include 'COMMON.INTERACT'
15 include 'COMMON.NAMES'
16 include 'COMMON.HOMRESTR'
17 include 'COMMON.HOMOLOGY'
22 c include 'include_unres/COMMON.VAR'
25 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
27 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
28 c & sigma_odl_temp(maxres,maxres,max_template)
30 character*24 model_ki_dist, model_ki_angle
31 character*500 controlcard
32 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
33 integer idomain(max_template,maxres)
39 c FP - Nov. 2014 Temporary specifications for new vars
41 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
42 double precision, dimension (max_template,maxres) :: rescore
43 double precision, dimension (max_template,maxres) :: rescore2
44 character*24 tpl_k_rescore
45 c -----------------------------------------------------------------
46 c Reading multiple PDB ref structures and calculation of retraints
47 c not using pre-computed ones stored in files model_ki_{dist,angle}
49 c -----------------------------------------------------------------
52 c Alternative: reading from input
53 call card_concat(controlcard,.true.)
54 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
55 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
56 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
57 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
58 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
59 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
60 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
61 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
62 call readi(controlcard,"IHSET",ihset,1)
63 if (homol_nset.gt.1)then
64 call card_concat(controlcard,.true.)
65 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
66 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
67 c write(iout,*) "iset homology_weight "
69 c write(iout,*) i,waga_homology(i)
72 iset=mod(kolor,homol_nset)+1
77 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
79 cd write (iout,*) "nnt",nnt," nct",nct
91 c Reading HM global scores (prob not required)
94 do k=1,constr_homology
98 c open (4,file="HMscore")
99 c do k=1,constr_homology
100 c read (4,*,end=521) hmscore_tmp
101 c hmscore(k)=hmscore_tmp ! Another transformation can be used
102 c write(*,*) "Model", k, ":", hmscore(k)
113 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
115 if (read_homol_frag) then
116 call read_klapaucjusz
119 do k=1,constr_homology
121 read(inp,'(a)') pdbfile
122 c Next stament causes error upon compilation (?)
123 c if(me.eq.king.or. .not. out1file)
124 c write (iout,'(2a)') 'PDB data will be read from file ',
125 c & pdbfile(:ilen(pdbfile))
126 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
127 & pdbfile(:ilen(pdbfile))
128 open(ipdbin,file=pdbfile,status='old',err=33)
130 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
131 & pdbfile(:ilen(pdbfile))
134 c print *,'Begin reading pdb data'
136 c Files containing res sim or local scores (former containing sigmas)
139 write(kic2,'(bz,i2.2)') k
141 tpl_k_rescore="template"//kic2//".sco"
145 call readpdb_template(k)
158 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
159 & (crefjlee(j,i+nres),j=1,3)
161 write (iout,*) "read_constr_homology: after reading pdb file"
166 c Distance restraints
169 C Copy the coordinates from reference coordinates (?)
173 c write (iout,*) "c(",j,i,") =",c(j,i)
177 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
179 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
180 open (ientin,file=tpl_k_rescore,status='old')
181 if (nnt.gt.1) rescore(k,1)=0.0d0
182 do irec=nnt,nct ! loop for reading res sim
184 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
187 idomain(k,i_tmp)=idomain_tmp
188 rescore(k,i_tmp)=rescore_tmp
189 rescore2(k,i_tmp)=rescore2_tmp
190 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
191 & i_tmp,rescore2_tmp,rescore_tmp,
195 read (ientin,*,end=1401) rescore_tmp
197 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
198 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
199 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
204 if (waga_dist.ne.0.0d0) then
212 distal=dsqrt(x12*x12+y12*y12+z12*z12)
213 c write (iout,*) k,i,j,distal,dist2_cut
215 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
216 & .and. distal.le.dist2_cut ) then
222 c write (iout,*) "k",k
223 c write (iout,*) "i",i," j",j," constr_homology",
231 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
233 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
234 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
235 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
237 if (odl(k,ii).le.dist_cut) then
238 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
241 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
242 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
244 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
245 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
249 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
259 c Theta, dihedral and SC retraints
261 if (waga_angle.gt.0.0d0) then
262 c open (ientin,file=tpl_k_sigma_dih,status='old')
263 c do irec=1,maxres-3 ! loop for reading sigma_dih
264 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
265 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
266 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
267 c & sigma_dih(k,i+nnt-1)
272 if (idomain(k,i).eq.0) then
276 dih(k,i)=phiref(i) ! right?
277 c read (ientin,*) sigma_dih(k,i) ! original variant
278 c write (iout,*) "dih(",k,i,") =",dih(k,i)
279 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
280 c & "rescore(",k,i-1,") =",rescore(k,i-1),
281 c & "rescore(",k,i-2,") =",rescore(k,i-2),
282 c & "rescore(",k,i-3,") =",rescore(k,i-3)
284 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
285 & rescore(k,i-2)+rescore(k,i-3))/4.0
286 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
287 c write (iout,*) "Raw sigmas for dihedral angle restraints"
288 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
289 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
290 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
291 c Instead of res sim other local measure of b/b str reliability possible
292 if (sigma_dih(k,i).ne.0)
293 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
294 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
299 if (waga_theta.gt.0.0d0) then
300 c open (ientin,file=tpl_k_sigma_theta,status='old')
301 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
302 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
303 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
304 c & sigma_theta(k,i+nnt-1)
309 do i = nnt+2,nct ! right? without parallel.
310 c do i = i=1,nres ! alternative for bounds acc to readpdb?
311 c do i=ithet_start,ithet_end ! with FG parallel.
312 if (idomain(k,i).eq.0) then
316 thetatpl(k,i)=thetaref(i)
317 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
318 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
319 c & "rescore(",k,i-1,") =",rescore(k,i-1),
320 c & "rescore(",k,i-2,") =",rescore(k,i-2)
321 c read (ientin,*) sigma_theta(k,i) ! 1st variant
322 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
323 & rescore(k,i-2))/3.0
324 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
325 if (sigma_theta(k,i).ne.0)
326 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
328 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
329 c rescore(k,i-2) ! right expression ?
330 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
334 if (waga_d.gt.0.0d0) then
335 c open (ientin,file=tpl_k_sigma_d,status='old')
336 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
337 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
338 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
339 c & sigma_d(k,i+nnt-1)
343 do i = nnt,nct ! right? without parallel.
344 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
345 c do i=loc_start,loc_end ! with FG parallel.
346 if (itype(i).eq.10) cycle
347 if (idomain(k,i).eq.0 ) then
354 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
355 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
356 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
357 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
358 sigma_d(k,i)=rescore(k,i) ! right expression ?
359 if (sigma_d(k,i).ne.0)
360 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
362 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
363 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
364 c read (ientin,*) sigma_d(k,i) ! 1st variant
369 c remove distance restraints not used in any model from the list
370 c shift data in all arrays
372 if (waga_dist.ne.0.0d0) then
378 if (ii_in_use(ii).eq.0.and.liiflag) then
382 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
383 & .not.liiflag.and.ii.eq.lim_odl) then
384 if (ii.eq.lim_odl) then
390 do ki=iistart,lim_odl-iishift
391 ires_homo(ki)=ires_homo(ki+iishift)
392 jres_homo(ki)=jres_homo(ki+iishift)
393 ii_in_use(ki)=ii_in_use(ki+iishift)
394 do k=1,constr_homology
395 odl(k,ki)=odl(k,ki+iishift)
396 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
397 l_homo(k,ki)=l_homo(k,ki+iishift)
401 lim_odl=lim_odl-iishift
407 endif ! .not. klapaucjusz
409 if (constr_homology.gt.0) call homology_partition
410 if (constr_homology.gt.0) call init_int_table
411 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
412 cd & "lim_xx=",lim_xx
413 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
414 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
418 if (.not.lprn) return
419 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
420 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
421 write (iout,*) "Distance restraints from templates"
423 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
424 & ii,ires_homo(ii),jres_homo(ii),
425 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
426 & ki=1,constr_homology)
428 write (iout,*) "Dihedral angle restraints from templates"
430 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
431 & (rad2deg*dih(ki,i),
432 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
434 write (iout,*) "Virtual-bond angle restraints from templates"
436 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
437 & (rad2deg*thetatpl(ki,i),
438 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
440 write (iout,*) "SC restraints from templates"
442 write(iout,'(i5,100(4f8.2,4x))') i,
443 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
444 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
447 c -----------------------------------------------------------------
450 c----------------------------------------------------------------------
451 subroutine read_klapaucjusz
454 include 'DIMENSIONS.ZSCOPT'
455 include 'DIMENSIONS.FREE'
459 include 'COMMON.SETUP'
460 include 'COMMON.CONTROL'
461 include 'COMMON.HOMOLOGY'
462 include 'COMMON.CHAIN'
463 include 'COMMON.IOUNITS'
465 include 'COMMON.INTERACT'
466 include 'COMMON.NAMES'
467 include 'COMMON.HOMRESTR'
468 character*256 fragfile
469 integer ninclust(maxclust),inclust(max_template,maxclust),
470 & nresclust(maxclust),iresclust(maxres,maxclust)
473 character*24 model_ki_dist, model_ki_angle
474 character*500 controlcard
475 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
476 integer idomain(max_template,maxres)
477 logical lprn /.true./
483 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
484 double precision, dimension (max_template,maxres) :: rescore
485 double precision, dimension (max_template,maxres) :: rescore2
486 character*24 tpl_k_rescore
493 call getenv("FRAGFILE",fragfile)
494 open(ientin,file=fragfile,status="old",err=10)
495 read(ientin,*) constr_homology,nclust
501 do k=1,constr_homology
502 read(ientin,'(a)') pdbfile
503 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
504 & pdbfile(:ilen(pdbfile))
505 open(ipdbin,file=pdbfile,status='old',err=33)
507 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
508 & pdbfile(:ilen(pdbfile))
512 call readpdb_template(k)
515 c chomo(j,i,k)=c(j,i)
525 read(ientin,*) ninclust(i),nresclust(i)
526 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
527 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
533 do ll = 1,ninclust(l)
541 idomain(k,iresclust(i,l)+1) = 1
543 idomain(k,iresclust(i,l)) = 1
547 c Distance restraints
550 C Copy the coordinates from reference coordinates (?)
554 c write (iout,*) "c(",j,i,") =",c(j,i)
557 call int_from_cart(.true.,.false.)
558 call sc_loc_geom(.false.)
563 if (waga_dist.ne.0.0d0) then
571 distal=dsqrt(x12*x12+y12*y12+z12*z12)
572 c write (iout,*) k,i,j,distal,dist2_cut
574 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
575 & .and. distal.le.dist2_cut ) then
581 c write (iout,*) "k",k
582 c write (iout,*) "i",i," j",j," constr_homology",
590 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
592 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
593 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
594 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
596 if (odl(k,ii).le.dist_cut) then
597 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
600 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
601 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
603 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
604 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
608 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
611 c l_homo(k,ii)=.false.
618 c Theta, dihedral and SC retraints
620 if (waga_angle.gt.0.0d0) then
622 if (idomain(k,i).eq.0) then
627 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
628 & rescore(k,i-2)+rescore(k,i-3))/4.0
629 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
630 c & " sigma_dihed",sigma_dih(k,i)
631 if (sigma_dih(k,i).ne.0)
632 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
637 if (waga_theta.gt.0.0d0) then
639 if (idomain(k,i).eq.0) then
640 c sigma_theta(k,i)=0.0
643 thetatpl(k,i)=thetaref(i)
644 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
645 & rescore(k,i-2))/3.0
646 if (sigma_theta(k,i).ne.0)
647 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
651 if (waga_d.gt.0.0d0) then
653 if (itype(i).eq.10) cycle
654 if (idomain(k,i).eq.0 ) then
661 sigma_d(k,i)=rescore(k,i)
662 if (sigma_d(k,i).ne.0)
663 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
664 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
670 c remove distance restraints not used in any model from the list
671 c shift data in all arrays
673 if (waga_dist.ne.0.0d0) then
679 if (ii_in_use(ii).eq.0.and.liiflag) then
683 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
684 & .not.liiflag.and.ii.eq.lim_odl) then
685 if (ii.eq.lim_odl) then
691 do ki=iistart,lim_odl-iishift
692 ires_homo(ki)=ires_homo(ki+iishift)
693 jres_homo(ki)=jres_homo(ki+iishift)
694 ii_in_use(ki)=ii_in_use(ki+iishift)
695 do k=1,constr_homology
696 odl(k,ki)=odl(k,ki+iishift)
697 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
698 l_homo(k,ki)=l_homo(k,ki+iishift)
702 lim_odl=lim_odl-iishift
708 write (iout,*) "ires_homo and jres_homo arrays, lim_odl",lim_odl
710 write (iout,*) i,ires_homo(i),jres_homo(i)
714 10 stop "Error in fragment file"