1 subroutine read_constr_homology
5 include 'DIMENSIONS.ZSCOPT'
6 include 'DIMENSIONS.FREE'
10 include 'COMMON.SETUP'
11 include 'COMMON.CONTROL'
12 include 'COMMON.CHAIN'
13 include 'COMMON.IOUNITS'
15 include 'COMMON.INTERACT'
16 include 'COMMON.NAMES'
17 include 'COMMON.HOMRESTR'
18 include 'COMMON.HOMOLOGY'
23 c include 'include_unres/COMMON.VAR'
26 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
28 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
29 c & sigma_odl_temp(maxres,maxres,max_template)
31 character*24 model_ki_dist, model_ki_angle
32 character*500 controlcard
33 integer ki,i,ii,ik,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,
34 & lim_theta,lim_xx,irec,iistart,iishift,i10,i01
35 double precision distal
36 integer idomain(max_template,maxres)
43 c FP - Nov. 2014 Temporary specifications for new vars
45 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
47 double precision, dimension (max_template,maxres) :: rescore
48 double precision, dimension (max_template,maxres) :: rescore2
49 double precision, dimension (max_template,maxres) :: rescore3
50 character*24 tpl_k_rescore
51 c -----------------------------------------------------------------
52 c Reading multiple PDB ref structures and calculation of retraints
53 c not using pre-computed ones stored in files model_ki_{dist,angle}
55 c -----------------------------------------------------------------
58 c Alternative: reading from input
59 call card_concat(controlcard,.true.)
60 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
61 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
62 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
63 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
64 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
65 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
66 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
67 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
68 call readi(controlcard,"IHSET",ihset,1)
69 if (homol_nset.gt.1)then
70 call card_concat(controlcard,.true.)
71 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
72 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
73 c write(iout,*) "iset homology_weight "
75 c write(iout,*) i,waga_homology(i)
78 iset=mod(kolor,homol_nset)+1
83 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
85 cd write (iout,*) "nnt",nnt," nct",nct
97 c Reading HM global scores (prob not required)
100 do k=1,constr_homology
104 c open (4,file="HMscore")
105 c do k=1,constr_homology
106 c read (4,*,end=521) hmscore_tmp
107 c hmscore(k)=hmscore_tmp ! Another transformation can be used
108 c write(*,*) "Model", k, ":", hmscore(k)
119 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
121 if (read_homol_frag) then
122 call read_klapaucjusz
125 do k=1,constr_homology
127 read(inp,'(a)') pdbfile
128 c Next stament causes error upon compilation (?)
129 c if(me.eq.king.or. .not. out1file)
130 c write (iout,'(2a)') 'PDB data will be read from file ',
131 c & pdbfile(:ilen(pdbfile))
132 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
133 & pdbfile(:ilen(pdbfile))
134 open(ipdbin,file=pdbfile,status='old',err=33)
136 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
137 & pdbfile(:ilen(pdbfile))
140 c print *,'Begin reading pdb data'
142 c Files containing res sim or local scores (former containing sigmas)
145 write(kic2,'(bz,i2.2)') k
147 tpl_k_rescore="template"//kic2//".sco"
152 call readpdb_template(k)
161 do i=1,2*nres_chomo(k)
168 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
169 & (crefjlee(j,i+nres),j=1,3)
171 write (iout,*) "read_constr_homology: after reading pdb file"
176 c Distance restraints
179 C Copy the coordinates from reference coordinates (?)
180 do i=1,2*nres_chomo(k)
183 c write (iout,*) "c(",j,i,") =",c(j,i)
187 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
189 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
190 open (ientin,file=tpl_k_rescore,status='old')
191 if (nnt.gt.1) rescore(k,1)=0.0d0
192 do irec=nnt,nct ! loop for reading res sim
194 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
195 & rescore3_tmp,idomain_tmp
197 idomain(k,i_tmp)=idomain_tmp
198 rescore(k,i_tmp)=rescore_tmp
199 rescore2(k,i_tmp)=rescore2_tmp
200 rescore3(k,i_tmp)=rescore3_tmp
201 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
202 & i_tmp,rescore2_tmp,rescore_tmp,
203 & rescore3_tmp,idomain_tmp
206 read (ientin,*,end=1401) rescore_tmp
208 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
209 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
210 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
215 if (waga_dist.ne.0.0d0) then
223 distal=dsqrt(x12*x12+y12*y12+z12*z12)
224 c write (iout,*) k,i,j,distal,dist2_cut
226 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
227 & .and. distal.le.dist2_cut ) then
233 c write (iout,*) "k",k
234 c write (iout,*) "i",i," j",j," constr_homology",
242 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
244 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
245 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
246 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
248 if (odl(k,ii).le.dist_cut) then
249 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
252 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
253 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
255 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
256 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
260 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
269 c write (iout,*) "Distance restraints set"
272 c Theta, dihedral and SC retraints
274 if (waga_angle.gt.0.0d0) then
275 c open (ientin,file=tpl_k_sigma_dih,status='old')
276 c do irec=1,maxres-3 ! loop for reading sigma_dih
277 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
278 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
279 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
280 c & sigma_dih(k,i+nnt-1)
285 if (idomain(k,i).eq.0) then
289 dih(k,i)=phiref(i) ! right?
290 c read (ientin,*) sigma_dih(k,i) ! original variant
291 c write (iout,*) "dih(",k,i,") =",dih(k,i)
292 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
293 c & "rescore(",k,i-1,") =",rescore(k,i-1),
294 c & "rescore(",k,i-2,") =",rescore(k,i-2),
295 c & "rescore(",k,i-3,") =",rescore(k,i-3)
297 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
298 & rescore(k,i-2)+rescore(k,i-3))/4.0
299 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
300 c write (iout,*) "Raw sigmas for dihedral angle restraints"
301 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
302 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
303 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
304 c Instead of res sim other local measure of b/b str reliability possible
305 if (sigma_dih(k,i).ne.0)
306 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
307 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
311 c write (iout,*) "Dihedral angle restraints set"
314 if (waga_theta.gt.0.0d0) then
315 c open (ientin,file=tpl_k_sigma_theta,status='old')
316 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
317 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
318 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
319 c & sigma_theta(k,i+nnt-1)
324 do i = nnt+2,nct ! right? without parallel.
325 c do i = i=1,nres ! alternative for bounds acc to readpdb?
326 c do i=ithet_start,ithet_end ! with FG parallel.
327 if (idomain(k,i).eq.0) then
331 thetatpl(k,i)=thetaref(i)
332 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
333 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
334 c & "rescore(",k,i-1,") =",rescore(k,i-1),
335 c & "rescore(",k,i-2,") =",rescore(k,i-2)
336 c read (ientin,*) sigma_theta(k,i) ! 1st variant
337 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
338 & rescore(k,i-2))/3.0
339 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
340 if (sigma_theta(k,i).ne.0)
341 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
343 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
344 c rescore(k,i-2) ! right expression ?
345 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
348 c write (iout,*) "Angle restraints set"
351 if (waga_d.gt.0.0d0) then
352 c open (ientin,file=tpl_k_sigma_d,status='old')
353 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
354 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
355 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
356 c & sigma_d(k,i+nnt-1)
360 do i = nnt,nct ! right? without parallel.
361 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
362 c do i=loc_start,loc_end ! with FG parallel.
363 if (itype(i).eq.10) cycle
364 if (idomain(k,i).eq.0 ) then
371 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
372 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
373 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
374 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
375 sigma_d(k,i)=rescore3(k,i) ! right expression ?
376 if (sigma_d(k,i).ne.0)
377 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
379 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
380 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
381 c read (ientin,*) sigma_d(k,i) ! 1st variant
385 c write (iout,*) "SC restraints set"
388 c remove distance restraints not used in any model from the list
389 c shift data in all arrays
391 c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
392 if (waga_dist.ne.0.0d0) then
399 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
400 c & .and. distal.le.dist2_cut ) then
401 c write (iout,*) "i",i," j",j," ii",ii
403 if (ii_in_use(ii).eq.0.and.liiflag.or.
404 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
411 if(i10.eq.lim_odl) i10=i10+1
413 ires_homo(iistart+ki)=ires_homo(ki+i01)
414 jres_homo(iistart+ki)=jres_homo(ki+i01)
415 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
416 do k=1,constr_homology
417 odl(k,iistart+ki)=odl(k,ki+i01)
418 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
419 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
422 iistart=iistart+i10-i01
425 if (ii_in_use(ii).ne.0.and..not.liiflag) then
433 c write (iout,*) "Removing distances completed"
435 endif ! .not. klapaucjusz
437 if (constr_homology.gt.0) call homology_partition
438 c write (iout,*) "After homology_partition"
440 if (constr_homology.gt.0) call init_int_table
441 c write (iout,*) "After init_int_table"
443 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
444 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
448 if (.not.out_template_restr) return
449 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
450 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
451 write (iout,*) "Distance restraints from templates"
453 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
454 & ii,ires_homo(ii),jres_homo(ii),
455 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
456 & ki=1,constr_homology)
458 write (iout,*) "Dihedral angle restraints from templates"
460 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
461 & (rad2deg*dih(ki,i),
462 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
464 write (iout,*) "Virtual-bond angle restraints from templates"
466 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
467 & (rad2deg*thetatpl(ki,i),
468 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
470 write (iout,*) "SC restraints from templates"
472 write(iout,'(i5,100(4f8.2,4x))') i,
473 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
474 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
477 c -----------------------------------------------------------------
480 c----------------------------------------------------------------------
481 subroutine read_klapaucjusz
484 include 'DIMENSIONS.ZSCOPT'
485 include 'DIMENSIONS.FREE'
489 include 'COMMON.SETUP'
490 include 'COMMON.CONTROL'
491 include 'COMMON.HOMOLOGY'
492 include 'COMMON.CHAIN'
493 include 'COMMON.IOUNITS'
495 include 'COMMON.INTERACT'
496 include 'COMMON.NAMES'
497 include 'COMMON.HOMRESTR'
498 character*256 fragfile
499 integer ninclust(maxclust),inclust(max_template,maxclust),
500 & nresclust(maxclust),iresclust(maxres,maxclust)
503 character*24 model_ki_dist, model_ki_angle
504 character*500 controlcard
505 integer i,ii,ik,ki,j,k,l,ll,ii_in_use(maxdim),i_tmp,idomain_tmp,
506 & nclust,iistart,iishift,lim_xx
508 integer idomain(max_template,maxres)
509 double precision distal
510 logical lprn /.true./
516 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
517 double precision, dimension (max_template,maxres) :: rescore
518 double precision, dimension (max_template,maxres) :: rescore2
519 character*24 tpl_k_rescore
526 call getenv("FRAGFILE",fragfile)
527 open(ientin,file=fragfile,status="old",err=10)
528 read(ientin,*) constr_homology,nclust
535 do k=1,constr_homology
536 read(ientin,'(a)') pdbfile
537 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
538 & pdbfile(:ilen(pdbfile))
539 open(ipdbin,file=pdbfile,status='old',err=33)
541 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
542 & pdbfile(:ilen(pdbfile))
547 call readpdb_template(k)
552 c chomo(j,i,k)=c(j,i)
562 read(ientin,*) ninclust(i),nresclust(i)
563 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
564 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
570 c write (iout,*) "CLUSTER",l," members:",ninclust(l)
571 do ll = 1,ninclust(l)
579 idomain(k,iresclust(i,l)+1) = 1
581 idomain(k,iresclust(i,l)) = 1
585 c Distance restraints
588 C Copy the coordinates from reference coordinates (?)
589 c write (iout,*) "k",k
595 c write (iout,*) "c(",j,i,") =",c(j,i)
598 c write (iout,*) "read_klapaucjusz: calling int_from_cart"
600 c write (iout,*) "idomain"
601 c write (iout,'(2i5)') (i,idomain(k,i),i=1,nres)
602 call int_from_cart(.true.,.false.)
603 call sc_loc_geom(.false.)
604 c write (iout,*) "en"
610 if (waga_dist.ne.0.0d0) then
618 distal=dsqrt(x12*x12+y12*y12+z12*z12)
619 c write (iout,*) k,i,j,distal,dist2_cut
621 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
622 & .and. distal.le.dist2_cut ) then
628 c write (iout,*) "k",k
629 c write (iout,*) "i",i," j",j," constr_homology",
637 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
639 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
640 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
641 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
643 if (odl(k,ii).le.dist_cut) then
644 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
647 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
648 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
650 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
651 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
655 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
658 c l_homo(k,ii)=.false.
665 c Theta, dihedral and SC retraints
667 if (waga_angle.gt.0.0d0) then
669 if (idomain(k,i).eq.0) then
674 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
675 & rescore(k,i-2)+rescore(k,i-3))/4.0
676 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
677 c & " sigma_dihed",sigma_dih(k,i)
678 if (sigma_dih(k,i).ne.0)
679 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
684 if (waga_theta.gt.0.0d0) then
686 if (idomain(k,i).eq.0) then
687 c sigma_theta(k,i)=0.0
690 thetatpl(k,i)=thetaref(i)
691 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
692 & rescore(k,i-2))/3.0
693 if (sigma_theta(k,i).ne.0)
694 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
698 if (waga_d.gt.0.0d0) then
700 if (itype(i).eq.10) cycle
701 if (idomain(k,i).eq.0 ) then
708 sigma_d(k,i)=rescore(k,i)
709 if (sigma_d(k,i).ne.0)
710 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
711 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
717 c remove distance restraints not used in any model from the list
718 c shift data in all arrays
720 if (waga_dist.ne.0.0d0) then
726 if (ii_in_use(ii).eq.0.and.liiflag) then
730 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
731 & .not.liiflag.and.ii.eq.lim_odl) then
732 if (ii.eq.lim_odl) then
738 do ki=iistart,lim_odl-iishift
739 ires_homo(ki)=ires_homo(ki+iishift)
740 jres_homo(ki)=jres_homo(ki+iishift)
741 ii_in_use(ki)=ii_in_use(ki+iishift)
742 do k=1,constr_homology
743 odl(k,ki)=odl(k,ki+iishift)
744 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
745 l_homo(k,ki)=l_homo(k,ki+iishift)
749 lim_odl=lim_odl-iishift
755 write (iout,*) "ires_homo and jres_homo arrays, lim_odl",lim_odl
757 write (iout,*) i,ires_homo(i),jres_homo(i)
761 10 stop "Error in fragment file"