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./
513 c FP - Nov. 2014 Temporary specifications for new vars
515 double precision rescore_tmp,x12,y12,z12
516 double precision, dimension (max_template,maxres) :: rescore
517 character*24 tpl_k_rescore
518 c -----------------------------------------------------------------
519 c Reading multiple PDB ref structures and calculation of retraints
520 c not using pre-computed ones stored in files model_ki_{dist,angle}
522 c -----------------------------------------------------------------
525 c Alternative: reading from input
526 call card_concat(controlcard,.true.)
527 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
528 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
529 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
530 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
531 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
533 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
534 if (homol_nset.gt.1)then
535 call card_concat(controlcard,.true.)
536 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
537 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
538 write(iout,*) "iset homology_weight "
540 c write(iout,*) i,waga_homology(i)
543 iset=mod(kolor,homol_nset)+1
548 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
550 cd write (iout,*) "nnt",nnt," nct",nct
562 c Reading HM global scores (prob not required)
564 c open (4,file="HMscore")
565 c do k=1,constr_homology
566 c read (4,*,end=521) hmscore_tmp
567 c hmscore(k)=hmscore_tmp ! Another transformation can be used
568 c write(*,*) "Model", k, ":", hmscore(k)
572 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
574 do k=1,constr_homology
576 read(inp,'(a)') pdbfile
577 c Next stament causes error upon compilation (?)
578 c if(me.eq.king.or. .not. out1file)
579 c write (iout,'(2a)') 'PDB data will be read from file ',
580 c & pdbfile(:ilen(pdbfile))
581 open(ipdbin,file=pdbfile,status='old',err=33)
583 33 write (iout,'(a)') 'Error opening PDB file.'
586 c print *,'Begin reading pdb data'
588 c Files containing res sim or local scores (former containing sigmas)
591 write(kic2,'(bz,i2.2)') k
593 tpl_k_rescore="template"//kic2//".sco"
594 c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
595 c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
596 c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
597 c tpl_k_sigma_d="template"//kic2//".sigma_d"
608 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
609 & (crefjlee(j,i+nres),j=1,3)
612 write (iout,*) "read_constr_homology: after reading pdb file"
616 c Distance restraints
619 C Copy the coordinates from reference coordinates (?)
623 c write (iout,*) "c(",j,i,") =",c(j,i)
627 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
629 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
630 open (ientin,file=tpl_k_rescore,status='old')
631 do irec=1,maxdim ! loop for reading res sim
633 rescore(k,irec)=0.0d0
636 read (ientin,*,end=1401) rescore_tmp
637 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
638 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
639 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
644 c open (ientin,file=tpl_k_sigma_odl,status='old')
645 c do irec=1,maxdim ! loop for reading sigma_odl
646 c read (ientin,*,end=1401) i, j,
647 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
648 c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
649 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
653 if (waga_dist.ne.0.0d0) then
655 do i = nnt,nct-2 ! right? without parallel.
656 do j=i+2,nct ! right?
657 c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
658 c do j=i+2,nres ! ibid
659 c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
660 c do j=i+2,nct ! ibid
662 c write (iout,*) "k",k
663 c write (iout,*) "i",i," j",j," constr_homology",
668 c Attempt to replace dist(i,j) by its definition in ...
673 distal=dsqrt(x12*x12+y12*y12+z12*z12)
676 c odl(k,ii)=dist(i,j)
677 c write (iout,*) "dist(",i,j,") =",dist(i,j)
678 c write (iout,*) "distal = ",distal
679 c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
680 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
681 c & "rescore(",k,j,") =",rescore(k,j)
683 c Calculation of sigma from res sim
685 c if (odl(k,ii).le.6.0d0) then
686 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
687 c Other functional forms possible depending on odl(k,ii), eg.
689 if (odl(k,ii).le.dist_cut) then
690 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
691 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
693 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
694 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
696 c Following expr replaced by a positive exp argument
697 c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
698 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
700 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
701 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
704 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
705 c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
707 c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
708 c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
710 c read (ientin,*) sigma_odl(k,ii) ! 1st variant
713 c if (constr_homology.gt.0) call homology_partition
716 c Theta, dihedral and SC retraints
718 if (waga_angle.gt.0.0d0) then
719 c open (ientin,file=tpl_k_sigma_dih,status='old')
720 c do irec=1,maxres-3 ! loop for reading sigma_dih
721 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
722 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
723 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
724 c & sigma_dih(k,i+nnt-1)
728 do i = nnt+3,nct ! right? without parallel.
729 c do i=1,nres ! alternative for bounds acc to readpdb?
730 c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
731 c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
732 dih(k,i)=phiref(i) ! right?
733 c read (ientin,*) sigma_dih(k,i) ! original variant
734 c write (iout,*) "dih(",k,i,") =",dih(k,i)
735 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
736 c & "rescore(",k,i-1,") =",rescore(k,i-1),
737 c & "rescore(",k,i-2,") =",rescore(k,i-2),
738 c & "rescore(",k,i-3,") =",rescore(k,i-3)
740 sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
741 & rescore(k,i-2)+rescore(k,i-3) ! right expression ?
743 c write (iout,*) "Raw sigmas for dihedral angle restraints"
744 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
745 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
746 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
747 c Instead of res sim other local measure of b/b str reliability possible
748 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
749 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
750 if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
751 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
755 if (waga_theta.gt.0.0d0) then
756 c open (ientin,file=tpl_k_sigma_theta,status='old')
757 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
758 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
759 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
760 c & sigma_theta(k,i+nnt-1)
765 do i = nnt+2,nct ! right? without parallel.
766 c do i = i=1,nres ! alternative for bounds acc to readpdb?
767 c do i=ithet_start,ithet_end ! with FG parallel.
768 thetatpl(k,i)=thetaref(i)
769 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
770 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
771 c & "rescore(",k,i-1,") =",rescore(k,i-1),
772 c & "rescore(",k,i-2,") =",rescore(k,i-2)
773 c read (ientin,*) sigma_theta(k,i) ! 1st variant
774 sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
775 & rescore(k,i-2) ! right expression ?
776 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
778 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
779 c rescore(k,i-2) ! right expression ?
780 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
781 if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
785 if (waga_d.gt.0.0d0) then
786 c open (ientin,file=tpl_k_sigma_d,status='old')
787 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
788 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
789 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
790 c & sigma_d(k,i+nnt-1)
795 do i = nnt,nct ! right? without parallel.
796 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
797 c do i=loc_start,loc_end ! with FG parallel.
798 if (itype(i).eq.10) goto 1 ! right?
802 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
803 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
804 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
805 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
806 sigma_d(k,i)=rescore(k,i) ! right expression ?
807 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
809 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
810 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
811 c read (ientin,*) sigma_d(k,i) ! 1st variant
812 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
818 if (waga_dist.ne.0.0d0) lim_odl=ii
819 if (constr_homology.gt.0) call homology_partition
820 if (constr_homology.gt.0) call init_int_table
821 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
822 cd & "lim_xx=",lim_xx
823 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
824 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
828 if (.not.lprn) return
829 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
830 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
831 write (iout,*) "Distance restraints from templates"
833 write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
834 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
836 write (iout,*) "Dihedral angle restraints from templates"
838 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
839 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
841 write (iout,*) "Virtual-bond angle restraints from templates"
843 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
844 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
846 write (iout,*) "SC restraints from templates"
848 write(iout,'(i5,10(4f8.2,4x))') i,
849 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
850 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
853 c -----------------------------------------------------------------
856 c----------------------------------------------------------------------