program wham_multparm
!      program WHAM_multparm
! Creation/update of the database of conformations
      use wham_data
      use io_wham
      use io_database
      use wham_calc
      use ene_calc
      use conform_compar
      use work_part
!
      use io_units
      use control_data, only:indpdb
#ifdef MPI
      use mpi_data
!      use mpi_
#endif
      use control, only:initialize
!el      use io_config, only:parmread
!
#ifndef ISNAN
      external proc_proc
#ifdef WINPGI
!MS$ATTRIBUTES C ::  proc_proc
#endif
#endif
!
!el#ifndef ISNAN
!el      external proc_proc
!el#endif
!el#ifdef WINPGI
!elcMS$ATTRIBUTES C ::  proc_proc
!el#endif
!      include "DIMENSIONS"
!      include "DIMENSIONS.ZSCOPT"
!      include "DIMENSIONS.FREE"
!      implicit none
#ifdef MPI
!      include "COMMON.MPI"
!      use mpi_data
      include "mpif.h"
      integer :: IERROR,ERRCODE
#endif
!      include "COMMON.IOUNITS"
!      include "COMMON.FREE"
!      include "COMMON.CONTROL"
!      include "COMMON.ALLPARM"
!      include "COMMON.PROT"
      real(kind=8) :: rr !,x(max_paropt)
      integer :: idumm
      integer :: i,ipar,islice

!el      run_wham=.true.
!#define WHAM_RUN
!      call alloc_wham_arrays
!write(iout,*) "after alloc wham"
#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
!el      if (nprocs.gt.MaxProcs+1) then
!el        write (2,*) "Error - too many processors",&
!el         nprocs,MaxProcs+1
!el        write (2,*) "Increase MaxProcs and recompile"
!el        call MPI_Finalize(IERROR)
!el        stop
!el      endif
#endif
! 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
!write(iout,*) "before init"
      call initialize
!write(iout,*)"after init"
      call openunits
!write(iout,*)"after open ui"
      call cinfo
!write(iout,*)"after cinfo"
      call read_general_data(*10)
!write(iout,*)"after read_gen"
      call flush(iout)
      call molread(*10)
!write(iout,*)"after molread"
      call flush(iout)
#ifdef MPI 
      write (iout,*) "Calling proc_groups"
      call proc_groups
      write (iout,*) "proc_groups exited"
      call flush(iout)
#endif
!el----------
      call alloc_wham_arrays
!el----------
      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)
!write(iout,*)"before proc_cont, define frag"
        call proc_cont
        call fragment_list
        if (constr_dist.gt.0) call read_dist_constr
      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
          write (iout,*) "call enecalc",islice,nslice
          call enecalc(islice,*10)
          write (iout,*) "enecalc OK"
          call flush(iout)
          call WHAMCALC(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"
#ifdef MPI
      call MPI_Finalize( IERROR )
#endif
      stop
      end program wham_multparm
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
#ifdef MPI
      subroutine proc_groups
! Split the processors into the Master and Workers group, if needed.
      use io_units
      use MPI_data
      use wham_data

      implicit none
!      include "DIMENSIONS"
!      include "DIMENSIONS.ZSCOPT"
!      include "DIMENSIONS.FREE"
!      include "COMMON.IOUNITS"
!      include "COMMON.MPI"
!      include "COMMON.FREE"
      include "mpif.h"
      integer :: n,chunk,i,j,ii,remainder
      integer :: kolorW,key,ierror,errcode
      logical :: lprint
      lprint=.true.
!
! Split the communicator if independent runs for different parameter
! sets will be performed.
!
      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
        kolorW = me/nprocs
        key = mod(me,nprocs)
        write (iout,*) "My old rank",me," kolor",kolorW," key",key
        call MPI_Comm_split(MPI_COMM_WORLD,kolorW,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=kolorW+1
        write (iout,*) "My parameter set is",myparm
        call flush(iout)
      else
        myparm=nparmset
      endif
      Me1 = Me
      Nprocs1 = Nprocs
      return
      end subroutine proc_groups
#endif
!------------------------------------------------------------------------------
#ifdef AIX
      subroutine flush(iu)
      call flush_(iu)
      return
      end subroutine flush
#endif
!-----------------------------------------------------------------------------
      subroutine promienie(*)

      use io_units
      use names
      use io_base, only:ucase
      use energy_data, only:sigma0,dsc,dsc_inv
      use wham_data
      use w_compar_data
      implicit none
!      include 'DIMENSIONS'
!      include 'COMMON.CONTROL'
!      include 'COMMON.INTERACT'
!      include 'COMMON.IOUNITS'
!      include 'COMMON.CONTPAR'
!      include 'COMMON.LOCAL'
      integer ::i,j
      real(kind=8) :: facont=1.569D0  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
      character(len=8) :: contfunc
      character(len=8) :: contfuncid(5)=reshape((/'GB      ',&
                      'DIST    ','CEN     ','ODC     ','SIG     '/),shape(contfuncid))
!el      character(len=8) ucase
      call getenv("CONTFUNC",contfunc)
      contfunc=ucase(contfunc)
      do icomparfunc=1,5
        if (contfunc.eq.contfuncid(icomparfunc)) goto 10
      enddo     
   10 continue
      write (iout,*) "Sidechain contact function is ",contfunc,&
        "icomparfunc",icomparfunc 
      do i=1,ntyp
        do j=1,ntyp
          if (icomparfunc.lt.3) then
            read(isidep1,*) chi_comp(i,j),chip_comp(i,j),sig_comp(i,j),&
             sc_cutoff(i,j)
          else if (icomparfunc.lt.5) then
            read(isidep1,*) sc_cutoff(i,j)
          else if (icomparfunc.eq.5) then
            sc_cutoff(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)*facont
          else
            write (iout,*) "Error - Unknown contact function"
            return 1
          endif
        enddo
      enddo
      close (isidep1)
      do i=1,ntyp1
        if (i.eq.10 .or. i.eq.ntyp1) then
          dsc_inv(i)=0.0d0
        else
          dsc_inv(i)=1.0d0/dsc(i)
        endif
      enddo
      return
      end subroutine promienie
!-----------------------------------------------------------------------------
      subroutine alloc_wham_arrays

      use names
      use geometry_data, only:nres
      use energy_data, only:maxcont
      use wham_data
      use w_compar_data
      integer :: i,j,k,l
!-------------------------
! COMMON.FREE
!      common /wham/
      allocate(stot(nslice)) !(maxslice)
      do i=1,nslice
        stot(i)=0
      enddo
      allocate(Kh(nQ,MaxR,MaxT_h,nParmSet),q0(nQ,MaxR,MaxT_h,nParmSet))!(MaxQ,MaxR,MaxT_h,max_parm)
      allocate(f(maxR,maxT_h,nParmSet)) !(maxR,maxT_h,max_parm)
      allocate(beta_h(maxT_h,nParmSet)) !(MaxT_h,max_parm)
      allocate(nR(maxT_h,nParmSet),nRR(maxT_h,nParmSet)) !(maxT_h,max_parm)
      allocate(snk(MaxR,MaxT_h,nParmSet,nSlice)) !(MaxR,MaxT_h,max_parm,MaxSlice)
!      do i=1,MaxR
!        do j=1,MaxT_h
!          do k=1,nParmSet
!            do l=1,nSlice
!              snk(i,j,k,l)=0
!            enddo
!          enddo
!        enddo
!      enddo

      allocate(totraj(maxR,nParmSet)) !(maxR,max_parm)

      allocate(nT_h(nParmSet))!(max_parm)
      allocate(replica(nParmSet))
      allocate(umbrella(nParmSet))
      allocate(read_iset(nParmSet))
!      allocate(nT_h(nParmSet))
!-------------------------
! COMMON.PROT
!      common /protein/
      allocate(ntot(nslice))  !(maxslice)
!      allocatable :: isampl        !(max_parm)
!-------------------------
! COMMON.PROTFILES
!      common /protfil/
      allocate(protfiles(maxfile_prot,2,MaxR,MaxT_h,nParmSet)) !(maxfile_prot,2,MaxR,MaxT_h,Max_Parm)
      allocate(nfile_bin(MaxR,MaxT_h,nParmSet))
      allocate(nfile_asc(MaxR,MaxT_h,nParmSet))
      allocate(nfile_cx(MaxR,MaxT_h,nParmSet))
      allocate(rec_start(MaxR,MaxT_h,nParmSet))
      allocate(rec_end(MaxR,MaxT_h,nParmSet)) !(MaxR,MaxT_h,Max_Parm)
!-------------------------
! COMMON.OBCINKA
!      common /obcinka/
      allocate(time_start_collect(maxR,MaxT_h,nParmSet))
      allocate(time_end_collect(maxR,MaxT_h,nParmSet)) !(maxR,MaxT_h,Max_Parm)
!-------------------------
! COMMON.CONTPAR
!      common /contpar/
      allocate(sig_comp(ntyp,ntyp),chi_comp(ntyp,ntyp),&
        chip_comp(ntyp,ntyp),sc_cutoff(ntyp,ntyp)) !(ntyp,ntyp)
!-------------------------
! COMMON.PEPTCONT
!      common /peptcont/
      allocate(icont_pept_ref(2,maxcont)) !(2,maxcont)
!      allocate(ncont_frag_ref()) !(mmaxfrag)
!      allocate(icont_frag_ref(2,maxcont)) !(2,maxcont,mmaxfrag)
      allocate(isec_ref(nres)) !(maxres)
!-------------------------
! COMMON.VAR
! Angles from experimental structure
!      common /varref/
      allocate(vbld_ref(nres),theta_ref(nres),&
        phi_ref(nres),alph_ref(nres),omeg_ref(nres)) !(maxres)
!-------------------------
      end subroutine alloc_wham_arrays
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------