source/wham/src-HCD/oligomer.F

   1       subroutine oligomer
   2       implicit none
   3       include "DIMENSIONS"
   4       include "COMMON.CHAIN"
   5       include "COMMON.INTERACT"
   6       include "COMMON.IOUNITS"
   7       integer i,j,k,l,ichain,jchain,ii,lshift(maxchain),lmin(maxchain),
   8      & nrestot
   9       double precision sumodl,sumodl_min,shift_coord(3,maxchain),
  10      & boxsize(3),cm(3,maxchain),ccm(3),ccm_prev(3)
  11       logical previous,change
  12       boxsize(1)=boxxsize
  13       boxsize(2)=boxysize
  14       boxsize(3)=boxzsize
  15       nrestot=0
  16       do ichain=1,nchain
  17 c        print *,"ichain",ichain
  18         nrestot=nrestot+chain_length(ichain)
  19         do j=1,3
  20           cm(j,ichain)=0.0d0
  21         enddo
  22         do i=chain_border(1,ichain),chain_border(2,ichain)
  23 c          print *,i,c(:,i)
  24           do j=1,3
  25             cm(j,ichain)=cm(j,ichain)+c(j,i)
  26           enddo
  27         enddo
  28         do j=1,3
  29           cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
  30         enddo
  31       enddo
  32 #ifdef DEBUG
  33       write (iout,*) "CM before shift"
  34       do i=1,nchain
  35         write (iout,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
  36       enddo
  37 #endif DEBUG
  38       do j=1,3
  39         lmin=0
  40         sumodl_min=1.0d10
  41         do i=0,3**nchain-1
  42           ii=i
  43           lshift(1)=0
  44           do ichain=2,nchain
  45             lshift(ichain)=mod(ii,3)-1
  46             ii=ii/3
  47           enddo
  48 #ifdef DEBUG
  49           write (iout,'(10i3)') (lshift(ichain),ichain=1,nchain)
  50 #endif
  51           sumodl=0.0d0
  52           do ichain=1,nchain
  53             do jchain=1,ichain-1
  54             sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
  55      &          -cm(j,jchain)-lshift(jchain)*boxsize(j))
  56             enddo
  57           enddo
  58 #ifdef DEBUG
  59           write(iout,*) "sumodl",sumodl," sumodl_min",sumodl_min
  60 #endif
  61           if (sumodl.lt.sumodl_min) then
  62             sumodl_min=sumodl
  63             lmin=lshift
  64           endif
  65         enddo
  66         do ichain=1,nchain
  67           cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
  68           shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
  69         enddo
  70       enddo
  71       do ichain=1,nchain
  72         do j=1,3
  73           ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
  74         enddo
  75       enddo
  76       do j=1,3
  77         ccm(j)=ccm(j)/nrestot
  78       enddo
  79 #ifdef DEBUG
  80       write (iout,'(a,3f10.5)') "ccm     ",(ccm(j),j=1,3)
  81       write (iout,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
  82       write (iout,'(a,3f10.5)') "diff    ",(ccm(j)-ccm_prev(j),j=1,3)
  83       write (iout,*) "shift_coord"
  84       do i=1,nchain
  85         write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
  86       enddo
  87 #endif;
  88 #ifdef PREVIOUS
  89 c
  90 c Shift the system by box dimensions so that its CM be closest to the
  91 c CM of the pevious snapshot
  92 c
  93       if (previous) then
  94         do j=1,3
  95           change=.true.
  96           do while (change)
  97           if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
  98             do ichain=1,nchain
  99               shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
 100             enddo
 101             ccm(j)=ccm(j)-boxsize(j)
 102           else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
 103             do ichain=1,nchain
 104               shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
 105             enddo
 106             ccm(j)=ccm(j)+boxsize(j)
 107           else
 108             change=.false.
 109           endif
 110         enddo
 111         enddo
 112         do j=1,3
 113           ccm_prev(j)=ccm(j)
 114         enddo
 115       else
 116         previous=.true.
 117         do j=1,3
 118           ccm_prev(j)=ccm(j)
 119         enddo
 120       endif
 121 #else
 122 c Simply shift the CM of the whole oligomer to the origin of the
 123 c coordinate system
 124       do ichain=1,nchain
 125         do j=1,3
 126           shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
 127         enddo
 128       enddo
 129 #endif
 130 #ifdef DEBUG
 131       write (iout,*) "shift_coord"
 132       do i=1,nchain
 133         write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
 134       enddo
 135 #endif
 136       do ichain=1,nchain
 137 c        write (iout,*) "shift_coord",shift_coord(:,ichain)
 138         do i=chain_border1(1,ichain),chain_border1(2,ichain)
 139           do j=1,3
 140             c(j,i)=c(j,i)+shift_coord(j,ichain)
 141             c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
 142           enddo
 143         enddo
 144       enddo
 145
 146       return
 147       end