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
527 logical lprn /.true./
530 c FP - Nov. 2014 Temporary specifications for new vars
532 double precision rescore_tmp,x12,y12,z12
533 double precision, dimension (max_template,maxres) :: rescore
534 character*24 tpl_k_rescore
535 c -----------------------------------------------------------------
536 c Reading multiple PDB ref structures and calculation of retraints
537 c not using pre-computed ones stored in files model_ki_{dist,angle}
539 c -----------------------------------------------------------------
542 c Alternative: reading from input
543 call card_concat(controlcard,.true.)
544 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
545 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
546 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
547 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
548 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
550 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
551 call readi(controlcard,"IHSET",ihset,1)
552 if (homol_nset.gt.1)then
553 call card_concat(controlcard,.true.)
554 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
555 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
556 write(iout,*) "iset homology_weight "
558 c write(iout,*) i,waga_homology(i)
561 iset=mod(kolor,homol_nset)+1
566 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
568 cd write (iout,*) "nnt",nnt," nct",nct
580 c Reading HM global scores (prob not required)
582 c open (4,file="HMscore")
583 c do k=1,constr_homology
584 c read (4,*,end=521) hmscore_tmp
585 c hmscore(k)=hmscore_tmp ! Another transformation can be used
586 c write(*,*) "Model", k, ":", hmscore(k)
590 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
592 do k=1,constr_homology
594 read(inp,'(a)') pdbfile
595 c Next stament causes error upon compilation (?)
596 c if(me.eq.king.or. .not. out1file)
597 c write (iout,'(2a)') 'PDB data will be read from file ',
598 c & pdbfile(:ilen(pdbfile))
599 open(ipdbin,file=pdbfile,status='old',err=33)
601 33 write (iout,'(a)') 'Error opening PDB file.'
604 c print *,'Begin reading pdb data'
606 c Files containing res sim or local scores (former containing sigmas)
609 write(kic2,'(bz,i2.2)') k
611 tpl_k_rescore="template"//kic2//".sco"
612 c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
613 c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
614 c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
615 c tpl_k_sigma_d="template"//kic2//".sigma_d"
626 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
627 & (crefjlee(j,i+nres),j=1,3)
630 write (iout,*) "read_constr_homology: after reading pdb file"
634 c Distance restraints
637 C Copy the coordinates from reference coordinates (?)
641 c write (iout,*) "c(",j,i,") =",c(j,i)
645 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
647 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
648 open (ientin,file=tpl_k_rescore,status='old')
649 do irec=1,maxdim ! loop for reading res sim
651 rescore(k,irec)=0.0d0
654 read (ientin,*,end=1401) rescore_tmp
655 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
656 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
657 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
662 c open (ientin,file=tpl_k_sigma_odl,status='old')
663 c do irec=1,maxdim ! loop for reading sigma_odl
664 c read (ientin,*,end=1401) i, j,
665 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
666 c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
667 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
671 if (waga_dist.ne.0.0d0) then
673 do i = nnt,nct-2 ! right? without parallel.
674 do j=i+2,nct ! right?
675 c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
676 c do j=i+2,nres ! ibid
677 c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
678 c do j=i+2,nct ! ibid
680 c write (iout,*) "k",k
681 c write (iout,*) "i",i," j",j," constr_homology",
686 c Attempt to replace dist(i,j) by its definition in ...
691 distal=dsqrt(x12*x12+y12*y12+z12*z12)
694 c odl(k,ii)=dist(i,j)
695 c write (iout,*) "dist(",i,j,") =",dist(i,j)
696 c write (iout,*) "distal = ",distal
697 c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
698 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
699 c & "rescore(",k,j,") =",rescore(k,j)
701 c Calculation of sigma from res sim
703 c if (odl(k,ii).le.6.0d0) then
704 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
705 c Other functional forms possible depending on odl(k,ii), eg.
707 if (odl(k,ii).le.dist_cut) then
708 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
709 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
710 c write (iout,*) "c sigma_odl",k,ii,sigma_odl(k,ii)
713 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
714 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
716 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
717 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
718 c write (iout,*) "d sigma_odl",k,ii,sigma_odl(k,ii),
719 c & odl(k,ii),dist_cut
721 c Following expr replaced by a positive exp argument
722 c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
723 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
725 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
726 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
729 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
730 c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
732 c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
733 c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
735 c read (ientin,*) sigma_odl(k,ii) ! 1st variant
738 c if (constr_homology.gt.0) call homology_partition
741 c Theta, dihedral and SC retraints
743 if (waga_angle.gt.0.0d0) then
744 c open (ientin,file=tpl_k_sigma_dih,status='old')
745 c do irec=1,maxres-3 ! loop for reading sigma_dih
746 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
747 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
748 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
749 c & sigma_dih(k,i+nnt-1)
753 do i = nnt+3,nct ! right? without parallel.
754 c do i=1,nres ! alternative for bounds acc to readpdb?
755 c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
756 c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
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) ! right expression ?
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)
775 if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
776 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
780 if (waga_theta.gt.0.0d0) then
781 c open (ientin,file=tpl_k_sigma_theta,status='old')
782 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
783 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
784 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
785 c & sigma_theta(k,i+nnt-1)
790 do i = nnt+2,nct ! right? without parallel.
791 c do i = i=1,nres ! alternative for bounds acc to readpdb?
792 c do i=ithet_start,ithet_end ! with FG parallel.
793 thetatpl(k,i)=thetaref(i)
794 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
795 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
796 c & "rescore(",k,i-1,") =",rescore(k,i-1),
797 c & "rescore(",k,i-2,") =",rescore(k,i-2)
798 c read (ientin,*) sigma_theta(k,i) ! 1st variant
799 sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
800 & rescore(k,i-2) ! right expression ?
801 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
803 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
804 c rescore(k,i-2) ! right expression ?
805 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
806 if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
810 if (waga_d.gt.0.0d0) then
811 c open (ientin,file=tpl_k_sigma_d,status='old')
812 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
813 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
814 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
815 c & sigma_d(k,i+nnt-1)
820 do i = nnt,nct ! right? without parallel.
821 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
822 c do i=loc_start,loc_end ! with FG parallel.
823 if (itype(i).eq.10) goto 1 ! right?
827 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
828 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
829 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
830 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
831 sigma_d(k,i)=rescore(k,i) ! right expression ?
832 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
834 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
835 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
836 c read (ientin,*) sigma_d(k,i) ! 1st variant
837 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
843 if (waga_dist.ne.0.0d0) lim_odl=ii
844 if (constr_homology.gt.0) call homology_partition
845 if (constr_homology.gt.0) call init_int_table
846 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
847 cd & "lim_xx=",lim_xx
848 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
849 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
853 if (.not.lprn) return
854 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
855 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
856 write (iout,*) "Distance restraints from templates"
858 write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
859 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
861 write (iout,*) "Dihedral angle restraints from templates"
863 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
864 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
866 write (iout,*) "Virtual-bond angle restraints from templates"
868 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
869 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
871 write (iout,*) "SC restraints from templates"
873 write(iout,'(i5,10(4f8.2,4x))') i,
874 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
875 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
883 c -----------------------------------------------------------------
886 c----------------------------------------------------------------------