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,
43 double precision, dimension (max_template,maxres) :: rescore
44 double precision, dimension (max_template,maxres) :: rescore2
45 double precision, dimension (max_template,maxres) :: rescore3
46 character*24 tpl_k_rescore
47 c -----------------------------------------------------------------
48 c Reading multiple PDB ref structures and calculation of retraints
49 c not using pre-computed ones stored in files model_ki_{dist,angle}
51 c -----------------------------------------------------------------
54 c Alternative: reading from input
55 call card_concat(controlcard,.true.)
56 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
57 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
58 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
59 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
60 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
61 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
62 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
63 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
64 call readi(controlcard,"IHSET",ihset,1)
65 if (homol_nset.gt.1)then
66 call card_concat(controlcard,.true.)
67 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
68 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
69 c write(iout,*) "iset homology_weight "
71 c write(iout,*) i,waga_homology(i)
74 iset=mod(kolor,homol_nset)+1
79 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
81 cd write (iout,*) "nnt",nnt," nct",nct
93 c Reading HM global scores (prob not required)
96 do k=1,constr_homology
100 c open (4,file="HMscore")
101 c do k=1,constr_homology
102 c read (4,*,end=521) hmscore_tmp
103 c hmscore(k)=hmscore_tmp ! Another transformation can be used
104 c write(*,*) "Model", k, ":", hmscore(k)
115 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
117 if (read_homol_frag) then
118 call read_klapaucjusz
121 do k=1,constr_homology
123 read(inp,'(a)') pdbfile
124 c Next stament causes error upon compilation (?)
125 c if(me.eq.king.or. .not. out1file)
126 c write (iout,'(2a)') 'PDB data will be read from file ',
127 c & pdbfile(:ilen(pdbfile))
128 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
129 & pdbfile(:ilen(pdbfile))
130 open(ipdbin,file=pdbfile,status='old',err=33)
132 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
133 & pdbfile(:ilen(pdbfile))
136 c print *,'Begin reading pdb data'
138 c Files containing res sim or local scores (former containing sigmas)
141 write(kic2,'(bz,i2.2)') k
143 tpl_k_rescore="template"//kic2//".sco"
147 call readpdb_template(k)
162 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
163 & (crefjlee(j,i+nres),j=1,3)
165 write (iout,*) "read_constr_homology: after reading pdb file"
170 c Distance restraints
173 C Copy the coordinates from reference coordinates (?)
177 c write (iout,*) "c(",j,i,") =",c(j,i)
181 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
183 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
184 open (ientin,file=tpl_k_rescore,status='old')
185 if (nnt.gt.1) rescore(k,1)=0.0d0
186 do irec=nnt,nct ! loop for reading res sim
188 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
189 & rescore3_tmp,idomain_tmp
191 idomain(k,i_tmp)=idomain_tmp
192 rescore(k,i_tmp)=rescore_tmp
193 rescore2(k,i_tmp)=rescore2_tmp
194 rescore3(k,i_tmp)=rescore3_tmp
195 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
196 & i_tmp,rescore2_tmp,rescore_tmp,
197 & rescore3_tmp,idomain_tmp
200 read (ientin,*,end=1401) rescore_tmp
202 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
203 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
204 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
209 if (waga_dist.ne.0.0d0) then
217 distal=dsqrt(x12*x12+y12*y12+z12*z12)
218 c write (iout,*) k,i,j,distal,dist2_cut
220 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
221 & .and. distal.le.dist2_cut ) then
227 c write (iout,*) "k",k
228 c write (iout,*) "i",i," j",j," constr_homology",
236 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
238 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
239 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
240 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
242 if (odl(k,ii).le.dist_cut) then
243 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
246 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
247 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
249 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
250 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
254 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
264 c Theta, dihedral and SC retraints
266 if (waga_angle.gt.0.0d0) then
267 c open (ientin,file=tpl_k_sigma_dih,status='old')
268 c do irec=1,maxres-3 ! loop for reading sigma_dih
269 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
270 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
271 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
272 c & sigma_dih(k,i+nnt-1)
277 if (idomain(k,i).eq.0) then
281 dih(k,i)=phiref(i) ! right?
282 c read (ientin,*) sigma_dih(k,i) ! original variant
283 c write (iout,*) "dih(",k,i,") =",dih(k,i)
284 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
285 c & "rescore(",k,i-1,") =",rescore(k,i-1),
286 c & "rescore(",k,i-2,") =",rescore(k,i-2),
287 c & "rescore(",k,i-3,") =",rescore(k,i-3)
289 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
290 & rescore(k,i-2)+rescore(k,i-3))/4.0
291 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
292 c write (iout,*) "Raw sigmas for dihedral angle restraints"
293 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
294 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
295 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
296 c Instead of res sim other local measure of b/b str reliability possible
297 if (sigma_dih(k,i).ne.0)
298 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
299 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
304 if (waga_theta.gt.0.0d0) then
305 c open (ientin,file=tpl_k_sigma_theta,status='old')
306 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
307 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
308 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
309 c & sigma_theta(k,i+nnt-1)
314 do i = nnt+2,nct ! right? without parallel.
315 c do i = i=1,nres ! alternative for bounds acc to readpdb?
316 c do i=ithet_start,ithet_end ! with FG parallel.
317 if (idomain(k,i).eq.0) then
321 thetatpl(k,i)=thetaref(i)
322 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
323 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
324 c & "rescore(",k,i-1,") =",rescore(k,i-1),
325 c & "rescore(",k,i-2,") =",rescore(k,i-2)
326 c read (ientin,*) sigma_theta(k,i) ! 1st variant
327 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
328 & rescore(k,i-2))/3.0
329 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
330 if (sigma_theta(k,i).ne.0)
331 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
333 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
334 c rescore(k,i-2) ! right expression ?
335 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
339 if (waga_d.gt.0.0d0) then
340 c open (ientin,file=tpl_k_sigma_d,status='old')
341 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
342 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
343 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
344 c & sigma_d(k,i+nnt-1)
348 do i = nnt,nct ! right? without parallel.
349 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
350 c do i=loc_start,loc_end ! with FG parallel.
351 if (itype(i).eq.10) cycle
352 if (idomain(k,i).eq.0 ) then
359 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
360 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
361 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
362 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
363 c sigma_d(k,i)=rescore(k,i) ! right expression ?
364 sigma_d(k,i)=rescore3(k,i) ! right expression ?
365 if (sigma_d(k,i).ne.0)
366 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
368 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
369 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
370 c read (ientin,*) sigma_d(k,i) ! 1st variant
375 c remove distance restraints not used in any model from the list
376 c shift data in all arrays
378 if (waga_dist.ne.0.0d0) then
384 if (ii_in_use(ii).eq.0.and.liiflag) then
388 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
389 & .not.liiflag.and.ii.eq.lim_odl) then
390 if (ii.eq.lim_odl) then
396 do ki=iistart,lim_odl-iishift
397 ires_homo(ki)=ires_homo(ki+iishift)
398 jres_homo(ki)=jres_homo(ki+iishift)
399 ii_in_use(ki)=ii_in_use(ki+iishift)
400 do k=1,constr_homology
401 odl(k,ki)=odl(k,ki+iishift)
402 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
403 l_homo(k,ki)=l_homo(k,ki+iishift)
407 lim_odl=lim_odl-iishift
413 endif ! .not. klapaucjusz
415 if (constr_homology.gt.0) call homology_partition
416 if (constr_homology.gt.0) call init_int_table
417 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
418 cd & "lim_xx=",lim_xx
419 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
420 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
424 if (.not.lprn) return
425 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
426 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
427 write (iout,*) "Distance restraints from templates"
429 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
430 & ii,ires_homo(ii),jres_homo(ii),
431 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
432 & ki=1,constr_homology)
434 write (iout,*) "Dihedral angle restraints from templates"
436 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
437 & (rad2deg*dih(ki,i),
438 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
440 write (iout,*) "Virtual-bond angle restraints from templates"
442 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
443 & (rad2deg*thetatpl(ki,i),
444 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
446 write (iout,*) "SC restraints from templates"
448 write(iout,'(i5,100(4f8.2,4x))') i,
449 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
450 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
453 c -----------------------------------------------------------------
456 c----------------------------------------------------------------------
457 subroutine read_klapaucjusz
460 include 'DIMENSIONS.ZSCOPT'
461 include 'DIMENSIONS.FREE'
465 include 'COMMON.SETUP'
466 include 'COMMON.CONTROL'
467 include 'COMMON.HOMOLOGY'
468 include 'COMMON.CHAIN'
469 include 'COMMON.IOUNITS'
471 include 'COMMON.INTERACT'
472 include 'COMMON.NAMES'
473 include 'COMMON.HOMRESTR'
474 character*256 fragfile
475 integer ninclust(maxclust),inclust(max_template,maxclust),
476 & nresclust(maxclust),iresclust(maxres,maxclust)
479 character*24 model_ki_dist, model_ki_angle
480 character*500 controlcard
481 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
482 integer idomain(max_template,maxres)
483 logical lprn /.true./
489 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
490 double precision, dimension (max_template,maxres) :: rescore
491 double precision, dimension (max_template,maxres) :: rescore2
492 character*24 tpl_k_rescore
499 call getenv("FRAGFILE",fragfile)
500 open(ientin,file=fragfile,status="old",err=10)
501 read(ientin,*) constr_homology,nclust
507 do k=1,constr_homology
508 read(ientin,'(a)') pdbfile
509 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
510 & pdbfile(:ilen(pdbfile))
511 open(ipdbin,file=pdbfile,status='old',err=33)
513 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
514 & pdbfile(:ilen(pdbfile))
518 call readpdb_template(k)
521 c chomo(j,i,k)=c(j,i)
531 read(ientin,*) ninclust(i),nresclust(i)
532 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
533 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
539 do ll = 1,ninclust(l)
547 idomain(k,iresclust(i,l)+1) = 1
549 idomain(k,iresclust(i,l)) = 1
553 c Distance restraints
556 C Copy the coordinates from reference coordinates (?)
557 c write (iout,*) "k",k
561 c write (iout,*) "c(",j,i,") =",c(j,i)
565 c write (iout,*) "read_klapaucjusz: calling int_from_cart"
566 call int_from_cart(.true.,.false.)
567 call sc_loc_geom(.false.)
568 c write (iout,*) "en"
573 if (waga_dist.ne.0.0d0) then
581 distal=dsqrt(x12*x12+y12*y12+z12*z12)
582 c write (iout,*) k,i,j,distal,dist2_cut
584 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
585 & .and. distal.le.dist2_cut ) then
591 c write (iout,*) "k",k
592 c write (iout,*) "i",i," j",j," constr_homology",
600 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
602 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
603 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
604 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
606 if (odl(k,ii).le.dist_cut) then
607 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
610 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
611 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
613 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
614 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
618 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
621 c l_homo(k,ii)=.false.
628 c Theta, dihedral and SC retraints
630 if (waga_angle.gt.0.0d0) then
632 if (idomain(k,i).eq.0) then
637 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
638 & rescore(k,i-2)+rescore(k,i-3))/4.0
639 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
640 c & " sigma_dihed",sigma_dih(k,i)
641 if (sigma_dih(k,i).ne.0)
642 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
647 if (waga_theta.gt.0.0d0) then
649 if (idomain(k,i).eq.0) then
650 c sigma_theta(k,i)=0.0
653 thetatpl(k,i)=thetaref(i)
654 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
655 & rescore(k,i-2))/3.0
656 if (sigma_theta(k,i).ne.0)
657 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
661 if (waga_d.gt.0.0d0) then
663 if (itype(i).eq.10) cycle
664 if (idomain(k,i).eq.0 ) then
671 sigma_d(k,i)=rescore(k,i)
672 if (sigma_d(k,i).ne.0)
673 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
674 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
680 c remove distance restraints not used in any model from the list
681 c shift data in all arrays
683 if (waga_dist.ne.0.0d0) then
689 if (ii_in_use(ii).eq.0.and.liiflag) then
693 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
694 & .not.liiflag.and.ii.eq.lim_odl) then
695 if (ii.eq.lim_odl) then
701 do ki=iistart,lim_odl-iishift
702 ires_homo(ki)=ires_homo(ki+iishift)
703 jres_homo(ki)=jres_homo(ki+iishift)
704 ii_in_use(ki)=ii_in_use(ki+iishift)
705 do k=1,constr_homology
706 odl(k,ki)=odl(k,ki+iishift)
707 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
708 l_homo(k,ki)=l_homo(k,ki+iishift)
712 lim_odl=lim_odl-iishift
718 write (iout,*) "ires_homo and jres_homo arrays, lim_odl",lim_odl
720 write (iout,*) i,ires_homo(i),jres_homo(i)
724 10 stop "Error in fragment file"