1 subroutine read_constr_homology
8 include 'COMMON.CONTROL'
9 include 'COMMON.HOMOLOGY'
10 include 'COMMON.CHAIN'
11 include 'COMMON.IOUNITS'
13 include 'COMMON.QRESTR'
15 include 'COMMON.INTERACT'
16 include 'COMMON.NAMES'
23 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
25 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
26 c & sigma_odl_temp(maxres,maxres,max_template)
28 character*24 model_ki_dist, model_ki_angle
29 character*500 controlcard
30 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
31 & ik,iistart,nres_temp
34 logical liiflag,lfirst
37 c FP - Nov. 2014 Temporary specifications for new vars
39 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
41 double precision, dimension (max_template,maxres) :: rescore
42 double precision, dimension (max_template,maxres) :: rescore2
43 double precision, dimension (max_template,maxres) :: rescore3
44 double precision distal
45 character*24 pdbfile,tpl_k_rescore
46 c -----------------------------------------------------------------
47 c Reading multiple PDB ref structures and calculation of retraints
48 c not using pre-computed ones stored in files model_ki_{dist,angle}
50 c -----------------------------------------------------------------
53 c Alternative: reading from input
54 call card_concat(controlcard)
55 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
56 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
57 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
58 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
59 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
60 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
61 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
62 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
63 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
64 if(.not.read2sigma.and.start_from_model) then
65 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
66 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
67 start_from_model=.false.
69 if(start_from_model .and. (me.eq.king .or. .not. out1file))
70 & write(iout,*) 'START_FROM_MODELS is ON'
71 if(start_from_model .and. rest) then
72 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
73 write(iout,*) 'START_FROM_MODELS is OFF'
74 write(iout,*) 'remove restart keyword from input'
77 if (homol_nset.gt.1)then
78 call card_concat(controlcard)
79 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
80 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
81 write(iout,*) "iset homology_weight "
83 write(iout,*) i,waga_homology(i)
86 iset=mod(kolor,homol_nset)+1
92 cd write (iout,*) "nnt",nnt," nct",nct
99 c write(iout,*) 'nnt=',nnt,'nct=',nct
102 do k=1,constr_homology
115 if (read_homol_frag) then
116 call read_klapaucjusz
119 do k=1,constr_homology
121 read(inp,'(a)') pdbfile
122 if(me.eq.king .or. .not. out1file)
123 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
124 & pdbfile(:ilen(pdbfile))
125 open(ipdbin,file=pdbfile,status='old',err=33)
127 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
128 & pdbfile(:ilen(pdbfile))
131 c print *,'Begin reading pdb data'
133 c Files containing res sim or local scores (former containing sigmas)
136 write(kic2,'(bz,i2.2)') k
138 tpl_k_rescore="template"//kic2//".sco"
143 call readpdb_template(k)
150 c Distance restraints
153 C Copy the coordinates from reference coordinates (?)
154 do i=1,2*nres_chomo(k)
157 c write (iout,*) "c(",j,i,") =",c(j,i)
161 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
163 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
164 open (ientin,file=tpl_k_rescore,status='old')
165 if (nnt.gt.1) rescore(k,1)=0.0d0
166 do irec=nnt,nct ! loop for reading res sim
168 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
169 & rescore3_tmp,idomain_tmp
171 idomain(k,i_tmp)=idomain_tmp
172 rescore(k,i_tmp)=rescore_tmp
173 rescore2(k,i_tmp)=rescore2_tmp
174 rescore3(k,i_tmp)=rescore3_tmp
175 if (.not. out1file .or. me.eq.king)
176 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
177 & i_tmp,rescore2_tmp,rescore_tmp,
178 & rescore3_tmp,idomain_tmp
181 read (ientin,*,end=1401) rescore_tmp
183 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
184 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
185 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
190 if (waga_dist.ne.0.0d0) then
198 distal=dsqrt(x12*x12+y12*y12+z12*z12)
199 c write (iout,*) k,i,j,distal,dist2_cut
201 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
202 & .and. distal.le.dist2_cut ) then
208 c write (iout,*) "k",k
209 c write (iout,*) "i",i," j",j," constr_homology",
217 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
219 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
220 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
221 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
223 if (odl(k,ii).le.dist_cut) then
224 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
227 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
228 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
230 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
231 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
235 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
244 c write (iout,*) "Distance restraints set"
247 c Theta, dihedral and SC retraints
249 if (waga_angle.gt.0.0d0) then
250 c open (ientin,file=tpl_k_sigma_dih,status='old')
251 c do irec=1,maxres-3 ! loop for reading sigma_dih
252 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
253 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
254 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
255 c & sigma_dih(k,i+nnt-1)
260 if (idomain(k,i).eq.0) then
264 dih(k,i)=phiref(i) ! right?
265 c read (ientin,*) sigma_dih(k,i) ! original variant
266 c write (iout,*) "dih(",k,i,") =",dih(k,i)
267 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
268 c & "rescore(",k,i-1,") =",rescore(k,i-1),
269 c & "rescore(",k,i-2,") =",rescore(k,i-2),
270 c & "rescore(",k,i-3,") =",rescore(k,i-3)
272 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
273 & rescore(k,i-2)+rescore(k,i-3))/4.0
274 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
275 c write (iout,*) "Raw sigmas for dihedral angle restraints"
276 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
277 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
278 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
279 c Instead of res sim other local measure of b/b str reliability possible
280 if (sigma_dih(k,i).ne.0)
281 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
282 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
286 c write (iout,*) "Dihedral angle restraints set"
289 if (waga_theta.gt.0.0d0) then
290 c open (ientin,file=tpl_k_sigma_theta,status='old')
291 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
292 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
293 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
294 c & sigma_theta(k,i+nnt-1)
299 do i = nnt+2,nct ! right? without parallel.
300 c do i = i=1,nres ! alternative for bounds acc to readpdb?
301 c do i=ithet_start,ithet_end ! with FG parallel.
302 if (idomain(k,i).eq.0) then
306 thetatpl(k,i)=thetaref(i)
307 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
308 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
309 c & "rescore(",k,i-1,") =",rescore(k,i-1),
310 c & "rescore(",k,i-2,") =",rescore(k,i-2)
311 c read (ientin,*) sigma_theta(k,i) ! 1st variant
312 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
313 & rescore(k,i-2))/3.0
314 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
315 if (sigma_theta(k,i).ne.0)
316 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
318 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
319 c rescore(k,i-2) ! right expression ?
320 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
323 c write (iout,*) "Angle restraints set"
326 if (waga_d.gt.0.0d0) then
327 c open (ientin,file=tpl_k_sigma_d,status='old')
328 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
329 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
330 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
331 c & sigma_d(k,i+nnt-1)
335 do i = nnt,nct ! right? without parallel.
336 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
337 c do i=loc_start,loc_end ! with FG parallel.
338 if (itype(i).eq.10) cycle
339 if (idomain(k,i).eq.0 ) then
346 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
347 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
348 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
349 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
350 sigma_d(k,i)=rescore3(k,i) ! right expression ?
351 if (sigma_d(k,i).ne.0)
352 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
354 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
355 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
356 c read (ientin,*) sigma_d(k,i) ! 1st variant
360 c write (iout,*) "SC restraints set"
363 c remove distance restraints not used in any model from the list
364 c shift data in all arrays
366 c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
367 if (waga_dist.ne.0.0d0) then
374 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
375 c & .and. distal.le.dist2_cut ) then
376 c write (iout,*) "i",i," j",j," ii",ii
378 if (ii_in_use(ii).eq.0.and.liiflag.or.
379 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
386 if(i10.eq.lim_odl) i10=i10+1
388 ires_homo(iistart+ki)=ires_homo(ki+i01)
389 jres_homo(iistart+ki)=jres_homo(ki+i01)
390 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
391 do k=1,constr_homology
392 odl(k,iistart+ki)=odl(k,ki+i01)
393 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
394 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
397 iistart=iistart+i10-i01
400 if (ii_in_use(ii).ne.0.and..not.liiflag) then
408 c write (iout,*) "Removing distances completed"
410 endif ! .not. klapaucjusz
412 if (constr_homology.gt.0) call homology_partition
413 c write (iout,*) "After homology_partition"
415 if (constr_homology.gt.0) call init_int_table
416 c write (iout,*) "After init_int_table"
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.out_template_restr) return
424 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
425 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
462 include 'COMMON.SETUP'
463 include 'COMMON.CONTROL'
464 include 'COMMON.HOMOLOGY'
465 include 'COMMON.CHAIN'
466 include 'COMMON.IOUNITS'
469 include 'COMMON.INTERACT'
470 include 'COMMON.NAMES'
471 character*256 fragfile
472 integer ninclust(maxclust),inclust(max_template,maxclust),
473 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
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 & ik,ll,ii,kk,iistart,iishift,lim_xx
480 double precision distal
481 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 pdbfile,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))
518 call readpdb_template(k)
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 (?)
559 c write (iout,*) "c(",j,i,") =",c(j,i)
562 call int_from_cart(.true.,.false.)
563 call sc_loc_geom(.false.)
569 if (waga_dist.ne.0.0d0) then
577 distal=dsqrt(x12*x12+y12*y12+z12*z12)
578 c write (iout,*) k,i,j,distal,dist2_cut
580 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
581 & .and. distal.le.dist2_cut ) then
587 c write (iout,*) "k",k
588 c write (iout,*) "i",i," j",j," constr_homology",
596 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
598 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
599 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
600 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
602 if (odl(k,ii).le.dist_cut) then
603 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
606 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
607 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
609 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
610 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
614 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
617 c l_homo(k,ii)=.false.
624 c Theta, dihedral and SC retraints
626 if (waga_angle.gt.0.0d0) then
628 if (idomain(k,i).eq.0) then
633 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
634 & rescore(k,i-2)+rescore(k,i-3))/4.0
635 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
636 c & " sigma_dihed",sigma_dih(k,i)
637 if (sigma_dih(k,i).ne.0)
638 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
643 if (waga_theta.gt.0.0d0) then
645 if (idomain(k,i).eq.0) then
646 c sigma_theta(k,i)=0.0
649 thetatpl(k,i)=thetaref(i)
650 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
651 & rescore(k,i-2))/3.0
652 if (sigma_theta(k,i).ne.0)
653 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
657 if (waga_d.gt.0.0d0) then
659 if (itype(i).eq.10) cycle
660 if (idomain(k,i).eq.0 ) then
667 sigma_d(k,i)=rescore(k,i)
668 if (sigma_d(k,i).ne.0)
669 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
670 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
676 c remove distance restraints not used in any model from the list
677 c shift data in all arrays
679 if (waga_dist.ne.0.0d0) then
685 if (ii_in_use(ii).eq.0.and.liiflag) then
689 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
690 & .not.liiflag.and.ii.eq.lim_odl) then
691 if (ii.eq.lim_odl) then
697 do ki=iistart,lim_odl-iishift
698 ires_homo(ki)=ires_homo(ki+iishift)
699 jres_homo(ki)=jres_homo(ki+iishift)
700 ii_in_use(ki)=ii_in_use(ki+iishift)
701 do k=1,constr_homology
702 odl(k,ki)=odl(k,ki+iishift)
703 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
704 l_homo(k,ki)=l_homo(k,ki+iishift)
708 lim_odl=lim_odl-iishift
715 10 stop "Error in fragment file"