cluster_wham src-M agrees with src for single chain

diff --git a/source/cluster/wham/src-M/COMMON.CLUSTER b/source/cluster/wham/src-M/COMMON.CLUSTER
index 4477d19..f1ad0fd 100644 (file)
--- a/source/cluster/wham/src-M/COMMON.CLUSTER
+++ b/source/cluster/wham/src-M/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-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index c7c2c2f..4551e68 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -3349,7 +3349,8 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.21) cycle
+        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &(itype(i).eq.ntyp1)) cycle
         dethetai=0.0d0
         dephii=0.0d0
         dephii1=0.0d0
@@ -3359,7 +3360,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.21) then
+        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -3373,13 +3374,13 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+          ityp1=ithetyp(itype(i-2))
           do k=1,nsingle
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.21) then
+        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -3394,7 +3395,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0

diff --git a/source/cluster/wham/src-M/main_clust.F b/source/cluster/wham/src-M/main_clust.F
index a2e4769..d33ce25 100644 (file)
--- a/source/cluster/wham/src-M/main_clust.F
+++ b/source/cluster/wham/src-M/main_clust.F
@@ -23,7 +23,7 @@ C
       logical printang(max_cut)
       integer printpdb(max_cut)
       integer printmol2(max_cut)
-      character*240 lineh
+      character*240 lineh,scrachdir2d
       REAL CRIT(maxconf),MEMBR(maxconf)
       REAL CRITVAL(maxconf-1)
       INTEGER IA(maxconf),IB(maxconf)
@@ -34,13 +34,14 @@ 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
+     & it,ncon_work,ind1,ilen,is,ie
       double precision t1,t2,tcpu,difconf
       real diss_(maxdist)
       
       double precision varia(maxvar)
       double precision hrtime,mintime,sectime
       logical eof
+      external ilen
 #ifdef MPI
       call MPI_Init( IERROR )
       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
@@ -80,6 +81,12 @@ c      write (iout,*) "after permut"
 c      call flush(iout)
       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
@@ -91,12 +98,21 @@ c      call flush(iout)
           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
@@ -112,7 +128,6 @@ c      call flush(iout)
 #ifdef MPI
       call work_partition(.true.,ncon)
 #endif
-
       call probabl(iT,ncon_work,ncon,*20)
 
       if (ncon_work.lt.2) then
@@ -174,7 +189,8 @@ c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
       if (me.eq.master) then
 #endif
-      open(80,file='/tmp/distance',form='unformatted')
+      scrachdir2d=scratchdir(:ilen(scratchdir))//'distance'
+      open(80,file=scrachdir2d,form='unformatted')
       do i=1,ndis
         write(80) diss(i)
       enddo
@@ -247,29 +263,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 +308,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),
@@ -310,17 +337,17 @@ C
 C
       close(icbase,status="delete")
 #ifdef MPI
-      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
 #endif
       stop '********** Program terminated normally.'
    20 write (iout,*) "Error reading coordinates"
 #ifdef MPI
-      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
 #endif
       stop
    30 write (iout,*) "Error reading reference structure"
 #ifdef MPI
-      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
 #endif
       stop
       end

diff --git a/source/cluster/wham/src-M/read_coords.F b/source/cluster/wham/src-M/read_coords.F
index c34aca4..f1e2c7b 100644 (file)
--- a/source/cluster/wham/src-M/read_coords.F
+++ b/source/cluster/wham/src-M/read_coords.F
@@ -212,6 +212,7 @@ c          call flush(iout)
                enddo
              enddo
           else
+            itmp=0
 #if (defined(AIX) && !defined(JUBL))
             call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
             if (iret.eq.0) goto 101

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index 4c15495..e33b6d6 100644 (file)
--- a/source/cluster/wham/src-M/readrtns.F
+++ b/source/cluster/wham/src-M/readrtns.F
@@ -38,7 +38,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,'SYM',symetr,1)
       write (iout,*) 'sym', symetr
       call readi(controlcard,'NSTART',nstart,0)
@@ -49,7 +50,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
@@ -206,6 +208,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)
 

diff --git a/source/cluster/wham/src-M/wrtclust.f b/source/cluster/wham/src-M/wrtclust.f
index f2f3eb7..3915ebc 100644 (file)
--- a/source/cluster/wham/src-M/wrtclust.f
+++ b/source/cluster/wham/src-M/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
@@ -191,7 +198,7 @@ c Write conformations of the family i to PDB files
 c            write (iout,*) i,ncon_out,nconf(i,ncon_out),
 c     &        totfree(nconf(i,ncon_out)),emin1,ecut
           enddo
-          write (iout,*) "ncon_out",ncon_out
+c          write (iout,*) "ncon_out",ncon_out
           call flush(iout)
           do j=1,nres
             tempfac(1,j)=5.0d0

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index af921d0..709a471 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -3399,7 +3399,8 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.21) cycle
+        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &(itype(i).eq.ntyp1)) cycle
         dethetai=0.0d0
         dephii=0.0d0
         dephii1=0.0d0
@@ -3409,7 +3410,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.21) then
+        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -3423,13 +3424,13 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+          ityp1=ithetyp(itype(i-2))
           do k=1,nsingle
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.21) then
+        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -3444,7 +3445,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0