program WHAM_multparm
c Creation/update of the database of conformations
      implicit none
#ifndef ISNAN
      external proc_proc
#endif
#ifdef WINPGI
cMS$ATTRIBUTES C ::  proc_proc
#endif
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "DIMENSIONS.FREE"
#ifdef MPI
      include "mpif.h"
      integer IERROR,ERRCODE
      include "COMMON.MPI"
#endif
      include "COMMON.IOUNITS"
      include "COMMON.FREE"
      include "COMMON.CONTROL"
      include "COMMON.ALLPARM"
      include "COMMON.PROT"
      double precision rr,x(max_paropt)
      integer idumm
      integer i,ipar,islice
#ifdef MPI
      call MPI_Init( IERROR )
      call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
      call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
      Master = 0
      if (ierror.gt.0) then
        write(iout,*) "SEVERE ERROR - Can't initialize MPI."
        call mpi_finalize(ierror)
        stop
      endif
      if (nprocs.gt.MaxProcs+1) then
        write (2,*) "Error - too many processors",
     &   nprocs,MaxProcs+1
        write (2,*) "Increase MaxProcs and recompile"
        call MPI_Finalize(IERROR)
        stop
      endif
#endif
c NaNQ initialization
#ifndef ISNAN
      i=-1
      rr=dacos(100.0d0)
#ifdef WINPGI
      idumm=proc_proc(rr,i)
#else
      call proc_proc(rr,i)
#endif
#endif
      call initialize
      call openunits
      call cinfo
      call read_general_data(*10)
      call flush(iout)
      call molread(*10)
      call flush(iout)
#ifdef MPI 
      write (iout,*) "Calling proc_groups"
      call proc_groups
      write (iout,*) "proc_groups exited"
      call flush(iout)
#endif
      do ipar=1,nParmSet
        write (iout,*) "Calling parmread",ipar
        call parmread(ipar,*10)
        if (.not.separate_parset) then
          call store_parm(ipar)
          write (iout,*) "Finished storing parameters",ipar
        else if (ipar.eq.myparm) then
          call store_parm(1)
          write (iout,*) "Finished storing parameters",ipar
        endif
        call flush(iout)
      enddo
      call read_efree(*10)
      write (iout,*) "Finished READ_EFREE"
      call flush(iout)
      call read_protein_data(*10)
      write (iout,*) "Finished READ_PROTEIN_DATA"
      call flush(iout)
      if (indpdb.gt.0) then
        call promienie
        call read_compar
        call read_ref_structure(*10)
        call proc_cont
        call fragment_list
      endif
      write (iout,*) "Begin read_database"
      call flush(iout)
      call read_database(*10)
      write (iout,*) "Finished read_database"
      call flush(iout)
      if (separate_parset) nparmset=1
      do islice=1,nslice
        if (ntot(islice).gt.0) then
#ifdef MPI 
          call work_partition(islice,.true.)
          write (iout,*) "work_partition OK"
          call flush(iout)
#endif
          call enecalc(islice,*10)
          write (iout,*) "enecalc OK"
          call flush(iout)
          call WHAM_CALC(islice,*10)
          write (iout,*) "wham_calc OK"
          call flush(iout)
          call write_dbase(islice,*10)
          write (iout,*) "write_dbase OK"
          call flush(iout)
          if (ensembles.gt.0) then
            call make_ensembles(islice,*10)
            write (iout,*) "make_ensembles OK"
            call flush(iout)
          endif
        endif
      enddo
#ifdef MPI
      call MPI_Finalize( IERROR )
#endif
      stop
   10 write (iout,*) "Error termination of the program"
      call MPI_Finalize( IERROR )
      stop
      end
c------------------------------------------------------------------------------
#ifdef MPI
      subroutine proc_groups
C Split the processors into the Master and Workers group, if needed.
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "DIMENSIONS.FREE"
      include "mpif.h"
      include "COMMON.IOUNITS"
      include "COMMON.MPI"
      include "COMMON.FREE"
      integer n,chunk,i,j,ii,remainder
      integer kolor,key,ierror,errcode
      logical lprint
      lprint=.true.
C 
C Split the communicator if independent runs for different parameter
C sets will be performed.
C
      if (nparmset.eq.1 .or. .not.separate_parset) then
        WHAM_COMM = MPI_COMM_WORLD
      else if (separate_parset) then
        if (nprocs.lt.nparmset) then
          write (iout,*) 
     & "*** Cannot split parameter sets for fewer processors than sets",
     &  nprocs,nparmset
          call MPI_Finalize(ierror)
          stop
        endif 
        write (iout,*) "nparmset",nparmset
        nprocs = nprocs/nparmset
        kolor = me/nprocs
        key = mod(me,nprocs)
        write (iout,*) "My old rank",me," kolor",kolor," key",key
        call MPI_Comm_split(MPI_COMM_WORLD,kolor,key,WHAM_COMM,ierror)
        call MPI_Comm_size(WHAM_COMM,nprocs,ierror)
        call MPI_Comm_rank(WHAM_COMM,me,ierror)
        write (iout,*) "My new rank",me," comm size",nprocs
        write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
     &   " WHAM_COMM",WHAM_COMM
        myparm=kolor+1
        write (iout,*) "My parameter set is",myparm
        call flush(iout)
      else
        myparm=nparmset
      endif
      Me1 = Me
      Nprocs1 = Nprocs
      return
      end
c------------------------------------------------------------------------------
      subroutine work_partition(islice,lprint)
c Split the conformations between processors
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "DIMENSIONS.FREE"
      include "mpif.h"
      include "COMMON.IOUNITS"
      include "COMMON.MPI"
      include "COMMON.PROT"
      integer islice
      integer n,chunk,i,j,ii,remainder
      integer kolor,key,ierror,errcode
      logical lprint
C
C Divide conformations between processors; the first and
C the last conformation to handle by ith processor is stored in 
C indstart(i) and indend(i), respectively.
C
C First try to assign equal number of conformations to each processor.
C
        n=ntot(islice)
        write (iout,*) "n=",n
        indstart(0)=1
        chunk = N/nprocs1
        scount(0) = chunk
c        print *,"i",0," indstart",indstart(0)," scount",
c     &     scount(0)
        do i=1,nprocs1-1
          indstart(i)=chunk+indstart(i-1) 
          scount(i)=scount(i-1)
c          print *,"i",i," indstart",indstart(i)," scount",
c     &     scount(i)
        enddo 
C
C Determine how many conformations remained yet unassigned.
C
        remainder=N-(indstart(nprocs1-1)
     &    +scount(nprocs1-1)-1)
c        print *,"remainder",remainder
C
C Assign the remainder conformations to consecutive processors, starting
C from the lowest rank; this continues until the list is exhausted.
C
        if (remainder .gt. 0) then 
          do i=1,remainder
            scount(i-1) = scount(i-1) + 1
            indstart(i) = indstart(i) + i
          enddo
          do i=remainder+1,nprocs1-1
            indstart(i) = indstart(i) + remainder
          enddo
        endif

        indstart(nprocs1)=N+1
        scount(nprocs1)=0

        do i=0,NProcs1
          indend(i)=indstart(i)+scount(i)-1
          idispl(i)=indstart(i)-1
        enddo

        N=0
        do i=0,Nprocs1-1
          N=N+indend(i)-indstart(i)+1
        enddo

c        print *,"N",n," NTOT",ntot(islice)
        if (N.ne.ntot(islice)) then
          write (iout,*) "!!! Checksum error on processor",me,
     &     " slice",islice
          call flush(iout)
          call MPI_Abort( MPI_COMM_WORLD, Ierror, Errcode )
        endif

      if (lprint) then
        write (iout,*) "Partition of work between processors"
          do i=0,nprocs1-1
            write (iout,'(a,i5,a,i7,a,i7,a,i7)')
     &        "Processor",i," indstart",indstart(i),
     &        " indend",indend(i)," count",scount(i)
          enddo
      endif
      return
      end
#endif
#ifdef AIX
      subroutine flush(iu)
      call flush_(iu)
      return
      end
#endif