cluster & wham update
[unres.git] / source / cluster / wham / src-HCD / read_constr_homology.F
index b188deb..6ae3ef4 100644 (file)
@@ -1,5 +1,5 @@
       subroutine read_constr_homology
-
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
@@ -27,16 +27,19 @@ 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, ii_in_use(maxdim),i_tmp,idomain_tmp
+      integer ki,i,ii,ik,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,
+     & lim_theta,lim_xx,irec,iistart,iishift,i10,i01
+      double precision distal
       integer idomain(max_template,maxres)
-      logical lprn /.true./
+      logical lfirst
       integer ilen
       external ilen
       logical liiflag
+      integer nres_temp
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
-      double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+      double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
      &    rescore3_tmp
       double precision, dimension (max_template,maxres) :: rescore
       double precision, dimension (max_template,maxres) :: rescore2
@@ -142,6 +145,7 @@ c
         tpl_k_rescore="template"//kic2//".sco"
 
         unres_pdb=.false.
+        nres_temp=nres
         if (read2sigma) then
           call readpdb_template(k)
           close(ipdbin)
@@ -149,15 +153,16 @@ c
           call readpdb(out_template_coord)
           close(ipdbin)
         endif
+        nres_chomo(k)=nres
+        nres=nres_temp
 
-c        call readpdb
         do i=1,2*nres
           do j=1,3
             crefjlee(j,i)=c(j,i)
           enddo
         enddo
 #ifdef DEBUG
-        do i=1,nres
+        do i=1,nres_chomo(k)
           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
      &      (crefjlee(j,i+nres),j=1,3)
         enddo
@@ -170,7 +175,7 @@ c     Distance restraints
 c
 c          ... --> odl(k,ii)
 C Copy the coordinates from reference coordinates (?)
-        do i=1,2*nres
+        do i=1,2*nres_chomo(k)
           do j=1,3
             c(j,i)=cref(j,i)
 c           write (iout,*) "c(",j,i,") =",c(j,i)
@@ -259,6 +264,8 @@ c    &                       constr_homology
           enddo
         lim_odl=ii
         endif
+c        write (iout,*) "Distance restraints set"
+c        call flush(iout)
 c
 c     Theta, dihedral and SC retraints
 c
@@ -294,11 +301,13 @@ 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
             if (sigma_dih(k,i).ne.0)
-     &      sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+     &       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 
         endif
+c        write (iout,*) "Dihedral angle restraints set"
+c        call flush(iout)
 
         if (waga_theta.gt.0.0d0) then
 c         open (ientin,file=tpl_k_sigma_theta,status='old')
@@ -334,6 +343,8 @@ c                             rescore(k,i-2) !  right expression ?
 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
           enddo
         endif
+c        write (iout,*) "Angle restraints set"
+c        call flush(iout)
 
         if (waga_d.gt.0.0d0) then
 c       open (ientin,file=tpl_k_sigma_d,status='old')
@@ -359,7 +370,6 @@ c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 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)
-c               sigma_d(k,i)=rescore(k,i) !  right expression ?
                sigma_d(k,i)=rescore3(k,i) !  right expression ?
                if (sigma_d(k,i).ne.0)
      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
@@ -370,57 +380,70 @@ c              read (ientin,*) sigma_d(k,i) ! 1st variant
           enddo
         endif
       enddo
+c      write (iout,*) "SC restraints set"
+c      call flush(iout)
 c
 c remove distance restraints not used in any model from the list
 c shift data in all arrays
 c
+c      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
       if (waga_dist.ne.0.0d0) then
         ii=0
         liiflag=.true.
+        lfirst=.true.
         do i=nnt,nct-2 
          do j=i+2,nct 
           ii=ii+1
-          if (ii_in_use(ii).eq.0.and.liiflag) then
+c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+c     &            .and. distal.le.dist2_cut ) then
+c          write (iout,*) "i",i," j",j," ii",ii
+c          call flush(iout)
+          if (ii_in_use(ii).eq.0.and.liiflag.or.
+     &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
             liiflag=.false.
-            iistart=ii
+            i10=ii
+            if (lfirst) then
+              lfirst=.false.
+              iistart=ii
+            else
+              if(i10.eq.lim_odl) i10=i10+1
+              do ki=0,i10-i01-1
+               ires_homo(iistart+ki)=ires_homo(ki+i01)
+               jres_homo(iistart+ki)=jres_homo(ki+i01)
+               ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+               do k=1,constr_homology
+                odl(k,iistart+ki)=odl(k,ki+i01)
+                sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+                l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+               enddo
+              enddo
+              iistart=iistart+i10-i01
+            endif
           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
+          if (ii_in_use(ii).ne.0.and..not.liiflag) then
+             i01=ii
              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
+        lim_odl=iistart-1
       endif
-
-      endif ! .not. klapaucjusz     
+c      write (iout,*) "Removing distances completed"
+c      call flush(iout)
+      endif ! .not. klapaucjusz
 
       if (constr_homology.gt.0) call homology_partition
+c      write (iout,*) "After homology_partition"
+c      call flush(iout)
       if (constr_homology.gt.0) call init_int_table
-cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
-cd     & "lim_xx=",lim_xx
-c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
-c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c      write (iout,*) "After init_int_table"
+c      call flush(iout)
+c      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 c
 c Print restraints
 c
-      if (.not.lprn) return
+      if (.not.out_template_restr) return
 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 c      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
        write (iout,*) "Distance restraints from templates"
@@ -477,7 +500,9 @@ c----------------------------------------------------------------------
       character*500 controlcard
       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
       integer idomain(max_template,maxres)
+      integer nres_temp
       logical lprn /.true./
+      logical lfirst
       integer ilen
       external ilen
       logical liiflag
@@ -513,7 +538,10 @@ c Read pdb files
         stop
   34    continue
         unres_pdb=.false.
+        nres_temp=nres
         call readpdb_template(k)
+        nres_chomo(k)=nres
+        nres=nres_temp
 c        do i=1,2*nres
 c          do j=1,3
 c            chomo(j,i,k)=c(j,i)
@@ -552,6 +580,8 @@ c     Distance restraints
 c
 c          ... --> odl(k,ii)
 C Copy the coordinates from reference coordinates (?)
+        nres_temp=nres
+        nres=nres_chomo(k)
         do i=1,2*nres
           do j=1,3
             c(j,i)=chomo(j,i,k)
@@ -564,6 +594,7 @@ c           write (iout,*) "c(",j,i,") =",c(j,i)
           thetaref(i)=theta(i)
           phiref(i)=phi(i)
         enddo
+        nres=nres_temp
         if (waga_dist.ne.0.0d0) then
           ii=0
           do i = nnt,nct-2