& with_dihed_constr,with_theta_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, with_dihed_constr,with_theta_constr,
& efree,iopt,nstart,nend,symetr,
& tor_mode,shield_mode,
& constr_dist,out1file,
- & constr_homology,homol_nset,read2sigma
+ & constr_homology,homol_nset,read2sigma,read_homol_frag
common /homol/ waga_homology(10),
& waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut,
& iset,ihset,l_homo(max_template,maxdim),
C Maximum number of bins in SAXS restraints
integer MaxSAXS
parameter (MaxSAXS=1000)
+c Maximum number of clusters of templates containing same fragments
+ integer maxclust
+ parameter(maxclust=1000)
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
call readi(controlcard,'NSAXS',nsaxs,0)
call readi(controlcard,'SAXS_MODE',saxs_mode,0)
call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
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
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,
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)=0.2d0
+ 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----------------------------------------------------------------------
logical refstr,pdbref,punch_dist,print_rms,caonly,verbose,
& merge_helices,bxfile,cxfile,histfile,entfile,zscfile,
& rmsrgymap,with_dihed_constr,check_conf,histout,out1file,
- & read2sigma,l_homo,with_theta_constr
+ & read2sigma,l_homo,with_theta_constr,read_homol_frag
common /cntrl/ iscode,indpdb,refstr,pdbref,outpdb,outmol2,
& punch_dist,print_rms,caonly,verbose,icomparfunc,pdbint,
& merge_helices,bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap,
& with_theta_constr,
& symetr,
& constr_dist,
- & constr_homology,out1file,homol_nset,read2sigma
+ & constr_homology,out1file,homol_nset,read2sigma,read_homol_frag
common /homol/ waga_homology(maxR),
& waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut,
& iset,ihset,l_homo(max_template,maxdim)
C Maximum number of bins in SAXS restraints
integer MaxSAXS
parameter (MaxSAXS=1000)
+c Maximum number of clusters of templates containing same fragments
+ integer maxclust
+ parameter(maxclust=1000)
+
enddo
c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+ if (read_homol_frag) then
+ call read_klapaucjusz
+ else
+
do k=1,constr_homology
read(inp,'(a)') pdbfile
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,
return
end
c----------------------------------------------------------------------
+ subroutine read_klapaucjusz
+
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
+#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)=0.2d0
+ 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----------------------------------------------------------------------
call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
write (iout,*) "with_homology_constr ",with_dihed_constr,
& " CONSTR_HOMOLOGY",constr_homology
+ read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
refstr = index(controlcard,'REFSTR').gt.0
pdbref = index(controlcard,'PDBREF').gt.0
dyn_ss=(index(controlcard,'DYN_SS').gt.0)
projects / unres.git / commitdiff