Adam's corrections

diff --git a/source/cluster/wham/src-HCD-5D/COMMON.CHAIN b/source/cluster/wham/src-HCD-5D/COMMON.CHAIN
index 9de64dd..2b481a5 100644 (file)
--- a/source/cluster/wham/src-HCD-5D/COMMON.CHAIN
+++ b/source/cluster/wham/src-HCD-5D/COMMON.CHAIN
@@ -1,6 +1,6 @@
       integer nres,nres0,nsup,nstart_sup,nend_sup,nstart_seq,
      & nchain,chain_border,chain_length,ireschain,npermchain,
-     & tabpermchain,ishift_pdb,iz_sc
+     & tabpermchain,ishift_pdb,iz_sc,nres_chomo
       double precision c,cref,crefjlee,cref_pdb,dc,xloc,xrot,dc_norm,
      & t,r,prod,rt,chomo
       common /chain/ c(3,maxres2+2),dc(3,maxres2),xloc(3,maxres),
@@ -18,4 +18,5 @@
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
       common /box/  boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
-      common /chomo_models/ chomo(3,maxres2+2,max_template)
+      common /chomo_models/ chomo(3,maxres2+2,max_template),
+     & nres_chomo(max_template)

diff --git a/source/cluster/wham/src-HCD-5D/energy_p_new.F b/source/cluster/wham/src-HCD-5D/energy_p_new.F
index 27e944b..5d07d5d 100644 (file)
--- a/source/cluster/wham/src-HCD-5D/energy_p_new.F
+++ b/source/cluster/wham/src-HCD-5D/energy_p_new.F
@@ -10203,16 +10203,21 @@ c
          enddo
          
 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) 
+         if (nexl.gt.0) then
+           min_odl=0.0d0
+         else
+           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
+           enddo
+         endif
+
 c        write (iout,* )"min_odl",min_odl
 #ifdef DEBUG
          write (iout,*) "ij dij",i,j,dij

diff --git a/source/cluster/wham/src-HCD-5D/read_constr_homology.F b/source/cluster/wham/src-HCD-5D/read_constr_homology.F
index b188deb..6ae3ef4 100644 (file)
--- a/source/cluster/wham/src-HCD-5D/read_constr_homology.F
+++ b/source/cluster/wham/src-HCD-5D/read_constr_homology.F
@@ -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 

diff --git a/source/cluster/wham/src-HCD-5D/readpdb.F b/source/cluster/wham/src-HCD-5D/readpdb.F
index a7bb52b..58c63e4 100644 (file)
--- a/source/cluster/wham/src-HCD-5D/readpdb.F
+++ b/source/cluster/wham/src-HCD-5D/readpdb.F
@@ -98,7 +98,7 @@ c          write (2,'(a)') card
 !              write (iout,*) "Calculating sidechain center iii",iii
               if (unres_pdb) then
                 do j=1,3
-                  dc(j,ires+nres)=sccor(j,iii)
+                  dc(j,ires_old)=sccor(j,iii)
                 enddo
               else
                 call sccenter(ires_old,iii,sccor)
@@ -347,17 +347,20 @@ c      character*5 atom
      & '     Gamma'
         endif
        endif
-      do i=1,nres-1
+      do i=2,nres
         iti=itype(i)
-        if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
+c        write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
+        if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
+     &    (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
-ctest          stop
+c          stop
         endif
-        vbld(i+1)=dist(i,i+1)
-        vbld_inv(i+1)=1.0d0/vbld(i+1)
-        if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
+        vbld(i)=dist(i-1,i)
+        vbld_inv(i)=1.0d0/vbld(i)
+        theta(i+1)=alpha(i-1,i,i+1)
         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
       enddo
+
 c      if (unres_pdb) then
 c        if (itype(1).eq.ntyp1) then
 c          theta(3)=90.0d0*deg2rad
@@ -806,7 +809,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
         endif
       endif
 C Calculate internal coordinates.
-      if (lprn) then
+      if (out_template_coord) then
       write (iout,'(/a)') 
      &  "Cartesian coordinates of the reference structure"
       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
@@ -818,6 +821,7 @@ C Calculate internal coordinates.
       enddo
       endif
 C Calculate internal coordinates.
+#ifdef DEBUG
        write (iout,'(a)') 
      &   "Backbone and SC coordinates as read from the PDB"
        do ires=1,nres
@@ -825,7 +829,8 @@ C Calculate internal coordinates.
      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
      &    (c(j,nres+ires),j=1,3)
        enddo
-      call int_from_cart(.true.,.false.)
+#endif
+      call int_from_cart(.true.,out_template_coord)
       call sc_loc_geom(.false.)
       do i=1,nres
         thetaref(i)=theta(i)

diff --git a/source/unres/src-HCD-5D/MREMD.F b/source/unres/src-HCD-5D/MREMD.F
index 38db8a8..9191402 100644 (file)
--- a/source/unres/src-HCD-5D/MREMD.F
+++ b/source/unres/src-HCD-5D/MREMD.F
@@ -62,6 +62,7 @@ cold      integer nup(0:maxprocs),ndown(0:maxprocs)
      & ene_i_iex,xxx,tmp,econstr_temp_iex,econstr_temp_i
       integer iran_num
       double precision ran_number
+      integer i_econstr/20/
 
 cdeb      imin_itime_old=0
       ntwx_cache=0
@@ -186,8 +187,9 @@ cd           write (*,*) me," After broadcast: file_exist",file_exist
              do i=nnt,nct
               if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
      &                   *dsqrt(gamsc(iabs(itype(i))))
-             enddo
+              enddo
            endif
+
        endif
 
        if(me.eq.king) then
@@ -250,20 +252,20 @@ cd           write (*,*) me," After broadcast: file_exist",file_exist
               write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
      &                    (ndowna(i,il),i=1,ndowna(0,il))
              enddo
+
              stdfp=dsqrt(2*Rb*t_bath/d_time)
              do i=1,ntyp
                stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
              enddo
              if (lang.gt.0 .and. .not.surfarea) then
-               do i=nnt,nct-1
-                 stdforcp(i)=stdfp*dsqrt(gamp)
-               enddo
+             do i=nnt,nct-1
+               stdforcp(i)=stdfp*dsqrt(gamp)
+             enddo
              do i=nnt,nct
               if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
      &                   *dsqrt(gamsc(iabs(itype(i))))
-               enddo
-             endif
-
+             enddo
+         endif
         ELSE IF (.not.(rest.and.file_exist)) THEN
          do il=1,remd_m(1)
           ifirst(il)=il
@@ -319,7 +321,7 @@ cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
           iset=i2set(me)
           if(me.eq.king.or..not.out1file) 
      &     write(iout,*) me,"iset=",iset,"t_bath=",t_bath
-          call flush(iout)
+c          call flush(iout)
        endif        
 c
        stdfp=dsqrt(2*Rb*t_bath/d_time)
@@ -721,7 +723,7 @@ ctime           call flush(iout)
                write(iout,*) "MIN ii_write=",ii_write
              endif
             endif
-           call flush(iout)
+c           call flush(iout)
         endif
         if (synflag) then
 c Update the time safety limiy
@@ -855,7 +857,7 @@ cd end
           time05=MPI_WTIME()
           if (me.eq.king .or. .not. out1file) then
             write(iout,*) 'REMD writing traj time=',time05-time04
-            call flush(iout)
+c            call flush(iout)
           endif
 
 
@@ -1089,12 +1091,12 @@ cd            write(iout,*) "i_dir=",i_dir
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned 
 c       on failed replica exchange attempt
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
                if (adaptive) then
-                 remd_ene(20,i)=remd_ene(n_ene+5,i)
-                 remd_ene(20,iex)=remd_ene(n_ene+6,iex)
+                 remd_ene(i_econstr,i)=remd_ene(n_ene+5,i)
+                 remd_ene(i_econstr,iex)=remd_ene(n_ene+6,iex)
                endif
             elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
 
@@ -1103,10 +1105,10 @@ c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
                i_iset1=i_iset+1
                i_mset1=iran_num(1,mset(i_iset1))
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
-               remd_ene(20,i)=remd_ene(n_ene+3,i)
-               remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
+               remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+               remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
 
             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
 
@@ -1115,14 +1117,14 @@ c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
                i_iset1=i_iset+1
                i_mset1=iran_num(1,mset(i_iset1))
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
                if (adaptive) then
-                 remd_ene(20,i)=remd_ene(n_ene+7,i)
-                 remd_ene(20,iex)=remd_ene(n_ene+8,iex)
+                 remd_ene(i_econstr,i)=remd_ene(n_ene+7,i)
+                 remd_ene(i_econstr,iex)=remd_ene(n_ene+8,iex)
                else
-                 remd_ene(20,i)=remd_ene(n_ene+3,i)
-                 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+                 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+                 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
                endif
             else
                goto 444 
@@ -1222,8 +1224,8 @@ cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
               else
                remd_ene(0,iex)=ene_iex_iex
                remd_ene(0,i)=ene_i_i
-               remd_ene(20,iex)=econstr_temp_iex
-               remd_ene(20,i)=econstr_temp_i
+               remd_ene(i_econstr,iex)=econstr_temp_iex
+               remd_ene(i_econstr,i)=econstr_temp_i
               endif
 
 cd      do il=1,nset
@@ -1259,7 +1261,7 @@ c-------------------------------------
               enddo
              endif
 
-             call flush(iout)
+c             call flush(iout)
 
 cd              write (iout,'(a6,100i4)') "ifirst",
 cd     &                    (ifirst(i),i=1,remd_m(1))
@@ -1305,7 +1307,6 @@ co     &    " rescaling weights with temperature",t_bath
          do i=1,ntyp
            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
          enddo
-
 c Compute the standard deviations of stochastic forces for Langevin dynamics
 c if the friction coefficients do not depend on surface area
          if (lang.gt.0 .and. .not.surfarea) then
@@ -1317,11 +1318,11 @@ c if the friction coefficients do not depend on surface area
      &                   *dsqrt(gamsc(iabs(itype(i))))
            enddo
          endif
-cde         write(iout,*) 'REMD after',me,t_bath
+cde           write(iout,*) 'REMD after',me,t_bath
            time08=MPI_WTIME()
            if (me.eq.king .or. .not. out1file) then
             write(iout,*) 'REMD exchange time=',time08-time00
-            call flush(iout)
+c            call flush(iout)
            endif
         endif
       enddo

diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F
index 5dd78ea..16c0f37 100644 (file)
--- a/source/unres/src-HCD-5D/readrtns_CSA.F
+++ b/source/unres/src-HCD-5D/readrtns_CSA.F
@@ -3559,13 +3559,12 @@ c----------------------------------------------------------------------
       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,
-     & ik,ll,ii,kk,iistart,lim_xx
+     & ik,ll,ii,kk,iistart,iishift,lim_xx
       double precision distal
       logical lprn /.true./
       integer ilen
       external ilen
-      logical liiflag,lfirst
-      integer i01,i10
+      logical liiflag
 c
 c
       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
@@ -3756,39 +3755,36 @@ c
       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.or.
-     &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
+          if (ii_in_use(ii).eq.0.and.liiflag) then
             liiflag=.false.
-            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
+            iistart=ii
           endif
-          if (ii_in_use(ii).ne.0.and..not.liiflag) then
-             i01=ii
+          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
-        lim_odl=iistart-1
       endif
 
       return

diff --git a/source/wham/src-HCD-5D/COMMON.CHAIN b/source/wham/src-HCD-5D/COMMON.CHAIN
index dfffc78..7b79a58 100644 (file)
--- a/source/wham/src-HCD-5D/COMMON.CHAIN
+++ b/source/wham/src-HCD-5D/COMMON.CHAIN
@@ -1,6 +1,6 @@
       integer nres,nres0,nsup,nstart_sup,nend_sup,nstart_seq,
      & ishift_pdb,chain_length,chain_border,chain_border1,ichanres,
-     & tabpermchain,nchain ,npermchain,ireschain,iz_sc
+     & tabpermchain,nchain ,npermchain,ireschain,iz_sc,nres_chomo
       double precision c,cref,crefjlee,dc,xloc,xrot,dc_norm,t,r,prod,rt,
      & rmssing,anatemp,chomo
       common /chain/ c(3,maxres2+2),dc(3,maxres2),xloc(3,maxres),
@@ -18,4 +18,5 @@
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
       common /box/  boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
-      common /chomo_models/ chomo(3,maxres2+2,max_template)
+      common /chomo_models/ chomo(3,maxres2+2,max_template),
+     & nres_chomo(max_template)

diff --git a/source/wham/src-HCD-5D/energy_p_new.F b/source/wham/src-HCD-5D/energy_p_new.F
index bd69774..6105156 100644 (file)
--- a/source/wham/src-HCD-5D/energy_p_new.F
+++ b/source/wham/src-HCD-5D/energy_p_new.F
@@ -4609,16 +4609,20 @@ c
          enddo
          
 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) 
+         if (nexl.gt.0) then
+           min_odl=0.0d0
+         else
+           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
+           enddo
+         endif
 c        write (iout,* )"min_odl",min_odl
 #ifdef DEBUG
          write (iout,*) "ij dij",i,j,dij

diff --git a/source/wham/src-HCD-5D/read_constr_homology.F b/source/wham/src-HCD-5D/read_constr_homology.F
index dece50f..fa81b80 100644 (file)
--- a/source/wham/src-HCD-5D/read_constr_homology.F
+++ b/source/wham/src-HCD-5D/read_constr_homology.F
@@ -1,5 +1,6 @@
       subroutine read_constr_homology
 
+      implicit none
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
       include 'DIMENSIONS.FREE'
@@ -29,12 +30,15 @@ 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
@@ -143,6 +147,7 @@ c
         tpl_k_rescore="template"//kic2//".sco"
 
         unres_pdb=.false.
+        nres_temp=nres
         if (read2sigma) then
           call readpdb_template(k)
           close(ipdbin)
@@ -150,27 +155,29 @@ c
           call readpdb
           close(ipdbin)
         endif
+        nres_chomo(k)=nres
+        nres=nres_temp
 
-c        call readpdb
-        do i=1,2*nres
+        do i=1,2*nres_chomo(k)
           do j=1,3
             crefjlee(j,i)=c(j,i)
           enddo
         enddo
 #ifdef DEBUG
-        do i=1,nres
+        do i=1,nres_chmo(k)
           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
      &      (crefjlee(j,i+nres),j=1,3)
         enddo
         write (iout,*) "read_constr_homology: after reading pdb file"
         call flush(iout)
 #endif
+
 c
 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 +266,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 +303,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 +345,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 +372,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 +382,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"
@@ -454,7 +479,7 @@ c -----------------------------------------------------------------
       end
 c----------------------------------------------------------------------
       subroutine read_klapaucjusz
-
+      implicit none
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
       include 'DIMENSIONS.FREE'
@@ -477,8 +502,11 @@ c----------------------------------------------------------------------
       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 i,ii,ik,ki,j,k,l,ll,ii_in_use(maxdim),i_tmp,idomain_tmp,
+     & nclust,iistart,iishift,lim_xx
+      integer nres_temp
       integer idomain(max_template,maxres)
+      double precision distal
       logical lprn /.true./
       integer ilen
       external ilen
@@ -498,6 +526,7 @@ c
       call getenv("FRAGFILE",fragfile) 
       open(ientin,file=fragfile,status="old",err=10)
       read(ientin,*) constr_homology,nclust
+      lim_xx=0
       l_homo = .false.
       sigma_theta=0.0
       sigma_d=0.0
@@ -514,7 +543,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)
@@ -535,6 +567,7 @@ c
 c Loop over clusters
 c
       do l=1,nclust
+c        write (iout,*) "CLUSTER",l," members:",ninclust(l)
         do ll = 1,ninclust(l)
         
         k = inclust(ll,l)
@@ -553,18 +586,27 @@ c     Distance restraints
 c
 c          ... --> odl(k,ii)
 C Copy the coordinates from reference coordinates (?)
+c        write (iout,*) "k",k
+        nres_temp=nres
+        nres=nres_chomo(k)
         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)
+c            write (iout,*) "c(",j,i,") =",c(j,i)
           enddo
         enddo
+c        write (iout,*) "read_klapaucjusz: calling int_from_cart"
+c        call cartprint
+c        write (iout,*) "idomain"
+c        write (iout,'(2i5)') (i,idomain(k,i),i=1,nres)
         call int_from_cart(.true.,.false.)
         call sc_loc_geom(.false.)
+c        write (iout,*) "en" 
         do i=1,nres
           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 

diff --git a/source/wham/src-HCD-5D/readpdb.F b/source/wham/src-HCD-5D/readpdb.F
index 83c63ab..6f4ba5f 100644 (file)
--- a/source/wham/src-HCD-5D/readpdb.F
+++ b/source/wham/src-HCD-5D/readpdb.F
@@ -27,10 +27,9 @@ C geometry.
       ibeg=1
       ishift1=0
       sccalc=.false.
-      ires=0
       do
         read (ipdbin,'(a80)',end=10) card
-c        write (iout,'(a)') card
+!       write (iout,'(a)') card
         if (card(:5).eq.'HELIX') then
           nhfrag=nhfrag+1
           lsecondary=.true.
@@ -59,7 +58,7 @@ c        write (iout,'(a)') card
           iterter(ires_old)=1
           ishift1=ishift1+1
           ibeg=2
-          write (iout,*) "Chain ended",ires,ishift,ires_old
+!          write (iout,*) "Chain ended",ires,ishift,ires_old
           if (unres_pdb) then
             do j=1,3
               dc(j,ires)=sccor(j,iii)
@@ -92,7 +91,7 @@ c          write (2,'(a)') card
 !              write (iout,*) "Calculating sidechain center iii",iii
               if (unres_pdb) then
                 do j=1,3
-                  dc(j,ires+nres)=sccor(j,iii)
+                  dc(j,ires_old)=sccor(j,iii)
                 enddo
               else
                 call sccenter(ires_old,iii,sccor)
@@ -120,7 +119,7 @@ c              write (iout,*) "BEG ires",ires
 ! Start a new chain
               ishift=-ires_old+ires-1 !!!!!
               ishift1=ishift1-1    !!!!!
-              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
+!              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
               ires=ires-ishift+ishift1
               ires_old=ires
               ibeg=0
@@ -326,7 +325,7 @@ c        write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
         if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
      &    (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
-          stop
+c          stop
         endif
         vbld(i)=dist(i-1,i)
         vbld_inv(i)=1.0d0/vbld(i)
@@ -805,8 +804,8 @@ C Calculate internal coordinates.
       endif
 C Calculate internal coordinates.
 c      call int_from_cart1(.false.)
-      call int_from_cart(.true.,.true.)
-      call sc_loc_geom(.true.)
+      call int_from_cart(.true.,out_template_coord)
+      call sc_loc_geom(.false.)
       do i=1,nres
         thetaref(i)=theta(i)
         phiref(i)=phi(i)