assign atoms to default group 0 when groups are not specified

Author: jayponder

source/cluster.f

@@ -24,7 +24,6 @@ subroutine cluster
       use inform
       use iounit
       use keys
-      use limits
       use molcul
       implicit none
       integer i,j,k
@@ -154,37 +153,67 @@ subroutine cluster
 c
 c     pack atoms of each group into a contiguous indexed list
 c
-      do i = 1, n
-         list(i) = grplist(i)
-      end do
-      call sort3 (n,list,kgrp)
+      if (use_group) then
+         do i = 1, n
+            list(i) = grplist(i)
+         end do
+         call sort3 (n,list,kgrp)
 c
 c     find the first and last atom in each of the groups
 c
-      k = list(1)
-      igrp(1,k) = 1
-      do i = 1, n
-         j = list(i)
-         if (j .ne. k) then
-            igrp(2,k) = i - 1
-            igrp(1,j) = i
-            k = j
-         end if
-         ngrp = max(j,ngrp)
-      end do
-      igrp(2,j) = n
+         k = list(1)
+         igrp(1,k) = 1
+         do i = 1, n
+            j = list(i)
+            if (j .ne. k) then
+               igrp(2,k) = i - 1
+               igrp(1,j) = i
+               k = j
+            end if
+            ngrp = max(j,ngrp)
+         end do
+         igrp(2,j) = n
 c
 c     sort the list of atoms in each group by atom number
 c
-      do i = 0, ngrp
-         size = igrp(2,i) - igrp(1,i) + 1
-         if (igrp(1,i) .ne. 0)  call sort (size,kgrp(igrp(1,i)))
-      end do
+         do i = 0, ngrp
+            size = igrp(2,i) - igrp(1,i) + 1
+            if (igrp(1,i) .ne. 0)  call sort (size,kgrp(igrp(1,i)))
+         end do
+c
+c     if groups are not used, put full system in default group
+c
+      else
+         ngrp = 0
+         igrp(1,0) = 1
+         igrp(2,0) = n
+         do i = 1, n
+            kgrp(i) = i
+            grplist(i) = 0
+         end do
+      end if
 c
 c     perform deallocation of some local arrays
 c
       deallocate (list)
 c
+c     compute the total mass of all atoms in each group
+c
+      do i = 0, ngrp
+         grpmass(i) = 0.0d0
+         do j = igrp(1,i), igrp(2,i)
+            grpmass(i) = grpmass(i) + mass(kgrp(j))
+         end do
+      end do
+c
+c     turn off bounds and replicas for intragroup calculations
+c
+      if (use_intra) then
+         use_bounds = .false.
+         use_replica = .false.
+         call cutoffs
+      end if
+c
 c     use only intragroup or intergroup interactions if selected
 c
       if (use_intra) then
@@ -204,7 +233,7 @@ subroutine cluster
          end do
       end if
 c
-c     disable consideration of interactions with any empty groups
+c     disable consideration of interactions with empty groups
 c
       do i = 0, ngrp
          size = igrp(2,i) - igrp(1,i) + 1
@@ -216,27 +245,10 @@ subroutine cluster
          end if
       end do
 c
-c     turn off bounds and replicas for intragroup calculations
-c
-      if (use_intra) then
-         use_bounds = .false.
-         use_replica = .false.
-         call cutoffs
-      end if
-c
-c     compute the total mass of all atoms in each group
-c
-      do i = 0, ngrp
-         grpmass(i) = 0.0d0
-         do j = igrp(1,i), igrp(2,i)
-            grpmass(i) = grpmass(i) + mass(kgrp(j))
-         end do
-      end do
-c
 c     output the final list of atoms in each group
 c
       if (use_group .and. debug) then
-         do i = 0, ngrp
+         do i = 1, ngrp
             size = igrp(2,i) - igrp(1,i) + 1
             if (size .ne. 0) then
                write (iout,50)  i

source/diffuse.f

@@ -37,6 +37,7 @@ program diffuse
       integer nframe,iframe
       integer iarc,start,stop
       integer step,skip,size
+      integer nave,nstart,nstop
       integer, allocatable :: list(:)
       integer, allocatable :: ntime(:)
       real*8 xmid,ymid,zmid
@@ -45,7 +46,8 @@ program diffuse
       real*8 xr,yr,zr,weigh
       real*8 tstep,dunits,delta
       real*8 xvalue,yvalue,zvalue
-      real*8 rvalue,dvalue,counts
+      real*8 rvalue,counts
+      real*8 dvalue,dave
       real*8, allocatable :: xmsd(:)
       real*8, allocatable :: ymsd(:)
       real*8, allocatable :: zmsd(:)
@@ -332,16 +334,25 @@ program diffuse
      &           7x,'R MSD',4x,'Diff Const',
      &        /,7x,'(ps)',9x,'(/2)',8x,'(/2)',8x,'(/2)',8x,'(/6)',
      &           5x,'(x 10^5)',/)
+      nstart = nint(dble(nframe-1)/3.0d0)
+      nstop = 2 * nstart
+      nave = nstop - nstart + 1
+      dave = 0.0d0
       do i = 1, nframe-1
          delta = tstep * dble(i)
          xvalue = xmsd(i) / 2.0d0
          yvalue = ymsd(i) / 2.0d0
          zvalue = zmsd(i) / 2.0d0
          rvalue = (xmsd(i) + ymsd(i) + zmsd(i)) / 6.0d0
          dvalue = rvalue / delta
+         if (i.ge.nstart .and. i.le.nstop)  dave = dave + dvalue/nave
          write (iout,190)  delta,xvalue,yvalue,zvalue,rvalue,dvalue
   190    format (f12.2,4f12.2,f12.4)
       end do
+      if (nave .ge. 8) then
+         write (iout,200)  dave
+  200    format (/,' Self-Diffusion Constant :',f12.4,' x 10^-5 cm^2/s')     
+      end if
 c
 c     perform deallocation of some local arrays
 c