From dc8aad4b659fef3940a18037dc98775344242ce3 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Thu, 19 May 2016 09:38:09 +0200 Subject: [PATCH] correction of dihedral homology restraints phi(i+3) --- source/cluster/wham/src-M/energy_p_new.F | 2 +- source/cluster/wham/src-M/initialize_p.F | 10 +++--- source/cluster/wham/src-M/readrtns.F | 41 +++++++++++++---------- source/unres/src_MD-M/energy_p_new_barrier.F | 10 ++++-- source/unres/src_MD-M/initialize_p.F | 10 +++--- source/unres/src_MD-M/readrtns_CSA.F | 46 ++++++++++++++------------ source/wham/src-M/energy_p_new.F | 2 +- source/wham/src-M/initialize_p.F | 10 +++--- source/wham/src-M/molread_zs.F | 41 +++++++++++++---------- 9 files changed, 96 insertions(+), 76 deletions(-) diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index 34807fe..4e4a386 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -4177,7 +4177,7 @@ c write (iout,*) idihconstr_start_homo,idihconstr_end_homo do i=idihconstr_start_homo,idihconstr_end_homo kat2=0.0d0 c betai=beta(i,i+1,i+2,i+3) - betai = phi(i+3) + betai = phi(i) c write (iout,*) "betai =",betai do k=1,constr_homology dih_diff(k)=pinorm(dih(k,i)-betai) diff --git a/source/cluster/wham/src-M/initialize_p.F b/source/cluster/wham/src-M/initialize_p.F index 8b01a36..1c83869 100644 --- a/source/cluster/wham/src-M/initialize_p.F +++ b/source/cluster/wham/src-M/initialize_p.F @@ -620,10 +620,10 @@ c include 'COMMON.SETUP' & " lim_dih",lim_dih #ifdef MPL call int_bounds(lim_odl,link_start_homo,link_end_homo) - call int_bounds(lim_dih-nnt+1,idihconstr_start_homo, + call int_bounds(lim_dih,idihconstr_start_homo, & idihconstr_end_homo) - idihconstr_start_homo=idihconstr_start_homo+nnt-1 - idihconstr_end_homo=idihconstr_end_homo+nnt-1 + idihconstr_start_homo=idihconstr_start_homo+nnt-1+3 + idihconstr_end_homo=idihconstr_end_homo+nnt-1+3 if (me.eq.king .or. .not. out1file) & write (iout,*) 'Processor',fg_rank,' CG group',kolor, & ' absolute rank',MyRank, @@ -634,8 +634,8 @@ c include 'COMMON.SETUP' #else link_start_homo=1 link_end_homo=lim_odl - idihconstr_start_homo=nnt - idihconstr_end_homo=lim_dih + idihconstr_start_homo=nnt+3 + idihconstr_end_homo=lim_dih+3 write (iout,*) & ' lim_odl',lim_odl,' link_start=',link_start_homo, & ' link_end',link_end_homo,' lim_dih',lim_dih, diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 1680dab..51ac51f 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -968,6 +968,7 @@ c====------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' + include 'COMMON.NAMES' include 'COMMON.HOMRESTR' c c For new homol impl @@ -992,9 +993,11 @@ c & sigma_odl_temp(maxres,maxres,max_template) c c FP - Nov. 2014 Temporary specifications for new vars c - double precision rescore_tmp,x12,y12,z12,rescore2_tmp + double precision rescore_tmp,x12,y12,z12,rescore2_tmp, + & rescore3_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 + double precision, dimension (max_template,maxres) :: rescore3 character*24 tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -1048,8 +1051,6 @@ cd call flush(iout) c c New c - lim_theta=0 - lim_xx=0 c c Reading HM global scores (prob not required) c @@ -1145,14 +1146,18 @@ 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 + do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, - & idomain_tmp + & rescore3_tmp,idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp + rescore3(k,i_tmp)=rescore3_tmp + write(iout,'(a7,i5,3f10.5,i5)') "rescore", + & i_tmp,rescore2_tmp,rescore_tmp, + & rescore3_tmp,idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -1252,7 +1257,8 @@ 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)) + if (sigma_dih(k,i).ne.0) + & 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) enddo lim_dih=nct-nnt-2 @@ -1284,14 +1290,14 @@ c read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 - sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + if (sigma_theta(k,i).ne.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) enddo endif - lim_theta=nct-nnt-1 if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') @@ -1317,15 +1323,14 @@ 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)) + sigma_d(k,i)=rescore3(k,i) ! right expression ? + if (sigma_d(k,i).ne.0) + & 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? enddo - lim_xx=nct-nnt+1 endif enddo c @@ -1374,17 +1379,19 @@ cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" - do i=nnt+3,lim_dih - write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i), + do i=nnt+3,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(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,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i), + do i=nnt+2,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(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 + do i=nnt,nct write(iout,'(i5,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index b3c4c22..5eca704 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -7060,11 +7060,12 @@ c write (iout,*) idihconstr_start_homo,idihconstr_end_homo do i=idihconstr_start_homo,idihconstr_end_homo kat2=0.0d0 c betai=beta(i,i+1,i+2,i+3) - betai = phi(i+3) + betai = phi(i) c write (iout,*) "betai =",betai do k=1,constr_homology dih_diff(k)=pinorm(dih(k,i)-betai) -c write (iout,*) "dih_diff(",k,") =",dih_diff(k) +cd write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k) +cd & ,sigma_dih(k,i) c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)= c & -(6.28318-dih_diff(i,k)) c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)= @@ -7172,6 +7173,9 @@ c c dtheta_i=theta(j)-thetaref(j,iref) c dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing theta_diff(k)=thetatpl(k,i)-theta(i) +cd write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k) +cd & ,sigma_theta(k,i) + c utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta? @@ -7242,7 +7246,7 @@ c Original sign inverted for calc of gradient dyy=-yytpl(k,i)+yytab(i) ! ibid y dzz=-zztpl(k,i)+zztab(i) ! ibid z c write(iout,*) "dxx, dyy, dzz" -c write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz +cd write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i) c usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index 4d92752..570600a 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -1487,10 +1487,10 @@ cd & " lim_dih",lim_dih #ifdef MPI if (me.eq.king .or. .not. out1file) write (iout,*) "MPI" call int_bounds(lim_odl,link_start_homo,link_end_homo) - call int_bounds(lim_dih-nnt+1,idihconstr_start_homo, + call int_bounds(lim_dih,idihconstr_start_homo, & idihconstr_end_homo) - idihconstr_start_homo=idihconstr_start_homo+nnt-1 - idihconstr_end_homo=idihconstr_end_homo+nnt-1 + idihconstr_start_homo=idihconstr_start_homo+nnt-1+3 + idihconstr_end_homo=idihconstr_end_homo+nnt-1+3 if (me.eq.king .or. .not. out1file) & write (iout,*) 'Processor',fg_rank,' CG group',kolor, & ' absolute rank',MyRank, @@ -1502,8 +1502,8 @@ cd & " lim_dih",lim_dih write (iout,*) "Not MPI" link_start_homo=1 link_end_homo=lim_odl - idihconstr_start_homo=nnt - idihconstr_end_homo=lim_dih + idihconstr_start_homo=nnt+3 + idihconstr_end_homo=lim_dih+3 write (iout,*) & ' lim_odl',lim_odl,' link_start=',link_start_homo, & ' link_end',link_end_homo,' lim_dih',lim_dih, diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index f944233..2a591cf 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -2542,6 +2542,7 @@ c------------------------------------------------------------------------------- include 'COMMON.MD' include 'COMMON.GEO' include 'COMMON.INTERACT' + include 'COMMON.NAMES' c c For new homol impl c @@ -2562,9 +2563,11 @@ c & sigma_odl_temp(maxres,maxres,max_template) c c FP - Nov. 2014 Temporary specifications for new vars c - double precision rescore_tmp,x12,y12,z12,rescore2_tmp + double precision rescore_tmp,x12,y12,z12,rescore2_tmp, + & rescore3_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 + double precision, dimension (max_template,maxres) :: rescore3 character*24 pdbfile,tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -2615,11 +2618,6 @@ cd call flush(iout) lim_odl=0 lim_dih=0 c -c New -c - lim_theta=0 - lim_xx=0 -c write(iout,*) 'nnt=',nnt,'nct=',nct c do i = nnt,nct @@ -2683,14 +2681,18 @@ 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 + do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, - & idomain_tmp + & rescore3_tmp,idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp + rescore3(k,i_tmp)=rescore3_tmp + write(iout,'(a7,i5,3f10.5,i5)') "rescore", + & i_tmp,rescore2_tmp,rescore_tmp, + & rescore3_tmp,idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -2790,7 +2792,8 @@ 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)) + if (sigma_dih(k,i).ne.0) + & 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) enddo lim_dih=nct-nnt-2 @@ -2822,14 +2825,14 @@ c read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 - sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + if (sigma_theta(k,i).ne.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) enddo endif - lim_theta=nct-nnt-1 if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') @@ -2855,15 +2858,14 @@ 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)) + sigma_d(k,i)=rescore3(k,i) ! right expression ? + if (sigma_d(k,i).ne.0) + & 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? enddo - lim_xx=nct-nnt+1 endif enddo c @@ -2894,8 +2896,6 @@ c 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, -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 @@ -2912,17 +2912,19 @@ cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" - do i=nnt+3,lim_dih - write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i), + do i=nnt+3,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(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,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i), + do i=nnt+2,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(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 + do i=nnt,nct write(iout,'(i5,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index d0775ad..a98fe50 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -3673,7 +3673,7 @@ c write (iout,*) idihconstr_start_homo,idihconstr_end_homo do i=idihconstr_start_homo,idihconstr_end_homo kat2=0.0d0 c betai=beta(i,i+1,i+2,i+3) - betai = phi(i+3) + betai = phi(i) c write (iout,*) "betai =",betai do k=1,constr_homology dih_diff(k)=pinorm(dih(k,i)-betai) diff --git a/source/wham/src-M/initialize_p.F b/source/wham/src-M/initialize_p.F index 459c312..3f5f299 100644 --- a/source/wham/src-M/initialize_p.F +++ b/source/wham/src-M/initialize_p.F @@ -632,10 +632,10 @@ c include 'COMMON.SETUP' & " lim_dih",lim_dih #ifdef MPL call int_bounds(lim_odl,link_start_homo,link_end_homo) - call int_bounds(lim_dih-nnt+1,idihconstr_start_homo, + call int_bounds(lim_dih,idihconstr_start_homo, & idihconstr_end_homo) - idihconstr_start_homo=idihconstr_start_homo+nnt-1 - idihconstr_end_homo=idihconstr_end_homo+nnt-1 + idihconstr_start_homo=idihconstr_start_homo+nnt-1+3 + idihconstr_end_homo=idihconstr_end_homo+nnt-1+3 if (me.eq.king .or. .not. out1file) & write (iout,*) 'Processor',fg_rank,' CG group',kolor, & ' absolute rank',MyRank, @@ -646,8 +646,8 @@ c include 'COMMON.SETUP' #else link_start_homo=1 link_end_homo=lim_odl - idihconstr_start_homo=nnt - idihconstr_end_homo=lim_dih + idihconstr_start_homo=nnt+3 + idihconstr_end_homo=lim_dih+3 write (iout,*) & ' lim_odl',lim_odl,' link_start=',link_start_homo, & ' link_end',link_end_homo,' lim_dih',lim_dih, diff --git a/source/wham/src-M/molread_zs.F b/source/wham/src-M/molread_zs.F index dee1afc..0c1a211 100644 --- a/source/wham/src-M/molread_zs.F +++ b/source/wham/src-M/molread_zs.F @@ -339,6 +339,7 @@ c====------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' + include 'include_unres/COMMON.NAMES' include 'COMMON.HOMRESTR' c c For new homol impl @@ -363,9 +364,11 @@ c & sigma_odl_temp(maxres,maxres,max_template) c c FP - Nov. 2014 Temporary specifications for new vars c - double precision rescore_tmp,x12,y12,z12,rescore2_tmp + double precision rescore_tmp,x12,y12,z12,rescore2_tmp, + & rescore3_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 + double precision, dimension (max_template,maxres) :: rescore3 character*24 tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -410,8 +413,6 @@ cd call flush(iout) c c New c - lim_theta=0 - lim_xx=0 c c Reading HM global scores (prob not required) c @@ -494,14 +495,18 @@ 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 + do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, - & idomain_tmp + & rescore3_tmp,idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp + rescore3(k,i_tmp)=rescore3_tmp + write(iout,'(a7,i5,3f10.5,i5)') "rescore", + & i_tmp,rescore2_tmp,rescore_tmp, + & rescore3_tmp,idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -601,7 +606,8 @@ 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)) + if (sigma_dih(k,i).ne.0) + & 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) enddo lim_dih=nct-nnt-2 @@ -633,14 +639,14 @@ c read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 - sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + if (sigma_theta(k,i).ne.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) enddo endif - lim_theta=nct-nnt-1 if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') @@ -666,15 +672,14 @@ 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)) + sigma_d(k,i)=rescore3(k,i) ! right expression ? + if (sigma_d(k,i).ne.0) + & 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? enddo - lim_xx=nct-nnt+1 endif enddo c @@ -723,17 +728,19 @@ cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" - do i=nnt+3,lim_dih - write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i), + do i=nnt+3,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(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,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i), + do i=nnt+2,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(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 + do i=nnt,nct write(iout,'(i5,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) -- 1.7.9.5