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 c write (iout,*) "Calling read_dist_constr"
379 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
381 call card_concat(controlcard,.true.)
382 call readi(controlcard,"NFRAG",nfrag_,0)
383 call readi(controlcard,"NPAIR",npair_,0)
384 call readi(controlcard,"NDIST",ndist_,0)
385 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
386 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
387 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
388 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
389 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
390 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
391 write (iout,*) "IFRAG"
393 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
395 write (iout,*) "IPAIR"
397 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
400 if (.not.refstr .and. nfrag_.gt.0) then
402 & "ERROR: no reference structure to compute distance restraints"
404 & "Restraints must be specified explicitly (NDIST=number)"
407 if (nfrag_.lt.2 .and. npair_.gt.0) then
408 write (iout,*) "ERROR: Less than 2 fragments specified",
409 & " but distance restraints between pairs requested"
414 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
415 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
416 & ifrag_(2,i)=nstart_sup+nsup-1
417 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
419 if (wfrag_(i).gt.0.0d0) then
420 do j=ifrag_(1,i),ifrag_(2,i)-1
422 write (iout,*) "j",j," k",k
424 if (constr_dist.eq.1) then
429 forcon(nhpb)=wfrag_(i)
430 else if (constr_dist.eq.2) then
431 if (ddjk.le.dist_cut) then
436 forcon(nhpb)=wfrag_(i)
443 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
445 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
446 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
452 if (wpair_(i).gt.0.0d0) then
460 do j=ifrag_(1,ii),ifrag_(2,ii)
461 do k=ifrag_(1,jj),ifrag_(2,jj)
465 forcon(nhpb)=wpair_(i)
467 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
468 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
474 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
475 & ibecarb(i),forcon(nhpb+1)
476 if (forcon(nhpb+1).gt.0.0d0) then
478 if (ibecarb(i).gt.0) then
482 if (dhpb(nhpb).eq.0.0d0)
483 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
487 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
488 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
496 c====-------------------------------------------------------------------
497 subroutine read_constr_homology
500 include 'DIMENSIONS.ZSCOPT'
501 include 'DIMENSIONS.FREE'
505 include 'COMMON.SETUP'
506 include 'COMMON.CONTROL'
507 include 'COMMON.CHAIN'
508 include 'COMMON.IOUNITS'
510 include 'COMMON.INTERACT'
511 include 'COMMON.HOMRESTR'
516 c include 'include_unres/COMMON.VAR'
519 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
521 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
522 c & sigma_odl_temp(maxres,maxres,max_template)
524 character*24 model_ki_dist, model_ki_angle
525 character*500 controlcard
526 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
527 integer idomain(max_template,maxres)
528 logical lprn /.true./
533 c FP - Nov. 2014 Temporary specifications for new vars
535 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
536 double precision, dimension (max_template,maxres) :: rescore
537 double precision, dimension (max_template,maxres) :: rescore2
538 character*24 tpl_k_rescore
539 c -----------------------------------------------------------------
540 c Reading multiple PDB ref structures and calculation of retraints
541 c not using pre-computed ones stored in files model_ki_{dist,angle}
543 c -----------------------------------------------------------------
546 c Alternative: reading from input
547 call card_concat(controlcard,.true.)
548 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
549 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
550 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
551 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
552 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
553 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
554 call readi(controlcard,"HOMOL_SET",homol_nset,1)
555 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
556 call readi(controlcard,"IHSET",ihset,1)
557 if (homol_nset.gt.1)then
558 call card_concat(controlcard,.true.)
559 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
560 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
561 write(iout,*) "iset homology_weight "
563 c write(iout,*) i,waga_homology(i)
566 iset=mod(kolor,homol_nset)+1
571 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
573 cd write (iout,*) "nnt",nnt," nct",nct
585 c Reading HM global scores (prob not required)
588 do k=1,constr_homology
592 c open (4,file="HMscore")
593 c do k=1,constr_homology
594 c read (4,*,end=521) hmscore_tmp
595 c hmscore(k)=hmscore_tmp ! Another transformation can be used
596 c write(*,*) "Model", k, ":", hmscore(k)
607 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
609 do k=1,constr_homology
611 read(inp,'(a)') pdbfile
612 c Next stament causes error upon compilation (?)
613 c if(me.eq.king.or. .not. out1file)
614 c write (iout,'(2a)') 'PDB data will be read from file ',
615 c & pdbfile(:ilen(pdbfile))
616 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
617 & pdbfile(:ilen(pdbfile))
618 open(ipdbin,file=pdbfile,status='old',err=33)
620 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
621 & pdbfile(:ilen(pdbfile))
624 c print *,'Begin reading pdb data'
626 c Files containing res sim or local scores (former containing sigmas)
629 write(kic2,'(bz,i2.2)') k
631 tpl_k_rescore="template"//kic2//".sco"
642 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
643 & (crefjlee(j,i+nres),j=1,3)
646 write (iout,*) "read_constr_homology: after reading pdb file"
650 c Distance restraints
653 C Copy the coordinates from reference coordinates (?)
657 c write (iout,*) "c(",j,i,") =",c(j,i)
661 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
663 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
664 open (ientin,file=tpl_k_rescore,status='old')
665 if (nnt.gt.1) rescore(k,1)=0.0d0
666 do irec=nnt,maxdim ! loop for reading res sim
668 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
671 idomain(k,i_tmp)=idomain_tmp
672 rescore(k,i_tmp)=rescore_tmp
673 rescore2(k,i_tmp)=rescore2_tmp
676 read (ientin,*,end=1401) rescore_tmp
678 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
679 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
680 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
685 if (waga_dist.ne.0.0d0) then
693 distal=dsqrt(x12*x12+y12*y12+z12*z12)
694 c write (iout,*) k,i,j,distal,dist2_cut
696 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
697 & .and. distal.le.dist2_cut ) then
703 c write (iout,*) "k",k
704 c write (iout,*) "i",i," j",j," constr_homology",
712 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
714 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
715 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
716 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
718 if (odl(k,ii).le.dist_cut) then
719 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
722 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
723 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
725 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
726 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
730 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
740 c Theta, dihedral and SC retraints
742 if (waga_angle.gt.0.0d0) then
743 c open (ientin,file=tpl_k_sigma_dih,status='old')
744 c do irec=1,maxres-3 ! loop for reading sigma_dih
745 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
746 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
747 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
748 c & sigma_dih(k,i+nnt-1)
753 if (idomain(k,i).eq.0) then
757 dih(k,i)=phiref(i) ! right?
758 c read (ientin,*) sigma_dih(k,i) ! original variant
759 c write (iout,*) "dih(",k,i,") =",dih(k,i)
760 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
761 c & "rescore(",k,i-1,") =",rescore(k,i-1),
762 c & "rescore(",k,i-2,") =",rescore(k,i-2),
763 c & "rescore(",k,i-3,") =",rescore(k,i-3)
765 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
766 & rescore(k,i-2)+rescore(k,i-3))/4.0
767 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
768 c write (iout,*) "Raw sigmas for dihedral angle restraints"
769 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
770 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
771 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
772 c Instead of res sim other local measure of b/b str reliability possible
773 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
774 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
779 if (waga_theta.gt.0.0d0) then
780 c open (ientin,file=tpl_k_sigma_theta,status='old')
781 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
782 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
783 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
784 c & sigma_theta(k,i+nnt-1)
789 do i = nnt+2,nct ! right? without parallel.
790 c do i = i=1,nres ! alternative for bounds acc to readpdb?
791 c do i=ithet_start,ithet_end ! with FG parallel.
792 if (idomain(k,i).eq.0) then
796 thetatpl(k,i)=thetaref(i)
797 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
798 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
799 c & "rescore(",k,i-1,") =",rescore(k,i-1),
800 c & "rescore(",k,i-2,") =",rescore(k,i-2)
801 c read (ientin,*) sigma_theta(k,i) ! 1st variant
802 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
803 & rescore(k,i-2))/3.0
804 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
805 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
807 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
808 c rescore(k,i-2) ! right expression ?
809 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
814 if (waga_d.gt.0.0d0) then
815 c open (ientin,file=tpl_k_sigma_d,status='old')
816 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
817 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
818 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
819 c & sigma_d(k,i+nnt-1)
823 do i = nnt,nct ! right? without parallel.
824 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
825 c do i=loc_start,loc_end ! with FG parallel.
826 if (itype(i).eq.10) cycle
827 if (idomain(k,i).eq.0 ) then
834 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
835 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
836 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
837 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
838 sigma_d(k,i)=rescore(k,i) ! right expression ?
839 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
841 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
842 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
843 c read (ientin,*) sigma_d(k,i) ! 1st variant
844 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
850 c remove distance restraints not used in any model from the list
851 c shift data in all arrays
853 if (waga_dist.ne.0.0d0) then
858 if (ii_in_use(ii).eq.0) then
860 ires_homo(ki)=ires_homo(ki+1)
861 jres_homo(ki)=jres_homo(ki+1)
862 ii_in_use(ki)=ii_in_use(ki+1)
863 do k=1,constr_homology
864 odl(k,ki)=odl(k,ki+1)
865 sigma_odl(k,ki)=sigma_odl(k,ki+1)
866 l_homo(k,ki)=l_homo(k,ki+1)
875 if (constr_homology.gt.0) call homology_partition
876 if (constr_homology.gt.0) call init_int_table
877 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
878 cd & "lim_xx=",lim_xx
879 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
880 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
884 if (.not.lprn) return
885 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
886 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
887 write (iout,*) "Distance restraints from templates"
889 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
890 & ii,ires_homo(ii),jres_homo(ii),
891 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
892 & ki=1,constr_homology)
894 write (iout,*) "Dihedral angle restraints from templates"
896 write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
897 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
899 write (iout,*) "Virtual-bond angle restraints from templates"
901 write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
902 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
904 write (iout,*) "SC restraints from templates"
906 write(iout,'(i5,100(4f8.2,4x))') i,
907 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
908 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
911 c -----------------------------------------------------------------
914 c----------------------------------------------------------------------