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"
169 c Distance restraints
172 C Copy the coordinates from reference coordinates (?)
176 c write (iout,*) "c(",j,i,") =",c(j,i)
180 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
182 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
183 open (ientin,file=tpl_k_rescore,status='old')
184 if (nnt.gt.1) rescore(k,1)=0.0d0
185 do irec=nnt,nct ! loop for reading res sim
187 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
188 & rescore3_tmp,idomain_tmp
190 idomain(k,i_tmp)=idomain_tmp
191 rescore(k,i_tmp)=rescore_tmp
192 rescore2(k,i_tmp)=rescore2_tmp
193 rescore3(k,i_tmp)=rescore3_tmp
194 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
195 & i_tmp,rescore2_tmp,rescore_tmp,
196 & rescore3_tmp,idomain_tmp
199 read (ientin,*,end=1401) rescore_tmp
201 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
202 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
203 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
208 if (waga_dist.ne.0.0d0) then
216 distal=dsqrt(x12*x12+y12*y12+z12*z12)
217 c write (iout,*) k,i,j,distal,dist2_cut
219 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
220 & .and. distal.le.dist2_cut ) then
226 c write (iout,*) "k",k
227 c write (iout,*) "i",i," j",j," constr_homology",
235 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
237 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
238 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
239 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
241 if (odl(k,ii).le.dist_cut) then
242 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
245 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
246 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
248 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
249 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
253 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
263 c Theta, dihedral and SC retraints
265 if (waga_angle.gt.0.0d0) then
266 c open (ientin,file=tpl_k_sigma_dih,status='old')
267 c do irec=1,maxres-3 ! loop for reading sigma_dih
268 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
269 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
270 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
271 c & sigma_dih(k,i+nnt-1)
276 if (idomain(k,i).eq.0) then
280 dih(k,i)=phiref(i) ! right?
281 c read (ientin,*) sigma_dih(k,i) ! original variant
282 c write (iout,*) "dih(",k,i,") =",dih(k,i)
283 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
284 c & "rescore(",k,i-1,") =",rescore(k,i-1),
285 c & "rescore(",k,i-2,") =",rescore(k,i-2),
286 c & "rescore(",k,i-3,") =",rescore(k,i-3)
288 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
289 & rescore(k,i-2)+rescore(k,i-3))/4.0
290 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
291 c write (iout,*) "Raw sigmas for dihedral angle restraints"
292 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
293 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
294 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
295 c Instead of res sim other local measure of b/b str reliability possible
296 if (sigma_dih(k,i).ne.0)
297 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
298 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
303 if (waga_theta.gt.0.0d0) then
304 c open (ientin,file=tpl_k_sigma_theta,status='old')
305 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
306 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
307 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
308 c & sigma_theta(k,i+nnt-1)
313 do i = nnt+2,nct ! right? without parallel.
314 c do i = i=1,nres ! alternative for bounds acc to readpdb?
315 c do i=ithet_start,ithet_end ! with FG parallel.
316 if (idomain(k,i).eq.0) then
320 thetatpl(k,i)=thetaref(i)
321 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
322 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
323 c & "rescore(",k,i-1,") =",rescore(k,i-1),
324 c & "rescore(",k,i-2,") =",rescore(k,i-2)
325 c read (ientin,*) sigma_theta(k,i) ! 1st variant
326 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
327 & rescore(k,i-2))/3.0
328 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
329 if (sigma_theta(k,i).ne.0)
330 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
332 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
333 c rescore(k,i-2) ! right expression ?
334 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
338 if (waga_d.gt.0.0d0) then
339 c open (ientin,file=tpl_k_sigma_d,status='old')
340 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
341 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
342 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
343 c & sigma_d(k,i+nnt-1)
347 do i = nnt,nct ! right? without parallel.
348 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
349 c do i=loc_start,loc_end ! with FG parallel.
350 if (itype(i).eq.10) cycle
351 if (idomain(k,i).eq.0 ) then
358 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
359 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
360 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
361 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
362 c sigma_d(k,i)=rescore(k,i) ! right expression ?
363 sigma_d(k,i)=rescore3(k,i) ! right expression ?
364 if (sigma_d(k,i).ne.0)
365 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
367 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
368 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
369 c read (ientin,*) sigma_d(k,i) ! 1st variant
374 c remove distance restraints not used in any model from the list
375 c shift data in all arrays
377 if (waga_dist.ne.0.0d0) then
383 if (ii_in_use(ii).eq.0.and.liiflag) then
387 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
388 & .not.liiflag.and.ii.eq.lim_odl) then
389 if (ii.eq.lim_odl) then
395 do ki=iistart,lim_odl-iishift
396 ires_homo(ki)=ires_homo(ki+iishift)
397 jres_homo(ki)=jres_homo(ki+iishift)
398 ii_in_use(ki)=ii_in_use(ki+iishift)
399 do k=1,constr_homology
400 odl(k,ki)=odl(k,ki+iishift)
401 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
402 l_homo(k,ki)=l_homo(k,ki+iishift)
406 lim_odl=lim_odl-iishift
412 endif ! .not. klapaucjusz
414 if (constr_homology.gt.0) call homology_partition
415 if (constr_homology.gt.0) call init_int_table
416 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
417 cd & "lim_xx=",lim_xx
418 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
419 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
423 if (.not.lprn) return
424 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
425 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
426 write (iout,*) "Distance restraints from templates"
428 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
429 & ii,ires_homo(ii),jres_homo(ii),
430 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
431 & ki=1,constr_homology)
433 write (iout,*) "Dihedral angle restraints from templates"
435 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
436 & (rad2deg*dih(ki,i),
437 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
439 write (iout,*) "Virtual-bond angle restraints from templates"
441 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
442 & (rad2deg*thetatpl(ki,i),
443 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
445 write (iout,*) "SC restraints from templates"
447 write(iout,'(i5,100(4f8.2,4x))') i,
448 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
449 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
452 c -----------------------------------------------------------------
455 c----------------------------------------------------------------------
456 subroutine read_klapaucjusz
459 include 'DIMENSIONS.ZSCOPT'
460 include 'DIMENSIONS.FREE'
464 include 'COMMON.SETUP'
465 include 'COMMON.CONTROL'
466 include 'COMMON.HOMOLOGY'
467 include 'COMMON.CHAIN'
468 include 'COMMON.IOUNITS'
470 include 'COMMON.INTERACT'
471 include 'COMMON.NAMES'
472 include 'COMMON.HOMRESTR'
473 character*256 fragfile
474 integer ninclust(maxclust),inclust(max_template,maxclust),
475 & nresclust(maxclust),iresclust(maxres,maxclust)
478 character*24 model_ki_dist, model_ki_angle
479 character*500 controlcard
480 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
481 integer idomain(max_template,maxres)
482 logical lprn /.true./
488 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
489 double precision, dimension (max_template,maxres) :: rescore
490 double precision, dimension (max_template,maxres) :: rescore2
491 character*24 tpl_k_rescore
498 call getenv("FRAGFILE",fragfile)
499 open(ientin,file=fragfile,status="old",err=10)
500 read(ientin,*) constr_homology,nclust
506 do k=1,constr_homology
507 read(ientin,'(a)') pdbfile
508 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
509 & pdbfile(:ilen(pdbfile))
510 open(ipdbin,file=pdbfile,status='old',err=33)
512 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
513 & pdbfile(:ilen(pdbfile))
517 call readpdb_template(k)
520 c chomo(j,i,k)=c(j,i)
530 read(ientin,*) ninclust(i),nresclust(i)
531 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
532 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
538 do ll = 1,ninclust(l)
546 idomain(k,iresclust(i,l)+1) = 1
548 idomain(k,iresclust(i,l)) = 1
552 c Distance restraints
555 C Copy the coordinates from reference coordinates (?)
559 c write (iout,*) "c(",j,i,") =",c(j,i)
562 call int_from_cart(.true.,.false.)
563 call sc_loc_geom(.false.)
568 if (waga_dist.ne.0.0d0) then
576 distal=dsqrt(x12*x12+y12*y12+z12*z12)
577 c write (iout,*) k,i,j,distal,dist2_cut
579 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
580 & .and. distal.le.dist2_cut ) then
586 c write (iout,*) "k",k
587 c write (iout,*) "i",i," j",j," constr_homology",
595 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
597 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
598 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
599 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
601 if (odl(k,ii).le.dist_cut) then
602 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
605 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
606 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
608 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
609 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
613 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
616 c l_homo(k,ii)=.false.
623 c Theta, dihedral and SC retraints
625 if (waga_angle.gt.0.0d0) then
627 if (idomain(k,i).eq.0) then
632 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
633 & rescore(k,i-2)+rescore(k,i-3))/4.0
634 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
635 c & " sigma_dihed",sigma_dih(k,i)
636 if (sigma_dih(k,i).ne.0)
637 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
642 if (waga_theta.gt.0.0d0) then
644 if (idomain(k,i).eq.0) then
645 c sigma_theta(k,i)=0.0
648 thetatpl(k,i)=thetaref(i)
649 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
650 & rescore(k,i-2))/3.0
651 if (sigma_theta(k,i).ne.0)
652 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
656 if (waga_d.gt.0.0d0) then
658 if (itype(i).eq.10) cycle
659 if (idomain(k,i).eq.0 ) then
666 sigma_d(k,i)=rescore(k,i)
667 if (sigma_d(k,i).ne.0)
668 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
669 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
675 c remove distance restraints not used in any model from the list
676 c shift data in all arrays
678 if (waga_dist.ne.0.0d0) then
684 if (ii_in_use(ii).eq.0.and.liiflag) then
688 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
689 & .not.liiflag.and.ii.eq.lim_odl) then
690 if (ii.eq.lim_odl) then
696 do ki=iistart,lim_odl-iishift
697 ires_homo(ki)=ires_homo(ki+iishift)
698 jres_homo(ki)=jres_homo(ki+iishift)
699 ii_in_use(ki)=ii_in_use(ki+iishift)
700 do k=1,constr_homology
701 odl(k,ki)=odl(k,ki+iishift)
702 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
703 l_homo(k,ki)=l_homo(k,ki+iishift)
707 lim_odl=lim_odl-iishift
713 write (iout,*) "ires_homo and jres_homo arrays, lim_odl",lim_odl
715 write (iout,*) i,ires_homo(i),jres_homo(i)
719 10 stop "Error in fragment file"