#include "copyright.h"
#include "../include/dprec.fh"
#include "../include/assert.fh"
#include "ncsu-config.h"

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ The Molecular Dynamics/NMR Refinement/Modeling Module of the AMBER
!-----------------------------------------------------------------------
!     --- SANDER ---
!-----------------------------------------------------------------------

subroutine sander()

   use state
#ifndef DISABLE_NCSU
   use ncsu_sander_hooks, only : &
      ncsu_on_sander_init => on_sander_init, &
      ncsu_on_sander_exit => on_sander_exit
#endif /* DISABLE_NCSU */

   use lmod_driver
   use constants, only : INV_AMBER_ELECTROSTATIC, plumed, plumedfile
   ! The main qmmm_struct contains all the QMMM variables and arrays
   use qmmm_module, only : qmmm_nml, qmmm_struct, deallocate_qmmm, qmmm_mpi, &
#ifdef MPI
         qmmm_mpi_setup, &
#endif
         qm2_struct, qmewald, qm_gb, qmmm_vsolv, qm2_params

   use qmmm_vsolv_module, only: qmmm_vsolv_store_parameters, new

   use qm2_extern_module, only: qm2_extern_finalize
   use sebomd_module, only : sebomd_obj, &
       sebomd_open_files, sebomd_close_files, &
#ifdef MPI
        sebomd_bcast_obj, &
#endif
        sebomd_setup

   use sebomd_arrays, only : init_sebomd_arrays, cleanup_sebomd_arrays

   use genborn
   use decomp, only : allocate_int_decomp, allocate_real_decomp, &
         deallocate_int_decomp, deallocate_real_decomp, &
#ifdef MPI
         synchronize_dec, build_dec_mask, decmask, &
#endif
         nat, nrs, jgroup, indx
   use fastwt
   use relax_mat
   use nmr, only: nmrrad, impnum
   use ew_recip, only: deallocate_m1m2m3,first_pme
   use parms
   use molecule, only : mol_info, n_iwrap_mask_atoms, iwrap_mask_atoms, &
                        allocate_molecule, deallocate_molecule
   use nblist, only:cutoffnb,skinnb,nblist_allocate,nblist_deallocate, &
         nblist_allreal,nblist_allint, num_calls_nblist, first_list_flag
   use stack
   use amoeba_runmd, only : AM_RUNMD_get_coords,AM_RUNMD
   use amoeba_mdin, only : beeman_integrator,iamoeba,am_nbead
   use amoeba_interface, only : AMOEBA_deallocate,AMOEBA_readparm

#ifdef RISMSANDER
   use sander_rism_interface, only: rismprm,rism_setparam, rism_init,rism_finalize 
#endif

#ifdef PUPIL_SUPPORT
   use pupildata
#endif /* PUPIL */

#ifdef APBS
   use apbs
#endif /* APBS */

#ifdef _XRAY
   use xray_interface_module, only: xray_init, xray_read_parm, xray_fini
#endif

#ifdef MPI /* SOFT CORE */
   use softcore, only: setup_sc, cleanup_sc, ifsc, extra_atoms, sc_sync_x, &
        summarize_ti_changes, sc_check_perturbed_molecules, ti_check_neutral, tishake
   use mbar, only: setup_mbar, cleanup_mbar, ifmbar
#endif

   ! for LIE calculations
   use linear_response, only: ilrt, setup_linear_response, &
                              cleanup_linear_response

#if defined(MPI)
   use evb_parm, only: xch_type
#  if defined(LES)
   use evb_pimd, only: evb_pimd_init, PE_slice, master_worldrank, jobs_per_node
#  endif
! REMD
   use remd, only : rem, mdloop, repnum, remd1d_setup, remd_exchange, &
                    remd_cleanup, hremd_exchange, ph_remd_exchange, &
                    multid_remd_setup, multid_remd_exchange
   use bintraj, only: setup_remd_indices
#else
#  define rem 0
#endif /* MPI */

   use pimd_vars, only: ipimd
   use neb_vars, only: ineb

   use trajenemod, only: trajene

!RCW+MJW CHARMM SUPPORT
   use charmm_mod, only : charmm_active, charmm_deallocate_arrays, &
                          charmm_filter_out_qm_atoms
   use ff11_mod, only : cmap_active, deallocate_cmap_arrays

   use memory_module, only: x, ix, ih, memory_init

! Self-Guided molecular/Langevin Dynamics (SGLD)
   use sgld, only : isgld,psgld
   
   use nbips, only: ipssys,ips

   use crg_reloc, only: ifcr, cr_backup_charge, cr_cleanup, cr_allocate, &
                        cr_read_input, cr_check_input, cr_print_info

   use emap,only: temap,pemap,qemap

   use file_io_dat
   use constantph, only : cnstph_finalize
   use barostats, only : mcbar_setup

!AMD
   use amd_mod
!scaledMD
   use scaledMD_mod

   use abfqmmm_module  ! lam81

   implicit none

   logical belly, erstop
   integer ier,ifind,jn,ncalls,xmin_iter
   character(len=4) itest
   logical ok
   logical newstyle
#  include "../include/memory.h"
#  include "nmr.h"
#  include "box.h"
#  include "../include/md.h"
#  include "extra.h"
#  include "tgtmd.h"
#  include "multitmd.h"
#  include "les.h"

#  include "parallel.h"
#ifdef MPI
  ! =========================== AMBER/MPI ===========================
#  ifdef MPI_DOUBLE_PRECISION
#    undef MPI_DOUBLE_PRECISION
#  endif
   include 'mpif.h'
#  ifdef CRAY_PVP
#    define MPI_DOUBLE_PRECISION MPI_REAL8
#  endif
#  ifdef MPI_BUFFER_SIZE
   integer*4 mpibuf(mpi_buffer_size)
#  endif
!  REMD: loop is the current exchange. runmd is called numexchg times. 
   integer loop

   integer nrank, istat
   _REAL_ ener(30),vir(4)
   integer ierr
   integer partner
! ========================= END AMBER/MPI =========================
#endif
#  include "ew_pme_recip.h"
#  include "ew_frc.h"
#  include "ew_erfc_spline.h"
#  include "ew_parallel.h"
#  include "ew_mpole.h"
#  include "ew_cntrl.h"
#  include "def_time.h"

   type(state_rec) ::  ene
   integer native,nr3,nr

   ! nmrcal vars
   _REAL_ f,enmr,devdis,devang,devtor,devplpt,devpln,devgendis,ag,bg,cg
   ! Updated 9/2007 by Matthew Seetin to enable plane-point and plane-plane
   ! restraints
   _REAL_ emtmd
   integer numphi,nhb

   ! runmin/trajene var
   _REAL_ carrms

   ! dipole moment stuff
   integer ngrp

   character(len=8) initial_date, setup_end_date, final_date
   character(len=10) initial_time, setup_end_time, final_time
   integer nstlim_total ! for final time printout

   _REAL_ time0, time1

   integer idiff,i,j,istop,index,ierror,itemp

   integer, dimension(:), allocatable :: ipairs
   integer :: n_force_calls
   logical qsetup
   logical :: do_list_update=.false.
   _REAL_ :: box_center(3)   ! lam81
#ifdef MPI_DEBUGGER
   integer, volatile :: release_debug

   ! So only the master master thread has release_debug = 0
   release_debug = worldrank

   ! Lock us into an infinite loop while release_debug == 0 on any thread (only
   ! the master here). This allows you to connect a debugger to any running
   ! process without having to 'race' program execution.  A debugger MUST be
   ! attached to the master thread (typically the thread with the lowest PID),
   ! and have release_debug set to NOT 0 (e.g., via "set release_debug=1").
   do
      if (release_debug .ne. 0) exit
   end do

   ! Prevent any other threads from progressing past this point until all
   ! threads you want to watch with a debugger are watched and all those
   ! threads are continued.
   call mpi_barrier(mpi_comm_world, ierr)
#endif

   ! ---- HERE BEGIN THE EXECUTABLE STATEMENTS ----

   ! Initialize the cpu timer. Needed for machines where returned cpu times
   ! are relative.
   call date_and_time( initial_date, initial_time )
   call wallclock( time0 )
   call init_timers()

   ! Initialize the printing of ongoing time and performance summaries.
!   call print_ongoing_time_summary(0,0,0.0d0,7)
   call print_ongoing_time_summary(0,0,0.0d0,7,time0)

   ! BPR - original location of PUPIL interface. I moved it further down
   ! because, if it's here, it can't print stuff; write(6,...) statements
   ! assume mdread1() has already been invoked. However, moving this down
   ! may break other things.

   ! ==== Flag to tell list builder to print size of list on first call =======
   first_list_flag = .true.
   ! ==== Flag to tell recip space routines to allocate on first call =======
   first_pme = .true.

   ! ==== Initialise first_call flags for QMMM ====
   qmmm_struct%qm_mm_first_call = .true.
   qmmm_struct%fock_first_call = .true.
   qmmm_struct%fock2_2atm_first_call = .true.
   qmmm_struct%qm2_allocate_e_repul_first_call = .true.
   qmmm_struct%qm2_calc_rij_eqns_first_call = .true.
   qmmm_struct%qm2_scf_first_call = .true.
   qmmm_struct%zero_link_charges_first_call = .true.
   qmmm_struct%adj_mm_link_pair_crd_first_call = .true.
   qmmm_struct%num_qmmm_calls = 0

#ifdef MPI
   ! Parallel initialization (setup is done in multisander.F90).

   ! Make PE 0 the master
   master = mytaskid == 0
   master_master = masterrank == 0

   if ( master .and. numtasks > MPI_MAX_PROCESSORS ) then
      write(0, '(a,i4,a,i4)') &
            'Error: the number of processors must not be greater than ', &
            MPI_MAX_PROCESSORS, ', but is ', numtasks
      call mexit(6,1)
   end if
#  ifdef MPI_BUFFER_SIZE
   call mpi_buffer_attach(mpibuf, mpi_buffer_size*4, ierr)
#  endif
#else /* not MPI follows */
   ! In the single-threaded version, the one process is master
   master = .true.
#endif /* MPI */

   erstop = .false.
   qsetup = .true.

   ! --- generic packing scheme ---

   nwdvar = 1
   native = 32
#ifdef ISTAR2

   ! --- Int*2 packing scheme ---

   nwdvar = 2
#endif  /*ISTAR2*/
   numpk = nwdvar
   nbit = native/numpk

   ! ----- Only the master node (only node when single-process)
   !       performs the initial setup and reading/writing -----

   call timer_start(TIME_TOTAL)

   call abfqmmm_init_param()                                  ! lam81

   do while ( (abfqmmm_param%qmstep <= abfqmmm_param%maxqmstep) &
        .or. (abfqmmm_param%maxqmstep == 0 .and. abfqmmm_param%system == 2) ) ! lam81

   masterwork: if (master) then

   if (abfqmmm_param%abfqmmm == 0) then ! lam81

      ! ---- first, initial reads to determine memory sizes:
      call mdread1()
      call amopen(8,parm,'O','F','R')
      call rdparm1(8)
      if (mtmd /= 'mtmd' .or. itgtmd == 2) call mtmdlx(natom)
      ! --- now, we can allocate memory:

      call locmem()
      write(6,'(/,a,5x,a)') '|','Memory Use     Allocated'
      write(6,'(a,5x,a,i14)') '|', 'Real      ', lastr
      write(6,'(a,5x,a,i14)') '|', 'Hollerith ', lasth
      write(6,'(a,5x,a,i14)') '|', 'Integer   ', lasti
      write(6,'(a,5x,a,i14)') '|', 'Max Pairs ', lastpr

      ! --- dynamic memory allocation:

      ! GMS: 
      ! Allocate space for module molecule
      ! in the master node
      mol_info%natom = natom
      mol_info%nres  = nres
      call allocate_molecule()

      ! Allocate all global arrays
      allocate( x(lastr), ix(lasti), ipairs(lastpr), ih(lasth), stat = ier )
      REQUIRE( ier == 0 )
      ix(1:lasti) = 0

      ! This sets up pointer arrays in MEMORY_MODULE to match array-offsets into
      ! the shared X, IX, and IH arrays. Eventually, LOCMEM code should be
      ! merged with MEMORY_MODULE to allocate individual allocatable arrays, but
      ! that will also require updating the MPI code to handle individual
      ! arrays.
      call memory_init()

      ! Allocate the parm arrays
      call allocate_parms()

      if ((igb /= 0 .and. igb /= 10 .and. ipb == 0) &
                    .or.hybridgb>0.or.icnstph.gt.1) &
         call allocate_gb( natom, ncopy )

      if( idecomp > 0 ) then
#ifdef MPI
         if (ifsc > 0) then
            call synchronize_dec(natom, nres)
         else
            nat = natom
            nrs = nres
         end if
#else
         nat = natom
         nrs = nres
#endif
         call allocate_int_decomp(natom, nres)
      else
         call allocate_int_decomp(1, 1)
      end if

      write(6,'(a,5x,a,i14)'  ) '|', 'nblistReal', nblist_allreal
      write(6,'(a,5x,a,i14)'  ) '|', 'nblist Int', nblist_allint
      write(6,'(a,5x,a,i14,a)') '|', '  Total   ', &
            (8*(lastr+lastrst+nblist_allreal)  &
            + 4*(lasth+lasti+lastpr+lastist+nblist_allint))/1024, &
            ' kbytes'

      ! --- finish reading the prmtop file and other user input:
      call rdparm2(x,ix,ih,ipairs,8)

      call AMOEBA_readparm(8,ntf,ntc,natom,x(lmass))! ntf,ntc get reset if amoeba prmtop
#ifdef _XRAY
      call xray_read_parm(8,6)
#endif

   end if  ! lam81

      if (qmmm_nml%ifqnt .or. abfqmmm_param%abfqmmm == 1) then                                            ! lam81
       if(abfqmmm_param%abfqmmm == 0) then                                                                ! lam81
        call sebomd_setup
        call read_qmmm_nm_and_alloc(igb, ih, ix, x, cut, use_pme, ntb, 0)                                 ! lam81
        if (qmmm_nml%qmtheory%SEBOMD) then
            ! don't do QM/MM
            qmmm_nml%ifqnt= .false.
            sebomd_obj%do_sebomd = .true.
        end if
       end if                                                                                             ! lam81
       if(qmmm_struct%abfqmmm == 1 .and. abfqmmm_param%abfqmmm == 0) then                                 ! lam81
        call abfqmmm_setup(natom,nres,ix(i02),ih(m04),ih(m02),x(lmass), &                                 ! lam81
                           nbonh,nbona,ix(iibh),ix(ijbh),ix(iiba),ix(ijba))                               ! lam81
        nr=natom                                                                                          ! lam81
        call AMOEBA_check_newstyle_inpcrd(inpcrd,newstyle)                                                ! lam81
        if (newstyle) then                                                                                ! lam81
         call AM_RUNMD_get_coords(natom,t,irest,ntb,x(lcrd),x(lvel))                                      ! lam81
        else                                                                                              ! lam81
         call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.FALSE.)                          ! lam81
        end if                                                                                            ! lam81
        abfqmmm_param%maxqmstep = nstlim                                                                  ! lam81
       end if                                                                                             ! lam81
       if(abfqmmm_param%abfqmmm == 1) then                                                                ! lam81
        if(abfqmmm_param%system == 1) then                                                                ! lam81
         call abfqmmm_update_qmatoms(x(lcrd))                                                             ! lam81
         if(abfqmmm_param%ntwpdb < 0) then                                                                ! lam81
          call abfqmmm_write_pdb(x(lcrd),ix(i70))                                                         ! lam81
          close(6)                                                                                        ! lam81
          call mexit(6,1)                                                                                 ! lam81
         end if                                                                                           ! lam81
        end if                                                                                            ! lma81
        call abfqmmm_select_system_qmatoms(natom)                                                         ! lam81
        if(qmmm_nml%ifqnt) then                                                                           ! lam81
         call read_qmmm_nm_and_alloc(igb,ih,ix,x,cut,use_pme,ntb,abfqmmm_param%qmstep, &                  ! lam81
                                     abfqmmm_param%isqm,abfqmmm_param%abfcharge)                          ! lam81
        end if                                                                                            ! lam81
       endif                                                                                              ! lam81
      end if                                                                                              ! lam81

      if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then                                 ! lam81

      call mdread2(x,ix,ih,ipairs)

      endif                                                                                               ! lam81

#if defined(RISMSANDER) 
      call rism_setparam(mdin,&
           commsander,&
           natom,ntypes,x(L15:L15+natom-1),&
           x(LMASS:LMASS+natom-1),cn1,cn2,&
           ix(i04:i04+ntypes**2-1), ix(i06:i06+natom-1))
#endif /*RISMSANDER*/

      if ( ifcr /= 0 ) then
         call cr_read_input(natom)
         call cr_check_input( ips )
         !call cr_print_info(6)
         call cr_backup_charge( x(l15), natom )
      end if

      ! --- alloc memory for decomp module that needs info from mdread2
      if (idecomp == 1 .or. idecomp == 2) then
         call allocate_real_decomp(nrs)
#ifdef MPI
         ! -- ti decomp
         if (ifsc > 0) then
            !  following lines don't really seem to make sense(?)
            ! partner = ieor(masterrank,1)
            ! if (nat == natom) then
            !    nrank = masterrank
            ! else
            !    nrank = partner
            ! end if
            call mpi_bcast(jgroup, nat, MPI_INTEGER, 0, commmaster, ierr)
         end if
#endif
      else if( idecomp == 3 .or. idecomp == 4 ) then
         call allocate_real_decomp(npdec*npdec)
      end if

      ! ----- EVALUATE SOME CONSTANTS FROM MDREAD SETTINGS -----

      nr = nrp
      nr3 = 3*nr
      belly = ibelly > 0

      ! ========================= PUPIL INTERFACE =========================
#ifdef PUPIL_SUPPORT

      ! I moved the PUPIL interface down here so that write() statements work
      ! as advertised. BPR 9/7/09

      ! Initialise the CORBA interface
      puperror = 0
      call fixport()
      call inicorbaintfcmd(puperror)
      if (puperror .ne. 0) then
         write(6,*) 'Error creating PUPIL CORBA interface.'
         call mexit(6,1)
      end if
      pupactive = .true.
      write(6,*) 'PUPIL CORBA interface initialized.'

      ! Allocation of memory and initialization
      pupStep  = 0
      puperror = 0
      allocate (qcell   (12     ),stat=puperror)
      allocate (pupmask (natom  ),stat=puperror)
      allocate (pupqlist(natom  ),stat=puperror)
      allocate (pupatm  (natom  ),stat=puperror)
      allocate (pupchg  (natom  ),stat=puperror)
      allocate (qfpup   (natom*3),stat=puperror)
      allocate (qcdata  (natom*9),stat=puperror)
      allocate (keyMM   (natom  ),stat=puperror)
      allocate (pupres  (nres   ),stat=puperror)
      allocate (keyres  (nres   ),stat=puperror)

      if (puperror /= 0) then
         write(6,*) 'Error allocating PUPIL interface memory.'
         call mexit(6,1)
      end if

      ! Initialise the "atomic numbers" and "quantum forces" vectors        
      pupqatoms = 0
      iresPup   = 1
      pupres(1) = 1
      do iPup=1,natom
         bs1  = (iPup-1)*3
         call get_atomic_number_pupil(ih(iPup+m06-1),x(lmass+iPup-1),pupatm(iPup))
         if (iresPup .lt. nres) then
            if (iPup .ge. ix(iresPup+i02)) then
               iresPup = iresPup + 1
               pupres(iresPup) = iPup
            end if
         end if
         write (strAux,"(A4,'.',A4)") trim(ih(iresPup+m02-1)),adjustl(ih(iPup+m04-1))
         keyres(iresPup) = trim(ih(iresPup+m02-1))
         keyMM(iPup)     = trim(strAux)

         ! Retrieve the initial charges
         pupchg(iPup) = x(L15+iPup-1)
         !write(6,*) 'Atom num.',iPup,'Label,Mass,Atomic Num.', keyMM(iPup),x(lmass+iPup-1),pupatm(iPup), 'Charge', pupchg(iPup)

         do jPup=1,3
            qfpup(bs1+jPup) = 0.0d0
         end do
      end do

      write(6,*) 'Got all atomic numbers.'

      ! Initialise the PUPIL cell
      do iPup=1,12
         qcell(iPup) = 0.0d0
      end do

      ! Submit the KeyMM particles and their respective atomic numbers to PUPIL
      puperror = 0
      call putatomtypes(natom,puperror,pupatm,keyMM)
      if (puperror .ne. 0) then
         write(6,*) 'Error sending MM atom types to PUPIL.'
         call mexit(6,1)
      end if

      ! Submit the Residue Pointer vector to PUPIL
      write(6,"(a20,1x,i6,3x,a17,1x,i6)") 'Number of residues =', nres, 'Number of atoms =', natom
      !do iPup=1,nres
      !   write(6,*) 'Residue  ',iPup,keyres(iPup),pupres(iPup)
      !end do
      puperror = 0
      call putresiduetypes(nres,puperror,pupres,keyres)
      if (puperror .ne. 0) then
         write(6,*) 'Error sending MM residue types to PUPIL.'
         call mexit(6,1)
      end if

      write(6,*) 'Sent system data to PUPIL.'
      write(*,*) 'PUPIL structure initialized.'
#endif
      ! ========================= PUPIL INTERFACE =========================

      ! --- seed the random number generator ---

      ! DAN ROE: Note master node only here
      if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then   ! lam81
#ifdef MPI
      if (rem == 0) then
         call amrset(ig)
      else
         ! carlos: set random seed different for different replicas
         ! but keep same seed for cpus in the same replica since 
         ! we want data from diff # cpus to match
         call amrset(ig + (17 * nodeid)) ! nodeid is in md.h and is repnum - 1
      end if
#else
      call amrset(ig)
#endif

      if (nbit < 32 .and. nr > 32767) then
         write(6, *) '  Too many atoms for 16 bit pairlist -'
         write(6, *) '    Recompile without ISTAR2'
         call mexit(6, 1)
      end if

      if (ntp > 0.and.iabs(ntb) /= 2) then
         write(6,*) 'Input of NTP/NTB inconsistent'
         call mexit(6, 1)
      end if
      end if ! lam81

      ! ----- READ COORDINATES AND VELOCITIES -----

      if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then   ! lam81 <<< QMSTEP=1 BLOCK >>>

      call timer_start(TIME_RDCRD)
#ifdef LES
      call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0les,.TRUE.,solvph)
#else
      call AMOEBA_check_newstyle_inpcrd(inpcrd,newstyle)
      if (newstyle) then
         call AM_RUNMD_get_coords(natom,t,irest,ntb,x(lcrd),x(lvel))
      else
         call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.TRUE.,solvph)
      end if
#endif
      if (iamoeba > 0) then
         natom = natom*am_nbead
         nrp   = nrp*am_nbead
         nr    = nr*am_nbead
         nr3   = nr3*am_nbead
         ncopy = am_nbead
      end if

      ! M-WJ
      !if( igb == 0 .and. induced == 1 ) call get_dips(x,nr)
! WJM  if is a polarizable model, reading input dipole information
      if (igb == 0 .and. ipb == 0 .and. induced > 0) call get_dips(x,nr)

#ifdef APBS
      ! APBS initialization
      if (mdin_apbs) then
         ! in: natom, coords, charge and radii (from prmtop)
         ! out: pb charges and pb radii (via apbs_vars module)
         call apbs_init(natom, x(lcrd), x(l15), x(l97))
      end if
#endif /* APBS */

#ifdef _XRAY
      call xray_init()
#endif

      ! ----- SET THE INITIAL VELOCITIES -----

      if (ntx <= 3) then
         call setvel(nr,x(lvel),x(lwinv),tempi,init,iscale,scalm)
         ! random numbers may have been "used up" in setting the intial
         ! velocities; re-set the generator so that all nodes are back in
         ! sync

         ! DAN ROE: Note master node only here
#ifdef MPI
         if (rem == 0) call amrset(ig)
#else
         call amrset(ig)
#endif

      end if
      if (belly) call bellyf(natom,ix(ibellygp),x(lvel))
      call timer_stop(TIME_RDCRD)

      if(abfqmmm_param%abfqmmm == 1 .and. ntb > 0) then                                 ! lam81
       call iwrap2(abfqmmm_param%n_user_qm, abfqmmm_param%user_qm, x(lcrd), box_center) ! lam81
      end if                                                                            ! lam81

      ! --- If we are reading NMR restraints/weight changes,
      !     read them now:

      if (nmropt >= 1) then
         call nmrcal(x(lcrd),f,ih(m04),ih(m02),ix(i02),x(lwinv),enmr, &
               devdis,devang,devtor,devplpt,devpln,devgendis,temp0,tautp,cut,ntb,x(lnmr01), &
               ix(inmr02),x(l95),5,6,rk,tk,pk,cn1,cn2, &
               ag,bg,cg,numbnd,numang,numphi,nimprp, &
               nhb,natom,natom,ntypes,nres,rad,wel,radhb, &
               welhb,rwell,isftrp,tgtrmsd,temp0les,-1,'READ')
         ! Updated 9/2007 by Matthew Seetin to enable plane-point and plane-plane restraints

         ! --- Determine how many of the torsional parameters
         !     are impropers
         call impnum(ix(i46),ix(i56),ix(i48),ix(i58),nphih,nphia, &
               0,nptra,nimprp)
      end if

      ! -- Set up info related to weight changes for the non-bonds:

      call nmrrad(rad,wel,cn1,cn2,ntypes,0,0.0d0)
      call decnvh(asol,bsol,nphb,radhb,welhb)

      if (iredir(4) > 0) call noeread(x,ix,ih)
      if (iredir(8) > 0) call alignread(natom, x(lcrd))
      if (iredir(9) > 0) call csaread

      end if                  ! lam81 <<< QMSTEP=1 BLOCK >>>

      !---------------------------------------------------------------
      ! --- Call FASTWAT, which will tag those bonds which are part
      !     of 3-point water molecules. Constraints will be effected
      !     for these waters using a fast analytic routine -- dap.

      call timer_start(TIME_FASTWT)

      call fastwat(ih(m04),nres,ix(i02),ih(m02), &
            nbonh,nbona,ix(iibh),ix(ijbh),ibelly,ix(ibellygp), &
            iwtnm,iowtnm,ihwtnm,jfastw,ix(iifstwt), &
            ix(iifstwr),ibgwat,ienwat,ibgion,ienion,iorwat, &
            6,natom)
      call timer_stop(TIME_FASTWT)

      call getwds(ih(m04), nres        , ix(i02)     , ih(m02)     , &
             nbonh       , nbona       , 0           , ix(iibh)    , &
             ix(ijbh)    , iwtnm       , iowtnm      , ihwtnm      , &
             jfastw      , ix(iicbh)   , req         , x(lwinv)    , &
             rbtarg      , ibelly      , ix(ibellygp), 6               )

      ! Assign link atoms between quantum mechanical and molecular mechanical
      ! atoms if quantum atoms are present.
      ! After assigning the link atoms, delete all connectivity between the
      ! QM atoms.
      if(qmmm_nml%ifqnt) then

         call identify_link_atoms(nbona,ix(iiba),ix(ijba))

         ! Variable QM solvent:
         ! Store the original bond parameters since we will need to rebuild
         ! the QM region (delete bonded terms etc) repeatedly
         if ( qmmm_nml%vsolv > 0 ) then
            call new(qmmm_vsolv, nbonh, nbona, ntheth, ntheta, nphih, nphia)
            call qmmm_vsolv_store_parameters(qmmm_vsolv, numbnd, &
                 ix(iibh), ix(ijbh), ix(iicbh), &
                 ix(iiba), ix(ijba), ix(iicba), &
                 ix(i24), ix(i26), ix(i28), ix(i30), &
                 ix(i32), ix(i34), ix(i36), ix(i38), &
                 ix(i40), ix(i42), ix(i44), ix(i46), ix(i48), &
                 ix(i50), ix(i52), ix(i54), ix(i56), ix(i58))
         end if

         if( abfqmmm_param%abfqmmm == 1 ) then                                                    ! lam81
          if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then                      ! lam81
            call abfqmmm_allocate_arrays_of_parameters(numbnd, nbonh, nbona, ntheth, ntheta, nphih, nphia) ! lam81
            call abfqmmm_store_parameters(ix(iibh), ix(ijbh), ix(iicbh), &                        ! lam81
                                          ix(iiba), ix(ijba), ix(iicba), &                        ! lam81
                                          ix(i24), ix(i26), ix(i28), ix(i30), &                   ! lam81
                                          ix(i32), ix(i34), ix(i36), ix(i38), &                   ! lam81
                                          ix(i40), ix(i42), ix(i44), ix(i46), ix(i48), &          ! lam81
                                          ix(i50), ix(i52), ix(i54), ix(i56), ix(i58), &          ! lam81
                                          x(l15), rk, req)                                        ! lam81
          else                                                                                    ! lam81
            call abfqmmm_set_parameters(numbnd, nbonh, nbona, ntheth, ntheta, nphih, nphia, &     ! lam81
                                        ix(iibh), ix(ijbh), ix(iicbh), &                          ! lam81
                                        ix(iiba), ix(ijba), ix(iicba), &                          ! lam81
                                        ix(i24), ix(i26), ix(i28), ix(i30), &                     ! lam81
                                        ix(i32), ix(i34), ix(i36), ix(i38), &                     ! lam81
                                        ix(i40), ix(i42), ix(i44), ix(i46), ix(i48), &            ! lam81
                                        ix(i50), ix(i52), ix(i54), ix(i56), ix(i58), &            ! lam81
                                        x(l15), rk, req)                                          ! lam81

            call init_extra_pts(ix(iibh),ix(ijbh),ix(iicbh), &                                    ! lam81
                                ix(iiba),ix(ijba),ix(iicba), &                                    ! lam81
                                ix(i24),ix(i26),ix(i28),ix(i30), &                                ! lam81
                                ix(i32),ix(i34),ix(i36),ix(i38), &                                ! lam81
                                ix(i40),ix(i42),ix(i44),ix(i46),ix(i48), &                        ! lam81
                                ix(i50),ix(i52),ix(i54),ix(i56),ix(i58), &                        ! lam81
                                ih(m06),ix,x,ix(i08),ix(i10), &                                   ! lam81
                                nspm,ix(i70),x(l75),tmass,tmassinv,x(lmass),x(lwinv),req)         ! lam81

          end if                                                                                  ! lam81
         end if                                                                                   ! lam81

         ! Remove bonds between QM atoms from list (Hydrogen)
         if (nbonh .gt. 0) call setbon(nbonh,ix(iibh),ix(ijbh),ix(iicbh), &
                ix(ibellygp))

         ! Remove bonds between QM atoms from list (Heavy)
         if (nbona .gt. 0) call setbon(nbona,ix(iiba),ix(ijba),ix(iicba), &
                ix(ibellygp))

         ! Remove angles between QM atoms from list (Hydrogen)
         if (ntheth .gt. 0) call setang(ntheth,ix(i24),ix(i26),ix(i28),ix(i30), &
                ix(ibellygp))

         ! Remove angles between QM atoms from list (Heavy)
         if (ntheta .gt. 0) call setang(ntheta,ix(i32),ix(i34),ix(i36),ix(i38),&
                ix(ibellygp))

         ! Remove dihedrals between QM atoms from list (Hydrogen)
         if (nphih .gt. 0) call setdih(nphih,ix(i40),ix(i42),ix(i44),ix(i46), &
               ix(i48), ix(ibellygp))

         ! Remove dihedrals between QM atoms from list (Heavy)
         if (nphia .gt. 0) call setdih(nphia,ix(i50),ix(i52),ix(i54),ix(i56), &
               ix(i58), ix(ibellygp))

         ! Remove CHARMM energy terms from QM region
         call charmm_filter_out_qm_atoms()

         ! Now we should work out the type of each quantum atom present. 
         ! This is used for our arrays of pre-computed parameters. It is 
         ! essentially a re-basing of the atomic numbers and is done to save 
         ! memory. Note: qm_assign_atom_types will allocate the qm_atom_type 
         ! array for us. Only the master calls this routine. All other 
         ! threads get this allocated and broadcast to them by the mpi setup 
         ! routine.
         call qm_assign_atom_types

         ! Set default QMMM MPI parameters - for single cpu operation.
         ! These will get overwritten by qmmm_mpi_setup if MPI is on.
         qmmm_mpi%commqmmm_master = master
         qmmm_mpi%numthreads = 1
         qmmm_mpi%mytaskid = 0
         qmmm_mpi%natom_start = 1
         qmmm_mpi%natom_end = natom
         qmmm_mpi%nquant_nlink_start = 1
         qmmm_mpi%nquant_nlink_end = qmmm_struct%nquant_nlink

         ! Now we know how many link atoms we can allocate the scf_mchg array...
         allocate(qm2_struct%scf_mchg(qmmm_struct%nquant_nlink), stat = ier)
         REQUIRE(ier == 0) ! Deallocated in deallocate qmmm

         ! We can also allocate ewald_memory
         if (qmmm_nml%qm_ewald > 0 ) then
            call allocate_qmewald(natom)
         end if
         if (qmmm_nml%qmgb == 2 ) then                     ! lam81
         call allocate_qmgb(qmmm_struct%nquant_nlink)
         end if                                            ! lam81

         allocate( qmmm_struct%dxyzqm(3, qmmm_struct%nquant_nlink), stat = ier )
         REQUIRE(ier == 0) ! Deallocated in deallocate qmmm

      else if(abfqmmm_param%abfqmmm == 1) then             ! lam81
      
       call abfqmmm_set_parameters(numbnd, nbonh, nbona, ntheth, ntheta, nphih, nphia, &     ! lam81
                                   ix(iibh), ix(ijbh), ix(iicbh), &                          ! lam81
                                   ix(iiba), ix(ijba), ix(iicba), &                          ! lam81
                                   ix(i24), ix(i26), ix(i28), ix(i30), &                     ! lam81
                                   ix(i32), ix(i34), ix(i36), ix(i38), &                     ! lam81
                                   ix(i40), ix(i42), ix(i44), ix(i46), ix(i48), &            ! lam81
                                   ix(i50), ix(i52), ix(i54), ix(i56), ix(i58), &            ! lam81
                                   x(l15), rk, req)                                          ! lam81

       call init_extra_pts(ix(iibh),ix(ijbh),ix(iicbh), &                                    ! lam81
                           ix(iiba),ix(ijba),ix(iicba), &                                    ! lam81
                           ix(i24),ix(i26),ix(i28),ix(i30), &                                ! lam81
                           ix(i32),ix(i34),ix(i36),ix(i38), &                                ! lam81
                           ix(i40),ix(i42),ix(i44),ix(i46),ix(i48), &                        ! lam81
                           ix(i50),ix(i52),ix(i54),ix(i56),ix(i58), &                        ! lam81
                           ih(m06),ix,x,ix(i08),ix(i10), &                                   ! lam81
                           nspm,ix(i70),x(l75),tmass,tmassinv,x(lmass),x(lwinv),req)         ! lam81

      end if !if (qmmm_nml%ifqnt)

      ! --- Open the data dumping files and position it depending
      !     on the type of run:

      if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81
#ifdef MPI
      ! adaptive QM/MM (qmmm_nml%vsolv=2) via multisander
      ! all groups have identical coords and velocities
      ! only master of first group needs to dump results
      if ( (qmmm_nml%vsolv < 2) .or. (worldrank == 0) ) then
#endif
         call open_dump_files
#ifdef MPI
      end if
#endif
      end if                                                              ! lam81
      if (master) call amflsh(6)

   ! --- end of master process setup ---
   end if masterwork ! (master)

#  if defined(RISMSANDER)
   call rism_init(commsander)
#  endif /* RISMSANDER */

#ifdef MPI

  call mpi_barrier(commsander,ierr)
   !     =========================== AMBER/MPI ===========================

   !     NOTE: in the current AMBER/MPI implementation, two means of
   !     running in parallel within sander are supported. The value
   !     of mpi_orig determines which approach is used.
   !     This is turned on when minimization (imin .ne. 0) is requested,
   !     and is otherwise off.

   !     When running the mpi_orig case, a variable notdone is now
   !     set by the master and determines when to exit the force()
   !     loop.  When the master has finished calling force, the
   !     master changes notdone to 0 and broadcasts the data one more
   !     time to signal end of the loop.  force() is modified so that
   !     in the mpi_orig case, an initial broadcast is done to receive
   !     the value from the master to decide whether to do the work or
   !     simply exit.

   !     ...set up initial data and send all needed data to other nodes,
   !     now that the master has it

   !     First, broadcast parameters in memory.h, so that all processors
   !     will know how much memory to allocate:

   call mpi_bcast(natom,BC_MEMORY,mpi_integer,0,commsander,ierr)
   ! -- ti decomp
   call mpi_bcast(idecomp,1,mpi_integer,0,commsander,ierr)
   call mpi_bcast(nat,1,mpi_integer,0,commsander,ierr)
   call mpi_bcast(nrs,1,mpi_integer,0,commsander,ierr)

   ! ----- Set up integer stack initial size --------------    
   call mpi_bcast(lastist,1,mpi_integer,0,commsander,ierr)
   call mpi_bcast(lastrst,1,mpi_integer,0,commsander,ierr)

   call mpi_bcast(qmmm_struct%abfqmmm,1,mpi_integer,0,commsander,ierr) ! lam81
   call mpi_bcast(abfqmmm_param%abfqmmm,1,mpi_integer,0,commsander,ierr) ! lam81
   if(abfqmmm_param%abfqmmm == 1) then   ! lam81
    call mpi_bcast(abfqmmm_param%natom,1,mpi_integer,0,commsander,ierr) ! lam81
    call mpi_bcast(abfqmmm_param%qmstep,1,mpi_integer,0,commsander,ierr) ! lam81
    call mpi_bcast(abfqmmm_param%maxqmstep,1,mpi_integer,0,commsander,ierr) ! lam81
    call mpi_bcast(abfqmmm_param%system,1,mpi_integer,0,commsander,ierr) ! lam81
   end if       ! lam81

   if(abfqmmm_param%abfqmmm == 1) then      ! lam81
    if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then  ! lam81
     call mpi_bcast(abfqmmm_param%r_core_in,1,mpi_double_precision,0,commsander,ierr) ! lam81
     call mpi_bcast(abfqmmm_param%r_core_out,1,mpi_double_precision,0,commsander,ierr) ! lam81
     call mpi_bcast(abfqmmm_param%r_qm_in,1,mpi_double_precision,0,commsander,ierr) ! lam81
     call mpi_bcast(abfqmmm_param%r_qm_out,1,mpi_double_precision,0,commsander,ierr) ! lam81
     call mpi_bcast(abfqmmm_param%r_buffer_in,1,mpi_double_precision,0,commsander,ierr) ! lam81
     call mpi_bcast(abfqmmm_param%r_buffer_out,1,mpi_double_precision,0,commsander,ierr) ! lam81
     call mpi_bcast(abfqmmm_param%mom_cons_type,1,mpi_integer,0,commsander,ierr) ! lam81
     call mpi_bcast(abfqmmm_param%mom_cons_region,1,mpi_integer,0,commsander,ierr) ! lam81
     if (.not. master) then     ! lam81
      allocate(abfqmmm_param%x(3*natom), stat=ier)   ! lam81
      REQUIRE(ier==0)          ! lam81
      allocate(abfqmmm_param%id(natom), stat=ier)   ! lam81
      REQUIRE(ier==0)          ! lam81
      allocate(abfqmmm_param%mass(natom), stat=ier) ! lam81
      REQUIRE(ier==0)          ! lam81
     end if                   ! lam81
     allocate(abfqmmm_param%v(3*natom+iscale), stat=ier)    ! lam81
     REQUIRE(ier==0)          ! lam81
     allocate(abfqmmm_param%f(3*natom+iscale), stat=ier)    ! lam81
     REQUIRE(ier==0)          ! lam81
     allocate(abfqmmm_param%f1(3*natom+iscale), stat=ier)   ! lam81
     REQUIRE(ier==0)          ! lam81
     allocate(abfqmmm_param%f2(3*natom+iscale), stat=ier)   ! lam81
     REQUIRE(ier==0)          ! lam81
    end if                    ! lam81
    call mpi_bcast(abfqmmm_param%x,3*natom,mpi_double_precision,0,commsander,ierr) ! lam81
    call mpi_bcast(abfqmmm_param%id,natom,mpi_integer,0,commsander,ierr) ! lam81
    call mpi_bcast(abfqmmm_param%mass,natom,mpi_double_precision,0,commsander,ierr) ! lam81
   end if     ! lam81

   if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then   ! lam81
   call stack_setup()
   else                    ! lam81
   call deallocate_stacks  ! lam81             
   call stack_setup()      ! lam81
   end if                  ! lam81

   call mpi_bcast(plumed,1,MPI_INTEGER,0,commsander,ierr)
   call mpi_bcast(plumedfile,256,MPI_CHARACTER,0,commsander,ierr)

   ! GMS: Broadcast parameters from module 'molecule'
   call mpi_bcast(mol_info%natom,1,mpi_integer,0,commsander,ierr)
   call mpi_bcast(mol_info%nres,1,mpi_integer,0,commsander,ierr)
   call mpi_barrier(commsander,ierr)

   ! ---allocate memory on the non-master nodes:
   if (.not. master) then

      if (abfqmmm_param%qmstep /= 1 .or. abfqmmm_param%system /= 1) then ! lam81
       deallocate(x,ix,ipairs,ih)                                        ! lam81
      end if                                                             ! lam81

      allocate( x(1:lastr), stat = ier )
      REQUIRE( ier == 0 )

      allocate( ix(1:lasti), stat = ier )
      REQUIRE( ier == 0 )

      allocate( ipairs(1:lastpr), stat = ier )
      REQUIRE( ier == 0 )

      allocate( ih(1:lasth), stat = ier )
      REQUIRE( ier == 0 )

      ! GMS: 
      ! Allocate space for molecule module 
      ! arrays in the other nodes
      if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81
      call allocate_molecule()
      else                                                                ! lam81
       call deallocate_molecule()                                         ! lam81
       call allocate_molecule()                                           ! lam81
      end if                                                              ! lam81

      ! --------------------------

      ! -- ti decomp
      if (idecomp > 0) then
         call allocate_int_decomp(natom, nres)
         if (idecomp == 1 .or. idecomp == 2) then
            call allocate_real_decomp(nrs)
         else if( idecomp == 3 .or. idecomp == 4 ) then
            call allocate_real_decomp(npdec*npdec)
         end if
      end if
   end if  ! ( .not.master )

   !GMS: Broadcast arrays from module 'molecule'
   call mpi_bcast(mol_info%natom_res,mol_info%nres,MPI_INTEGER,0,commsander,ierr)
   call mpi_bcast(mol_info%atom_to_resid_map,mol_info%natom,MPI_INTEGER,0,commsander,ierr)
   call mpi_bcast(mol_info%atom_mass,mol_info%natom, MPI_DOUBLE_PRECISION, 0, commsander, ier)

   if(idecomp == 1 .or. idecomp == 2) then
      call mpi_bcast(jgroup,nat,MPI_INTEGER,0,commsander,ierr)
   end if

   if(icfe == 0 .and. (idecomp ==3 .or. idecomp == 4)) then
      call mpi_bcast(jgroup,natom,MPI_INTEGER,0,commsander,ierr)
      call mpi_bcast(indx,nres,MPI_INTEGER,0,commsander,ierr)
   end if

   call startup_groups(ierr)
   call startup(x,ix,ih)

!  +---------------------------------------------------------------------------+
!  |  Broadcast EVB/PIMD inputs/parameters to all PEs                          |
!  +:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::+
!  |  Note: The masters have all required EVB/PIMD inputs/parameters via call  |
!  |  to mdread2 ( evb_input, evb_init, evb_pimd_init ).  For EVB/PIMD, all    |
!  |  PEs need the inputs/parameters ... so we perform this initialization     |
!  |  again for all PEs besides the masters.  The alternative is to use        |
!  |  MPI_BCAST.                                                               |
!  +---------------------------------------------------------------------------+

   call mpi_bcast( ievb , 1, MPI_INTEGER, 0, commworld, ierr )
   call mpi_bcast( ipimd, 1, MPI_INTEGER, 0, commworld, ierr )
  
   ! KFW
   if (ievb /= 0) then
      call mpi_bcast(evbin, MAX_FN_LEN, MPI_CHARACTER, 0, commworld, ierr)
      if (.not. master) then
         call evb_input
         call evb_init
      end if
   end if
  ! KFW


#  if defined(LES)
   call mpi_bcast ( ncopy, 1, MPI_INTEGER, 0, commworld, ierr )
   call mpi_bcast ( cnum(1:natom), natom, MPI_INTEGER, 0, commworld, ierr )
   call mpi_bcast ( evbin, MAX_FN_LEN, MPI_CHARACTER, 0, commworld, ierr )
#  endif /* LES */
   if (ievb /= 0) then
#  if defined(LES)
      !call mpi_bcast ( mastersize, 1, MPI_INTEGER, 0, commworld, ierr )
      !call mpi_bcast ( jobs_per_node, 1, MPI_INTEGER, 0, commworld, ierr )
      !call mpi_bcast ( nsize, 1, MPI_INTEGER, 0, commworld, ierr )
      if (ipimd > 0 .and. .not. master) then
         call evb_input
         call evb_init
         !call evb_pimd_init
      end if
      call evb_pimd_init
      !call mpi_bcast ( master_worldrank, mastersize, MPI_INTEGER, 0, commworld, ierr )
      !call mpi_bcast ( PE_slice, jobs_per_node*nsize, MPI_INTEGER, 0, commworld, ierr )
#  endif /* LES */
      call evb_bcast
      call evb_alloc
   end if

!  +---------------------------------------------------------------------------+
!  |  Obtain B vector for Schlegel's distributed Gaussian method               |
!  +---------------------------------------------------------------------------+

   if (trim(adjustl(xch_type)) == "dist_gauss") call schlegel_dg

   if (ifsc /= 0) then
      ! multi-CPU minimization does not work with soft core !
      if (imin > 0 .and. numtasks > 1) then
         call sander_bomb('imin > 0 and numtasks > 1', &
                     'TI minimizations cannot be performed with > 2 CPUs','')
      end if
      call setup_sc(natom, nres, ih(m04), ih(m06), &
          ix(i02), ih(m02), x(lcrd), ntypes, clambda, nstlim)
      if (ntp > 0 .and. master) then
         ! Check which molecules are perturbed in NPT runs 
         call sc_check_perturbed_molecules(nspm, ix(i70))
      end if
      ! -- ti decomp
      if (idecomp > 0) then
         if (sanderrank == 0) call build_dec_mask(natom)
         call mpi_bcast(decmask, natom, MPI_INTEGER, 0, commsander, ierr)
      end if
      ! Make sure all common atoms have the same v (that of V0) in TI runs
      if ( ifsc /= 2 ) then
         if (master) call sc_sync_x(x(lvel),nr3)
         !call mpi_barrier(commsander,ierr)
         if (numtasks > 1) then
            call mpi_bcast(nr3,1,MPI_INTEGER,0,commsander,ierr)
            call mpi_bcast(x(lvel),nr3,MPI_DOUBLE_PRECISION,0,commsander,ierr)               
         end if
      end if
      if (tishake /= 0) call setnoshake_sc(ix,ntc,num_noshake,master)
   else
      extra_atoms=0      
   end if
   if (ifmbar /= 0) then
      call setup_mbar(nstlim)
   end if

   if (.not. master .and. igb == 0 .and. ipb == 0) then
        if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81
         call nblist_allocate(natom,ntypes,num_direct,numtasks)
        else                                                     ! lam81
         call nblist_deallocate                                  ! lam81
         call nblist_allocate(natom,ntypes,num_direct,numtasks)  ! lam81
        end if                                                   ! lam81
   end if

   ! --- Allocate memory for GB on the non-master nodes:

   if( .not.master ) then

      if ((igb /= 0 .and. igb /= 10 .and. ipb == 0) .or. hybridgb>0 .or. icnstph.gt.1) then ! lam81
         if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81
          call allocate_gb( natom, ncopy )
         else                                                ! lam81
          call deallocate_gb                                 ! lam81
          call allocate_gb( natom, ncopy )                   ! lam81
         end if                                              ! lam81
      end if                                                 ! lam81

   end if  ! ( .not.master )

   nr = nrp
   nr3 = 3*nr
   belly = ibelly > 0

   !  ---Do setup for QMMM in parallel if ifqnt is on - this is essentially 
   !     everything in the QMMM namelist and a few other bits and pieces.
   !  ---Note, currently only master node knows if qmmm_nml%ifqnt is 
   !     on so we need to broadcast this first and then make decisions based 
   !     on this.

   call mpi_bcast(qmmm_nml%ifqnt, 1, mpi_logical, 0, commsander, ierror)

   if (qmmm_nml%ifqnt) then
      ! Broadcast all of the stuff in qmmm_nml and allocate the relevant 
      ! arrays on all processors. This will also set up information for
      ! openmp on the master processor if it is in use.
      call qmmm_mpi_setup( master, natom )
      if (qmmm_nml%qm_ewald > 0 .and. .not. master) then
         call allocate_qmewald(natom)
      end if
      if (qmmm_nml%qmgb==2 .and. .not. master) then
         call allocate_qmgb(qmmm_struct%nquant_nlink)
      end if

      if (.not. master) then
        allocate( qmmm_struct%dxyzqm(3, qmmm_struct%nquant_nlink), stat = ier )
        REQUIRE(ier == 0) !Deallocated in deallocate qmmm
      end if

   end if

   !    --- END QMMM MPI SETUP ---

   ! DAN ROE: Note: all nodes are calling this. amrset(ig) has already been called
   !          by the master node (twice if initial velocities are set, ntx<=3).
   ! DAN ROE: REMD: Now need to call amrset on all child threads, masters have called
   !          it above before initial coord read 
   if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then   ! lam81
   if(rem == 0) then
      call amrset(ig+1) 
   else
      if (.not. master) call amrset(ig + 17 * nodeid)
   endif

   if (nmropt >= 1) &
         call nmrcal(x(lcrd),f,ih(m04),ih(m02),ix(i02),x(lwinv),enmr, &
         devdis,devang,devtor,devplpt,devpln,devgendis,temp0,tautp,cut,ntb,x(lnmr01), &
         ix(inmr02),x(l95),5,6,rk,tk,pk,cn1,cn2, &
         ag,bg,cg,numbnd,numang,numphi,nimprp, &
         nhb,natom,natom,ntypes,nres,rad,wel,radhb, &
         welhb,rwell,isftrp,tgtrmsd,temp0les,-1,'MPI ')
   end if ! lam81
   ! Updated 9/2007 by Matthew Seetin to enable plane-point and plane-plane restraints

   call mpi_bcast(lmtmd01,1,mpi_integer,0, commsander,ierr)
   call mpi_bcast(imtmd02,1,mpi_integer,0, commsander,ierr)

   if (itgtmd == 2) call mtmdcall(emtmd,x(lmtmd01),ix(imtmd02),x(lcrd),x(lforce),&
                            ih(m04),ih(m02),ix(i02),ih(m06),x(lmass),natom,nres,'MPI ')


   ! ---------------- Check system is neutral and print warning message ------
   ! ---------------- adjust charges for roundoff error.                ------
   if( igb == 0 .and. ipb == 0 .and. iyammp == 0 ) then
     if ( icfe == 0 ) then
       call check_neutral(x(l15),natom)
     else
       call ti_check_neutral(x(l15),natom)
     end if
   end if

   ! tell all nodes if this is a SEBOMD run
   call mpi_bcast(sebomd_obj%do_sebomd, 1, MPI_LOGICAL, 0, commsander, ierr)
   if (sebomd_obj%do_sebomd) then
     ! transfer SEBOMD info to the nodes
     call sebomd_bcast_obj()
     ! open necessary files
     if (master) then
        call sebomd_open_files
     end if
     ! initialize SEBOMD arrays
     call init_sebomd_arrays(natom)
   end if

   ! ---------------- Old parallel for minimization ----------------------

   if (imin /= 0) then
      mpi_orig = .true.
      notdone = 1
   else
      mpi_orig = .false.
   end if

   if (mpi_orig .and. .not. master) then

      ! ...all nodes only do the force calculations (JV)
      ! Minimisation so only only master gets past the loop below
      ! hence need to zero QM charges on non-master threads here.

      if (qmmm_nml%ifqnt) then
         ! Apply charge correction if required.
         if (qmmm_nml%adjust_q>0) then
            call qmmm_adjust_q(qmmm_nml%adjust_q, natom, qmmm_struct%nquant, qmmm_struct%nquant_nlink, &
                  qmmm_struct%nlink, x(L15), &
                  qmmm_struct%iqmatoms, qmmm_nml%qmcharge, qmmm_struct%atom_mask, &
                  qmmm_struct%mm_link_mask, master,x(LCRD), qmmm_nml%vsolv)
         end if
         ! At this point we can also fill the qmmm_struct%scaled_mm_charges
         ! array - we only need to do this once as the charges are constant
         ! during a run. Having a separate array of scaled charges saves us
         ! having to do it on every qmmm routine call. Do this BEFORE zeroing
         ! the QM charges since that routine take care of these values as well.
         do i = 1, natom
            qmmm_struct%scaled_mm_charges(i) = x(L15+(i-1)) * INV_AMBER_ELECTROSTATIC &
             * qmmm_nml%chg_lambda ! charge scaling factor for FEP
         end do
         ! Zero out the charges on the quantum mechanical atoms
         call qm_zero_charges(x(L15),qmmm_struct%scaled_mm_charges,.true.)
         if (qmmm_struct%nlink > 0 ) then
            ! We need to exclude all electrostatic
            ! interactions with MM link pairs, both QM-MM and MM-MM. Do this by
            ! zeroing the MM link pair charges in the main charge array.
            ! These charges are stored in qmmm_struct%mm_link_pair_resp_charges in case
            ! they are later needed.
            call qm_zero_mm_link_pair_main_chg(qmmm_struct%nlink,qmmm_struct%link_pairs,x(L15), &
                  qmmm_struct%scaled_mm_charges,.true.)
         end if

      end if

      if (igb == 7 .or. igb == 8) call igb7_init(natom, ncopy, x(l97)) !x(l97) is rborn()
      !Hai Nguyen: add igb == 8 here

      if ( ifcr /= 0 ) call cr_allocate( master, natom )

      n_force_calls = 0
      do while( notdone == 1 )
         n_force_calls = n_force_calls+1
         call force(x,ix,ih,ipairs,x(lcrd),x(lforce),ene,vir, &
                x(l96), x(l97), x(l98), x(l99), qsetup, &
                do_list_update,n_force_calls)
      end do

      ! Deallocate and return
      goto 999

   end if

   ! ----------------------------------------------------------------------

   if (master) then
      if(abfqmmm_param%abfqmmm /=1) write(6, '(a,i4,a,/)') '|  Running AMBER/MPI version on ',numtasks, ' nodes' ! lam81
      !BTREE is selected by default if cpu is a power of two.
      !The number of processes is required to be a power of two for Btree
      !Print a warning about inefficiency with cpus not being a power of 2.
      if (numtasks > 1 .and. logtwo(numtasks) <= 0) then 
       if(abfqmmm_param%abfqmmm /=1) write(6, '(a,i4,a,/)') '|  WARNING: The number of processors is not a power of 2' ! lam81
       if(abfqmmm_param%abfqmmm /=1) write(6, '(a,i4,a,/)') '|           this may be inefficient on some systems.' ! lam81
      end if
   end if
   if (master .and. numgroup > 1) write(6, '(a,i4,a,i4,a,i4,a)') &
         '|  MULTISANDER: ', numgroup, ' groups. ', &
         numtasks, ' processors out of ', worldsize, ' total.'
   if (master) call amflsh(6)

   ! ========================= END AMBER/MPI =========================
#else
   if(abfqmmm_param%abfqmmm == 1) then      ! lam81
    if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then  ! lam81
     allocate(abfqmmm_param%v(3*natom+iscale), stat=ier)    ! lam81
     REQUIRE(ier==0)          ! lam81
     allocate(abfqmmm_param%f(3*natom+iscale), stat=ier)    ! lam81
     REQUIRE(ier==0)          ! lam81
     allocate(abfqmmm_param%f1(3*natom+iscale), stat=ier)   ! lam81
     REQUIRE(ier==0)          ! lam81
     allocate(abfqmmm_param%f2(3*natom+iscale), stat=ier)   ! lam81
     REQUIRE(ier==0)          ! lam81
    end if                    ! lam81
   end if     ! lam81

   !   debug needs to copy charges at start and they can't change later
   !   ---------------- Check system is neutral and print warning message ------
   !   ---------------- adjust charges for roundoff error.                ------
   if( igb == 0 .and. ipb == 0 .and. iyammp == 0 ) call check_neutral(x(l15),natom)

   if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then   ! lam81
   call amrset(ig+1)
   end if ! lam81

   if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then   ! lam81
   call stack_setup()
   else                    ! lam81
   call deallocate_stacks  ! lam81             
   call stack_setup()      ! lam81
   end if                  ! lam81

   if (sebomd_obj%do_sebomd) then
     ! open necessary files
     call sebomd_open_files
     ! initialize SEBOMD arrays
     call init_sebomd_arrays(natom)
   endif

#endif /* MPI */

#ifdef OPENMP

   ! If -openmp was specified to configure_amber then -DOPENMP is defined and the 
   ! threaded version of MKL will have been linked in. It is important here that
   ! we set the default number of openmp threads for MKL to be 1 to stop conflicts
   ! with threaded vectorization routines when running in parallel etc.
   ! Individual calls to MKL from routines that know what they are doing - e.g.
   ! QMMM calls to diagonalizers etc can increase this limit as long as they
   ! put it back afterwards.
   call omp_set_num_threads(1)

   ! If we are using openmp for matrix diagonalization print some information.
   if (qmmm_nml%ifqnt .and. master) call qm_print_omp_info()
#endif

   ! allocate memory for crg relocation
   if (ifcr /= 0) call cr_allocate( master, natom )

   ! initialize LIE module if used
   if ( ilrt /= 0 ) then
      call setup_linear_response(natom,nres,ih(m04),ih(m06),ix(i02),ih(m02),x(lcrd),x(l15), &
                                 ntypes, ix(i04), ix(i06), cn1, cn2, master)
   end if

   call date_and_time( setup_end_date, setup_end_time )

   ! Initialize the printing of ongoing time and performance summaries.
   ! We do this quite late here to avoid including all the setup time.
   if (master) call print_ongoing_time_summary(0,0,0.0d0,7)

   ! ----------------------------------------------------------------------
   ! Now do the dynamics or minimization.
   ! ----------------------------------------------------------------------

   if (igb == 7 .or. igb == 8 ) call igb7_init(natom, ncopy, x(l97)) !x(l97) is rborn()
     !Hai Nguyen: add igb ==8 here

   if (qmmm_nml%ifqnt) then
      ! Apply charge correction if required.
      if (qmmm_nml%adjust_q>0) then
         call qmmm_adjust_q(qmmm_nml%adjust_q, natom, qmmm_struct%nquant, qmmm_struct%nquant_nlink, &
               qmmm_struct%nlink, x(L15), &
               qmmm_struct%iqmatoms, qmmm_nml%qmcharge, qmmm_struct%atom_mask, &
               qmmm_struct%mm_link_mask, master,x(LCRD), qmmm_nml%vsolv)
      end if
      ! At this point we can also fill the qmmm_struct%scaled_mm_charges
      ! array - we only need to do this once as the charges are constant
      ! during a run. Having a separate array of scaled charges saves us
      ! having to do it on every qmmm routine call. Do this BEFORE zeroing
      ! the QM charges since that routine take care of these values as well.
      do i = 1, natom
         qmmm_struct%scaled_mm_charges(i) = x(L15+(i-1)) * INV_AMBER_ELECTROSTATIC &
                                           * qmmm_nml%chg_lambda ! charge scaling factor for FEP
      end do

      ! Zeroing of QM charges MUST be done AFTER call to check_neutral.
      ! Zero out the charges on the quantum mechanical atoms.
      call qm_zero_charges(x(L15),qmmm_struct%scaled_mm_charges,.true.)

      if (qmmm_struct%nlink > 0) then
         ! We need to exclude all electrostatic
         ! interactions with MM link pairs, both QM-MM and MM-MM. Do this by
         ! zeroing the MM link pair charges in the main charge array.
         ! These charges are stored in qmmm_struct%mm_link_pair_resp_charges in case
         ! they are later needed.
         call qm_zero_mm_link_pair_main_chg(qmmm_struct%nlink,qmmm_struct%link_pairs,x(L15), &
                                            qmmm_struct%scaled_mm_charges,.true.)
      end if

   end if

   ! Use the debugf namelist to activate
   call debug_frc(x,ix,ih,ipairs,x(lcrd),x(lforce),cn1,cn2,qsetup)

  ! Prepare for SGLD simulation
   if (isgld > 0) call psgld(natom,x(lmass),x(lvel), rem)

  ! Prepare for EMAP constraints
   if (temap) call pemap(dt,temp0,x,ix,ih)

      ! Prepare for Isotropic periodic sum of nonbonded interaction
   if (ips .gt. 0) call ipssys(natom,ntypes,ntb,x(l15), &
            cut,cn1,cn2,ix(i04),ix(i06),x(lcrd))

   if (master .and. (.not. qmmm_nml%ifqnt) .and. (abfqmmm_param%abfqmmm /= 1)) &   ! lam81
         write(6,'(/80(''-'')/,''   4.  RESULTS'',/80(''-'')/)')

   ! Set up the MC barostat if requested
   if (ntp > 0 .and. barostat == 2) call mcbar_setup(ig)

   ! Input flag imin determines the type of calculation: MD, minimization, ...
   select case (imin)
   case (0)
      ! --- Dynamics:

      call timer_start(TIME_RUNMD)
      ! Set up AMD 
      if (iamd .gt. 0) then
        call amd_setup(ntwx)
      endif

      ! Set up scaledMD 
      if (scaledMD .gt. 0) then
        call scaledMD_setup(ntwx)
      endif

      if (ipimd > 0) then
         call pimd_init(natom,x(lmass),x(lwinv),x(lvel),ipimd)
      end if

      if(ineb>0) call neb_init()

#ifdef MPI
      ! ----===== REMD =====----
      ! If this is not a REMD run, runmd is called only once. 
      ! If this is a REMD run, runmd is called 0 to numexchg times,
      ! where the 0th runmd is just for getting initial PEs (no dynamics).
      if (rem == 0) then
         ! Not a REMD run. runmd will be called once.
         loop = 0
      else if (rem > 0) then
         ! This is a REMD run. runmd will be called numexchg times.
         loop = numexchg

         ! Set up temptable, open remlog, etc.
         call remd1d_setup(numexchg, hybridgb, numwatkeep, &
                           temp0, mxvar, natom, ig, solvph)
      else
         ! Multi-D REMD run
         loop = numexchg
         call multid_remd_setup(numexchg, hybridgb, numwatkeep, &
                                temp0, mxvar, natom, ig, solvph, irest)
         ! Now set up REMD indices for traj/restart writes. Only do this on
         ! master since only master handles writes.
         if (master) call setup_remd_indices 
      end if ! Replica run setup

      ! Loop over REMD exchanges
      do mdloop = 0, loop

         ! ----===== REMD EXCHANGE HANDLING =====----
         ! Note: mdloop==0 is just used to get initial energies for the
         !        first exchange. 
         if (rem < 0 .and. mdloop > 0) then
            call multid_remd_exchange(x, ix, ih, ipairs, qsetup, &
                                      do_list_update, temp0, solvph)
         else if ((rem == 1 .or. rem == 2) .and. mdloop > 0) then
            call remd_exchange(1, rem, x(lcrd), x(lvel), x(lmass), &
                               nr3, natom, nr, temp0)
         else if (rem == 3 .and. mdloop > 0) then
            call hremd_exchange(1, x, ix, ih, ipairs, qsetup, do_list_update)
            ! force was called inside hremd_exchange, so call nmrdcp
            ! to decrement the NMR counter, since this should not count
            ! as a real step. This is OK, since the counter got
            ! incremented at the _very_ end of nmrcal, so we haven't already
            ! printed an unwanted value (JMS 2/12)
            if (nmropt /= 0) call nmrdcp
         else if (rem == 4 .and. mdloop > 0) then
            call ph_remd_exchange(1, solvph)
         end if ! rem>0 and mdloop>0 

         ! ----===== END REMD EXCHANGE HANDLING =====----

#ifdef VERBOSE_REMD
         if (rem > 0 .and. mdloop .eq. 0 .and. master) &
            write (6,'(a,i4)') "REMD: Getting initial energy for replica ",repnum
#endif /* VERBOSE_REMD */
#endif  /* MPI */

#ifndef DISABLE_NCSU
         if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then  ! lam81
            call ncsu_on_sander_init(ih, x(lmass), x(lcrd), rem)
         end if                ! lam81
#endif /* DISABLE_NCSU */

         if (beeman_integrator == 1) then
            call AM_RUNMD(ix,ih,ipairs, &
                  x(lwinv),x(lmass),x, &
                  x(lcrd),x(lvel),x(lforce),qsetup)
         else
            call runmd(x,ix,ih,ipairs, &
                  x(lcrd),x(lwinv),x(lmass),x(lforce), &
                  x(lvel),x(lvel2),x(l45),x(lcrdr), &
                  x(l50),x(l95),ix(i70),x(l75), &
                  erstop,qsetup)
         end if !beeman_integrator == 1
   
#ifndef DISABLE_NCSU
         if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then  ! lam81
         call ncsu_on_sander_exit()
         end if                ! lam81
#endif /* DISABLE_NCSU */

    ! ----===== END REMD =====----
#ifdef MPI
      end do ! Loop over REMD exchanges (mdloop)

      ! Cleanup REMD files.
      if (rem /= 0) call remd_cleanup()

#endif

      call timer_stop(TIME_RUNMD)

      if (master) call amflsh(6)

      if (erstop) then
         ! This error condition stems from subroutine shake;
         ! furthermore, it seems that erstop can never be true since shake
         ! can never return with its third last argument, niter, equal to 0.
         ! SRB, Sep 24, 2003
         if (master) then
            write(6, *) 'FATAL ERROR'
         end if
         call mexit(6,1)
      end if

   case (1)

      !--- Minimization:

      ! Input flag ntmin determines the method of minimization
      select case (ntmin)
      case (0, 1, 2)
         call runmin(x,ix,ih,ipairs,x(lcrd),x(lforce),x(lvel), &
               ix(iibh),ix(ijbh),x(l50),x(lwinv),ix(ibellygp), &
               x(l95),ene, carrms, qsetup)
         ! If a conventional minimisation is being done,
         ! the restart file is written inside the runmin routine.
      case (LMOD_NTMIN_XMIN)
         write(6,'(a)') '  LMOD XMIN Minimization.'
         write(6,'(a)') ''
         write(6,'(a)') '  Note: Owing to the behaviour of the XMIN algorithm,'
         write(6,'(a)') '        coordinates in the trajectory and intermediate'
         write(6,'(a)') '        restart files will not match up with energies'
         write(6,'(a)') '        in the mdout and mdinfo files. The final energy'
         write(6,'(a)') '        and final coordinates do match.'
         write(6,'(a)') ''
         xmin_iter = 0
         call run_xmin( x, ix, ih, ipairs, &
               x(lcrd), x(lforce), ene, qsetup, xmin_iter, ntpr )
         if (master) call minrit(0,nrp,ntxo,x(lcrd))  ! Write the restart file
      case (LMOD_NTMIN_LMOD)
         write(6,'(a)') '  LMOD LMOD Minimization.'
         write(6,'(a)') ''
         write(6,'(a)') '  Note: Owing to the behaviour of the XMIN algorithm,'
         write(6,'(a)') '        coordinates in the trajectory and intermediate'
         write(6,'(a)') '        restart files will not match up with energies'
         write(6,'(a)') '        in the mdout and mdinfo files. The final energy'
         write(6,'(a)') '        and final coordinates do match.'
         write(6,'(a)') ''
         call run_lmod( x, ix, ih, ipairs, &
               x(lcrd), x(lforce), ene, qsetup )
         if (master) call minrit(0,nrp,ntxo,x(lcrd))  ! Write the restart file
      case default
         ! invalid ntmin
         ! ntmin input validation occurs in mdread.f
         ASSERT( .false. )
      end select
      

   case (5)
      ! ---carlos modified for reading trajectories (trajene option)

      write (6,*) "POST-PROCESSING OF TRAJECTORY ENERGIES"

      ! ---read trajectories and calculate energies for each frame

      call trajene(x,ix,ih,ipairs,ene,ok,qsetup)

      if (.not. ok) then
         write (6,*) 'error in trajene()'
         call mexit(6,1)
      end if

   case default
      ! invalid imin
      ! imin input validation should be transferred to mdread.f
      write(6,'(/2x,a,i3,a)') 'Error: Invalid IMIN (',imin,').'
      ASSERT( .false. )
   end select

#ifdef MPI /* SOFT CORE */
   if (master) then
      if (icfe /=0 .and. ifsc == 1) call summarize_ti_changes(natom,resat)
   end if
#endif

   ! finish up EMAP 
   if (temap) call qemap()

   if (abfqmmm_param%abfqmmm /= 1) then                 ! lam81
    exit                                                ! lam81
   else                                                 ! lam81
#ifdef MPI /* lam81 */
    call mpi_barrier(commsander,ierr)                   ! lam81
#endif /* lam81 */
    if (abfqmmm_param%qmstep /= abfqmmm_param%maxqmstep .or. abfqmmm_param%system /= 2) then  ! lam81
      if(qmmm_nml%ifqnt) call deallocate_qmmm(qmmm_nml, qmmm_struct, qmmm_vsolv, qm2_params)  ! lam81
    end if                                                                                    ! lam81
    call deallocate_m1m2m3()                            ! lam81
    if (abfqmmm_param%system == 2 .and. master) then    ! lam81
      call abfqmmm_write_idrst()                        ! lam81
      call abfqmmm_write_pdb(x(lcrd),ix(i70))           ! lam81
    end if                                              ! lam81
    call abfqmmm_next_step()                            ! lam81
   end if                                               ! lam81

   end do ! lam81

   if(abfqmmm_param%abfqmmm == 1) then          ! lam81
    deallocate(abfqmmm_param%id, stat=ier)      ! lam81
    REQUIRE( ier == 0 )                         ! lam81 
    deallocate(abfqmmm_param%v, stat=ier)       ! lam81
    REQUIRE( ier == 0 )                         ! lam81
    deallocate(abfqmmm_param%f, stat=ier)       ! lam81
    REQUIRE( ier == 0 )                         ! lam81
    deallocate(abfqmmm_param%f1, stat=ier)      ! lam81
    REQUIRE( ier == 0 )                         ! lam81
    deallocate(abfqmmm_param%f2, stat=ier)      ! lam81
    REQUIRE( ier == 0 )                         ! lam81
    if(master) then                             ! lam81
     deallocate(abfqmmm_param%isqm, stat=ier)   ! lam81
     REQUIRE( ier == 0 )                        ! lam81
    endif                                       ! lam81
   end if                                       ! lam81

   ! -- calc time spent running vs setup
   call timer_stop(TIME_TOTAL)
   call wallclock( time1 )
   call date_and_time( final_date, final_time )
   call profile_time( time1 - time0, num_calls_nblist, profile_mpi)

#ifdef MPI
   ! =========================== AMBER/MPI ===========================

   ! Set and broadcast notdone in mpi_orig case to inform
   ! other nodes that we are finished calling force(). (tec3)

   if (mpi_orig) then
      notdone = 0
      call mpi_bcast(notdone,1,mpi_integer,0, commsander,ierr)
   end if

   ! ========================= END AMBER/MPI =========================
#endif

! ========================= PUPIL INTERFACE =========================
#ifdef PUPIL_SUPPORT
   ! Finalize Corba Interface
   puperror = 0
   call killcorbaintfc(puperror)
   if (puperror /= 0) then
      write(6,*) 'Error ending PUPIL CORBA interface.'
   end if
   write(6,'(a)') 'PUPIL CORBA interface finalized.'
   pupactive = .false.
#endif
! ========================= PUPIL INTERFACE =========================

#ifdef _XRAY
   !(out_lun,residue_pointer,residue_label,atom_name,coor,num_bonds,ibond,jbond)
   call xray_fini()
#endif

   call amflsh(6)
   if (master) then
#ifdef MPI
      ! adaptive QM/MM (qmmm_nml%vsolv=2) via multisander
      ! all groups have identical coords and velocities
      ! only master of first group needs to dump results
      if ( (qmmm_nml%vsolv < 2) .or. (worldrank == 0) ) then
#endif
         call close_dump_files
#ifdef MPI
      end if
#endif
      
      if (icnstph /= 0) &
         call cnstph_finalize
      ! --- write out final times, taking REMD into account
#ifdef MPI
      if (rem .ne. 0) then
         nstlim_total = nstlim * numexchg
      else
         nstlim_total = nstlim
      end if
#else
      nstlim_total = nstlim
#endif
      if (imin == 0) &
         call print_ongoing_time_summary(nstlim_total,nstlim_total,dt,6)

      write(6,'(12(a))') '|           Job began  at ', initial_time(1:2), &
            ':', initial_time(3:4), ':', initial_time(5:10), '  on ',&
            initial_date(5:6), '/', initial_date(7:8), '/', initial_date(1:4)
      write(6,'(12(a))') '|           Setup done at ', setup_end_time(1:2),  &
            ':', setup_end_time(3:4), ':', setup_end_time(5:10), '  on ', &
            setup_end_date(5:6), '/',setup_end_date(7:8),'/',setup_end_date(1:4)
      write(6,'(12(a))') '|           Run   done at ', final_time(1:2),  &
            ':', final_time(3:4), ':', final_time(5:10), '  on ', &
            final_date(5:6), '/', final_date(7:8), '/', final_date(1:4)
      call nwallclock( ncalls )
      write(6, '(''|'',5x,''wallclock() was called'',I8,'' times'')') ncalls
      call amflsh(6)

      if (iesp > 0) then
         call esp(natom,x(lcrd),x(linddip))
      end if
   end if
   call amflsh(6)

   ! --- dynamic memory deallocation:
   999 continue 

   if(qmmm_nml%ifqnt .and. qmmm_nml%qmtheory%EXTERN .and. master) then
     call qm2_extern_finalize()
   endif

   if (qmmm_nml%ifqnt .and. .not. qmmm_struct%qm_mm_first_call) then   
      ! If first_call is still true, this thread never really
      ! called the QMMM routine. E.g. more threads than PIMD replicates
      call deallocate_qmmm(qmmm_nml, qmmm_struct, qmmm_vsolv, qm2_params)
   end if

   if (ipimd > 0) call pimd_finalize(ipimd)

   if (ineb > 0) call neb_finalize()

   if (idecomp > 0) then
      call deallocate_real_decomp()
      call deallocate_int_decomp()
   end if
   if (master .and. idecomp == 0) call deallocate_int_decomp()

#ifdef MPI /* SOFT CORE */
   if (ifsc /= 0) call cleanup_sc()
   if (ifmbar /= 0) call cleanup_mbar()
#endif

   ! finalize LIE module if initiated above
   if ( ilrt /= 0 ) then
      call cleanup_linear_response(master)
   end if

#ifdef RISMSANDER
   call rism_finalize()
#endif

   if ( ifcr /= 0 ) call cr_cleanup() 

   if (sebomd_obj%do_sebomd) then
     if (master) then
       call sebomd_close_files
     end if
     call cleanup_sebomd_arrays
   end if

   if (master .and. iwrap == 2) then
      deallocate(iwrap_mask_atoms, stat=ier)
      REQUIRE(ier == 0)
   end if
   call nblist_deallocate()
   call deallocate_stacks()
   if ((igb /= 0 .and. igb /= 10 .and. ipb == 0) .or. hybridgb > 0 .or. icnstph > 1) then
      call deallocate_gb()
   end if
   if (master) then  
      if(igb == 10 .or. ipb /= 0) then
         call pb_free()
      end if
   end if
   deallocate(ih, stat = ier)
   REQUIRE(ier == 0)
   deallocate(ipairs, stat = ier)
   REQUIRE(ier == 0)
   deallocate(ix, stat = ier)
   REQUIRE(ier == 0)
   deallocate(x, stat = ier)
   REQUIRE(ier == 0)
   if(ntb > 0 .and. ifbox == 1 .and. ew_type == 0 .and. mpoltype == 0) &
         call deallocate_m1m2m3()
   call AMOEBA_deallocate
   ! GMS: -- Module molecule --
   call deallocate_molecule()
   ! --------------------------

   if (charmm_active) call charmm_deallocate_arrays()
   if (cmap_active) call deallocate_cmap_arrays()

   if (master.and.mdout /= 'stdout') close(6)

   return

end subroutine sander

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Calculate the ElectroStatic Potential
subroutine esp(natom,x,mom_ind)

   ! routine to calculate the ESP due to the induced moments (only)
   ! at the same spatial points as the reference QM.
   use constants, only : zero, BOHRS_TO_A, INV_AMBER_ELECTROSTATIC
   use file_io_dat
   implicit none
   integer natom
   _REAL_  x(3,*)
   _REAL_  mom_ind(3,*)

#  include "ew_mpole.h"

   integer dat_unit, new_unit, minus_new_unit
   parameter(dat_unit=30, new_unit=31, minus_new_unit=33)

   integer inat, nesp, idum
   _REAL_  xin, yin, zin
   integer jn, kn
   _REAL_  esp_qm, xb_esp, yb_esp, zb_esp
   _REAL_  x_esp, y_esp, z_esp
   _REAL_  e_x, e_y, e_z, e_q, esp_new
   _REAL_  dist, dist3
   integer iptr

   call amopen(dat_unit,"esp.dat",'O','F','R')
   call amopen(new_unit,"esp.induced",owrite,'F','W')
   call amopen(minus_new_unit,"esp.qm-induced",owrite,'F','W')
   read (dat_unit,'(3i5)')inat,nesp,idum
   write(6,'(t2,''inat = '',i5)')inat
   write(6,'(t2,''nesp = '',i5)')nesp

   write(new_unit,'(2i5)')inat,nesp
   write(minus_new_unit,'(2i5)')inat,nesp

   if (inat /= natom) then
      write(6,'(t2,''natom mismatch with esp file'')')
      call mexit(6,1)
   end if

   do jn = 1,inat
      read (dat_unit,'(17x,3e16.0)')xin,yin,zin
      write(new_unit,'(17x,3e16.7)')xin,yin,zin
      write(minus_new_unit,'(17x,3e16.7)')xin,yin,zin
   end do

   do jn = 1,nesp
      e_x = zero
      e_y = zero
      e_z = zero
      e_q = zero
      read(dat_unit,'(1x,4e16.0)')esp_qm,xb_esp,yb_esp,zb_esp
      x_esp = xb_esp * BOHRS_TO_A
      y_esp = yb_esp * BOHRS_TO_A
      z_esp = zb_esp * BOHRS_TO_A

      do kn = 1,natom
         dist = (sqrt((x(1,kn)-x_esp)**2 + &
               (x(2,kn)-y_esp)**2 + &
               (x(3,kn)-z_esp)**2))
         dist3 = dist**3
         e_x = e_x - mom_ind(1,kn   )*(x(1,kn)-x_esp)/dist3
         e_y = e_y - mom_ind(2,kn   )*(x(2,kn)-y_esp)/dist3
         e_z = e_z - mom_ind(3,kn   )*(x(3,kn)-z_esp)/dist3
      end do

      e_x = e_x * BOHRS_TO_A * INV_AMBER_ELECTROSTATIC
      e_y = e_y * BOHRS_TO_A * INV_AMBER_ELECTROSTATIC
      e_z = e_z * BOHRS_TO_A * INV_AMBER_ELECTROSTATIC
      e_q = e_q * BOHRS_TO_A * INV_AMBER_ELECTROSTATIC
      esp_new = e_x + e_y + e_z

      write(new_unit,      '(1x,4e16.7)')esp_new, &
            xb_esp,yb_esp,zb_esp
      write(minus_new_unit,'(1x,4e16.7)')esp_qm-esp_new, &
            xb_esp,yb_esp,zb_esp
   end do

   close(dat_unit)
   close(new_unit)
   close(minus_new_unit)

   return

end subroutine esp

