From: Cezary Czaplewski Date: Thu, 5 May 2016 17:09:44 +0000 (+0200) Subject: read2sigma in wham and cluster_wham X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=42251e287a50083c0337e4d78b357b370f832666;p=unres.git read2sigma in wham and cluster_wham some other changes for dfa and energy to match unres --- diff --git a/source/cluster/wham/src/COMMON.CONTROL b/source/cluster/wham/src/COMMON.CONTROL index fe7c63a..efed33f 100644 --- a/source/cluster/wham/src/COMMON.CONTROL +++ b/source/cluster/wham/src/COMMON.CONTROL @@ -3,16 +3,19 @@ & constr_homology,homol_nset,iset real*8 waga_homology real*8 waga_dist, waga_angle,waga_theta, waga_d, dist_cut + & ,dist2_cut logical refstr,pdbref,punch_dist,print_dist,caonly,lside, & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx, & with_dihed_constr,out1file,print_homology_restraints, & print_contact_map,print_homology_models + & ,read2sigma,l_homo common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2, & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int, & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,constr_dist, & with_dihed_constr, constr_homology,homol_nset,out1file, - & print_contact_map + & print_contact_map,read2sigma,l_homo(max_template,maxdim) common /cntrlr/ waga_homology(10), - & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,iset, + & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,dist2_cut, + & iset, & print_homology_restraints,print_homology_models diff --git a/source/cluster/wham/src/COMMON.DFA b/source/cluster/wham/src/COMMON.DFA index 3c74bc8..c6add4f 100644 --- a/source/cluster/wham/src/COMMON.DFA +++ b/source/cluster/wham/src/COMMON.DFA @@ -94,3 +94,8 @@ c & ,DFAEXP & WSHET(MAXRES,MAXRES), EDFABET, & CK(4),SCK(4),S1(4),S2(4) c & ,DFAEXP(15001), + + DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/ + DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/ + DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/ + DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/ diff --git a/source/cluster/wham/src/DIMENSIONS b/source/cluster/wham/src/DIMENSIONS index d9f0c81..d8cdbf4 100644 --- a/source/cluster/wham/src/DIMENSIONS +++ b/source/cluster/wham/src/DIMENSIONS @@ -9,7 +9,7 @@ C Max. number of processors. parameter (maxprocs=16) C Max. number of AA residues integer maxres,maxres2 - parameter (maxres=600) + parameter (maxres=800) C Appr. max. number of interaction sites parameter (maxres2=2*maxres) C Max. number of variables diff --git a/source/cluster/wham/src/dfa.F b/source/cluster/wham/src/dfa.F index 922f933..00cf12c 100644 --- a/source/cluster/wham/src/dfa.F +++ b/source/cluster/wham/src/dfa.F @@ -1,14 +1,3 @@ - block data dfadata - include 'DIMENSIONS' - include 'COMMON.INTERACT' - include 'COMMON.DFA' - DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/ - DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/ - DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/ - DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/ - end -c------------------------------------------------------------------------ - subroutine init_dfa_vars include 'DIMENSIONS' @@ -841,9 +830,9 @@ c write(*,*) n1atom, j, dist c write(*,*) 'case1:',tmp_n - t1dx=t1dx+0.0d0 - t1dy=t1dy+0.0d0 - t1dz=t1dz+0.0d0 +cc t1dx=t1dx+0.0d0 +cc t1dy=t1dy+0.0d0 +cc t1dz=t1dz+0.0d0 t2dx(j)=0.0d0 t2dy(j)=0.0d0 t2dz(j)=0.0d0 @@ -887,9 +876,9 @@ c neighbor atoms tmp_n=tmp_n+1.0d0 c write(*,*) 'case4:',tmp_n - t1dx=t1dx+0.0d0 - t1dy=t1dy+0.0d0 - t1dz=t1dz+0.0d0 +cc t1dx=t1dx+0.0d0 +cc t1dy=t1dy+0.0d0 +cc t1dz=t1dz+0.0d0 t2dx(j)=0.0d0 t2dy(j)=0.0d0 t2dz(j)=0.0d0 @@ -1046,8 +1035,8 @@ C End of sheet variables double precision enebet enebet=0.0d0 - bx=0.0d0;by=0.0d0;bz=0.0d0 - shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0 +c bx=0.0d0;by=0.0d0;bz=0.0d0 +c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0 gdfab=0.0d0 @@ -2130,6 +2119,9 @@ c------------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc****************************************************************************** c real*8 dfaexp(15001) @@ -2163,6 +2155,9 @@ cc********************************************************************** real*8 uum, uup real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2 + real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca) + common /shetf/ shetfx,shetfy,shetfz + common /sheca/ bx,by,bz common /phys1/ inb,nmax,iselect common /kyori2/ dis @@ -2181,6 +2176,7 @@ cc********************************************************************** common /sinsin/ sth real*8 r_pair_mat(maxca,maxca) + real*8 e_gcont,fprim_gcont,de_gcont ci integer istrand(maxca,maxca) ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca) ci common /shetest/ istrand,istrand_p,istrand_m @@ -2201,6 +2197,12 @@ c do i=1,inb-7 do j=i+4,inb-3 + + if (dis(i,j).lt.dfa_cutoff) then + call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + + ip=i+1 ipp=i+2 jp=j+1 @@ -2415,7 +2417,30 @@ c uum=vbetam1(i,j)+vbetam2(i,j) vbet(i,j)=uup+uum vbetp=vbetp+uup vbetm=vbetm+uum - vbeta=vbeta+vbet(i,j) + vbeta=vbeta+vbet(i,j)*e_gcont + + + if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then +c gradient correction from gcont + de_gcont=vbet(i,j)*fprim_gcont/dis(i,j) + shetfx(i)=shetfx(i) + de_gcont*rx(i,j) + shetfy(i)=shetfy(i) + de_gcont*ry(i,j) + shetfz(i)=shetfz(i) + de_gcont*rz(i,j) + + shetfx(j)=shetfx(j) - de_gcont*rx(i,j) + shetfy(j)=shetfy(j) - de_gcont*ry(i,j) + shetfz(j)=shetfz(j) - de_gcont*rz(i,j) + +c energy correction from gcont + vbet(i,j)=vbet(i,j)*e_gcont + vbetap(i,j)=vbetap(i,j)*e_gcont + vbetap1(i,j)=vbetap1(i,j)*e_gcont + vbetap2(i,j)=vbetap2(i,j)*e_gcont + vbetam(i,j)=vbetam(i,j)*e_gcont + vbetam1(i,j)=vbetam1(i,j)*e_gcont + vbetam2(i,j)=vbetam2(i,j)*e_gcont + endif + ci elseif(istrand(i,j).eq.0)then ci vbet(i,j)=0 @@ -2434,7 +2459,16 @@ c write(*,*) 'vbetp:',vbetp c write(*,*) 'vbetm:',vbetm c write(*,*) 'vbet(i,j):',vbet(i,j) c stop - + + else + vbetap(i,j)=0 + vbetap1(i,j)=0 + vbetap2(i,j)=0 + vbetam(i,j)=0 + vbetam1(i,j)=0 + vbetam2(i,j)=0 + vbet(i,j)=0 + endif enddo enddo @@ -2458,6 +2492,9 @@ c------------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbet(maxca,maxca) real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) @@ -2516,6 +2553,7 @@ c local variables real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2 + real*8 e_gcont,fprim_gcont C-------------------------------------------------------------------------------- do i=4,inb-4 im3=i-3 @@ -2755,6 +2793,9 @@ c---------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -2795,11 +2836,17 @@ c local variables real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b + real*8 e_gcont,fprim_gcont c******************************************************************************** do i=3,inb-5 imm=i-2 im=i-1 do j=i+2,inb-3 + + if (dis(imm,j).lt.dfa_cutoff) then + call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + jp=j+1 jpp=j+2 @@ -2923,7 +2970,8 @@ ci & .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j) ! shefz(i,5)=shefz(i,5) ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j) - + + endif ci endif enddo @@ -2936,6 +2984,9 @@ c--------------------------------------------------------------------------c implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -2976,11 +3027,17 @@ C local variables real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b + real*8 e_gcont,fprim_gcont C******************************************************************************** do i=2,inb-6 ip=i+1 im=i-1 do j=i+3,inb-3 + + if (dis(im,j).lt.dfa_cutoff) then + call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + jp=j+1 jpp=j+2 @@ -3102,7 +3159,8 @@ ci & .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j) ! shefz(i,6)=shefz(i,6) ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j) - + + endif ci endif enddo @@ -3115,6 +3173,9 @@ c----------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -3155,12 +3216,18 @@ C local variables real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8 + real*8 e_gcont,fprim_gcont C******************************************************************************** do j=7,inb-1 jm=j-1 jmm=j-2 do i=1,j-6 + + if (dis(i,jmm).lt.dfa_cutoff) then + call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + ip=i+1 ipp=i+2 @@ -3284,6 +3351,7 @@ ci & .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then ! shefz(j,11)=shefz(j,11) ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm) + endif ci endif enddo @@ -3296,6 +3364,9 @@ c----------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -3335,11 +3406,17 @@ C local variables real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8 + real*8 e_gcont,fprim_gcont !c*************************************************************************c do j=6,inb-2 jp=j+1 jm=j-1 do i=1,j-5 + + if (dis(i,jm).lt.dfa_cutoff) then + call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + ip=i+1 ipp=i+2 @@ -3456,6 +3533,8 @@ ci & .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm) $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm) + endif + ci endif ENDDO diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F index 8ce7e5b..3b1c095 100644 --- a/source/cluster/wham/src/energy_p_new.F +++ b/source/cluster/wham/src/energy_p_new.F @@ -3127,6 +3127,7 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d dij=dist(i,j) c write (iout,*) "dij(",i,j,") =",dij do k=1,constr_homology + if(.not.l_homo(k,ii)) cycle distance(k)=odl(k,ii)-dij c write (iout,*) "distance(",k,") =",distance(k) c @@ -3146,7 +3147,17 @@ c endif enddo - min_odl=minval(distancek) +c min_odl=minval(distancek) + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + & min_odl=distancek(kk) + enddo c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij @@ -3159,6 +3170,7 @@ c write (iout,* )"min_odl",min_odl c Nie wiem po co to liczycie jeszcze raz! c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ c & (2*(sigma_odl(i,j,k))**2)) + if(.not.l_homo(k,ii)) cycle if (waga_dist.ge.0.0d0) then c c For Gaussian-type Urestr @@ -3209,6 +3221,7 @@ c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) c & *waga_dist)+min_odl c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist c + if(.not.l_homo(k,ii)) cycle if (waga_dist.ge.0.0d0) then c For Gaussian-type Urestr c @@ -4018,7 +4031,7 @@ C coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-3).ne.ntyp1) then + if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -5133,6 +5146,7 @@ c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor esccor=0.0D0 do i=itau_start,itau_end esccor_ii=0.0D0 + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) phii=phi(i) diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F index ac18e17..16fdbeb 100644 --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@ -932,14 +932,18 @@ 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 + integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp + integer idomain(max_template,maxres) + integer ilen + external ilen logical lprn logical unres_pdb 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 tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -963,8 +967,9 @@ c Alternative: reading from input 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 reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) call readi(controlcard,"HOMOL_NSET",homol_nset,1) + read2sigma=(index(controlcard,'READ2SIGMA').gt.0) if (homol_nset.gt.1)then call readi(controlcard,"ISET",iset,1) call card_concat(controlcard) @@ -997,6 +1002,11 @@ c c c Reading HM global scores (prob not required) c + do i = nnt,nct + do k=1,constr_homology + idomain(k,i)=0 + enddo + enddo c open (4,file="HMscore") c do k=1,constr_homology c read (4,*,end=521) hmscore_tmp @@ -1005,6 +1015,13 @@ c write(*,*) "Model", k, ":", hmscore(k) c enddo c521 continue + ii=0 + do i = nnt,nct-2 + do j=i+2,nct + ii=ii+1 + ii_in_use(ii)=0 + enddo + enddo c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d write (iout,*) "CONSTR_HOMOLOGY",constr_homology @@ -1016,9 +1033,12 @@ 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)) + write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file', + & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 - 33 write (iout,'(a)') 'Error opening PDB file.' + 33 write (iout,'(a,5x,a)') 'Error opening PDB file', + & pdbfile(:ilen(pdbfile)) stop 34 continue c print *,'Begin reading pdb data' @@ -1029,10 +1049,6 @@ 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 @@ -1077,93 +1093,79 @@ 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 + 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)=rescore_tmp + rescore2(k,i_tmp)=rescore2_tmp + 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 + 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 + endif 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) + 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 + do i = nnt,nct-2 + do j=i+2,nct + + 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) +c write (iout,*) k,i,j,distal,dist2_cut + + if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 + & .and. distal.le.dist2_cut ) 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 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) - else + 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) + else #ifdef OLDSIGMA - sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else - sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #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 + endif + sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) + else + ii=ii+1 + l_homo(k,ii)=.false. 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 + lim_odl=ii endif c c Theta, dihedral and SC retraints @@ -1178,10 +1180,11 @@ 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. + do i = nnt+3,nct + 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) @@ -1191,8 +1194,8 @@ 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))/4.0 ! right expression ? -c + & rescore(k,i-2)+rescore(k,i-3))/4.0 +c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 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)* @@ -1200,9 +1203,8 @@ 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 + lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then @@ -1218,6 +1220,10 @@ 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. + 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), @@ -1225,15 +1231,16 @@ 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))/3.0 ! right expression ? + & 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)) 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') @@ -1243,12 +1250,15 @@ c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use 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) @@ -1263,12 +1273,36 @@ 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, @@ -1283,22 +1317,24 @@ 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) + write(iout,'(3i5,100(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 - write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i), + write (iout,'(i5,100(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), + write (iout,'(i5,100(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, + 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) enddo @@ -1306,5 +1342,4 @@ cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d c ----------------------------------------------------------------- return end - - +c---------------------------------------------------------------------- diff --git a/source/wham/src/COMMON.CONTROL b/source/wham/src/COMMON.CONTROL index 209c22a..4b021d0 100644 --- a/source/wham/src/COMMON.CONTROL +++ b/source/wham/src/COMMON.CONTROL @@ -2,14 +2,17 @@ & ensembles,constr_dist,constr_homology,homol_nset, & iset,ihset real*8 waga_homology - real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut + real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut, + & dist2_cut logical refstr,pdbref,punch_dist,print_rms,caonly,verbose, & merge_helices,bxfile,cxfile,histfile,entfile,zscfile, - & rmsrgymap,with_dihed_constr,check_conf,histout,out1file + & rmsrgymap,with_dihed_constr,check_conf,histout,out1file, + & read2sigma,l_homo common /cntrl/ iscode,indpdb,refstr,pdbref,outpdb,outmol2, & punch_dist,print_rms,caonly,verbose,icomparfunc,pdbint, & merge_helices,bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap, & ensembles,with_dihed_constr,check_conf,histout,constr_dist, - & constr_homology,out1file,homol_nset + & constr_homology,out1file,homol_nset,read2sigma common /homol/ waga_homology(maxR), - & waga_dist,waga_angle,waga_theta,waga_d,dist_cut,iset,ihset + & waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut, + & iset,ihset,l_homo(max_template,maxdim) diff --git a/source/wham/src/COMMON.DFA b/source/wham/src/COMMON.DFA index 3c74bc8..c6add4f 100644 --- a/source/wham/src/COMMON.DFA +++ b/source/wham/src/COMMON.DFA @@ -94,3 +94,8 @@ c & ,DFAEXP & WSHET(MAXRES,MAXRES), EDFABET, & CK(4),SCK(4),S1(4),S2(4) c & ,DFAEXP(15001), + + DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/ + DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/ + DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/ + DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/ diff --git a/source/wham/src/DIMENSIONS b/source/wham/src/DIMENSIONS index 50ccc95..d6f4b71 100644 --- a/source/wham/src/DIMENSIONS +++ b/source/wham/src/DIMENSIONS @@ -13,7 +13,7 @@ C Max. number of coarse-grain processors c parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres - parameter (maxres=600) + parameter (maxres=800) c parameter (maxres=400) C Appr. max. number of interaction sites integer maxres2 diff --git a/source/wham/src/dfa.F b/source/wham/src/dfa.F index 922f933..00cf12c 100644 --- a/source/wham/src/dfa.F +++ b/source/wham/src/dfa.F @@ -1,14 +1,3 @@ - block data dfadata - include 'DIMENSIONS' - include 'COMMON.INTERACT' - include 'COMMON.DFA' - DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/ - DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/ - DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/ - DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/ - end -c------------------------------------------------------------------------ - subroutine init_dfa_vars include 'DIMENSIONS' @@ -841,9 +830,9 @@ c write(*,*) n1atom, j, dist c write(*,*) 'case1:',tmp_n - t1dx=t1dx+0.0d0 - t1dy=t1dy+0.0d0 - t1dz=t1dz+0.0d0 +cc t1dx=t1dx+0.0d0 +cc t1dy=t1dy+0.0d0 +cc t1dz=t1dz+0.0d0 t2dx(j)=0.0d0 t2dy(j)=0.0d0 t2dz(j)=0.0d0 @@ -887,9 +876,9 @@ c neighbor atoms tmp_n=tmp_n+1.0d0 c write(*,*) 'case4:',tmp_n - t1dx=t1dx+0.0d0 - t1dy=t1dy+0.0d0 - t1dz=t1dz+0.0d0 +cc t1dx=t1dx+0.0d0 +cc t1dy=t1dy+0.0d0 +cc t1dz=t1dz+0.0d0 t2dx(j)=0.0d0 t2dy(j)=0.0d0 t2dz(j)=0.0d0 @@ -1046,8 +1035,8 @@ C End of sheet variables double precision enebet enebet=0.0d0 - bx=0.0d0;by=0.0d0;bz=0.0d0 - shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0 +c bx=0.0d0;by=0.0d0;bz=0.0d0 +c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0 gdfab=0.0d0 @@ -2130,6 +2119,9 @@ c------------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc****************************************************************************** c real*8 dfaexp(15001) @@ -2163,6 +2155,9 @@ cc********************************************************************** real*8 uum, uup real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2 + real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca) + common /shetf/ shetfx,shetfy,shetfz + common /sheca/ bx,by,bz common /phys1/ inb,nmax,iselect common /kyori2/ dis @@ -2181,6 +2176,7 @@ cc********************************************************************** common /sinsin/ sth real*8 r_pair_mat(maxca,maxca) + real*8 e_gcont,fprim_gcont,de_gcont ci integer istrand(maxca,maxca) ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca) ci common /shetest/ istrand,istrand_p,istrand_m @@ -2201,6 +2197,12 @@ c do i=1,inb-7 do j=i+4,inb-3 + + if (dis(i,j).lt.dfa_cutoff) then + call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + + ip=i+1 ipp=i+2 jp=j+1 @@ -2415,7 +2417,30 @@ c uum=vbetam1(i,j)+vbetam2(i,j) vbet(i,j)=uup+uum vbetp=vbetp+uup vbetm=vbetm+uum - vbeta=vbeta+vbet(i,j) + vbeta=vbeta+vbet(i,j)*e_gcont + + + if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then +c gradient correction from gcont + de_gcont=vbet(i,j)*fprim_gcont/dis(i,j) + shetfx(i)=shetfx(i) + de_gcont*rx(i,j) + shetfy(i)=shetfy(i) + de_gcont*ry(i,j) + shetfz(i)=shetfz(i) + de_gcont*rz(i,j) + + shetfx(j)=shetfx(j) - de_gcont*rx(i,j) + shetfy(j)=shetfy(j) - de_gcont*ry(i,j) + shetfz(j)=shetfz(j) - de_gcont*rz(i,j) + +c energy correction from gcont + vbet(i,j)=vbet(i,j)*e_gcont + vbetap(i,j)=vbetap(i,j)*e_gcont + vbetap1(i,j)=vbetap1(i,j)*e_gcont + vbetap2(i,j)=vbetap2(i,j)*e_gcont + vbetam(i,j)=vbetam(i,j)*e_gcont + vbetam1(i,j)=vbetam1(i,j)*e_gcont + vbetam2(i,j)=vbetam2(i,j)*e_gcont + endif + ci elseif(istrand(i,j).eq.0)then ci vbet(i,j)=0 @@ -2434,7 +2459,16 @@ c write(*,*) 'vbetp:',vbetp c write(*,*) 'vbetm:',vbetm c write(*,*) 'vbet(i,j):',vbet(i,j) c stop - + + else + vbetap(i,j)=0 + vbetap1(i,j)=0 + vbetap2(i,j)=0 + vbetam(i,j)=0 + vbetam1(i,j)=0 + vbetam2(i,j)=0 + vbet(i,j)=0 + endif enddo enddo @@ -2458,6 +2492,9 @@ c------------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbet(maxca,maxca) real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) @@ -2516,6 +2553,7 @@ c local variables real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2 + real*8 e_gcont,fprim_gcont C-------------------------------------------------------------------------------- do i=4,inb-4 im3=i-3 @@ -2755,6 +2793,9 @@ c---------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -2795,11 +2836,17 @@ c local variables real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b + real*8 e_gcont,fprim_gcont c******************************************************************************** do i=3,inb-5 imm=i-2 im=i-1 do j=i+2,inb-3 + + if (dis(imm,j).lt.dfa_cutoff) then + call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + jp=j+1 jpp=j+2 @@ -2923,7 +2970,8 @@ ci & .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j) ! shefz(i,5)=shefz(i,5) ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j) - + + endif ci endif enddo @@ -2936,6 +2984,9 @@ c--------------------------------------------------------------------------c implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -2976,11 +3027,17 @@ C local variables real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b + real*8 e_gcont,fprim_gcont C******************************************************************************** do i=2,inb-6 ip=i+1 im=i-1 do j=i+3,inb-3 + + if (dis(im,j).lt.dfa_cutoff) then + call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + jp=j+1 jpp=j+2 @@ -3102,7 +3159,8 @@ ci & .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j) ! shefz(i,6)=shefz(i,6) ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j) - + + endif ci endif enddo @@ -3115,6 +3173,9 @@ c----------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -3155,12 +3216,18 @@ C local variables real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8 + real*8 e_gcont,fprim_gcont C******************************************************************************** do j=7,inb-1 jm=j-1 jmm=j-2 do i=1,j-6 + + if (dis(i,jmm).lt.dfa_cutoff) then + call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + ip=i+1 ipp=i+2 @@ -3284,6 +3351,7 @@ ci & .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then ! shefz(j,11)=shefz(j,11) ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm) + endif ci endif enddo @@ -3296,6 +3364,9 @@ c----------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -3335,11 +3406,17 @@ C local variables real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8 + real*8 e_gcont,fprim_gcont !c*************************************************************************c do j=6,inb-2 jp=j+1 jm=j-1 do i=1,j-5 + + if (dis(i,jm).lt.dfa_cutoff) then + call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + ip=i+1 ipp=i+2 @@ -3456,6 +3533,8 @@ ci & .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm) $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm) + endif + ci endif ENDDO diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F index a3e2b12..13267da 100644 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@ -3201,6 +3201,7 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d dij=dist(i,j) c write (iout,*) "dij(",i,j,") =",dij do k=1,constr_homology + if(.not.l_homo(k,ii)) cycle distance(k)=odl(k,ii)-dij c write (iout,*) "distance(",k,") =",distance(k) c @@ -3220,7 +3221,17 @@ c endif enddo - min_odl=minval(distancek) +c min_odl=minval(distancek) + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + & min_odl=distancek(kk) + enddo c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij @@ -3233,6 +3244,7 @@ c write (iout,* )"min_odl",min_odl c Nie wiem po co to liczycie jeszcze raz! c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ c & (2*(sigma_odl(i,j,k))**2)) + if(.not.l_homo(k,ii)) cycle if (waga_dist.ge.0.0d0) then c c For Gaussian-type Urestr @@ -3283,6 +3295,7 @@ c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) c & *waga_dist)+min_odl c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist c + if(.not.l_homo(k,ii)) cycle if (waga_dist.ge.0.0d0) then c For Gaussian-type Urestr c diff --git a/source/wham/src/make_ensemble1.F b/source/wham/src/make_ensemble1.F index e7e68de..fd548fc 100644 --- a/source/wham/src/make_ensemble1.F +++ b/source/wham/src/make_ensemble1.F @@ -25,7 +25,8 @@ double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors, & escloc,ehomology_constr, & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, - & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt + & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, + & edfadis,edfator,edfanei,edfabet integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist double precision qfree,sumprob,eini,efree,rmsdev character*80 bxname @@ -163,6 +164,10 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft esccor=enetb(19,i,iparm) edihcnstr=enetb(20,i,iparm) ehomology_constr=enetb(22,i,iparm) + edfadis=enetb(23,i,iparm) + edfator=enetb(24,i,iparm) + edfanei=enetb(25,i,iparm) + edfabet=enetb(26,i,iparm) if (homol_nset.gt.1) & ehomology_constr=waga_homology(homol_nset)*ehomology_constr #ifdef SPLITELE @@ -175,6 +180,8 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor & +wbond*estr+ehomology_constr + & +wdfa_dist*edfadis + & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet #else etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) @@ -185,6 +192,8 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor & +wbond*estr+ehomology_constr + & +wdfa_dist*edfadis + & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet #endif #ifdef MPI Fdimless_(i)= diff --git a/source/wham/src/molread_zs.F b/source/wham/src/molread_zs.F index 68d0723..b11e3fb 100644 --- a/source/wham/src/molread_zs.F +++ b/source/wham/src/molread_zs.F @@ -523,14 +523,18 @@ 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 + integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp + integer idomain(max_template,maxres) logical lprn /.true./ + integer ilen + external ilen logical unres_pdb 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 tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -546,8 +550,9 @@ c Alternative: reading from input 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 reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) + call readi(controlcard,"HOMOL_SET",homol_nset,1) + read2sigma=(index(controlcard,'READ2SIGMA').gt.0) call readi(controlcard,"IHSET",ihset,1) if (homol_nset.gt.1)then call card_concat(controlcard,.true.) @@ -579,6 +584,11 @@ c c c Reading HM global scores (prob not required) c + do i = nnt,nct + do k=1,constr_homology + idomain(k,i)=0 + enddo + enddo c open (4,file="HMscore") c do k=1,constr_homology c read (4,*,end=521) hmscore_tmp @@ -587,6 +597,13 @@ c write(*,*) "Model", k, ":", hmscore(k) c enddo c521 continue + ii=0 + do i = nnt,nct-2 + do j=i+2,nct + ii=ii+1 + ii_in_use(ii)=0 + enddo + enddo c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d do k=1,constr_homology @@ -596,9 +613,12 @@ 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)) + write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file', + & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 - 33 write (iout,'(a)') 'Error opening PDB file.' + 33 write (iout,'(a,5x,a)') 'Error opening PDB file', + & pdbfile(:ilen(pdbfile)) stop 34 continue c print *,'Begin reading pdb data' @@ -609,10 +629,6 @@ 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 @@ -646,96 +662,79 @@ 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 + 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)=rescore_tmp + rescore2(k,i_tmp)=rescore2_tmp + 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) - 1301 continue + 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) -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) + 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 + do i = nnt,nct-2 + do j=i+2,nct + + 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) +c write (iout,*) k,i,j,distal,dist2_cut + + if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 + & .and. distal.le.dist2_cut ) 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 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 + 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) + else #ifdef OLDSIGMA - sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else - sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & 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 + endif + sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) + else + ii=ii+1 + l_homo(k,ii)=.false. 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 + lim_odl=ii endif c c Theta, dihedral and SC retraints @@ -750,10 +749,11 @@ 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. + do i = nnt+3,nct + 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) @@ -764,7 +764,7 @@ 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))/4.0 -c +c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 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)* @@ -772,9 +772,8 @@ 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 + lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then @@ -790,6 +789,10 @@ 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. + 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), @@ -798,14 +801,15 @@ 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))/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)) 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') @@ -815,12 +819,15 @@ c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use 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) @@ -835,12 +842,36 @@ 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, @@ -855,31 +886,28 @@ 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) + write(iout,'(3i5,100(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 - write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i), + write (iout,'(i5,100(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), + write (iout,'(i5,100(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, + 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) enddo endif -#ifdef AIX - call flush_(iout) -#else - call flush(iout) -#endif c ----------------------------------------------------------------- return end