cluster & wham Adam's PBC corrections

diff --git a/source/cluster/wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos b/source/cluster/wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos
index d1e64a7..4f7f61f 100644 (file)
--- a/source/cluster/wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos
+++ b/source/cluster/wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos
@@ -1,7 +1,7 @@
 #INSTALL_DIR = /opt/cray/mpt/7.3.2/gni/mpich-intel/15.0
 FC = ftn
-#OPT = -O3 -ip -mcmodel=medium -shared-intel -dynamic
-OPT = -CB -g -mcmodel=medium -shared-intel -dynamic
+OPT = -O3 -ip -mcmodel=medium -shared-intel -dynamic
+#OPT = -CB -g -mcmodel=medium -shared-intel -dynamic
 FFLAGS =  ${OPT} -c -I. -Iinclude_unres -I$(INSTALL_DIR)/include
 LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a
 

diff --git a/source/cluster/wham/src-HCD-5D/geomout.F b/source/cluster/wham/src-HCD-5D/geomout.F
index 4ef656f..be43686 100644 (file)
--- a/source/cluster/wham/src-HCD-5D/geomout.F
+++ b/source/cluster/wham/src-HCD-5D/geomout.F
@@ -17,12 +17,15 @@
       iatom=0
       ichain=1
       ires=0
+      iti_prev=0
       do i=nnt,nct
         iti=itype(i)
         if (iti.eq.ntyp1) then
-          ichain=ichain+1
           ires=0
-          write (ipdb,'(a)') 'TER'
+          if (iti_prev.ne.ntyp1) then
+            write (ipdb,'(a)') 'TER'
+            ichain=ichain+1
+          endif
         else
         ires=ires+1
         iatom=iatom+1
@@ -35,6 +38,7 @@
      &      ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
         endif
         endif
+        iti_prev=iti
       enddo
       write (ipdb,'(a)') 'TER'
       do i=nnt,nct-1

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 cac0d06..b188deb 100644 (file)
--- a/source/cluster/wham/src-HCD-5D/read_constr_homology.F
+++ b/source/cluster/wham/src-HCD-5D/read_constr_homology.F
@@ -144,8 +144,10 @@ c
         unres_pdb=.false.
         if (read2sigma) then
           call readpdb_template(k)
+          close(ipdbin)
         else
           call readpdb(out_template_coord)
+          close(ipdbin)
         endif
 
 c        call readpdb

diff --git a/source/cluster/wham/src-HCD-5D/sizesclu.dat b/source/cluster/wham/src-HCD-5D/sizesclu.dat
index 7d0d666..cb8572c 100644 (file)
--- a/source/cluster/wham/src-HCD-5D/sizesclu.dat
+++ b/source/cluster/wham/src-HCD-5D/sizesclu.dat
@@ -5,7 +5,7 @@
 * Max. number of conformations in the data set.
 *
       integer maxconf,maxstr_proc
-      PARAMETER (MAXCONF=8000)
+      PARAMETER (MAXCONF=12000)
       parameter (maxstr_proc=maxconf/2)
 *
 * Max. number of "distances" between conformations.

diff --git a/source/unres/src-HCD-5D/readpdb.F b/source/unres/src-HCD-5D/readpdb.F
index 1190bb1..c56b1df 100644 (file)
--- a/source/unres/src-HCD-5D/readpdb.F
+++ b/source/unres/src-HCD-5D/readpdb.F
@@ -30,6 +30,7 @@ C geometry.
       lsecondary=.false.
       nhfrag=0
       nbfrag=0
+      ires=0
       do
         read (ipdbin,'(a80)',end=10) card
         if (card(:5).eq.'HELIX') then

diff --git a/source/wham/src-HCD-5D/COMMON.CHAIN b/source/wham/src-HCD-5D/COMMON.CHAIN
index 7369baa..dfffc78 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,ichanres,tabpermchain,
-     & nchain ,npermchain,ireschain,iz_sc
+     & ishift_pdb,chain_length,chain_border,chain_border1,ichanres,
+     & tabpermchain,nchain ,npermchain,ireschain,iz_sc
       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),
@@ -11,7 +11,8 @@
      & rmssing,anatemp,iz_sc,nsup,
      & nstart_sup,nend_sup,chain_length(maxchain),npermchain,
      & ireschain(maxres),tabpermchain(maxchain,maxperm),
-     & chain_border(2,maxchain),nchain,nstart_seq,ishift_pdb
+     & chain_border(2,maxchain),chain_border1(2,maxchain),nchain,
+     & nstart_seq,ishift_pdb
       double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,
      &  sssgrad,
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick

diff --git a/source/wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos b/source/wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos
index 5280f0a..e667382 100644 (file)
--- a/source/wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos
+++ b/source/wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos
@@ -1,9 +1,9 @@
 BIN = ~/bin
 FC = ftn
-#OPT = -mcmodel=medium -shared-intel -O3 -dynamic
+OPT = -mcmodel=medium -shared-intel -O3 -dynamic
 #OPT = -O3 -intel-static -mcmodel=medium 
 #OPT = -O3 -ip -w 
-OPT = -g -CB -mcmodel=medium -shared-intel -dynamic
+#OPT = -g -CB -mcmodel=medium -shared-intel -dynamic
 FFLAGS = ${OPT} -c -I. -I./include_unres -I$(INSTALL_DIR)/include
 LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a
 

diff --git a/source/wham/src-HCD-5D/enecalc1.F b/source/wham/src-HCD-5D/enecalc1.F
index f037ae8..0040e37 100644 (file)
--- a/source/wham/src-HCD-5D/enecalc1.F
+++ b/source/wham/src-HCD-5D/enecalc1.F
@@ -213,8 +213,8 @@ c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
      &                            ((c(l,k+nres),l=1,3),k=nnt,nct)
 c              call intout
-              call pdbout(indstart(me1)+iii,
-     & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0)
+c              call pdbout(indstart(me1)+iii,
+c     & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0)
               call enerprint(energia(0),fT)
               errmsg_count=errmsg_count+1
               if (errmsg_count.gt.maxerrmsg_count) 

diff --git a/source/wham/src-HCD-5D/molread_zs.F b/source/wham/src-HCD-5D/molread_zs.F
index d7f586d..878e4dd 100644 (file)
--- a/source/wham/src-HCD-5D/molread_zs.F
+++ b/source/wham/src-HCD-5D/molread_zs.F
@@ -73,10 +73,18 @@ C Convert sequence to numeric code
       nct=nres
       call seq2chains(nres,itype,nchain,chain_length,chain_border,
      &  ireschain)
+      chain_border1(1,1)=1
+      chain_border1(2,1)=chain_border(2,1)+1
+      do i=2,nchain-1
+        chain_border1(1,i)=chain_border(1,i)-1
+        chain_border1(2,i)=chain_border(2,i)+1
+      enddo
+      chain_border1(1,nchain)=chain_border(1,nchain)-1
+      chain_border1(2,nchain)=nres
       write(iout,*) "nres",nres," nchain",nchain
       do i=1,nchain
         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
-     &    chain_border(2,i)
+     &    chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
       enddo
       call chain_symmetry(nchain,nres,itype,chain_border,
      &    chain_length,npermchain,tabpermchain)

diff --git a/source/wham/src-HCD-5D/oligomer.F b/source/wham/src-HCD-5D/oligomer.F
index 34b7be0..c63ea47 100644 (file)
--- a/source/wham/src-HCD-5D/oligomer.F
+++ b/source/wham/src-HCD-5D/oligomer.F
@@ -4,73 +4,142 @@
       include "COMMON.CHAIN"
       include "COMMON.INTERACT"
       include "COMMON.IOUNITS"
-      integer i,ii,ipi,ipj,ipmin,j,jmin,k,ix,iy,iz,
-     &  ixmin,iymin,izmin,ir_start,ir_end
-      integer iper(maxchain),iaux
-      double precision dchain,dchainmin,cmchain(3,20)
-      cmchain=0.0d0
-      do i=1,nchain
-        ii=0
-        do j=chain_border(1,i),chain_border(2,i)
-          if (itype(j).eq.ntyp1) cycle
-          ii=ii+1
-          do k=1,3
-            cmchain(k,i)=cmchain(k,i)+c(k,j)
+      integer i,j,k,l,ichain,jchain,ii,lshift(maxchain),lmin(maxchain),
+     & nrestot
+      double precision sumodl,sumodl_min,shift_coord(3,maxchain),
+     & boxsize(3),cm(3,maxchain),ccm(3),ccm_prev(3)
+      logical previous,change
+      boxsize(1)=boxxsize
+      boxsize(2)=boxysize
+      boxsize(3)=boxzsize
+      nrestot=0
+      do ichain=1,nchain
+c        print *,"ichain",ichain
+        nrestot=nrestot+chain_length(ichain)
+        do j=1,3
+          cm(j,ichain)=0.0d0
+        enddo
+        do i=chain_border(1,ichain),chain_border(2,ichain)
+c          print *,i,c(:,i)
+          do j=1,3
+            cm(j,ichain)=cm(j,ichain)+c(j,i)
           enddo
         enddo
-        do k=1,3
-          cmchain(k,i)=cmchain(k,i)/ii
+        do j=1,3
+          cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
         enddo
       enddo
+#ifdef DEBUG
+      write (iout,*) "CM before shift"
       do i=1,nchain
-        iper(i)=i
+        write (iout,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
       enddo
-      do i=1,nchain-1
-        dchainmin=1.0d10
-        do j=i+1,nchain
-          ipi=iper(i)
-          ipj=iper(j)
-          do ix=-1,1
-            do iy=-1,1
-              do iz=-1,1
-                dchain=(cmchain(1,ipj)-cmchain(1,ipi)+ix*boxxsize)**2+
-     &                 (cmchain(2,ipj)-cmchain(2,ipi)+iy*boxysize)**2+
-     &                 (cmchain(3,ipj)-cmchain(3,ipi)+iz*boxzsize)**2
-c                write (iout,*) "i",i," ipi",ipi," j",j," ipj",ipj," d",
-c     &               dsqrt(dchain)," dmin",dsqrt(dchainmin)," jmin",jmin
-                if (dchain.lt.dchainmin) then
-                  dchainmin=dchain
-                  ixmin=ix
-                  iymin=iy
-                  izmin=iz
-                  jmin=j
-                endif
-              enddo
+#endif DEBUG
+      do j=1,3
+        lmin=0
+        sumodl_min=1.0d10
+        do i=0,3**nchain-1
+          ii=i
+          do ichain=1,nchain
+            lshift(ichain)=mod(ii,3)-1
+            ii=ii/3
+          enddo
+#ifdef DEBUG
+          write (iout,'(10i3)') (lshift(ichain),ichain=1,nchain)
+#endif
+          sumodl=0.0d0
+          do ichain=1,nchain
+            do jchain=1,ichain-1
+            sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
+     &          -cm(j,jchain)-lshift(jchain)*boxsize(j))
             enddo
           enddo
+#ifdef DEBUG
+          write(iout,*) "sumodl",sumodl," sumodl_min",sumodl_min
+#endif
+          if (sumodl.lt.sumodl_min) then
+            sumodl_min=sumodl
+            lmin=lshift
+          endif
         enddo
-        if (ixmin.eq.0 .and. iymin.eq.0 .and. izmin.eq.0) cycle
-        ipj=iper(jmin)
-        cmchain(1,ipj)=cmchain(1,ipj)+ixmin*boxxsize
-        cmchain(2,ipj)=cmchain(2,ipj)+iymin*boxysize
-        cmchain(3,ipj)=cmchain(3,ipj)+izmin*boxzsize
-        ir_start=chain_border(1,ipj)
-        if (ir_start.gt.1) ir_start=ir_start-1 
-        ir_end=chain_border(2,ipj)
-        if (ir_end.lt.nres) ir_end=ir_end+1
-        do k=ir_start,ir_end
-          c(1,k)=c(1,k)+ixmin*boxxsize
-          c(2,k)=c(2,k)+iymin*boxysize
-          c(3,k)=c(3,k)+izmin*boxzsize
-          c(1,k+nres)=c(1,k+nres)+ixmin*boxxsize
-          c(2,k+nres)=c(2,k+nres)+iymin*boxysize
-          c(3,k+nres)=c(3,k+nres)+izmin*boxzsize
+        do ichain=1,nchain
+          cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
+          shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
+        enddo
+      enddo
+      do ichain=1,nchain
+        do j=1,3
+          ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
         enddo
-c        write (iout,*) "jmin",jmin," ipj",ipj,
-c     &   " ixmin",ixmin," iymin",iymin," izmin",izmin
-        iaux=iper(i+1)
-        iper(i+1)=iper(jmin)
-        iper(jmin)=iaux
       enddo
+      do j=1,3
+        ccm(j)=ccm(j)/nrestot
+      enddo
+#ifdef DEBUG
+      write (iout,'(a,3f10.5)') "ccm     ",(ccm(j),j=1,3)
+      write (iout,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
+      write (iout,'(a,3f10.5)') "diff    ",(ccm(j)-ccm_prev(j),j=1,3)
+      write (iout,*) "shift_coord"
+      do i=1,nchain
+        write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+      enddo
+#endif;
+#ifdef PREVIOUS
+c
+c Shift the system by box dimensions so that its CM be closest to the
+c CM of the pevious snapshot
+c 
+      if (previous) then
+        do j=1,3
+          change=.true.
+          do while (change)
+          if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
+            do ichain=1,nchain
+              shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
+            enddo
+            ccm(j)=ccm(j)-boxsize(j)
+          else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
+            do ichain=1,nchain
+              shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
+            enddo
+            ccm(j)=ccm(j)+boxsize(j)
+          else
+            change=.false.
+          endif
+        enddo
+        enddo
+        do j=1,3
+          ccm_prev(j)=ccm(j)
+        enddo
+      else
+        previous=.true.
+        do j=1,3
+          ccm_prev(j)=ccm(j)
+        enddo
+      endif
+#else
+c Simply shift the CM of the whole oligomer to the origin of the
+c coordinate system 
+      do ichain=1,nchain
+        do j=1,3
+          shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
+        enddo
+      enddo 
+#endif
+#ifdef DEBUG
+      write (iout,*) "shift_coord"
+      do i=1,nchain
+        write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+      enddo
+#endif
+      do ichain=1,nchain
+        do i=chain_border1(1,ichain),chain_border1(2,ichain)
+          do j=1,3
+            c(j,i)=c(j,i)+shift_coord(j,ichain)
+            c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
+          enddo
+        enddo
+      enddo 
+
       return
       end

diff --git a/source/wham/src-HCD-5D/read_constr_homology.F b/source/wham/src-HCD-5D/read_constr_homology.F
index ab9901d..dece50f 100644 (file)
--- a/source/wham/src-HCD-5D/read_constr_homology.F
+++ b/source/wham/src-HCD-5D/read_constr_homology.F
@@ -145,8 +145,10 @@ c
         unres_pdb=.false.
         if (read2sigma) then
           call readpdb_template(k)
+          close(ipdbin)
         else
           call readpdb
+          close(ipdbin)
         endif
 
 c        call readpdb
@@ -163,7 +165,6 @@ c        call readpdb
         write (iout,*) "read_constr_homology: after reading pdb file"
         call flush(iout)
 #endif
-
 c
 c     Distance restraints
 c

diff --git a/source/wham/src-HCD-5D/readpdb.F b/source/wham/src-HCD-5D/readpdb.F
index 8a85b8c..83c63ab 100644 (file)
--- a/source/wham/src-HCD-5D/readpdb.F
+++ b/source/wham/src-HCD-5D/readpdb.F
@@ -27,9 +27,10 @@ C geometry.
       ibeg=1
       ishift1=0
       sccalc=.false.
+      ires=0
       do
         read (ipdbin,'(a80)',end=10) card
-!       write (iout,'(a)') card
+c        write (iout,'(a)') card
         if (card(:5).eq.'HELIX') then
           nhfrag=nhfrag+1
           lsecondary=.true.
@@ -58,7 +59,7 @@ C geometry.
           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)
@@ -119,7 +120,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