READ_HOMOL_FRAG
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Mon, 26 Mar 2018 01:08:24 +0000 (03:08 +0200)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Mon, 26 Mar 2018 01:08:24 +0000 (03:08 +0200)
source/unres/src_MD/COMMON.CONTROL
source/unres/src_MD/DIMENSIONS
source/unres/src_MD/readrtns.F

index eeb46be..4dbaa2a 100644 (file)
@@ -8,7 +8,7 @@
      &                 sideadd,lsecondary,read_cart,unres_pdb,
      &                 vdisulf,searchsc,lmuca,dccart,extconf,out1file,
      &                 gnorm_check,gradout,split_ene,read2sigma
-     &                ,start_from_model
+     &                ,start_from_model,read_homol_frag
       common /cntrl/ checkgrad_inc,
      & modecalc,iscode,indpdb,indback,indphi,iranconf,
      & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,outx,
@@ -16,7 +16,7 @@
      & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb
      & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file,
      & constr_dist,gnorm_check,gradout,split_ene,constr_homology,
-     & homol_nset,read2sigma,start_from_model
+     & homol_nset,read2sigma,start_from_model,read_homol_frag
       common /homol/ waga_homology(maxprocs/20),
      & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,dist2_cut
 C... minim = .true. means DO minimization.
index 117b9fd..9d4b545 100644 (file)
@@ -140,3 +140,6 @@ C Nose-Hoover chain - chain length and order of Yoshida algorithm
 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)
index 93ed3a4..7fa2ea1 100644 (file)
@@ -98,6 +98,7 @@ C Set up the time limit (caution! The time must be input in minutes!)
       read_cart=index(controlcard,'READ_CART').gt.0
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
+      read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
@@ -2825,6 +2826,10 @@ c
         enddo
       enddo
 
+      if (read_homol_frag) then
+        call read_klapaucjusz
+      else
+
       do k=1,constr_homology
 
         read(inp,'(a)') pdbfile
@@ -3097,6 +3102,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,
@@ -3138,6 +3146,258 @@ cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 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.MD'
+      include 'COMMON.GEO'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      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
+      logical lprn /.true./
+      integer ilen
+      external ilen
+      logical 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 pdbfile,tpl_k_rescore
+
+c
+c For new homol impl
+c
+      include 'COMMON.VAR'
+c
+      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_template(k)
+        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----------------------------------------------------------------------
 
 #ifdef WINIFL