! modified.or.of $NWCHEM_TOP/src/nwchem.F
! originally NWChem version 6.3, updated to match v6.6 
! you.lu@stfc.ac.uk
! April 2017
! YL: refactered April 2019 to the modern.or.of v6.8
! YL 12/01/2021: changed the version number to 7.0.2

      integer function chemshell_nwchem(natoms,
     &                                  nbqs,
     &                                  coords,
     &                                  charges,
     &                                  bq_coords,
     &                                  bq_charges,
     &                                  bq_grads,
     &                                  energy,
     &                                  gradients,
     &                                  rtdb,
     &                                  debug) result(ierror)

          use ChemShellNWChemModule  , only : chemsh_nwchem_input, chemsh_nwchem_output, fragmovecs, NWChemObj
          use ChemShellParallelModule, only : chemsh_parallel

          implicit none

! YL: use Fortran directives to switch to fixed-format F77 legacies
!YL: DEC$ NOFREEFORM
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "rtdb.fh"
!#include "tcgmsg.fh"
#include "pstat.fh"
#include "util.fh"
#include "inp.fh"
#include "bgj_common.fh"
#include "stdio.fh"
#include "bq.fh"
#include "bq_params.fh"
!YL: DEC$ FREEFORM


          integer     , intent(in)  :: natoms, nbqs
          real(kind=8), intent(in)  :: coords(3,natoms)
          real(kind=8), intent(in)  :: charges(natoms)
          real(kind=8), intent(in)  :: bq_coords(3,nbqs)
          real(kind=8), intent(in)  :: bq_charges(nbqs)
          real(kind=8), intent(out) :: bq_grads(3,nbqs)
          real(kind=8), intent(out) :: energy
          real(kind=8), intent(out) :: gradients(3,natoms)
          integer     , intent(out) :: rtdb
          integer     , intent(in)  :: debug

          integer :: stack, heap, global
          logical :: status
#ifdef CRAY_T3D
          integer oldact, fsigctl
#endif
#ifdef PSCALE
          integer io_code
#else
          integer*4 io_code
#endif

          integer           :: i, ios, icounter
          character(len=32) :: key, theory
          logical           :: otmp
          logical, external :: bq_init, bq_end

          ! YL 15/12/2023: excited state
          integer :: nroots
          real(kind=8), allocatable :: excitation_energies(:)
!ENDCHEMSH


!  $Id$

!  ======================================================================================================
!> \mainpage Northwest Computational Chemistry Package (NWChem) @NWCHEM_VERSION@
!>
!> NWChem is an open-source computational chemistry package distributed under the terms of
!> the Educational Community License (ECL) 2.0
!>
!> This software and its documentation were developed at the EMSL at Pacific Northwest National Laboratory,
!> a multiprogram national laboratory, operated for the U.S. Department of Energy by Battelle under
!> Contract Number DE-AC05-76RL01830. Support for this work was provided by the Department of Energy
!> Office of Biological and Environmental Research, Office of Basic Energy Science, and the Office of
!> Advanced Scientific Computing.
!>
!> Licensed under the Educational Community License, Version 2.0 (the "License"); you may
!> not use this file except in compliance with the License. You may obtain a copy of the
!> License at <a href="http://www.osedu.org/licenses/ECL-2.0">http://www.osedu.org/licenses/ECL-2.0</a>.
!>
!> Unless required by applicable law or agreed to in writing, software distributed under the
!> License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
!> either express or implied. See the License for the specific language governing
!> permissions and limitations under the License.
!>
!> Further information, including user documentation and forums, may be found at
!> <a href="http://www.nwchem-sw.org/">http://www.nwchem-sw.org/</a>. Alternatively,
!> the paper
!>
!> * M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, T.P. Straatsma, H.J.J. Van Dam,
!>   D. Wang, J. Nieplocha, E. Apra, T.L. Windus, W.A. de Jong (2010)<br>
!>   "NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations"<br>
!>   <I>Computer Physics Communications</I>, <b>181</b>, 1477–1489, DOI: <a href="http://dx.doi.org/10.1016/j.cpc.2010.04.018">10.1016/j.cpc.2010.04.018</a>
!>
!> provides details on the codes capabilities.
!>
!> Copyright (c) 1994-2018 Pacific Northwest National Laboratory, Battelle Memorial Institute
!>
!> Environmental Molecular Sciences Laboratory (EMSL)<br>
!> Pacific Northwest National Laboratory<br>
!> Richland, WA 99352

!  ======================================================================================================

          character(len=NW_MAX_PATH_LEN) :: input_filename, output_filename
          character(len=NW_MAX_PATH_LEN) :: rtdb_name = ""
          double precision               :: total_wall, total_cpu
#ifdef USE_OFFLOAD
          integer  ppnout
          logical, external :: offload_master
#endif

          logical           :: ostartup, ocontinue, orestart
          logical, external :: input_parse
          integer, external :: set_geometry
! YL TODO: OPENMP
#if defined(USE_OPENMP)
          integer, external :: omp_get_max_threads
#endif

          integer           :: nprescfs

          luout = 2016
          ierror = 1

          if(chemsh_parallel%is_root()) then
              write(input_filename , '(a)') chemsh_nwchem_input
              write(output_filename, '(a)') chemsh_nwchem_output
          endif

          ! initialise the output file
          if(chemsh_parallel%is_root().and.output_filename.ne.'stdout') then
              open(unit=luout, status='unknown', form='formatted', file=output_filename)
          ! reopen the output file to avoid conflict with other modules (e.g., DL_FIND)
      !    elseif(chemsh_parallel%is_root().and.output_filename.ne.'stdout') then
      !        open(unit=luout, status='unknown', form='formatted', position='append', file=output_filename)
          endif

          ! MXINIT is needed by PeIGS and PFFT to initialize
          ! the communication fabric they use.
          call mxinit

          ! initialize timers so they are relative to job start
          total_wall = -util_wallsec()
          total_cpu  = -util_cpusec()

          ! only process 0 opens the input file
          if(chemsh_parallel%is_root()) then   
              open(luin, file=input_filename, status='old', form='formatted', iostat=ios)
              if(ios.ne.0) then
                  call errquit('nwchem: failed to open the input file', 0, INPUT_ERR)
              endif
          else
#if !(defined(KSR) || defined(IBM) || defined(FUJITSU_SOLARIS) ||defined(FUJITSU_VPP) ||defined(FUJITSU_VPP64))
          close(luin)
#endif
          endif

          ! Touch OpenMP here so that any runtime initialization happens up-front.
          ! This ensures that any printout that the OpenMP runtime generates,
          ! such as affinity information, appears at the top of the output file.
          ! Otherwise, it might not appear until e.g. the CCSD module, at which
          ! point it will pollute the output file in an undesirable way.
          ! Do not move this in front of GA/MPI/TCGMSG initialization, since the
          ! OpenMP runtime may inherit affinity information from MPI that is only
          ! determined during MPI initialization.
          ! Format definition is outside of preprocessor protection to ensure the
          ! label is not accidentally reused, since that will not be caught by
          ! testing that does not enable OpenMP.
99    format(2x,'NWChem w/ OpenMP: maximum threads = ',i2)
#if defined(USE_OPENMP)
          !$omp parallel
          !$omp master
          if(chemsh_parallel%is_root()) then
              write(luout,99) omp_get_max_threads()
          endif
          !$omp end master
          !$omp end parallel
#endif

          call rtdb_init()

#if defined(DECOSF) || defined(LINUXALPHA)
          call dec_fpe  ! To avoid underflow problems on Alpha in Texas
#endif
#ifdef CRAY_T3D
          ! this as a temporary fix for SIGFPE in Texas that does not seem
          ! to affect the final results
          oldact = fsigctl('IGNORE','SIGFPE',0)
#endif
      
          ! initialize pstat
          if(.not.pstat_init(20, 1, ' ')) then
              call errquit('nwchem: pstat_init failed', 0, UNKNOWN_ERR)
          endif

          ! YL: we had to patch the original input_file_info() which contained a bug of using nodeid() instead of ga_nodeid()
          !     fixed version released as NWChem 6.7
          call input_file_info(input_filename, rtdb_name, ostartup, ocontinue)

          ! now are ready to summarize the environment
          call nwchem_banner(input_filename, rtdb_name, ostartup, ocontinue)

          ! actually open the database and store the file prefix 
          ! note that only process 0 has the database name ... that is OK.
          if(ostartup) then
              if(.not.rtdb_open(rtdb_name, 'empty', rtdb)) then
                  call errquit('start: rtdb_open empty failed ', 0, RTDB_ERR)
              endif
          else
              if(.not.rtdb_open(rtdb_name, 'old', rtdb)) then
                  call errquit('start: rtdb_open old failed ', 0, RTDB_ERR)
              endif
          endif

          ! initialize nxtask
          call nxtask_init(rtdb)

          ! BGJ
          bgj_rtdb = rtdb

          if(ostartup.or.ocontinue) then
            orestart = .false.
          else
            orestart = .true.
          endif

          call util_set_rtdb_state(rtdb,ostartup,ocontinue,orestart)

          call util_file_info_rtdb(rtdb) ! Save file info for restart

          call movecs_ecce_print_on()

          call geom_hnd_parallel(.true.)

          call perfm_start()

#ifdef USE_OFFLOAD
          call util_getppn(ppnout)
          if(ppnout.eq.0) then
              call errquit('util_getppn failed ', 0, UERR)
          endif
          if(chemsh_parallel%is_root()) then
              write(luout, *) chemsh_parallel%node, ' ppn ', ppnout
          endif
          if(offload_master()) then
              if(chemsh_parallel%node.lt.ppnout) then
                  write(luout, *) chemsh_parallel%node, ' offload_master '
              endif
              call util_mic_set_affinity()
          endif
          ierror = chemsh_parallel%sync()
#endif

          if(orestart.or.ocontinue) then
              call nw_print_restart_info(rtdb)
          endif

          ! if continue then go right to task stored on rtdb do not further parse
          ! input.  if input is required then user should have used restart
          if(ocontinue) then
              call task(rtdb)
          endif

          nprescfs = 0

          ! parse input data, shove into the database and execute tasks
          ! YL 01/05/2019: reduced the previous itask-based procedure to enable
          !                multiple tasks in a single input, so we can pre-calculate
          !                fragment MOs at SCF level
          do while(input_parse(rtdb))

              call util_print_rtdb_load(rtdb,' ') ! High level print

              ! initialise BQ handling (taken from task/task.F)
              if(nbqs.gt.0) then
                  if(.not.bq_init(rtdb)) then
                      call errquit(" >>> ChemShell/NWChem ERROR: failed to initialise BQs", 0, CALC_ERR)
                  endif
              endif

              if(util_print('tcgmsg', print_never)) then
                  call setdbg(1)
              else
                  call setdbg(0)
              endif

#ifdef CATAMOUNT
              call util_allocga()
#endif

              ierror = set_geometry(rtdb, coords, bq_coords, bq_charges, natoms, nbqs, nprescfs, debug)

              call task(rtdb)

              ! YL 01/05/2019: reactivate BQs
              if(nbqs.gt.0) then
                  call bq_grads_direct(rtdb, natoms, nbqs, bq_grads, debug)
              endif

              ! YL 01/05/2019: there may be methods better than this
              if(rtdb_get(rtdb, 'scf:energy', MT_DBL, 1, energy)) then
                  if(chemsh_parallel%is_root()) then
                      print *
                      print *, "ChemShell/NWChem: HF Energy =", energy, "a.u."
                      print *
                  endif
                  if(rtdb_get(rtdb, 'scf:gradient', MT_DBL, 3*natoms, gradients)) then
      !                print *, "### SCF gradients calculated"
                  endif
              else
                  if(rtdb_get(rtdb, 'dft:energy', MT_DBL, 1, energy)) then
                      if(chemsh_parallel%is_root()) then
                          print *
                          print *, "ChemShell/NWChem: DFT Energy =", energy, "a.u."
                          print *
                      endif
                      if(.not.rtdb_get(rtdb, 'dft:gradient', MT_DBL, 3*natoms, gradients)) then
      !                    print *, "### DFT gradients calculated"
                      endif
                  else
                      call errquit(" >>> ChemShell/NWChem ERROR: no energy", 0, CALC_ERR)
                  endif
              endif

          enddo

          ! end bq handling
          if(nbqs.ne.0) then
              if(.not.bq_end(rtdb)) then
                  call errquit(" >>> ChemShell/NWChem ERROR: failed to finalise BQs", 0, CALC_ERR)
              endif
          endif

          ! YL 15/12/2023: TD-DFT excited states
          if(rtdb_get(rtdb,'tddft:nroots',MT_INT,1,nroots)) then
              allocate(excitation_energies(nroots))
              if(rtdb_get(rtdb,'tddft:energy-all',MT_DBL,nroots,excitation_energies)) then
                  NWChemObj%excitation_energies = excitation_energies + energy
              else
                  call errquit(" >>> ChemShell/NWChem ERROR: Failed to retrieve TD-DFT energies", 0, CALC_ERR)
              endif
              ! KS - 20/02/24: in case that no gradients are requested for TD-DFT
              if(.not.rtdb_get(rtdb,'tddft:gradient',MT_DBL,3*natoms, gradients)) then
              !    call errquit(" >>> ChemShell/NWChem ERROR: Failed to retrieve TD-DFT gradients", 0, CALC_ERR)
              endif
              deallocate(excitation_energies)
          endif

          ! close the RTDB
          call util_print_rtdb_load(rtdb, ' ') ! high level print

          if(util_print('rtdbvalues', print_debug)) then
              if(.not.rtdb_print(rtdb, .true.)) then
                  call errquit('control: rtdb_print failed', 0, RTDB_ERR)
              endif
          elseif(util_print('rtdb', print_high)) then
              if(.not.rtdb_print(rtdb, .false.)) then
                  call errquit('control: rtdb_print failed', 0, RTDB_ERR)
              endif
          endif

          if(.not.rtdb_close(rtdb, 'keep')) then
               if(chemsh_parallel%is_root()) then
                   print *, ">>> ChemShell/NWChem WARNING: nwchem rtdb_close failed!"
               endif
          endif

          if(util_print('rtdb', print_high).or.util_print('rtdbvalues', print_high)) then
              call rtdb_print_usage  ! Called after closing so memory leaks apparent
          endif

          ! tidy up pstat
          if(.not.pstat_terminate()) then
              call errquit('nwchem: pstat_terminate failed', 0, UNKNOWN_ERR)
          endif

          ! print memory and other info
          ierror = chemsh_parallel%sync()

          if(chemsh_parallel%is_root()) then
              if(util_print('ga summary', print_default)) then
                  ierror = chemsh_parallel%summarise(int(0, kind=4))
              endif
              if(util_print('ga stats', print_default)) then
                  ierror = chemsh_parallel%print()
                  write(luout, *)
              endif
              if(util_print('ma summary', print_default).and.debug.gt.3) then
                  call MA_SUMMARIZE_ALLOCATED_BLOCKS()
              endif
              if(util_print('ma stats', print_high).and.debug.gt.3) then
                  call MA_PRINT_STATS(.true.)
              elseif(util_print('ma stats', print_default).and.debug.gt.3) then
                  call MA_PRINT_STATS(.false.)
              endif
              if(util_print('version', print_debug)) then
                  call util_version
              endif
              if(util_print('citation', print_none)) then
                  call util_cite()
              endif
          endif

          call perfm_end()
          total_wall = total_wall + util_wallsec()
          total_cpu  = total_cpu  + util_cpusec()

          if(chemsh_parallel%is_root()) then
              if(util_print('total time', print_none)) then
                  write(luout, 1, iostat=io_code) total_cpu, total_wall
1     format(/' Total times  cpu:', f11.1, 's     wall:', f11.1, 's')
                  if(io_code.ne.0.and.util_print('total time', print_high))  then
                      write(luout,*) ' iostat ',io_code
                  endif
                  call util_flush(luout)
              endif
          endif
          call ecce_print1('all tasks cpu time' , MT_DBL,total_cpu , 1)
          call ecce_print1('all tasks wall time', MT_DBL,total_wall, 1)

          ierror = chemsh_parallel%sync()

          ! close the ECCE' output file if there is one
          call ecce_print_file_close()

          ! YL 05/08/2019: it wasn't closed and caused serious problems
          if(chemsh_parallel%is_root().and.output_filename.ne.'stdout') then
              close(unit=luin)
              close(unit=luout)
          endif

          ierror = 0

      endfunction chemshell_nwchem
      subroutine nwchem_banner(input_filename, rtdb_name, ostartup, ocontinue)

          use ChemShellParallelModule, only : chemsh_parallel

          implicit none

!YL: DEC$ NOFREEFORM
#include "global.fh"
#include "inp.fh"
#include "mafdecls.fh"
#include "util.fh"
#include "stdio.fh"
!YL: DEC$ FREEFORM

          ! print a banner at the top of the output to identify the origin
          ! of the code along with hostname, no. of processors, date, ...

          character*(*) input_filename, rtdb_name
          logical ostartup
          logical ocontinue

          real(kind=8)     , parameter   :: ndbls2MB   = 8.0/1024.0/1024.0
          real(kind=8)     , parameter   :: Byte2MB    = 1.0/1024.0/1024.0
          integer          , parameter   :: ndbls2Byte = 8
          character(len=80)              :: hostname, compiled, batchid
          character(len=NW_MAX_PATH_LEN) :: executable, srcdir, file_prefix
          character(len=NW_MAX_PATH_LEN) :: branch, nwchem_rev, ga_rev
          character(len=26)              :: date
          character(len=2048)            :: strbuff
#if defined(CRAY)
          integer ilen
#endif
          integer heap, stack, global, global_bytes, total, nprocs, ierror
          logical status
#ifdef CATAMOUNT
          integer istatus,setvbuf3f
          external setvbuf3f
#endif
          logical, external :: util_scalapack_info
          character(len=3)  :: cstatus, fstatus
          character(len=10) :: cstart
          character(len=40) :: GStatus

#ifdef CATAMOUNT
          istatus = setvbuf3f(luout, 0, 4096)
          istatus = setvbuf3f(0, 0, 4096)
#endif

          if(chemsh_parallel%is_root()) then

              call util_print_centered(luout, '                                                          ', 40, .true.)
              write(luout,*)
              write(luout,*)
              call util_print_centered(luout, 'Northwest Computational Chemistry Package (NWChem) @NWCHEM_VERSION@'  , 40, .true.)
              write(luout,*)
              write(luout,*)
              call util_print_centered(luout, 'Environmental Molecular Sciences Laboratory'               , 40, .false.)
              call util_print_centered(luout, 'Pacific Northwest National Laboratory'                     , 40, .false.)
              call util_print_centered(luout, 'Richland, WA 99352'                                        , 40, .false.)
              write(luout,*)
              call util_print_centered(luout, 'Copyright (c) 1994-2020'                                   , 40, .false.)
              call util_print_centered(luout, 'Pacific Northwest National Laboratory'                     , 40, .false.)
              call util_print_centered(luout, 'Battelle Memorial Institute'                               , 40, .false.)
              write(luout,*)
              call util_print_centered(luout, 'NWChem is an open-source computational chemistry package'  , 40, .false.)
              call util_print_centered(luout, 'distributed under the terms of the'                        , 40, .false.)
              call util_print_centered(luout, 'Educational Community License (ECL) 2.0'                   , 40, .false.)
              call util_print_centered(luout, 'A copy of the license is included with this distribution'  , 40, .false.)
              call util_print_centered(luout, 'in the LICENSE.TXT file'                                   , 40, .false.)
              write(luout,*)
              call util_legal()

              call util_getarg(0, executable)
              call util_hostname(hostname)
              call util_date(date)

              nprocs = chemsh_parallel%nnodes_wg

              if(ostartup) then
                  cstart = 'startup'
              else
                  if(ocontinue) then
                      cstart = 'continue'
                  else
                      cstart = 'restart'
                  endif
              endif

              call util_file_prefix(' ', file_prefix)

              compiled = ""
              call util_nwchem_srcdir(srcdir)
              srcdir = srcdir(1:max(1,inp_strlen(srcdir)-4))

              branch = "ChemShell/NWChem Interface"
              call util_nwchem_version(nwchem_rev)

              call util_ga_version(ga_rev)

              ! see if we can get the batch system id.  We understand NQS and
              ! loadLeveler at present
              BatchID = ' '
              call util_getenv('QSUB_REQID', batchid)
              if(inp_strlen(batchid).eq.0) then
                  call util_getenv('LOADL_JOB_NAME', batchid)
                  if ( inp_strlen(batchid) .eq. 0 ) then
                      call util_getenv('LSB_JOBID', batchid)
                  endif
              endif

              call util_print_centered(luout, 'Job information', 17, .true.)

 999          format(/4x,"hostname        = ",a/ 
     &                4x,"program         = ",a/ 
     &                4x,"date            = ",a/ 
     &                4x,"compiled        = ",a/ 
     &                4x,"source          = ",a/ 
     &                4x,"nwchem branch   = ",a/ 
     &                4x,"nwchem revision = ",a/ 
     &                4x,"ga revision     = ",a/ 
     &                4x,"use scalapack   = ",l1/
     &                4x,"input           = ",a/ 
     &                4x,"prefix          = ",a/ 
     &                4x,"data base       = ",a/ 
     &                4x,"status          = ",a/ 
     &                4x,"nprocs          = ",i8/
     &                4x,"time left       = ",i6,"s")
              write(luout,999) hostname(1:inp_strlen(hostname)),
     &                         executable(1:inp_strlen(executable)),
     &                         date(1:inp_strlen(date)),
     &                         compiled(1:inp_strlen(compiled)),
     &                         srcdir(1:inp_strlen(srcdir)),
     &                         branch(1:inp_strlen(branch)),
     &                         nwchem_rev(1:inp_strlen(nwchem_rev)),
     &                         ga_rev(1:inp_strlen(ga_rev)),
     &                         util_scalapack_info(),
     &                         input_filename(1:inp_strlen(input_filename)),
     &                         file_prefix(1:inp_strlen(file_prefix)),
     &                         rtdb_name(1:inp_strlen(rtdb_name)),
     &                         cstart(1:inp_strlen(cstart)),
     &                         nprocs,
     &                         util_batch_job_time_remaining()

              ! print version information to ecce.out file
              compiled(inp_strlen(compiled)+1:inp_strlen(compiled)+25) = ' Version @NWCHEM_VERSION@'
              call ecce_print_version(compiled(1:inp_strlen(compiled)))

              if(inp_strlen(batchid).gt.0) then
                  strbuff = '(4x,"batch job id    = ",a/)'
                  write(luout,strbuff) batchid(:inp_strlen(batchid))
              else
                  write(luout,'(//)') 
              endif

              ! determine, without altering, setting of memory verification and hardfail
              status = MA_SET_AUTO_VERIFY(.true.)
              if(status) then
                  cstatus = 'yes'
              else
                  cstatus = 'no '
              endif
              status = MA_SET_AUTO_VERIFY(.true.)

              status = MA_SET_HARD_FAIL(.true.)
              if(status) then
                  fstatus = 'yes'
              else
                  fstatus = 'no '
              endif
              status = MA_SET_HARD_FAIL(status)

              ! memory size information
              ! note:  subtleties of MA & friends will likely result in actual
              ! heap and stack sizes _slightly_ larger than requested at initialization,
              ! and other slight "inconsistencies". it has to do with the overhead associated
              ! with an allocation

              ! available in MA (in number of doubles)
              heap  = MA_INQUIRE_HEAP(MT_DBL)
              stack = MA_INQUIRE_STACK(MT_DBL)

              ! (in Byte)
              global_bytes = GA_MEMORY_AVAIL()
              ! if GA cannot determine a limit on the available memory, it
              ! will return a negative number, which MA_SIZEOF does not
              ! appreciate. this should not happen in NWChem, but we might
              ! as well play it safe.
              if(global_bytes.ge.0) then
                  ! in number of doubles
                  global = MA_SIZEOF(MT_BYTE, global_bytes, MT_DBL)
              else
                  global = 0
              endif

              if(GA_USES_MA()) then
                  total = heap + stack
                  gstatus = '(using local memory)'
              else
                  total = heap + stack + global
                  gstatus = '(using shared memory)'
              endif

              call util_print_centered(luout, 'Memory information', 19, .true.)

              strbuff = '(/4x,"heap     = ", i10," doubles = ",f8.3," MB",/,
     &                     4x,"stack    = ", i10," doubles = ",f8.3," MB",/,
     &                     4x,"global   = ", i10," doubles = ",f8.3," MB",
     &                     1x,a/,
     &                     4x,"total    = ", i10," doubles = ",f8.3," MB",/,
     &                     4x,"verify   = ", a3/,
     &                     4x,"hardfail = ", a3//)'
              write(luout,strbuff)
     &              heap,  ((dble(heap)  *ndbls2MB)),
     &              stack, ((dble(stack) *ndbls2MB)),
     &              global,((dble(global)*ndbls2MB)) ,
     &              gstatus(1:inp_strlen(gstatus)),
     &              total, ((dble(total) *ndbls2MB)),
     &              cstatus,
     &              fstatus

              call util_print_centered(luout, 'Directory information', 20, .true.)
              write(luout, *)

          endif

          ! below here back to parallel execution
          call util_file_print_dirs()
          if(chemsh_parallel%is_root()) then
              write(luout, *)
          endif

      endsubroutine nwchem_banner
      subroutine nw_print_restart_info(rtdb)

          use ChemShellParallelModule, only : chemsh_parallel

          implicit none

!YL: DEC$ NOFREEFORM
#include "global.fh"
#include "errquit.fh"
#include "inp.fh"
#include "mafdecls.fh"
#include "rtdb.fh"
#include "stdio.fh"
#include "util.fh"
!YL: DEC$ FREEFORM
          integer :: ierror, rtdb

          ! summarize the status of a restart calculation
          logical :: mode, task_qmmm, ignore, lstatus
          character*30 operation, status
          character*32 theory

          if(chemsh_parallel%is_root()) then
              mode = rtdb_parallel(.false.)
              if(.not.rtdb_get(rtdb,'task:QMMM', mt_log, 1, task_qmmm)) then
                  task_qmmm = .false.
              endif
              if(.not.rtdb_cget(rtdb, 'task:operation', 1, operation)) then
                  operation = ' '
              endif
              if(.not.rtdb_cget(rtdb,'task:theory',1,theory)) then
                  theory = ' '
              endif
              if(.not.rtdb_get(rtdb, 'task:ignore', mt_log, 1, ignore)) then
                  ignore = .false.
              endif
              if(rtdb_get(rtdb, 'task:status', mt_log, 1, lstatus)) then
                  if(lstatus) then
                      status = 'ok'
                  else
                      status = 'fail'
                  endif
              else
                  status = 'unknown'
              endif

              if(theory.ne.' '.or.operation.ne.' '.or.task_qmmm) then
                  call util_print_centered(luout,'Previous task information', 22, .true.)
                  write(luout,1) theory, operation, status, task_qmmm, ignore
1     format(/,
     &       '    Theory    = ', a,/,
     &       '    Operation = ', a,/,
     &       '    Status    = ', a,/,
     &       '    Qmmm      = ', l1,/,
     &       '    Ignore    = ', l1,/)
              endif
              call geom_print_known_geoms(rtdb)
              call bas_print_known_bases(rtdb)
              call util_flush(luout)
              mode = rtdb_parallel(mode)
          endif

          ierror = chemsh_parallel%sync()

      endsubroutine nw_print_restart_info
      ! directly access bq gradients instead of through force.dat 
      ! cf bq/bq_force.F
      subroutine bq_grads_direct(rtdb, natoms, nbqs, bq_grads, debug)

          use ChemShellParallelModule, only : chemsh_parallel

          implicit none

!YL: DEC$ NOFREEFORM
#include "rtdb.fh"
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "bq_params.fh"
#include "stdio.fh"
!YL: DEC$ FREEFORM

          integer     , intent(in)  :: rtdb
          integer     , intent(in)  :: natoms
          integer     , intent(in)  :: nbqs
          real(kind=8), intent(out) :: bq_grads(3,nbqs)
          integer     , intent(in)  :: debug

          character(len=32)  :: pname
          integer            :: bq_handle, bq_ncent
          logical            :: ostatus

          integer            :: i, j
          integer            :: i_cbq, i_qbq
          integer            :: h_f
          integer            :: i_f
          character(len=255) :: fname

          logical, external  :: bq_init, bq_end
          logical, external  :: bq_ncenter, bq_index_coord, bq_index_charge, bq_destroy, bq_get_active
          integer            :: ibq, ioff

          pname = "bq_grads_direct"

          ! YL 01/05/2019: after task(rtdb) is run BQs becomes inactive
          if(.not.bq_init(rtdb)) then
      !        if(debug.ge.5.and.chemsh_parallel%is_root()) then
      !            print *, " >>> ChemShell/NWChem bq_grads_direct: cannot reactivate BQs, quitting..."
      !        endif
              return
          endif

          call bq_force_status(rtdb, ostatus)

          if(.not.ostatus) then
      !        if(debug.ge.5.and.chemsh_parallel%is_root()) then
      !            print *, " >>> ChemShell/NWChem bq_grads_direct: bq_force_status = .false., quitting..."
      !        endif
              return
          endif

          if(.not.rtdb_cget(rtdb, "bq:force_filename", 1, fname)) then
              call util_file_prefix("_nwchem.bq_grads", fname)
          endif

          ! get bq information
          if(.not.bq_get_active(bq_handle)) then
      !        if(debug.ge.5.and.chemsh_parallel%is_root()) then
      !            print *, " >>> ChemShell/NWChem bq_grads_direct: no active BQs, quitting..."
      !        endif
              return
          endif

          if(.not.bq_ncenter(bq_handle, bq_ncent)) then
              call errquit(pname//':no bq centers', 0, 0)
          endif

          if(bq_ncent.ne.nbqs) then
              call errquit(pname//':nbqs /= bq_ncent!', 0, 0)
          endif

          if(.not.bq_index_coord(bq_handle, i_cbq)) then
              call errquit(pname//':no bq coords', 0, 0)
          endif

          if(.not.bq_index_charge(bq_handle, i_qbq)) then
              call errquit(pname//':no bq coords', 0, 0)
          endif

          if(.not.MA_PUSH_GET(MT_DBL, 3*nbqs, 'bq f', h_f, i_f)) then
              call errquit(pname//': Failed to allocate memory for c', 3*nbqs, MA_ERR)
          endif

          call qmmm_bq_forces(rtdb,
     &                        nbqs,
     &                        DBL_MB(i_cbq),
     &                        DBL_MB(i_qbq),
     &                        DBL_MB(i_f))

          ! return gradients
          bq_grads = reshape(DBL_MB(i_f:i_f+3*nbqs-1), (/3,nbqs/))

          if(.not.MA_POP_STACK(h_f)) then
              call errquit(pname//': Failed to deallocate stack ', 0, MA_ERR)
          endif

          if(.not.bq_end(rtdb)) then
              print *, " >>> ChemShell/NWChem bq_grads_direct: bq_end(rtdb) failed."
          endif

      endsubroutine bq_grads_direct
      integer function set_geometry(rtdb, coords, bq_coords, bq_charges, natoms, nbqs, nprescfs, debug) result(ierror)

          use ChemShellNWChemModule  , only : fragmovecs
          use ChemShellParallelModule, only : chemsh_parallel

          implicit none

!YL: DEC$ NOFREEFORM
#include "errquit.fh"
#include "rtdb.fh"
#include "pstat.fh"
#include "util.fh"
#include "inp.fh"
#include "stdio.fh"
#include "geom.fh"
#include "bq.fh"
#include "bq_params.fh"
!YL: DEC$ FREEFORM

          integer     , intent(inout) :: rtdb
          real(kind=8), intent(in)    :: coords(3,natoms)
          real(kind=8), intent(in)    :: bq_coords(3,nbqs)
          real(kind=8), intent(in)    :: bq_charges(nbqs)
          integer     , intent(in)    :: natoms
          integer     , intent(in)    :: nbqs
          integer     , intent(inout) :: nprescfs
          integer     , intent(in)    :: debug

          real(kind=8)                :: twodimdbl(3,natoms+nbqs), onedimdbl(natoms+nbqs)
          character(len=16)           :: tags(natoms+nbqs)
          integer                     :: geom_handle, bq_handle, ncentres

          if(chemsh_parallel%is_root().and.debug.ge.3) then
              print *
              print *, ">>> ChemShell/NWChem: setting up geometry..."
              print *
          endif

          ! create `geom` object
          if(.not.geom_create(geom_handle, 'geometry')) then
              ierror = -1
              return
          endif

          ! load geometry from RTDB
          if(.not.geom_rtdb_load(rtdb, geom_handle, 'geometry')) then
              ierror = -2
              return
          endif

          ! get geometry information from RTDB
          tags = ''
          if(.not.geom_cart_get(geom_handle, ncentres, tags, twodimdbl, onedimdbl)) then
              ierror = -3
              return
          endif

          ! note number of centres no longer includes BQs
          ! YL 01/05/2019: this was to check if natoms are consistent, but since we need
          !                to support pre-SCF now we don't throw an error when they
          !                are not identical
          if(ncentres.ne.natoms.and.chemsh_parallel%is_root()) then
              ! YL 10/04/2019: fragmovecs is a switch from ChemShellNWChemModule
              if(fragmovecs) then
                  nprescfs = nprescfs + 1
                  if(ncentres.gt.1) then
                      write(*,'(A,I0,A,I0,A)') " ChemShell/NWChem: processing fragment MO No. ",
     &                                         nprescfs, " of ", ncentres, " atoms"
                  else
                      write(*,'(A,I0,A,I0,A)') " ChemShell/NWChem: processing fragment MO No. ",
     &                                         nprescfs, " of ", ncentres, " atom"
                  endif
              else
                  print *, ">>> ChemShell/NWChem: ncentres =", ncentres
                  print *, "                      natoms   =", natoms
                  ierror = -10
                  return
              endif
          endif 

          twodimdbl(:,1       :natoms)      = coords
          twodimdbl(:,natoms+1:natoms+nbqs) = bq_coords
          ! set geom
          if(.not.geom_cart_set(geom_handle, ncentres, tags, twodimdbl, onedimdbl)) then
              ierror = -4
              return
          endif

          ! save to RTDB
          if(.not.geom_rtdb_store(rtdb, geom_handle, 'geometry')) then
              ierror = -5
              return
          endif

          ! delete
          if(.not.geom_destroy(geom_handle)) then
              ierror = -6
              return
          endif

          ! we need to duplicate the functionality above
          ! for BQs. so we will find the active handle, set new values 
          ! and then replace the old values in the RTDB with this information.
          ! This has been copied from nwchem-6.3/src/bq/bq_input.F
          ! YL: in PCS we pass in bq_coords and bq_charges from the Python side, so we don't use bq_get()
          if(nbqs.gt.0) then

              if(.not.bq_get_active(bq_handle)) then
                  call errquit('bq_get_active: no active bq handle', 0, 0)
              endif

              ! set BQs for rtdb
              if(.not.bq_set(bq_handle, nbqs, bq_charges, reshape(bq_coords, (/nbqs*3/)))) then
                  call errquit('bq_set failed', 0, 0)
              endif

              call bq_print_info(bq_handle)

              ! overwrite the stored BQ points
              if(.not.bq_rtdb_store(rtdb, bq_handle)) then
                  call errquit('bq_store failed', 0, 0)
              endif

          endif

          ierror = 0

      endfunction set_geometry
