Adam's changes to wham and cluster following previous commit

XG.. and ..GX sequence is automatically corrected

also in cluster_wham it is possible to select number of clusters
instead of cutoff

diff --git a/source/cluster/wham/src/COMMON.CLUSTER b/source/cluster/wham/src/COMMON.CLUSTER
index 4477d19..f1ad0fd 100644 (file)
--- a/source/cluster/wham/src/COMMON.CLUSTER
+++ b/source/cluster/wham/src/COMMON.CLUSTER
@@ -4,11 +4,11 @@
       real*4 diss,allcart
       double precision enetb,entfac,totfree,energy,rmstb
       integer ncut,ngr,licz,nconf,iass,icc,mult,list_conf,
-     &  nss_all,ihpb_all,jhpb_all,iass_tot,iscore,nprop
+     &  nss_all,ihpb_all,jhpb_all,iass_tot,iscore,nprop,nclust
       common /clu/ diss(maxdist),energy(0:maxconf),
      &  enetb(max_ene,maxstr_proc),ecut,
      &  entfac(maxconf),totfree(0:maxconf),totfree_gr(maxgr),
-     &  rcutoff(max_cut+1),ncut,min_var,tree,plot_tree,lgrp
+     &  rcutoff(max_cut+1),ncut,nclust,min_var,tree,plot_tree,lgrp
       common /clu1/ ngr,licz(maxgr),nconf(maxgr,maxingr),iass(maxgr),
      &  iass_tot(maxgr,max_cut),list_conf(maxconf)
       common /alles/ allcart(3,maxres2,maxstr_proc),rmstb(maxconf),

diff --git a/source/cluster/wham/src/COMMON.CONTROL b/source/cluster/wham/src/COMMON.CONTROL
index 2b9b9e5..fe7c63a 100644 (file)
--- a/source/cluster/wham/src/COMMON.CONTROL
+++ b/source/cluster/wham/src/COMMON.CONTROL
@@ -5,11 +5,14 @@
       real*8 waga_dist, waga_angle,waga_theta, waga_d, dist_cut
       logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
      & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx,
-     & with_dihed_constr,out1file
+     & with_dihed_constr,out1file,print_homology_restraints,
+     & print_contact_map,print_homology_models
       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,
-     & with_dihed_constr, constr_homology,homol_nset,out1file
+     & with_dihed_constr, constr_homology,homol_nset,out1file,
+     & print_contact_map
       common /cntrlr/ waga_homology(10),
-     & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,iset
+     & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,iset,
+     & print_homology_restraints,print_homology_models
 

diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F
index a5ba7fb..8ce7e5b 100644 (file)
--- a/source/cluster/wham/src/energy_p_new.F
+++ b/source/cluster/wham/src/energy_p_new.F
@@ -5166,11 +5166,9 @@ c   3 = SC...Ca...Ca...SCi
           cosphi=dcos(j*tauangle(intertyp,i))
           sinphi=dsin(j*tauangle(intertyp,i))
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
-#define DEBUG
 #ifdef DEBUG
           esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
 #endif
-#undef DEBUG
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci

diff --git a/source/cluster/wham/src/main_clust.F b/source/cluster/wham/src/main_clust.F
index c31847f..980871c 100644 (file)
--- a/source/cluster/wham/src/main_clust.F
+++ b/source/cluster/wham/src/main_clust.F
@@ -34,7 +34,7 @@ C
       DIMENSION NN(maxconf),DISNN(maxconf)
       LOGICAL FLAG(maxconf)
       integer i,j,k,l,m,n,len,lev,idum,ii,ind,ioffset,jj,icut,ncon,
-     & it,ncon_work,ind1,ilen
+     & it,ncon_work,ind1,ilen,is,ie
       double precision t1,t2,tcpu,difconf
       real diss_(maxdist)
       
@@ -76,6 +76,12 @@ c      if (refstr) call read_ref_structure(*30)
 
       print *,'MAIN: nnt=',nnt,' nct=',nct
 
+      if (nclust.gt.0) then
+        PRINTANG(1)=.TRUE.
+        PRINTPDB(1)=outpdb
+        printmol2(1)=outmol2
+        ncut=0
+      else
       DO I=1,NCUT
         PRINTANG(I)=.FALSE.
         PRINTPDB(I)=0
@@ -87,12 +93,21 @@ c      if (refstr) call read_ref_structure(*30)
           printmol2(i)=outmol2
         ENDIF
       ENDDO
+      endif
+      if (ncut.gt.0) then
       write (iout,*) 'Number of cutoffs:',NCUT
       write (iout,*) 'Cutoff values:'
       DO ICUT=1,NCUT
         WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
      &    printpdb(icut),printmol2(icut)
       ENDDO
+      else if (nclust.gt.0) then
+      write (iout,'("Number of clusters requested",i5)') nclust
+      else
+      if (me.eq.Master)
+     & write (iout,*) "ERROR: Either nclust or ncut must be >0"
+      stop
+      endif
       DO I=1,NRES-3  
         MULT(I)=1
       ENDDO
@@ -247,29 +262,39 @@ C
       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 
+c 3/3/16 AL: added explicit number of cluters
+      if (nclust.gt.0) then
+        is=nclust-1
+        ie=nclust-1
+        icut=1
+      else
+        is=1
+        ie=lev-1
+      endif
       do i=1,maxgr
         licz(i)=0
       enddo
       icut=1
-      i=1
-      NGR=i+1
+      i=is
+      NGR=is+1
       do j=1,n
         licz(iclass(j,i))=licz(iclass(j,i))+1
         nconf(iclass(j,i),licz(iclass(j,i)))=j
 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
       enddo        
-      do i=1,lev-1
-
+c      do i=1,lev-1
+      do i=is,ie
          idum=lev-i
          DO L=1,LEV
             IF (HEIGHT(L).EQ.IDUM) GOTO 190
          ENDDO
  190     IDUM=L
-         write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
-     &    " icut",icut," cutoff",rcutoff(icut)
-         IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
-          WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
+c         write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
+c     &    " icut",icut," cutoff",rcutoff(icut)
+         IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
+         if (nclust.le.0) 
+     &    WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
           write (iout,'(a,f8.2)') 'Maximum distance found:',
      &              CRITVAL(IDUM)
           CALL SRTCLUST(ICUT,ncon_work,iT)
@@ -282,9 +307,10 @@ c     &    nconf(iclass(j,i),licz(iclass(j,i)))
          do l=1,maxgr
           licz(l)=0
          enddo
+         ii=i-is+1
          do j=1,n
-          licz(iclass(j,i))=licz(iclass(j,i))+1
-          nconf(iclass(j,i),licz(iclass(j,i)))=j
+          licz(iclass(j,ii))=licz(iclass(j,ii))+1
+          nconf(iclass(j,ii),licz(iclass(j,ii)))=j
 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 cd          print *,j,iclass(j,i),

diff --git a/source/cluster/wham/src/probabl.F b/source/cluster/wham/src/probabl.F
index 2f34f61..62da77f 100644 (file)
--- a/source/cluster/wham/src/probabl.F
+++ b/source/cluster/wham/src/probabl.F
@@ -197,16 +197,20 @@ cc        if (wcorr6.eq.0) ecorr6=0.0d0
 #ifdef MPI
         Fdimless_(i)=beta_h(ib)*etot+entfac(ii)
         totfree_(i)=etot
+#ifdef DEBUG
+        write (iout,*) "etrop", i,ii,ib,
+     &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
+     &   entfac(ii),Fdimless_(i)
+#endif
 #else
         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
         totfree(i)=etot
-#endif
 #ifdef DEBUG
-        
         write (iout,*) "etrop", i,ii,ib,
      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
      &   entfac(ii),Fdimless(i)
 #endif
+#endif
       enddo   ! i
 #ifdef MPI
       call MPI_Gatherv(Fdimless_(1),scount(me),
@@ -226,6 +230,7 @@ cc        if (wcorr6.eq.0) ecorr6=0.0d0
           entfac(i)=entfac_(i)
         enddo
 #endif
+
 #ifdef DEBUG
         write (iout,*) "The FDIMLESS array before sorting"
         do i=1,ncon
@@ -251,14 +256,12 @@ c        write (iout,*) "qfree",qfree
         sumprob=0.0
         do i=1,min0(ncon,maxstr_proc)-1 
           sumprob=sumprob+dexp(dble(-fdimless(i)+fdimless(1)))/qfree 
-c#define DEBUG
 #ifdef DEBUG
           write (iout,*) 'i=',i,ib,beta_h(ib),
      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
      &     totfree(list_conf(i)),
      &     -entfac(list_conf(i)),fdimless(i),sumprob
 #endif
-c#undef DEBUG
           if (sumprob.gt.prob_limit) goto 122
 c          if (sumprob.gt.1.00d0) goto 122
           nlist=nlist+1

diff --git a/source/cluster/wham/src/readpdb.F b/source/cluster/wham/src/readpdb.F
index 8351702..604b840 100644 (file)
--- a/source/cluster/wham/src/readpdb.F
+++ b/source/cluster/wham/src/readpdb.F
@@ -223,7 +223,7 @@ C Calculate internal coordinates.
       enddo
       endif
 C Calculate internal coordinates.
-      if(me.eq.king.or..not.out1file)then
+      if(print_homology_models.and.(me.eq.king.or..not.out1file))then
        write (iout,'(a)') 
      &   "Backbone and SC coordinates as read from the PDB"
        do ires=1,nres

diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F
index be830b5..037eb1f 100644 (file)
--- a/source/cluster/wham/src/readrtns.F
+++ b/source/cluster/wham/src/readrtns.F
@@ -37,7 +37,8 @@ C
       min_var=(index(controlcard,'MINVAR').gt.0)
       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
-      call readi(controlcard,'NCUT',ncut,1)
+      call readi(controlcard,'NCUT',ncut,0)
+      call readi(controlcard,'NCLUST',nclust,5)
       call readi(controlcard,'NSTART',nstart,0)
       call readi(controlcard,'NEND',nend,0)
       call reada(controlcard,'ECUT',ecut,10.0d0)
@@ -46,7 +47,8 @@ C
       lgrp=(index(controlcard,'LGRP').gt.0)
       caonly=(index(controlcard,'CA_ONLY').gt.0)
       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
-      call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
+      if (ncut.gt.0) 
+     & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
       call readi(controlcard,'IOPT',iopt,2) 
       lside = index(controlcard,"SIDE").gt.0
       efree = index(controlcard,"EFREE").gt.0
@@ -69,6 +71,11 @@ C
       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
       write (iout,*) "with_homology_constr ",with_dihed_constr,
      & " CONSTR_HOMOLOGY",constr_homology
+      print_homology_restraints=
+     & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
+      print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
+      print_homology_models=
+     & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
 
 #ifdef AIX
       call flush_(iout)
@@ -268,6 +275,19 @@ C Convert sequence to numeric code
       do i=1,nres
         itype(i)=rescode(i,sequence(i),iscode)
       enddo
+      if (itype(2).eq.10) then
+        write (iout,*)
+     &   "Glycine is the first full residue, initial dummy deleted"
+        do i=1,nres
+          itype(i)=itype(i+1)
+        enddo
+        nres=nres-1
+      endif
+      if (itype(nres).eq.10) then
+        write (iout,*)
+     &   "Glycine is the last full residue, terminal dummy deleted"
+        nres=nres-1
+      endif
       print *,nres
       print '(20i4)',(itype(i),i=1,nres)
 
@@ -362,7 +382,7 @@ C     both routines should be in dfa.F file!!
       endif
 
       if (constr_homology.gt.0) then
-        call read_constr_homology
+        call read_constr_homology(print_homology_restraints)
       endif
 
 c      if (pdbref) then
@@ -432,7 +452,7 @@ c      endif
             enddo
           enddo
         endif
-        call contact(.true.,ncont_ref,icont_ref)
+        call contact(print_contact_map,ncont_ref,icont_ref)
       endif
 c Read distance restraints
       if (constr_dist.gt.0) then
@@ -885,7 +905,7 @@ c        call flush(iout)
       end
 
 c====-------------------------------------------------------------------
-      subroutine read_constr_homology
+      subroutine read_constr_homology(lprn)
 
       include 'DIMENSIONS'
 #ifdef MPI
@@ -913,7 +933,7 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
       integer ki, i, j, k, l
-      logical lprn /.true./
+      logical lprn
       logical unres_pdb
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
@@ -991,7 +1011,7 @@ c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
       do k=1,constr_homology
 
         read(inp,'(a)') pdbfile
-        write (iout,*) "k ",k," pdbfile ",pdbfile
+c        write (iout,*) "k ",k," pdbfile ",pdbfile
 c  Next stament causes error upon compilation (?)
 c       if(me.eq.king.or. .not. out1file)
 c         write (iout,'(2a)') 'PDB data will be read from file ',
@@ -1026,7 +1046,6 @@ c       tpl_k_sigma_d="template"//kic2//".sigma_d"
           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
      &      (crefjlee(j,i+nres),j=1,3)
         enddo
-#endif
         write (iout,*) "READ HOMOLOGY INFO"
         write (iout,*) "read_constr_homology x: after reading pdb file"
         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
@@ -1035,6 +1054,7 @@ c       tpl_k_sigma_d="template"//kic2//".sigma_d"
         write (iout,*) "waga_theta",waga_theta
         write (iout,*) "waga_d",waga_d
         write (iout,*) "dist_cut",dist_cut
+#endif
 #ifdef AIX
       call flush_(iout)
 #else

diff --git a/source/cluster/wham/src/wrtclust.f b/source/cluster/wham/src/wrtclust.f
index 6012b05..ec43c77 100644 (file)
--- a/source/cluster/wham/src/wrtclust.f
+++ b/source/cluster/wham/src/wrtclust.f
@@ -84,12 +84,16 @@ C 12/8/93 Estimation of "diameters" of the subsequent families.
       ave_dim=0.0
       amax_dim=0.0
 c      write (iout,*) "ecut",ecut
+      emin=totfree(nconf(igr,1))
+c      write (2,*) "emin",emin," ecut",ecut
       do i=2,licz(igr)
         ii=nconf(igr,i)
+c        write (2,*) " igr",igr," i",i," ii",ii," totfree",totfree(ii),
+c     &     " emin",emin," diff",totfree(ii)-emin," ecut",ecut
         if (totfree(ii)-emin .gt. ecut) goto 10
         do j=1,i-1
           jj=nconf(igr,j)
-          if (jj.eq.1) exit
+c          if (jj.eq.1) exit
           if (ii.lt.jj) then
             ind=ioffset(ncon,ii,jj)
           else
@@ -112,9 +116,12 @@ c     &      list_conf(jj),curr_dist
      & '; average distance in the family:',ave_dim 
       rmsave(igr)=0.0d0
       qpart=0.0d0
+      emin=totfree(nconf(igr,1))
       do i=1,licz(igr)
         icon=nconf(igr,i)
-        boltz=dexp(-totfree(icon))
+        boltz=dexp(-totfree(icon)+emin)
+c        write (2,*) "igr",igr," i",i," icon",icon," totfree",
+c     &     totfree(icon)," emin",emin," boltz",boltz," rms",rmstb(icon)
         rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
         qpart=qpart+boltz
       enddo

diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F
index 8a1f9c7..81abea5 100644 (file)
--- a/source/wham/src/energy_p_new.F
+++ b/source/wham/src/energy_p_new.F
@@ -5261,11 +5261,9 @@ c   3 = SC...Ca...Ca...SCi
           cosphi=dcos(j*tauangle(intertyp,i))
           sinphi=dsin(j*tauangle(intertyp,i))
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
-#define DEBUG
 #ifdef DEBUG
           esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
 #endif
-#undef DEBUG
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci

diff --git a/source/wham/src/molread_zs.F b/source/wham/src/molread_zs.F
index a98f47c..9c999cc 100644 (file)
--- a/source/wham/src/molread_zs.F
+++ b/source/wham/src/molread_zs.F
@@ -57,6 +57,19 @@ C Convert sequence to numeric code
       do i=1,nres
         itype(i)=rescode(i,sequence(i),iscode)
       enddo
+      if (itype(2).eq.10) then
+        write (iout,*)
+     &   "Glycine is the first full residue, initial dummy deleted"
+        do i=1,nres
+          itype(i)=itype(i+1)
+        enddo
+        nres=nres-1
+      endif
+      if (itype(nres).eq.10) then
+        write (iout,*)
+     &   "Glycine is the last full residue, terminal dummy deleted"
+        nres=nres-1
+      endif
       write (iout,*) "Numeric code:"
       write (iout,'(20i4)') (itype(i),i=1,nres)
       do i=1,nres-1