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 write (iout,*) "Numeric code:"
61 write (iout,'(20i4)') (itype(i),i=1,nres)
64 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
66 if (itype(i).eq.21) then
70 else if (itype(i+1).ne.20) then
72 else if (itype(i).ne.20) then
81 if (with_dihed_constr) then
83 read (inp,*) ndih_constr
84 if (ndih_constr.gt.0) then
86 write (iout,*) 'FTORS',ftors
87 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
89 & 'There are',ndih_constr,' constraints on phi angles.'
91 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
94 phi0(i)=deg2rad*phi0(i)
95 drange(i)=deg2rad*drange(i)
103 if (itype(1).eq.21) nnt=2
104 if (itype(nres).eq.21) nct=nct-1
105 write(iout,*) 'NNT=',NNT,' NCT=',NCT
107 C Juyong:READ init_vars
108 C Initialize variables!
109 C Juyong:READ read_info
110 C READ fragment information!!
111 C both routines should be in dfa.F file!!
113 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
114 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
115 write (iout,*) "Calling init_dfa_vars"
118 write (iout,*) 'init_dfa_vars finished!'
121 write (iout,*) 'read_dfa_info finished!'
125 c Read distance restraints
126 if (constr_dist.gt.0) then
127 if (refstr) call read_ref_structure(*11)
128 call read_dist_constr
132 if (constr_homology.gt.0) then
133 c write (iout,*) "About to call read_constr_homology"
135 call read_constr_homology
136 c write (iout,*) "Exit read_constr_homology"
138 if (indpdb.gt.0 .or. pdbref) then
142 cref(j,i)=crefjlee(j,i)
147 write (iout,*) "Array C"
149 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
150 & (c(j,i+nres),j=1,3)
152 write (iout,*) "Array Cref"
154 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
155 & (cref(j,i+nres),j=1,3)
159 call int_from_cart1(.false.)
160 call sc_loc_geom(.false.)
164 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
168 dc(j,i)=c(j,i+1)-c(j,i)
169 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
174 dc(j,i+nres)=c(j,i+nres)-c(j,i)
175 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
187 write (iout,'(/a,i3,a)') 'The chain contains',ns,
188 & ' disulfide-bridging cysteines.'
189 write (iout,'(20i4)') (iss(i),i=1,ns)
190 write (iout,'(/a/)') 'Pre-formed links are:'
196 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
197 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
198 & dhpb(i),ebr,forcon(i)
203 11 stop "Error reading reference structure"
205 c-----------------------------------------------------------------------------
206 logical function seq_comp(itypea,itypeb,length)
208 integer length,itypea(length),itypeb(length)
211 if (itypea(i).ne.itypeb(i)) then
219 c-----------------------------------------------------------------------------
220 subroutine read_bridge
221 C Read information about disulfide bridges.
222 implicit real*8 (a-h,o-z)
224 include 'DIMENSIONS.ZSCOPT'
225 include 'COMMON.IOUNITS'
228 include 'COMMON.INTERACT'
229 include 'COMMON.NAMES'
230 include 'COMMON.CHAIN'
231 include 'COMMON.FFIELD'
232 include 'COMMON.SBRIDGE'
233 C Read bridging residues.
234 read (inp,*) ns,(iss(i),i=1,ns)
236 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
237 C Check whether the specified bridging residues are cystines.
239 if (itype(iss(i)).ne.1) then
240 write (iout,'(2a,i3,a)')
241 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
242 & ' can form a disulfide bridge?!!!'
243 write (*,'(2a,i3,a)')
244 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
245 & ' can form a disulfide bridge?!!!'
249 C Read preformed bridges.
251 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
252 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
255 C Check if the residues involved in bridges are in the specified list of
259 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
260 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
261 write (iout,'(a,i3,a)') 'Disulfide pair',i,
262 & ' contains residues present in other pairs.'
263 write (*,'(a,i3,a)') 'Disulfide pair',i,
264 & ' contains residues present in other pairs.'
269 if (ihpb(i).eq.iss(j)) goto 10
271 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
274 if (jhpb(i).eq.iss(j)) goto 20
276 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
287 if (ns.gt.0.and.dyn_ss) then
288 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
291 C /06/28/2013 Adasko: nss number of full SS bonds
294 forcon(i-nss)=forcon(i)
301 dyn_ss_mask(iss(i))=.true.
302 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
303 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
308 c------------------------------------------------------------------------------
309 subroutine read_angles(kanal,iscor,energ,iprot,*)
310 implicit real*8 (a-h,o-z)
312 include 'DIMENSIONS.ZSCOPT'
313 include 'COMMON.INTERACT'
314 include 'COMMON.SBRIDGE'
317 include 'COMMON.CHAIN'
318 include 'COMMON.IOUNITS'
320 read(kanal,'(a80)',end=10,err=10) lineh
321 read(lineh(:5),*,err=8) ic
322 read(lineh(6:),*,err=8) energ
325 print *,'error, assuming e=1d10',lineh
329 read(lineh(18:),*,end=10,err=10) nss
331 read (lineh(20:),*,end=10,err=10)
332 & (IHPB(I),JHPB(I),I=1,NSS),iscor
334 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
335 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
338 c print *,"energy",energ," iscor",iscor
339 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
340 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
341 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
342 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
344 theta(i)=deg2rad*theta(i)
345 phi(i)=deg2rad*phi(i)
346 alph(i)=deg2rad*alph(i)
347 omeg(i)=deg2rad*omeg(i)
352 c-------------------------------------------------------------------------------
353 subroutine read_dist_constr
354 implicit real*8 (a-h,o-z)
356 include 'DIMENSIONS.ZSCOPT'
357 include 'DIMENSIONS.FREE'
358 include 'COMMON.CONTROL'
359 include 'COMMON.CHAIN'
360 include 'COMMON.IOUNITS'
361 include 'COMMON.SBRIDGE'
362 integer ifrag_(2,100),ipair_(2,100)
363 double precision wfrag_(100),wpair_(100)
364 character*500 controlcard
365 c write (iout,*) "Calling read_dist_constr"
366 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
368 call card_concat(controlcard,.true.)
369 call readi(controlcard,"NFRAG",nfrag_,0)
370 call readi(controlcard,"NPAIR",npair_,0)
371 call readi(controlcard,"NDIST",ndist_,0)
372 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
373 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
374 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
375 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
376 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
377 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
378 write (iout,*) "IFRAG"
380 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
382 write (iout,*) "IPAIR"
384 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
387 if (.not.refstr .and. nfrag_.gt.0) then
389 & "ERROR: no reference structure to compute distance restraints"
391 & "Restraints must be specified explicitly (NDIST=number)"
394 if (nfrag_.lt.2 .and. npair_.gt.0) then
395 write (iout,*) "ERROR: Less than 2 fragments specified",
396 & " but distance restraints between pairs requested"
401 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
402 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
403 & ifrag_(2,i)=nstart_sup+nsup-1
404 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
406 if (wfrag_(i).gt.0.0d0) then
407 do j=ifrag_(1,i),ifrag_(2,i)-1
409 write (iout,*) "j",j," k",k
411 if (constr_dist.eq.1) then
416 forcon(nhpb)=wfrag_(i)
417 else if (constr_dist.eq.2) then
418 if (ddjk.le.dist_cut) then
423 forcon(nhpb)=wfrag_(i)
430 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
432 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
433 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
439 if (wpair_(i).gt.0.0d0) then
447 do j=ifrag_(1,ii),ifrag_(2,ii)
448 do k=ifrag_(1,jj),ifrag_(2,jj)
452 forcon(nhpb)=wpair_(i)
454 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
455 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
461 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
462 & ibecarb(i),forcon(nhpb+1)
463 if (forcon(nhpb+1).gt.0.0d0) then
465 if (ibecarb(i).gt.0) then
469 if (dhpb(nhpb).eq.0.0d0)
470 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
474 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
475 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
483 c====-------------------------------------------------------------------
484 subroutine read_constr_homology
487 include 'DIMENSIONS.ZSCOPT'
488 include 'DIMENSIONS.FREE'
492 include 'COMMON.SETUP'
493 include 'COMMON.CONTROL'
494 include 'COMMON.CHAIN'
495 include 'COMMON.IOUNITS'
497 include 'COMMON.INTERACT'
498 include 'COMMON.HOMRESTR'
503 c include 'include_unres/COMMON.VAR'
506 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
508 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
509 c & sigma_odl_temp(maxres,maxres,max_template)
511 character*24 model_ki_dist, model_ki_angle
512 character*500 controlcard
513 integer ki, i, j, k, l
514 logical lprn /.true./
517 c FP - Nov. 2014 Temporary specifications for new vars
519 double precision rescore_tmp,x12,y12,z12
520 double precision, dimension (max_template,maxres) :: rescore
521 character*24 tpl_k_rescore
522 c -----------------------------------------------------------------
523 c Reading multiple PDB ref structures and calculation of retraints
524 c not using pre-computed ones stored in files model_ki_{dist,angle}
526 c -----------------------------------------------------------------
529 c Alternative: reading from input
530 call card_concat(controlcard,.true.)
531 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
532 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
533 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
534 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
535 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
537 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
538 call readi(controlcard,"IHSET",ihset,1)
539 if (homol_nset.gt.1)then
540 call card_concat(controlcard,.true.)
541 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
542 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
543 write(iout,*) "iset homology_weight "
545 c write(iout,*) i,waga_homology(i)
548 iset=mod(kolor,homol_nset)+1
553 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
555 cd write (iout,*) "nnt",nnt," nct",nct
567 c Reading HM global scores (prob not required)
569 c open (4,file="HMscore")
570 c do k=1,constr_homology
571 c read (4,*,end=521) hmscore_tmp
572 c hmscore(k)=hmscore_tmp ! Another transformation can be used
573 c write(*,*) "Model", k, ":", hmscore(k)
577 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
579 do k=1,constr_homology
581 read(inp,'(a)') pdbfile
582 c Next stament causes error upon compilation (?)
583 c if(me.eq.king.or. .not. out1file)
584 c write (iout,'(2a)') 'PDB data will be read from file ',
585 c & pdbfile(:ilen(pdbfile))
586 open(ipdbin,file=pdbfile,status='old',err=33)
588 33 write (iout,'(a)') 'Error opening PDB file.'
591 c print *,'Begin reading pdb data'
593 c Files containing res sim or local scores (former containing sigmas)
596 write(kic2,'(bz,i2.2)') k
598 tpl_k_rescore="template"//kic2//".sco"
599 c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
600 c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
601 c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
602 c tpl_k_sigma_d="template"//kic2//".sigma_d"
613 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
614 & (crefjlee(j,i+nres),j=1,3)
617 write (iout,*) "read_constr_homology: after reading pdb file"
621 c Distance restraints
624 C Copy the coordinates from reference coordinates (?)
628 c write (iout,*) "c(",j,i,") =",c(j,i)
632 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
634 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
635 open (ientin,file=tpl_k_rescore,status='old')
636 do irec=1,maxdim ! loop for reading res sim
638 rescore(k,irec)=0.0d0
641 read (ientin,*,end=1401) rescore_tmp
642 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
643 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
644 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
649 c open (ientin,file=tpl_k_sigma_odl,status='old')
650 c do irec=1,maxdim ! loop for reading sigma_odl
651 c read (ientin,*,end=1401) i, j,
652 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
653 c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
654 c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
658 if (waga_dist.ne.0.0d0) then
660 do i = nnt,nct-2 ! right? without parallel.
661 do j=i+2,nct ! right?
662 c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
663 c do j=i+2,nres ! ibid
664 c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
665 c do j=i+2,nct ! ibid
667 c write (iout,*) "k",k
668 c write (iout,*) "i",i," j",j," constr_homology",
673 c Attempt to replace dist(i,j) by its definition in ...
678 distal=dsqrt(x12*x12+y12*y12+z12*z12)
681 c odl(k,ii)=dist(i,j)
682 c write (iout,*) "dist(",i,j,") =",dist(i,j)
683 c write (iout,*) "distal = ",distal
684 c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
685 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
686 c & "rescore(",k,j,") =",rescore(k,j)
688 c Calculation of sigma from res sim
690 c if (odl(k,ii).le.6.0d0) then
691 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
692 c Other functional forms possible depending on odl(k,ii), eg.
694 if (odl(k,ii).le.dist_cut) then
695 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
696 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
697 c write (iout,*) "c sigma_odl",k,ii,sigma_odl(k,ii)
700 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
701 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
703 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
704 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
705 c write (iout,*) "d sigma_odl",k,ii,sigma_odl(k,ii),
706 c & odl(k,ii),dist_cut
708 c Following expr replaced by a positive exp argument
709 c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
710 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
712 c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
713 c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
716 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
717 c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
719 c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
720 c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
722 c read (ientin,*) sigma_odl(k,ii) ! 1st variant
725 c if (constr_homology.gt.0) call homology_partition
728 c Theta, dihedral and SC retraints
730 if (waga_angle.gt.0.0d0) then
731 c open (ientin,file=tpl_k_sigma_dih,status='old')
732 c do irec=1,maxres-3 ! loop for reading sigma_dih
733 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
734 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
735 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
736 c & sigma_dih(k,i+nnt-1)
740 do i = nnt+3,nct ! right? without parallel.
741 c do i=1,nres ! alternative for bounds acc to readpdb?
742 c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
743 c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
744 dih(k,i)=phiref(i) ! right?
745 c read (ientin,*) sigma_dih(k,i) ! original variant
746 c write (iout,*) "dih(",k,i,") =",dih(k,i)
747 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
748 c & "rescore(",k,i-1,") =",rescore(k,i-1),
749 c & "rescore(",k,i-2,") =",rescore(k,i-2),
750 c & "rescore(",k,i-3,") =",rescore(k,i-3)
752 sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
753 & rescore(k,i-2)+rescore(k,i-3) ! right expression ?
755 c write (iout,*) "Raw sigmas for dihedral angle restraints"
756 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
757 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
758 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
759 c Instead of res sim other local measure of b/b str reliability possible
760 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
761 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
762 if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
763 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
767 if (waga_theta.gt.0.0d0) then
768 c open (ientin,file=tpl_k_sigma_theta,status='old')
769 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
770 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
771 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
772 c & sigma_theta(k,i+nnt-1)
777 do i = nnt+2,nct ! right? without parallel.
778 c do i = i=1,nres ! alternative for bounds acc to readpdb?
779 c do i=ithet_start,ithet_end ! with FG parallel.
780 thetatpl(k,i)=thetaref(i)
781 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
782 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
783 c & "rescore(",k,i-1,") =",rescore(k,i-1),
784 c & "rescore(",k,i-2,") =",rescore(k,i-2)
785 c read (ientin,*) sigma_theta(k,i) ! 1st variant
786 sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
787 & rescore(k,i-2) ! right expression ?
788 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
790 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
791 c rescore(k,i-2) ! right expression ?
792 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
793 if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
797 if (waga_d.gt.0.0d0) then
798 c open (ientin,file=tpl_k_sigma_d,status='old')
799 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
800 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
801 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
802 c & sigma_d(k,i+nnt-1)
807 do i = nnt,nct ! right? without parallel.
808 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
809 c do i=loc_start,loc_end ! with FG parallel.
810 if (itype(i).eq.10) goto 1 ! right?
814 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
815 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
816 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
817 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
818 sigma_d(k,i)=rescore(k,i) ! right expression ?
819 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
821 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
822 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
823 c read (ientin,*) sigma_d(k,i) ! 1st variant
824 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
830 if (waga_dist.ne.0.0d0) lim_odl=ii
831 if (constr_homology.gt.0) call homology_partition
832 if (constr_homology.gt.0) call init_int_table
833 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
834 cd & "lim_xx=",lim_xx
835 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
836 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
840 if (.not.lprn) return
841 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
842 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
843 write (iout,*) "Distance restraints from templates"
845 write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
846 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
848 write (iout,*) "Dihedral angle restraints from templates"
850 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
851 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
853 write (iout,*) "Virtual-bond angle restraints from templates"
855 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
856 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
858 write (iout,*) "SC restraints from templates"
860 write(iout,'(i5,10(4f8.2,4x))') i,
861 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
862 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
870 c -----------------------------------------------------------------
873 c----------------------------------------------------------------------