program ENE_OPT
c Optimize the UNRES energy function by minimization of a quartic target
c function or by the VMC method.
      implicit none
#ifndef ISNAN
      external proc_proc
#endif
#ifdef WINPGI
cMS$ATTRIBUTES C ::  proc_proc
#endif
      include "DIMENSIONS.ZSCOPT"
#ifdef MPI
      include "mpif.h"
      integer IERROR,ERRCODE,kolor,key
      include "COMMON.MPI"
#endif
      include "COMMON.IOUNITS"
      include "COMMON.OPTIM"
      integer nvarr,iparm
      double precision rr,x(max_paropt)
      integer idumm
      integer i
c      print *,"Starting..."
#ifdef MPI
c      print *,"Initializing MPI..."
      call MPI_Init( IERROR )
      call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
      call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
c      print *,"Finished initializing MPI..."
      Master = 0
      Master1 = 1
c      print *,"Me",me," Master",master," Ierror",ierror
      if (ierror.gt.0) then
        write(iout,*) "SEVERE ERROR - Can't initialize MPI."
        call mpi_finalize(ierror)
        stop
      endif
#endif
#ifndef ISNAN
c NaNQ initialization
      i=-1
      rr=dacos(100.0d0)
#ifdef WINPGI
      idumm=proc_proc(rr,i)
#else
      call proc_proc(rr,i)
#endif
#endif
      call initialize
c      print *,"calling openunits"
      call openunits
c      print *,"openunits called"
      call read_general_data(*10)
      write (iout,'(80(1h-)/10x,
     & "Maximum likelihood optimization of UNRES energy function",
     & " v. 05/10/16"/80(1h-))')
      call flush(iout)
      call cinfo
c      call promienie(*10)
      write (iout,*) "Finished READ_GENERAL_DATA"
      call flush(iout)
      do iparm=1,nparmset
        call parmread(iparm,*10)
      enddo
      call read_pmf_data(*10)
      write (iout,*) "Finished parmread"
      call flush(iout)
      call read_optim_parm(*10)
      call print_general_data(*10)
      call read_protein_data(*10)
      write (iout,*) "Finished READ_PROTEIN_DATA"
      call flush(iout)
      call read_database(*10)
      write (iout,*) "Finished READ_DATABASE"
      call flush(iout)
#ifdef MPI 
c      write (iout,*) Me,' calling PROC_GROUPS'
      call proc_groups
c      write (iout,*) Me,' calling WORK_PARTITION_MAP'
c      call work_partition_map(nvarr)
#endif
#ifdef NEWCORR
      call pmfread(*10)
#endif
      call proc_data(nvarr,x,*10)
      call read_thermal
#ifdef MPI
      if (me.eq.Master) then
#endif
      call maxlikopt(nvarr,x)
#ifdef MPI
      call jebadelko(nvarr)
      else
        call jebadelko(nvarr)
      endif
      call bilans
      call MPI_Finalize( IERROR )
#else
      call bilans
#endif
      stop
   10 write (iout,*) "Error termination of the program"
      call MPI_Finalize( IERROR )
      stop
      end
c------------------------------------------------------------------------------
      subroutine bilans
      implicit none
      integer ntasks 
      parameter (ntasks=11)
      include "DIMENSIONS.ZSCOPT"
#ifdef MPI
      include "mpif.h"
      include "COMMON.MPI"
      integer IERROR
      double precision ttask_all(ntasks)
      integer nctask_all(ntasks)
#endif
      include "COMMON.IOUNITS"
      integer i
      double precision ttask
      integer nctask
      common /timing/ ttask(ntasks),nctask(ntasks)
      character*16 tname(ntasks) /"function","gradient",9*''/

      write (iout,*)
      write (iout,'(80(1h-))')
#ifdef MPI
      write (iout,*) "Routine call info from the processor ",me," ..."
#else
      write (iout,*) "Routine call info ..."
#endif
      write (iout,*)
      write (iout,'(65(1h-))')
      write (iout,100) "task","   # calls","     total time",
     &  "  ave.time/call"
      write (iout,'(65(1h-))')
      do i=1,ntasks
        write (iout,200) tname(i),nctask(i),ttask(i),
     &  ttask(i)/(nctask(i)+1.0d-10)
      enddo
      write (iout,*) 
#ifdef MPI
      call MPI_Reduce(ttask(1),ttask_all(1),ntasks, 
     &  MPI_DOUBLE_PRECISION, MPI_SUM, Master, WHAM_COMM,IERROR)
      call MPI_Reduce(nctask(1),nctask_all(1),ntasks, 
     &  MPI_INTEGER, MPI_SUM, Master, WHAM_COMM,IERROR)
      if (Me.eq.Master) then
      write (iout,'(80(1h-))')
      write (iout,*) "Routine call info from all processors ..."
      write (iout,*)
      write (iout,'(65(1h-))')
      write (iout,100) "task","   # calls","     total time",
     &  "  ave.time/call"
      write (iout,'(65(1h-))')
      do i=1,ntasks
        write (iout,200) tname(i),nctask_all(i),ttask_all(i),
     &    ttask_all(i)/(nctask_all(i)+1.0d-10)
      enddo
      write (iout,*) 
      endif
#endif
      return
  100 format(a,t21,a,t31,a,t46,a)
  200 format(a,t21,i10,f15.2,f15.8)
      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 "mpif.h"
      include "COMMON.IOUNITS"
      include "COMMON.MPI"
      include "COMMON.VMCPAR"
      integer n,chunk,iprot,i,j,ii,remainder
      integer kolor,key,ierror,errcode
      logical lprint
      lprint=.true.
C No secondary structure constraints. 
      Comm1 = WHAM_COMM
      Me1 = Me
      Nprocs1 = Nprocs
      return
      end
c-------------------------------------------------------------------------------
      subroutine work_partition(lprint)
c Split the conformations between processors
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "mpif.h"
      include "COMMON.CLASSES"
      include "COMMON.IOUNITS"
      include "COMMON.MPI"
      include "COMMON.VMCPAR"
      integer n,chunk,iprot,i,j,ii,remainder
      integer kolor,key,ierror,errcode
      logical lprint
C
C Divide conformations between processors; for each proteins the first and
C the last conformation to handle by ith processor is stored in 
C indstart(i,iprot) and indend(i,iprot), respectively.
C
C First try to assign equal number of conformations to each processor.
C
      do iprot=1,nprot
        n=ntot_work(iprot)
        write (iout,*) "Protein",iprot," n=",n
        indstart(0,iprot)=1
        chunk = N/nprocs1
        scount(0,iprot) = chunk
c        print *,"i",0," indstart",indstart(0,iprot)," scount",
c     &     scount(0,iprot)
        do i=1,nprocs1-1
          indstart(i,iprot)=chunk+indstart(i-1,iprot) 
          scount(i,iprot)=scount(i-1,iprot)
c          print *,"i",i," indstart",indstart(i,iprot)," scount",
c     &     scount(i,iprot)
        enddo 
C
C Determine how many conformations remained yet unassigned.
C
        remainder=N-(indstart(nprocs1-1,iprot)
     &    +scount(nprocs1-1,iprot)-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,iprot) = scount(i-1,iprot) + 1
            indstart(i,iprot) = indstart(i,iprot) + i
          enddo
          do i=remainder+1,nprocs1-1
            indstart(i,iprot) = indstart(i,iprot) + remainder
          enddo
        endif

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

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

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

c        print *,"N",n," NTOT",ntot_work(iprot)
        if (N.ne.ntot_work(iprot)) then
          write (iout,*) "!!! Checksum error on processor",me
          call flush(iout)
          call MPI_Abort( WHAM_COMM, Ierror, Errcode )
        endif

        ii=0
        do i=1,ntot_work(iprot)
          if (i.ge.indstart(me1,iprot) .and. i.le.indend(me1,iprot)) 
     &    then
            ii=ii+1
            i2ii(i,iprot)=ii
          else
            i2ii(i,iprot)=0
          endif
c          write (iout,*) "i",i," iprot",iprot," i2ii",i2ii(i,iprot)
        enddo

      enddo ! iprot
      if (lprint) then
        write (iout,*) "Partition of work between processors"
        do iprot=1,nprot
          write (iout,*) "Protein",iprot
#ifdef DEBUG
          write (iout,*) "The i2ii array"
          do j=1,ntot_work(iprot)
            write (iout,*) j,i2ii(j,iprot)
          enddo
#endif
          do i=0,nprocs1-1
            write (iout,'(a,i5,a,i7,a,i7,a,i7)')
     &        "Processor",i," indstart",indstart(i,iprot),
     &        " indend",indend(i,iprot)," count",scount(i,iprot)
          enddo
        enddo
      endif
      return
      end
c------------------------------------------------------------------------------
      subroutine jebadelko(nvarr)
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "mpif.h"
      include "COMMON.IOUNITS"
      include "COMMON.MPI"
      include "COMMON.VMCPAR"
      include "COMMON.CLASSES"
      include "COMMON.OPTIM"
      include "COMMON.WEIGHTS"
      include "COMMON.WEIGHTDER"
      include "COMMON.ENERGIES"
      include "COMMON.TIME1"
      include "COMMON.PROTNAME"
      include "COMMON.PROTFILES"
      include "COMMON.TORSION"
      include "COMMON.COMPAR"
      integer What, TAG, IERROR, status(MPI_STATUS_SIZE), istop, iprot, 
     &  nvarr, nf, errcode, ider
      integer count
      double precision x(max_paropt), g(max_paropt), viol
      integer uiparm,i,j
      integer ibatch,ib
      external fdum
      double precision rdum,rdif
      double precision tcpu,t1,t1w,t1_ini,t1w_ini
      logical lprn,op
      integer iun
      lprn=.false.
      write(iout,*) "Processor",me,me1," called JEBADELKO"
      call flush(iout)
      if (me.eq.Master) then
        istop=0
        call func1(nvarr,istop,x)
      else
#ifdef DEBUG
      write (iout,*) "ELOWEST at slave  starting JEBADELKO"
      do iprot=1,nprot
        do ibatch=1,natlike(iprot)+2
          do ib=1,nbeta(ibatch,iprot)
            write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
     &         " elowest",elowest(ib,ibatch,iprot)
          enddo
        enddo
      enddo
#endif
        istop=1
        t1w_ini = MPI_WTIME()
        t1_ini = tcpu()
        do while (istop.ne.0)
#ifdef DEBUG
      write (iout,*) "ELOWEST at slave  calling FUNC1 from JBADELKO"
      do iprot=1,nprot
        do ibatch=1,natlike(iprot)+2
          do ib=1,nbeta(ibatch,iprot)
            write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
     &         " elowest",elowest(ib,ibatch,iprot)
          enddo
        enddo
      enddo
#endif
          call func1(nvarr,istop,x)
c          write (iout,*) "slave: after func1"
c          call flush(iout)
#ifdef NEWCORR
          if (istop.eq.1 .and. mod_fourier(nloctyp).gt.0) then
            rdum = rdif(nvarr,x,g,ider)
c            write (iout,*) "slave: after rdif"
c            call flush(iout)
          endif
#endif
        enddo
        t1w = mpi_wtime() - t1w_ini
        t1 = tcpu() - t1_ini
        write (iout,*)
        write (iout,*) "CPU time",t1," wall clock time",t1w
        call flush(iout)
      endif
#ifdef DEBUG
      write (iout,*)
      write (iout,*) "Energies of processor",me
      do j=1,nprot
        write (iout,*)
        write (iout,*) "Protein ",protname(iprot)
        do i=1,ntot_work(j)
          if (i.ge.indstart(me1,j).and.i.le.indend(me1,j)) then
            write(iout,*)i,e_total(i,j),rmstb(i,j),iscore(i,0,j)
          endif
        enddo
      enddo
#endif
      do iprot=1,nprot
        write (iout,*) "Deleting scratchfile",bprotfiles(iprot)
        inquire (file=bprotfiles(iprot),number=iun,opened=op)
        write (iout,*) "unit",iun," icbase",icbase
        if (.not. op) then
          open (icbase,file=bprotfiles(iprot),status="old",
     &    form="unformatted",access="direct",recl=lenrec(iprot))
          close(icbase,status="delete")
        else
          close(iun,status="delete")
        endif
        call flush(iout)
        if (.not.mod_other_params) then
        write (iout,*) "Deleting scratchfile",benefiles(iprot)
        inquire (file=benefiles(iprot),number=iun,opened=op)
        write (iout,*) "unit",iun," ientout",icbase
        if (.not.op) then
          open (ientout,file=benefiles(iprot),status="old",
     &    form="unformatted",access="direct",recl=lenrec_ene(iprot))
          close(ientout,status="delete")
        else
          close (iun,status="delete")
        endif
        endif
        call flush(iout)
      enddo
      write (iout,*) "Processor",me,"leaves JEBADELKO"
      call flush(iout)
      return
      end
#endif