#include "copyright.h"
#include "../include/dprec.fh"

#ifndef PBSA
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Emit the final minimization report 
!-----------------------------------------------------------------------
!     --- REPORT_MIN_RESULTS ---
!-----------------------------------------------------------------------
! Find the maximum component of the gradient (grdmax),
! emit this and other gradient details (printe),
! optionally do pseudocontact shift constraints,
! optionally decompose energies for mm_pbsa,
! and emit nmr related information (nmrptx).

subroutine report_min_results( nstep, gradient_rms, coordinates, &
      forces, energies, igraph, xx, ix, ih )

   use state
   use decomp, only : checkdec, printdec, printpdec
   use file_io_dat
   implicit none

   integer, intent(in)             :: nstep
   _REAL_,  intent(in)             :: gradient_rms
   _REAL_,  intent(in)             :: coordinates(*)
   _REAL_,  intent(in)             :: forces(*)
   type(state_rec), intent(in)     :: energies
   character(len=4), intent(in)    :: igraph(*)    ! atom name map
   _REAL_,  intent(inout)          :: xx(*)        ! real dynamic memory
   integer, intent(inout)          :: ix(*)        ! integer dynamic memory
   character(len=4), intent(inout) :: ih(*)        ! hollerith dynamic memory

#  include "box.h"
#  include "extra.h"
#  include "md.h"
#  include "memory.h"
#  include "nmr.h"
#  include "tgtmd.h"
#  include "multitmd.h"

   character(len=4) :: atom_name_of_gmax
   integer          :: atom_number_of_gmax
   _REAL_           :: gradient_max
   _REAL_ emtmd

   if (master) then
      write(6, '(/ /20x,a,/)' ) 'FINAL RESULTS'
      call grdmax( forces, gradient_max, atom_number_of_gmax )
      atom_name_of_gmax = igraph(atom_number_of_gmax)
      if (imin /= 5) rewind(MDINFO_UNIT)
      call printe( nstep, gradient_rms, gradient_max, energies, &
             atom_number_of_gmax, atom_name_of_gmax )
      if (idecomp > 0) call checkdec(idecomp)
      if (idecomp == 1 .or. idecomp == 2) call printdec(ix)
      if (idecomp == 3 .or. idecomp == 4) call printpdec(ix)
      if (nmropt > 0) then
         if (iredir(7) /= 0) call pcshift(-1,coordinates,forces)
         call nmrptx(6)
         call nmrptx(MDINFO_UNIT)
         call ndvptx(coordinates,forces,ih(m04),ih(m02),ix(i02),nres, &
               xx(l95),natom,xx(lwinv),ntb,xx(lnmr01),ix(inmr02),6)
      end if
      if (itgtmd == 2) then
        emtmd = 0.0d0
        call mtmdcall(emtmd,xx(lmtmd01),ix(imtmd02),coordinates,forces,ih(m04),ih(m02),ix(i02),&
                                      ih(m06),xx(lmass),natom,nres,'PRNT')
      end if
      if (imin /= 5) call amflsh(MDINFO_UNIT)
   end if

   return
end subroutine report_min_results


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Emit a minimization progress report 
!-----------------------------------------------------------------------
!     --- REPORT_MIN_PROGRESS ---
!-----------------------------------------------------------------------
! Find the maximum component of the gradient (grdmax),
! emit this and other gradient details (printe),
! and emit nmr related information (nmrptx).

subroutine report_min_progress( nstep, gradient_rms, forces, energies, &
                                igraph, charge )

   use state
   use file_io_dat
   use crg_reloc, only : ifcr, crprintcharges, cr_print_charge
   implicit none

   integer, intent(in)    :: nstep
   _REAL_,  intent(in)    :: gradient_rms
   _REAL_,  intent(in)    :: forces(*)
   type(state_rec), intent(in)    :: energies
   character(len=4), intent(in)    :: igraph(*)    ! atom name map
   _REAL_, intent(in)     :: charge(*)

#  include "extra.h"
#  include "md.h"
#  include "nmr.h"

   character(len=4) :: atom_name_of_gmax
   integer          :: atom_number_of_gmax
   _REAL_           :: gradient_max

   if (master) then
      call grdmax( forces, gradient_max, atom_number_of_gmax )
      atom_name_of_gmax = igraph(atom_number_of_gmax)
      if (imin /= 5) rewind(MDINFO_UNIT)
      call printe( nstep, gradient_rms, gradient_max, energies, &
             atom_number_of_gmax, atom_name_of_gmax )
      if (nmropt > 0) then
         call nmrptx(6)
         call nmrptx(MDINFO_UNIT)
      end if
      if ( ifcr > 0 .and. crprintcharges > 0 ) then
         call cr_print_charge( charge, nstep ) 
      end if
      if (imin /= 5) call amflsh(MDINFO_UNIT)
   end if

   return
end subroutine report_min_progress
#endif /*ifndef PBSA*/

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Compute the maximum gradient component and the corresponding atom
!-----------------------------------------------------------------------
!     --- GRDMAX ---
!-----------------------------------------------------------------------

subroutine grdmax( gradient, gradient_max, atom_number_of_gmax )

   use constants, only : zero
   implicit none
   _REAL_,  intent(in)  :: gradient(*)
   !     This is actually two-dimensional (3,natoms), but to enable
   !     vectorization on IA32 SSE platforms they are treated as
   !     one-dimensional; this may also improve software pipelining !

   _REAL_,  intent(out) :: gradient_max
   integer, intent(out) :: atom_number_of_gmax

#  include "memory.h"

   integer :: i
   _REAL_  :: gi

   gradient_max = ZERO
   atom_number_of_gmax = 1
   do i = 1,3*natom
      gi = abs(gradient(i))
      if (gi > gradient_max) then
         gradient_max = gi
         atom_number_of_gmax = i
      end if
   end do
   atom_number_of_gmax = (atom_number_of_gmax - 1)/3 + 1

   return
end subroutine grdmax


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Print status of a minimization calculation.
!-----------------------------------------------------------------------
!     --- PRINTE ---
!-----------------------------------------------------------------------
! Output the step number, the gradient rms, the max gradient component,
! and the atom label and atom number of the max gradient component.

subroutine printe( nstep, gradient_rms, gradient_max, ene, &
      atom_number_of_gmax, atom_name_of_gmax )
   
#ifdef APBS
   use file_io_dat
#endif
#ifdef RISMSANDER
   use sander_rism_interface, only : rismprm, RISM_NONE, RISM_FULL, RISM_INTERP,&
        rism_calc_type, rism_thermo_print
#endif
   use qmmm_module, only : qmmm_nml, qmmm_struct
   use cns_xref
   use state ! Access to energy_rec
   use charmm_mod, only : charmm_active
   use crg_reloc, only : ifcr
   use emap,only : temap
   use ff11_mod, only : cmap_active

   implicit none
   integer, intent(in)          :: nstep
   _REAL_,  intent(in)          :: gradient_rms
   _REAL_,  intent(in)          :: gradient_max
   type(state_rec), intent(in)  :: ene
   integer, intent(in)          :: atom_number_of_gmax
   character(len=4), intent(in) :: atom_name_of_gmax

#  include "md.h"
#  include "ew_mpole.h"
#  include "ew_cntrl.h"
#  include "tgtmd.h"

   _REAL_  epot,enonb,enele,ehbond,ebond,eangle,edihed,enb14,eel14,egb,epb
   _REAL_  econst,epolar,aveper,aveind,avetot,esurf,edisp,diprms,dipiter
   _REAL_  dipole_temp,escf,dvdl,enemap
#ifdef RISMSANDER
   _REAL_ :: erism
   _REAL_ :: pot_array(potential_energy_rec_len)
#endif /*RISMSANDER*/
   _REAL_ :: ect

   !     ----- Extract Energies. -----

   epot        = ene%pot%tot
   enonb       = ene%pot%vdw
   enele       = ene%pot%elec
   ehbond      = ene%pot%hbond
   egb         = ene%pot%gb 
   epb         = ene%pot%pb
   ebond       = ene%pot%bond
   eangle      = ene%pot%angle
   edihed      = ene%pot%dihedral
   enb14       = ene%pot%vdw_14
   eel14       = ene%pot%elec_14
   econst      = ene%pot%constraint
   epolar      = ene%pot%polar
   dvdl        = ene%pot%dvdl
   aveper      = ene%aveper
   aveind      = ene%aveind
   avetot      = ene%avetot
   esurf       = ene%pot%surf
   diprms      = ene%diprms
   dipiter     = ene%dipiter
   dipole_temp = ene%dipole_temp
   
   escf    = ene%pot%scf
   edisp   = ene%pot%disp
   enemap   = ene%pot%emap
#ifdef RISMSANDER
   erism   = ene%pot%rism
#endif /*RISMSANDER*/
   ect = ene%pot%ct


   

   write(6,9018)
   write(6,9028) nstep, epot, gradient_rms, gradient_max, &
         atom_name_of_gmax, atom_number_of_gmax
   write(6,9038) ebond,eangle,edihed
! CHARMM SPECIFIC ENERGY TERMS
   if (charmm_active) write(6,9039) ene%pot%angle_ub, &
                                    ene%pot%imp,      &
                                    ene%pot%cmap

#ifdef RISMSANDER
   if( igb == 0 .and. ipb == 0 .and. rismprm%irism == 0) then
#else
   if( igb == 0 .and. ipb == 0 ) then
#endif
      write(6,9048) enonb,enele,ehbond
   else if ( igb == 10 .or. ipb /= 0 ) then
      write(6,9050) enonb,enele,epb
#ifdef APBS
   else if ( igb == 6 .and. mdin_apbs ) then
      write(6,9050) enonb,enele,epb
#endif /* APBS */
#ifdef RISMSANDER
   else if(rismprm%irism == 1 )then
      write(6,9051) enonb,enele,erism
#endif
   else
      write(6,9049) enonb,enele,egb
   end if
   write(6,9058) enb14,eel14,econst

   if ( ifcr /= 0 ) write(6,9099) ect

!  wxw: EMAP energy
   if (temap) then
      write (6,9062) enemap
   endif

   if (qmmm_nml%ifqnt) then
     !write the SCF energy
     if (qmmm_nml%qmtheory%PM3) then
        if (qmmm_nml%qmmm_int == 3) then
           write(6,9090) escf ! PM3-MM*
        else if (qmmm_nml%qmmm_int == 4) then
           write(6,9096) escf ! PM3/MMX2
        else
           write(6,9080) escf
        end if
     else if (qmmm_nml%qmtheory%AM1) then
        write(6,9081) escf
     else if (qmmm_nml%qmtheory%AM1D) then
        write(6,9981) escf
     else if (qmmm_nml%qmtheory%MNDO) then
        write(6,9082) escf
     else if (qmmm_nml%qmtheory%MNDOD) then
        write(6,9982) escf
     else if (qmmm_nml%qmtheory%PDDGPM3) then
        write(6,9083) escf
     else if (qmmm_nml%qmtheory%PDDGMNDO) then
        write(6,9084) escf
     else if (qmmm_nml%qmtheory%PM3CARB1) then
        write(6,9085) escf
     else if (qmmm_nml%qmtheory%DFTB) then
        write(6,9086) escf
     else if (qmmm_nml%qmtheory%RM1) then
        write(6,9087) escf
     else if (qmmm_nml%qmtheory%PDDGPM3_08) then
        write(6,9088) escf
     else if (qmmm_nml%qmtheory%PM6) then
        write(6,9089) escf
     else if (qmmm_nml%qmtheory%PM3ZNB) then
        write(6,9091) escf
     else if (qmmm_nml%qmtheory%EXTERN) then
        write(6,9092) escf
     else if (qmmm_nml%qmtheory%PM3MAIS) then
        write(6,9093) escf
     else
        write(6,'(" ERROR - UNKNOWN QM THEORY")')
     end if
   end if
#ifdef PUPIL_SUPPORT
   write(6,9900) escf
#endif
   if( gbsa > 0 ) write(6,9077) esurf
   if (igb == 10 .or. ipb /= 0) write(6,9074) esurf,edisp
   call write_cns_xref_min_energies( ene )
#ifdef APBS
   if (igb == 6 .and. mdin_apbs ) write(6,9069) esurf
#endif /* APBS */
      if (cmap_active .and. ipol > 0 ) then
          write(6,9066) epolar, ene%pot%cmap
      else
          if (cmap_active) write(6,9067) ene%pot%cmap
          if (ipol > 0) write(6,9068) epolar
      end if
   if (econst /= 0.0) write(6,9078) epot-econst
   if ( dvdl /= 0.d0) write(6,9100) dvdl
   if (induced > 0.and.indmeth < 3) write(6,9190)diprms,dipiter
   if (induced > 0.and.indmeth == 3) write(6,9191)diprms, &
         dipole_temp
   if (itgtmd == 1) then
      write (6,'(a,f8.3)') "Current RMSD from reference: ",rmsdvalue
      write (6,'(a,f8.3)') "Current target RMSD:         ",tgtrmsd
   end if

   call amflsh(6)

   !     ----- SEND IT TO THE INFO FILE -----

   if (imin /= 5) then
      write(7,9018)
      write(7,9028) nstep, epot, gradient_rms, gradient_max, &
            atom_name_of_gmax, atom_number_of_gmax
      write(7,9038) ebond,eangle,edihed
! CHARMM SPECIFIC ENERGY TERMS
   if (charmm_active) write(7,9039) ene%pot%angle_ub, &
                                    ene%pot%imp,      &
                                    ene%pot%cmap


#ifdef RISMSANDER
      if( igb == 0 .and. ipb == 0 .and. rismprm%irism == 0) then
#else
      if( igb == 0 .and. ipb == 0 ) then
#endif
         write(7,9048) enonb,enele,ehbond
      else if ( igb == 10 .or. ipb /= 0 ) then
         write(7,9050) enonb,enele,epb
#ifdef APBS
      else if ( igb == 6 .and. mdin_apbs ) then
         write(7,9050) enonb,enele,epb
#endif /* APBS */
#ifdef RISMSANDER
      else if ( rismprm%irism == 1 ) then
         write(7,9051) enonb,enele,erism
#endif
      else
         write(7,9049) enonb,enele,egb
      end if
      write(7,9058) enb14,eel14,econst

!  wxw: EMAP energy
   if (temap) then
      write (7,9062) enemap
   endif

      if (qmmm_nml%ifqnt) then
        !write the SCF energy
        if (qmmm_nml%qmtheory%PM3) then
          if (qmmm_nml%qmmm_int == 3) then
             write(7,9090) escf ! PM3-MM*
          else if (qmmm_nml%qmmm_int == 4) then
             write(7,9096) escf ! PM3/MMX2
          else
             write(7,9080) escf
          end if
        else if (qmmm_nml%qmtheory%AM1) then
           write(7,9081) escf
        else if (qmmm_nml%qmtheory%AM1D) then
           write(7,9981) escf
        else if (qmmm_nml%qmtheory%MNDO) then
           write(7,9082) escf
        else if (qmmm_nml%qmtheory%MNDOD) then
           write(7,9982) escf           
        else if (qmmm_nml%qmtheory%PDDGPM3) then
           write(7,9083) escf
        else if (qmmm_nml%qmtheory%PDDGMNDO) then
           write(7,9084) escf
        else if (qmmm_nml%qmtheory%PM3CARB1) then
           write(7,9085) escf
        else if (qmmm_nml%qmtheory%DFTB) then
           write(7,9086) escf
        else if (qmmm_nml%qmtheory%RM1) then
           write(7,9087) escf
        else if (qmmm_nml%qmtheory%PDDGPM3_08) then
           write(7,9088) escf
        else if (qmmm_nml%qmtheory%PM6) then
           write(7,9089) escf
        else if (qmmm_nml%qmtheory%PM3ZNB) then
           write(7,9091) escf
        else if (qmmm_nml%qmtheory%EXTERN) then
           write(7,9092) escf
        else if (qmmm_nml%qmtheory%PM3MAIS) then
           write(7,9093) escf
        else
           write(7,'(" ERROR - UNKNOWN QM THEORY")')
        end if
      end if
#ifdef PUPIL_SUPPORT
      write(7,9900) escf
#endif
      if( gbsa > 0 ) write(7,9077) esurf
      if ( igb == 10 .or. ipb /= 0 ) write(7,9074) esurf,edisp
#ifdef APBS
      if (igb == 6 .and. mdin_apbs ) write(7,9069) esurf
#endif /* APBS */
! FF11 CMAP SPECIFIC ENERGY TERMS
      if (cmap_active .and. epolar /= 0.0 ) then
          write(7,9066) epolar, ene%pot%cmap
      else
          if (cmap_active) write(7,9067) ene%pot%cmap
          if (epolar /= 0.0) write(7,9068) epolar
      end if
      if (econst /= 0.0) write(7,9078) epot-econst
      if ( dvdl /= 0.d0) write(7,9100) dvdl
      if (induced > 0.and.indmeth < 3) write(7,9190)diprms,dipiter
      if (induced > 0.and.indmeth == 3) write(7,9191)diprms, &
            dipole_temp
   end if

#ifdef RISMSANDER
   if(rismprm%irism==1 .and. rismprm%write_thermo==1)then
      if(rism_calc_type(nstep) == RISM_FULL)&
           call rism_thermo_print(.false.,transfer(ene%pot,pot_array))
   end if
#endif /*RISMSANDER*/

   9018 format (/ /,3x,'NSTEP',7x,'ENERGY',10x,'RMS',12x,'GMAX',9x, &
         'NAME',4x,'NUMBER')
   9028 format(1x,i6,2x,3(2x,1pe13.4),5x,a4,2x,i7,/)
   9038 format (1x,'BOND    = ',f13.4,2x,'ANGLE   = ',f13.4,2x, &
         'DIHED      = ',f13.4)
   9039 format (1x, 'UB      = ', f13.4, 2x, 'IMP     = ', f13.4, 2x, &
            'CMAP       = ', f13.4)
   9048 format (1x,'VDWAALS = ',f13.4,2x,'EEL     = ',f13.4,2x, &
         'HBOND      = ',f13.4)
   9049 format (1x,'VDWAALS = ',f13.4,2x,'EEL     = ',f13.4,2x, &
         'EGB        = ',f13.4)
   9050 format (1x,'VDWAALS = ',f13.4,2x,'EEL     = ',f13.4,2x, &
         'EPB        = ',f13.4)

#ifdef RISMSANDER
   9051 format (1x,'VDWAALS = ',f13.4,2x,'EEL     = ',f13.4,2x, &
         'ERISM      = ',f13.4)
#endif

   9058 format (1x,'1-4 VDW = ',f13.4,2x,'1-4 EEL = ',f13.4,2x, &
         'RESTRAINT  = ',f13.4)
   9062 format (1x,'EMAP    = ',f14.4)
   9066 format (1x,'EPOLAR  = ',f13.4,2x,'CMAP    = ',f13.4)
   9067 format (1x,'CMAP    = ',f13.4)
   9068 format (1x,'EPOLAR  = ',f13.4)
#ifdef APBS
   9069 format (1x,'ENPOLAR = ',f13.4)
#endif /* APBS */
   9074 format (1x,'ECAVITY = ',f13.4,2x,'EDISPER = ',f13.4)
   9077 format (1x,'ESURF   = ',f13.4)
   9078 format (1x,'EAMBER  = ',f13.4)
   9080 format (1x,'PM3ESCF =',f14.4)
   9081 format (1x,'AM1ESCF =',f14.4)
   9981 format (1x,'AM1DESCF =',f14.4)
   9082 format (1x,'MNDOESCF=',f14.4)
   9982 format (1x,'MNDODESCF=',f14.4)
   9083 format (1x,'PDDGPM3-ESCF=',f14.4)
   9084 format (1x,'PDDGMNDO-ESCF=',f14.4)
   9085 format (1x,'PM3CARB1-ESCF=',f14.4)
   9086 format (1x,'DFTBESCF=',f14.4)
   9087 format (1x,'RM1ESCF =',f14.4)
   9088 format (1x,'PDDGPM3_08-ESCF=',f14.4)
   9089 format (1x,'PM6ESCF =',f14.4)
   9090 format (1x,'PM3MMXESCF =',f14.4)
   9091 format (1x,'PM3ZNBESCF =',f14.4)
   9092 format (1x,'EXTERNESCF =',f14.4)
   9093 format (1x,'PM3MAISESCF =',f14.4)
   9099 format (1x,'ECRG    =',f14.4)
   9096 format (1x,'PM3MMX2ESCF =',f14.4)
   9190 format(1x,'Dipole convergence: rms = ',e10.3,' iters = ',f6.2)
   9191 format(1x,'Dipole convergence: rms = ',e10.3, &
         ' temperature = ',f6.2)
   9100 format (1x,'DV/DL  = ',f14.4)
#ifdef PUPIL_SUPPORT
   9900 format (1x,'PUPESCF =',f14.4)
#endif
   return
end subroutine printe