subroutine nip
use Numeric_Kinds_Module
use CommonInfo_Module
USE Convrsns_Module
use APH_Module
use NIP_Module
use Masses_Module
implicit none

!======================================================================================
! Calculate the magnitude of a negative imaginary potential (NIP)
!       Reference: Bear & Neuhasuer,  J. Chem. Phys. 90, 4351 (1989)
!
! Note that this routine calculates a real number for the NIP.
! Need to multiply by -i before applying to wave function.
!
! Varaibles:
!   order = order of the NIP function
!   rmax  = APH rho value of outer NIP boundary (the rho grid outer boundary)
!   rmin  = APH rho value of inner NIP boundary
!   ufac  = amplitude of the NIP function
!   vnip  = array of NIP values
!
!=========================================================================================

!=========================================================================================
!     I N T E R N A L S
INTEGER  :: irho
REAL(dp) :: rmin,rmax
!=========================================================================================

 ALLOCATE(vnip(nrho))

 rmin=rho_val(ind_nip_rho)
 rmax=rho_val(nrho)
 ufac=ufac/autoev
 vnip=0.d0

 DO irho=ind_nip_rho,nrho
   vnip(irho)=ufac*((rho_val(irho)-rmin)**order)/((rmax-rmin)**order)
 ENDDO

 WRITE(Out_Unit,31) 'NIP Info:'
 WRITE(Out_Unit,32) 'order',order
 WRITE(Out_Unit,33) 'ufac ',ufac,ufac*autoev
 31 FORMAT(/1X,a,/1X,50('-'))
 32 FORMAT(1X,a,10('.'),I1)
 33 FORMAT(1X,a,10('.'),F7.4,' au',3X,F6.3,' eV')

End Subroutine nip