5 implicit real*8 (a-h,o-z)
7 include 'DIMENSIONS.ZSCOPT'
8 include 'DIMENSIONS.FREE'
9 include 'COMMON.IOUNITS'
12 c include 'include_unres/COMMON.VAR'
13 include 'COMMON.INTERACT'
14 include 'COMMON.LOCAL'
15 include 'COMMON.NAMES'
16 include 'COMMON.CHAIN'
17 include 'COMMON.FFIELD'
18 include 'COMMON.SBRIDGE'
19 include 'COMMON.TORCNSTR'
20 include 'COMMON.CONTROL'
21 character*4 sequence(maxres)
23 double precision x(maxvar)
24 character*320 controlcard,ucase
25 dimension itype_pdb(maxres)
27 call card_concat(controlcard,.true.)
28 call reada(controlcard,'SCAL14',scal14,0.4d0)
29 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
30 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
31 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
33 call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
34 call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
35 call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
36 call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
37 write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
38 & " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
39 r0_corr=cutoff_corr-delt_corr
40 call readi(controlcard,"NRES",nres,0)
41 iscode=index(controlcard,"ONE_LETTER")
43 write (iout,*) "Error: no residues in molecule"
46 if (nres.gt.maxres) then
47 write (iout,*) "Error: too many residues",nres,maxres
49 write(iout,*) 'nres=',nres
50 C Read sequence of the protein
52 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
54 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
56 C Convert sequence to numeric code
58 itype(i)=rescode(i,sequence(i),iscode)
60 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
62 & "Glycine is the first full residue, initial dummy deleted"
68 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
70 & "Glycine is the last full residue, terminal dummy deleted"
73 write (iout,*) "Numeric code:"
74 write (iout,'(20i4)') (itype(i),i=1,nres)
77 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
79 if (itype(i).eq.21) then
83 else if (itype(i+1).ne.20) then
85 else if (itype(i).ne.20) then
94 if (with_dihed_constr) then
96 read (inp,*) ndih_constr
97 if (ndih_constr.gt.0) then
99 write (iout,*) 'FTORS',ftors
100 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
102 & 'There are',ndih_constr,' constraints on phi angles.'
104 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
107 phi0(i)=deg2rad*phi0(i)
108 drange(i)=deg2rad*drange(i)
116 if (itype(1).eq.21) nnt=2
117 if (itype(nres).eq.21) nct=nct-1
118 write(iout,*) 'NNT=',NNT,' NCT=',NCT
120 C Juyong:READ init_vars
121 C Initialize variables!
122 C Juyong:READ read_info
123 C READ fragment information!!
124 C both routines should be in dfa.F file!!
126 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
127 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
128 write (iout,*) "Calling init_dfa_vars"
131 write (iout,*) 'init_dfa_vars finished!'
134 write (iout,*) 'read_dfa_info finished!'
138 c Read distance restraints
139 if (constr_dist.gt.0) then
140 if (refstr) call read_ref_structure(*11)
141 call read_dist_constr
145 if (constr_homology.gt.0) then
146 c write (iout,*) "About to call read_constr_homology"
148 call read_constr_homology
149 c write (iout,*) "Exit read_constr_homology"
151 if (indpdb.gt.0 .or. pdbref) then
155 cref(j,i)=crefjlee(j,i)
160 write (iout,*) "Array C"
162 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
163 & (c(j,i+nres),j=1,3)
165 write (iout,*) "Array Cref"
167 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
168 & (cref(j,i+nres),j=1,3)
172 call int_from_cart1(.false.)
173 call sc_loc_geom(.false.)
177 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
181 dc(j,i)=c(j,i+1)-c(j,i)
182 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
187 dc(j,i+nres)=c(j,i+nres)-c(j,i)
188 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
200 write (iout,'(/a,i3,a)') 'The chain contains',ns,
201 & ' disulfide-bridging cysteines.'
202 write (iout,'(20i4)') (iss(i),i=1,ns)
203 write (iout,'(/a/)') 'Pre-formed links are:'
209 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
210 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
211 & dhpb(i),ebr,forcon(i)
216 11 stop "Error reading reference structure"
218 c-----------------------------------------------------------------------------
219 logical function seq_comp(itypea,itypeb,length)
221 integer length,itypea(length),itypeb(length)
224 if (itypea(i).ne.itypeb(i)) then
232 c-----------------------------------------------------------------------------
233 subroutine read_bridge
234 C Read information about disulfide bridges.
235 implicit real*8 (a-h,o-z)
237 include 'DIMENSIONS.ZSCOPT'
238 include 'COMMON.IOUNITS'
241 include 'COMMON.INTERACT'
242 include 'COMMON.NAMES'
243 include 'COMMON.CHAIN'
244 include 'COMMON.FFIELD'
245 include 'COMMON.SBRIDGE'
246 C Read bridging residues.
247 read (inp,*) ns,(iss(i),i=1,ns)
249 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
250 C Check whether the specified bridging residues are cystines.
252 if (itype(iss(i)).ne.1) then
253 write (iout,'(2a,i3,a)')
254 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
255 & ' can form a disulfide bridge?!!!'
256 write (*,'(2a,i3,a)')
257 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
258 & ' can form a disulfide bridge?!!!'
262 C Read preformed bridges.
264 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
265 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
268 C Check if the residues involved in bridges are in the specified list of
272 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
273 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
274 write (iout,'(a,i3,a)') 'Disulfide pair',i,
275 & ' contains residues present in other pairs.'
276 write (*,'(a,i3,a)') 'Disulfide pair',i,
277 & ' contains residues present in other pairs.'
282 if (ihpb(i).eq.iss(j)) goto 10
284 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
287 if (jhpb(i).eq.iss(j)) goto 20
289 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
300 if (ns.gt.0.and.dyn_ss) then
301 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
304 C /06/28/2013 Adasko: nss number of full SS bonds
307 forcon(i-nss)=forcon(i)
314 dyn_ss_mask(iss(i))=.true.
315 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
316 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
321 c------------------------------------------------------------------------------
322 subroutine read_angles(kanal,iscor,energ,iprot,*)
323 implicit real*8 (a-h,o-z)
325 include 'DIMENSIONS.ZSCOPT'
326 include 'COMMON.INTERACT'
327 include 'COMMON.SBRIDGE'
330 include 'COMMON.CHAIN'
331 include 'COMMON.IOUNITS'
333 read(kanal,'(a80)',end=10,err=10) lineh
334 read(lineh(:5),*,err=8) ic
335 read(lineh(6:),*,err=8) energ
338 print *,'error, assuming e=1d10',lineh
342 read(lineh(18:),*,end=10,err=10) nss
344 read (lineh(20:),*,end=10,err=10)
345 & (IHPB(I),JHPB(I),I=1,NSS),iscor
347 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
348 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
351 c print *,"energy",energ," iscor",iscor
352 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
353 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
354 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
355 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
357 theta(i)=deg2rad*theta(i)
358 phi(i)=deg2rad*phi(i)
359 alph(i)=deg2rad*alph(i)
360 omeg(i)=deg2rad*omeg(i)
365 c-------------------------------------------------------------------------------
366 subroutine read_dist_constr
367 implicit real*8 (a-h,o-z)
369 include 'DIMENSIONS.ZSCOPT'
370 include 'DIMENSIONS.FREE'
371 include 'COMMON.CONTROL'
372 include 'COMMON.CHAIN'
373 include 'COMMON.IOUNITS'
374 include 'COMMON.SBRIDGE'
375 integer ifrag_(2,100),ipair_(2,100)
376 double precision wfrag_(100),wpair_(100)
377 character*500 controlcard
378 logical normalize,next
380 double precision xlink(4,0:4) /
382 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
383 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
384 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
385 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
386 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
387 write (iout,*) "Calling read_dist_constr"
388 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
394 call card_concat(controlcard,.true.)
395 next = index(controlcard,"NEXT").gt.0
396 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
397 write (iout,*) "restr_type",restr_type
398 call readi(controlcard,"NFRAG",nfrag_,0)
399 call readi(controlcard,"NFRAG",nfrag_,0)
400 call readi(controlcard,"NPAIR",npair_,0)
401 call readi(controlcard,"NDIST",ndist_,0)
402 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
403 if (restr_type.eq.10)
404 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
405 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
406 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
407 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
408 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
409 normalize = index(controlcard,"NORMALIZE").gt.0
410 write (iout,*) "WBOLTZD",wboltzd
411 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
412 c write (iout,*) "IFRAG"
414 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
416 c write (iout,*) "IPAIR"
418 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
420 if (nfrag_.gt.0) then
422 read(inp,'(a)') pdbfile
424 & "Distance restraints will be constructed from structure ",pdbfile
425 open(ipdbin,file=pdbfile,status='old',err=11)
431 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
432 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
433 & ifrag_(2,i)=nstart_sup+nsup-1
434 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
436 if (wfrag_(i).eq.0.0d0) cycle
437 do j=ifrag_(1,i),ifrag_(2,i)-1
439 c write (iout,*) "j",j," k",k
441 if (restr_type.eq.1) then
447 forcon(nhpb)=wfrag_(i)
448 else if (constr_dist.eq.2) then
449 if (ddjk.le.dist_cut) then
455 forcon(nhpb)=wfrag_(i)
457 else if (restr_type.eq.3) then
463 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
465 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
466 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
471 if (wpair_(i).eq.0.0d0) cycle
479 do j=ifrag_(1,ii),ifrag_(2,ii)
480 do k=ifrag_(1,jj),ifrag_(2,jj)
481 if (restr_type.eq.1) then
487 forcon(nhpb)=wfrag_(i)
488 else if (constr_dist.eq.2) then
489 if (ddjk.le.dist_cut) then
495 forcon(nhpb)=wfrag_(i)
497 else if (restr_type.eq.3) then
503 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
505 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
506 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
512 write (iout,*) "Distance restraints as read from input"
514 if (restr_type.eq.11) then
515 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
516 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
517 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
518 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
521 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
522 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
523 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
524 if (ibecarb(nhpb).gt.0) then
525 ihpb(nhpb)=ihpb(nhpb)+nres
526 jhpb(nhpb)=jhpb(nhpb)+nres
528 else if (constr_dist.eq.10) then
529 c Cross-lonk Markov-like potential
530 call card_concat(controlcard,.true.)
531 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
532 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
534 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
535 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
536 if (index(controlcard,"ZL").gt.0) then
538 else if (index(controlcard,"ADH").gt.0) then
540 else if (index(controlcard,"PDH").gt.0) then
542 else if (index(controlcard,"DSS").gt.0) then
547 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
548 & xlink(1,link_type))
549 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
550 & xlink(2,link_type))
551 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
552 & xlink(3,link_type))
553 call reada(controlcard,"SIGMA",forcon(nhpb+1),
554 & xlink(4,link_type))
555 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
556 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
557 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
558 if (forcon(nhpb+1).le.0.0d0 .or.
559 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
562 if (ibecarb(nhpb).gt.0) then
563 ihpb(nhpb)=ihpb(nhpb)+nres
564 jhpb(nhpb)=jhpb(nhpb)+nres
566 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
567 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
568 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
572 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
573 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
574 if (forcon(nhpb+1).gt.0.0d0) then
576 if (dhpb1(nhpb).eq.0.0d0) then
581 if (ibecarb(nhpb).gt.0) then
582 ihpb(nhpb)=ihpb(nhpb)+nres
583 jhpb(nhpb)=jhpb(nhpb)+nres
585 if (dhpb(nhpb).eq.0.0d0)
586 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
588 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
589 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
591 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
592 C if (forcon(nhpb+1).gt.0.0d0) then
594 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
602 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
603 & fordepthmax=fordepth(i)
606 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
609 if (nhpb.gt.nss) then
610 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
611 & "The following",nhpb-nss,
612 & " distance restraints have been imposed:",
613 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
616 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
617 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
625 11 write (iout,*)"read_dist_restr: error reading reference structure"
628 c====-------------------------------------------------------------------
629 subroutine read_constr_homology
632 include 'DIMENSIONS.ZSCOPT'
633 include 'DIMENSIONS.FREE'
637 include 'COMMON.SETUP'
638 include 'COMMON.CONTROL'
639 include 'COMMON.CHAIN'
640 include 'COMMON.IOUNITS'
642 include 'COMMON.INTERACT'
643 include 'COMMON.NAMES'
644 include 'COMMON.HOMRESTR'
649 c include 'include_unres/COMMON.VAR'
652 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
654 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
655 c & sigma_odl_temp(maxres,maxres,max_template)
657 character*24 model_ki_dist, model_ki_angle
658 character*500 controlcard
659 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
660 integer idomain(max_template,maxres)
661 logical lprn /.true./
664 logical unres_pdb,liiflag
666 c FP - Nov. 2014 Temporary specifications for new vars
668 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
669 double precision, dimension (max_template,maxres) :: rescore
670 double precision, dimension (max_template,maxres) :: rescore2
671 character*24 tpl_k_rescore
672 c -----------------------------------------------------------------
673 c Reading multiple PDB ref structures and calculation of retraints
674 c not using pre-computed ones stored in files model_ki_{dist,angle}
676 c -----------------------------------------------------------------
679 c Alternative: reading from input
680 call card_concat(controlcard,.true.)
681 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
682 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
683 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
684 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
685 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
686 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
687 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
688 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
689 call readi(controlcard,"IHSET",ihset,1)
690 if (homol_nset.gt.1)then
691 call card_concat(controlcard,.true.)
692 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
693 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
694 write(iout,*) "iset homology_weight "
696 c write(iout,*) i,waga_homology(i)
699 iset=mod(kolor,homol_nset)+1
704 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
706 cd write (iout,*) "nnt",nnt," nct",nct
718 c Reading HM global scores (prob not required)
721 do k=1,constr_homology
725 c open (4,file="HMscore")
726 c do k=1,constr_homology
727 c read (4,*,end=521) hmscore_tmp
728 c hmscore(k)=hmscore_tmp ! Another transformation can be used
729 c write(*,*) "Model", k, ":", hmscore(k)
740 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
742 if (read_homol_frag) then
743 call read_klapaucjusz
746 do k=1,constr_homology
748 read(inp,'(a)') pdbfile
749 c Next stament causes error upon compilation (?)
750 c if(me.eq.king.or. .not. out1file)
751 c write (iout,'(2a)') 'PDB data will be read from file ',
752 c & pdbfile(:ilen(pdbfile))
753 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
754 & pdbfile(:ilen(pdbfile))
755 open(ipdbin,file=pdbfile,status='old',err=33)
757 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
758 & pdbfile(:ilen(pdbfile))
761 c print *,'Begin reading pdb data'
763 c Files containing res sim or local scores (former containing sigmas)
766 write(kic2,'(bz,i2.2)') k
768 tpl_k_rescore="template"//kic2//".sco"
779 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
780 & (crefjlee(j,i+nres),j=1,3)
783 write (iout,*) "read_constr_homology: after reading pdb file"
787 c Distance restraints
790 C Copy the coordinates from reference coordinates (?)
794 c write (iout,*) "c(",j,i,") =",c(j,i)
798 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
800 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
801 open (ientin,file=tpl_k_rescore,status='old')
802 if (nnt.gt.1) rescore(k,1)=0.0d0
803 do irec=nnt,nct ! loop for reading res sim
805 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
808 idomain(k,i_tmp)=idomain_tmp
809 rescore(k,i_tmp)=rescore_tmp
810 rescore2(k,i_tmp)=rescore2_tmp
811 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
812 & i_tmp,rescore2_tmp,rescore_tmp,
816 read (ientin,*,end=1401) rescore_tmp
818 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
819 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
820 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
825 if (waga_dist.ne.0.0d0) then
833 distal=dsqrt(x12*x12+y12*y12+z12*z12)
834 c write (iout,*) k,i,j,distal,dist2_cut
836 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
837 & .and. distal.le.dist2_cut ) then
843 c write (iout,*) "k",k
844 c write (iout,*) "i",i," j",j," constr_homology",
852 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
854 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
855 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
856 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
858 if (odl(k,ii).le.dist_cut) then
859 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
862 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
863 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
865 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
866 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
870 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
880 c Theta, dihedral and SC retraints
882 if (waga_angle.gt.0.0d0) then
883 c open (ientin,file=tpl_k_sigma_dih,status='old')
884 c do irec=1,maxres-3 ! loop for reading sigma_dih
885 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
886 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
887 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
888 c & sigma_dih(k,i+nnt-1)
893 if (idomain(k,i).eq.0) then
897 dih(k,i)=phiref(i) ! right?
898 c read (ientin,*) sigma_dih(k,i) ! original variant
899 c write (iout,*) "dih(",k,i,") =",dih(k,i)
900 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
901 c & "rescore(",k,i-1,") =",rescore(k,i-1),
902 c & "rescore(",k,i-2,") =",rescore(k,i-2),
903 c & "rescore(",k,i-3,") =",rescore(k,i-3)
905 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
906 & rescore(k,i-2)+rescore(k,i-3))/4.0
907 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
908 c write (iout,*) "Raw sigmas for dihedral angle restraints"
909 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
910 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
911 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
912 c Instead of res sim other local measure of b/b str reliability possible
913 if (sigma_dih(k,i).ne.0)
914 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
915 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
920 if (waga_theta.gt.0.0d0) then
921 c open (ientin,file=tpl_k_sigma_theta,status='old')
922 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
923 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
924 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
925 c & sigma_theta(k,i+nnt-1)
930 do i = nnt+2,nct ! right? without parallel.
931 c do i = i=1,nres ! alternative for bounds acc to readpdb?
932 c do i=ithet_start,ithet_end ! with FG parallel.
933 if (idomain(k,i).eq.0) then
937 thetatpl(k,i)=thetaref(i)
938 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
939 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
940 c & "rescore(",k,i-1,") =",rescore(k,i-1),
941 c & "rescore(",k,i-2,") =",rescore(k,i-2)
942 c read (ientin,*) sigma_theta(k,i) ! 1st variant
943 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
944 & rescore(k,i-2))/3.0
945 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
946 if (sigma_theta(k,i).ne.0)
947 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
949 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
950 c rescore(k,i-2) ! right expression ?
951 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
955 if (waga_d.gt.0.0d0) then
956 c open (ientin,file=tpl_k_sigma_d,status='old')
957 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
958 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
959 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
960 c & sigma_d(k,i+nnt-1)
964 do i = nnt,nct ! right? without parallel.
965 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
966 c do i=loc_start,loc_end ! with FG parallel.
967 if (itype(i).eq.10) cycle
968 if (idomain(k,i).eq.0 ) then
975 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
976 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
977 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
978 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
979 sigma_d(k,i)=rescore(k,i) ! right expression ?
980 if (sigma_d(k,i).ne.0)
981 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
983 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
984 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
985 c read (ientin,*) sigma_d(k,i) ! 1st variant
990 c remove distance restraints not used in any model from the list
991 c shift data in all arrays
993 if (waga_dist.ne.0.0d0) then
999 if (ii_in_use(ii).eq.0.and.liiflag) then
1003 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1004 & .not.liiflag.and.ii.eq.lim_odl) then
1005 if (ii.eq.lim_odl) then
1006 iishift=ii-iistart+1
1011 do ki=iistart,lim_odl-iishift
1012 ires_homo(ki)=ires_homo(ki+iishift)
1013 jres_homo(ki)=jres_homo(ki+iishift)
1014 ii_in_use(ki)=ii_in_use(ki+iishift)
1015 do k=1,constr_homology
1016 odl(k,ki)=odl(k,ki+iishift)
1017 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1018 l_homo(k,ki)=l_homo(k,ki+iishift)
1022 lim_odl=lim_odl-iishift
1028 endif ! .not. klapaucjusz
1030 if (constr_homology.gt.0) call homology_partition
1031 if (constr_homology.gt.0) call init_int_table
1032 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1033 cd & "lim_xx=",lim_xx
1034 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1035 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1039 if (.not.lprn) return
1040 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1041 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1042 write (iout,*) "Distance restraints from templates"
1044 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1045 & ii,ires_homo(ii),jres_homo(ii),
1046 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1047 & ki=1,constr_homology)
1049 write (iout,*) "Dihedral angle restraints from templates"
1051 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1052 & (rad2deg*dih(ki,i),
1053 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1055 write (iout,*) "Virtual-bond angle restraints from templates"
1057 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1058 & (rad2deg*thetatpl(ki,i),
1059 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1061 write (iout,*) "SC restraints from templates"
1063 write(iout,'(i5,100(4f8.2,4x))') i,
1064 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1065 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1068 c -----------------------------------------------------------------
1071 c----------------------------------------------------------------------
1072 subroutine read_klapaucjusz
1074 include 'DIMENSIONS'
1075 include 'DIMENSIONS.ZSCOPT'
1076 include 'DIMENSIONS.FREE'
1080 include 'COMMON.SETUP'
1081 include 'COMMON.CONTROL'
1082 include 'COMMON.CHAIN'
1083 include 'COMMON.IOUNITS'
1084 include 'COMMON.GEO'
1085 include 'COMMON.INTERACT'
1086 include 'COMMON.NAMES'
1087 include 'COMMON.HOMRESTR'
1088 character*256 fragfile
1089 integer ninclust(maxclust),inclust(max_template,maxclust),
1090 & nresclust(maxclust),iresclust(maxres,maxclust)
1093 character*24 model_ki_dist, model_ki_angle
1094 character*500 controlcard
1095 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1096 integer idomain(max_template,maxres)
1097 logical lprn /.true./
1100 logical unres_pdb,liiflag
1103 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1104 double precision, dimension (max_template,maxres) :: rescore
1105 double precision, dimension (max_template,maxres) :: rescore2
1106 character*24 tpl_k_rescore
1109 c For new homol impl
1111 include 'COMMON.VAR'
1113 double precision chomo(3,maxres2+2,max_template)
1114 call getenv("FRAGFILE",fragfile)
1115 open(ientin,file=fragfile,status="old",err=10)
1116 read(ientin,*) constr_homology,nclust
1122 do k=1,constr_homology
1123 read(ientin,'(a)') pdbfile
1124 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1125 & pdbfile(:ilen(pdbfile))
1126 open(ipdbin,file=pdbfile,status='old',err=33)
1128 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
1129 & pdbfile(:ilen(pdbfile))
1146 read(ientin,*) ninclust(i),nresclust(i)
1147 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1148 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1151 c Loop over clusters
1154 do ll = 1,ninclust(l)
1162 idomain(k,iresclust(i,l)+1) = 1
1164 idomain(k,iresclust(i,l)) = 1
1168 c Distance restraints
1171 C Copy the coordinates from reference coordinates (?)
1175 c write (iout,*) "c(",j,i,") =",c(j,i)
1178 call int_from_cart(.true.,.false.)
1179 call sc_loc_geom(.false.)
1181 thetaref(i)=theta(i)
1184 if (waga_dist.ne.0.0d0) then
1192 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1193 c write (iout,*) k,i,j,distal,dist2_cut
1195 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1196 & .and. distal.le.dist2_cut ) then
1202 c write (iout,*) "k",k
1203 c write (iout,*) "i",i," j",j," constr_homology",
1208 if (read2sigma) then
1211 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1213 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1214 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1215 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1217 if (odl(k,ii).le.dist_cut) then
1218 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1221 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1222 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1224 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1225 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1229 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1232 c l_homo(k,ii)=.false.
1239 c Theta, dihedral and SC retraints
1241 if (waga_angle.gt.0.0d0) then
1243 if (idomain(k,i).eq.0) then
1244 c sigma_dih(k,i)=0.0
1248 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1249 & rescore(k,i-2)+rescore(k,i-3))/4.0
1250 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1251 c & " sigma_dihed",sigma_dih(k,i)
1252 if (sigma_dih(k,i).ne.0)
1253 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1258 if (waga_theta.gt.0.0d0) then
1260 if (idomain(k,i).eq.0) then
1261 c sigma_theta(k,i)=0.0
1264 thetatpl(k,i)=thetaref(i)
1265 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1266 & rescore(k,i-2))/3.0
1267 if (sigma_theta(k,i).ne.0)
1268 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1272 if (waga_d.gt.0.0d0) then
1274 if (itype(i).eq.10) cycle
1275 if (idomain(k,i).eq.0 ) then
1282 sigma_d(k,i)=rescore(k,i)
1283 if (sigma_d(k,i).ne.0)
1284 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1285 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1291 c remove distance restraints not used in any model from the list
1292 c shift data in all arrays
1294 if (waga_dist.ne.0.0d0) then
1300 if (ii_in_use(ii).eq.0.and.liiflag) then
1304 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1305 & .not.liiflag.and.ii.eq.lim_odl) then
1306 if (ii.eq.lim_odl) then
1307 iishift=ii-iistart+1
1312 do ki=iistart,lim_odl-iishift
1313 ires_homo(ki)=ires_homo(ki+iishift)
1314 jres_homo(ki)=jres_homo(ki+iishift)
1315 ii_in_use(ki)=ii_in_use(ki+iishift)
1316 do k=1,constr_homology
1317 odl(k,ki)=odl(k,ki+iishift)
1318 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1319 l_homo(k,ki)=l_homo(k,ki+iishift)
1323 lim_odl=lim_odl-iishift
1330 10 stop "Error infragment file"
1332 c----------------------------------------------------------------------