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           do ichain=1,nchain
  44             lshift(ichain)=mod(ii,3)-1
  45             ii=ii/3
  46           enddo
  47 #ifdef DEBUG
  48           write (iout,'(10i3)') (lshift(ichain),ichain=1,nchain)
  49 #endif
  50           sumodl=0.0d0
  51           do ichain=1,nchain
  52             do jchain=1,ichain-1
  53             sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
  54      &          -cm(j,jchain)-lshift(jchain)*boxsize(j))
  55             enddo
  56           enddo
  57 #ifdef DEBUG
  58           write(iout,*) "sumodl",sumodl," sumodl_min",sumodl_min
  59 #endif
  60           if (sumodl.lt.sumodl_min) then
  61             sumodl_min=sumodl
  62             lmin=lshift
  63           endif
  64         enddo
  65         do ichain=1,nchain
  66           cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
  67           shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
  68         enddo
  69       enddo
  70       do ichain=1,nchain
  71         do j=1,3
  72           ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
  73         enddo
  74       enddo
  75       do j=1,3
  76         ccm(j)=ccm(j)/nrestot
  77       enddo
  78 #ifdef DEBUG
  79       write (iout,'(a,3f10.5)') "ccm     ",(ccm(j),j=1,3)
  80       write (iout,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
  81       write (iout,'(a,3f10.5)') "diff    ",(ccm(j)-ccm_prev(j),j=1,3)
  82       write (iout,*) "shift_coord"
  83       do i=1,nchain
  84         write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
  85       enddo
  86 #endif;
  87 #ifdef PREVIOUS
  88 c
  89 c Shift the system by box dimensions so that its CM be closest to the
  90 c CM of the pevious snapshot
  91 c
  92       if (previous) then
  93         do j=1,3
  94           change=.true.
  95           do while (change)
  96           if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
  97             do ichain=1,nchain
  98               shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
  99             enddo
 100             ccm(j)=ccm(j)-boxsize(j)
 101           else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
 102             do ichain=1,nchain
 103               shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
 104             enddo
 105             ccm(j)=ccm(j)+boxsize(j)
 106           else
 107             change=.false.
 108           endif
 109         enddo
 110         enddo
 111         do j=1,3
 112           ccm_prev(j)=ccm(j)
 113         enddo
 114       else
 115         previous=.true.
 116         do j=1,3
 117           ccm_prev(j)=ccm(j)
 118         enddo
 119       endif
 120 #else
 121 c Simply shift the CM of the whole oligomer to the origin of the
 122 c coordinate system
 123       do ichain=1,nchain
 124         do j=1,3
 125           shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
 126         enddo
 127       enddo
 128 #endif
 129 #ifdef DEBUG
 130       write (iout,*) "shift_coord"
 131       do i=1,nchain
 132         write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
 133       enddo
 134 #endif
 135       do ichain=1,nchain
 136         do i=chain_border1(1,ichain),chain_border1(2,ichain)
 137           do j=1,3
 138             c(j,i)=c(j,i)+shift_coord(j,ichain)
 139             c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
 140           enddo
 141         enddo
 142       enddo
 143
 144       return
 145       end