implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.VAR'
+c include 'include_unres/COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.LOCAL'
include 'COMMON.NAMES'
call reada(controlcard,'SCALSCP',scalscp,1.0d0)
call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
+C Bartek
+ call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
+ call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
+ call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
+ call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
+ write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
+ & " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
r0_corr=cutoff_corr-delt_corr
call readi(controlcard,"NRES",nres,0)
iscode=index(controlcard,"ONE_LETTER")
do i=1,nres
itype(i)=rescode(i,sequence(i),iscode)
enddo
+ if (itype(2).eq.10) then
+ write (iout,*)
+ & "Glycine is the first full residue, initial dummy deleted"
+ do i=1,nres
+ itype(i)=itype(i+1)
+ enddo
+ nres=nres-1
+ endif
+ if (itype(nres).eq.10) then
+ write (iout,*)
+ & "Glycine is the last full residue, terminal dummy deleted"
+ nres=nres-1
+ endif
write (iout,*) "Numeric code:"
write (iout,'(20i4)') (itype(i),i=1,nres)
do i=1,nres-1
if (itype(1).eq.21) nnt=2
if (itype(nres).eq.21) nct=nct-1
write(iout,*) 'NNT=',NNT,' NCT=',NCT
+
+C Juyong:READ init_vars
+C Initialize variables!
+C Juyong:READ read_info
+C READ fragment information!!
+C both routines should be in dfa.F file!!
+
+ if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
+ & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
+ write (iout,*) "Calling init_dfa_vars"
+ call flush(iout)
+ call init_dfa_vars
+ write (iout,*) 'init_dfa_vars finished!'
+ call flush(iout)
+ call read_dfa_info
+ write (iout,*) 'read_dfa_info finished!'
+ call flush(iout)
+ endif
+
+c Read distance restraints
+ if (constr_dist.gt.0) then
+ if (refstr) call read_ref_structure(*11)
+ call read_dist_constr
+ call hpb_partition
+ endif
+
+ if (constr_homology.gt.0) then
+c write (iout,*) "About to call read_constr_homology"
+c call flush(iout)
+ call read_constr_homology
+c write (iout,*) "Exit read_constr_homology"
+c call flush(iout)
+ if (indpdb.gt.0 .or. pdbref) then
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=crefjlee(j,i)
+ cref(j,i)=crefjlee(j,i)
+ enddo
+ enddo
+ endif
+#ifdef DEBUG
+ write (iout,*) "Array C"
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
+ & (c(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Array Cref"
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
+ & (cref(j,i+nres),j=1,3)
+ enddo
+#endif
+#ifdef DEBUG
+ call int_from_cart1(.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
+ enddo
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ enddo
+#endif
+ else
+ homol_nset=0
+ endif
+
+
call setup_var
call init_int_table
if (ns.gt.0) then
endif
write (iout,'(a)')
return
+ 11 stop "Error reading reference structure"
end
c-----------------------------------------------------------------------------
logical function seq_comp(itypea,itypeb,length)
enddo
endif
endif
+ if (ns.gt.0.and.dyn_ss) then
+C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
+C the bond
+ do i=nss+1,nhpb
+C /06/28/2013 Adasko: nss number of full SS bonds
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ call hpb_partition
+ do i=1,ns
+ dyn_ss_mask(iss(i))=.true.
+C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
+c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
+ enddo
+ endif
return
end
c------------------------------------------------------------------------------
return
10 return1
end
+c-------------------------------------------------------------------------------
+ subroutine read_dist_constr
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ integer ifrag_(2,100),ipair_(2,100)
+ double precision wfrag_(100),wpair_(100)
+ character*500 controlcard
+c write (iout,*) "Calling read_dist_constr"
+c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+c call flush(iout)
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,"NFRAG",nfrag_,0)
+ call readi(controlcard,"NPAIR",npair_,0)
+ call readi(controlcard,"NDIST",ndist_,0)
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+ call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+ call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+ call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+ call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+ write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+ write (iout,*) "IFRAG"
+ do i=1,nfrag_
+ write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ enddo
+ write (iout,*) "IPAIR"
+ do i=1,npair_
+ write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+ enddo
+ call flush(iout)
+ if (.not.refstr .and. nfrag_.gt.0) then
+ write (iout,*)
+ & "ERROR: no reference structure to compute distance restraints"
+ write (iout,*)
+ & "Restraints must be specified explicitly (NDIST=number)"
+ stop
+ endif
+ if (nfrag_.lt.2 .and. npair_.gt.0) then
+ write (iout,*) "ERROR: Less than 2 fragments specified",
+ & " but distance restraints between pairs requested"
+ stop
+ endif
+ call flush(iout)
+ do i=1,nfrag_
+ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+ if (ifrag_(2,i).gt.nstart_sup+nsup-1)
+ & ifrag_(2,i)=nstart_sup+nsup-1
+c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ call flush(iout)
+ if (wfrag_(i).gt.0.0d0) then
+ do j=ifrag_(1,i),ifrag_(2,i)-1
+ do k=j+1,ifrag_(2,i)
+ write (iout,*) "j",j," k",k
+ ddjk=dist(j,k)
+ if (constr_dist.eq.1) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ else if (constr_dist.eq.2) then
+ if (ddjk.le.dist_cut) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,npair_
+ if (wpair_(i).gt.0.0d0) then
+ ii = ipair_(1,i)
+ jj = ipair_(2,i)
+ if (ii.gt.jj) then
+ itemp=ii
+ ii=jj
+ jj=itemp
+ endif
+ do j=ifrag_(1,ii),ifrag_(2,ii)
+ do k=ifrag_(1,jj),ifrag_(2,jj)
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ forcon(nhpb)=wpair_(i)
+ dhpb(nhpb)=dist(j,k)
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,ndist_
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+ & ibecarb(i),forcon(nhpb+1)
+ if (forcon(nhpb+1).gt.0.0d0) then
+ nhpb=nhpb+1
+ if (ibecarb(i).gt.0) then
+ ihpb(i)=ihpb(i)+nres
+ jhpb(i)=jhpb(i)+nres
+ endif
+ if (dhpb(nhpb).eq.0.0d0)
+ & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ endif
+ enddo
+ do i=1,nhpb
+ write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
+ & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
+ enddo
+ call flush(iout)
+ return
+ end
+
+
+
+c====-------------------------------------------------------------------
+ subroutine read_constr_homology
+
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.INTERACT'
+ include 'COMMON.HOMRESTR'
+c
+c For new homol impl
+c
+ include 'COMMON.VAR'
+c include 'include_unres/COMMON.VAR'
+c
+
+c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
+c & dist_cut
+c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+c & sigma_odl_temp(maxres,maxres,max_template)
+ character*2 kic2
+ character*24 model_ki_dist, model_ki_angle
+ character*500 controlcard
+ integer ki, i, j, k, l
+ logical lprn /.true./
+ logical unres_pdb
+c
+c FP - Nov. 2014 Temporary specifications for new vars
+c
+ double precision rescore_tmp,x12,y12,z12
+ double precision, dimension (max_template,maxres) :: rescore
+ character*24 tpl_k_rescore
+c -----------------------------------------------------------------
+c Reading multiple PDB ref structures and calculation of retraints
+c not using pre-computed ones stored in files model_ki_{dist,angle}
+c FP (Nov., 2014)
+c -----------------------------------------------------------------
+c
+c
+c Alternative: reading from input
+ call card_concat(controlcard,.true.)
+ call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
+ call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
+ call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
+ call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
+
+ call readi(controlcard,"HOMOL_NSET",homol_nset,1)
+ call readi(controlcard,"IHSET",ihset,1)
+ if (homol_nset.gt.1)then
+ call card_concat(controlcard,.true.)
+ read(controlcard,*) (waga_homology(i),i=1,homol_nset)
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "iset homology_weight "
+c do i=1,homol_nset
+c write(iout,*) i,waga_homology(i)
+c enddo
+ endif
+ iset=mod(kolor,homol_nset)+1
+ else
+ iset=1
+ waga_homology(1)=1.0
+ endif
+c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
+
+cd write (iout,*) "nnt",nnt," nct",nct
+cd call flush(iout)
+
+
+ lim_odl=0
+ lim_dih=0
+c
+c New
+c
+ lim_theta=0
+ lim_xx=0
+c
+c Reading HM global scores (prob not required)
+c
+c open (4,file="HMscore")
+c do k=1,constr_homology
+c read (4,*,end=521) hmscore_tmp
+c hmscore(k)=hmscore_tmp ! Another transformation can be used
+c write(*,*) "Model", k, ":", hmscore(k)
+c enddo
+c521 continue
+
+c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+
+ do k=1,constr_homology
+
+ read(inp,'(a)') pdbfile
+c Next stament causes error upon compilation (?)
+c if(me.eq.king.or. .not. out1file)
+c write (iout,'(2a)') 'PDB data will be read from file ',
+c & pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a)') 'Error opening PDB file.'
+ stop
+ 34 continue
+c print *,'Begin reading pdb data'
+c
+c Files containing res sim or local scores (former containing sigmas)
+c
+
+ write(kic2,'(bz,i2.2)') k
+
+ tpl_k_rescore="template"//kic2//".sco"
+c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
+c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
+c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
+c tpl_k_sigma_d="template"//kic2//".sigma_d"
+
+ unres_pdb=.false.
+ call readpdb
+ do i=1,2*nres
+ do j=1,3
+ crefjlee(j,i)=c(j,i)
+ enddo
+ enddo
+#ifdef DEBUG
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
+ & (crefjlee(j,i+nres),j=1,3)
+ enddo
+#endif
+ write (iout,*) "read_constr_homology: after reading pdb file"
+ call flush(iout)
+
+c
+c Distance restraints
+c
+c ... --> odl(k,ii)
+C Copy the coordinates from reference coordinates (?)
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=cref(j,i)
+c write (iout,*) "c(",j,i,") =",c(j,i)
+ enddo
+ enddo
+c
+c From read_dist_constr (commented out 25/11/2014 <-> res sim)
+c
+c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
+ open (ientin,file=tpl_k_rescore,status='old')
+ do irec=1,maxdim ! loop for reading res sim
+ if (irec.eq.1) then
+ rescore(k,irec)=0.0d0
+ goto 1301
+ endif
+ read (ientin,*,end=1401) rescore_tmp
+c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
+ rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
+c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
+ 1301 continue
+ enddo
+ 1401 continue
+ close (ientin)
+c open (ientin,file=tpl_k_sigma_odl,status='old')
+c do irec=1,maxdim ! loop for reading sigma_odl
+c read (ientin,*,end=1401) i, j,
+c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
+c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
+c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
+c enddo
+c 1401 continue
+c close (ientin)
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i = nnt,nct-2 ! right? without parallel.
+ do j=i+2,nct ! right?
+c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
+c do j=i+2,nres ! ibid
+c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
+c do j=i+2,nct ! ibid
+ ii=ii+1
+c write (iout,*) "k",k
+c write (iout,*) "i",i," j",j," constr_homology",
+c & constr_homology
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+c
+c Attempt to replace dist(i,j) by its definition in ...
+c
+ x12=c(1,i)-c(1,j)
+ y12=c(2,i)-c(2,j)
+ z12=c(3,i)-c(3,j)
+ distal=dsqrt(x12*x12+y12*y12+z12*z12)
+ odl(k,ii)=distal
+c
+c odl(k,ii)=dist(i,j)
+c write (iout,*) "dist(",i,j,") =",dist(i,j)
+c write (iout,*) "distal = ",distal
+c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,j,") =",rescore(k,j)
+c
+c Calculation of sigma from res sim
+c
+c if (odl(k,ii).le.6.0d0) then
+c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
+c Other functional forms possible depending on odl(k,ii), eg.
+c
+ if (odl(k,ii).le.dist_cut) then
+ sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
+c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
+c write (iout,*) "c sigma_odl",k,ii,sigma_odl(k,ii)
+ else
+#ifdef OLDSIGMA
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
+#else
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+c write (iout,*) "d sigma_odl",k,ii,sigma_odl(k,ii),
+c & odl(k,ii),dist_cut
+#endif
+c Following expr replaced by a positive exp argument
+c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+
+c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
+c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+ endif
+c
+ sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
+c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
+c
+c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
+c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
+ enddo
+c read (ientin,*) sigma_odl(k,ii) ! 1st variant
+ enddo
+c lim_odl=ii
+c if (constr_homology.gt.0) call homology_partition
+ endif
+c
+c Theta, dihedral and SC retraints
+c
+ if (waga_angle.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_dih,status='old')
+c do irec=1,maxres-3 ! loop for reading sigma_dih
+c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
+c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
+c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_dih(k,i+nnt-1)
+c enddo
+c1402 continue
+c close (ientin)
+ do i = nnt+3,nct ! right? without parallel.
+c do i=1,nres ! alternative for bounds acc to readpdb?
+c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
+c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
+ dih(k,i)=phiref(i) ! right?
+c read (ientin,*) sigma_dih(k,i) ! original variant
+c write (iout,*) "dih(",k,i,") =",dih(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,i-1,") =",rescore(k,i-1),
+c & "rescore(",k,i-2,") =",rescore(k,i-2),
+c & "rescore(",k,i-3,") =",rescore(k,i-3)
+
+ sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)+rescore(k,i-3) ! right expression ?
+c
+c write (iout,*) "Raw sigmas for dihedral angle restraints"
+c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
+c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
+c Instead of res sim other local measure of b/b str reliability possible
+ sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
+ if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
+c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
+ enddo
+ endif
+
+ if (waga_theta.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_theta,status='old')
+c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
+c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
+c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_theta(k,i+nnt-1)
+c enddo
+c1403 continue
+c close (ientin)
+
+ do i = nnt+2,nct ! right? without parallel.
+c do i = i=1,nres ! alternative for bounds acc to readpdb?
+c do i=ithet_start,ithet_end ! with FG parallel.
+ thetatpl(k,i)=thetaref(i)
+c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,i-1,") =",rescore(k,i-1),
+c & "rescore(",k,i-2,") =",rescore(k,i-2)
+c read (ientin,*) sigma_theta(k,i) ! 1st variant
+ sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2) ! right expression ?
+ sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+
+c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+c rescore(k,i-2) ! right expression ?
+c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
+ if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
+ enddo
+ endif
+
+ if (waga_d.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_d,status='old')
+c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
+c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
+c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_d(k,i+nnt-1)
+c enddo
+c1404 continue
+ close (ientin)
+
+ do i = nnt,nct ! right? without parallel.
+c do i=2,nres-1 ! alternative for bounds acc to readpdb?
+c do i=loc_start,loc_end ! with FG parallel.
+ if (itype(i).eq.10) goto 1 ! right?
+ xxtpl(k,i)=xxref(i)
+ yytpl(k,i)=yyref(i)
+ zztpl(k,i)=zzref(i)
+c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
+c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
+c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
+ sigma_d(k,i)=rescore(k,i) ! right expression ?
+ sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+
+c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
+c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
+c read (ientin,*) sigma_d(k,i) ! 1st variant
+ if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
+ 1 continue
+ enddo
+ endif
+ close(ientin)
+ enddo
+ if (waga_dist.ne.0.0d0) lim_odl=ii
+ if (constr_homology.gt.0) call homology_partition
+ if (constr_homology.gt.0) call init_int_table
+cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
+cd & "lim_xx=",lim_xx
+c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c
+c Print restraints
+c
+ if (.not.lprn) return
+cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write (iout,*) "Distance restraints from templates"
+ do ii=1,lim_odl
+ write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
+ & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "Dihedral angle restraints from templates"
+ do i=nnt+3,lim_dih
+ write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+ & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "Virtual-bond angle restraints from templates"
+ do i=nnt+2,lim_theta
+ write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+ & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "SC restraints from templates"
+ do i=nnt,lim_xx
+ write(iout,'(i5,10(4f8.2,4x))') i,
+ & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
+ & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
+ enddo
+ endif
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+c -----------------------------------------------------------------
+ return
+ end
+c----------------------------------------------------------------------
projects / unres.git / blobdiff