5 implicit real*8 (a-h,o-z)
7 include 'DIMENSIONS.ZSCOPT'
8 include 'COMMON.IOUNITS'
11 c include 'include_unres/COMMON.VAR'
12 include 'COMMON.INTERACT'
13 include 'COMMON.LOCAL'
14 include 'COMMON.NAMES'
15 include 'COMMON.CHAIN'
16 include 'COMMON.FFIELD'
17 include 'COMMON.SBRIDGE'
18 include 'COMMON.TORCNSTR'
19 include 'COMMON.CONTROL'
20 character*4 sequence(maxres)
22 double precision x(maxvar)
23 character*320 controlcard,ucase
24 dimension itype_pdb(maxres)
26 call card_concat(controlcard,.true.)
27 call reada(controlcard,'SCAL14',scal14,0.4d0)
28 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
29 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
30 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
32 call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
33 call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
34 call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
35 call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
36 write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
37 & " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
38 r0_corr=cutoff_corr-delt_corr
39 call readi(controlcard,"NRES",nres,0)
40 iscode=index(controlcard,"ONE_LETTER")
42 write (iout,*) "Error: no residues in molecule"
45 if (nres.gt.maxres) then
46 write (iout,*) "Error: too many residues",nres,maxres
48 write(iout,*) 'nres=',nres
49 C Read sequence of the protein
51 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
53 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
55 C Convert sequence to numeric code
57 itype(i)=rescode(i,sequence(i),iscode)
59 write (iout,*) "Numeric code:"
60 write (iout,'(20i4)') (itype(i),i=1,nres)
63 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
65 if (itype(i).eq.21) then
69 else if (itype(i+1).ne.20) then
71 else if (itype(i).ne.20) then
80 if (with_dihed_constr) then
82 read (inp,*) ndih_constr
83 if (ndih_constr.gt.0) then
85 write (iout,*) 'FTORS',ftors
86 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
88 & 'There are',ndih_constr,' constraints on phi angles.'
90 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
93 phi0(i)=deg2rad*phi0(i)
94 drange(i)=deg2rad*drange(i)
102 if (itype(1).eq.21) nnt=2
103 if (itype(nres).eq.21) nct=nct-1
104 write(iout,*) 'NNT=',NNT,' NCT=',NCT
106 C Juyong:READ init_vars
107 C Initialize variables!
108 C Juyong:READ read_info
109 C READ fragment information!!
110 C both routines should be in dfa.F file!!
112 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
113 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
114 write (iout,*) "Calling init_dfa_vars"
117 write (iout,*) 'init_dfa_vars finished!'
120 write (iout,*) 'read_dfa_info finished!'
124 c Read distance restraints
125 if (constr_dist.gt.0) then
126 if (refstr) call read_ref_structure(*11)
127 call read_dist_constr
131 if (constr_homology.gt.0) then
132 c write (iout,*) "About to call read_constr_homology"
134 call read_constr_homology
135 c write (iout,*) "Exit read_constr_homology"
137 if (indpdb.gt.0 .or. pdbref) then
141 cref(j,i)=crefjlee(j,i)
146 write (iout,*) "Array C"
148 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
149 & (c(j,i+nres),j=1,3)
151 write (iout,*) "Array Cref"
153 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
154 & (cref(j,i+nres),j=1,3)
158 call int_from_cart1(.false.)
159 call sc_loc_geom(.false.)
163 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
167 dc(j,i)=c(j,i+1)-c(j,i)
168 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
173 dc(j,i+nres)=c(j,i+nres)-c(j,i)
174 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
186 write (iout,'(/a,i3,a)') 'The chain contains',ns,
187 & ' disulfide-bridging cysteines.'
188 write (iout,'(20i4)') (iss(i),i=1,ns)
189 write (iout,'(/a/)') 'Pre-formed links are:'
195 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
196 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
197 & dhpb(i),ebr,forcon(i)
202 11 stop "Error reading reference structure"
204 c-----------------------------------------------------------------------------
205 logical function seq_comp(itypea,itypeb,length)
207 integer length,itypea(length),itypeb(length)
210 if (itypea(i).ne.itypeb(i)) then
218 c-----------------------------------------------------------------------------
219 subroutine read_bridge
220 C Read information about disulfide bridges.
221 implicit real*8 (a-h,o-z)
223 include 'DIMENSIONS.ZSCOPT'
224 include 'COMMON.IOUNITS'
227 include 'COMMON.INTERACT'
228 include 'COMMON.NAMES'
229 include 'COMMON.CHAIN'
230 include 'COMMON.FFIELD'
231 include 'COMMON.SBRIDGE'
232 C Read bridging residues.
233 read (inp,*) ns,(iss(i),i=1,ns)
235 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
236 C Check whether the specified bridging residues are cystines.
238 if (itype(iss(i)).ne.1) then
239 write (iout,'(2a,i3,a)')
240 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
241 & ' can form a disulfide bridge?!!!'
242 write (*,'(2a,i3,a)')
243 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
244 & ' can form a disulfide bridge?!!!'
248 C Read preformed bridges.
250 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
251 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
254 C Check if the residues involved in bridges are in the specified list of
258 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
259 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
260 write (iout,'(a,i3,a)') 'Disulfide pair',i,
261 & ' contains residues present in other pairs.'
262 write (*,'(a,i3,a)') 'Disulfide pair',i,
263 & ' contains residues present in other pairs.'
268 if (ihpb(i).eq.iss(j)) goto 10
270 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
273 if (jhpb(i).eq.iss(j)) goto 20
275 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
286 if (ns.gt.0.and.dyn_ss) then
287 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
290 C /06/28/2013 Adasko: nss number of full SS bonds
293 forcon(i-nss)=forcon(i)
300 dyn_ss_mask(iss(i))=.true.
301 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
302 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
307 c------------------------------------------------------------------------------
308 subroutine read_angles(kanal,iscor,energ,iprot,*)
309 implicit real*8 (a-h,o-z)
311 include 'DIMENSIONS.ZSCOPT'
312 include 'COMMON.INTERACT'
313 include 'COMMON.SBRIDGE'
316 include 'COMMON.CHAIN'
317 include 'COMMON.IOUNITS'
319 read(kanal,'(a80)',end=10,err=10) lineh
320 read(lineh(:5),*,err=8) ic
321 read(lineh(6:),*,err=8) energ
324 print *,'error, assuming e=1d10',lineh
328 read(lineh(18:),*,end=10,err=10) nss
330 read (lineh(20:),*,end=10,err=10)
331 & (IHPB(I),JHPB(I),I=1,NSS),iscor
333 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
334 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
337 c print *,"energy",energ," iscor",iscor
338 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
339 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
340 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
341 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
343 theta(i)=deg2rad*theta(i)
344 phi(i)=deg2rad*phi(i)
345 alph(i)=deg2rad*alph(i)
346 omeg(i)=deg2rad*omeg(i)
351 c-------------------------------------------------------------------------------
352 subroutine read_dist_constr
353 implicit real*8 (a-h,o-z)
355 include 'DIMENSIONS.ZSCOPT'
356 include 'COMMON.CONTROL'
357 include 'COMMON.CHAIN'
358 include 'COMMON.IOUNITS'
359 include 'COMMON.SBRIDGE'
360 integer ifrag_(2,100),ipair_(2,100)
361 double precision wfrag_(100),wpair_(100)
362 character*500 controlcard
363 c write (iout,*) "Calling read_dist_constr"
364 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
366 call card_concat(controlcard,.true.)
367 call readi(controlcard,"NFRAG",nfrag_,0)
368 call readi(controlcard,"NPAIR",npair_,0)
369 call readi(controlcard,"NDIST",ndist_,0)
370 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
371 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
372 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
373 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
374 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
375 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
376 write (iout,*) "IFRAG"
378 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
380 write (iout,*) "IPAIR"
382 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
385 if (.not.refstr .and. nfrag_.gt.0) then
387 & "ERROR: no reference structure to compute distance restraints"
389 & "Restraints must be specified explicitly (NDIST=number)"
392 if (nfrag_.lt.2 .and. npair_.gt.0) then
393 write (iout,*) "ERROR: Less than 2 fragments specified",
394 & " but distance restraints between pairs requested"
399 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
400 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
401 & ifrag_(2,i)=nstart_sup+nsup-1
402 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
404 if (wfrag_(i).gt.0.0d0) then
405 do j=ifrag_(1,i),ifrag_(2,i)-1
407 write (iout,*) "j",j," k",k
409 if (constr_dist.eq.1) then
414 forcon(nhpb)=wfrag_(i)
415 else if (constr_dist.eq.2) then
416 if (ddjk.le.dist_cut) then
421 forcon(nhpb)=wfrag_(i)
428 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
430 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
431 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
437 if (wpair_(i).gt.0.0d0) then
445 do j=ifrag_(1,ii),ifrag_(2,ii)
446 do k=ifrag_(1,jj),ifrag_(2,jj)
450 forcon(nhpb)=wpair_(i)
452 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
453 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
459 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
460 & ibecarb(i),forcon(nhpb+1)
461 if (forcon(nhpb+1).gt.0.0d0) then
463 if (ibecarb(i).gt.0) then
467 if (dhpb(nhpb).eq.0.0d0)
468 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
472 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
473 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
481 c====-------------------------------------------------------------------
482 subroutine read_constr_homology
485 include 'DIMENSIONS.ZSCOPT'
489 include 'COMMON.SETUP'
490 include 'COMMON.CONTROL'
491 include 'COMMON.CHAIN'
492 include 'COMMON.IOUNITS'
494 include 'COMMON.INTERACT'
495 include 'COMMON.HOMRESTR'
500 c include 'include_unres/COMMON.VAR'
503 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
505 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
506 c & sigma_odl_temp(maxres,maxres,max_template)
508 character*24 model_ki_dist, model_ki_angle
509 character*500 controlcard
510 integer ki, i, j, k, l
511 logical lprn /.true./
514 c FP - Nov. 2014 Temporary specifications for new vars
516 double precision rescore_tmp,x12,y12,z12
517 double precision, dimension (max_template,maxres) :: rescore
518 character*24 tpl_k_rescore
519 c -----------------------------------------------------------------
520 c Reading multiple PDB ref structures and calculation of retraints
521 c not using pre-computed ones stored in files model_ki_{dist,angle}
523 c -----------------------------------------------------------------
526 c Alternative: reading from input
527 call card_concat(controlcard,.true.)
528 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
529 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
530 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
531 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
532 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
534 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
535 if (homol_nset.gt.1)then
536 call card_concat(controlcard,.true.)
537 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
538 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
539 write(iout,*) "iset homology_weight "
541 c write(iout,*) i,waga_homology(i)
544 iset=mod(kolor,homol_nset)+1
549 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
551 cd write (iout,*) "nnt",nnt," nct",nct
563 c Reading HM global scores (prob not required)
565 c open (4,file="HMscore")
566 c do k=1,constr_homology
567 c read (4,*,end=521) hmscore_tmp
568 c hmscore(k)=hmscore_tmp ! Another transformation can be used
569 c write(*,*) "Model", k, ":", hmscore(k)
573 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
575 do k=1,constr_homology
577 read(inp,'(a)') pdbfile
578 c Next stament causes error upon compilation (?)
579 c if(me.eq.king.or. .not. out1file)
580 c write (iout,'(2a)') 'PDB data will be read from file ',
581 c & pdbfile(:ilen(pdbfile))
582 open(ipdbin,file=pdbfile,status='old',err=33)
584 33 write (iout,'(a)') 'Error opening PDB file.'
587 c print *,'Begin reading pdb data'
589 c Files containing res sim or local scores (former containing sigmas)
592 write(kic2,'(bz,i2.2)') k
594 tpl_k_rescore="template"//kic2//".sco"
595 c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
596 c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
597 c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
598 c tpl_k_sigma_d="template"//kic2//".sigma_d"
609 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
610 & (crefjlee(j,i+nres),j=1,3)
613 write (iout,*) "read_constr_homology: after reading pdb file"
617 c Distance restraints
620 C Copy the coordinates from reference coordinates (?)
624 c write (iout,*) "c(",j,i,") =",c(j,i)
628 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
630 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
631 open (ientin,file=tpl_k_rescore,status='old')
632 do irec=1,maxdim ! loop for reading res sim
634 rescore(k,irec)=0.0d0
637 read (ientin,*,end=1401) rescore_tmp
638 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
639 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
640 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
645 c open (ientin,file=tpl_k_sigma_odl,status='old')
646 c do irec=1,maxdim ! loop for reading sigma_odl
647 c read (ientin,*,end=1401) i, j,
648 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
649 c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
650 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
654 if (waga_dist.ne.0.0d0) then
656 do i = nnt,nct-2 ! right? without parallel.
657 do j=i+2,nct ! right?
658 c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
659 c do j=i+2,nres ! ibid
660 c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
661 c do j=i+2,nct ! ibid
663 c write (iout,*) "k",k
664 c write (iout,*) "i",i," j",j," constr_homology",
669 c Attempt to replace dist(i,j) by its definition in ...
674 distal=dsqrt(x12*x12+y12*y12+z12*z12)
677 c odl(k,ii)=dist(i,j)
678 c write (iout,*) "dist(",i,j,") =",dist(i,j)
679 c write (iout,*) "distal = ",distal
680 c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
681 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
682 c & "rescore(",k,j,") =",rescore(k,j)
684 c Calculation of sigma from res sim
686 c if (odl(k,ii).le.6.0d0) then
687 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
688 c Other functional forms possible depending on odl(k,ii), eg.
690 if (odl(k,ii).le.dist_cut) then
691 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
692 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
695 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
696 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
698 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
699 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
701 c Following expr replaced by a positive exp argument
702 c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
703 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
705 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
706 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
709 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
710 c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
712 c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
713 c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
715 c read (ientin,*) sigma_odl(k,ii) ! 1st variant
718 c if (constr_homology.gt.0) call homology_partition
721 c Theta, dihedral and SC retraints
723 if (waga_angle.gt.0.0d0) then
724 c open (ientin,file=tpl_k_sigma_dih,status='old')
725 c do irec=1,maxres-3 ! loop for reading sigma_dih
726 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
727 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
728 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
729 c & sigma_dih(k,i+nnt-1)
733 do i = nnt+3,nct ! right? without parallel.
734 c do i=1,nres ! alternative for bounds acc to readpdb?
735 c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
736 c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
737 dih(k,i)=phiref(i) ! right?
738 c read (ientin,*) sigma_dih(k,i) ! original variant
739 c write (iout,*) "dih(",k,i,") =",dih(k,i)
740 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
741 c & "rescore(",k,i-1,") =",rescore(k,i-1),
742 c & "rescore(",k,i-2,") =",rescore(k,i-2),
743 c & "rescore(",k,i-3,") =",rescore(k,i-3)
745 sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
746 & rescore(k,i-2)+rescore(k,i-3) ! right expression ?
748 c write (iout,*) "Raw sigmas for dihedral angle restraints"
749 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
750 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
751 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
752 c Instead of res sim other local measure of b/b str reliability possible
753 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
754 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
755 if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
756 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
760 if (waga_theta.gt.0.0d0) then
761 c open (ientin,file=tpl_k_sigma_theta,status='old')
762 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
763 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
764 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
765 c & sigma_theta(k,i+nnt-1)
770 do i = nnt+2,nct ! right? without parallel.
771 c do i = i=1,nres ! alternative for bounds acc to readpdb?
772 c do i=ithet_start,ithet_end ! with FG parallel.
773 thetatpl(k,i)=thetaref(i)
774 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
775 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
776 c & "rescore(",k,i-1,") =",rescore(k,i-1),
777 c & "rescore(",k,i-2,") =",rescore(k,i-2)
778 c read (ientin,*) sigma_theta(k,i) ! 1st variant
779 sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
780 & rescore(k,i-2) ! right expression ?
781 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
783 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
784 c rescore(k,i-2) ! right expression ?
785 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
786 if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
790 if (waga_d.gt.0.0d0) then
791 c open (ientin,file=tpl_k_sigma_d,status='old')
792 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
793 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
794 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
795 c & sigma_d(k,i+nnt-1)
800 do i = nnt,nct ! right? without parallel.
801 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
802 c do i=loc_start,loc_end ! with FG parallel.
803 if (itype(i).eq.10) goto 1 ! right?
807 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
808 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
809 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
810 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
811 sigma_d(k,i)=rescore(k,i) ! right expression ?
812 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
814 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
815 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
816 c read (ientin,*) sigma_d(k,i) ! 1st variant
817 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
823 if (waga_dist.ne.0.0d0) lim_odl=ii
824 if (constr_homology.gt.0) call homology_partition
825 if (constr_homology.gt.0) call init_int_table
826 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
827 cd & "lim_xx=",lim_xx
828 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
829 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
833 if (.not.lprn) return
834 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
835 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
836 write (iout,*) "Distance restraints from templates"
838 write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
839 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
841 write (iout,*) "Dihedral angle restraints from templates"
843 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
844 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
846 write (iout,*) "Virtual-bond angle restraints from templates"
848 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
849 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
851 write (iout,*) "SC restraints from templates"
853 write(iout,'(i5,10(4f8.2,4x))') i,
854 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
855 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
858 c -----------------------------------------------------------------
861 c----------------------------------------------------------------------