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)
160 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
161 & (crefjlee(j,i+nres),j=1,3)
163 write (iout,*) "read_constr_homology: after reading pdb file"
168 c Distance restraints
171 C Copy the coordinates from reference coordinates (?)
175 c write (iout,*) "c(",j,i,") =",c(j,i)
179 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
181 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
182 open (ientin,file=tpl_k_rescore,status='old')
183 if (nnt.gt.1) rescore(k,1)=0.0d0
184 do irec=nnt,nct ! loop for reading res sim
186 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
187 & rescore3_tmp,idomain_tmp
189 idomain(k,i_tmp)=idomain_tmp
190 rescore(k,i_tmp)=rescore_tmp
191 rescore2(k,i_tmp)=rescore2_tmp
192 rescore3(k,i_tmp)=rescore3_tmp
193 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
194 & i_tmp,rescore2_tmp,rescore_tmp,
195 & rescore3_tmp,idomain_tmp
198 read (ientin,*,end=1401) rescore_tmp
200 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
201 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
202 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
207 if (waga_dist.ne.0.0d0) then
215 distal=dsqrt(x12*x12+y12*y12+z12*z12)
216 c write (iout,*) k,i,j,distal,dist2_cut
218 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
219 & .and. distal.le.dist2_cut ) then
225 c write (iout,*) "k",k
226 c write (iout,*) "i",i," j",j," constr_homology",
234 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
236 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
237 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
238 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
240 if (odl(k,ii).le.dist_cut) then
241 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
244 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
245 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
247 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
248 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
252 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
262 c Theta, dihedral and SC retraints
264 if (waga_angle.gt.0.0d0) then
265 c open (ientin,file=tpl_k_sigma_dih,status='old')
266 c do irec=1,maxres-3 ! loop for reading sigma_dih
267 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
268 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
269 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
270 c & sigma_dih(k,i+nnt-1)
275 if (idomain(k,i).eq.0) then
279 dih(k,i)=phiref(i) ! right?
280 c read (ientin,*) sigma_dih(k,i) ! original variant
281 c write (iout,*) "dih(",k,i,") =",dih(k,i)
282 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
283 c & "rescore(",k,i-1,") =",rescore(k,i-1),
284 c & "rescore(",k,i-2,") =",rescore(k,i-2),
285 c & "rescore(",k,i-3,") =",rescore(k,i-3)
287 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
288 & rescore(k,i-2)+rescore(k,i-3))/4.0
289 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
290 c write (iout,*) "Raw sigmas for dihedral angle restraints"
291 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
292 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
293 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
294 c Instead of res sim other local measure of b/b str reliability possible
295 if (sigma_dih(k,i).ne.0)
296 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
297 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
302 if (waga_theta.gt.0.0d0) then
303 c open (ientin,file=tpl_k_sigma_theta,status='old')
304 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
305 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
306 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
307 c & sigma_theta(k,i+nnt-1)
312 do i = nnt+2,nct ! right? without parallel.
313 c do i = i=1,nres ! alternative for bounds acc to readpdb?
314 c do i=ithet_start,ithet_end ! with FG parallel.
315 if (idomain(k,i).eq.0) then
319 thetatpl(k,i)=thetaref(i)
320 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
321 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
322 c & "rescore(",k,i-1,") =",rescore(k,i-1),
323 c & "rescore(",k,i-2,") =",rescore(k,i-2)
324 c read (ientin,*) sigma_theta(k,i) ! 1st variant
325 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
326 & rescore(k,i-2))/3.0
327 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
328 if (sigma_theta(k,i).ne.0)
329 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
331 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
332 c rescore(k,i-2) ! right expression ?
333 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
337 if (waga_d.gt.0.0d0) then
338 c open (ientin,file=tpl_k_sigma_d,status='old')
339 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
340 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
341 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
342 c & sigma_d(k,i+nnt-1)
346 do i = nnt,nct ! right? without parallel.
347 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
348 c do i=loc_start,loc_end ! with FG parallel.
349 if (itype(i).eq.10) cycle
350 if (idomain(k,i).eq.0 ) then
357 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
358 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
359 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
360 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
361 sigma_d(k,i)=rescore3(k,i) ! right expression ?
362 if (sigma_d(k,i).ne.0)
363 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
365 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
366 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
367 c read (ientin,*) sigma_d(k,i) ! 1st variant
372 c remove distance restraints not used in any model from the list
373 c shift data in all arrays
375 if (waga_dist.ne.0.0d0) then
381 if (ii_in_use(ii).eq.0.and.liiflag) then
385 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
386 & .not.liiflag.and.ii.eq.lim_odl) then
387 if (ii.eq.lim_odl) then
393 do ki=iistart,lim_odl-iishift
394 ires_homo(ki)=ires_homo(ki+iishift)
395 jres_homo(ki)=jres_homo(ki+iishift)
396 ii_in_use(ki)=ii_in_use(ki+iishift)
397 do k=1,constr_homology
398 odl(k,ki)=odl(k,ki+iishift)
399 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
400 l_homo(k,ki)=l_homo(k,ki+iishift)
404 lim_odl=lim_odl-iishift
410 endif ! .not. klapaucjusz
412 if (constr_homology.gt.0) call homology_partition
413 if (constr_homology.gt.0) call init_int_table
414 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
415 cd & "lim_xx=",lim_xx
416 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
417 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
421 if (.not.lprn) return
422 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
423 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
424 write (iout,*) "Distance restraints from templates"
426 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
427 & ii,ires_homo(ii),jres_homo(ii),
428 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
429 & ki=1,constr_homology)
431 write (iout,*) "Dihedral angle restraints from templates"
433 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
434 & (rad2deg*dih(ki,i),
435 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
437 write (iout,*) "Virtual-bond angle restraints from templates"
439 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
440 & (rad2deg*thetatpl(ki,i),
441 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
443 write (iout,*) "SC restraints from templates"
445 write(iout,'(i5,100(4f8.2,4x))') i,
446 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
447 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
450 c -----------------------------------------------------------------
453 c----------------------------------------------------------------------
454 subroutine read_klapaucjusz
457 include 'DIMENSIONS.ZSCOPT'
458 include 'DIMENSIONS.FREE'
462 include 'COMMON.SETUP'
463 include 'COMMON.CONTROL'
464 include 'COMMON.HOMOLOGY'
465 include 'COMMON.CHAIN'
466 include 'COMMON.IOUNITS'
468 include 'COMMON.INTERACT'
469 include 'COMMON.NAMES'
470 include 'COMMON.HOMRESTR'
471 character*256 fragfile
472 integer ninclust(maxclust),inclust(max_template,maxclust),
473 & nresclust(maxclust),iresclust(maxres,maxclust)
476 character*24 model_ki_dist, model_ki_angle
477 character*500 controlcard
478 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
479 integer idomain(max_template,maxres)
480 logical lprn /.true./
486 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
487 double precision, dimension (max_template,maxres) :: rescore
488 double precision, dimension (max_template,maxres) :: rescore2
489 character*24 tpl_k_rescore
496 call getenv("FRAGFILE",fragfile)
497 open(ientin,file=fragfile,status="old",err=10)
498 read(ientin,*) constr_homology,nclust
504 do k=1,constr_homology
505 read(ientin,'(a)') pdbfile
506 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
507 & pdbfile(:ilen(pdbfile))
508 open(ipdbin,file=pdbfile,status='old',err=33)
510 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
511 & pdbfile(:ilen(pdbfile))
515 call readpdb_template(k)
518 c chomo(j,i,k)=c(j,i)
528 read(ientin,*) ninclust(i),nresclust(i)
529 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
530 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
536 do ll = 1,ninclust(l)
544 idomain(k,iresclust(i,l)+1) = 1
546 idomain(k,iresclust(i,l)) = 1
550 c Distance restraints
553 C Copy the coordinates from reference coordinates (?)
557 c write (iout,*) "c(",j,i,") =",c(j,i)
560 call int_from_cart(.true.,.false.)
561 call sc_loc_geom(.false.)
566 if (waga_dist.ne.0.0d0) then
574 distal=dsqrt(x12*x12+y12*y12+z12*z12)
575 c write (iout,*) k,i,j,distal,dist2_cut
577 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
578 & .and. distal.le.dist2_cut ) then
584 c write (iout,*) "k",k
585 c write (iout,*) "i",i," j",j," constr_homology",
593 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
595 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
596 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
597 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
599 if (odl(k,ii).le.dist_cut) then
600 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
603 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
604 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
606 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
607 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
611 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
614 c l_homo(k,ii)=.false.
621 c Theta, dihedral and SC retraints
623 if (waga_angle.gt.0.0d0) then
625 if (idomain(k,i).eq.0) then
630 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
631 & rescore(k,i-2)+rescore(k,i-3))/4.0
632 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
633 c & " sigma_dihed",sigma_dih(k,i)
634 if (sigma_dih(k,i).ne.0)
635 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
640 if (waga_theta.gt.0.0d0) then
642 if (idomain(k,i).eq.0) then
643 c sigma_theta(k,i)=0.0
646 thetatpl(k,i)=thetaref(i)
647 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
648 & rescore(k,i-2))/3.0
649 if (sigma_theta(k,i).ne.0)
650 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
654 if (waga_d.gt.0.0d0) then
656 if (itype(i).eq.10) cycle
657 if (idomain(k,i).eq.0 ) then
664 sigma_d(k,i)=rescore(k,i)
665 if (sigma_d(k,i).ne.0)
666 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
667 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
673 c remove distance restraints not used in any model from the list
674 c shift data in all arrays
676 if (waga_dist.ne.0.0d0) then
682 if (ii_in_use(ii).eq.0.and.liiflag) then
686 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
687 & .not.liiflag.and.ii.eq.lim_odl) then
688 if (ii.eq.lim_odl) then
694 do ki=iistart,lim_odl-iishift
695 ires_homo(ki)=ires_homo(ki+iishift)
696 jres_homo(ki)=jres_homo(ki+iishift)
697 ii_in_use(ki)=ii_in_use(ki+iishift)
698 do k=1,constr_homology
699 odl(k,ki)=odl(k,ki+iishift)
700 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
701 l_homo(k,ki)=l_homo(k,ki+iishift)
705 lim_odl=lim_odl-iishift
711 write (iout,*) "ires_homo and jres_homo arrays, lim_odl",lim_odl
713 write (iout,*) i,ires_homo(i),jres_homo(i)
717 10 stop "Error in fragment file"