#include "copyright.h"
#include "../include/dprec.fh"
#include "ncsu-config.h"
#include "../include/assert.fh"
#ifndef PBSA
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Open input files and read cntrl namelist.
subroutine mdread1()

   use file_io_dat
   use lmod_driver, only : read_lmod_namelist
   use qmmm_module, only : qmmm_nml,qmmm_struct, qm_gb
   use constants, only : RETIRED_INPUT_OPTION, zero, one, two, three, seven, &
                         eight, NO_INPUT_VALUE_FLOAT, NO_INPUT_VALUE, &
                         plumed, plumedfile
   use constantph, only : mccycles
   use amoeba_mdin, only: AMOEBA_read_mdin, iamoeba
   use nose_hoover_module, only: nchain  ! APJ
   use lscivr_vars, only: ilscivr, icorf_lsc
   use pimd_vars, only: ipimd,itimass
   use neb_vars, only: ineb
   use cmd_vars, only: restart_cmd, eq_cmd, adiab_param
   use stack, only: lastist,lastrst
   use nmr, only: echoin
   use crg_reloc, only: ifcr, cropt, crcut, crskin, crin, crprintcharges
   use sgld, only : isgld, isgsta,isgend,fixcom, &
                    tsgavg,sgft,sgff,sgfd,tempsg,treflf,tsgavp
   use amd_mod, only: iamd,iamdlag,EthreshD,alphaD,EthreshP,alphaP, &
        w_amd,EthreshD_w,alphaD_w,EthreshP_w,alphaP_w
   use scaledMD_mod, only: scaledMD,scaledMD_lambda
   use nbips, only: ips,teips,tvips,teaips,tvaips,raips,mipsx,mipsy,mipsz, &
                    mipso,gridips,dvbips
   use emap,only: temap,gammamap
#ifdef DSSP
   use dssp, only: idssp
#endif /* DSSP */

   use emil_mod,          only : emil_do_calc
   use mdin_emil_dat_mod, only : error_hdr

#if !defined(DISABLE_NCSU) && defined(MPI)
   use ncsu_sander_hooks, only : ncsu_on_mdread1 => on_mdread1
#endif /* ! DISABLE_NCSU && MPI */
#ifdef _XRAY
   use xray_interface_module, only: xray_active, xray_read_mdin
#endif /* _XRAY */
#ifdef MPI /* SOFT CORE */
   use softcore, only : scalpha,scbeta,ifsc,scmask,logdvdl,dvdl_norest,dynlmb, &
                        sceeorder, tishake, emil_sc
   use mbar, only : ifmbar, bar_intervall, bar_l_min, bar_l_max, bar_l_incr
   use remd, only  : rem
#endif /* MPI */
   ! Parameter for LIE module
   use linear_response, only: ilrt, lrt_interval, lrtmask
#ifdef RISMSANDER
   use sander_rism_interface, only: rismprm,xvvfile, guvfile, huvfile, cuvfile,&
        uuvfile, asympfile, quvFile, chgDistFile
#endif /*RISMSANDER*/
#ifdef APBS
   use apbs
#endif /* APBS */
   use sebomd_module, only: read_sebomd_namelist, sebomd_namelist_default
   implicit none
#  include "box.h"
#  include "def_time.h"
#  include "ew_cntrl.h"
#  include "ew_pme_recip.h"
#  include "../include/md.h"
#  include "../include/memory.h"
#  include "mmtsb.h"
#  include "nmr.h"
#  include "tgtmd.h"
#  include "multitmd.h"
#  include "ew_erfc_spline.h"
#ifdef LES
#  include "les.h"
#else
   _REAL_      temp0les
#endif

   character(len=4) watdef(4),watnam,owtnm,hwtnm1,hwtnm2

   _REAL_      dele
   integer     ierr
   integer     ifind
   integer     imcdo
   integer     itotst
   integer     jn
   integer     inerr
   logical     mdin_cntrl, mdin_lmod, mdin_qmmm  ! true if namelists are in mdin
   logical     mdin_sebomd
   integer :: ifqnt    ! local here --> put into qmmm_nml%ifqnt after read here
   integer     mxgrp
   integer     iemap
   character(len=8) date
   character(len=10) time
   character(len=512) :: char_tmp_512
   _REAL_      dtemp  ! retired 
   _REAL_      dxm  ! retired 
   _REAL_      heat  ! retired 
   _REAL_      timlim ! retired

#ifdef RISMSANDER
   integer irism
#endif /*RISMSANDER*/

   namelist /cntrl/ irest,ibelly, &
         ntx,ntxo,ntcx,ig,tempi, &
         ntb,ntt,nchain,temp0,tautp, &
         ntp,pres0,comp,taup,barostat,mcbarint, &
         nscm,nstlim,t,dt, &
         ntc,ntcc,nconp,tol,ntf,ntn,nsnb, &
         cut,dielc, &
         ntpr,ntwx,ntwv,ntwe,ntwf,ntave,ntpp,ioutfm, &
         ntr,nrc,ntrx,taur,nmropt, &
         ivcap,cutcap,xcap,ycap,zcap,fcap, &
         xlorth,ylorth,zlorth,xorth,yorth,zorth,forth, &
         imin,drms,dele,dx0, &
         pencut,ipnlty,iscale,scalm,noeskp, &
         maxcyc,ncyc,ntmin,vlimit, &
         mxsub,ipol,jfastw,watnam,owtnm,hwtnm1,hwtnm2, iesp, &
         skmin, skmax, vv,vfac, tmode, ips, &
         mipsx,mipsy,mipsz,mipso,gridips,raips,dvbips, &
         iamd,iamdlag,EthreshD,alphaD,EthreshP,alphaP, &
         w_amd,EthreshD_w,alphaD_w,EthreshP_w,alphaP_w, &
         scaledMD,scaledMD_lambda, &
         iemap,gammamap, &
         isgld,isgsta,isgend,fixcom,tsgavg,sgft,sgff,sgfd,tempsg,treflf,tsgavp,&
         jar, iamoeba, &
         numexchg, repcrd, numwatkeep, hybridgb, &
         ntwprt,tausw, &
         ntwr,iyammp,imcdo, &
         igb,alpb,Arad,rgbmax,saltcon,offset,gbsa,vrand, &
         surften,iwrap,nrespa,nrespai,gamma_ln,extdiel,intdiel, &
         cut_inner,icfe,clambda,klambda, rbornstat,lastrst,lastist,  &
         itgtmd,tgtrmsd,tgtmdfrc,tgtfitmask,tgtrmsmask, dec_verbose, &
         idecomp,temp0les,restraintmask,restraint_wt,bellymask, &
         noshakemask,crgmask, iwrap_mask,mmtsb_switch,mmtsb_iterations, &
         rdt,icnstph,solvph,ntcnstph,ntrelax, mccycles, &
         ifqnt,ievb, ipimd, itimass, ineb,profile_mpi, ilscivr, icorf_lsc, &
         ipb, inp, plumed,plumedfile, &
         gbneckscale, & 
         gbalphaH,gbbetaH,gbgammaH, &
         gbalphaC,gbbetaC,gbgammaC, &
         gbalphaN,gbbetaN,gbgammaN, &
         gbalphaOS,gbbetaOS,gbgammaOS, &
         gbalphaP,gbbetaP,gbgammaP, &
         Sh,Sc,Sn,So,Ss,Sp, &
         lj1264, &
         ifcr, cropt, crcut, crskin, crin, crprintcharges, &
         csurften, ninterface, gamma_ten, &
#ifdef MPI /* SOFT CORE */
         scalpha, scbeta, ifsc, scmask, logdvdl, dvdl_norest, dynlmb, &
         sceeorder, &
         ifmbar, bar_intervall, bar_l_min, bar_l_max, bar_l_incr, tishake, &
         emil_sc, &
#endif
         ilrt, lrt_interval, lrtmask, &
#ifdef DSSP
         idssp, &
#endif
#ifdef RISMSANDER
         irism,&
#endif /*RISMSANDER*/
         emil_do_calc, & 
         restart_cmd, eq_cmd, adiab_param,  &
         vdwmodel, & ! mjhsieh - the model used for van der Waals
         dtemp, heat, timlim  !all retired 

   ! Define default water residue name and the names of water oxygen & hydrogens
   
   data watdef/'WAT ','O   ','H1  ','H2  '/
   
   !     ----- READ THE CONTROL DATA AND OPEN DIFFERENT FILES -----
   
   if (mdout /= "stdout" ) &
         call amopen(6,mdout,owrite,'F','W')
   call amopen(5,mdin,'O','F','R')
   write(6,9308)
   call date_and_time( DATE=date, TIME=time )
   write(6,'(12(a))') '| Run on ', date(5:6), '/', date(7:8), '/',  &
        date(1:4), ' at ', time(1:2), ':', time(3:4), ':', time(5:6)

   ! Write the path of the current executable and working directory
   call get_command_argument(0, char_tmp_512)
   write(6,'(/,a,a)') '|   Executable path: ', trim(char_tmp_512)
   call getcwd(char_tmp_512)
   write(6,'(a,a)') '| Working directory: ', trim(char_tmp_512)
! Write the hostname if we can get it from environment variable
! Note: get_environment_variable is part of the F2003 standard but seems
!       to be supported by GNU, Intel, IBM and Portland (2010+) compilers
   call get_environment_variable("HOSTNAME", char_tmp_512, inerr)
   if (inerr .eq. 0) then
     write(6,'(a,a,/)') '|          Hostname: Unknown'
   else
     write(6,'(a,a,/)') '|          Hostname: ', trim(char_tmp_512)
   end if

   if (owrite /= 'N') write(6, '(2x,a)') '[-O]verwriting output'

   ! Echo the file assignments to the user:
   
   write(6,9700) 'MDIN'   ,mdin(1:70)  , 'MDOUT' ,mdout(1:70) , &
         'INPCRD' ,inpcrd(1:70), 'PARM'  ,parm(1:70)  , &
         'RESTRT',restrt(1:70) , 'REFC'  ,refc(1:70)  , &
         'MDVEL' ,mdvel(1:70)  , 'MDFRC' ,mdfrc(1:70) , &
         'MDEN'   ,mden(1:70)  , &
         'MDCRD' ,mdcrd(1:70)  , 'MDINFO' ,mdinfo(1:70), &
         'MTMD'  ,mtmd(1:70)   , 'INPDIP', inpdip(1:70), &
         'RSTDIP', rstdip(1:70), 'INPTRAJ', inptraj(1:70)
#  ifdef MPI
   write(6,9702) 'REMLOG',     trim(remlog), &
                 'REMTYPE',    trim(remtype), &
                 'REMSTRIP',   trim(remstripcoord), &
                 'SAVEENE',    trim(saveenefile), &
                 'CLUSTERINF', trim(clusterinfofile), &
                 'RESERVOIR',  trim(reservoirname), &
                 'REMDDIM',    trim(remd_dimension_file)
#  endif
#ifdef RISMSANDER
   if(len_trim(xvvfile) > 0)&
        write(6,9701) 'Xvv',trim(xvvfile)
   if(len_trim(guvfile) > 0)&
        write(6,9701) 'Guv',trim(Guvfile)
   if(len_trim(huvfile) > 0)&
        write(6,9701) 'Huv',trim(Huvfile)
   if(len_trim(cuvfile) > 0)&
        write(6,9701) 'Cuv',trim(Cuvfile)
   if(len_trim(uuvfile) > 0)&
        write(6,9701) 'Uuv',trim(Uuvfile)
   if(len_trim(asympfile) > 0)&
        write(6,9701) 'Asymptotics',trim(asympfile)
   if(len_trim(quvfile) > 0)&
        write(6,9701) 'Quv',trim(Quvfile)
   if(len_trim(chgDistfile) > 0)&
        write(6,9701) 'ChgDist',trim(chgDistfile)
#endif /*RISMSANDER*/

   ! Echo the input file to the user:
   call echoin(5,6)
   !     ----- READ DATA CHARACTERIZING THE MD-RUN -----
   read(5,'(a80)') title
   !       ----read input in namelist format, first setting up defaults
   
   dtemp = RETIRED_INPUT_OPTION
   dxm   = RETIRED_INPUT_OPTION
   heat  = RETIRED_INPUT_OPTION
   timlim = RETIRED_INPUT_OPTION
   irest = 0
   ibelly = 0
   ipol = RETIRED_INPUT_OPTION
   iesp = 0
   ntx = 1
   ntxo = NO_INPUT_VALUE
   ig = 71277
   tempi = ZERO
   ntb = NO_INPUT_VALUE
   ntt = 0
   nchain = 1
   temp0 = 300.0d0
   plumed = 0
   plumedfile = 'plumed.dat'

#ifdef LES
   ! alternate temp for LES copies, if negative then use single bath
   ! single bath not the same as 2 baths with same target T
   temp0les = -ONE
   rdt = ZERO
#endif
   ipimd =0
   itimass = 0   ! Default = no TI w.r.t. mass.
   ineb  =0

   tautp = ONE
   ntp = 0
   barostat = 1
   mcbarint = 100
   pres0 = ONE
   comp = 44.6d0
   taup = ONE
   npscal = 1
   nscm = 1000
   nstlim = 1
   t = ZERO
   dt = 0.001d0
   ntc = 1
   tol = 0.00001
   ntf = 1
   nsnb = 25
   cut =  NO_INPUT_VALUE_FLOAT
   dielc = ONE
   ntpr = 50
   ntwr = 500
   ntwx = 0
   ntwv = 0
   ntwf = 0
   ntwe = 0
   ipb = 0
   inp = 2

#ifdef RISMSANDER
   irism = 0
#endif /*RISMSANDER*/

   ntave = 0
   ioutfm = 0
   ntr = 0
   ntrx = 1
   ivcap = 0
   natcap = 0
   fcap = 1.5d0
   cutcap = 0.0d0
   xcap = 0.0d0
   ycap = 0.0d0
   zcap = 0.0d0
   forth = 1.5d0
   xlorth = -1.0d0
   ylorth = -1.0d0
   zlorth = -1.0d0
   xorth = 47114711.0d0
   yorth = 47114711.0d0
   zorth = 47114711.0d0
   numexchg = 0
   repcrd   = 1
   lj1264 = 0

   profile_mpi = 0 !whether to write profile_mpi timing file - default = 0 (NO).

   ! number of waters to keep for hybrid model,
   ! numwatkeep: the number of closest
   ! waters to keep. close is defined as close to non-water.
   ! for simulations with ions, ions should be stripped too
   ! or at least ignored in the "closest" calculation. this
   ! is not currently done.

   ! if it stays at -1 then we keep all waters
   ! 0 would mean to strip them all

    numwatkeep=-1

   ! hybridgb: gb model to use with hybrid REMD.
   hybridgb=0
 
   ! carlos targeted MD, like ntr
   
   itgtmd=0
   tgtrmsd=0.
   tgtmdfrc=0.
   tgtfitmask=''
   tgtrmsmask=''

   pencut = 0.1d0
   taumet = 0.0001d0
   omega = 500.0d0
   ipnlty = 1
   scalm = 100.0d0
   iscale = 0
   noeskp = 1
   nmropt = 0
   jar = 0
   tausw = 0.1d0
   imin = 0
   isftrp = 0
   rwell = ONE
   maxcyc = 1
   ncyc = 10
   ntmin = 1
   dx0 = 0.01d0
   drms = 1.0d-4
   vlimit = 20.0d0
   mxsub = 1
   jfastw = 0
   watnam = '    '
   owtnm =  '    '
   hwtnm1 = '    '
   hwtnm2 = '    '
   ntwprt = 0
   igb = 0
   alpb = 0
   Arad = 15.0d0
   rgbmax = 25.d0
   saltcon = ZERO

   !  default offset depends on igb value, and users need to
   !  be able to modify it, so we need to set a dummy value. if it's still the
   !  dummy after we read the namelist, we set the default based on igb. if not,
   !  we leave it at what the user set.
   !  best solution would be to create a GB namelist.
   offset = -999999.d0 
   gbneckscale = -999999.d0 

   iyammp = 0
   imcdo = -1
   gbsa = 0
   vrand=1000
   surften = 0.005d0
   iwrap = 0
   nrespa = 1
   nrespai = 1
   irespa = 1
   gamma_ln = ZERO
   extdiel = 78.5d0
   intdiel = ONE
   gbgamma = ZERO
   gbbeta = ZERO
   gbalpha = ONE

   !Hai Nguyen: set default parameters for igb = 8
   ! NOTE THAT NONE OF THESE ARE USED UNLESS IGB=8, SO USERS SHOULD NOT EVEN SET
   ! THEM
   gbalphaH = 0.788440d0
   gbbetaH = 0.798699d0
   gbgammaH = 0.437334d0
   gbalphaC = 0.733756d0
   gbbetaC = 0.506378d0
   gbgammaC = 0.205844d0
   gbalphaN = 0.503364d0
   gbbetaN = 0.316828d0
   gbgammaN = 0.192915d0
   gbalphaOS = 0.867814d0
   gbbetaOS = 0.876635d0
   gbgammaOS = 0.387882d0
   gbalphaP = 1.0d0    !P parameters are not optimized yet
   gbbetaP = 0.8d0     !P parameters are not optimized yet
   gbgammaP = 4.85d0   !P parameters are not optimized yet
   !scaling parameters below will only be used for igb=8. 
   ! the actual code does not use these variables, it uses X(l96)
   ! if igb=8, we will use these to set the X(l96) array.
   Sh = 1.425952d0
   Sc = 1.058554d0
   Sn = 0.733599d0
   So = 1.061039d0
   Ss = -0.703469d0
   Sp = 0.5d0          !P parameters are not optimized yet
   !End Hai Nguyen
         
   iconstreff = 0
   cut_inner = EIGHT
   icfe = 0
   clambda = ZERO
   klambda = 1
   ievb = 0
   rbornstat = 0
   idecomp = 0
   ! added a flag to control output of BDC/SDC synonymous with MMPBSA.py's
   ! version of the same variable.
   dec_verbose = 3 
   lastrst = 1
   lastist = 1
   restraintmask=''
   restraint_wt = ZERO
   bellymask=''
   noshakemask=''
   iwrap_mask=''  ! GMS: mask to wrap around if iwrap == 2
   crgmask=''
   mmtsb_switch = mmtsb_off ! MMTSB Replica Exchange Off by Default
   mmtsb_iterations = 100   ! MMTSB Replica Exchange Frequency in Iterations

   icnstph = 0
   solvph = SEVEN
   ntcnstph = 10
   ntrelax = 500 ! how long to let waters relax
   mccycles = 1  ! How many cycles of Monte Carlo steps to run
   skmin = 50 !used by neb calculation
   skmax = 100 !used by neb calculation
   vv = 0 !velocity verlet -- off if vv/=1
   vfac = 0 !velocity verlet scaling factor, 0 by default
   tmode = 1 !default tangent mode for NEB calculation

   ifqnt = NO_INPUT_VALUE

   ifcr = 0 ! no charge relocation
   cropt = 0 ! 1-4 EEL is calculated with the original charges
   crcut = 3.0
   crskin = 2.0 
   crin = ''
   crprintcharges = 0

   ips = 0    ! no isotropic periodic sum
   raips=-1.0d0   ! automatically determined
   mipsx=-1   ! number of grids in x direction, <0 for automatically determined
   mipsy=-1   ! number of grids in y direction, <0 for automatically determined
   mipsz=-1   ! number of grids in z direction, <0 for automatically determined
   mipso=4    ! default 4th order b-spline
   gridips=2   ! grid size. used to determine grid number if not defined
   dvbips=1.0d-8   ! Volume change tolerance. aips will be done when change more than dvbips

   iamd = 0 ! No accelerated MD used
   iamdlag = 0 !frequency of boosting in steps
   EthreshD = 0.d0
   alphaD = 0.d0
   EthreshP = 0.d0
   alphaP = 0.d0
   w_amd = 0 ! windowed amd
   EthreshD_w = 0.d0
   alphaD_w = 0.d0
   EthreshP_w = 0.d0
   alphaP_w = 0.d0
 
   scaledMD = 0 ! No scaled MD used
   scaledMD_lambda = 0.d0

   iemap=0     ! no emap constraint
   gammamap=1     ! default friction constant for map motion, 1/ps
   isgld = 0   ! no self-guiding
   isgsta=1    ! Begining index of SGLD range
   isgend=0    ! Ending index of SGLD range
   fixcom=-1    ! fix center of mass in SGLD simulation
   tsgavg=0.2d0    !  Local averaging time of SGLD simulation
   sgft=-1.0d3      !  Guiding factor of SGLD simulation
   sgff=-1.0d3      !  Guiding factor of SGLD simulation
   sgfd=-1.0d3      !  Guiding factor of SGLD simulation
   tempsg=0.0d0    !  Guiding temperature of SGLD simulation
   treflf=0.0d0    !  Reference low frequency temperature of SGLD simulation
   tsgavp=2.0d0    !  Convergency time of SGLD simulation

   !     Check to see if "cntrl" namelist has been defined.
   mdin_cntrl=.false.
   mdin_qmmm = .false.
   mdin_ewald=.false.
   mdin_pb=.false.
#ifdef APBS
   mdin_apbs = .false.
#endif /* APBS */
   mdin_lmod=.false.
   mdin_amoeba=.false.
   mdin_sebomd=.false.
   iamoeba = 0
#ifdef MPI /* SOFT CORE */
   scalpha=0.5
   scbeta=12.0
   sceeorder=2
   ifsc=0
   logdvdl=0
   dvdl_norest=0
   dynlmb=0.0
   ifmbar=0
   bar_intervall=100
   bar_l_min=0.1
   bar_l_max=0.9
   bar_l_incr=0.1
   tishake = 0
   emil_sc = 0
#endif
   ilrt = 0
   lrt_interval = 50
   lrtmask=''
#ifdef DSSP
   idssp = 0
#endif
   emil_do_calc = 0
   
!  Constant Surface Tension
   csurften = 0      !constant surface tension off (valid options are 0,1,2,3)
   gamma_ten = 0.0d0 !0.0 dyne/cm - default used in charmm. Ignored for csurften=0
   ninterface = 2   !Number of interfaces in the surface tension (Must be greater than 2)

   call nmlsrc('cntrl',5,ifind)
   if (ifind /= 0) mdin_cntrl=.true.

   call nmlsrc('ewald',5,ifind)
   if (ifind /= 0) mdin_ewald=.true.

   call nmlsrc('pb',5,ifind)
   if (ifind /= 0) mdin_pb=.true.
   
   call nmlsrc('qmmm', 5, ifind)
   if (ifind /= 0) mdin_qmmm = .true.

#ifdef APBS
   call nmlsrc('apbs',5,ifind)
   if (ifind /= 0) mdin_apbs=.true.
#endif /* APBS */

   call nmlsrc('lmod',5,ifind)
   if (ifind /= 0) mdin_lmod=.true.

   call nmlsrc('amoeba',5,ifind)
   if (ifind /= 0) mdin_amoeba=.true.

   call nmlsrc('sebomd',5,ifind)
   if (ifind /= 0) mdin_sebomd=.true.

#ifdef _XRAY
   call nmlsrc('xray',5,ifind)
   xray_active = (ifind /= 0)
#endif

   rewind 5
   if ( mdin_cntrl ) then
      read(5,nml=cntrl,err=999)
   else
      write(6, '(1x,a,/)') 'Could not find cntrl namelist'
      call mexit(6,1)
   end if

   if ( igb == 10 .and. ipb == 0 ) ipb = 2
   if ( igb == 0  .and. ipb /= 0 ) igb = 10

   if (plumed.eq.1) then
     write(6, '(1x,a,/)') 'PLUMED is on'
     write(6, '(1x,a,a,/)') 'PLUMEDfile is ',plumedfile
   endif
   
   if (ifqnt == NO_INPUT_VALUE) then
      ifqnt = 0 ! default value
      if (mdin_qmmm) then
         write(6, '(1x,a,/)') &
            '| WARNING qmmm namelist found, but ifqnt was not set! QMMM NOT &
            &active.'
      end if
   end if

   ! Now that we've read the input file, set up the defaults for variables
   ! whose values depend on other input values (ntb, cut)
   if (ntb == NO_INPUT_VALUE) then
      if (ntp > 0) then
         ntb = 2
      else if (igb > 0) then
         ntb = 0
      else
         ntb = 1
      end if
   end if

   if (ntxo == NO_INPUT_VALUE) then
#ifdef MPI
      if (rem < 0) then
         ntxo = 2
      else
         ntxo = 1
      end if
#else
      ntxo = 1
#endif
   end if

   if (cut == NO_INPUT_VALUE_FLOAT) then
      if (igb == 0) then
         cut = EIGHT
      else
         cut = 9999.d0
      end if
   end if

#ifdef RISMSANDER
   !force igb=6 to get vacuum electrostatics.  This must be done ASAP to ensure SANDER's
   !electrostatics are initialized properly
   rismprm%irism=irism
   if(irism/=0) then
      write(6,'(a)') "|3D-RISM Forcing igb=6"
      igb=6
   end if
#endif /*RISMSANDER*/   


   if (ifqnt>0) then
     qmmm_nml%ifqnt = .true.
     if (saltcon /= 0.0d0) then
       qm_gb%saltcon_on = .true.
     else
       qm_gb%saltcon_on = .false.
     end if
     if (alpb == 1) then
       qm_gb%alpb_on = .true.
     else
       qm_gb%alpb_on = .false.
     end if
     if (igb == 10 .or. ipb /= 0) then
       write(6, '(1x,a,/)') 'QMMM is not compatible with Poisson Boltzmann (igb=10 or ipb/=0).'
       call mexit(6,1)
     end if
   else
     qmmm_nml%ifqnt = .false.
   end if

   if ( mdin_lmod ) then
      rewind 5
      call read_lmod_namelist()
   end if

   !--------------------------------------------------------------------
   !     --- vars have been read ---
   !--------------------------------------------------------------------
   
   write(6,9309)
   
   ! emit warnings for retired cntrl namelist variables

   if ( dtemp /= RETIRED_INPUT_OPTION ) then
      write(6,'(/,a,/,a,/,a)') 'Warning: dtemp has been retired.', &
            '  Check the Retired Namelist Variables Appendix in the manual.'
   end if
   if ( dxm /= RETIRED_INPUT_OPTION ) then
      write(6,'(/,a,/,a,/,a)') 'Warning: dxm has been retired.', &
            '  Check the Retired Namelist Variables Appendix in the manual.'
            ! '  The step length will be unlimited.'
   end if
   if ( heat /= RETIRED_INPUT_OPTION ) then
      write(6,'(/,a,/,a,/,a)') 'Warning: heat has been retired.', &
            '  Check the Retired Namelist Variables Appendix in the manual.'
   end if
   
   if ( timlim /= RETIRED_INPUT_OPTION ) then
      write(6,'(/,a,/,a,/,a)') 'Warning: timlim has been retired.', &
            '  Check the Retired Namelist Variables Appendix in the manual.'
   end if

! Constant surface tension valid options
   if (csurften > 0) then
      if (csurften < 0 .or. csurften > 3) then
         write(6,'(/2x,a)') &
         'Invalid csurften value. csurften must be between 0 and 3'
         call mexit(6,1)
      end if
      if (ntb /= 2) then
         write(6,'(/2x,a)') &
         'ntb invalid. ntb must be 2 for constant surface tension.'
         call mexit(6,1)
      end if
      if (ntp < 2) then
         write(6,'(/2x,a)') &
         'ntp invalid. ntp must be 2 or 3 for constant surface tension.'
         call mexit(6,1)
      end if
      if (ninterface < 2) then
         write(6,'(/2x,a)') &
         'ninterface must be greater than 2 for constant surface tension.'
         call mexit(6,1)
      end if
   
      if (iamoeba > 0 ) then
         write(6,'(/2x,a)') &
         'Constant Surface Tension is not compatible with Amoeba Runs.'
         call mexit(6,1)
      end if
   
      if (ipimd > 0 ) then
         write(6,'(/2x,a)') &
         'Constant Surface Tension is not compatible with PIMD Runs.'
         call mexit(6,1)
      end if

   end if

! MC Barostat valid options. Some of these may work, but disable them until they
! are fully tested.
   
   if (ntp > 0 .and. barostat == 2) then
      inerr = 0
      if (ievb /= 0) then
         write(6, '(/2x,a)') 'AMOEBA is not compatible with the MC Barostat'
         inerr = 1
      end if
      if (ipimd /= 0) then
         write(6, '(/2x,a)') 'PIMD is not compatible with the MC Barostat'
         inerr = 1
      end if
      if (icfe /= 0) then
         write(6, '(/2x,a)') 'TI is not compatible with the MC Barostat'
         inerr = 1
      end if
#ifdef LES
      write(6, '(/2x,a)') 'LES is not compatible with the MC Barostat'
      inerr = 1
#endif
      ! Any others? Hopefully most or all of the above can be made compatible.
      if (inerr == 1) &
         call mexit(6, 1)
   end if

   call printflags()

   !--------------------------------------------------------------------
   ! If user has requested ewald electrostatics, read some more input
   !--------------------------------------------------------------------
   
   if( igb == 0 .and. ipb == 0 ) call load_ewald_info(parm,inpcrd,ntp)

   !--------------------------------------------------------------------
   ! parameters for IPS and for SGLD:
   ! ips=1  3D IPS for electrostatic and Lennard-Jones potentials
   ! ips=2  3D IPS for electrostatic potential only
   ! ips=3  3D IPS for  Lennard-Jones potential only
   ! ips=4  3D IPS/DFFT for electrostatic and Lennard-Jones potentials
   ! ips=5  3D IPS/DFFT for electrostatic potential only
   ! ips=6  3D IPS/DFFT for Lennard-Jones potential only
   !--------------------------------------------------------------------
   
   teips=.false.
   tvips=.false.
   teaips=.false.
   tvaips=.false.
   if((ips-4)*(ips-6) == 0 )tvaips =.true.
   if ( (ips-4)*(ips-5) == 0 )teaips =.true.
   if( tvaips.OR.( (ips -1)*(ips-3) == 0 ))tvips =.true.
   if( teaips.OR.((ips -1)*(ips-2) == 0 ))teips =.true.
   if( teips ) then
      use_pme = 0
      eedmeth = 6
   end if
   if( tvips ) then
      vdwmeth = 2
      if(use_pme/=0.and.tvaips)then
        mipsx=nfft1   ! number of grids in x direction, <0 for automatically determined
        mipsy=nfft2   ! number of grids in y direction, <0 for automatically determined
        mipsz=nfft3   ! number of grids in z direction, <0 for automatically determined
        mipso=order    ! default 6th order b-spline
      endif
   end if
   temap=iemap>0
   ishake = 0
   if (ntc > 1) ishake = 1
   
   !--------------------------------------------------------------------
   ! Set up some parameters for AMD simulations:
   ! AMD initialization
   ! iamd=0 no boost is used, 1 boost on the total energy, 
   ! 2 boost on the dohedrals, 3 boost on dihedrals and total energy
   !--------------------------------------------------------------------
   if(iamd.gt.0)then
      if(iamd.eq.1)then !only total potential energy will be boosted
         EthreshD=0.d0
         alphaD=0.d0
      else if(iamd.eq.2)then !only dihedral energy will be boosted
         EthreshP=0.d0
         alphaP=0.d0
      endif
      if(w_amd.gt.0)then
         if(iamd.eq.1)then !only total potential energy will be boosted
            EthreshD_w=0.d0
            alphaD_w=0.d0
         else if(iamd.eq.2)then !only dihedral energy will be boosted
            EthreshP_w=0.d0
            alphaP_w=0.d0
         endif
         write(6,'(a,i3)')'| Using Windowed Accelerated MD (wAMD) &
                          &LOWERING BARRIERS to enhance sampling w_amd =', w_amd
         write(6,'(a,2f22.12)')'| AMD boost to total energy: EthreshP,alphaP',&
                           EthreshP, alphaP
         write(6,'(a,2f22.12)')'| AMD boost to dihedrals: EthreshD,alphaD',&
                           EthreshD,alphaD
         write(6,'(a,2f22.12)')'| AMD extra parameters boost to total energy: &
                           &EthreshP_w,alphaP_w', EthreshP_w, alphaP_w
         write(6,'(a,2f22.12)')'| AMD extra parameters boost to dihedrals: &
                           &EthreshD_w,alphaD_w', EthreshD_w, alphaD_w
      else
         write(6,'(a,i3)')'| Using Accelerated MD (AMD) RASING VALLEYS to &
                           &enhance sampling iamd =',iamd
         write(6,'(a,2f22.12)')'| AMD boost to total energy: EthreshP,alphaP', &
                           EthreshP, alphaP
         write(6,'(a,2f22.12)')'| AMD boost to dihedrals: EthreshD,alphaD', &
                           EthreshD, alphaD
      endif
   endif
 

   !--------------------------------------------------------------------
   ! Set up some parameters for scaledMD simulations:
   ! scaledMD initialization
   ! scaledMD=0 no scaling is used, 1 scale the potential energy 
   !--------------------------------------------------------------------
   if(scaledMD.gt.0)then
      write(6,'(a,i3)')'| Using Scaled MD to enhance sampling scaledMD =',&
                        scaledMD
      write(6,'(a,f22.12)')'| scaledMD scaling factor lambda: ',scaledMD_lambda
   endif
   



   !--------------------------------------------------------------------
   ! Set up some parameters for GB simulations:
   !--------------------------------------------------------------------
   !Hai Nguyen: update offset = 0.09d0 for igb /= 8
   !I add this step because I want to use different offset value as default value
   !for igb = 8
   
   if ( igb == 8 ) then
     if (offset == -999999.d0) then
        offset = 0.195141d0  !set to default for igb=8
     end if 
     if (gbneckscale == -999999.d0) then
        gbneckscale = 0.826836d0  
     end if 
   else 
      ! not igb=8, use old defaults
      if (offset == -999999.d0) then 
         offset = 0.09d0
      end if 
      if (gbneckscale == -999999.d0) then
         gbneckscale = 0.361825d0
      end if 
   endif
   
   if( igb == 2 .or. hybridgb == 2 ) then
      !       --- use our best guesses for Onufriev/Case GB  (GB^OBC I)
      
      gbgamma = 2.90912499999d0  ! (the "99999" to force roundoff on print)
      gbbeta = ZERO
      gbalpha = 0.8d0
   end if

   if( igb == 5 .or. hybridgb == 5 ) then
      
      !       --- use our second best guesses for Onufriev/Case GB (GB^OBC II)
      
      gbgamma = 4.850d0
      gbbeta = 0.8d0
      gbalpha = ONE
   end if
   
   if( igb == 7 ) then
      
      !       --- use parameters for Mongan et al. CFA GBNECK
      
      gbgamma = 2.50798245d0
      gbbeta = 1.90792938d0
      gbalpha = 1.09511284d0
   end if
   
   !--------------------------------------------------------------------
   ! If user has requested PB electrostatics, read some more input
   !--------------------------------------------------------------------

   if ( igb == 10 .or. ipb /= 0 ) then
#ifdef MPI
      write(6,'(a)') "PBSA currently doesn't work with MPI inside SANDER."
      call mexit(6,1)
#endif /*MPI*/
      call pb_read
   end if

#ifdef APBS
   if ( mdin_apbs ) then
      call apbs_read
   end if
#endif /* APBS */

#ifdef _XRAY
   call xray_read_mdin(mdin_lun=5)
#endif

   call sebomd_namelist_default
   if (mdin_sebomd) then
     rewind 5
     call read_sebomd_namelist
   endif

   if( iamoeba == 1 ) then
      if( mdin_amoeba ) then
         call AMOEBA_read_mdin(5)
      else
        write(6,*) ' iamoeba is set but the &amoeba namelist was not found'
        call mexit(6,1)
      end if
   end if
   
   ! -------------------------------------------------------------------
   ! If the user has requested NMR restraints, do a cursory read of the
   ! restraints file(s) now to determine the amount of memory necessary
   ! for these restraints:
   ! -------------------------------------------------------------------
   
   if (jar == 1 ) nmropt = 1
   intreq = 0
   irlreq = 0
   if (nmropt > 0) then
      mxgrp = 0
      itotst = 1
      
      ! Set ITOTST to 0 if IMIN equals 1 (i.e. if minimization, not dynamics)
      ! This will cause any "time-averaged" requests to be over-ridden.
      
      if (imin == 1) then 
         itotst = 0
      end if 
      !         CALL AMOPEN(31,NMR,'O','F','R')
      call restlx(5,itotst,mxgrp,dt,6,ierr)
      !         CLOSE(31)
   end if

   ! -------------------------------------------------------------------
   ! If EMIL was requested, make sure it was compiled in, and validate.
   ! -------------------------------------------------------------------
  if( emil_do_calc .gt. 0 ) then
#ifdef EMIL 
    if( ntc .ne. 1 ) then
       write (6, '(a,a)') error_hdr, 'emil_do_calc == 1,'
       write (6, '(a)') '      and ntc != 1.'
       write (6, '(a)') '      Current thinking is that SHAKE and '
       write (6, '(a)') '      EMIL do not mix well, consider setting ntc = 1, ntf = 1,'
       write (6, '(a)') '      and dt = 0.001.'
       call mexit(6,1)
    end if
#else
    write (6, '(a,a)') error_hdr, 'emil_do_calc = 1,'
    write (6, '(a)') '      but AMBER was compiled with EMIL switched out.'
    write (6, '(a)') '      Run $AMBERHOME/configure --help for more info.'
    call mexit(6,1)
#endif
  end if


   ! Set the definition of the water molecule. The default definition is in
   ! WATDEF(4).
   
   read(watdef(1),'(A4)') iwtnm
   read(watdef(2),'(A4)') iowtnm
   read(watdef(3),'(A4)') ihwtnm(1)
   read(watdef(4),'(A4)') ihwtnm(2)
   if (watnam /= '    ') read(watnam,'(A4)') iwtnm
   if (owtnm /= '    ') read(owtnm, '(A4)') iowtnm
   if (hwtnm1 /= '    ') read(hwtnm1,'(A4)') ihwtnm(1)
   if (hwtnm2 /= '    ') read(hwtnm2,'(A4)') ihwtnm(2)

#if !defined(DISABLE_NCSU) && defined(MPI)
   call ncsu_on_mdread1()
#endif

   return
   999 continue   ! bad cntrl read
   write(6,*) 'error in reading namelist cntrl'
   call mexit(6,1)
   
   ! --- input file polar opts read err trapping:
   
   9308 format(/10x,55('-'),/10x, &
         'Amber 14 SANDER                              2014', &
         /10x,55('-')/)
   9309 format(/80('-')/'   1.  RESOURCE   USE: ',/80('-')/)
   9700 format(/,'File Assignments:',/,15('|',a6,': ',a,/))
   9701 format('|',a6,': ',a)
   9702 format(7('|',a10,': ',a,/))
end subroutine mdread1 
#endif /*ifndef PBSA*/

#ifndef PBSA
!======================================================================
!          MDREAD2
!======================================================================

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Initialize to defaults and print the inputable variables.
subroutine mdread2(x,ix,ih,ipairs)
   
   use molecule, only: n_iwrap_mask_atoms, iwrap_mask_atoms
   use lmod_driver, only : LMOD_NTMIN_LMOD, LMOD_NTMIN_XMIN, write_lmod_namelist
   use decomp, only : jgroup, indx, irespw
   use findmask
   use qmmm_module, only : qmmm_nml,qmmm_struct, qm2_struct, qm2_rij_eqns, qm_gb, qmmm_vsolv
   use qmmm_vsolv_module, only : print
   use constants, only : ZERO, ONE, TWO
   use parms, only: req
   use nbips, only: ips
   use amd_mod, only: iamd,iamdlag,EthreshD,alphaD,EthreshP,alphaP, &
        w_amd,EthreshD_w,alphaD_w,EthreshP_w,alphaP_w
   use scaledMD_mod, only: scaledMD
   use nblist, only: a,b,c,alpha,beta,gamma,nbflag,skinnb,sphere,nbtell,cutoffnb
   use amoeba_mdin, only : iamoeba,beeman_integrator
   use amoeba_runmd, only : AM_RUNMD_get_coords
   use nose_hoover_module, only: nchain  ! APJ
   use pimd_vars, only: ipimd, itimass
   use neb_vars, only: last_neb_atom, ineb
   use cmd_vars, only: restart_cmd, eq_cmd, adiab_param
   use constantph, only: cnstphread, cnstph_zero, cph_igb, mccycles
   use file_io_dat
   use sander_lib, only: upper

   use emap,only: temap,emap_options

   use qmmm_module, only: get_atomic_number

   use sebomd_module, only: sebomd_obj, sebomd_write_info, sebomd_write_options

#ifdef APBS
   use apbs
#endif /* APBS */
#ifdef EMIL
   use emil_mod,          only : emil_do_calc
   use mdin_emil_dat_mod, only : error_hdr, init_emil_dat
#endif /* EMIL */
#ifdef _XRAY
   use xray_interface_module, only: xray_write_options
#endif /* _XRAY */
#if defined(LES) && defined(MPI)
   use evb_pimd, only: evb_pimd_init
#endif /* LES && MPI */
#ifdef MPI /* SOFT CORE */
   use softcore, only : ifsc, scmask, scalpha, scbeta, dvdl_norest, &
                        sceeorder, logdvdl, dynlmb, tishake
   use mbar, only : ifmbar, bar_intervall, bar_l_min, bar_l_max, bar_l_incr
! REMD
   use remd, only : rem, rremd
#endif
   use crg_reloc, only: ifcr, cropt, crcut, crskin, crprintcharges
   use linear_response, only: ilrt, lrt_interval, lrtmask
! SGLD
   use sgld, only : isgld

   implicit none
   _REAL_ x(*)
#  include "../include/memory.h"
   integer ix(lasti),ipairs(*)
   character(len=4) ih(*)
   integer nbond
   integer atom1,atom2
   integer ntmp
   logical belly,konst
   character(len=1) atsymb,atsymb2
   character(len=2) atype
   integer atomicnumber, hybridization
   integer ngrp,inerr,nr,iaci,ir,i,mxresat,j
   integer noshakegp( natom ), natnos
   integer iwrap_maskgp( natom ) , ier
   logical errFlag
   integer crggp( natom )
   _REAL_ dummy,rvdw,dcharge,emtmd
   logical newstyle
   
#ifdef MPI
   !     =========================== AMBER/MPI ===========================
#ifdef MPI_DOUBLE_PRECISION
#  undef MPI_DOUBLE_PRECISION
#endif
   include 'mpif.h'
#  include "parallel.h"
   integer ist(MPI_STATUS_SIZE), partner, ierr, nbonh_c, num_noshake_c
   integer nquant_c, noshake_overlap_c
#ifdef CRAY_PVP
#  define MPI_DOUBLE_PRECISION MPI_REAL8
#endif
   !     ========================= END AMBER/MPI =========================
#endif
#  include "../include/md.h"
#  include "box.h"
#  include "mmtsb.h"
#  include "nmr.h"
#  include "extra_pts.h"
#  include "ew_cntrl.h"
#  include "ew_pme_recip.h"
#  include "ew_erfc_spline.h"
#  include "ew_mpole.h"
#  include "ew_legal.h"
#  include "def_time.h"
#  include "tgtmd.h"
#  include "multitmd.h"
#ifdef LES
#  include "les.h"
#endif
   
   ! -------------------------------------------------------------------
   !      --- set up resat array, containing string identifying
   !          residue for each atom
   ! -------------------------------------------------------------------
   
   mxresat = min( natom, matom )
   ir = 0
   do i=1,mxresat
      if (i >= ix(ir+i02)) ir=ir+1
      write(resat(i),'(a4,1x,a4,i4)') ih(m04+i-1), &
            ih(ir+m02-1),ir
      !                                    ---null terminator:
      resat(i)(14:14) = char(0)
   end do
   close(unit=8)
   
   ! -------------------------------------------------------------------
   !     ----- SET THE DEFAULT VALUES FOR SOME VARIABLES -----
   ! -------------------------------------------------------------------
   
   nrp = natom
   
   if (ifbox == 1) write(6, '(/5x,''BOX TYPE: RECTILINEAR'')')
   if (ifbox == 2) write(6, '(/5x,''BOX TYPE: TRUNCATED OCTAHEDRON'')')
   if (ifbox == 3) write(6, '(/5x,''BOX TYPE: GENERAL'')')
   
   ! For 0< ipimd <= 3, no removal of COM motion
   if (ipimd>0.and.ipimd<=3) then
      nscm = 0
      ndfmin = 0
   endif

   if (ntr.eq.1) &
      nscm = 0

   nsolut =  nrp
   if ( nscm > 0 .and. ntb == 0 ) then
      ndfmin = 6   ! both translation and rotation com motion removed
      if (nsolut == 1) ndfmin = 3
      if (nsolut == 2) ndfmin = 5
   else if ( nscm > 0 ) then
      ndfmin = 3    ! just translation com will be removed
   else
      ndfmin = 0
   end if
   if (ibelly > 0) then   ! No COM Motion Removal, ever.
      nscm = 0
      ndfmin = 0

      !Do not allow ntt=3 with ibelly=1 - this does not make sense and will cause
      !issues.
      if (ntt==3) then
        call sander_bomb("mdread2","ibelly=1 with ntt=3 is not a valid option.", &
                         "Either use a different thermostat or avoid using ibelly.")
      end if
   end if
   if(nscm <= 0) nscm = 0
   if(gamma_ln > 0.0d0)ndfmin=0  ! No COM motion removal for LD simulation
   if(ntt == 4)ndfmin=0  ! No COM motion removal for Nose'-Hoover simulation
   init = 3
   if (irest > 0) init = 4
   if (dielc <= ZERO ) dielc = ONE
   if (tautp <= ZERO ) tautp = 0.2d0
   if (taup <= ZERO ) taup = 0.2d0
   
   !     ----- RESET THE CAP IF NEEDED -----
   
   ! ivcap == 0: Cap will be in effect if it is in the prmtop file (ifcap = 1)

   if(ivcap == 1) then
      ! Cap will be in effect even if not in prmtop file
      !   requires additional information in sander.in file as in the case of ivcap == 3, 4, or 5
      ifcap = 2
   else if(ivcap == 2) then
      ! Inactivate cap
      ifcap = 0
   else if(ivcap == 3) then
      ! Sphere -> not yet implemented
      ifcap = 3
   else if(ivcap == 4) then
      ! Orthorhombus
      ifcap = 4
   else if(ivcap == 5) then
      ! Shell of waters around solute
      ifcap = 5
   end if

   !Support for random seed using time of day in microsecs
   if( ig==-1 ) then
     !Turn on NO_NTT3_SYNC when ig=-1. This means the code no
     !longer synchronized the random numbers between streams when
     !running in parallel giving better scaling.
     no_ntt3_sync = 1
     call microsec(ig)
#ifdef MPI
     write (6, '(a,i8,a)') "Note: ig = -1. Setting random seed to ", ig ," based on wallclock &
                               &time in microseconds"
     write (6, '(a)') "      and disabling the synchronization of random &
                               &numbers between tasks"
     write (6, '(a)') "      to improve performance."
#else 
     write (6, '(a,i8,a)') "Note: ig = -1. Setting random seed to ", ig ," based on wallclock &
                                &time in microseconds."
#endif
   end if

#ifdef MPI
   !For some runs using multisander where the coordinates of the two 'replicas' need to be identical,
   !for example TI, it is critical that the random number stream is synchronized between all replicas.
   !Only use the IG value from worldrank=0. Ok to broadcast between the various sander masters since
   !they all call mdread2.
   !Also needs to be synchronized for adaptive QM/MM (qmmm_nml%vsolv > 1)
   if ( (icfe > 0) .or. (qmmm_nml%vsolv > 1) ) then
      ! no_ntt3_sync currently does not work with softcore TI simulations 
      ! see sc_lngdyn in softcore.F90
      ! AWG: I think I also need to set no_ntt3_sync=0 ???
      if ( (ifsc > 0) .or. (qmmm_nml%vsolv > 1) ) no_ntt3_sync = 0 
      call mpi_bcast(ig, 1, MPI_INTEGER, 0, commmaster, ierr)
   end if
#endif
 
   ! -------------------------------------------------------------------
   !     ----- PRINT DATA CHARACTERIZING THE RUN -----
   ! -------------------------------------------------------------------
   
   nr = nrp
   write(6,9328)
   write(6,9008) title
   write(6,'(/a)') 'General flags:'
   write(6,'(5x,2(a,i8))') 'imin    =',imin,', nmropt  =',nmropt

   ! Error Checking for REMD
#ifdef MPI
   if (rem/=0) then
      ! Make sure that the number of replicas is even so
      !  that they all have partners in the exchange step
      if (mod(numgroup,2).ne.0) then
         write (6,'(a)') "==================================="
         write (6,'(a)') "REMD requires an even # of replicas"
         write (6,'(a)') "==================================="
         call mexit (6,1)
      endif

      write(6,'(/a)') 'Replica exchange'
      write(6,'(5x,4(a,i8))') 'numexchg=',numexchg,', rem=',rem

      ! REPCRD option temporarily disabled
      if (repcrd == 0) write(6,'(a)') &
         "REMD WARNING: repcrd disabled. Only replica &
         &trajectories/output can be written."

      ! Check for correct number of exchanges
      if (numexchg <= 0) then
         write(6,'(a)') "REMD ERROR: numexchg must be > 0, "
         call mexit(6,1)
      endif

      ! RXSGLD
      if (isgld > 0 .and. rem < 0) then
         write(6, '(a)') 'Multi-D REMD and replica-exchange SGLD are not &
                         &supported yet!'
         call mexit(6, 1)
      end if

      ! Hybrid GB
      if (numwatkeep >= 0) then
         write(6,'(5x,4(a,i8))') 'numwatkeep=',numwatkeep,', hybridgb=',hybridgb
         ! Check that user specified GB model for hybrid REMD
         if (hybridgb /= 2 .and. hybridgb /= 5 .and. hybridgb /= 1) then
            write(6,'(a)') "HYBRID REMD ERROR: hybridgb must be 1, 2, or 5."
            call mexit(6,1)
         endif
      else
         !Check that user did not specify GB model if no hybrid run.
         if (hybridgb /= 0) then
            write(6,'(a)') &
            "HYBRID REMD ERROR: numwatkeep must be >= 0 when hybridgb is set."
            call mexit(6,1)
         endif
      endif
      ! RREMD
      if (rremd>0) then
         write(6,'(5x,4(a,i8))') "rremd=",rremd
      endif

      ! ntp > 0 not allowed for remd 
      if (ntp > 0) then
         write(6,'(a,i1)') "ERROR: ntp > 0 not allowed for rem > 0, ntp=", ntp
         call mexit(6,1)
      endif

      ! M-REMD (rem < 0) requires netcdf output.
      if (rem < 0 .and. ioutfm .ne. 1) then
         write(6,'(a)') "ERROR: Multi-D REMD (rem < 0) requires NetCDF &
                        &trajectories (ioutfm=1)"
         call mexit(6,1)
      endif

#  ifdef LES
      ! DAN ROE: Temporarily disable LES REMD until it is verified with new
      !           REMD code
      if (rem==2) then
         write (6,*) "******* LES REM (rem==2) temporarily disabled. Stop. *******"
         call mexit(6,1)
      endif

      if (rem==2 .and. igb/=1) then
         write (6,*) ' partial REM (rem=2) only works with igb=1'
         call mexit(6,1)
      endif
#  else
      if (rem==2) then
         write(6,*) '*******  For rem ==2, partial REM'
         write(6,*) 'use sander.LES with topology created by addles'
         call mexit(6,1)
      endif
#  endif

   endif ! rem>0

#else /* _NO_ MPI */
   ! Check if user set numexchg with no MPI
   if (numexchg>0) write(6,'(a)') &
      "WARNING: numexchg > 0 - for REMD run please recompile sander for &
      &parallel runs."

   ! Check if user set numwatkeep or hybridgb with no MPI - not sensible.
   if (numwatkeep>=0) write(6,'(a)') &
      "WARNING: numwatkeep >= 0 - for hybrid REMD run please recompile &
      &sander for parallel runs."

   if (hybridgb>0) write(6,'(a)') &
      "WARNING: hybridgb > 0 - for hybrid REMD run please recompile &
      &sander for parallel runs."
#endif /* MPI */
   ! End error checking for REMD

   write(6,'(/a)') 'Nature and format of input:'
   write(6,'(5x,4(a,i8))') 'ntx     =',ntx,', irest   =',irest, &
         ', ntrx    =',ntrx

   write(6,'(/a)') 'Nature and format of output:'
   write(6,'(5x,4(a,i8))') 'ntxo    =',ntxo,', ntpr    =',ntpr, &
         ', ntrx    =',ntrx,', ntwr    =',ntwr
   write(6,'(5x,5(a,i8))') 'iwrap   =',iwrap,', ntwx    =',ntwx, &
         ', ntwv    =',ntwv,', ntwe    =',ntwe
   write(6,'(5x,3(a,i8),a,i7)') 'ioutfm  =',ioutfm, &
         ', ntwprt  =',ntwprt, &
         ', idecomp =',idecomp,', rbornstat=',rbornstat
   if (ntwf > 0) &
      write(6,'(5x, a,i8)') 'ntwf    =',ntwf
   write(6,'(/a)') 'Potential function:'
   write(6,'(5x,5(a,i8))') 'ntf     =',ntf,', ntb     =',ntb, &
         ', igb     =',igb,', nsnb    =',nsnb
   write(6,'(5x,3(a,i8))') 'ipol    =',ipol,', gbsa    =',gbsa, &
         ', iesp    =',iesp
   write(6,'(5x,3(a,f10.5))') 'dielc   =',dielc, &
         ', cut     =',cut,', intdiel =',intdiel

   ! charge relocation
   if ( ifcr /= 0 ) then
      write(6,'(/a)') 'Charge relocation:'
      write(6,'(5x,2(a,i8))') 'cropt   =', cropt, &
                                ', crprintcharges=', crprintcharges
      write(6,'(5x,2(a,f10.5))') 'crcut   =', crcut, ', crskin  =', crskin
   end if

   if (( igb /= 0 .and. igb /= 10 .and. ipb == 0 .and. igb /= 8) &
                                  .or.hybridgb>0.or.icnstph>1) then
      write(6,'(5x,3(a,f10.5))') 'saltcon =',saltcon, &
            ', offset  =',offset,', gbalpha= ',gbalpha
      write(6,'(5x,3(a,f10.5))') 'gbbeta  =',gbbeta, &
            ', gbgamma =',gbgamma,', surften =',surften
      write(6,'(5x,3(a,f10.5))') 'rdt     =',rdt, ', rgbmax  =',rgbmax, &
            '  extdiel =',extdiel
      write(6,'(5x,3(a,i8))') 'alpb  = ',alpb
   end if
   
   !Hai Nguyen: print output for igb=8
   if ( igb == 8 ) then
      write(6,'(5x,3(a,f10.5))') 'saltcon =',saltcon, &
            ', offset  =',offset,', surften =',surften
      write(6,'(5x,3(a,f10.5))') 'rdt     =',rdt, ', rgbmax  =',rgbmax, &
            '  extdiel =',extdiel
      write(6,'(5x,3(a,i8))') 'alpb  = ',alpb
      write(6,'(5x,3(a,f10.5))') 'gbalphaH  =',gbalphaH, &
            ', gbbetaH   =',gbbetaH,',  gbgammaH  = ',gbgammaH
      write(6,'(5x,3(a,f10.5))') 'gbalphaC  =',gbalphaC, &
            ', gbbetaC   =',gbbetaC,',  gbgammaC  = ',gbgammaC
      write(6,'(5x,3(a,f10.5))') 'gbalphaN  =',gbalphaN, &
            ', gbbetaN   =',gbbetaN,',  gbgammaN  = ',gbgammaN
      write(6,'(5x,3(a,f10.5))') 'gbalphaOS =',gbalphaOS, &
            ', gbbetaOS  =',gbbetaOS,',  gbgammaOS = ',gbgammaOS
      write(6,'(5x,3(a,f10.5))') 'gbalphaP  =',gbalphaP, &
            ', gbbetaP   =',gbbetaP,',  gbgammaP  = ',gbgammaP
   end if


   if( alpb /= 0 ) then
      write(6,'(5x,3(a,f10.5))') 'Arad =', Arad
   end if    

   write(6,'(/a)') 'Frozen or restrained atoms:'
   write(6,'(5x,4(a,i8))') 'ibelly  =',ibelly,', ntr     =',ntr
   if( ntr == 1 ) write(6,'(5x,a,f10.5)') 'restraint_wt =', restraint_wt

   if( imin /= 0 ) then
      if( ipimd > 0 ) then
         write(6,'(/a)') 'pimd cannot be used in energy minimization'
         stop
      end if

      write(6,'(/a)') 'Energy minimization:'
      ! print inputable variables applicable to all minimization methods.
      write(6,'(5x,4(a,i8))') 'maxcyc  =',maxcyc,', ncyc    =',ncyc, &
            ', ntmin   =',ntmin
      write(6,'(5x,2(a,f10.5))') 'dx0     =',dx0, ', drms    =',drms

      ! Input flag ntmin determines the method of minimization
      select case ( ntmin )
      case ( 0, 1, 2 )
         ! no specific output
      case ( LMOD_NTMIN_XMIN, LMOD_NTMIN_LMOD )
         call write_lmod_namelist( )
      case default
         ! invalid ntmin
         write(6,'(/2x,a,i3,a)') 'Error: Invalid NTMIN (',ntmin,').'
         stop
      end select
   else
      write(6,'(/a)') 'Molecular dynamics:'
      write(6,'(5x,4(a,i10))') 'nstlim  =',nstlim,', nscm    =',nscm, &
            ', nrespa  =',nrespa
      write(6,'(5x,3(a,f10.5))') 't       =',t, &
            ', dt      =',dt,', vlimit  =',vlimit
      
      if ( ntt == 0 .and. tempi > 0.0d0 .and. irest == 0 ) then
         write(6,'(/a)') 'Initial temperature generation:'
         write(6,'(5x,a,i8)') 'ig      =',ig
         write(6,'(5x,a,f10.5)') 'tempi   =',tempi
      else if( ntt == 1 ) then
         write(6,'(/a)') 'Berendsen (weak-coupling) temperature regulation:'
         write(6,'(5x,3(a,f10.5))') 'temp0   =',temp0, &
               ', tempi   =',tempi,', tautp   =', tautp
#ifdef LES
         write(6,'(5x,3(a,f10.5))') 'temp0LES   =',temp0les
#endif
      else if( ntt == 2 ) then
         write(6,'(/a)') 'Anderson (strong collision) temperature regulation:'
         write(6,'(5x,4(a,i8))') 'ig      =',ig, ', vrand   =',vrand
         write(6,'(5x,3(a,f10.5))') 'temp0   =',temp0, ', tempi   =',tempi
      else if( ntt == 3 ) then
         write(6,'(/a)') 'Langevin dynamics temperature regulation:'
         write(6,'(5x,4(a,i8))') 'ig      =',ig
         write(6,'(5x,3(a,f10.5))') 'temp0   =',temp0, &
               ', tempi   =',tempi,', gamma_ln=', gamma_ln
      else if( ntt == 4 ) then
         write(6,'(/a)') 'Nose-Hoover chains'
         write(6,'(5x,(a,f10.5))') 'gamma_ln=', gamma_ln
         write(6,'(5x,(a,i8))') 'number of oscillators=', nchain
      else if( ntt == 5 ) then                                       ! APJ
         write(6,'(/a)') 'Nose-Hoover chains Langevin'               ! APJ
         write(6,'(5x,4(a,i8))') 'ig      =',ig                      ! APJ
         write(6,'(5x,(a,f10.5))') 'gamma_ln=', gamma_ln             ! APJ
         write(6,'(5x,(a,i8))') 'number of oscillators=', nchain     ! APJ
      else if( ntt == 6 ) then                                       ! APJ
         write(6,'(/a)') 'Adaptive Langevin temperature regulation:' ! APJ
         write(6,'(5x,4(a,i8))') 'ig      =',ig                      ! APJ
         write(6,'(5x,3(a,f10.5))') 'temp0   =',temp0, &             ! APJ
               ', tempi   =',tempi,', gamma_ln=', gamma_ln           ! APJ
      else if( ntt == 7 ) then                                       ! APJ
         write(6,'(/a)') 'Adaptive Nose-Hoover chains'               ! APJ
         write(6,'(5x,(a,f10.5))') 'gamma_ln=', gamma_ln             ! APJ
         write(6,'(5x,(a,i8))') 'number of oscillators=', nchain     ! APJ
      else if( ntt == 8 ) then                                       ! APJ
         write(6,'(/a)') 'Adaptive Nose-Hoover chains Langevin'      ! APJ
         write(6,'(5x,4(a,i8))') 'ig      =',ig                      ! APJ
         write(6,'(5x,(a,f10.5))') 'gamma_ln=', gamma_ln             ! APJ
         write(6,'(5x,(a,i8))') 'number of oscillators=', nchain     ! APJ
      end if

      if( ntp /= 0 ) then
         write(6,'(/a)') 'Pressure regulation:'
         write(6,'(5x,4(a,i8))') 'ntp     =',ntp
         write(6,'(5x,3(a,f10.5))') 'pres0   =',pres0, &
               ', comp    =',comp,', taup    =',taup
         if (barostat == 2) then
            write(6, '(5x,a)') 'Monte-Carlo Barostat:'
            write(6, '(5x,a,i8)') 'mcbarint  =', mcbarint
         end if
      end if

      if (csurften /= 0) then
         write(6,'(/a)') 'Constant surface tension:'
         write(6,'(5x,a,i8)') 'csurften  =', csurften
         write(6,'(5x,a,f10.5,a,i8)') 'gamma_ten =', gamma_ten, ' ninterface =', ninterface
      end if

   end if

   if( ntc /= 1 ) then
      write(6,'(/a)') 'SHAKE:'
      write(6,'(5x,4(a,i8))') 'ntc     =',ntc,', jfastw  =',jfastw
      write(6,'(5x,3(a,f10.5))') 'tol     =',tol
   end if

   if( ifcap == 1 .or. ifcap == 2 .or. ifcap == 3 ) then
      write(6,'(/a)') 'Water cap:'
      write(6,'(5x,2(a,i8))') 'ivcap   =',ivcap,', natcap  =',natcap
      write(6,'(5x,2(a,f10.5))') 'fcap    =',fcap, ', cutcap  =',cutcap
      write(6,'(5x,3(a,f10.5))') 'xcap    =',xcap, ', ycap    =',ycap,    &
                                 ', zcap    =',zcap
   else if( ifcap == 4 ) then
      write(6,'(/a)') 'Orthorhombus:'
      write(6,'(5x,1(a,i8))')    'ivcap   =',ivcap
      write(6,'(5x,1(a,f10.5))') 'forth   =',forth
      write(6,'(5x,3(a,f10.5))') 'xlorth  =',xlorth,', ylorth  =',ylorth, &
                                 ', zlorth  =',zlorth
      write(6,'(5x,3(a,f10.5))') 'xorth   =',xorth, ', yorth   =',yorth,  &
                                 ', zorth   =',zorth
   else if( ifcap == 5 ) then
      write(6,'(/a)') 'Water shell:'
      write(6,'(5x,(a,i8,a,f10.5))') 'ivcap   =',ivcap,', cutcap  =',cutcap
   endif

   if( nmropt > 0 ) then
      write(6,'(/a)') 'NMR refinement options:'
      write(6,'(5x,4(a,i8))')'iscale  =',iscale,', noeskp  =',noeskp, &
            ', ipnlty  =',ipnlty,', mxsub   =',mxsub
      write(6,'(5x,3(a,f10.5))') 'scalm   =',scalm, &
            ', pencut  =',pencut,', tausw   =',tausw
   end if

   if( numextra > 0 ) then
      write(6,'(/a)') 'Extra-points options:'
      write(6,'(5x,4(a,i8))') 'frameon =',frameon, &
            ', chngmask=',chngmask
   end if

   if( ipol > 0 ) then
      write(6,'(/a)') 'Polarizable options:'
      write(6,'(5x,4(a,i8))') 'indmeth =',indmeth, &
            ', maxiter =',maxiter,', irstdip =',irstdip, &
            ', scaldip =',scaldip
      write(6,'(5x,3(a,f10.5))') &
            'diptau  =',diptau,', dipmass =',dipmass
     if ( ipol > 1 ) then
      write(6,'(5x,3(a,f10.5))') &
            'Default Thole coefficient = ',dipdamp
     end if
   end if

#ifdef MPI /* SOFT CORE */
   if( icfe /= 0 .or. ifsc/=0) then
      write(6,'(/a)') 'Free energy options:'
      write(6,'(5x,a,i8,a,i8,a,i8)')       'icfe    =', icfe   , ', ifsc    =', ifsc,    ', klambda =', klambda
      write(6,'(5x,a,f8.4,a,f8.4,a,f8.4)') 'clambda =', clambda, ', scalpha =', scalpha, ', scbeta  =', scbeta
      write(6,'(5x,a,i8,a,i8)')            'sceeorder =', sceeorder, ' dvdl_norest =', dvdl_norest
      write(6,'(5x,a,f8.4,a,i8)')          'dynlmb =', dynlmb, ' logdvdl =', logdvdl
   end if
   if ( ifmbar /= 0 ) then
      write (6,'(/a)') 'FEP MBAR options:'
      write(6,'(5x,a,i8,a,i8)')       'ifmbar  =', ifmbar, ',  bar_intervall = ', bar_intervall
      write(6,'(5x,a,f6.4,a,f6.4,a,f6.4)')  'bar_l_min =', bar_l_min, ',  bar_l_max =', bar_l_max, ',  bar_l_incr =', bar_l_incr
   end if
#endif

   if (ilrt /= 0) then
      write (6,*)
      write (6,'(a,i4,a,i4)') ' Linear Response Theory: ilrt =', ilrt, ' lrt_interval =', lrt_interval
      write (6,*)
   end if

   ! Options for TI w.r.t. mass.
   select case (itimass)
   case (0)     ! Default: no TI wrt. mass.
   case (1,2)   ! 1 = use virial est., 2 = use thermodynamic est.
      write(6,'(/a)') 'Isotope effects (thermodynamic integration w.r.t. mass):'
      write(6,'(5x,4(a,i8))') 'itimass =',itimass      
      write(6,'(5x,3(a,f10.5))') 'clambda =',clambda
      if (icfe /= 0) then 
         write(6,'(/2x,a,i2,a,i2,a)') 'Error: Cannot do TI w.r.t. both potential (icfe =', &
            icfe, ') and mass (itimass =', itimass, ').'
         stop
      endif
      if (ipimd == 0) then 
         write(6,'(/2x,a)') 'Error (IPIMD=0): TI w.r.t. mass requires a PIMD run.'
         stop       
      endif
   case default ! Invalid itimass
      write(6,'(/2x,a,i2,a)') 'Error: Invalid ITIMASS (', itimass, ' ).'
      stop
   end select

!KFW
!   call mpi_bcast ( ievb, 1, MPI_INTEGER, 0, commworld, ierr )
!   call mpi_barrier ( commworld, ierr )

!   if( ievb == 1 ) then
!KFW  write(6,'(/a)') 'EVB options:'
!KFW  write(6,'(5x,3(a,f10.5))') 'V11     =',v11,', V22    =', v22, &
!KFW         ', V12     =', v12
!kfw  write(6,'(5x,3(a,f10.5))') 'kevb    =',kevb,', evbt   =', evbt
!   end if

   if( itgtmd == 1 ) then
      write(6,'(/a)') 'Targeted molecular dynamics:'
      write(6,'(5x,3(a,f10.5))') 'tgtrmsd =',tgtrmsd, &
            ', tgtmdfrc=',tgtmdfrc
   end if
 
   if( icnstph /= 0) then
      write(6, '(/a)') 'Constant pH options:'
      if ( icnstph .ne. 1 ) &
         write(6, '(5x,a,i8)') 'icnstph =', icnstph
      write(6, '(5x,a,i8)') 'ntcnstph =', ntcnstph
      write(6, '(5x,a,f10.5)') 'solvph =', solvph
      if ( icnstph .ne. 1 ) &
         write(6,'(5x,2(a,i8))') 'ntrelax =', ntrelax, ' mccycles =', mccycles
   end if
   if( icnstph /= 2) then
      ntrelax = 0 ! needed for proper behavior of timing
   end if

   if( ntb > 0 ) then
      write(6,'(/a)') 'Ewald parameters:'
      write(6,'(5x,4(a,i8))') 'verbose =',verbose, &
            ', ew_type =',ew_type,', nbflag  =',nbflag, &
            ', use_pme =',use_pme
      write(6,'(5x,4(a,i8))') 'vdwmeth =',vdwmeth, &
            ', eedmeth =',eedmeth,', netfrc  =',netfrc
      write(6, 9002) a, b, c
      write(6, 9003) alpha, beta, gamma
      write(6, 9004) nfft1, nfft2, nfft3
      write(6, 9006) cutoffnb, dsum_tol
      write(6, 9007) ew_coeff
      write(6, 9005) order
      9002 format (5x,'Box X =',f9.3,3x,'Box Y =',f9.3,3x,'Box Z =',f9.3)
      9003 format (5x,'Alpha =',f9.3,3x,'Beta  =',f9.3,3x,'Gamma =',f9.3)
      9004 format (5x,'NFFT1 =',i5  ,7x,'NFFT2 =',i5  ,7x,'NFFT3 =',i5)
      9005 format (5x,'Interpolation order =',i5)
      9006 format (5x,'Cutoff=',f9.3,3x,'Tol   =',e9.3)
      9007 format (5x,'Ewald Coefficient =',f9.5)
   end if

   if( mmtsb_switch /= mmtsb_off ) then
      call mmtsb_print_banner()
      call mmtsb_init( temp0, clambda )
   end if

!---- QMMM Options ----

   if( qmmm_nml%ifqnt ) then
      write(6, '(/a)') 'QMMM options:'
      write(6, '(5x,"        ifqnt = True       nquant = ",i8)') &
                 qmmm_struct%nquant
      write(6, '(5x,"         qmgb = ",i8,"  qmcharge = ",i8,"   adjust_q = ",i8)') &
                 qmmm_nml%qmgb, qmmm_nml%qmcharge, qmmm_nml%adjust_q
      write(6, '(5x,"         spin = ",i8,"     qmcut = ",f8.4, "    qmshake = ",i8)') qmmm_nml%spin, &
                 qmmm_nml%qmcut, qmmm_nml%qmshake
      write(6, '(5x,"     qmmm_int = ",i8)') qmmm_nml%qmmm_int
      write(6, '(5x,"lnk_atomic_no = ",i8,"   lnk_dis = ",f8.4," lnk_method = ",i8)') &
                 qmmm_nml%lnk_atomic_no,qmmm_nml%lnk_dis, qmmm_nml%lnk_method
      if ( qmmm_nml%qmtheory%PM3 ) then 
         write(6, '(5x,"     qm_theory =     PM3")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%AM1 ) then
         write(6, '(5x,"     qm_theory =     AM1")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%AM1D ) then
         write(6, '(5x,"     qm_theory =   AM1/d")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%MNDO ) then
         write(6, '(5x,"     qm_theory =    MNDO")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%MNDOD ) then
         write(6, '(5x,"     qm_theory =  MNDO/d")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%PDDGPM3 ) then
         write(6, '(5x,"     qm_theory = PDDGPM3")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%PDDGMNDO ) then
         write(6, '(5x,"     qm_theory =PDDGMNDO")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%PM3CARB1 ) then
         write(6, '(5x,"     qm_theory =PM3CARB1")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%DFTB ) then
         write(6, '(5x,"     qm_theory =    DFTB")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%RM1 ) then
         write(6, '(5x,"     qm_theory =     RM1")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%PDDGPM3_08 ) then
         write(6, '(5x,"     qm_theory = PDDGPM3_08")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%PM6 ) then
         write(6, '(5x,"     qm_theory =     PM6")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%PM3ZNB ) then
         write(6, '(5x,"     qm_theory = PM3/ZnB")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%PM3MAIS ) then
         write(6, '(5x,"     qm_theory = PM3-MAIS")',ADVANCE='NO') 
      else if ( qmmm_nml%qmtheory%EXTERN ) then
         write(6, '(5x,"     qm_theory =     EXTERN")',ADVANCE='NO') 
      else
         write(6, '(5x,"     qm_theory = UNKNOWN!")',ADVANCE='NO') 
      end if
      write (6, '(" verbosity = ",i8)') qmmm_nml%verbosity
      if (qmmm_nml%qmqm_analyt) then
        write(6, '(5x,"       qmqmdx = Analytical")')
      else
        write(6, '(5x,"       qmqmdx = Numerical")')
      end if
      !AWG: if EXTERN in use, skip printing of options that do not apply
      EXTERN: if ( .not. qmmm_nml%qmtheory%EXTERN ) then
         if (qmmm_nml%tight_p_conv) then
            write(6, '(5x," tight_p_conv = True (converge density to SCFCRT)")')
         else
            write(6, '(5x," tight_p_conv = False (converge density to 0.05xSqrt[SCFCRT])")')
         end if
         write(6, '(5x,"      scfconv = ",e9.3,"  itrmax = ",i8)') qmmm_nml%scfconv, qmmm_nml%itrmax
         if (qmmm_nml%printcharges) then
            write(6, '(5x," printcharges = True ")',ADVANCE='NO')
         else
            write(6, '(5x," printcharges = False")',ADVANCE='NO')
         end if
         select case (qmmm_nml%printdipole)
            case (1)
               write(6, '(5x," printdipole = QM   ")',ADVANCE='NO')
            case (2)
               write(6, '(5x," printdipole = QM+MM")',ADVANCE='NO')
            case default
               write(6, '(5x," printdipole = False")',ADVANCE='NO')
         end select 
         if (qmmm_nml%peptide_corr) then
            write(6, '(5x," peptide_corr = True")')
         else
            write(6, '(5x," peptide_corr = False")')
         end if
         if (qmmm_nml%qmmmrij_incore) then
            write(6, '(4x,"qmmmrij_incore = True ")')
         else
            write(6, '(4x,"qmmmrij_incore = False")')
         end if
         if (qmmm_nml%qmqm_erep_incore) then
            write(6, '(2x,"qmqm_erep_incore = True ")')
         else
            write(6, '(2x,"qmqm_erep_incore = False")')
         end if
         if (qmmm_nml%allow_pseudo_diag) then
            write(6, '(7x,"pseudo_diag = True ")',ADVANCE='NO')
            write(6, '("pseudo_diag_criteria = ",f8.4)') qmmm_nml%pseudo_diag_criteria
         else
            write(6, '(7x,"pseudo_diag = False")')
         end if
         write(6,   '(6x,"diag_routine = ",i8)') qmmm_nml%diag_routine
      end if EXTERN
      !If ntb=0 or use_pme =0 then we can't do qm_ewald so overide what the user may
      !have put in the namelist and set the value to false.
      if (qmmm_nml%qm_ewald>0) then
        if (qmmm_nml%qm_pme) then
          write(6, '(10x,"qm_ewald = ",i8, " qm_pme = True ")') qmmm_nml%qm_ewald
        else
          write(6, '(10x,"qm_ewald = ",i8, " qm_pme = False ")') qmmm_nml%qm_ewald
        end if
          write(6, '(10x,"  kmaxqx = ",i4," kmaxqy = ",i4," kmaxqz = ",i4," ksqmaxq = ",i4)') &
                  qmmm_nml%kmaxqx, qmmm_nml%kmaxqy, qmmm_nml%kmaxqz, qmmm_nml%ksqmaxq
      else
        write(6, '(10x,"qm_ewald = ",i8, " qm_pme = False ")') qmmm_nml%qm_ewald
      end if
      !Print the fock matrix prediction params if it is in use.
      if (qmmm_nml%fock_predict>0) then
          write(6, '(6x,"fock_predict = ",i4)') qmmm_nml%fock_predict
          write(6, '(6x,"    fockp_d1 = ",f8.4," fockp_d2 = ",f8.4)') qmmm_nml%fockp_d1, qmmm_nml%fockp_d2
          write(6, '(6x,"    fockp_d2 = ",f8.4," fockp_d4 = ",f8.4)') qmmm_nml%fockp_d3, qmmm_nml%fockp_d4
      end if

      if (qmmm_nml%qmmm_switch) then
         write(6, '(7x,"qmmm_switch = True",3x,"r_switch_lo =",f8.4,3x,"r_switch_hi =",f8.4)') &
               & qmmm_nml%r_switch_lo, qmmm_nml%r_switch_hi
      !else
      !   write(6, '(7x,"qmmm_switch = False")') 
      end if

      if (qmmm_nml%printdipole==2) then
         write(6, '("|",2x,"INFO: To compute MM dipole WAT residues will be stripped")') 
      end if 
   end if

   if (qmmm_nml%vsolv > 0) then
      call print(qmmm_vsolv)
   end if

!---- SEBOMD options ----

   if (sebomd_obj%do_sebomd) then
     call sebomd_write_info()
     call sebomd_write_options()
   endif

#ifdef _XRAY
!---- XRAY Options ----
   call xray_write_options()
#endif

! ----EMAP Options-----
   if(temap)call emap_options(5)
! ---------------------

#ifdef EMIL
! ----EMIL Options-----
   if(emil_do_calc.gt.0)call init_emil_dat(5, 6)
! ---------------------
#endif


#ifdef MPI
! --- MPI TIMING OPTIONS ---
      write(6, '(/a)') '| MPI Timing options:'
      write(6, '("|",5x," profile_mpi = ",i8)') profile_mpi
! Sanity check for profile_mpi
      call int_legal_range('profile_mpi',profile_mpi,0,1)
! --------------------------
#endif
 
   cut = cut*cut
   cut_inner = cut_inner*cut_inner
   
   !------------------------------------------------------------------------
   ! If user has requested generalized born electrostatics, set up variables
   !------------------------------------------------------------------------
   
   if( igb == 0 .and. gbsa > 0 ) then
      write(0,*) 'GB/SA calculation is performed only when igb>0'
      call mexit( 6,1 )
   end if
   if( gbsa == 2 .and. &
       ((imin == 0 .and. nstlim > 1) .or. &
        (imin == 1 .and. maxcyc > 1)) ) then
      write(0,*) 'GBSA=2 only works for single point energy calc'
      call mexit( 6,1 )
   end if
#ifdef APBS
   if( igb /= 0 .and. igb /= 10 .and. ipb == 0 .and. .not. mdin_apbs) then
#else
   if (( igb /= 0 .and. igb /= 10 .and. ipb == 0 ).or.hybridgb>0.or.icnstph>1) then
#endif /* APBS */
#ifdef LES
      write(6,*) 'igb=1,5,7 are working with LES, no SA term included'
#endif
      ! igb7 uses special S_x screening params.
      ! overwrite the tinker values read from the prmtop
      if (igb == 7) then
         do i=1,natom
            if(ix(i100) .eq. 1) then
               atomicnumber = ix(i100+i)
            else
               call get_atomic_number(ih(m04+i-1), x(lmass+i-1), &
                  atomicnumber, errFlag)
#ifdef LES
               if(errFlag) then
                  call get_atomic_number(ih(m04+i-1), &
                     x(lmass+i-1)*lesfac((lestyp(i)-1)*nlesty+lestyp(i)), &
                     atomicnumber, errFlag)
               end if
#endif
            end if

            if (atomicnumber .eq. 6) then
               x(l96+i-1) = 4.84353823306d-1
            else if (atomicnumber .eq. 1) then 
               x(l96+i-1) = 1.09085413633d0
            else if (atomicnumber .eq. 7) then
               x(l96+i-1) = 7.00147318409d-1
            else if (atomicnumber .eq. 8) then
               x(l96+i-1) = 1.06557401132d0
            else if (atomicnumber .eq. 16) then
               x(l96+i-1) = 6.02256336067d-1
            else if (atomicnumber .eq. 15) then
               x(l96+i-1) = 5d-1
            else
               x(l96+i-1) = 5d-1
            end if
         end do
      end if
      
      ! Hai Nguyen: changing S_x screening params for igb = 8
      ! overwrite the tinker values read from the prmtop
      
      if (igb == 8) then
         do i=1,natom
            if (ix(i100) .eq. 1) then
               atomicnumber = ix(i100+i)
            else
               call get_atomic_number(ih(m04+i-1), x(lmass+i-1), &
                  atomicnumber, errFlag)
#ifdef LES
               if(errFlag) then
                  call get_atomic_number(ih(m04+i-1), &
                     x(lmass+i-1)*lesfac((lestyp(i)-1)*nlesty+lestyp(i)), &
                     atomicnumber, errFlag)
               end if
#endif
            end if

            if (atomicnumber .eq. 6) then
               x(l96+i-1) = Sc
            else if (atomicnumber .eq. 1) then
               x(l96+i-1) = Sh
            else if (atomicnumber .eq. 7) then
               x(l96+i-1) = Sn
            else if (atomicnumber .eq. 8) then
               x(l96+i-1) = So
            else if (atomicnumber .eq. 16) then
               x(l96+i-1) = Ss 
            else if (atomicnumber .eq. 15) then
               x(l96+i-1) = Sp ! Hai Nguyen: We still don't have an optimized Sp parameter
            else
               !for atom type Cl,Br,...
               !These parameters are also not optimized.
               x(l96+i-1) = 5d-1
            end if
         end do
      end if


      ! Hai Nguyen: set up for igb == 2, 5, 7, 8
      ! Put gb parameters in arrays 
      if ( igb == 2 .or. igb == 5 .or. igb == 7 .or. &
           hybridgb == 2 .or. hybridgb == 5) then
        do i=1,natom
            x(l2402+i-1) = gbalpha
            x(l2403+i-1) = gbbeta
            x(l2404+i-1) = gbgamma
        end do
      end if
      
      !Hai Nguyen: IGB = 8
        if ( igb == 8 ) then
            do i=1,natom
                if (ix(i100) .eq. 1) then
                   atomicnumber = ix(i100+i)
                else
                   call get_atomic_number(ih(m04+i-1), x(lmass+i-1), &
                      atomicnumber, errFlag)
#ifdef LES
                   if(errFlag) then
                      call get_atomic_number(ih(m04+i-1), &
                         x(lmass+i-1)*lesfac((lestyp(i)-1)*nlesty+lestyp(i)), &
                         atomicnumber, errFlag)
                   end if
#endif
                end if

                if (atomicnumber .eq. 1) then
                        x(l2402+i-1) = gbalphaH
                        x(l2403+i-1) = gbbetaH
                        x(l2404+i-1) = gbgammaH
                else if (atomicnumber .eq. 6) then
                        x(l2402+i-1) = gbalphaC
                        x(l2403+i-1) = gbbetaC
                        x(l2404+i-1) = gbgammaC
                else if (atomicnumber .eq. 7) then
                        x(l2402+i-1) = gbalphaN
                        x(l2403+i-1) = gbbetaN
                        x(l2404+i-1) = gbgammaN
                else if (atomicnumber .eq. 8) then
                        x(l2402+i-1) = gbalphaOS
                        x(l2403+i-1) = gbbetaOS
                        x(l2404+i-1) = gbgammaOS
                else if (atomicnumber .eq. 16) then
                        x(l2402+i-1) = gbalphaOS
                        x(l2403+i-1) = gbbetaOS
                        x(l2404+i-1) = gbgammaOS 
                else if (atomicnumber .eq. 15) then
                        x(l2402+i-1) = gbalphaP
                        x(l2403+i-1) = gbbetaP
                        x(l2404+i-1) = gbgammaP
                else
                        !use GBOBC set for other atom types
                        x(l2402+i-1) = 1.0d0 
                        x(l2403+i-1) = 0.8d0 
                        x(l2404+i-1) = 4.85d0 
                end if
            end do
        end if ! end Hai Nguyen section
        
     !       put fs(i)*(rborn(i) - offset) into the "fs" array
      
      fsmax = 0.d0
      do i=1,natom
         x(l96-1+i) = x(l96-1+i)*( x(l97-1+i) - offset )
         fsmax = max( fsmax, x(l96-1+i) )
         if (rbornstat == 1) then
            x(l186-1+i) = 0.d0
            x(l187-1+i) = 999.d0
            x(l188-1+i) = 0.d0
            x(l189-1+i) = 0.d0
         end if
      end do
      
      !     ---------------------------------------------------------------------
      !       ---get Debye-Huckel kappa (A**-1) from salt concentration (M), assuming:
      !         T = 298.15, epsext=78.5,
      
      kappa = sqrt( 0.10806d0 * saltcon )
      
      !       ---scale kappa by 0.73 to account(?) for lack of ion exclusions:
      
      kappa = 0.73d0* kappa

      !Set kappa for qmmm if needed
      qm_gb%kappa = kappa
      !     ---------------------------------------------------------------------
      
      if ( gbsa == 1 ) then
         
         !     --- assign parameters for calculating SASA according to the
         !         LCPO method ---
         
         do i=1,natom
            ix(i80+i-1)=0
         end do
         
         !         --- get the number of bonded neighbors for each atom:
         
         do i=1,nbona
            atom1=ix(iiba+i-1)/3+1
            atom2=ix(ijba+i-1)/3+1
            ix(i80+atom1-1)=ix(i80+atom1-1)+1
            ix(i80+atom2-1)=ix(i80+atom2-1)+1
         end do

         do i=1,natom
            ix(i80-i)=ix(i80+i-1)
         end do

         do i=1,nbonh
            atom1=ix(iibh+i-1)/3+1
            atom2=ix(ijbh+i-1)/3+1
            ix(i80-atom1)=ix(i80-atom1)+1
            ix(i80-atom2)=ix(i80-atom2)+1
         end do

         
         !         --- construct parameters for SA calculation; note that the
         !             radii stored in L165 are augmented by 1.4 Ang.
         
         do i=1,natom
            write(atype,'(a2)') ih(m06+i-1)
            call upper(atype)
            if (ix(i100) .eq. 1) then
               atomicnumber = ix(i100+i)
            else
               call get_atomic_number(ih(m04+i-1), x(lmass+i-1), &
                  atomicnumber, errFlag)
#ifdef LES
               if(errFlag) then
                  call get_atomic_number(ih(m04+i-1), &
                     x(lmass+i-1)*lesfac((lestyp(i)-1)*nlesty+lestyp(i)), &
                     atomicnumber, errFlag)
               end if
#endif
            end if
            hybridization = ix(i80-i)
            nbond=ix(i80+i-1)
            if (atomicnumber .eq. 6) then
               if (hybridization .eq. 4) then
                  if (nbond == 1) then
                     x(l165-1+i) = 1.70d0 + 1.4d0
                     x(l170-1+i) = 0.77887d0
                     x(l175-1+i) = -0.28063d0
                     x(l180-1+i) = -0.0012968d0
                     x(l185-1+i) = 0.00039328d0
                  else if (nbond == 2) then
                     x(l165-1+i) = 1.70d0 + 1.4d0
                     x(l170-1+i) = 0.56482d0
                     x(l175-1+i) = -0.19608d0
                     x(l180-1+i) = -0.0010219d0
                     x(l185-1+i) = 0.0002658d0
                  else if (nbond == 3) then
                     x(l165-1+i) = 1.70d0 + 1.4d0
                     x(l170-1+i) = 0.23348d0
                     x(l175-1+i) = -0.072627d0
                     x(l180-1+i) = -0.00020079d0
                     x(l185-1+i) = 0.00007967d0
                  else if (nbond == 4) then
                     x(l165-1+i) = 1.70d0 + 1.4d0
                     x(l170-1+i) = 0.00000d0
                     x(l175-1+i) = 0.00000d0
                     x(l180-1+i) = 0.00000d0
                     x(l185-1+i) = 0.00000d0
                  else
                     write(6,*) 'Unusual nbond for CT:', i, nbond, &
                        ' Using default carbon LCPO parameters'
                     x(l165-1+i) = 1.70d0 + 1.4d0
                     x(l170-1+i) = 0.77887d0
                     x(l175-1+i) = -0.28063d0
                     x(l180-1+i) = -0.0012968d0
                     x(l185-1+i) = 0.00039328d0
                  end if
               else
                  if (nbond == 2) then
                     x(l165-1+i) = 1.70d0 + 1.4d0
                     x(l170-1+i) = 0.51245d0
                     x(l175-1+i) = -0.15966d0
                     x(l180-1+i) = -0.00019781d0
                     x(l185-1+i) = 0.00016392d0
                  else if (nbond == 3) then
                     x(l165-1+i) = 1.70d0 + 1.4d0
                     x(l170-1+i) = 0.070344d0
                     x(l175-1+i) = -0.019015d0
                     x(l180-1+i) = -0.000022009d0
                     x(l185-1+i) = 0.000016875d0
                  else
                     write(6,*) 'Unusual nbond for C :', i, nbond, &
                       ' Using default carbon LCPO parameters'
                     x(l165-1+i) = 1.70d0 + 1.4d0
                     x(l170-1+i) = 0.77887d0
                     x(l175-1+i) = -0.28063d0
                     x(l180-1+i) = -0.0012968d0
                     x(l185-1+i) = 0.00039328d0
                  end if
               end if
            else if (atomicnumber .eq. 8) then
               if (atype == 'O ') then
                  x(l165-1+i) = 1.60d0 + 1.4d0
                  x(l170-1+i) = 0.68563d0
                  x(l175-1+i) = -0.1868d0
                  x(l180-1+i) = -0.00135573d0
                  x(l185-1+i) = 0.00023743d0
               else if (atype == 'O2') then
                  x(l165-1+i) = 1.60d0 + 1.4d0
                  x(l170-1+i) = 0.88857d0
                  x(l175-1+i) = -0.33421d0
                  x(l180-1+i) = -0.0018683d0
                  x(l185-1+i) = 0.00049372d0
               else
                  if (nbond == 1) then
                     x(l165-1+i) = 1.60d0 + 1.4d0
                     x(l170-1+i) = 0.77914d0
                     x(l175-1+i) = -0.25262d0
                     x(l180-1+i) = -0.0016056d0
                     x(l185-1+i) = 0.00035071d0
                  else if (nbond == 2) then
                     x(l165-1+i) = 1.60d0 + 1.4d0
                     x(l170-1+i) = 0.49392d0
                     x(l175-1+i) = -0.16038d0
                     x(l180-1+i) = -0.00015512d0
                     x(l185-1+i) = 0.00016453d0
                  else
                     write(6,*) 'Unusual nbond for O:', i, nbond, &
                        ' Using default oxygen LCPO parameters'
                     x(l165-1+i) = 1.60d0 + 1.4d0
                     x(l170-1+i) = 0.77914d0
                     x(l175-1+i) = -0.25262d0
                     x(l180-1+i) = -0.0016056d0
                     x(l185-1+i) = 0.00035071d0
                  end if
               end if
            else if(atomicnumber .eq. 7) then
               if (atype == 'N3') then
                  if (nbond == 1) then
                     x(l165-1+i) = 1.65d0 + 1.4d0
                     x(l170-1+i) = 0.078602d0
                     x(l175-1+i) = -0.29198d0
                     x(l180-1+i) = -0.0006537d0
                     x(l185-1+i) = 0.00036247d0
                  else if (nbond == 2) then
                     x(l165-1+i) = 1.65d0 + 1.4d0
                     x(l170-1+i) = 0.22599d0
                     x(l175-1+i) = -0.036648d0
                     x(l180-1+i) = -0.0012297d0
                     x(l185-1+i) = 0.000080038d0
                  else if (nbond == 3) then
                     x(l165-1+i) = 1.65d0 + 1.4d0
                     x(l170-1+i) = 0.051481d0
                     x(l175-1+i) = -0.012603d0
                     x(l180-1+i) = -0.00032006d0
                     x(l185-1+i) = 0.000024774d0
                  else
                     write(6,*) 'Unusual nbond for N3:', i, nbond, &
                        ' Using default nitrogen LCPO parameters'
                     x(l165-1+i) = 1.65d0 + 1.4d0
                     x(l170-1+i) = 0.078602d0
                     x(l175-1+i) = -0.29198d0
                     x(l180-1+i) = -0.0006537d0
                     x(l185-1+i) = 0.00036247d0
                  end if
               else
                  if (nbond == 1) then
                     x(l165-1+i) = 1.65d0 + 1.4d0
                     x(l170-1+i) = 0.73511d0
                     x(l175-1+i) = -0.22116d0
                     x(l180-1+i) = -0.00089148d0
                     x(l185-1+i) = 0.0002523d0
                  else if (nbond == 2) then
                     x(l165-1+i) = 1.65d0 + 1.4d0
                     x(l170-1+i) = 0.41102d0
                     x(l175-1+i) = -0.12254d0
                     x(l180-1+i) = -0.000075448d0
                     x(l185-1+i) = 0.00011804d0
                  else if (nbond == 3) then
                     x(l165-1+i) = 1.65d0 + 1.4d0
                     x(l170-1+i) = 0.062577d0
                     x(l175-1+i) = -0.017874d0
                     x(l180-1+i) = -0.00008312d0
                     x(l185-1+i) = 0.000019849d0
                  else
                     write(6,*) 'Unusual nbond for N:', i, nbond, &
                        ' Using default nitrogen LCPO parameters'
                     x(l165-1+i) = 1.65d0 + 1.4d0
                     x(l170-1+i) = 0.078602d0
                     x(l175-1+i) = -0.29198d0
                     x(l180-1+i) = -0.0006537d0
                     x(l185-1+i) = 0.00036247d0
                  end if
               end if
            else if(atomicnumber .eq. 16) then
               if (atype == 'SH') then
                  x(l165-1+i) = 1.90d0 + 1.4d0
                  x(l170-1+i) = 0.7722d0
                  x(l175-1+i) = -0.26393d0
                  x(l180-1+i) = 0.0010629d0
                  x(l185-1+i) = 0.0002179d0
               else
                  x(l165-1+i) = 1.90d0 + 1.4d0
                  x(l170-1+i) = 0.54581d0
                  x(l175-1+i) = -0.19477d0
                  x(l180-1+i) = -0.0012873d0
                  x(l185-1+i) = 0.00029247d0
               end if
            else if (atomicnumber .eq. 15) then
               if (nbond == 3) then
                  x(l165-1+i) = 1.90d0 + 1.4d0
                  x(l170-1+i) = 0.3865d0
                  x(l175-1+i) = -0.18249d0
                  x(l180-1+i) = -0.0036598d0
                  x(l185-1+i) = 0.0004264d0
               else if (nbond == 4) then
                  x(l165-1+i) = 1.90d0 + 1.4d0
                  x(l170-1+i) = 0.03873d0
                  x(l175-1+i) = -0.0089339d0
                  x(l180-1+i) = 0.0000083582d0
                  x(l185-1+i) = 0.0000030381d0
               else
                  write(6,*) 'Unusual nbond for P:', i, nbond, &
                     ' Using default phosphorus LCPO parameters'
                  x(l165-1+i) = 1.90d0 + 1.4d0
                  x(l170-1+i) = 0.3865d0
                  x(l175-1+i) = -0.18249d0
                  x(l180-1+i) = -0.0036598d0
                  x(l185-1+i) = 0.0004264d0
               end if
            else if (atype(1:1) == 'Z') then
               x(l165-1+i) = 0.00000d0 + 1.4d0
               x(l170-1+i) = 0.00000d0
               x(l175-1+i) = 0.00000d0
               x(l180-1+i) = 0.00000d0
               x(l185-1+i) = 0.00000d0
            else if (atomicnumber .eq. 1) then
               x(l165-1+i) = 0.00000d0 + 1.4d0
               x(l170-1+i) = 0.00000d0
               x(l175-1+i) = 0.00000d0
               x(l180-1+i) = 0.00000d0
               x(l185-1+i) = 0.00000d0
            else if (atype == 'MG') then
               !  Mg radius = 0.99A: ref. 21 in J. Chem. Phys. 1997, 107, 5422
               !  Mg radius = 1.18A: ref. 30 in J. Chem. Phys. 1997, 107, 5422
               !  Mg radius = 1.45A: Aqvist 1992
               x(l165-1+i) = 1.18d0 + 1.4d0
               !  The following values were taken from O.sp3 with two bonded 
               !  neighbors -> O has the smallest van der Waals radius 
               ! compared to all other elements which had been parametrized
               x(l170-1+i) = 0.49392d0
               x(l175-1+i) = -0.16038d0
               x(l180-1+i) = -0.00015512d0
               x(l185-1+i) = 0.00016453d0
            else if (atype == 'F') then
               x(l165-1+i) = 1.47d0 + 1.4d0
               x(l170-1+i) = 0.68563d0
               x(l175-1+i) = -0.1868d0
               x(l180-1+i) = -0.00135573d0
               x(l185-1+i) = 0.00023743d0 
            else
               ! write( 0,* ) 'bad atom type: ',atype
               ! call mexit( 6,1 )
               x(l165-1+i) = 1.70 + 1.4;
               x(l170-1+i) = 0.51245;
               x(l175-1+i) = -0.15966;
               x(l180-1+i) = -0.00019781;
               x(l185-1+i) = 0.00016392;
               write(6,'(a,a)') 'Using carbon SA parms for atom type', atype 
            end if
         end do  !  i=1,natom
         !
      else if ( gbsa == 2 ) then

         !     --- assign parameters for calculating SASA according to the
         !         ICOSA method; the radii are augmented by 1.4 A ---

         do i=1,natom
            write(atype,'(a2)') ih(m06+i-1)
            if(ix(i100) .eq. 1) then
               atomicnumber = ix(i100+i)
            else
               call get_atomic_number(ih(m04+i-1), x(lmass+i-1), &
                  atomicnumber, errFlag)
#ifdef LES
               if(errFlag) then
                  call get_atomic_number(ih(m04+i-1), &
                     x(lmass+i-1)*lesfac((lestyp(i)-1)*nlesty+lestyp(i)), &
                     atomicnumber, errFlag)
               end if
#endif
            end if
            if (atomicnumber .eq. 7) then
               x(L165-1+i) = 1.55d0 + 1.4d0
            else if (atomicnumber .eq. 6) then
               x(L165-1+i) = 1.70d0 + 1.4d0
            else if (atomicnumber .eq. 1) then
               x(L165-1+i) = 1.20d0 + 1.4d0
            else if (atomicnumber .eq. 8) then
               x(L165-1+i) = 1.50d0 + 1.4d0
            else if (atomicnumber .eq. 15) then
               x(L165-1+i) = 1.80d0 + 1.4d0
            else if (atomicnumber .eq. 16) then
               x(L165-1+i) = 1.80d0 + 1.4d0
            else if (atomicnumber .eq. 17) then
               ! Cl radius
               x(L165-1+i) = 1.70d0 + 1.4d0
            else if (atomicnumber .eq. 9) then
               ! F radius
               x(L165-1+i) = 1.50d0 + 1.4d0
            else if (atomicnumber .eq. 35) then
               ! Br radius
               ! Bondi, J. Phys. Chem. 1964, 68, 441.
               x(L165-1+i) = 1.85d0 + 1.4d0
            else if (atomicnumber .eq. 20) then
               ! Ca radius
               ! Calculated from Aqvist, J. Phys. Chem. 1990, 94, 8021.
               x(L165-1+i) = 1.33d0 + 1.4d0
            else if (atomicnumber .eq. 11) then
               ! Na radius
               ! Calculated from Aqvist, J. Phys. Chem. 1990, 94, 8021.
               x(L165-1+i) = 1.87d0 + 1.4d0
            else if (atomicnumber .eq. 30) then
               ! Zn radius
               ! Hoops, Anderson, Merz, JACS 1991, 113, 8262.
               x(L165-1+i) = 1.10d0 + 1.4d0
            else if (atomicnumber .eq. 12) then
               !  Mg radius = 0.99A: ref. 21 in J. Chem. Phys. 1997, 107, 5422
               !  Mg radius = 1.18A: ref. 30 in J. Chem. Phys. 1997, 107, 5422
               !  Mg radius = 1.45A: Aqvist 1992
               x(L165-1+i) = 1.18d0 + 1.4d0
            else
               write( 0,* ) 'bad atom type: ',atype
               call mexit( 6,1 )
            end if

            ! dummy LCPO values:
            x(L170-1+i) = 0.0d0
            x(L175-1+i) = 0.0d0
            x(L180-1+i) = 0.0d0
            x(L185-1+i) = 0.0d0
            !  write(6,*) i,' ',atype,x(L165-1+i)
         end do  !  i=1,natom

      end if ! ( gbsa == 1 )
      
   end if  ! ( igb /= 0 .and. igb /= 10 .and. ipb == 0 )

   !-----------------------------------
   ! If a LRT calculation is requested, 
   ! setup the icosa-SASA parameters 
   ! this code is copied from above
   !-----------------------------------

   if ( ilrt /= 0 ) then
      !     --- assign parameters for calculating SASA according to the
      !         ICOSA method; the radii are augmented by 1.4 A ---
      do i=1,natom
         write(atype,'(a2)') ih(m06+i-1)
         if(ix(i100) .eq. 1) then
            atomicnumber = ix(i100+i)
         else
            call get_atomic_number(ih(m04+i-1), x(lmass+i-1), &
               atomicnumber, errFlag)
#ifdef LES
            if(errFlag) then
               call get_atomic_number(ih(m04+i-1), &
                  x(lmass+i-1)*lesfac((lestyp(i)-1)*nlesty+lestyp(i)), &
                  atomicnumber, errFlag)
            end if
#endif
         end if
         if (atomicnumber .eq. 7) then
            x(L165-1+i) = 1.55d0 + 1.4d0
         else if (atomicnumber .eq. 6) then
            x(L165-1+i) = 1.70d0 + 1.4d0
         else if (atomicnumber .eq. 1 .or. &
              ! added for lone pairs
              atype == 'EP') then
            x(L165-1+i) = 1.20d0 + 1.4d0
         else if (atomicnumber .eq. 8) then
            x(L165-1+i) = 1.50d0 + 1.4d0
         else if (atomicnumber .eq. 15) then
            x(L165-1+i) = 1.80d0 + 1.4d0
         else if (atomicnumber .eq. 16) then
            x(L165-1+i) = 1.80d0 + 1.4d0
         else if (atomicnumber .eq. 12) then
            !             Mg radius = 0.99A: ref. 21 in J. Chem. Phys. 1997, 107, 5422
            !             Mg radius = 1.18A: ref. 30 in J. Chem. Phys. 1997, 107, 5422
            !             Mg radius = 1.45A: Aqvist 1992
            x(L165-1+i) = 1.18d0 + 1.4d0
         else
            write( 0,* ) 'bad atom type: ',atype,' cannot perform SASA calculation'
            call mexit( 6,1 )
         end if  !  atype(1:1) == 'N'
         x(L170-1+i) = 0.0d0
         x(L175-1+i) = 0.0d0
         x(L180-1+i) = 0.0d0
         x(L185-1+i) = 0.0d0
         !write(6,*) i,' ',atype,x(L165-1+i)
      end do  !  i=1,natom
      
   end if ! ( ilrt /= 0 )
   
   !------------------------------------------------------------------------
   ! If user has requested Poisson-Boltzmann electrostatics, set up variables
   !------------------------------------------------------------------------
  
   if ( igb == 10 .or. ipb /= 0 ) then
      call pb_init(ifcap,natom,nres,ntypes,nbonh,nbona,ix(i02),ix(i04),ix(i06),ix(i08),ix(i10),&
                   ix(iibh),ix(ijbh),ix(iiba),ix(ijba),ix(ibellygp),ih(m02),ih(m04),ih(m06),x(l15),x(l97))
   end if  ! ( igb == 10 .or. ipb /= 0 ) 

   if (icnstph /= 0) then
      !  Initialize all constant pH data to 0 and read it in
      call cnstph_zero()
      call cnstphread(x(l15))

      !     Fill proposed charges array from current charges
      do i=1,natom
         x(l190-1+i) = x(l15-1+i)
      end do

      !  If we're doing explicit CpH, fill gbv* arrays
      if ( icnstph .gt. 1 .and. cph_igb == 2 .or. cph_igb == 5) then
        do i=1,natom
            x(l2402+i-1) = gbalpha
            x(l2403+i-1) = gbbeta
            x(l2404+i-1) = gbgamma
        end do
      end if
   end if

!  +---------------------------------------------------------------+
!  |  Read EVB input file                                          |
!  +---------------------------------------------------------------+

   if( ievb /= 0 ) then
#ifdef MPI
      call evb_input
      call evb_init
#  if defined(LES)
!KFW  call evb_pimd_init
#  endif
#else
      write(6,'(/2x,a)') 'Setting ievb>0 requires compilation with MPI'
      call mexit(6,1)
#endif
   endif

   if( iyammp /= 0 ) write( 6, '(a)' ) '  Using yammp non-bonded potential'
   
   ! -------------------------------------------------------------------
   !
   ! -- add check to see if the space in nmr.h is likely to be
   !     too small for this run:
   !     [Note: this check does *not* indicate if MXTAU and MXP are
   !      too small.  There is no easy way to ensure this, since
   !      the experimental intensities are read in a namelist
   !      command: if too many intensities are input, the read
   !      statment may cause a coredump before returning control
   !      to the main program.  Be careful.  sigh....]
   
   if (natom > matom .and. nmropt > 1) then
      write(6,*) 'WARNING: MATOM in nmr.h is smaller than the ', &
            natom,' atoms in this molecule.'
      write(6,*) 'Printout of NMR violations may be compromised.'
   end if
   
   ! -------------------------------------------------------------------
   !     --- checks on bogus data ---
   ! -------------------------------------------------------------------
   
   inerr = 0
      
   if( icfe < 0 .or. icfe > 1 ) then
      write(6,*) 'icfe must be 0 or 1 (icfe=2 is no longer supported)'
      inerr = 1
   end if
   if( icfe /= 0 .and. numgroup /= 2 ) then
      write(6,*) 'numgroup must be 2 if icfe is set'
      inerr = 1
   end if
   if (ievb>0) then
#ifdef MPI
!KFW  if( numgroup /= 2 ) then
!KFW     write(6,*) 'numgroup must be 2 if ievb is set'
!KFW     inerr = 1
!KFW  end if
#else
      write(6,'(/2x,a)') 'Setting ievb>0 requires compilation with MPI'
      inerr = 1
#endif
   end if
#ifdef PUPIL_SUPPORT
   ! BPR: PUPIL does not work with GB (or, I suppose, PB) for the
   ! time being. It is known to either crash or produce bogus
   ! results.
   if (igb > 0 .or. ipb /= 0) then
      write(6,'(a)') 'Cannot use implicit solvation (GB or PB) with PUPIL'
      inerr = 1
   end if
#endif /*PUPIL_SUPPORT*/
   if( (igb > 0 .or. ipb /= 0) .and. numextra > 0) then
      write(6,'(a)') 'Cannot use igb>0 with extra-point force fields'
      inerr = 1
   end if

!AMD validation
   if(iamd.gt.0)then
     if (EthreshD .eq. 0.d0 .and. alphaD .eq. 0.d0 .and. EthreshP .eq. 0.d0 .and. alphaP .eq. 0.d0) then
     write(6,'(a,i3)')'| AMD error all main parameters are 0.0 for Accelerated MD (AMD) or Windowed Accelerated MD (wAMD) '
     inerr = 1
     endif
    if(w_amd.gt.0)then
     if (EthreshD_w .eq. 0.d0 .and. alphaD_w .eq. 0.d0 .and. EthreshP_w .eq. 0.d0 .and. alphaP_w .eq. 0.d0) then
     write(6,'(a,i3)')'| AMD error all extra parameters are 0.0 for Windowed Accelerated MD (wAMD) LOWERING BARRIERS'
     inerr = 1
     endif
    endif
   endif
   if (iamd .gt. 0 .and. iamoeba ==1) then
      write(6,*)'amoeba is incompatible with AMD for now'
      inerr=1
   end if

   if (ips < 0 .or. ips > 6) then
      write(6,'(/2x,a,i3,a)') 'IPS (',ips,') must be between 0 and 6'
      inerr = 1
   end if
   if (ips /= 0 .and. ipol > 0 ) then
      write(6,'(/2x,a)') 'IPS and IPOL are inconsistent options'
      inerr = 1
   endif
   if (ips /= 0 .and. lj1264 /= 0) then
      write(6, '(/2x,a)') 'IPS and the LJ 12-6-4 potential are incompatible'
      inerr = 1
   endif
   if ( (igb > 0 .or. ipb /= 0) .and. ips > 0 ) then
      write(6,'(/2x,a,i3,a,i3,a)') 'IGB (',igb,') and ips (',ips, &
          ') cannot both be turned on'
      inerr = 1
   end if
   if (igb /= 0 .and. igb /= 1 .and. igb /= 2 .and. igb /= 5 &
         .and. igb /= 6 .and. igb /= 7 .and. igb /= 8 .and. igb /= 10) then
      write(6,'(/2x,a,i3,a)') 'IGB (',igb,') must be 0,1,2,5,6,7,8 or 10.'
      inerr = 1
   end if
   if (alpb /= 0 .and. alpb /= 1 )  then
      write(6,'(/2x,a,i3,a)') 'ALPB (',alpb,') must be 0 or 1.'
      inerr = 1
   end if
   if (alpb /= 0 .and. igb /= 1 .and. igb /= 2 .and. igb /= 5 .and. igb /=7 )  then
      write(6,'(/2x,a,i3,a)') 'IGB (',igb,') must be 1,2,5, or 7 if ALPB > 0.'
      inerr = 1
   end if

   if (jar < 0 .or. jar > 1) then
      write(6,'(/2x,a,i3,a)') 'JAR (',jar,') must be 0 or 1'
      inerr = 1
   end if

#ifdef LES
   if( igb /= 0 .and. igb /= 1 .and. igb /= 5 .and. igb /=7 ) then
      write(6,'(/,a)') 'Error: LES is only compatible with IGB > 0,1,5,7'
      inerr = 1
   end if
   if( alpb /= 0) then
      write(6,'(/,a)') 'Error: LES is not compatible with ALPB'
      inerr = 1
   end if
   if( gbsa > 0 ) then
      write(6,'(/,a)') 'Error: LES is not compatible with GBSA > 0'
      inerr = 1
   end if
   if( qmmm_nml%ifqnt ) then
      write(6,'(/,a)') 'Error: LES is not compatible with QM/MM'
      inerr = 1
   end if
   if( ipol > 0 ) then
      write(6,'(/,a)') 'Error: LES is not compatible with IPOL > 0'
      inerr = 1
   end if
   if (temp0les >= 0.d0 .and. iscale > 0 ) then
      write (6,'(/,a)') 'Error: iscale cannot be used with temp0les'
      inerr = 1
   end if
#endif
   if (irest /= 0 .and. irest /= 1) then
      write(6,'(/2x,a,i3,a)') 'IREST (',irest,') must be 0 or 1.'
      inerr = 1
   end if
   if (ibelly /= 0 .and. ibelly /= 1) then
      write(6,'(/2x,a,i3,a)') 'IBELLY (',ibelly,') must be 0 or 1.'
      inerr = 1
   end if
   if (imin < 0) then
      write(6,'(/2x,a,i3,a)') 'IMIN (',imin,') must be >= 0.'
      inerr = 1
   end if
   if (imin == 5) then
      if (ifbox /= 0 .and. ntb == 2) then
         write(6,'(/2x,a)') 'WARNING: IMIN=5 with changing periodic boundaries (NTB=2) can result in'
         write(6,'(/2x,a)') '         odd energies being calculated. Use with caution.'
      endif
#ifdef MPI
      if (sandersize > 1 .and. ntb == 2) then
         write(6,'(/2x,a)') 'ERROR: IMIN=5 and NTB=2 cannot be run with multiple processors.'
         inerr = 1
      endif
#endif
   end if


   if (iscale > mxvar) then
      write(6,9501) iscale,mxvar
      9501 format('ERROR: ISCALE (',i5,') exceeds MXVAR (',i5, &
            '). See nmr.h')
      inerr = 1
   end if
   if (ntx < 1 .or. ntx > 7) then
      write(6,'(/2x,a,i3,a)') 'NTX (',ntx,') must be in 1..7'
      inerr = 1
   end if
   
   if (ntb /= 0 .and. ntb /= 1 .and. ntb /= 2) then
      write(6,'(/2x,a,i3,a)') 'NTB (',ntb,') must be 0, 1 or 2.'
      inerr = 1
   end if
   if (ntb == 0 .and. iwrap > 0) then
      write(6,'(/2x,a)') 'Error: IWRAP > 0 cannot be used without a periodic box.'
      inerr = 1
   end if
   
   if (ntt < 0 .or. ntt > 8) then                                      ! APJ
      write(6,'(/2x,a,i3,a)') 'NTT (',ntt,') must be between 0 and 8.' ! APJ
      inerr = 1
   end if
   if (ntt == 1 .and. tautp < dt) then
      write(6, '(/2x,a,f6.2,a)') 'TAUTP (',tautp,') < DT (step size)'
      inerr = 1
   end if
   if( ntt < 3 .or. ntt > 8 ) then  ! APJ
      if( gamma_ln > 0.d0 ) then
         write(6,'(a)') 'ntt must be 3 to 8 if gamma_ln > 0' ! APJ
         inerr = 1
      end if
   end if

   if ( ntt==3 .or. ntt==6 ) nchain = 0 !APJ: Langevin, Adaptive-Langevin must have chain set to zero. ! APJ

   if (ntt == 3 .or. ntt == 4) then
      if ( ntb == 2) then
        !Require gamma_ln > 0.0d0 for ntt=3 and ntb=2 - strange things happen
        !if you run NPT with NTT=3 and gamma_ln = 0.
        if ( gamma_ln <= 0.0d0 ) then
          write(6,'(a)') 'gamma_ln must be > 0 for ntt=3 .or. 4 with ntb=2.'
          inerr = 1
        end if
      end if
   end if

   if (ntp /= 0 .and. ntp /= 1 .and. ntp /= 2 .and. ntp /= 3) then
      write(6,'(/2x,a,i3,a)') 'NTP (',ntp,') must be 0, 1, 2, or 3.'
      inerr = 1
   end if
   if (ntp == 3 .and. csurften < 1) then
      write(6,'(/2x,a)') 'csurften must be greater than 0 for ntp=3.'
      inerr = 1
   end if
   if (ntp > 0 .and. taup < dt .and. barostat == 1) then
      write(6, '(/2x,a,f6.2,a)') 'TAUP (',taup,') < DT (step size)'
      inerr = 1
   end if
   if (npscal < 0 .or. npscal > 1) then
      write(6,'(/2x,a,i3,a)') 'NPSCAL (',npscal,') must be 0 or 1.'
      inerr = 1
   end if
   if (ntp > 0) then
      if (barostat /= 1 .and. barostat /= 2) then
         write(6, '(/2x,a,i3,a)') 'BAROSTAT (', barostat, ') must be 1 or 2'
         inerr = 1
      end if
      if (barostat == 2) then
         if (mcbarint <= 0) then
            write(6, '(/2x,a,i3,a)') 'MCBARINT (',mcbarint,') must be positive'
            inerr = 1
         end if
         if (mcbarint >= nstlim) then
            write(6, '(a)') 'WARNING: mcbarint is greater than the number of &
                     &steps. This is effectively constant volume.'
         end if
      end if
   end if
   
   if (ntc < 1 .or. ntc > 4) then
      write(6,'(/2x,a,i3,a)') 'NTC (',ntc,') must be 1,2,3 or 4.'
      inerr = 1
   end if
   if (jfastw < 0 .or. jfastw > 4) then
      write(6,'(/2x,a,i3,a)') 'JFASTW (',jfastw,') must be 0->4.'
      inerr = 1
   end if
   
   if (ntf < 1 .or. ntf > 8) then
      write(6,'(/2x,a,i3,a)') 'NTF (',ntf,') must be in 1..8.'
      inerr = 1
   end if
   
   if (ioutfm /= 0 .and. ioutfm /= 1) then
      write(6,'(/2x,a,i3,a)') 'IOUTFM (',ioutfm,') must be 0 or 1.'
      inerr = 1
   end if
   
   if (ntpr < 0) then
      write(6,'(/2x,a,i3,a)') 'NTPR (',ntpr,') must be >= 0.'
      inerr = 1
   end if
   if (ntwx < 0) then
      write(6,'(/2x,a,i3,a)') 'NTWX (',ntwx,') must be >= 0.'
      inerr = 1
   end if
   if (ntwv < -1) then
      write(6,'(/2x,a,i3,a)') 'NTWV (',ntwv,') must be >= -1.'
      inerr = 1
   end if
   if (ntwf < -1) then
      write(6, '(/2x,a,i3,a)') 'NTWF (',ntwf,') must be >= -1.'
      inerr = 1
   end if
   if (ntwv == -1 .and. ioutfm /= 1) then
      write (6, '(/2x,a)') 'IOUTFM must be 1 for NTWV == -1.'
      inerr = 1
   end if
   if (ntwf == -1 .and. ioutfm /= 1) then
      write(6, '(/2x,a)') 'IOUTFM must be 1 for NTWF == -1.'
      inerr = 1
   end if
   if (ntwv == -1 .and. ntwx == 0) then
      write (6, '(/2x,a)') 'NTWX must be > 0 for NTWV == -1.'
      inerr = 1
   end if
   if (ntwe < 0) then
      write(6,'(/2x,a,i3,a)') 'NTWE (',ntwe,') must be >= 0.'
      inerr = 1
   end if
   if (ntave < 0) then
      write(6,'(/2x,a,i3,a)') 'NTAVE (',ntave,') must be >= 0.'
      inerr = 1
   end if
   if (ntr /= 0 .and. ntr /= 1) then
      write(6,'(/2x,a,i3,a)') 'NTR (',ntr,') must be 0 or 1.'
      inerr = 1
   end if
   if (ntrx /= 0 .and. ntrx /= 1) then
      write(6,'(/2x,a,i3,a)') 'NTRX (',ntrx,') must be 0 or 1.'
      inerr = 1
   end if
   if (nmropt < 0 .or. nmropt > 2) then
      write(6,'(/2x,a,i3,a)') 'NMROPT (',nmropt,') must be in 0..2.'
      inerr = 1
   end if
   
   if (idecomp < 0 .or. idecomp > 4) then
      write(6,'(/2x,a)') 'IDECOMP must be 0..4'
      inerr = 1
   end if
      
   ! check settings related to ivcap

   if(ivcap == 3 .or. ivcap == 4) then
      write(6,'(/2x,a)') 'IVCAP == 3 and IVCAP == 4 currently not implemented'
      inerr = 1
   endif
   if (ivcap < 0 .and. ivcap > 5) then
      write(6,'(/2x,a)') 'IVCAP must be 0 ... 5'
      inerr = 1
   end if
   if ((ivcap == 1 .or. ivcap == 5) .and. igb /= 10 .and. ipb == 0) then
      write(6,'(/2x,a)') 'IVCAP == 1,5 only works with Poisson Boltzmann (igb=10 or ipb/=0)'
      inerr = 1
   end if
   if((ivcap == 1 .or. ivcap == 3 .or. ivcap == 5 ) .and. cutcap <= 0.0d0) then
      write(6,'(/2x,a)') 'For IVCAP == 1,3, or 5, cutcap must be > 0.0'
      inerr = 1
   endif
   if (ivcap == 4 .and. &
         (xlorth < ZERO .or. ylorth < ZERO .or. zlorth < ZERO .or. &
!  give magic numbers a name  srb aug 2007 !
          xorth > 47114710.0d0 .or. &
!  give magic numbers a name !
          yorth > 47114710.0d0 .or. &
!  give magic numbers a name !
          zorth > 47114710.0d0)) then
      write(6,'(/2x,a)') &
      'For IVCAP == 4, xlorth, ylorth, zlorth, xorth, yorth, zorth must be set'
      inerr = 1
   end if
   if ((ivcap == 3 .or. ivcap == 4) .and. ibelly == 0) then
      write(6,'(/2x,a,a)') &
         'For IVCAP == 3 or 4, ibelly must be 1 and all atoms', &
         '  not in the spherical or orthorhombic region must be set NOT moving'
      inerr = 1
   end if
   if (ivcap == 5 .and. (imin /= 1 .or. maxcyc > 1)) then
      write(6,'(/2x,a,a)') &
         'IVCAP == 5 only works for single-point energy calculation'
      inerr = 1
   end if

   ! check if ifbox variable from prmtop file matches actual angles:

   if ( igb == 0 .and. ipb == 0 .and. ntb /= 0 ) then
      if ( ifbox == 1 ) then
         if ( abs(alpha - 90.0d0) > 1.d-5 .or. &
           abs(beta  - 90.0d0) > 1.d-5 .or. &
           abs(gamma - 90.0d0) > 1.d-5 ) then
           ifbox =3
           write(6,'(a)') '     Setting ifbox to 3 for non-orthogonal unit cell'
         end if
      end if

      if ( ifbox == 2 ) then
         if ( abs(alpha - 109.4712190d0) > 1.d-5 .or. &
              abs(beta  - 109.4712190d0) > 1.d-5 .or. &
              abs(gamma - 109.4712190d0) > 1.d-5 ) then
              write(6,'(/2x,a)') &
              'Error: ifbox=2 in prmtop but angles are not correct'
              inerr = 1
         end if
      end if
   end if

   ! checks for targeted MD
   if (itgtmd /= 0 .and. itgtmd > 2) then
      write(6,'(/2x,a,i3,a)') 'ITGTMD (',itgtmd,') must be 0, 1 or 2.'
      inerr = 1
   end if
   if (itgtmd == 1 .and. ntr == 1) then
      if (len_trim(tgtfitmask) > 0 .or. len_trim(tgtrmsmask) <= 0) then
         write(6,'(/2x,a)') 'ITGTMD: tgtrmsmask (and not tgtfitmask) ' //  &
                            'should be specified if NTR=1'
         inerr = 1
      end if
   end if
   ! skip this test until fallback to rgroup() is supported
   !if (itgtmd == 1 .and. ntr == 0) then
   !   if (len_trim(tgtfitmask) == 0 .and. len_trim(tgtrmsmask) == 0) then
   !      write(6,'(/2x,a)')  &
   !        'ITGTMD: both tgtfitmask and tgtrmsmask should be specified if NTR=0'
   !      inerr = 1
   !   end if
   !end if
   
   !     -- consistency checking
   
   if (imin > 0.and.nrespa > 1)  then
      write(6,'(/2x,a)') 'For minimization, set nrespa,nrespai=1'
      inerr = 1
   end if
   if (ntp > 0 .and. nrespa > 1) then
      write(6,'(/2x,a)') 'nrespa must be 1 if ntp>0'
      inerr = 1
   end if
   if  (ntx < 4.and.init /= 3)  then
      write(6,'(/2x,a)') 'NTX / IREST inconsistency'
      inerr = 1
   end if
   if (ntb == 2 .and. ntp == 0) then
      write(6,'(/2x,a)') 'NTB set but no NTP option (must be 1 or 2)'
      inerr = 1
   end if
   if (ntp /= 0 .and. ntb /= 2) then
      write(6,'(/,a,a)')' NTP > 0 but not constant pressure P.B.C.', &
            ' (NTB = 2) must be used'
      inerr = 1
   end if
   if (ntb /= 0 .and. ifbox == 0 .and. ntp /= 0) then
      write(6,'(/,a)') ' (NTB /= 0 and NTP /= 0) but IFBOX == 0'
      write(6,'(/,a)') ' This combination is not supported'
      inerr = 1
   end if
   if (ntb /= 0 .and. &
         ( box(1) < 1.d0  .or. &
         box(2) < 1.d0  .or. &
         box(3) < 1.d0 ) ) then
      write(6,'(/,a,3f10.3)') ' BOX is too small: ',box(1),box(2),box(3)
      inerr = 1
   else if (ntb /= 0 .and. &
         (sqrt(cut) >= box(1)*0.5d0 .or. &
         sqrt(cut) >= box(2)*0.5d0 .or. &
         sqrt(cut) >= box(3)*0.5d0) ) then
      write(6,'(/,a)') ' CUT must be < half smallest box dimension'
      inerr = 1
   end if
   if (ntb /= 0 .and. (igb > 0 .or. ipb /= 0) ) then
      write(6,'(/,a)') ' igb>0 is only compatible with ntb=0'
      inerr = 1
   end if
#ifdef APBS
   if ( ntb == 0 .and. sqrt(cut) < 8.05 .and. igb /= 10 .and. ipb == 0 .and. &
      .not. mdin_apbs) then
#else
   if ( ntb == 0 .and. sqrt(cut) < 8.05 .and. igb /= 10 .and. ipb == 0 ) then
#endif /* APBS */
      write(6,'(/,a,f8.2)') ' unreasonably small cut for non-periodic run: ', &
         sqrt(cut)
      inerr = 1
   end if
   if ( rgbmax < 5.d0*fsmax ) then
      write(6,'(/,a,f8.2)') ' rgbmax must be at least ', 5.d0*fsmax
      inerr = 1
   end if
   if (icfe /= 0 .and. indmeth == 3 ) then
      write(6,'(/,a)') ' indmeth=3 cannot be used with icfe>0'
      inerr = 1
   end if
   if (icfe /= 0 .and. ibelly /= 0 ) then
      write(6,'(/,a)') ' ibelly cannot be used with icfe'
      inerr = 1
   end if
#ifdef MPI /* SOFT CORE */
   if (icfe /= 0 .and. dvdl_norest /= 0 ) then
      write(6,'(/,a)') 'dvdl_norest must == 0!'
      write(6,'(/,a)') 'The dvdl_norest option is deprecated.'
      write(6,'(/,a)') 'Restraint energies are seperated, &
        &and do not contribute to dvdl.'
      inerr = 1
   end if
#endif
   ! Modification done by Ilyas Yildirim
   if (icfe == 1 .and. (klambda < 1 .or. klambda > 6)) then
     write(6,'(/,a)') ' klambda must be between 1 and 6'
     inerr = 1
   end if
   ! End of modification done by Ilyas Yildirim                                       

   if (clambda < 0.d0 .or. clambda > 1.d0 ) then
      write(6,'(/,a)') ' clambda must be between 0 and 1'
      inerr = 1
   end if

   if (icfe /= 0 .and. (idecomp == 3 .or. idecomp == 4)) then
      write(6,'(/,a)') ' Pairwise decomposition for thermodynamic integration not implemented'
      inerr = 1
   end if
   if (icfe /= 0 .and. idecomp /= 0 .and. ipol > 0) then
      write(6,'(/,a)') ' IPOL is incompatible with IDECOMP and ICFE'
      inerr = 1
   end if
 
#ifdef MPI /* SOFT CORE */
   if (ifsc /= 0) then
      if (ifsc == 2) then
         write (6,'(/,a)') 'The ifsc=2 option is no longer supported. &
               &Internal energies of the soft core region.'
         write (6,'(a)') 'are now handled implicitly and setting ifsc=2 &
               &is no longer needed'
         inerr = 1
      end if
      if (icfe /= 1 .and. ifsc==1) then
         write (6,'(/,a)') ' Softcore potential requires a standard TI run, set icfe to 1'
         inerr = 1
      end if
      if ( igb > 0 .or. ipb /= 0 ) then
         write (6,'(/,a)') ' Softcore potential is incompatible with GB (for now)'
         inerr = 1
      end if
      if ( ntf > 1 ) then
         write (6,'(/,a)') ' Softcore potentials require ntf=1 because SHAKE &
               &constraints on some bonds might be removed'
         inerr = 1
      end if
      if (clambda > 0.995 .or. clambda < 0.005) then
         write (6,'(/,a)') ' Softcore potentials cannot be used with clambda < 0.005 or > 0.995'
         inerr = 1
      end if
      if (klambda /= 1) then
         write (6,'(/,a)') ' Softcore potential requires linear mixing, set klambda to 1'
         inerr = 1
      end if
      if (imin == 1 .and. ntmin /= 2) then
         write (6,'(/,a)') ' Minimizations with ifsc=1 require the steepest descent algorithm.'
         write (6,'(/,a)') ' Set ntmin to 2 and restart'
         inerr = 1
      end if
   end if
   if (ifmbar /= 0) then
      if (icfe /= 1) then
         write (6,'(/,a)') ' MBAR requires a standard TI run, set icfe to 1'
         inerr = 1
      end if
      if (ifsc /=0 .and. (bar_l_max > 0.995 .or. bar_l_min < 0.005) ) then
         write (6,'(/,a)') ' Softcore potentials cannot be used with &
               &bar_l_min < 0.005 or bar_l_max > 0.995'
         inerr = 1
      end if
      if (klambda /= 1) then
         write (6,'(/,a)') ' MBAR requires linear mixing, set klambda to 1'
         inerr = 1
      end if
   end if
#endif

   if (ilrt /= 0 .and. lrtmask == '') then
      write (6,'(a)') 'Linear Response Theory activated, but lrtmask is not set'
      inerr = 1
   end if
  
   if (idecomp > 0 .and. (ntr > 0 .or. ibelly > 0)) then
      write(6,'(/,a)') 'IDECOMP is not compatible with NTR or IBELLY'
      inerr = 1
   end if
   if (icnstph /= 0) then

      if ( icnstph < 0 ) then
         write(6, '(/,a)') 'icnstph must be greater than 0'
         inerr = 1
      end if
      if ( igb == 0 .and. ipb == 0 .and. icnstph == 1 ) then
         write(6, '(/,a)') 'Constant pH using icnstph = 1 requires &
                           &GB implicit solvent'
         inerr = 1
      end if
      if ( ntb .eq. 0 .and. icnstph .gt. 1 ) then
         write(6, '(/,a)') 'Constant pH using icnstph = 2 requires &
                           &periodic boundary conditions'
         inerr = 1
      end if
      if (icfe /= 0) then
         write(6, '(/,a)') &
         'Constant pH and thermodynamic integration are incompatable'
         inerr = 1
      end if

      if (ntcnstph <= 0) then
         write(6, '(/,a)') 'ntcnstph must be a positive integer.'
         inerr = 1
      end if

      if (icnstph > 1 .and. mccycles <= 0) then
         write(6, '(/,a)') 'mccycles must be a positive integer.'
         inerr = 1
      end if

   end if ! icnstph

#ifdef noVIRIAL
   if( ntp > 0 .and. barostat == 1 ) then
      write(6,'(/,a)') 'Error: Berendsen barostat is incompatible with noVIRIAL'
      inerr = 1
   end if
#endif
   
   !-----------------------------------------------------
   !     ----sanity checks for Ewald
   !-----------------------------------------------------
   
   if( igb == 0 .and. ipb == 0 ) then
      call float_legal_range('skinnb: (nonbond list skin) ', &
            skinnb,skinlo,skinhi)
      
      !  --- Will check on sanity of settings after coords are read in
      !      and the extent of the system is determined.
      
      if(periodic == 1)then
         call float_legal_range('skinnb+cutoffnb: (nonbond list cut) ', &
               skinnb+cutoffnb,zero,sphere)
      end if
      if (ntb==0 .and. use_pme/=0) then
         write(6,'(/,a)') &
         'Using PME with a non-periodic simulation does not make sense. Set either ntb>0 of use_pme=0.'
         inerr = 1
      end if
      call float_legal_range('a: (unit cell size) ',a,boxlo,boxhi)
      call float_legal_range('b: (unit cell size) ',b,boxlo,boxhi)
      call float_legal_range('c: (unit cell size) ',c,boxlo,boxhi)
      call float_legal_range('alpha: (unit cell angle) ', &
            alpha,anglo,anghi)
      call float_legal_range('beta: (unit cell angle) ', &
            beta,anglo,anghi)
      call float_legal_range('gamma: (unit cell angle) ', &
            gamma,anglo,anghi)
      call int_legal_range('order: (interpolation order) ', &
            order,orderlo,orderhi)
      call opt_legal_range('verbose: ',verbose,0,4)
      call opt_legal_range('netfrc: ',netfrc,0,1)
      call opt_legal_range('nbflag: ',nbflag,0,1)
      call opt_legal_range('nbtell: ',nbtell,0,2)
      call opt_legal_range('ew_type: ',ew_type,0,1)
      call opt_legal_range('vdwmeth: ',vdwmeth,0,2)
      call opt_legal_range('eedmeth: ',eedmeth,1,6)
      call opt_legal_range('ee_type: ',ee_type,1,2)
      call opt_legal_range('maxiter: ',maxiter,1,50)
      call opt_legal_range('indmeth: ',indmeth,0,3)
      call opt_legal_range('fix_quad: ',fix_quad,0,1)
      call float_legal_range('eedtbdns: (erfc table density) ', &
            eedtbdns,denslo,denshi)
   end if  ! ( igb == 0 .and. ipb == 0 )

   if( ntb==2 .and. ipimd==1) then
      write(6,*) 'primitive PIMD is incompatible with NTP ensemble'
      inerr=1
   endif

   if( ntb==2 .and. ipimd==3 ) then
      write(6,*) 'CMD is incompatible with NTP ensemble'
      inerr=1
   endif

   if( ipimd==3 .and. adiab_param>=1.0 ) then
      write(6,*) 'For CMD adiab_param must be <=1'
      inerr=1
   endif

   if( ntb==2 .and. ipimd==4 ) then
      write(6,*) 'RPMD is incompatible with NTP ensemble'
      inerr=1
   endif

   if( ntt/=0 .and. ipimd==4 ) then
      write(6,*) 'RPMD is incompatible with NVT ensemble'
      inerr=1
   endif

   if( ntt/=4 .and. ipimd==2 ) then
      write(6,*) 'NMPIMD requires Nose-Hoover chains (ntt=4)'
      inerr=1
   endif

   if( ntt/=4 .and. ipimd==3 ) then
      write(6,*) 'CMD requires Nose-Hoover chains (ntt=4)'
      inerr=1
   endif

#if !defined(MPI)
   if( ineb > 0 ) then
      write(6,*) 'NEB requires MPI'
      inerr=1
   endif
#endif

#ifdef LES
   if( ineb > 0 ) then
      write(6,*) 'NEB no longer works with LES: use multiple groups and a groupfile instead'
      inerr=1
   endif
#endif

   if( ineb > 0 .and. ipimd > 0 ) then
      write(6,*) 'ineb>0 and ipimd>0 are incompatible options'
      inerr=1
   endif

   if( iamoeba == 1 )then
#ifdef LES
      write(6,*)'amoeba is incompatible with LES'
      inerr=1
#endif

      if( ntc > 1 ) then
         write(6,*) 'SHAKE (ntc>1) and amoeba are incompatible options'
         inerr=1
      end if
      if( ntp > 1 .and. beeman_integrator > 0 ) then
         write(6,*) 'ntp>1 is not consistent with the beeman integrator'
         inerr=1
      end if
      if ( igb /= 0 ) then
         write (6,'(/,a)') 'amoeba (iamoeba=1) is incompatible with GB (igb>0)'
         inerr=1
      end if
      if ( ipb /= 0 ) then
         write (6,'(/,a)') 'amoeba (iamoeba=1) is incompatible with PB (ipb>0)'
         inerr=1
      end if
#ifdef MPI
      if (numtasks > 1) then
         write(6, '(/,a)') 'amoeba (iamoeba=1) cannot be run in parallel in &
                           &sander'
         inerr=1
      end if
#endif

   end if

   ! ---WARNINGS:
   
   if ( ibelly == 1 .and. igb == 0 .and. ipb == 0 .and. ntb /= 0 ) then
      write(6,'(/,a,/,a,/,a)') &
            'Warning: Although EWALD will work with belly', &
            '(for equilibration), it is not strictly correct!'
   end if
   
   if (inerr == 1) then
      write(6, '(/,a)') ' *** input error(s)'
      call mexit(6,1)
   end if
   
   ! Load the restrained atoms (ntr=1) or the belly atoms (ibelly=1)
   ! or atoms for targeted md (itgtmd=1). Selections are read from
   ! &cntrl variables or, if these are not defined, it falls back to
   ! the old group input format.
   
   if(mtmd /= 'mtmd') then
     itgtmd = 2
     ntr = 0
     emtmd = 0.0d0
   end if
   konst = ntr > 0
   dotgtmd = itgtmd > 0
   belly = .false.
   natc = 0
   ngrp = 0
   natbel = 0
   nattgtfit = 0  ! number of atoms for tgtmd fitting (=overlap)
   nattgtrms = 0  ! number of atoms for tgtmd rmsd calculation
   nrc = 0
   if(konst.or.dotgtmd) then
      
      ! inserted here to fix the bug that coords are not available
      ! yet when distance based selection (<,>) is requested
#ifdef LES
      call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0les,.FALSE.,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
         if( irest == 1 .and. beeman_integrator > 0 ) then
            write(6,*) 'Cannot do a beeman_integrator restart with old-style coordinates'
            call mexit(6,1)
         end if
         call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.FALSE.,solvph)
      endif
#endif
      
      if(itgtmd == 2) then
        call mtmdcall(emtmd,x(lmtmd01),ix(imtmd02),x(lcrd),x(lforce),ih(m04),ih(m02),ix(i02),&
                    ih(m06),x(lmass),natom,nres,'READ')  
      else
        ! DRR - Open and close calls are now in rdrest
        !if (ntrx <= 0) then
        !  call amopen(10,refc,'O','U','R')
        !else
        !  call amopen(10,refc,'O','F','R')
        !end if
        ! these messages should be written after "5.  REFERENCE..." ?
        if (konst) write(6,9408)
        if (dotgtmd) write(6,9409)
  
        call rdrest(natom,ntrx,refc,x(lcrdr))
        !close(10)



      
        ! VH - tgtmd change: preferably call atommask() instead of rgroup()
        if (konst) then
          if( len_trim(restraintmask) <= 0 ) then
              call rgroup(natom,natc,nres,ngrp,ix(i02),ih(m02), &
                    ih(m04),ih(m06),ih(m08),ix(icnstrgp),jgroup,indx,irespw,npdec, &
                    x(l60),x(lcrdr),konst,dotgtmd,belly,idecomp,5,.true.)
          else
              call atommask( natom, nres, 0, ih(m04), ih(m06), &
                ix(i02), ih(m02), x(lcrd), restraintmask, ix(icnstrgp) )
  
              ! for now, emulate the "GATHER ALL THE CONSTRAINED ATOMS TOGETHER"
              ! section of rgroup(); later, the various masks should be done
              ! differently, i.e. without the "gather", as in the following:
              !     x(l60:l60+natom-1) = restraint_wt
              !     natc = sum(ix(icnstrgp:icnstrgp+natom-1))
  
              natc = 0
              do i=1,natom
                if( ix(icnstrgp-1+i) <= 0 ) cycle
                natc = natc + 1
                ix(icnstrgp-1+natc) = i
                x(l60-1+natc) = restraint_wt
              end do
              write(6,'(a,a,a,i5,a)') '     Mask ', &
              restraintmask(1:len_trim(restraintmask)), ' matches ',natc,' atoms'
          end if
        end if
        nrc = natc
        
        if (itgtmd == 1) then
          if (len_trim(tgtfitmask) <= 0 .and. len_trim(tgtrmsmask) <= 0) then
              ! the following if-endif can be deleted when we stop
              ! supporting rgroup()
              if (konst) then
                ! cannot do both ntr and tgtmd together using old group format
                write(6,'(/2x,a)') 'NTR must be 0 for targeted MD (TGTMD=1)'
                call mexit(6,1)
              else  ! the following only for backward compatibility
                call rgroup(natom,natc,nres,ngrp,ix(i02),ih(m02), &
                      ih(m04),ih(m06),ih(m08),ix(icnstrgp), &
                      jgroup,indx,irespw,npdec, &
                      x(l60),x(lcrdr),konst,dotgtmd,belly,idecomp,5,.true.)
                ! tgtmd atoms are now stored in nattgt, igroup -> icnstrgp
                nattgtfit = natc
                nattgtrms = natc
                do i=1,nattgtfit
                    ix(itgtfitgp-1+i) = ix(icnstrgp-1+i)
                    ix(itgtrmsgp-1+i) = ix(icnstrgp-1+i)
                end do
              end if
          else
              if (ntr == 0) then  ! read tgtfitmask only if ntr=1 
                ! read in atom group for tgtmd fitting (=overlap region)
                call atommask( natom, nres, 0, ih(m04), ih(m06), &
                    ix(i02), ih(m02), x(lcrd), tgtfitmask, ix(itgtfitgp) )
                ! see comments above (for ntr) for the following reduction cycle
                nattgtfit = 0
                do i=1,natom
                  if( ix(itgtfitgp-1+i) <= 0 ) cycle
                  nattgtfit = nattgtfit + 1
                  ix(itgtfitgp-1+nattgtfit) = i
                end do
                write(6,'(a,a,a,i5,a)')  &
                '     Mask "', tgtfitmask(1:len_trim(tgtfitmask)-1),  &
                '" matches ',nattgtfit,' atoms'
              end if
              ! read in atom group for tgtmd rmsd calculation
              call atommask( natom, nres, 0, ih(m04), ih(m06), &
                ix(i02), ih(m02), x(lcrd), tgtrmsmask, ix(itgtrmsgp) )
              nattgtrms = 0
              do i=1,natom
                if( ix(itgtrmsgp-1+i) <= 0 ) cycle
                nattgtrms = nattgtrms + 1
                ix(itgtrmsgp-1+nattgtrms) = i
              end do
              write(6,'(a,a,a,i5,a)')  &
              '     Mask "', tgtrmsmask(1:len_trim(tgtrmsmask)-1),  &
              '" matches ',nattgtrms,' atoms'
          end if
        end if
      
      end if ! (itgtmd == 2)

   end if  ! (konst.or.dotgtmd)

   if (ineb>0) then
      ! carlos: read in fitmask and rmsmask info for NEB, just as done for tgtmd
      ! init last_neb_atom, which is used to determine the limits for the
      ! communication of neighbor coordinates (to reduce size for explicit
      ! water)
      last_neb_atom = 0

      if (ntr /= 0) then
        write(6,'(/2x,a)') 'cannot use NEB with ntr restraints'
! CARLOS: WHY NOT? SHOULD BE OK.
! potential for user error is restrained region overlaps with NEB region. would
! blow up.
        call mexit(6,1)
      else
         ! read in atom group for fitting (=overlap region)
         call atommask( natom, nres, 0, ih(m04), ih(m06), &
            ix(i02), ih(m02), x(lcrd), tgtfitmask, ix(itgtfitgp) )
         ! see comments above (for ntr) for the following reduction cycle
         nattgtfit = 0
         do i=1,natom
           if( ix(itgtfitgp-1+i) <= 0 ) cycle
           nattgtfit = nattgtfit + 1
           ix(itgtfitgp-1+nattgtfit) = i
           if (i.gt.last_neb_atom) last_neb_atom = i
         end do
         write(6,'(a)') "The following selection will be used for NEB structure fitting"
         write(6,'(a,a,a,i5,a)')  &
         '     Mask "', tgtfitmask(1:len_trim(tgtfitmask)-1),  &
         '" matches ',nattgtfit,' atoms'
      end if
      ! read in atom group for tgtmd rmsd calculation
      call atommask( natom, nres, 0, ih(m04), ih(m06), &
         ix(i02), ih(m02), x(lcrd), tgtrmsmask, ix(itgtrmsgp) )
      nattgtrms = 0
      do i=1,natom
        if( ix(itgtrmsgp-1+i) <= 0 ) cycle
        nattgtrms = nattgtrms + 1
        ix(itgtrmsgp-1+nattgtrms) = i
        if (i.gt.last_neb_atom) last_neb_atom = i
      end do
      write(6,'(a)') "The following selection will be used for NEB force application"
      write(6,'(a,a,a,i5,a)')  &
      '     Mask "', tgtrmsmask(1:len_trim(tgtrmsmask)-1),  &
      '" matches ',nattgtrms,' atoms'
      write(6,'(/2x,a,i6)') "Last atom in NEB fitmask or rmsmask is ",last_neb_atom

      if (nattgtrms<=0 .or. nattgtfit <= 0) then
        write(6,'(/2x,a)') 'NEB requires use of tgtfitmask and tgtrmsmask'
        call mexit(6,1)
      endif
   endif
   
   ! dotgtmd may be false here even if doing tgtmd
   ! this is so belly info is read properly? following existing KONST code
   
   dotgtmd=.false.
   konst = .false.
   belly = ibelly > 0
   ngrp = 0
   if(belly) then
      ! inserted here to fix the bug that coords are not available
      ! yet when distance based selection (<,>) is requested
#ifdef LES
      call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0les,.FALSE.,solvph)
#else
      call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.FALSE.,solvph)
#endif
      write(6,9418)
      if( len_trim(bellymask) <= 0 ) then
         call rgroup(natom,natbel,nres,ngrp,ix(i02),ih(m02), &
            ih(m04),ih(m06),ih(m08),ix(ibellygp), &
            jgroup,indx,irespw,npdec, &
            x(l60),x(lcrdr),konst,dotgtmd,belly,idecomp,5,.true.)
      else
         call atommask( natom, nres, 0, ih(m04), ih(m06), &
            ix(i02), ih(m02), x(lcrd), bellymask, ix(ibellygp) )
         natbel = sum(ix(ibellygp:ibellygp+natom-1))
         write(6,'(a,a,a,i5,a)') '     Mask ', &
            bellymask(1:len_trim(bellymask)), ' matches ',natbel,' atoms'
      end if
   end if
   call setvar(ix,belly)

   !  see if the user has input a noshakemask string, and process it:
   natnos = 0
   if( len_trim(noshakemask) > 0 ) then
      call atommask( natom, nres, 0, ih(m04), ih(m06), &
         ix(i02), ih(m02), x(lcrd), noshakemask, noshakegp )
      natnos = sum(noshakegp(1:natom))
      write(6,*)
      write(6,'(a,a,a,i5,a)') 'Noshake mask ', &
         noshakemask(1:len_trim(noshakemask)), ' matches ',natnos,' atoms'
      call setnoshake(ix,noshakegp,ntc,num_noshake)
      if( ntf > 1 ) then
         write(6,'(a)') '   Setting ntf to 1'
         ntf = 1
      end if
   end if
   
   ! GMS ------------------------------------
   !  Check for 'iwrap_mask', and process it.
   ! ----------------------------------------
   n_iwrap_mask_atoms = 0
   if( len_trim(iwrap_mask) > 0 .and. iwrap == 2) then
      call atommask( natom, nres, 0, ih(m04), ih(m06), &
         ix(i02), ih(m02), x(lcrd), iwrap_mask, iwrap_maskgp )
      ! iwrap_maskgp is a natom long integer array, with elements:
      ! 0 --> atom is not in iwrap_mask
      ! 1 --> atom is in iwrap_mask
      n_iwrap_mask_atoms = sum(iwrap_maskgp(1:natom))
      write(6,*)
      write(6,'(a,a,a,i5,a)') 'Wrap mask ', &
         iwrap_mask(1:len_trim(iwrap_mask)), ' matches ',n_iwrap_mask_atoms,' atoms:'
      ! Set an array to store the atom numbers of the atoms
      ! in the iwrap_mask
      allocate(iwrap_mask_atoms(n_iwrap_mask_atoms), stat=ier)
      REQUIRE(ier == 0)

      j = 0
      do i=1,natom
        if( iwrap_maskgp(i)>0 ) then
          j = j+1
          iwrap_mask_atoms(j) = i
        end if
      end do
      
      write(6,'(10i5)') (iwrap_mask_atoms(i),i=1,n_iwrap_mask_atoms)
   end if

#ifdef MPI /* SOFT CORE */
   ! lower charges if a crgmask is set
   if ( len_trim(crgmask) > 0 ) then
      call atommask( natom, nres, 0, ih(m04), ih(m06), &
         ix(i02), ih(m02), x(lcrd), crgmask, crggp )
      write(6,'(a,a,a,i5,a)') 'Zero-Charge Mask ',crgmask(1:len_trim(crgmask)), ' matches ',sum(crggp(1:natom)),' atoms'
      call remove_charges(crggp, natom, x(l15))
   end if
#endif

   konst = .false.
   belly = .false.
   if(idecomp > 0) then
      write(6,9428)
      call rgroup(natom,ntmp,nres,ngrp,ix(i02),ih(m02), &
            ih(m04),ih(m06),ih(m08),ix(ibellygp), &
            jgroup,indx,irespw,npdec, &
            x(l60),x(lcrdr),konst,dotgtmd,belly,idecomp,5,.true.)
   end if
   
   if( ibelly > 0 .and. (igb > 0 .or. ipb /= 0) ) then
      
      !          ---here, the only allowable belly has just the first
      !             NATBEL atoms in the moving part.  Check to see that this
      !             requirement is satisfied:
      
      do i=natbel+1,natom
         if( ix(ibellygp+i-1) /= 0 ) then
            write(6,*) 'When igb>0, the moving part must be at the'
            write(6,*) '   start of the molecule.  This does not seem'
            write(6,*) '   to be the case here.'
            write(6,*) 'natbel,i,igroup(i) = ' &
                  ,natbel,i,ix(ibellygp+i-1)
            call mexit(6,1)
         end if
      end do
   end if
   

   !     ----- CALCULATE THE SQUARE OF THE BOND PARAMETERS FOR SHAKE
   !           THE PARAMETERS ARE PUT SEQUENTIALLY IN THE ARRAY CONP -----
   
   do i=1,nbonh + nbona + nbper
      j = ix(iicbh+i-1)
      x(l50+i-1) = req(j)**2
   end do
   
#ifdef MPI
      if( icfe /= 0 ) then

      !  use the masses of the prmtop file for the first group for both groups:
      !  [only the master nodes communicate here, since non-master nodes
      !   have not yet allocated space]
      ! This leads to problems for dual topology runs, and is therefore skipped
      ! if ifsc is set to one, the masses from both prmtop files are used
         if (ifsc == 0) then
            call mpi_bcast(x(lmass),natom,MPI_DOUBLE_PRECISION,0,commmaster,ierr)
            call mpi_bcast(x(lwinv),natom,MPI_DOUBLE_PRECISION,0,commmaster,ierr)            
            call mpi_bcast(x(l75),natom,MPI_DOUBLE_PRECISION,0,commmaster,ierr)
         end if
         tmass = sum(x(lmass:lmass+natom-1))
         tmassinv = 1.d0/tmass

      !  next, do a minimal sanity check that the SHAKE parameters are
      !  consistent on the two processors:

         ! For Softcore this might be allowed
         ! Put a better check here later
         if( ntc == 2 .and. ifsc == 0) then
            partner = ieor(masterrank,1)
            call mpi_sendrecv( nbonh, 1, MPI_INTEGER, partner, 5, &
                               nbonh_c, 1, MPI_INTEGER, partner, 5, &
                               commmaster, ist, ierr )
            call mpi_sendrecv( num_noshake, 1, MPI_INTEGER, partner, 5, &
                               num_noshake_c, 1, MPI_INTEGER, partner, 5, &
                               commmaster, ist, ierr )
            if (qmmm_nml%ifqnt .and. qmmm_nml%qmshake == 0) then
               ! qtw - if qmshake=0, we need to check the QM atoms
               call mpi_sendrecv( &
                       qmmm_struct%nquant, 1, MPI_INTEGER, partner, 5, &
                       nquant_c, 1, MPI_INTEGER, partner, 5, &
                       commmaster, ist, ierr)
               call mpi_sendrecv( &
                       qmmm_struct%noshake_overlap, 1, MPI_INTEGER,partner, 5, &
                       noshake_overlap_c, 1, MPI_INTEGER, partner, 5, &
                       commmaster, ist, ierr)
               if ( (qmmm_struct%nquant-qmmm_struct%noshake_overlap) /= &
                    (nquant_c-noshake_overlap_c) ) then
                  call sander_bomb('mdread2', &
                       'QMMM: NOSHAKE lists are not match in two groups!', &
                       'try noshakemask in cntrl to match the noshake list')
               end if
            else if( nbonh - num_noshake /= nbonh_c - num_noshake_c ) then
               write(6,*) 'SHAKE lists are not compatible in the two groups!'
               call mexit(6,1)
            end if
         else if( ntc == 3 ) then
            write(6,*) 'ntc = 3 is not compatible with icfe>0'
            call mexit(6,1)
         end if

      end if
#endif

   if ( iamoeba /= 1 )then
      if( igb == 0 .and. ipb == 0 ) &
         call init_extra_pts( &
         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), &
         ih(m06),ix,x,ix(i08),ix(i10), &
         nspm,ix(i70),x(l75),tmass,tmassinv,x(lmass),x(lwinv),req)
   endif
  
   !  DEBUG input; force checking
   call load_debug(5)

   return
   ! -------------------------------------------------------------------------
   ! Standard format statements:
   
   9328 format(/80('-')/,'   2.  CONTROL  DATA  FOR  THE  RUN',/80('-')/)
   9408 format(/4x,'LOADING THE CONSTRAINED ATOMS AS GROUPS',/)
   9409 format(/4x,'LOADING THE TARGETED MD ATOMS AS GROUPS',/)
   9418 format(/4x,'LOADING THE BELLY ATOMS AS GROUPS',/)
   9428 format(/4x,'LOADING THE DECOMP ATOMS AS GROUPS',/)
   9008 format(a80)
end subroutine mdread2 


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Emit defined preprocessor names, ie, flags.
subroutine printflags()

   implicit none
   integer     max_line_length
   parameter ( max_line_length = 80 )

   character(len=max_line_length) line  ! output string of active flags
   integer n                            ! len(line)
   
   line = '| Flags:'
   n = 8
   
#ifdef ISTAR2
   call printflags2(' ISTAR2',7,n,line,.false.)
#endif
#ifdef MPI
   call printflags2(' MPI',4,n,line,.false.)
# ifdef USE_MPI_IN_PLACE
   call printflags2(' USE_MPI_IN_PLACE',17,n,line,.false.)
# endif
#endif
#ifdef LES
   call printflags2(' LES',4,n,line,.false.)
#endif
#ifdef NMODE
   call printflags2(' NMODE',6,n,line,.false.)
#endif
#ifdef HAS_10_12
   call printflags2(' HAS_10_12',10,n,line,.false.)
#endif
#ifdef DNA_SHIFT
   call printflags2(' DNA_SHIFT',10,n,line,.false.)
#endif
#ifdef MMTSB
   call printflags2(' MMTSB',6,n,line,.false.)
#endif

#ifdef noVIRIAL
   call printflags2(' noVIRIAL',9,n,line,.false.)
#endif
   
   call printflags2(' ',1,n,line,.true.)
   return
end subroutine printflags 


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Primitive pre-Fortran90 implementation of printflags.
subroutine printflags2(flag,flag_len,line_len,line,last)

   implicit none
   integer     max_line_length
   parameter ( max_line_length = 80 )

   character(*) flag                ! flag name with blank prefix, intent(in)
   integer flag_len                 ! len(flag), intent(in)
   integer line_len                 ! len(line), intent(inout)
   character(len=max_line_length) line ! intent(inout)
   logical last                     ! is this the last flag ?, intent(in)

   if (line_len + flag_len > max_line_length) then
      write( 6,'(a)') line
      ! begin another line
      line = '| Flags:'
      line_len=8
   end if
   line=line(1:line_len) // flag(1:flag_len)
   line_len=line_len+flag_len
   if(last)write( 6,'(a)') line
   return
end subroutine printflags2 

!-------------------------------------------------
!     --- FLOAT_LEGAL_RANGE ---
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Check the range of a float; abort on illegal values.
subroutine float_legal_range(string,param,lo,hi)
   implicit none
   _REAL_ param,lo,hi
   character(len=*)string

   if ( param < lo .or. param > hi )then
      write(6,59)
      write(6,60)string,param
      write(6,61)
      write(6,62)lo,hi
      write(6,63)
      call mexit(6,1)
   end if
   59 format(/,1x,'Ewald PARAMETER RANGE CHECKING: ')
   60 format(1x,'parameter ',a,' has value ',e12.5)
   61 format(1x,'This is outside the legal range')
   62 format(1x,'Lower limit: ',e12.5,' Upper limit: ',e12.5)
   63 format(1x,'Check ew_legal.h')
   return
end subroutine float_legal_range 

!-------------------------------------------------
!     --- INT_LEGAL_RANGE ---
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Check the range of an integer; abort on illegal values.
subroutine int_legal_range(string,param,lo,hi)
   implicit none
   integer param,lo,hi
   character(len=*)string

   if ( param < lo .or. param > hi )then
      write(6,59)
      write(6,60)string,param
      write(6,61)
      write(6,62)lo,hi
      write(6,63)
      call mexit(6,1)
   end if
   59 format(/,1x,'PARAMETER RANGE CHECKING: ')
   60 format(1x,'parameter ',a,' has value ',i8)
   61 format(1x,'This is outside the legal range')
   62 format(1x,'Lower limit: ',i8,' Upper limit: ',i8)
   63 format(1x,'The limits may be adjustable; search in the .h files ')
   return
end subroutine int_legal_range 

!-------------------------------------------------
!     --- OPT_LEGAL_RANGE ---
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Check the range of an integer option; abort on illegal values.
subroutine opt_legal_range(string,param,lo,hi)
   implicit none
   integer param,lo,hi
   character(len=*)string

   if ( param < lo .or. param > hi )then
      write(6,59)
      write(6,60)string,param
      write(6,61)
      write(6,62)lo,hi
      write(6,63)
      call mexit(6,1)
   end if
   59 format(/,1x,'Ewald OPTION CHECKING: ')
   60 format(1x,'option ',a,' has value ',i5)
   61 format(1x,'This is outside the legal range')
   62 format(1x,'Lower limit: ',i5,' Upper limit: ',i5)
   63 format(1x,'Check the manual')
   return
end subroutine opt_legal_range 
#endif /*ifndef PBSA*/

!-------------------------------------------------
!     --- SANDER_BOMB ---
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Print an error message and quit
subroutine sander_bomb(routine,string1,string2)
   implicit none
   character(len=*) routine,string1,string2

   write(6, '(1x,2a)') &
         'SANDER BOMB in subroutine ', routine
   write(6, '(1x,a)') string1
   write(6, '(1x,a)') string2
   call mexit(6,1)
end subroutine sander_bomb
!-------------------------------------------------

!-------------------------------------------------
!     --- remove_charges ---
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Zero charges on some atoms
subroutine remove_charges(crggp,natom,charge)
  use constants, only: INV_AMBER_ELECTROSTATIC
  implicit none
  integer natom, crggp(*),i
  _REAL_ charge(*), charge_removed
  
  charge_removed = 0.d0
  do i=1,natom
     if (crggp(i)==1) then
        charge_removed = charge_removed + charge(i) * INV_AMBER_ELECTROSTATIC
        write (6,'(a,f12.4,a,i5)') 'Removing charge of ', charge(i) * INV_AMBER_ELECTROSTATIC,' from atom ',i
        charge(i)=0
     end if
  end do
  write(6, '(a,f12.4,a)') 'Total charge of ',charge_removed,' removed'
  RETURN
end subroutine remove_charges
!-------------------------------------------------
