From eeac8f8e7683a0ff43598f9b7dee32636646c83b Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Mon, 26 Mar 2018 12:12:50 +0200 Subject: [PATCH] READ_HOMOL_FRAG cluster_wham --- source/cluster/wham/src/COMMON.CONTROL | 4 +- source/cluster/wham/src/DIMENSIONS | 3 + source/cluster/wham/src/readrtns.F | 266 ++++++++++++++++++++++++++++++++ 3 files changed, 271 insertions(+), 2 deletions(-) diff --git a/source/cluster/wham/src/COMMON.CONTROL b/source/cluster/wham/src/COMMON.CONTROL index 6978419..52db8c8 100644 --- a/source/cluster/wham/src/COMMON.CONTROL +++ b/source/cluster/wham/src/COMMON.CONTROL @@ -8,7 +8,7 @@ & 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 + & ,read2sigma,l_homo,read_homol_frag 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, @@ -17,5 +17,5 @@ common /cntrlr/ waga_homology(20), & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,dist2_cut, & iset, - & print_homology_restraints,print_homology_models + & print_homology_restraints,print_homology_models,read_homol_frag diff --git a/source/cluster/wham/src/DIMENSIONS b/source/cluster/wham/src/DIMENSIONS index d8cdbf4..554a1fe 100644 --- a/source/cluster/wham/src/DIMENSIONS +++ b/source/cluster/wham/src/DIMENSIONS @@ -68,3 +68,6 @@ C Maximum number of terms in SC bond-stretching potential C Maximum number of templates in homology-modeling restraints integer max_template parameter(max_template=25) +c Maximum number of clusters of templates containing same fragments + integer maxclust + parameter(maxclust=1000) diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F index fac2508..4115cca 100644 --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@ -76,6 +76,7 @@ C print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0 print_homology_models= & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0 + read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0 #ifdef AIX call flush_(iout) @@ -1032,6 +1033,9 @@ c521 continue enddo c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d + if (read_homol_frag) then + call read_klapaucjusz + else write (iout,*) "CONSTR_HOMOLOGY",constr_homology do k=1,constr_homology @@ -1327,6 +1331,9 @@ c enddo enddo endif + + endif ! .not. klapaucjusz + 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, @@ -1369,3 +1376,262 @@ c ----------------------------------------------------------------- return end c---------------------------------------------------------------------- + subroutine read_klapaucjusz + + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.HOMRESTR' + character*256 fragfile + integer ninclust(maxclust),inclust(max_template,maxclust), + & nresclust(maxclust),iresclust(maxres,maxclust) + + character*2 kic2 + character*24 model_ki_dist, model_ki_angle + character*500 controlcard + 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,liiflag +c +c + 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 For new homol impl +c + include 'COMMON.VAR' +c + double precision chomo(3,maxres2+2,max_template) + call getenv("FRAGFILE",fragfile) + open(ientin,file=fragfile,status="old",err=10) + read(ientin,*) constr_homology,nclust + l_homo = .false. + sigma_theta=0.0 + sigma_d=0.0 + sigma_dih=0.0 +c Read pdb files + do k=1,constr_homology + read(ientin,'(a)') pdbfile + write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', + & pdbfile(:ilen(pdbfile)) + open(ipdbin,file=pdbfile,status='old',err=33) + goto 34 + 33 write (iout,'(a,5x,a)') 'Error opening PDB file', + & pdbfile(:ilen(pdbfile)) + stop + 34 continue + unres_pdb=.false. + call readpdb + do i=1,2*nres + do j=1,3 + chomo(j,i,k)=c(j,i) + enddo + enddo + do i=1,nres + rescore(k,i)=1.0d0 + rescore2(k,i)=1.0d0 + enddo + enddo +c Read clusters + do i=1,nclust + read(ientin,*) ninclust(i),nresclust(i) + read(ientin,*) (inclust(k,i),k=1,ninclust(i)) + read(ientin,*) (iresclust(k,i),k=1,nresclust(i)) + enddo +c +c Loop over clusters +c + do l=1,nclust + do ll = 1,ninclust(l) + + k = inclust(ll,l) + do i=1,nres + idomain(k,i)=0 + enddo + do i=1,nresclust(l) + if (nnt.gt.1) then + idomain(k,iresclust(i,l)+1) = 1 + else + idomain(k,iresclust(i,l)) = 1 + endif + enddo +c +c Distance restraints +c +c ... --> odl(k,ii) +C Copy the coordinates from reference coordinates (?) + do i=1,2*nres + do j=1,3 + c(j,i)=chomo(j,i,k) +c write (iout,*) "c(",j,i,") =",c(j,i) + enddo + enddo + call int_from_cart(.true.,.false.) + call sc_loc_geom(.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + if (waga_dist.ne.0.0d0) then + ii=0 + 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 + odl(k,ii)=distal + 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))* + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) +#else + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) +#endif + endif + endif + sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) + else + ii=ii+1 +c l_homo(k,ii)=.false. + endif + enddo + enddo + lim_odl=ii + endif +c +c Theta, dihedral and SC retraints +c + if (waga_angle.gt.0.0d0) then + do i = nnt+3,nct + if (idomain(k,i).eq.0) then +c sigma_dih(k,i)=0.0 + cycle + endif + dih(k,i)=phiref(i) + sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2)+rescore(k,i-3))/4.0 +c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), +c & " sigma_dihed",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)) + enddo + lim_dih=nct-nnt-2 + endif + + if (waga_theta.gt.0.0d0) then + do i = nnt+2,nct + if (idomain(k,i).eq.0) then +c sigma_theta(k,i)=0.0 + cycle + endif + thetatpl(k,i)=thetaref(i) + sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2))/3.0 + if (sigma_theta(k,i).ne.0) + & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + enddo + endif + + if (waga_d.gt.0.0d0) then + do i = nnt,nct + if (itype(i).eq.10) cycle + if (idomain(k,i).eq.0 ) then +c sigma_d(k,i)=0.0 + cycle + endif + xxtpl(k,i)=xxref(i) + yytpl(k,i)=yyref(i) + zztpl(k,i)=zzref(i) + sigma_d(k,i)=rescore(k,i) + if (sigma_d(k,i).ne.0) + & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) + if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 + enddo + endif + enddo ! l + enddo ! ll +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 + liiflag=.true. + do i=nnt,nct-2 + do j=i+2,nct + ii=ii+1 + if (ii_in_use(ii).eq.0.and.liiflag) then + liiflag=.false. + iistart=ii + endif + if (ii_in_use(ii).ne.0.and..not.liiflag.or. + & .not.liiflag.and.ii.eq.lim_odl) then + if (ii.eq.lim_odl) then + iishift=ii-iistart+1 + else + iishift=ii-iistart + endif + liiflag=.true. + do ki=iistart,lim_odl-iishift + ires_homo(ki)=ires_homo(ki+iishift) + jres_homo(ki)=jres_homo(ki+iishift) + ii_in_use(ki)=ii_in_use(ki+iishift) + do k=1,constr_homology + odl(k,ki)=odl(k,ki+iishift) + sigma_odl(k,ki)=sigma_odl(k,ki+iishift) + l_homo(k,ki)=l_homo(k,ki+iishift) + enddo + enddo + ii=ii-iishift + lim_odl=lim_odl-iishift + endif + enddo + enddo + endif + + return + 10 stop "Error infragment file" + end +c---------------------------------------------------------------------- -- 1.7.9.5