source/wham/src-M/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 nchain,i,ii,ipi,ipj,ipmin,j,jmin,k,ix,iy,iz,
   8      &  ixmin,iymin,izmin
   9       logical newchain
  10       integer ichain(2,20),iper(20),iaux
  11       double precision dchain,dchainmin,cmchain(3,20)
  12       nchain=1
  13       newchain=.false.
  14       ichain(1,nchain)=1
  15       do i=2,nres
  16         if (itype(i).eq.ntyp1) then
  17           if (.not.newchain) then
  18             ichain(2,nchain)=i
  19             if (i.lt.nres) nchain=nchain+1
  20             newchain=.true.
  21           else
  22             newchain=.false.
  23             ichain(1,nchain)=i
  24           endif
  25         endif
  26       enddo
  27 #ifdef DEBUG
  28       write (iout,*) "Chains"
  29       do i=1,nchain
  30         write (iout,*) i,ichain(1,i),ichain(2,i)
  31       enddo
  32 #endif
  33       cmchain=0.0d0
  34       do i=1,nchain
  35         ii=0
  36         do j=ichain(1,i),ichain(2,i)
  37           if (itype(j).eq.ntyp1) cycle
  38           ii=ii+1
  39           do k=1,3
  40             cmchain(k,i)=cmchain(k,i)+c(k,j)
  41           enddo
  42         enddo
  43         do k=1,3
  44           cmchain(k,i)=cmchain(k,i)/ii
  45         enddo
  46       enddo
  47       do i=1,nchain
  48         iper(i)=i
  49       enddo
  50       do i=1,nchain-1
  51         dchainmin=1.0d10
  52         do j=i+1,nchain
  53           ipi=iper(i)
  54           ipj=iper(j)
  55           do ix=-1,1
  56             do iy=-1,1
  57               do iz=-1,1
  58                 dchain=(cmchain(1,ipj)-cmchain(1,ipi)+ix*boxxsize)**2+
  59      &                 (cmchain(2,ipj)-cmchain(2,ipi)+iy*boxysize)**2+
  60      &                 (cmchain(3,ipj)-cmchain(3,ipi)+iz*boxzsize)**2
  61 c                write (iout,*) "i",i," ipi",ipi," j",j," ipj",ipj," d",
  62 c     &               dsqrt(dchain)," dmin",dsqrt(dchainmin)," jmin",jmin
  63                 if (dchain.lt.dchainmin) then
  64                   dchainmin=dchain
  65                   ixmin=ix
  66                   iymin=iy
  67                   izmin=iz
  68                   jmin=j
  69                 endif
  70               enddo
  71             enddo
  72           enddo
  73         enddo
  74         if (ixmin.eq.0 .and. iymin.eq.0 .and. izmin.eq.0) cycle
  75         ipj=iper(jmin)
  76         cmchain(1,ipj)=cmchain(1,ipj)+ixmin*boxxsize
  77         cmchain(2,ipj)=cmchain(2,ipj)+iymin*boxysize
  78         cmchain(3,ipj)=cmchain(3,ipj)+izmin*boxzsize
  79         do k=ichain(1,ipj),ichain(2,ipj)
  80           c(1,k)=c(1,k)+ixmin*boxxsize
  81           c(2,k)=c(2,k)+iymin*boxysize
  82           c(3,k)=c(3,k)+izmin*boxzsize
  83           c(1,k+nres)=c(1,k+nres)+ixmin*boxxsize
  84           c(2,k+nres)=c(2,k+nres)+iymin*boxysize
  85           c(3,k+nres)=c(3,k+nres)+izmin*boxzsize
  86         enddo
  87 c        write (iout,*) "jmin",jmin," ipj",ipj,
  88 c     &   " ixmin",ixmin," iymin",iymin," izmin",izmin
  89         iaux=iper(i+1)
  90         iper(i+1)=iper(jmin)
  91         iper(jmin)=iaux
  92       enddo
  93       return
  94       end