From: Cezary Czaplewski Date: Wed, 8 Jun 2016 12:01:29 +0000 (+0200) Subject: single chain correction of dihedral homology restraints phi(i+3) X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=175f9b9db89039e8aaea547c80e7fcb8ab979803;p=unres.git single chain correction of dihedral homology restraints phi(i+3) --- diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F index 3b1c095..7ee0e64 100644 --- a/source/cluster/wham/src/energy_p_new.F +++ b/source/cluster/wham/src/energy_p_new.F @@ -3312,7 +3312,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/initialize_p.F b/source/cluster/wham/src/initialize_p.F index 3fe87fc..e8eac93 100644 --- a/source/cluster/wham/src/initialize_p.F +++ b/source/cluster/wham/src/initialize_p.F @@ -590,10 +590,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, @@ -604,8 +604,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/readrtns.F b/source/cluster/wham/src/readrtns.F index 16fdbeb..34692c5 100644 --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@ -917,6 +917,7 @@ c====------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' + include 'COMMON.NAMES' include 'COMMON.HOMRESTR' c c For new homol impl @@ -995,11 +996,6 @@ 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 do i = nnt,nct @@ -1094,7 +1090,7 @@ 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 @@ -1102,6 +1098,9 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp + write(iout,'(a7,i5,2f10.5,i5)') "rescore", + & i_tmp,rescore2_tmp,rescore_tmp, + & idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -1201,7 +1200,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 @@ -1233,14 +1233,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') @@ -1267,14 +1267,14 @@ 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)) + 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 @@ -1323,17 +1323,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/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index a91710d..115cb57 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -6234,7 +6234,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/unres/src_MD/initialize_p.F b/source/unres/src_MD/initialize_p.F index 369e6bc..53a90bf 100644 --- a/source/unres/src_MD/initialize_p.F +++ b/source/unres/src_MD/initialize_p.F @@ -1415,10 +1415,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, @@ -1430,8 +1430,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/readrtns.F b/source/unres/src_MD/readrtns.F index 0b6f1ae..b51386c 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -2721,6 +2721,7 @@ c------------------------------------------------------------------------------- include 'COMMON.MD' include 'COMMON.GEO' include 'COMMON.INTERACT' + include 'COMMON.NAMES' c c For new homol impl c @@ -2794,11 +2795,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 @@ -2862,7 +2858,7 @@ 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 @@ -2870,6 +2866,9 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp + write(iout,'(a7,i5,2f10.5,i5)') "rescore", + & i_tmp,rescore2_tmp,rescore_tmp, + & idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -2969,7 +2968,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 @@ -3001,14 +3001,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') @@ -3035,14 +3035,14 @@ 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)) + 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 @@ -3091,17 +3091,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/energy_p_new.F b/source/wham/src/energy_p_new.F index 13267da..cf4b2e2 100644 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@ -3386,7 +3386,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/initialize_p.F b/source/wham/src/initialize_p.F index c1c1950..11bf4ab 100644 --- a/source/wham/src/initialize_p.F +++ b/source/wham/src/initialize_p.F @@ -596,10 +596,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, @@ -610,8 +610,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/molread_zs.F b/source/wham/src/molread_zs.F index b11e3fb..579d85a 100644 --- a/source/wham/src/molread_zs.F +++ b/source/wham/src/molread_zs.F @@ -508,6 +508,7 @@ c====------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' + include 'COMMON.NAMES' include 'COMMON.HOMRESTR' c c For new homol impl @@ -663,7 +664,7 @@ 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 @@ -671,6 +672,9 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp + write(iout,'(a7,i5,2f10.5,i5)') "rescore", + & i_tmp,rescore2_tmp,rescore_tmp, + & idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -770,7 +774,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 @@ -802,14 +807,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') @@ -836,14 +841,13 @@ 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)) + 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 @@ -892,17 +896,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)