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.NAMES'
512 include 'COMMON.HOMRESTR'
517 c include 'include_unres/COMMON.VAR'
520 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
522 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
523 c & sigma_odl_temp(maxres,maxres,max_template)
525 character*24 model_ki_dist, model_ki_angle
526 character*500 controlcard
527 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
528 integer idomain(max_template,maxres)
529 logical lprn /.true./
532 logical unres_pdb,liiflag
534 c FP - Nov. 2014 Temporary specifications for new vars
536 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
537 double precision, dimension (max_template,maxres) :: rescore
538 double precision, dimension (max_template,maxres) :: rescore2
539 character*24 tpl_k_rescore
540 c -----------------------------------------------------------------
541 c Reading multiple PDB ref structures and calculation of retraints
542 c not using pre-computed ones stored in files model_ki_{dist,angle}
544 c -----------------------------------------------------------------
547 c Alternative: reading from input
548 call card_concat(controlcard,.true.)
549 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
550 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
551 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
552 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
553 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
554 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
555 call readi(controlcard,"HOMOL_SET",homol_nset,1)
556 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
557 call readi(controlcard,"IHSET",ihset,1)
558 if (homol_nset.gt.1)then
559 call card_concat(controlcard,.true.)
560 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
561 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
562 write(iout,*) "iset homology_weight "
564 c write(iout,*) i,waga_homology(i)
567 iset=mod(kolor,homol_nset)+1
572 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
574 cd write (iout,*) "nnt",nnt," nct",nct
586 c Reading HM global scores (prob not required)
589 do k=1,constr_homology
593 c open (4,file="HMscore")
594 c do k=1,constr_homology
595 c read (4,*,end=521) hmscore_tmp
596 c hmscore(k)=hmscore_tmp ! Another transformation can be used
597 c write(*,*) "Model", k, ":", hmscore(k)
608 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
610 do k=1,constr_homology
612 read(inp,'(a)') pdbfile
613 c Next stament causes error upon compilation (?)
614 c if(me.eq.king.or. .not. out1file)
615 c write (iout,'(2a)') 'PDB data will be read from file ',
616 c & pdbfile(:ilen(pdbfile))
617 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
618 & pdbfile(:ilen(pdbfile))
619 open(ipdbin,file=pdbfile,status='old',err=33)
621 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
622 & pdbfile(:ilen(pdbfile))
625 c print *,'Begin reading pdb data'
627 c Files containing res sim or local scores (former containing sigmas)
630 write(kic2,'(bz,i2.2)') k
632 tpl_k_rescore="template"//kic2//".sco"
643 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
644 & (crefjlee(j,i+nres),j=1,3)
647 write (iout,*) "read_constr_homology: after reading pdb file"
651 c Distance restraints
654 C Copy the coordinates from reference coordinates (?)
658 c write (iout,*) "c(",j,i,") =",c(j,i)
662 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
664 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
665 open (ientin,file=tpl_k_rescore,status='old')
666 if (nnt.gt.1) rescore(k,1)=0.0d0
667 do irec=nnt,nct ! loop for reading res sim
669 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
672 idomain(k,i_tmp)=idomain_tmp
673 rescore(k,i_tmp)=rescore_tmp
674 rescore2(k,i_tmp)=rescore2_tmp
675 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
676 & i_tmp,rescore2_tmp,rescore_tmp,
680 read (ientin,*,end=1401) rescore_tmp
682 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
683 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
684 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
689 if (waga_dist.ne.0.0d0) then
697 distal=dsqrt(x12*x12+y12*y12+z12*z12)
698 c write (iout,*) k,i,j,distal,dist2_cut
700 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
701 & .and. distal.le.dist2_cut ) then
707 c write (iout,*) "k",k
708 c write (iout,*) "i",i," j",j," constr_homology",
716 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
718 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
719 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
720 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
722 if (odl(k,ii).le.dist_cut) then
723 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
726 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
727 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
729 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
730 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
734 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
744 c Theta, dihedral and SC retraints
746 if (waga_angle.gt.0.0d0) then
747 c open (ientin,file=tpl_k_sigma_dih,status='old')
748 c do irec=1,maxres-3 ! loop for reading sigma_dih
749 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
750 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
751 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
752 c & sigma_dih(k,i+nnt-1)
757 if (idomain(k,i).eq.0) then
761 dih(k,i)=phiref(i) ! right?
762 c read (ientin,*) sigma_dih(k,i) ! original variant
763 c write (iout,*) "dih(",k,i,") =",dih(k,i)
764 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
765 c & "rescore(",k,i-1,") =",rescore(k,i-1),
766 c & "rescore(",k,i-2,") =",rescore(k,i-2),
767 c & "rescore(",k,i-3,") =",rescore(k,i-3)
769 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
770 & rescore(k,i-2)+rescore(k,i-3))/4.0
771 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
772 c write (iout,*) "Raw sigmas for dihedral angle restraints"
773 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
774 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
775 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
776 c Instead of res sim other local measure of b/b str reliability possible
777 if (sigma_dih(k,i).ne.0)
778 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
779 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
784 if (waga_theta.gt.0.0d0) then
785 c open (ientin,file=tpl_k_sigma_theta,status='old')
786 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
787 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
788 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
789 c & sigma_theta(k,i+nnt-1)
794 do i = nnt+2,nct ! right? without parallel.
795 c do i = i=1,nres ! alternative for bounds acc to readpdb?
796 c do i=ithet_start,ithet_end ! with FG parallel.
797 if (idomain(k,i).eq.0) then
801 thetatpl(k,i)=thetaref(i)
802 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
803 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
804 c & "rescore(",k,i-1,") =",rescore(k,i-1),
805 c & "rescore(",k,i-2,") =",rescore(k,i-2)
806 c read (ientin,*) sigma_theta(k,i) ! 1st variant
807 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
808 & rescore(k,i-2))/3.0
809 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
810 if (sigma_theta(k,i).ne.0)
811 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
813 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
814 c rescore(k,i-2) ! right expression ?
815 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
819 if (waga_d.gt.0.0d0) then
820 c open (ientin,file=tpl_k_sigma_d,status='old')
821 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
822 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
823 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
824 c & sigma_d(k,i+nnt-1)
828 do i = nnt,nct ! right? without parallel.
829 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
830 c do i=loc_start,loc_end ! with FG parallel.
831 if (itype(i).eq.10) cycle
832 if (idomain(k,i).eq.0 ) then
839 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
840 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
841 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
842 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
843 sigma_d(k,i)=rescore(k,i) ! right expression ?
844 if (sigma_d(k,i).ne.0)
845 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
847 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
848 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
849 c read (ientin,*) sigma_d(k,i) ! 1st variant
854 c remove distance restraints not used in any model from the list
855 c shift data in all arrays
857 if (waga_dist.ne.0.0d0) then
863 if (ii_in_use(ii).eq.0.and.liiflag) then
867 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
868 & .not.liiflag.and.ii.eq.lim_odl) then
869 if (ii.eq.lim_odl) then
875 do ki=iistart,lim_odl-iishift
876 ires_homo(ki)=ires_homo(ki+iishift)
877 jres_homo(ki)=jres_homo(ki+iishift)
878 ii_in_use(ki)=ii_in_use(ki+iishift)
879 do k=1,constr_homology
880 odl(k,ki)=odl(k,ki+iishift)
881 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
882 l_homo(k,ki)=l_homo(k,ki+iishift)
886 lim_odl=lim_odl-iishift
891 if (constr_homology.gt.0) call homology_partition
892 if (constr_homology.gt.0) call init_int_table
893 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
894 cd & "lim_xx=",lim_xx
895 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
896 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
900 if (.not.lprn) return
901 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
902 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
903 write (iout,*) "Distance restraints from templates"
905 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
906 & ii,ires_homo(ii),jres_homo(ii),
907 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
908 & ki=1,constr_homology)
910 write (iout,*) "Dihedral angle restraints from templates"
912 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
913 & (rad2deg*dih(ki,i),
914 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
916 write (iout,*) "Virtual-bond angle restraints from templates"
918 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
919 & (rad2deg*thetatpl(ki,i),
920 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
922 write (iout,*) "SC restraints from templates"
924 write(iout,'(i5,100(4f8.2,4x))') i,
925 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
926 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
929 c -----------------------------------------------------------------
932 c----------------------------------------------------------------------