X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-HCD%2Foligomer.F;h=c63ea4704c059ef6f81a685b49e9e156783c9bf8;hb=df61eb1a0ad598d7c6c5b4037dee4bec200d5bc2;hp=34b7be0bbcaabeb3b03739d01928bac6bd123519;hpb=5aac27f7c6317fa0975184363614b42716034716;p=unres.git

diff --git a/source/wham/src-HCD/oligomer.F b/source/wham/src-HCD/oligomer.F
index 34b7be0..c63ea47 100644
--- a/source/wham/src-HCD/oligomer.F
+++ b/source/wham/src-HCD/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