C READ fragment information!!
C both routines should be in dfa.F file!!
+#ifdef DFA
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
call init_dfa_vars
call read_dfa_info
print*, 'read_dfa_info finished!'
endif
+#endif
C
if (pdbref) then
if(me.eq.king.or..not.out1file)
include 'COMMON.IOUNITS'
include 'COMMON.MD'
include 'COMMON.CONTROL'
+ integer iset1
read(inp,*) nset,nfrag,npair,nfrag_back
if(me.eq.king.or..not.out1file)
& write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
& " nfrag_back",nfrag_back
- do iset=1,nset
- read(inp,*) mset(iset)
+ do iset1=1,nset
+ read(inp,*) mset(iset1)
do i=1,nfrag
- read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
- & qinfrag(i,iset)
+ read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
+ & qinfrag(i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
- & ifrag(2,i,iset), qinfrag(i,iset)
+ & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
+ & ifrag(2,i,iset1), qinfrag(i,iset1)
enddo
do i=1,npair
- read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
- & qinpair(i,iset)
+ read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
+ & qinpair(i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
- & ipair(2,i,iset), qinpair(i,iset)
+ & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
+ & ipair(2,i,iset1), qinpair(i,iset1)
enddo
do i=1,nfrag_back
- read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
- & wfrag_back(3,i,iset),
- & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
+ & wfrag_back(3,i,iset1),
+ & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
if(me.eq.king.or..not.out1file)
- & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
- & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ & write(iout,*) "A",i,wfrag_back(1,i,iset1),
+ & wfrag_back(2,i,iset1),
+ & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
+ & ifrag_back(2,i,iset1)
enddo
enddo
return
character*2 kic2
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
- integer ki, i, j, k, l
+ integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
logical lprn /.true./
c
c FP - Nov. 2014 Temporary specifications for new vars
c
- double precision rescore_tmp,x12,y12,z12
+ double precision rescore_tmp,x12,y12,z12,rescore2_tmp
double precision, dimension (max_template,maxres) :: rescore
+ double precision, dimension (max_template,maxres) :: rescore2
character*24 pdbfile,tpl_k_rescore
c -----------------------------------------------------------------
c Reading multiple PDB ref structures and calculation of retraints
call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
call readi(controlcard,"HOMOL_NSET",homol_nset,1)
+ read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
if (homol_nset.gt.1)then
call card_concat(controlcard)
read(controlcard,*) (waga_homology(i),i=1,homol_nset)
lim_theta=0
lim_xx=0
c
-c Reading HM global scores (prob not required)
+ write(iout,*) 'nnt=',nnt,'nct=',nct
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 i = nnt,nct
+ do k=1,constr_homology
+ idomain(k,i)=0
+ enddo
+ enddo
+
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+ ii_in_use(ii)=0
+ enddo
+ enddo
+
do k=1,constr_homology
+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
+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')
+ if (nnt.gt.1) rescore(k,1)=0.0d0
+ do irec=nnt,maxdim ! loop for reading res sim
+ if (read2sigma) then
+ read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
+ & idomain_tmp
+ i_tmp=i_tmp+nnt-1
+ idomain(k,i_tmp)=idomain_tmp
+ rescore(k,i_tmp)=0.5d0*(rescore_tmp+0.5d0)
+ rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0)
+ else
+ idomain(k,irec)=1
+ 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)
+ endif
+ enddo
+ 1401 continue
+ close (ientin)
+
read(inp,'(a)') pdbfile
c Next stament causes error upon compilation (?)
c if(me.eq.king.or. .not. out1file)
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
-c
-c Distance restraints
+ if (read2sigma) then
+ call readpdb_template(k)
+ else
+ call readpdb
+ endif
c
-c ... --> odl(k,ii)
C Copy the coordinates from reference coordinates (?)
do i=1,2*nres
do j=1,3
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.
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
+
+ if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then
+
ii=ii+1
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+
c write (iout,*) "k",k
c write (iout,*) "i",i," j",j," constr_homology",
c & constr_homology
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
+ if (read2sigma) then
+ sigma_odl(k,ii)=0
+ do ik=i,j
+ sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
+ enddo
+ sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
+ if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
+ & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+ else
+ 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)
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)
+#endif
+ endif
+ 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
+ else
+ ii=ii+1
+ l_homo(k,ii)=.false.
+ endif
enddo
-c read (ientin,*) sigma_odl(k,ii) ! 1st variant
enddo
-c lim_odl=ii
-c if (constr_homology.gt.0) call homology_partition
+ lim_odl=ii
endif
c
c Theta, dihedral and SC retraints
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.
+ if (idomain(k,i).eq.0) then
+ sigma_dih(k,i)=0.0
+ cycle
+ endif
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 & "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 ?
+ sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)+rescore(k,i-3))/4.0
c
c write (iout,*) "Raw sigmas for dihedral angle restraints"
c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
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
+ lim_dih=nct-nnt-2
endif
if (waga_theta.gt.0.0d0) then
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.
+ if (idomain(k,i).eq.0) then
+ sigma_theta(k,i)=0.0
+ cycle
+ endif
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)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)/3.0
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
+ lim_theta=nct-nnt-1
if (waga_d.gt.0.0d0) then
c open (ientin,file=tpl_k_sigma_d,status='old')
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?
+ if (itype(i).eq.10) cycle
+ if (idomain(k,i).eq.0 ) then
+ sigma_d(k,i)=0.0
+ cycle
+ endif
xxtpl(k,i)=xxref(i)
yytpl(k,i)=yyref(i)
zztpl(k,i)=zzref(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
+ lim_xx=nct-nnt+1
endif
- close(ientin)
enddo
- if (waga_dist.ne.0.0d0) lim_odl=ii
+c
+c remove distance restraints not used in any model from the list
+c shift data in all arrays
+c
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i=nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+ if (ii_in_use(ii).eq.0) then
+ do ki=ii,lim_odl-1
+ ires_homo(ki)=ires_homo(ki+1)
+ jres_homo(ki)=jres_homo(ki+1)
+ ii_in_use(ki)=ii_in_use(ki+1)
+ do k=1,constr_homology
+ odl(k,ki)=odl(k,ki+1)
+ sigma_odl(k,ki)=sigma_odl(k,ki+1)
+ l_homo(k,ki)=l_homo(k,ki+1)
+ enddo
+ enddo
+ ii=ii-1
+ lim_odl=lim_odl-1
+ endif
+ enddo
+ enddo
+ endif
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,
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)
+ write(iout,'(3i5,10(2f8.2,1x,l1,4x))')
+ & ii,ires_homo(ii),jres_homo(ii),
+ & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
+ & ki=1,constr_homology)
enddo
write (iout,*) "Dihedral angle restraints from templates"
do i=nnt+3,lim_dih