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 c sigma_d(k,i)=rescore(k,i) ! right expression ?
362 sigma_d(k,i)=rescore3(k,i) ! right expression ?
363 if (sigma_d(k,i).ne.0)
364 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
366 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
367 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
368 c read (ientin,*) sigma_d(k,i) ! 1st variant
373 c remove distance restraints not used in any model from the list
374 c shift data in all arrays
376 if (waga_dist.ne.0.0d0) then
382 if (ii_in_use(ii).eq.0.and.liiflag) then
386 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
387 & .not.liiflag.and.ii.eq.lim_odl) then
388 if (ii.eq.lim_odl) then
394 do ki=iistart,lim_odl-iishift
395 ires_homo(ki)=ires_homo(ki+iishift)
396 jres_homo(ki)=jres_homo(ki+iishift)
397 ii_in_use(ki)=ii_in_use(ki+iishift)
398 do k=1,constr_homology
399 odl(k,ki)=odl(k,ki+iishift)
400 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
401 l_homo(k,ki)=l_homo(k,ki+iishift)
405 lim_odl=lim_odl-iishift
411 endif ! .not. klapaucjusz
413 if (constr_homology.gt.0) call homology_partition
414 if (constr_homology.gt.0) call init_int_table
415 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
416 cd & "lim_xx=",lim_xx
417 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
418 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
422 if (.not.lprn) return
423 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
424 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
425 write (iout,*) "Distance restraints from templates"
427 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
428 & ii,ires_homo(ii),jres_homo(ii),
429 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
430 & ki=1,constr_homology)
432 write (iout,*) "Dihedral angle restraints from templates"
434 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
435 & (rad2deg*dih(ki,i),
436 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
438 write (iout,*) "Virtual-bond angle restraints from templates"
440 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
441 & (rad2deg*thetatpl(ki,i),
442 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
444 write (iout,*) "SC restraints from templates"
446 write(iout,'(i5,100(4f8.2,4x))') i,
447 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
448 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
451 c -----------------------------------------------------------------
454 c----------------------------------------------------------------------
455 subroutine read_klapaucjusz
458 include 'DIMENSIONS.ZSCOPT'
459 include 'DIMENSIONS.FREE'
463 include 'COMMON.SETUP'
464 include 'COMMON.CONTROL'
465 include 'COMMON.HOMOLOGY'
466 include 'COMMON.CHAIN'
467 include 'COMMON.IOUNITS'
469 include 'COMMON.INTERACT'
470 include 'COMMON.NAMES'
471 include 'COMMON.HOMRESTR'
472 character*256 fragfile
473 integer ninclust(maxclust),inclust(max_template,maxclust),
474 & nresclust(maxclust),iresclust(maxres,maxclust)
477 character*24 model_ki_dist, model_ki_angle
478 character*500 controlcard
479 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
480 integer idomain(max_template,maxres)
481 logical lprn /.true./
487 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
488 double precision, dimension (max_template,maxres) :: rescore
489 double precision, dimension (max_template,maxres) :: rescore2
490 character*24 tpl_k_rescore
497 call getenv("FRAGFILE",fragfile)
498 open(ientin,file=fragfile,status="old",err=10)
499 read(ientin,*) constr_homology,nclust
505 do k=1,constr_homology
506 read(ientin,'(a)') pdbfile
507 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
508 & pdbfile(:ilen(pdbfile))
509 open(ipdbin,file=pdbfile,status='old',err=33)
511 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
512 & pdbfile(:ilen(pdbfile))
516 call readpdb_template(k)
519 c chomo(j,i,k)=c(j,i)
529 read(ientin,*) ninclust(i),nresclust(i)
530 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
531 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
537 do ll = 1,ninclust(l)
545 idomain(k,iresclust(i,l)+1) = 1
547 idomain(k,iresclust(i,l)) = 1
551 c Distance restraints
554 C Copy the coordinates from reference coordinates (?)
558 c write (iout,*) "c(",j,i,") =",c(j,i)
561 call int_from_cart(.true.,.false.)
562 call sc_loc_geom(.false.)
567 if (waga_dist.ne.0.0d0) then
575 distal=dsqrt(x12*x12+y12*y12+z12*z12)
576 c write (iout,*) k,i,j,distal,dist2_cut
578 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
579 & .and. distal.le.dist2_cut ) then
585 c write (iout,*) "k",k
586 c write (iout,*) "i",i," j",j," constr_homology",
594 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
596 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
597 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
598 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
600 if (odl(k,ii).le.dist_cut) then
601 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
604 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
605 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
607 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
608 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
612 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
615 c l_homo(k,ii)=.false.
622 c Theta, dihedral and SC retraints
624 if (waga_angle.gt.0.0d0) then
626 if (idomain(k,i).eq.0) then
631 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
632 & rescore(k,i-2)+rescore(k,i-3))/4.0
633 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
634 c & " sigma_dihed",sigma_dih(k,i)
635 if (sigma_dih(k,i).ne.0)
636 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
641 if (waga_theta.gt.0.0d0) then
643 if (idomain(k,i).eq.0) then
644 c sigma_theta(k,i)=0.0
647 thetatpl(k,i)=thetaref(i)
648 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
649 & rescore(k,i-2))/3.0
650 if (sigma_theta(k,i).ne.0)
651 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
655 if (waga_d.gt.0.0d0) then
657 if (itype(i).eq.10) cycle
658 if (idomain(k,i).eq.0 ) then
665 sigma_d(k,i)=rescore(k,i)
666 if (sigma_d(k,i).ne.0)
667 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
668 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
674 c remove distance restraints not used in any model from the list
675 c shift data in all arrays
677 if (waga_dist.ne.0.0d0) then
683 if (ii_in_use(ii).eq.0.and.liiflag) then
687 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
688 & .not.liiflag.and.ii.eq.lim_odl) then
689 if (ii.eq.lim_odl) then
695 do ki=iistart,lim_odl-iishift
696 ires_homo(ki)=ires_homo(ki+iishift)
697 jres_homo(ki)=jres_homo(ki+iishift)
698 ii_in_use(ki)=ii_in_use(ki+iishift)
699 do k=1,constr_homology
700 odl(k,ki)=odl(k,ki+iishift)
701 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
702 l_homo(k,ki)=l_homo(k,ki+iishift)
706 lim_odl=lim_odl-iishift
712 write (iout,*) "ires_homo and jres_homo arrays, lim_odl",lim_odl
714 write (iout,*) i,ires_homo(i),jres_homo(i)
718 10 stop "Error in fragment file"