!  DATE (LAST UPDATE): 14/07/23        STATUS: ONGOING                 |
!  AUTHOR: JOY DUTTA, PRECURSOR'S AUTHOR: BIJIT MUKHERJEE              |
!-----------------------------------------------------------------------
!       THIS SUBROUTINE SET THE BASIS SET (I.E. QUANTUM NUMBERS,       |
!       COUPLING MATRIX ELEMENTS...) FOR THE INTERACTION TYPES:        |
!         RIGID-ROTOR + RIGID-ROTOR IN STATIC ELECTRIC FIELDS          |
!                      USING DRESSED BASIS SET                         |
!                     |n1 mn1> |n2 mn2> |L ML>                         |
!               WHERE n & mn ARE DIATOM ROTATION AND ITS               |
!        PROJECTION QUANTUM NUMBERS ON SPACE-QUANTIZATION AXIS         |
!               THIS SETS A DIAGONAL ATOMIC HAMILTONIAN.               |      
!       USER NEEDS TO PROVIDE REFERENCE ELECTRIC FIELD AT WHICH        |
!             THE DRESSED BASIS SET WILL BE CONSTRUCTED.               |
!-----------------------------------------------------------------------
!       IMPORTANT OPTIONS:                                             |
!       1. TO USE FULL PAIR BASIS FUNCTIONS - SET NPAIR = 0 (DEFAULT). |
!       2. TO USE RESTRICTED PAIR BASIS FUNCTIONS - SET NPAIR > 0.     |
!          'NPAIR' WILL DEFINE THE NUMBER OF PAIR FUNCTIONS TO BE USED.|
!          USER NEEDS TO PROVIDE QUANTUM NUMBERS FOR SUCH STATES IN    |
!          INPUT 'IPAIR' SEPARATED BY COMMAS.                          |
!       3. TO INCLUDE 2ND ORDER VAN VLECK PERTURBATIVE CORRECTION      |
!          TO A GIVEN ROTOR-ROTOR INTERACTION TERM, USER NEEDS TO      |
!          GIVE AN EXTRA LAMBDA TERM WITH THE FOLLOWING INDICES:       |  
!          LM1, LM2, LM,                                               |
!          WHERE LM1 & LM2 DEFINE THE PARTICULAR INTERACTION AND LM    |
!          MUST BE SET TO A NEGATIVE INTEGER. FOR EXAMPLE, IF THIS     |
!          CORRECTION IS NEEDED FOR DIPOLE-DIPOLE INTERACTION, THEN    |
!          USER SHOULD GIVE LM1 = 1, LM2 = 1, LM = -1.                  |
!          NOTE THAT THIS OPTION WORKS ONLY FOR RESTRICED BASIS SET.   |      
!-----------------------------------------------------------------------
      MODULE ROTORS_MW_MOD

      IMPLICIT NONE

      ! Quantities available to all routines in the basis-set suite
      INTEGER, PARAMETER        :: NQNS = 4, NLABVS = 3

      ! NQNS=5 is added by Joy Dutta

      ! Anything smaller than ZERO will be considered to be 0
      DOUBLE PRECISION, PARAMETER :: ZERO = 1.D-16

      LOGICAL :: LRSTRCT

      DOUBLE PRECISION :: FLDREF, FROT, FSTARK, DIPOLE, FOMGA, FDLTA,
     &                    ETRNS

      DOUBLE PRECISION :: OMEGA

      INTEGER          :: NMAX, NMAXSELC, NMONBASIS, NBLMN, NPAIR, NQN1,
     &                    LMIN, LMAX, LSTEP, MPAIR, POLZN, NPHMX, 
     &                    ITRNS(2), NPLVL

      INTEGER, ALLOCATABLE :: IPAIR(:), IQN(:,:), JQN(:,:), NDIM(:),
     &                        MNIDX(:), NNIDX(:,:), IMLST(:,:), JPLVL(:)

      INTEGER, ALLOCATABLE :: NMNIDXMONBASIS(:,:)
      INTEGER, ALLOCATABLE :: IDX_SELC_BIG(:)

      INTEGER, ALLOCATABLE :: NTHR(:),MNTHR(:) ! Added by Joy for
                                               ! THRSH9

      INTEGER, ALLOCATABLE :: QNBTL(:,:),QMNBTL(:,:),QNPHBTL(:,:)

      INTEGER, ALLOCATABLE :: QMTTL(:)

      INTEGER, ALLOCATABLE :: MTTHR(:)

      INTEGER :: NMONBASISSELC

      INTEGER :: NMONPHBASIS 

      DOUBLE PRECISION, ALLOCATABLE :: EVALREF(:),!, SLOPE(:),
     &                                 EVECREF(:,:,:)

      DOUBLE PRECISION, ALLOCATABLE :: EVALDRSDREF(:),
     &                                 EVECDRSDREF(:,:,:)

      INTEGER ::           INT_NC,INT_MNC,INT_NPHC,
     &                     INT_NR,INT_MNR,INT_NPHR

      INTEGER ::           LAMPH 

      DOUBLE PRECISION :: EAC_D_IN_MHZ, EAC_D_IN_INV_CM


      END MODULE ROTORS_MW_MOD
!=========================================================================
      SUBROUTINE BAS9IN(PRTP, IBOUND, IPRINT)

      USE POTENTIAL, ONLY: NCONST, NRSQ, VCONST

      USE EFVS, ONLY: MAPEFV, EFVNAM, EFVUNT, NEFV  

      USE BASIS_DATA, ONLY: JHALF

      USE ROTORS_MW_MOD, ONLY: NMAX, NMAXSELC,FLDREF, IPAIR, LMIN, LMAX,
     &                      LSTEP, NPAIR, FROT, FSTARK, DIPOLE,
     &              LRSTRCT, OMEGA, FOMGA, FDLTA, POLZN, NPHMX,
     &                          ITRNS, INT_NC,INT_MNC,INT_NPHC,
     &                          INT_NR,INT_MNR,INT_NPHR,LAMPH

      USE PHYSICAL_CONSTANTS, ONLY: HZ_IN_INV_CM, HARTREE_IN_INV_CM,
     &                              HARTREE_IN_SI, BOHR_IN_SI,
     &                              ELEMENTARY_CHARGE_IN_SI,
     &                              DEBYE_IN_SI

      IMPLICIT NONE

      CHARACTER(32), INTENT(OUT)   :: PRTP

      INTEGER,       INTENT(INOUT) :: IBOUND

      INTEGER,       INTENT(IN)    :: IPRINT

      INTEGER,       PARAMETER     :: NMAX_ALLOWED = 10,
     &                                NPAIR_ALLOWED = 5000


      DOUBLE PRECISION             :: BROT, DEL_OMG, DELTA


      DOUBLE PRECISION, PARAMETER  :: EFLD_SI_TO_AU = BOHR_IN_SI*
     &                                ELEMENTARY_CHARGE_IN_SI/
     &                                HARTREE_IN_SI,
     &                                DEBYE_IN_AU = DEBYE_IN_SI/
     &                                ELEMENTARY_CHARGE_IN_SI/
     &                                BOHR_IN_SI

      NAMELIST / BASIS9 /  NMAX, NMAXSELC, NPHMX, ITRNS, DELTA, OMEGA, 
     &                     DEL_OMG, POLZN, FLDREF,
     &                     DIPOLE, BROT, LMIN, LMAX, LSTEP, 
     &                     INT_NC,INT_MNC,INT_NPHC, 
     &                     INT_NR,INT_MNR,INT_NPHR,LAMPH 

      ! Allocate arrays which will store pair state quantum numbers
      IF (.NOT. ALLOCATED(IPAIR)) ALLOCATE(IPAIR(NPAIR_ALLOWED*4))

      ! Assign an arbitrary number to fill up the allocated space
      IPAIR = -999

      ! Set default values of variables in namelist basis9
      NMAX         = 3       ! Maximum rotational quan. no.
      NMAXSELC     = 1       ! Maximum rotational quan. no. for calc
      NPAIR        = 0       ! No. of restricted basis functions. 0 for full.
      FLDREF       = 0.D0    ! Reference electric field (kV/cm)
      DIPOLE       = 0.D0    ! Dipole momenet (Debye)
      BROT         = 0.D0    ! Rotational constant (MHz)
      LMIN         = 0       ! Minimum end-over-end angular momentum
      LMAX         = 1       ! Maximum end-over-end angular momentum
      LSTEP        = 1       ! Difference between adjacent L quan. nos
      LRSTRCT      = .FALSE. ! Flag for restricted basis set
      OMEGA        = 0.D0    ! MW Rabi frequency (2pi * MHz)
      DEL_OMG      = 0.0D0   ! Ratio of detuning and Rabi frequency
      DELTA        = 0.D0    ! MW detuning frequency (2pi * MHz)
      POLZN        = 1       ! Polarization of MW field
                             ! +/-1: Circular, 0: linear along z
      NPHMX        = 2       ! Max number of photons
      ITRNS(1)     = 0       ! ITRNS array defines the two monomer 
      ITRNS(2)     = 1       ! dressed states for the laser transition
      INT_NC  =  0           ! 
      INT_MNC =  0           !
      INT_NPHC=  0           ! Info of states of 
      INT_NR  =  1           ! expt. Rabi frequency
      INT_MNR =  1           !
      INT_NPHR=  -1          !
      LAMPH   =  1           ! Angular Momentum of Photon

      ! Read quantities to define the basis set
      READ(5,BASIS9)

      ! Notes: DEL_OMG, INT_NC, INT_MNC, INT_NPHC, INT_NR, INT_MNR
      !        ,INT_NPHR, LAMPH are added by Joy Dutta

      ! Brief Description of the interaction type
      PRTP = 'ROTOR + ROTOR IN ELECTRIC FIELD'

      ! The matrix elements of L^2 operator are diagonal values L(L+1)      
      IBOUND = 0

      ! No. of terms contributing to the asymptotic Hamiltonian
      ! Rotation + DC Stark + photon energy + molecule-field interaction
      NCONST = 1 ! Rotation+DC Stark coming from dressed monomer
                 ! energies

      ! NCONSTth term: Diagonal contribution to H_int for electric
      ! field other than EFLDREF
!     NCONST = NCONST+1

      ! +ve only if L^2 is non-diagonal
      NRSQ = 0

      ! This version uses JTOT for MTOT (projection of total J)
      JHALF = 0

      ! Number of external field variables: static electric field
      NEFV = 1

      ! Non-linear mapping between EFVs and VCONST. See EFV9 routine.
      MAPEFV = -1

      ! Name and unit of efvs
      EFVNAM(1) = 'ELECTRIC Z FIELD'
      EFVUNT(1) = ' kV/cm '

      IF (NMAX .GT. NMAX_ALLOWED) THEN
         WRITE(6,100) NMAX_ALLOWED
 100     FORMAT(/' STOP: THIS VERSION USES NMAX UPTO ',I3)
         STOP
      ENDIF

      IF (NPAIR .GT. NPAIR_ALLOWED) STOP "INCREASE ALLOCATION SPACE"//
     &    " OF IPAIR ARRAY AND RECOMPILE."

      WRITE(15,*)
      WRITE(15,*) 'NPAIR, NPAIR_ALLOWED:'
      WRITE(15,*) NPAIR, NPAIR_ALLOWED
      WRITE(15,*)


      IF (NPAIR.NE.0) LRSTRCT = .TRUE.

      ! Important prefactors
      ! FROT (cm-1):   MHz*(1.e6 Hz/MHz)*(cm-1/Hz)
      ! FSTARK (1/kV): (1.e3 V/kV)*(1.e2 m-1/cm-1)*
      !                (efld in au/Vm-1)*debye*
      !                (debye in au/debye)*(cm-1/hartree)
      !                [hartree = (efld in au)*(debye in au)]
      ! FDLTA (cm^-1)
      ! FOMGA (cm^-1)

      DELTA=DEL_OMG*ABS(OMEGA)       ! This is done by Joy Dutta
 
!     DELTA=43.475D0

      FROT   = BROT*HZ_IN_INV_CM*1.D6
      FSTARK = -1.D5*EFLD_SI_TO_AU*DIPOLE*DEBYE_IN_AU*HARTREE_IN_INV_CM
      FDLTA = DELTA*HZ_IN_INV_CM*1.D6
      FOMGA = OMEGA*HZ_IN_INV_CM*1.D6
      print*,"fdlta,fomga=",fdlta,fomga

      ! The following print statements are written by Joy Dutta
      WRITE(6,*) 
      WRITE(6,*) "==========================================="
      WRITE(6,*)
      WRITE(6,*) "Ratio of detuning and Rabi freq:", DEL_OMG
      WRITE(6,*) 
      WRITE(6,*) "==========================================="
      WRITE(6,*) 
      WRITE(6,*) "Delta obtained from DEL_OMG*OMEGA:", DELTA 
      WRITE(6,*) 
      WRITE(6,*) "==========================================="
      WRITE(6,*) "Maximum rotational quantum number for photon"
      WRITE(6,*) "dressing" 
      WRITE(6,*)  
      WRITE(6,*) 'NMAXSELC:', NMAXSELC 
      WRITE(6,*)  
      WRITE(6,*) "==========================================="
      WRITE(6,*) "Info of states of given exp. Rabi frequency" 
      WRITE(6,*) 
      WRITE(6,*) "States in the column:" 
      WRITE(6,*) 
      WRITE(6,*) "INT_NC:",  INT_NC 
      WRITE(6,*) "INT_MNC:",  INT_MNC
      WRITE(6,*) "INT_NPHC:",  INT_NPHC
      WRITE(6,*) 
      WRITE(6,*) "States in the row:" 
      WRITE(6,*) 
      WRITE(6,*) "INT_NR:",  INT_NR 
      WRITE(6,*) "INT_MNR:",  INT_MNR
      WRITE(6,*) "INT_NPHR:",  INT_NPHR
      WRITE(6,*) 
      WRITE(6,*) "==========================================="
      WRITE(6,*) "Angular momentum of photon" 
      WRITE(6,*) 
      WRITE(6,*) "LAMPH:", LAMPH 
      WRITE(6,*) 
      WRITE(6,*) "NCONST:", NCONST 

      ! Define VCONST terms in cm-1
      VCONST(1) = 1.D0

      RETURN
      END SUBROUTINE BAS9IN
!=========================================================================
      SUBROUTINE SET9(LEVIN, EIN, NSTATE, JSTATE, NQN, QNAME, NBLOCK,
     &                NLABV, IPRINT)

      USE BASIS_DATA, ONLY: NLEVEL, JLEVEL, ELEVEL

      USE ROTORS_MW_MOD, ONLY: NMAX, FLDREF, NMONBASIS, NQNS,
     &                         NQN1, NPAIR, IPAIR, IQN, NLABVS,
     &                         EVALREF, NPHMX, POLZN, ETRNS, FDLTA,
     &                         NPLVL, JPLVL, ITRNS

      USE PHYSICAL_CONSTANTS, ONLY: HZ_IN_INV_CM

      IMPLICIT NONE

      INTEGER, INTENT(OUT)      :: NQN, NBLOCK, NLABV, NSTATE,
     &                             JSTATE(*)

      CHARACTER(8), INTENT(OUT) :: QNAME(10)

      LOGICAL, INTENT(IN)       :: LEVIN, EIN

      INTEGER, INTENT(IN)       :: IPRINT

      INTEGER                   :: I, II, K, INDX

      INTEGER                   :: IA, IA1, IA2, N, MN, NPH

      INTEGER                   :: N1, MN1, N2, MN2

      DOUBLE PRECISION          :: TITIME1, TITIME2, TETIME

      CALL GCLOCK(TITIME1)

      ! LEVIN must be false otherwise program expects explicit jlevel
      ! array supplied by user from &basis, which is not the case here.
      IF (LEVIN) THEN
         WRITE(6,99)
 99      FORMAT(/"***ERROR IN SET9: LEVIN IS TRUE, MUST BE FALSE")
         STOP
      ENDIF

      ! EIN must be fasle otherwise program expects explicit elevel
      ! array supplied by user from &basis, which is not the case here.
      IF (EIN) THEN
         WRITE(6,100)
 100     FORMAT(/"***ERROR IN SET9: EIN IS TRUE, MUST BE FALSE")
         STOP
      ENDIF

      ! Print out some general information      
      WRITE(6,101)
 101  FORMAT(/"  BASIS SET IS SYMMETRIZED PRODUCT OF INTERNAL DRESSED",
     &        " STATES OF DIATOMS AND PHOTON STATES, WHICH IN GENERAL",
     &        " READS |n1 mn1> |n2 mn2> |N> |L ML>")

      WRITE(6,102) NMAX
 102  FORMAT(/'  BASIS SET USES MONOMER DIATOM NMAX =',I3)

      WRITE(6,103) FLDREF
 103  FORMAT(/'  DRESSED BASIS SET WILL BE CONSTRUCTED AT', 
     &        ' REFERENCE ELECTRIC FIELD =', F8.4, ' kV/cm')


      ! The following subroutine is wriiten by Joy Dutta
      CALL EAC_DIP_DETER


      ! Total monomer basis functions      
      NMONBASIS = (NMAX+1)*(NMAX+1)
      
      CALL GCLOCK(TITIME2)

      ! Sets up mapping between a single quantum number to the set of
      ! monomer quantum numbers for a given monomer state
      CALL MONOMER_MAPPING

!     CALL MONOMER_ENERGIES



      CALL GCLOCK(TETIME)
      WRITE(6,*)
      WRITE(6,*)' TIME IN MONOMER_MAPPING:',TETIME-TITIME2

      IF (IPRINT .GE. 6) WRITE(6,104) NMONBASIS
 104  FORMAT(/'  NUMBER OF MONOMER BASIS FUNCTIONS =',I5)
      
      WRITE(6,105)
 105  FORMAT(/"  ASSIGNING A NON-NEGATIVE INTEGER QUANTUM NUMBER",
     &        " 'a' TO REPRESENT EACH MONOMER STATE")

      WRITE(6,106)
 106  FORMAT(/"  MAPPING OF MONOMER QUANTUM NUMBERS:",
     &       /"  ----------------------------- ",
     &       /"       a       n       mn       ",
     &       /"  ----------------------------- ")

      DO IA = 1, NMONBASIS
         N  = IQN(IA,1)
         MN = IQN(IA,2)
         WRITE(6,'(5I8)') IA-1,N,MN
      ENDDO

      ! Number of symmetry blocks for each jtot: bosonic & fermionic
      NBLOCK = 2

      ! Number of indices needed to describe each term
      NLABV = NLABVS 
      
      ! One greater than the number of quantum labels      
      NQN  = NQNS
      NQN1 = NQN-1

      ! Sanity check on IPAIR array
      DO I = 1, NPAIR*2
         IF (IPAIR(I) .LT. 0) 
     &   STOP "'INCOMPATIBLE INPUTS FOR NPAIR AND IPAIR IN &BASIS9"//
     &        " OR ILLEGAL VALUE IN IPAIR'"
      ENDDO

      CALL GCLOCK(TITIME2)

      ! Sets up dressed monomer and pair eigenvalues and eigenvectors 
      CALL BASIS_BUILDER


      CALL GCLOCK(TETIME)
      WRITE(6,*)
      WRITE(6,*)' TIME IN BASIS_BUILDER:',TETIME-TITIME2

      QNAME(1) = '    a1  '
      QNAME(2) = '    a2  '
      QNAME(3) = '     N  '

      WRITE(15,*)'====================================================' 
      WRITE(15,*)
      WRITE(15,*)'List of pair levels: Sl no. Mon1 idx, N, MN, Mon2 idx'
      WRITE(15,*)'N, MN, NPH'
      WRITE(15,*)
      WRITE(15,*)
      WRITE(15,*) '==================================================='
      WRITE(15,*)
      WRITE(15,*) 'Checking NPAIR Value:', NPAIR
      WRITE(15,*)
      WRITE(15,*) 'Size of IPAIR array:', SIZE(IPAIR)
      WRITE(15,*) '==================================================='


      K = 0
      NLEVEL = 0
      ! NPAIR and IPAIR get updated from BASIS_BUILDER
      DO II = 1, NPAIR
         IA1 = IPAIR(K+1)
         IA2 = IPAIR(K+2)
         N1  = IQN(IA1+1,1)
         MN1 = IQN(IA1+1,2)
         N2  = IQN(IA2+1,1)
         MN2 = IQN(IA2+1,2)
         DO NPH = 0, -NPHMX, -1
!        DO NPH =  -NPHMX, NPHMX
            IF (N1+N2+NPH.NE.2) CYCLE    ! For N=1<===>N=2 Transition: Joy Dutta 
!           IF ((N1+N2+NPH) .GT. 1) CYCLE
!           IF (IA1.NE.ITRNS(1).AND.IA2.NE.ITRNS(1).AND.
!    &          MN1.NE.POLZN.AND.MN2.NE.POLZN) CYCLE
!           IF (N1.NE.ITRNS(1).AND.N2.NE.ITRNS(2).AND.
!    &          MN1.NE.(POLZN-1).AND.MN2.NE.POLZN) CYCLE
!           IF (((N1+N2) /= 1) ) CYCLE 
!           IF (((N1+N2) == 2).AND.((MN1/=1).AND.(MN2/=1))) CYCLE 
!           IF (((N1+N2) == 2).AND.((MN1+MN2)/=2)) CYCLE 
            JLEVEL(NQN1*NLEVEL+1) = IA1 
            JLEVEL(NQN1*NLEVEL+2) = IA2 
            JLEVEL(NQN1*NLEVEL+3) = NPH
            NLEVEL = NLEVEL+1
         WRITE(15,'(8I3)') NLEVEL, IA1+1, N1, MN1, IA2+1, N2, MN2, NPH 
         ENDDO
         K = K+2
      ENDDO

      NSTATE = NLEVEL

      WRITE(15,*)
      WRITE(15,*)'====================================================' 
      WRITE(15,*)
      WRITE(15,*) 'Total No. of Pair level:', NSTATE
      WRITE(15,*)
      WRITE(15,*)'====================================================' 

      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) "==================================================="
      WRITE(15,*)  "Pair State and Pair Level relations" 
      WRITE(15,*) "==================================================="

      II = 0
      DO I = 1, NSTATE
         II = II+1
         INDX = NQN1*(I-1)
         DO K = 0, NQN1-1
            JSTATE(I+K*NSTATE) = JLEVEL(INDX+K+1)
            WRITE(15,*) I+K*NSTATE, INDX+K+1 
         ENDDO
         JSTATE(I+NSTATE*NQN1) = II
      ENDDO
 

      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) "==================================================="
      WRITE(15,*) 
      WRITE(15,*) 




      WRITE(15,*) "==================================================="
      WRITE(15,*) 
      WRITE(15,*) "**** Printing the pair state details: ****" 
      WRITE(15,*) 
      WRITE(15,*) "Pair state serial numbers, Corresponding monomer 1" 
      WRITE(15,*) "monomer 2 serial numbers, values of monomer quantum"
      WRITE(15,*) "numbers"
      WRITE(15,*) 
      WRITE(15,*) "HZ_IN_INV_CM", HZ_IN_INV_CM 
      WRITE(15,*) 
      WRITE(15,*) "==================================================="
      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) 

      DO II = 1, NSTATE
         IA1 = JLEVEL(NQN1*(II-1)+1)
         IA2 = JLEVEL(NQN1*(II-1)+2)
         NPH = JLEVEL(NQN1*(II-1)+3)
         N1  = IQN(IA1+1,1)
         MN1 = IQN(IA1+1,2)
         N2  = IQN(IA2+1,1)
         MN2 = IQN(IA2+1,2)
         ELEVEL(II) = EVALREF(IA1+1)+EVALREF(IA2+1)
     &              + (ETRNS+FDLTA)*NPH
      
      WRITE(15,'(6I3, f18.6)') II, N1, MN1, N2, MN2, NPH, 
     &            (ELEVEL(II)/(HZ_IN_INV_CM*1.D6))  
      ENDDO
      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) "==================================================="
!     STOP 123

!     NPLVL  = NLEVEL
!     IF(.NOT. ALLOCATED(JPLVL)) ALLOCATE(JPLVL(NQN1*NPLVL))
!     JPLVL = -999

!     JPLVL = JLEVEL(1:NQN1*NPLVL)
!     JPLVL = JLEVEL(NQN1*NPLVL)

!     ! Purpose of the following array finished hence deallocate
!     IF (ALLOCATED(IPAIR)) DEALLOCATE(IPAIR)

      CALL GCLOCK(TETIME)
      WRITE(6,*)
      WRITE(6,*)' TIME IN SET9:',TETIME-TITIME1

!     WRITE(20,*) JPLVL
!     WRITE(20,*) 
!     WRITE(20,*)
!     DO II=1,NSTATE
!        JPLVL(NQN1*(II-1)+1) = JLEVEL(NQN1*(II-1)+1)
!        JPLVL(NQN1*(II-1)+2) = JLEVEL(NQN1*(II-1)+2)
!        JPLVL(NQN1*(II-1)+3) = JLEVEL(NQN1*(II-1)+3)
!     ENDDO

!     WRITE(20,*) JPLVL
!     WRITE(20,*) 
!     WRITE(20,*) 

!     DO II=1,NSTATE
!       IA1=JPLVL(NQN1*(II-1)+1)
!       IA2=JPLVL(NQN1*(II-1)+2)
!       NPH=JPLVL(NQN1*(II-1)+3)
!        N1  = IQN(IA1+1,1)
!        MN1 = IQN(IA1+1,2)
!        N2  = IQN(IA2+1,1)
!        MN2 = IQN(IA2+1,2)
!     WRITE(20,'(6I3)') II, N1, MN1, N2, MN2, NPH
!     ENDDO
      RETURN
      END SUBROUTINE SET9
!=========================================================================
      SUBROUTINE BASE9(LCOUNT, N, JTOT, IBLOCK, JSTATE, NSTATE, NQN,
     &                 JSINDX, L, IPRINT)

      ! Contructs basis set using pair basis functions (dressed) and L
      ! quantum number for the current values of JTOT and IBLOCK

      USE ROTORS_MW_MOD, ONLY: LMIN, LMAX, LSTEP, IQN, POLZN

      IMPLICIT NONE

      INTEGER, INTENT(INOUT) :: N

      INTEGER, INTENT(OUT)   :: JSINDX(N), L(N)

      LOGICAL, INTENT(IN)    :: LCOUNT

      INTEGER, INTENT(IN)    :: JTOT, IBLOCK, NSTATE, NQN,
     &                          JSTATE(NSTATE,NQN), IPRINT

      INTEGER :: IFUNC, ISTATE, IA1, IA2, LL, MTOT, NPH

      INTEGER :: MN1, MN2

      INTEGER :: NN1, NN2         ! Done by J.D 

      DOUBLE PRECISION :: TITIME, TETIME

      CALL GCLOCK(TITIME)

      MTOT = JTOT


         WRITE(15,*) 
         WRITE(15,*) 
         WRITE(15,*) 
         WRITE(15,*) "==============================================="
         WRITE(15,*) 
         WRITE(15,*) "_______________________________________________" 
         WRITE(15,*) 
         WRITE(15,*) "IBLOCK:", IBLOCK, "LCOUNT:", LCOUNT 
         WRITE(15,*) 
         WRITE(15,*) 
         WRITE(15,*) "_______________________________________________" 
         WRITE(15,*) 
         WRITE(15,*) 

      IFUNC = 0
      DO ISTATE = 1, NSTATE
         IA1 = JSTATE(ISTATE,1)
         IA2 = JSTATE(ISTATE,2)
         NPH = JSTATE(ISTATE,3)
         ! Access explicit quantum numbers from IQN array
         MN1 = IQN(IA1+1,2)
         MN2 = IQN(IA2+1,2)

         NN1 = IQN(IA1+1,1) ! Done by J.D
         NN2 = IQN(IA2+1,1) !
      

         WRITE(15,*) 
         WRITE(15,'(8I3)')  ISTATE, IA1+1, NN1, MN1, IA2+1, NN2, MN2,
     &                      NPH
         WRITE(15,*) 

         DO LL = LMIN, LMAX, LSTEP
            IF (ABS(MTOT-MN1-MN2-NPH*POLZN) .GT. LL) CYCLE
            IF (IBLOCK .EQ. 1) THEN
               ! Fermionic symmetry: identical states must have odd L
               IF ((IA1.EQ.IA2) .AND. (MOD(LL,2).EQ.0)) CYCLE
               IFUNC = IFUNC+1
               IF (.NOT. LCOUNT) THEN
                  JSINDX(IFUNC) = ISTATE
                  L(IFUNC) = LL
               ENDIF
            ELSEIF (IBLOCK .EQ. 2) THEN
               ! Bosonic symmetry: identical states must have even L
               IF ((IA1.EQ.IA2) .AND. (MOD(LL,2).NE.0)) CYCLE
               IFUNC = IFUNC+1
               IF (.NOT. LCOUNT) THEN
                  JSINDX(IFUNC) = ISTATE
                  L(IFUNC) = LL
!              WRITE(16,'(2I3)') L(IFUNC), MTOT-MN1-MN2-NPH*POLZN
               ENDIF
            ELSE
               STOP "'WRONG IBLOCK PASSED INTO BASE9'"
            ENDIF
         ENDDO
!        WRITE(16,*) 
      ENDDO

      N = IFUNC

!     STOP 123

      CALL GCLOCK(TETIME)
      WRITE(6,*)
      WRITE(6,*)' TIME IN BASE9:',TETIME-TITIME

      RETURN
      END SUBROUTINE BASE9
!=========================================================================
      SUBROUTINE CPL9(N, IBLOCK, NHAM, LAM, MXLAM, NSTATE, JSTATE,
     &                JSINDX, L, JTOT, VL, IV, CENT, DGVL, IBOUND,
     &                IEXCH, IPRINT)

      ! Sets up the matrix elements for the potential and other
      ! operators between the basis functions set up in base9 for the
      ! current values of JTOT & IBLOCK

!     USE ROTORS_MW_MOD, ONLY: SLOPE, NQNS, NLABVS, IQN, LMIN, LMAX,
!    &                         LSTEP, MPAIR, IMLST, EVALREF, LRSTRCT,
!    &                         ZERO, POLZN, ETRNS, FDLTA

      USE ROTORS_MW_MOD, ONLY: NQNS, NLABVS, IQN, LMIN, LMAX,
     &                         LSTEP, MPAIR, IMLST, EVALREF, LRSTRCT,
     &                         ZERO, POLZN, ETRNS, FDLTA

      USE POTENTIAL, ONLY: NDGVL, NVLBLK

      USE PHYSICAL_CONSTANTS, ONLY: HZ_IN_INV_CM

      IMPLICIT NONE

      DOUBLE PRECISION, INTENT(OUT) :: VL(NVLBLK,N*(N+1)/2), CENT(N),
     &                                 DGVL(N,NDGVL)

      INTEGER, INTENT(OUT)          :: IV(NVLBLK,N*(N+1)/2)

      INTEGER, INTENT(IN)           :: N, IBLOCK, NHAM, MXLAM,
     &                                 NSTATE, JSTATE(NSTATE,NQNS),
     &                                 JSINDX(N), L(N), JTOT, IBOUND,
     &                                 IEXCH, IPRINT, LAM(*)

      INTEGER :: IRC, ICOL, IROW, IHAM, LM1, LM2, LM, LMSUM, MTOT, IEX

      INTEGER :: II, LC, JA1C, JA2C, N1C, MN1C, NPHC

      INTEGER :: LR, JA1R, JA2R, N1R, MN1R, NPHR

      INTEGER :: N2C, MN2C, N2R, MN2R

      INTEGER :: IETAC, IETAR, IETAK, LMX

      INTEGER :: NPHDIFF

      INTEGER :: N1K, MN1K, N2K, MN2K, LK, IINT, I1K, I2K

      DOUBLE PRECISION :: TERM

      DOUBLE PRECISION :: TERMLOC, TERMLOC1, TERMLOC2
      
      DOUBLE PRECISION :: TERMLOC3,TERMLOC4

      DOUBLE PRECISION :: ENRC, ENRK, ENRR, EDIFC, EDIFR

      DOUBLE PRECISION :: P1, P2C, P2K, P3CR, P3CK, P3KR
 
      DOUBLE PRECISION :: ANORMC, ANORMR, ANORMK

      DOUBLE PRECISION :: HVVCK, HVVKR 

      DOUBLE PRECISION :: T1TIME, T2TIME

      LOGICAL :: LCOLDIF, LROWDIF, LINTDIF

      LOGICAL :: DDMX=.FALSE.

      DOUBLE PRECISION, PARAMETER :: SQRTH = SQRT(0.5D0)

      ! Added by Joy Dutta
      INTEGER :: ITEST
      !


      ! Functions
      INTEGER(KIND=8)  :: FACTORIAL

      DOUBLE PRECISION :: THREEJ

      WRITE(6,*)
      IF (IBLOCK .EQ. 1) 
     &  WRITE(6,*)" THIS SYMMETRY BLOCK CONSIDERS IDENTICAL FERMIONS"

      IF (IBLOCK .EQ. 2) 
     &  WRITE(6,*)" THIS SYMMETRY BLOCK CONSIDERS IDENTICAL BOSONS"

      MTOT = JTOT

      IEX = 1
      IF (IBLOCK .EQ. 1) IEX = -1 

      ! Calculate the coupling matrices and store them in the VL array.
      VL = 0.D0
      TERMLOC=0.0D0
      TERMLOC1=0.0D0
      TERMLOC2=0.0D0
      TERMLOC3=0.0D0
      TERMLOC4=0.0D0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! Added by Joy Dutta
      WRITE(16,*) 
      WRITE(16,*) 
      WRITE(16,*)"========================================="
      WRITE(16,*) 
      WRITE(16,*) "Total Number of Basis function, N=", N
      WRITE(16,*) 
      WRITE(16,*) "MTOT:", MTOT 
      WRITE(16,*) 
      WRITE(16,*) "POLZN:", POLZN 
      WRITE(16,*) 
      WRITE(16,*) "NHAM:", NHAM 
      WRITE(16,*) 
      WRITE(16,*) "NDGVL:", NDGVL, "NVLBLK:", NVLBLK 
      WRITE(16,*) 
      WRITE(16,*) "Cond. for evaluating dipole-dipole matrix:", DDMX 
      WRITE(16,*)"========================================="
      WRITE(16,*) 
      WRITE(16,*) 

      WRITE(16,*)"========================================="
      WRITE(16,*)
      WRITE(16,*)"Details of Basis set function:"
      WRITE(16,*)
      WRITE(16,*) "SL No., MON1-INDX, N1, MN1"
      WRITE(16,*) "MON2-INDX, N2, MN2, NPHC, L, ML"
      WRITE(16,*)
      WRITE(16,*)
      WRITE(16,*)"========================================="

      DO ITEST=1,N
      JA1C = JSTATE(JSINDX(ITEST),1)+1
              N1C  = IQN(JA1C,1)
              MN1C = IQN(JA1C,2)
      JA2C = JSTATE(JSINDX(ITEST),2)+1
              N2C  = IQN(JA2C,1)
              MN2C = IQN(JA2C,2)
      NPHC = JSTATE(JSINDX(ITEST),3)
      LC   = L(ITEST)

      WRITE(16,'(11I3)') ITEST, JA1C, N1C, MN1C,
     &   JA2C, N2C, MN2C, NPHC, LC, MTOT-MN1C-MN2C-POLZN*NPHC,
     &   MN1C+MN2C+POLZN*NPHC 
      ENDDO 


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!





      ! IHAM: Outermost loop over different independent blocks of the
      ! potential
      ! First loop up to MXLAM terms each of which are defined by a set
      ! of LM1,LM2,LM
      DO IHAM = 1, MXLAM

         CALL GCLOCK(T1TIME)

         II = (IHAM-1)*NLABVS

         ! Get the lambda values to identify the potential term
         LM1 = LAM(II+1)
         LM2 = LAM(II+2)
         LM  = LAM(II+3)

         LMX = MAX(LM1,LM2)

         IF (LM1.LT.0 .OR. LM2.LT.0)
     &   STOP 'STOP: CPL9: ILLEGAL VALUES IN LAMBDA ARRAY OF &POTL'

         ! Isotropic interaction term
         IF (LM1.EQ.0 .AND. LM2.EQ.0 .AND. LM.EQ.0) THEN
            DO ICOL = 1, N
               IRC = ICOL+ICOL*(ICOL-1)/2
               VL(IHAM,IRC) = 1.D0
            ENDDO
         ENDIF

         ! Any multipole-multipole potential term
         IF (LM1+LM2.NE.0 .AND. LM.GE.ABS(LM1-LM2) .AND.
     &       LM.LE.LM1+LM2) THEN

           LMSUM = LM1+LM2

           P1 = SQRT(DBLE(FACTORIAL(2*LM1+2*LM2+1))/
     &               DBLE(FACTORIAL(2*LM1)*FACTORIAL(2*LM2)))

           ! Loop over explicit basis set
           ! Column => Bra basis functions


!          IF(DDMX) THEN   ! Condition for DD-matrix evaluation by J.D

           WRITE(16,*)"========================================="
           WRITE(16,*)
           WRITE(16,*) "Cond. for DD-matrix calcuation:", DDMX
           WRITE(16,*)
           WRITE(16,*)"========================================="

           DO ICOL = 1, N
              
              LC   = L(ICOL)
              JA1C = JSTATE(JSINDX(ICOL),1)+1
              JA2C = JSTATE(JSINDX(ICOL),2)+1
              NPHC = JSTATE(JSINDX(ICOL),3)
              LCOLDIF = .TRUE.
              ANORMC  = SQRTH
              IF (JA1C .EQ. JA2C) THEN
                LCOLDIF = .FALSE.
                ANORMC  = 1.D0
              ENDIF

              IETAC = IEX
              IF (MOD(LC,2).NE.0) IETAC = -IEX

              ! Get actual quantum numbers from IQN array
              N1C  = IQN(JA1C,1)
              MN1C = IQN(JA1C,2)

              N2C  = IQN(JA2C,1)
              MN2C = IQN(JA2C,2)

              P2C = P1*SQRT(DBLE(2*LC+1))

           ! Row => Ket basis functions
           DO IROW = 1, ICOL

              LR   = L(IROW)
              ! 3-j restrictions
              IF (LC+LR .EQ. 0) CYCLE
              IF (ABS(LC-LR) .GT. LMSUM) CYCLE

              ! Parity restriction
              IF (LM1.EQ.LM2 .AND. MOD(LC+LR,2).NE.0) CYCLE

              JA1R = JSTATE(JSINDX(IROW),1)+1
              JA2R = JSTATE(JSINDX(IROW),2)+1
              NPHR = JSTATE(JSINDX(IROW),3)

              ! Dipole-dipole matrix elements are diagonal in photon
              ! states
              IF (NPHC .NE. NPHR) CYCLE

              ! Set normalization constant according to whether the pair
              ! function is a repeated one or non-repeated one.
              LROWDIF = .TRUE.
              ANORMR  = SQRTH
              IF (JA1R .EQ. JA2R) THEN
                LROWDIF = .FALSE.
                ANORMR  = 1.D0
              ENDIF

              IETAR = IEX
              IF (MOD(LR,2).NE.0) IETAR = -IEX

              N1R  = IQN(JA1R,1)
              MN1R = IQN(JA1R,2)

              N2R  = IQN(JA2R,1)
              MN2R = IQN(JA2R,2)

              P3CR = P2C*SQRT(DBLE(2*LR+1))*THREEJ(LC,LMSUM,LR)
              IF (ABS(P3CR) .LT. ZERO) CYCLE

              IRC = IROW+ICOL*(ICOL-1)/2

              ! Calculates multipole-multipole matrix element taking
              ! into effect proper exchange symmetries
              CALL VMMDR(N1C,MN1C,N2C,MN2C,N1R,MN1R,N2R,MN2R,
     &             LCOLDIF,ANORMC,IETAC,LROWDIF,ANORMR,IETAR,
     &             LC,LR,LM1,LM2,LMSUM,MTOT-NPHC*POLZN,VL(IHAM,IRC)) 

              VL(IHAM,IRC) = P3CR*VL(IHAM,IRC)

              IF(ICOL /= IROW) THEN
              IF((LC==0 .AND. LR==2) .AND. (NPHC==NPHR)) THEN
!             IF((N1C==0).AND.(MN1C==0).AND.(N2C==1).AND.(MN2C==1)) THEN
!             IF((N1R==1).AND.(MN1R==1).AND.(N2R==0).AND.(MN2R==0)) THEN
              WRITE(18,'(8I3,F12.6)') N1C, MN1C, N2C, MN2C, N1R, MN1R, 
     &        N2R, MN2R, VL(IHAM,IRC)  
!             ENDIF 
!             ENDIF 
              ENDIF 
              ENDIF 

           ENDDO ! IROW loop
           ENDDO ! ICOL loop

         ! -------------------------------------------------------------------------
         ! Potential 1st order term evaluation ends here
         ! -------------------------------------------------------------------------
         ELSEIF (LM.LT.0 .AND. LRSTRCT) THEN
         ! Evaluation of 2nd order Van-Vleck terms
         ! LM = -1: 2nd order Van-Vleck in multipole-multipole

           LMSUM = LM1+LM2

           P1 = SQRT(DBLE(FACTORIAL(2*LM1+2*LM2+1))/
     &               DBLE(FACTORIAL(2*LM1)*FACTORIAL(2*LM2)))

           ! Van Vleck: loop over Class-I functions: explicit basis set
           ! Column => Bra basis functions
           DO ICOL = 1, N
              
              LC   = L(ICOL)
              JA1C = JSTATE(JSINDX(ICOL),1)+1
              JA2C = JSTATE(JSINDX(ICOL),2)+1
              NPHC = JSTATE(JSINDX(ICOL),3)
              LCOLDIF = .TRUE.
              ANORMC  = SQRTH
              IF (JA1C .EQ. JA2C) THEN
                LCOLDIF = .FALSE.
                ANORMC  = 1.D0
              ENDIF

              IETAC = IEX
              IF (MOD(LC,2).NE.0) IETAC = -IEX

              ! Get actual quantum numbers from IQN array
              N1C  = IQN(JA1C,1)
              MN1C = IQN(JA1C,2)

              N2C  = IQN(JA2C,1)
              MN2C = IQN(JA2C,2)

              P2C = P1*SQRT(DBLE(2*LC+1))

              ! Get the Stark and rotation energy
              ENRC = EVALREF(JA1C)+EVALREF(JA2C)

           ! Row => Ket basis functions
           DO IROW = 1, ICOL

              LR   = L(IROW)

              ! Parity restriction
              IF (LM1.EQ.LM2 .AND. MOD(LC+LR,2).NE.0) CYCLE

              ! 3-j restriction
              IF (ABS(LC-LR).GT.2*LMSUM) CYCLE
              
              JA1R = JSTATE(JSINDX(IROW),1)+1
              JA2R = JSTATE(JSINDX(IROW),2)+1
              NPHR = JSTATE(JSINDX(IROW),3)

              ! VV matrix elements are diagonal in photon states
              IF (NPHC .NE. NPHR) CYCLE

              LROWDIF = .TRUE.
              ANORMR  = SQRTH
              IF (JA1R .EQ. JA2R) THEN
                LROWDIF = .FALSE.
                ANORMR  = 1.D0
              ENDIF
              
              IETAR = IEX
              IF (MOD(LR,2).NE.0) IETAR = -IEX
              
              N1R  = IQN(JA1R,1)
              MN1R = IQN(JA1R,2)
              
              N2R  = IQN(JA2R,1)
              MN2R = IQN(JA2R,2)
              
              ! Get the Stark and rotation energy
              ENRR = EVALREF(JA1R)+EVALREF(JA2R)
              
              IRC = IROW+ICOL*(ICOL-1)/2
              
           ! Now start the sum over intermediate left over basis
           ! functions
           ! Van Vleck: loop over Class-II functions: excluded spin-free
           ! pair levels
           DO IINT = 1, MPAIR

              I1K = IMLST(IINT,1)
              I2K = IMLST(IINT,2)
              
              N1K  = IQN(I1K,1)
              MN1K = IQN(I1K,2)

              N2K  = IQN(I2K,1)
              MN2K = IQN(I2K,2)

              ! Exploit some of 3j-symbol restrictions to skip
              ! intermediate pair functions
              IF (MN1K.EQ.MN2K .AND. (ABS(MN1K-MN1C).GT.LMX.OR.
     &           ABS(MN1K-MN2C).GT.LMX)) CYCLE

              IF (MN1C.EQ.MN2C .AND. (ABS(MN1C-MN1K).GT.LMX.OR.
     &           ABS(MN1C-MN2K).GT.LMX)) CYCLE

              IF (ABS(MN1K-MN1C).GT.LMX .AND. 
     &            ABS(MN1K-MN2C).GT.LMX) CYCLE
              IF (ABS(MN2K-MN1C).GT.LMX .AND. 
     &            ABS(MN2K-MN2C).GT.LMX) CYCLE

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

              ! Get the Stark and rotation energy
              ENRK = EVALREF(I1K)+EVALREF(I2K)

              EDIFC = ENRC-ENRK
              EDIFR = ENRR-ENRK
              IF(ABS(EDIFC).LT.ZERO .OR. ABS(EDIFR).LT.ZERO) THEN
                WRITE(*,*) '  STOP: CPL9: DEGENERATE STATES IN CLASS-I',
     &                     ' AND II BASIS SETS IN VAN VLECK'
                STOP
              ENDIF

              LINTDIF = .TRUE.
              ANORMK  = SQRTH
              IF (I1K .EQ. I2K) THEN
                LINTDIF = .FALSE.
                ANORMK  = 1.D0
              ENDIF

           ! Now loop over L functions for the intermediate pair
           ! functions
           DO LK = ABS(LC-LMSUM), LC+LMSUM, LSTEP
              IF (LK.LT.LMIN .OR. LK.GT.LMAX) CYCLE
              IF (ABS(MTOT-MN1K-MN2K-NPHC*POLZN).GT.LK) CYCLE

              ! Parity restriction
              IF (LM1.EQ.LM2 .AND. MOD(LC+LK,2).NE.0) CYCLE

              ! 3-j restrictions for multipole-multipole Van-Vleck term
              IF (LM.EQ.-1 .AND. LK+LR.EQ.0) CYCLE
              IF (LM.EQ.-1 .AND. ABS(LK-LR).GT.LMSUM) CYCLE

              IETAK = IEX
              IF (MOD(LK,2).NE.0) IETAK = -IEX

              P3CK = P2C*SQRT(DBLE(2*LK+1))*THREEJ(LC,LMSUM,LK)
              IF (ABS(P3CK) .LT. ZERO) CYCLE

              P2K = P1*SQRT(DBLE(2*LK+1))

              P3KR = P2K*SQRT(DBLE(2*LR+1))*THREEJ(LK,LMSUM,LR)
              IF (ABS(P3KR) .LT. ZERO) CYCLE

              ! Calculate multipole-multipole matrix element between bra
              ! and intermediate basis functions
              CALL VMMDR(N1C,MN1C,N2C,MN2C,N1K,MN1K,N2K,MN2K,
     &                   LCOLDIF,ANORMC,IETAC,LINTDIF,ANORMK,IETAK,
     &                   LC,LK,LM1,LM2,LMSUM,MTOT-NPHC*POLZN,HVVCK)

              IF (ABS(HVVCK) .LT. ZERO) CYCLE

              HVVCK = P3CK*HVVCK
              HVVKR = 0.D0

              CALL VMMDR(N1K,MN1K,N2K,MN2K,N1R,MN1R,N2R,MN2R,
     &                   LINTDIF,ANORMK,IETAK,LROWDIF,ANORMR,IETAR,
     &                   LK,LR,LM1,LM2,LMSUM,MTOT-NPHC*POLZN,HVVKR)

              HVVKR = P3KR*HVVKR

              TERM = 0.5D0*HVVCK*HVVKR*((1.D0/EDIFC)+(1.D0/EDIFR))
              
              VL(IHAM,IRC) = VL(IHAM,IRC)+TERM

           ENDDO ! LK loop
           ENDDO ! IINT (Class-II pair functions) loop ends here

           ENDDO ! Ket N states loop ends here
           ENDDO ! Bra N states: Class-I function explicit basis set loop 

         ! -------------------------------------------------------------------------
         ! Evaluation of Van-Vleck terms ends here
         ! -------------------------------------------------------------------------

!          ENDIF ! Condition of DD-matrix evaluation by J.D

         ENDIF ! multipole-multipole potential term


         CALL GCLOCK(T2TIME)
         WRITE(6,*)
         WRITE(6,*)' CPL9 TIME:',T2TIME-T1TIME,' IN POTENTIAL TERM:',
     &             IHAM

      ENDDO ! IHAM loop

      !--------------------------------------------------------------
      ! Now build VL matrix elements for asymptotic Hamiltonian terms
      !--------------------------------------------------------------

      CALL GCLOCK(T1TIME)

       WRITE(17,*) "================================================" 
       WRITE(17,*)   
       WRITE(17,*)   
       WRITE(17,*) "Info of Basis sl no., quantum numbers and energies" 
       WRITE(17,*)   
       WRITE(17,*) "Energies in MHz" 
       WRITE(17,*)   
       WRITE(17,*) "================================================"
       WRITE(17,*)   
       WRITE(17,*)   
      ! N basis: Outermost loop over explicit basis set
      ! Bra basis functions 
      DO ICOL = 1, N
         
         LC   = L(ICOL)
         JA1C = JSTATE(JSINDX(ICOL),1)+1
         JA2C = JSTATE(JSINDX(ICOL),2)+1
         NPHC = JSTATE(JSINDX(ICOL),3)

         LCOLDIF = .TRUE.
         ANORMC  = SQRTH
         IF (JA1C .EQ. JA2C) THEN
           LCOLDIF = .FALSE.
           ANORMC  = 1.D0
         ENDIF

         IETAC = IEX
         IF (MOD(LC,2).NE.0) IETAC = -IEX

         ! Get actual quantum numbers from IQN array
         N1C  = IQN(JA1C,1)
         MN1C = IQN(JA1C,2)

         N2C  = IQN(JA2C,1)
         MN2C = IQN(JA2C,2)

      DO IROW = 1, ICOL
         
         LR   = L(IROW)

         ! Asymptotic matrix elements are diagonal in L
         IF (LR .NE. LC) CYCLE

         JA1R = JSTATE(JSINDX(IROW),1)+1
         JA2R = JSTATE(JSINDX(IROW),2)+1
         NPHR = JSTATE(JSINDX(IROW),3)

         LROWDIF = .TRUE.
         ANORMR  = SQRTH
         IF (JA1R .EQ. JA2R) THEN
           LROWDIF = .FALSE.
           ANORMR  = 1.D0
         ENDIF

         IETAR = IEX
         IF (MOD(LR,2).NE.0) IETAR = -IEX

         N1R  = IQN(JA1R,1)
         MN1R = IQN(JA1R,2)

         N2R  = IQN(JA2R,1)
         MN2R = IQN(JA2R,2)

         NPHDIFF=NPHR-NPHC    ! Defined by J.D

         ! Total projection of pair state should be conserved
         IF (MN1C+MN2C+POLZN*NPHC.NE.MN1R+MN2R+POLZN*NPHR) CYCLE

         IRC = IROW+ICOL*(ICOL-1)/2

         ! Diagonal terms
         IF (ICOL .EQ. IROW) THEN

           ! Get the Stark and rotation energy
           VL(MXLAM+1,IRC) = EVALREF(JA1C)+EVALREF(JA2C)
           ! Photon energy
           TERMLOC=(ETRNS+FDLTA)*NPHC
           VL(MXLAM+1,IRC) = VL(MXLAM+1,IRC)+TERMLOC

         WRITE(17,'(7I3,f18.6)') ICOL, N1C, MN1C, N2C, MN2C, NPHC, LC 
     &              ,(VL(MXLAM+1,IRC)/(HZ_IN_INV_CM*1.0D6))  

         ENDIF

         IF (ABS(NPHC-NPHR).NE.1) CYCLE
                 
         ! Hm: Monomer Hamiltonian
         ! <1C|Hm|1R>*krondel(2C,2R)
!        IF (JA2C.EQ.JA2R .AND. ABS(MN1C-MN1R).LE.POLZN) THEN
         ! This modified by J.D
         IF (JA2C.EQ.JA2R .AND. ((MN1C-MN1R) == POLZN*NPHDIFF)) THEN 
            CALL VMFDR(N1C,MN1C,NPHC,N1R,MN1R,NPHR,TERM)
            TERMLOC1=TERM*ANORMC*ANORMR*2.0D0
            VL(MXLAM+1,IRC) = VL(MXLAM+1,IRC)+TERMLOC1
         ENDIF

         IF (LCOLDIF) THEN
            ! <2C|Hm|1R>*krondel(1C,2R)
!           IF (JA1C.EQ.JA2R .AND. ABS(MN2C-MN1R).LE.POLZN) THEN
         ! This modified by J.D
         IF (JA1C.EQ.JA2R .AND. ((MN2C-MN1R) == POLZN*NPHDIFF)) THEN 
               CALL VMFDR(N2C,MN2C,NPHC,N1R,MN1R,NPHR,TERM)
               TERMLOC2=TERM*IETAC*ANORMC*ANORMR*2.0D0
               VL(MXLAM+1,IRC) = VL(MXLAM+1,IRC)
     &                         + TERMLOC2 
         ENDIF
         ENDIF ! LCOLDIF IF block

         IF (LROWDIF) THEN
            ! <1C|Hm|2R>*krondel(2C,1R)
!           IF (JA2C.EQ.JA1R .AND. ABS(MN1C-MN2R).LE.POLZN) THEN
         ! This modified by J.D
         IF (JA2C.EQ.JA1R .AND. ((MN1C-MN2R) == POLZN*NPHDIFF)) THEN
               CALL VMFDR(N1C,MN1C,NPHC,N2R,MN2R,NPHR,TERM)
               TERMLOC3=TERM*IETAC*ANORMC*ANORMR*2.0D0
               VL(MXLAM+1,IRC) = VL(MXLAM+1,IRC)
     &                         + TERMLOC3
         ENDIF
         ENDIF ! LROWDIF IF block

         IF (LCOLDIF .AND. LROWDIF) THEN
            ! <2C|Hm|2R>*krondel(1C,1R)
!           IF (JA1C.EQ.JA1R .AND. ABS(MN2C-MN2R).LE.POLZN) THEN
         ! This modified by J.D
         IF (JA1C.EQ.JA1R .AND. ((MN2C-MN2R) == POLZN*NPHDIFF)) THEN
               CALL VMFDR(N2C,MN2C,NPHC,N2R,MN2R,NPHR,TERM)
               TERMLOC4=TERM*ANORMC*ANORMR*2.0D0
               VL(MXLAM+1,IRC) = VL(MXLAM+1,IRC)
     &                         + TERMLOC4
         ENDIF

         ENDIF ! (LCOLDIF .AND. LROWDIF) IF block

      ENDDO ! IROW loop
      ENDDO ! ICOL loop: N basis loop ends here

      ! Diagonal contribution to H_int for electric field other than
      ! EFLDREF
!     DO ICOL = 1, N
!        IRC = ICOL+ICOL*(ICOL-1)/2
!        VL(NVLBLK,IRC) = SLOPE(JSINDX(ICOL))
!     ENDDO

      CALL GCLOCK(T2TIME)
      WRITE(6,*)
      WRITE(6,*)' CPL9 TIME:',T2TIME-T1TIME,' IN REST OF THE POTENTIAL',
     &          ' TERMS'

      ! Check whether any particular potential block has got all
      ! zeros or not.
      DO IHAM = 1, NVLBLK
         IF (ANY(ABS(VL(IHAM,:)) .GT. ZERO)) CYCLE
         WRITE(6,110) JTOT, IHAM
      ENDDO

  110 FORMAT('  * * * NOTE.  FOR JTOT =',I4,',  ALL COUPLING',
     &  ' COEFFICIENTS ARE 0.0 FOR POTENTIAL SYMMETRY',I4)

!     STOP 123

      RETURN
      END SUBROUTINE CPL9
!=========================================================================
      SUBROUTINE POTIN9(ITYPP, LAM, MXLAM, NPTS, NDIM, XPT, XWT, MXPT,
     &                  IVMIN, IVMAX, L1MAX, L2MAX, MXLMB, XFN, MX,
     &                  IXFAC)

      IMPLICIT NONE

      INTEGER, INTENT(INOUT)          :: ITYPP, MXLAM

      INTEGER, INTENT(OUT)            :: LAM(*)

      INTEGER, INTENT(IN)             :: IVMIN, IVMAX, L1MAX, L2MAX

      INTEGER, INTENT(INOUT)          :: NDIM, NPTS(NDIM), IXFAC, MX

      DOUBLE PRECISION, INTENT(INOUT) :: XFN(*)

      DOUBLE PRECISION, INTENT(OUT)   :: XPT(MXPT,NDIM), XWT(MXPT,NDIM)

      INTEGER, INTENT(IN)             :: MXLMB, MXPT

      ! Main function of potin9 is to set itypp to one of the precoded
      ! values understood by potenl (1 to 8) and do nothing else.
      ! Here this rotor-rotor interaction potential expansion to be used.
      ITYPP = 3

      RETURN
      END SUBROUTINE POTIN9
!=========================================================================
      SUBROUTINE EFV9(IFVARY)

      USE EFVS, ONLY: EFV

      USE POTENTIAL, ONLY: NCONST, VCONST

      USE ROTORS_MW_MOD, ONLY: FLDREF, FSTARK

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: IFVARY

!     VCONST(NCONST) = FSTARK*(EFV(1) - FLDREF)

      RETURN
      END SUBROUTINE EFV9
C=========================================================================
!     SUBROUTINE THRSH9(IREF, MONQN, NQN, EREF, IPRINT)

!     USE ROTORS_MW_MOD, ONLY: NPLVL, JPLVL, IQN, POLZN, EVALREF,
!    &                ETRNS, FDLTA, NMNIDXMONBASIS, IDX_SELC_BIG,
!    &                NMONBASISSELC,NTHR,MNTHR,NPHMX,MTTHR,
!    &                EVALDRSDREF

!     USE PHYSICAL_CONSTANTS, ONLY: HZ_IN_INV_CM

!     IMPLICIT NONE

!     DOUBLE PRECISION, INTENT(OUT) :: EREF

!     INTEGER, INTENT(IN)           :: IREF, MONQN(NQN), NQN, IPRINT

!     INTEGER :: IA1T, IA2T, NPHT, NQN1
!     
!     ! Counter and other local variables 
!     INTEGER :: IIC,IIC1
!     ! Counter and other local variables 

!     INTEGER :: ICOL, IROW, IA1C, IA2C, NPHC, IA1R, IA2R, NPHR

!     INTEGER :: N1C, MN1C, N2C, MN2C, N1R, MN1R, N2R, MN2R

!     INTEGER :: ITREF, ITREF2, IC, IR

!     INTEGER :: IMONT1, IMONT2, IMONT3, IMONSL

!     INTEGER :: MTMIN,MTMAX

!     DOUBLE PRECISION :: SQRTH, ANORMC, ANORMR, TERM

!     DOUBLE PRECISION :: T1TIME, T2TIME

!     DOUBLE PRECISION, ALLOCATABLE :: HAM(:,:), EVEC(:,:), EVAL(:),
!    &                                 DIAG(:), HAM2(:,:)

!     INTEGER, ALLOCATABLE :: INDX(:)

!     INTEGER, ALLOCATABLE :: IPHT(:)
!     

!     INTEGER :: NBLMT

!     LOGICAL :: LFG, LCOLDIF, LROWDIF

!     LOGICAL :: CFAC=.FALSE.

!     NQN1 = NQN-1

!     CALL GCLOCK(T1TIME)

!     IF (MONQN(1).EQ.-99999) THEN
!       WRITE(6,*) ' *** ERROR - THRSH9 CALLED WITH MONQN UNSET'
!       STOP
!     ENDIF
!

!     EREF=EVALDRSDREF(MONQN(1))+EVALDRSDREF(MONQN(2))

!     WRITE(33,*)
!     WRITE(33,*)
!     WRITE(33,*)"===================================="
!     WRITE(33,*) "Threshold Energy: |0,0>|0>+|0,0>|0>"
!     WRITE(33,*) "|0,0>|0> ==> Dressed monomer state"
!     WRITE(33,*)"===================================="
!     WRITE(33,*)
!     WRITE(33,*)"EREF in Cm^-1:", EREF
!     WRITE(33,*)
!     WRITE(33,*)"EREF in MHz:", (EREF/(HZ_IN_INV_CM*1.0D6))

!     ! THIS PART IS INEFFECTIVE
!     IF(CFAC) THEN

!     ! This encapsulated part is written by Joy Dutta
!     WRITE(19,*) "MONQN:"
!     WRITE(19,*) 
!     WRITE(19,*) MONQN 
!     WRITE(19,*) 
!     WRITE(19,*) 
!     WRITE(19,*) "Sl no. of monomers in MONQN"
!     WRITE(19,*) "i.e., |0,0>|0,0> in big basis:" 
!     WRITE(19,*) 

!     WRITE(25,*)
!     WRITE(25,*)
!     WRITE(25,*) "TEST FOR STORED N and MN for ONLY" 
!     WRITE(25,*) "MICROWAVE FIELD CALCULATIONS"
!     WRITE(25,*)

!     DO ICOL=1,NMONBASISSELC
!     WRITE(25,*) ICOL,NTHR(ICOL),MNTHR(ICOL),IDX_SELC_BIG(ICOL)
!     ENDDO 

!     IF(.NOT. ALLOCATED(IPHT)) THEN
!     ALLOCATE(IPHT(NPHMX))
!     IPHT=-999
!     ENDIF

!     WRITE(25,*)
!     WRITE(25,*) "========================="
!     WRITE(25,*)    "Photon quantum no. "
!     WRITE(25,*) "========================="
!     WRITE(25,*)
!     IIC=0
!     DO ICOL=NPHMX,1,-1
!     IIC=IIC+1
!     IPHT(IIC)=ICOL-NPHMX
!     WRITE(25,*) IIC,IPHT(IIC)
!     ENDDO

!     MTMIN=MINVAL(MNTHR)+MINVAL(POLZN*IPHT)
!     MTMAX=MAXVAL(MNTHR)+MAXVAL(POLZN*IPHT)


!     WRITE(25,*)
!     WRITE(25,*) "========================="
!     WRITE(25,*) "MIN & MAX VAL of MT ARRAY"
!     WRITE(25,*) "========================="
!     WRITE(25,*)
!     WRITE(25,*) MTMIN, MTMAX


!     NBLMT=0
!     DO IIC1=MTMIN,MTMAX
!     NBLMT=NBLMT+1
!     ENDDO


!     WRITE(25,*)
!     WRITE(25,*) "========================="
!     WRITE(25,*) "TOTAL NUM OF MT BLOCK"
!     WRITE(25,*) "========================="
!     WRITE(25,*)
!     WRITE(25,*) NBLMT


!     IF(.NOT. ALLOCATED(MTTHR)) THEN
!     ALLOCATE(MTTHR(NBLMT))
!     MTTHR=-999
!     ENDIF

!     IIC=0
!     DO IIC1=MTMIN,MTMAX
!     IIC=IIC+1
!     MTTHR(IIC)=IIC1
!     ENDDO

!     WRITE(25,*)
!     WRITE(25,*) "========================="
!     WRITE(25,*) " MT VALUE FOR EACH BLOCK"
!     WRITE(25,*) "========================="
!     WRITE(25,*)
!     WRITE(25,*) (MTTHR(IIC),IIC=1,NBLMT)
!     
!     ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     ! THIS PART IS INEFFECTIVE
!     IF(CFAC) THEN

!     IMONT1=NMNIDXMONBASIS(MONQN(1),MONQN(2)) 
!     IMONT2=NMNIDXMONBASIS(MONQN(3),MONQN(4))
!     IMONT3=MONQN(5)

!     WRITE(19,*) IMONT1, IMONT2, IMONT3 
!     WRITE(19,*) 
!     

!     N1C  = IQN(IMONT1,1)
!     MN1C = IQN(IMONT1,2)
!     
!     N2C  = IQN(IMONT2,1)
!     MN2C = IQN(IMONT2,2)
!     
!     WRITE(19,*) "Revival of quantum numbers:" 
!     WRITE(19,*) "N1C, MN1C, N2C, MN2C:"
!     WRITE(19,*) 
!     WRITE(19,*) N1C, MN1C, N2C, MN2C
!     WRITE(19,*) 
!     WRITE(19,*) "Total number of pair state:" 
!     WRITE(19,*) 
!     WRITE(19,*) "NPLVL:", NPLVL 
!     WRITE(19,*) 
!     WRITE(19,*) 
!     ! This encapsulated part is written by Joy Dutta

!     DO ICOL=1,NPLVL

!        IA1C = JPLVL(NQN1*(ICOL-1)+1)
!        IA2C = JPLVL(NQN1*(ICOL-1)+2)
!        NPHC = JPLVL(NQN1*(ICOL-1)+3)

!        N1C  = IQN(IA1C+1,1)
!        MN1C = IQN(IA1C+1,2)

!        N2C  = IQN(IA2C+1,1)
!        MN2C = IQN(IA2C+1,2)

!       WRITE(19,'(7I3)') ICOL,IA1C+1,N1C,MN1C,IA2C+1,N2C,MN2C
         
!     IF((((IA1C+1).EQ.IMONT1).AND.
!    & ((IA2C+1).EQ.IMONT2)).OR.
!    & (((IA1C+1).EQ.IMONT2).AND.
!    & ((IA2C+1).EQ.IMONT1))) THEN
!     IF(((IA1C+1).EQ.IMONT1).AND.
!    & ((IA2C+1).EQ.IMONT2)
!    & .AND.(IMONT3.EQ.NPHC)) THEN
!       IMONSL=ICOL
!     WRITE(19,*) "SL no. of |0,0>|0,0>|0>:", IMONSL
!     WRITE(19,*) 
!     ENDIF

!     ENDDO

!     ! Evaluation of diagonal energies

!     ALLOCATE(HAM(NPLVL,NPLVL))
!     HAM = 0.D0

!     DO ICOL = 1, NPLVL

!        IA1C = JPLVL(NQN1*(ICOL-1)+1)
!        IA2C = JPLVL(NQN1*(ICOL-1)+2)
!        NPHC = JPLVL(NQN1*(ICOL-1)+3)
!        LCOLDIF = .TRUE.
!        ANORMC  = SQRTH
!        IF (IA1C .EQ. IA2C) THEN
!          LCOLDIF = .FALSE.
!          ANORMC  = 1.D0
!        ENDIF
!        ! Get actual quantum numbers from IQN array
!        N1C  = IQN(IA1C+1,1)
!        MN1C = IQN(IA1C+1,2)

!        N2C  = IQN(IA2C+1,1)
!        MN2C = IQN(IA2C+1,2)

!       DO IROW = ICOL, NPLVL

!        IA1R = JPLVL(NQN1*(IROW-1)+1)
!        IA2R = JPLVL(NQN1*(IROW-1)+2)
!        NPHR = JPLVL(NQN1*(IROW-1)+3)

!        LROWDIF = .TRUE.
!        ANORMR  = SQRTH
!        IF (IA1R .EQ. IA2R) THEN
!          LROWDIF = .FALSE.
!          ANORMR  = 1.D0
!        ENDIF

!        ! Get actual quantum numbers from IQN array
!        N1R  = IQN(IA1R+1,1)
!        MN1R = IQN(IA1R+1,2)

!        N2R  = IQN(IA2R+1,1)
!        MN2R = IQN(IA2R+1,2)

!        ! Total projection of pair state should be conserved
!        IF (MN1C+MN2C+POLZN*NPHC.NE.MN1R+MN2R+POLZN*NPHR) CYCLE

!        ! Diagonal terms
!        IF (ICOL .EQ. IROW) THEN
!          ! Get the Stark + rotation + photon energy
!          HAM(IROW,ICOL) = EVALREF(IA1C+1)+EVALREF(IA2C+1)
!    &                    + (ETRNS+FDLTA)*NPHC

!       WRITE(21,'(6I3,f12.6)') ICOL, N1C, MN1C, N2C, MN2C, NPHC,
!    &    (HAM(IROW,ICOL)/(HZ_IN_INV_CM*1.0D6))             

!        ENDIF

!       ENDDO  ! IROW LOOP
!      ENDDO   ! ICOL LOOP

!     ! Evaluation of diagonal energies


!     Printing and checking the diagonal energies
!     Printing and checking the diagonal energies

!     ALLOCATE(DIAG(NPLVL),INDX(NPLVL))
!     DIAG = 0.D0
!     INDX = -999
!     FORALL (ICOL=1:NPLVL) DIAG(ICOL) = HAM(ICOL,ICOL)
!     FORALL (ICOL=1:NPLVL) INDX(ICOL) = ICOL

!     ! Do a bubble sort on the diagonal elements of the Hamiltonian
!     CALL BUBBLE_SORT(NPLVL,DIAG,INDX)

!     WRITE(21,*) 
!     WRITE(21,*) "Printing INDX:" 
!     WRITE(21,*) 
!     WRITE(21,*) (INDX(ICOL),ICOL=1,NPLVL) 
!     WRITE(21,*) 
!     WRITE(21,*) 
!     WRITE(21,*) "Printing Sorted energy element and index"
!     WRITE(21,*) 
!     WRITE(21,*) "New Sl, New Energy, Old Sl, Old Energy" 
!     WRITE(21,*) 
!     WRITE(21,*) 

!     DO ICOL=1,NPLVL

!     WRITE(21,'(I5,1x,f12.6,1x,I5,1x,f12.6)') ICOL, 
!    &   (DIAG(ICOL)/(HZ_IN_INV_CM*1.0D6)),
!    &   INDX(ICOL), (HAM(ICOL,ICOL)/(HZ_IN_INV_CM*1.0D6))
!     ENDDO

!     DO ICOL=1,NPLVL
!        IA1C = JPLVL(NQN1*(INDX(ICOL)-1)+1)                        
!        IA2C = JPLVL(NQN1*(INDX(ICOL)-1)+2)                       
!        NPHC = JPLVL(NQN1*(INDX(ICOL)-1)+3)                      
!                                                          
!        N1C  = IQN(IA1C+1,1)                            
!        MN1C = IQN(IA1C+1,2)                           
!                                                      
!        N2C  = IQN(IA2C+1,1)                         
!        MN2C = IQN(IA2C+1,2)                        
!                                                   
!       WRITE(19,'(7I3)') ICOL,IA1C+1,N1C,MN1C,IA2C+1,N2C,MN2C                              
!                                                  
!     IF((((IA1C+1).EQ.IMONT1).AND.                                                         
!    & ((IA2C+1).EQ.IMONT2)).OR.                                                            
!    & (((IA1C+1).EQ.IMONT2).AND.                                                           
!    & ((IA2C+1).EQ.IMONT1))) THEN                                                          
!     IF(((IA1C+1).EQ.IMONT1).AND.                
!    & ((IA2C+1).EQ.IMONT2)                      
!    & .AND.(IMONT3.EQ.NPHC)) THEN              
!       IMONSL=INDX(ICOL)                            
!     WRITE(21,*)                                         
!     WRITE(21,*) "New SL no:", ICOL,
!    & "Old SL no.:", IMONSL       
!     WRITE(21,*)                                         
!     ENDIF                                           

!     ENDDO

!     ENDIF ! CFAC CONDITION 



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     Inactive part
!     IF(CFAC) THEN  ! This part is now inactive

!     IA1T = MIN(MONQN(1),MONQN(2))
!     IA2T = MAX(MONQN(1),MONQN(2))
!     NPHT = MONQN(3)

!     IF (IA1T.LT.0 .OR. IA2T.LT.0)
!    & STOP 'THRSH9: IILEGAL VALUE(S) PASSED IN MONQN'

!     NQN1 = NQN-1

!     ALLOCATE(HAM(NPLVL,NPLVL))
!     HAM = 0.D0

!     SQRTH = SQRT(0.5D0)
!     LFG = .FALSE.

!     DO ICOL = 1, NPLVL

!        IA1C = JPLVL(NQN1*(ICOL-1)+1)
!        IA2C = JPLVL(NQN1*(ICOL-1)+2)
!        NPHC = JPLVL(NQN1*(ICOL-1)+3)

!        IF(IA1C.EQ.IA1T.AND.IA2C.EQ.IA2T.AND.
!    &      NPHC.EQ.NPHT) THEN
!           LFG = .TRUE.
!           ITREF = ICOL
!        ENDIF

!        LCOLDIF = .TRUE.
!        ANORMC  = SQRTH
!        IF (IA1C .EQ. IA2C) THEN
!          LCOLDIF = .FALSE.
!          ANORMC  = 1.D0
!        ENDIF

!        ! Get actual quantum numbers from IQN array
!        N1C  = IQN(IA1C+1,1)
!        MN1C = IQN(IA1C+1,2)

!        N2C  = IQN(IA2C+1,1)
!        MN2C = IQN(IA2C+1,2)

!     DO IROW = ICOL, NPLVL

!        IA1R = JPLVL(NQN1*(IROW-1)+1)
!        IA2R = JPLVL(NQN1*(IROW-1)+2)
!        NPHR = JPLVL(NQN1*(IROW-1)+3)

!        LROWDIF = .TRUE.
!        ANORMR  = SQRTH
!        IF (IA1R .EQ. IA2R) THEN
!          LROWDIF = .FALSE.
!          ANORMR  = 1.D0
!        ENDIF

!        ! Get actual quantum numbers from IQN array
!        N1R  = IQN(IA1R+1,1)
!        MN1R = IQN(IA1R+1,2)

!        N2R  = IQN(IA2R+1,1)
!        MN2R = IQN(IA2R+1,2)

!        ! Total projection of pair state should be conserved
!        IF (MN1C+MN2C+POLZN*NPHC.NE.MN1R+MN2R+POLZN*NPHR) CYCLE

!        ! Diagonal terms
!        IF (ICOL .EQ. IROW) THEN
!          ! Get the Stark + rotation + photon energy
!          HAM(IROW,ICOL) = EVALREF(IA1C+1)+EVALREF(IA2C+1)
!    &                    + (ETRNS+FDLTA)*NPHC
!        ENDIF

!        IF (ABS(NPHC-NPHR).NE.1) CYCLE
!                
!        ! Hm: Monomer Hamiltonian
!        ! <1C|Hm|1R>*krondel(2C,2R)
!        IF (IA2C.EQ.IA2R .AND. ABS(MN1C-MN1R).LE.POLZN) THEN
!           CALL VMFDR(N1C,MN1C,NPHC,N1R,MN1R,NPHR,TERM)
!           HAM(IROW,ICOL) = TERM*ANORMC*ANORMR*2.D0
!        ENDIF

!        IF (LCOLDIF) THEN
!           ! <2C|Hm|1R>*krondel(1C,2R)
!           IF (IA1C.EQ.IA2R .AND. ABS(MN2C-MN1R).LE.POLZN) THEN
!              CALL VMFDR(N2C,MN2C,NPHC,N1R,MN1R,NPHR,TERM)
!              HAM(IROW,ICOL) = HAM(IROW,ICOL)+TERM*ANORMC*ANORMR*2.D0
!           ENDIF
!        ENDIF

!        IF (LROWDIF) THEN
!           ! <1C|Hm|2R>*krondel(2C,1R)
!           IF (IA2C.EQ.IA1R .AND. ABS(MN1C-MN2R).LE.POLZN) THEN
!              CALL VMFDR(N1C,MN1C,NPHC,N2R,MN2R,NPHR,TERM)
!              HAM(IROW,ICOL) = HAM(IROW,ICOL)+TERM*ANORMC*ANORMR*2.D0
!           ENDIF
!        ENDIF ! (LCOLDIF .OR. LROWDIF) IF block

!        IF (LCOLDIF .AND. LROWDIF) THEN
!           ! <2C|Hm|2R>*krondel(1C,1R)
!           IF (IA1C.EQ.IA1R .AND. ABS(MN2C-MN2R).LE.POLZN) THEN
!              CALL VMFDR(N2C,MN2C,NPHC,N2R,MN2R,NPHR,TERM)
!              HAM(IROW,ICOL) = HAM(IROW,ICOL)+TERM*ANORMC*ANORMR*2.D0
!           ENDIF

!        ENDIF ! (LCOLDIF .AND. LROWDIF) IF block

!     ENDDO  ! IROW loop
!     ENDDO  ! ICOL loop

!     IF (.NOT. LFG) THEN
!       WRITE(6,*)'STOP: THRSHRSTRCT: SUPPLIED THRESHOLD IN MONQN DOES',
!    &            'NOT BELONG TO ONE OF THE PAIR THRESHOLDS.'
!       STOP
!     ENDIF

!     ALLOCATE(DIAG(NPLVL),INDX(NPLVL))
!     DIAG = 0.D0
!     INDX = -999

!     DO ICOL = 1, NPLVL
!        INDX(ICOL) = ICOL
!        DO IROW = 1,ICOL-1
!          HAM(IROW,ICOL) = HAM(ICOL,IROW)
!        ENDDO
         !write(*,'(7e16.6)')(ham(irow,icol),irow=1,nplvl)
!     ENDDO

!     ! Sort diagonal elements of the Hamiltonian in increasing order
!     ! so that the eigenvalues have one-to-one correspondence to the
!     ! basis indices
!     FORALL (ICOL=1:NPLVL) DIAG(ICOL) = HAM(ICOL,ICOL)

!     ! Do a bubble sort on the diagonal elements of the Hamiltonian
!     CALL BUBBLE_SORT(NPLVL,DIAG,INDX)

!     ! Construct the re-ordered Hamiltonian
!     ALLOCATE(HAM2(NPLVL,NPLVL))
!     HAM2 = 0.D0

!     DO ICOL = 1, NPLVL
!        IC = INDX(ICOL)
!        IF (IC .EQ. ITREF) ITREF2 = ICOL
!        DO IROW = 1, NPLVL
!           IR = INDX(IROW)
!           HAM2(IROW,ICOL) = HAM(IR,IC)
!        ENDDO
!        !write(*,'(7e16.6)')(ham2(irow,icol),irow=1,nplvl)
!     ENDDO

!     ALLOCATE(EVEC(NPLVL,NPLVL))
!     ALLOCATE(EVAL(NPLVL))

!     EVEC = 0.D0
!     EVAL = 0.D0

!     ! Diagonalize the re-ordered Hamiltonian
!     CALL DIAGVC(HAM2,NPLVL,NPLVL,EVAL,EVEC)

!     ! The quantum numbers of the threshold are poorly defined if the
      ! eigen vectors are mixed more than 50%.
!     IF ((EVEC(ITREF2,ITREF2)**2).LT.0.50) THEN
!        WRITE(6,*) "WARNING: IN THRSH9 SUPPLIED MONQN Q. NOS. ARE NO",
!    &              "LONGER GOOD TO DESCRIBE THE INCOMING THRESHOLD."
!     ENDIF 

!     EREF = EVAL(ITREF2)

      !print*,"itref",itref
!     !do icol=1,nplvl
!     !write(*,*)eval(icol)
!     !enddo

!     DEALLOCATE(HAM,HAM2)
!     DEALLOCATE(EVEC,EVAL)
!     DEALLOCATE(DIAG,INDX)

!     ENDIF

!     RETURN
!     END SUBROUTINE THRSH9
!=========================================================================
      SUBROUTINE EAC_DIP_DETER
      USE ROTORS_MW_MOD, ONLY:  POLZN,INT_NC,INT_MNC,INT_NPHC,
     &        INT_NR,INT_MNR,INT_NPHR,LAMPH,POLZN,OMEGA,FOMGA,
     &        EAC_D_IN_MHZ,EAC_D_IN_INV_CM 
      IMPLICIT NONE

      DOUBLE PRECISION :: THREEJ, THRJ, RES1, RES2, RES3, RES4

      WRITE(15,*) "==========================================="
      WRITE(15,*) "Polarization:", POLZN  
      WRITE(15,*) "==========================================="
      WRITE(15,*) "Info of states of given exp. Rabi frequency"
      WRITE(15,*)
      WRITE(15,*) "States in the column:"
      WRITE(15,*)
      WRITE(15,*) "INT_NC:",  INT_NC
      WRITE(15,*) "INT_MNC:",  INT_MNC
      WRITE(15,*) "INT_NPHC:",  INT_NPHC
      WRITE(15,*)
      WRITE(15,*) "States in the row:"
      WRITE(15,*)
      WRITE(15,*) "INT_NR:",  INT_NR
      WRITE(15,*) "INT_MNR:",  INT_MNR
      WRITE(15,*) "INT_NPHR:",  INT_NPHR
      WRITE(15,*)
      WRITE(15,*) "==========================================="
      WRITE(15,*) "Angular momentum of photon"
      WRITE(15,*)
      WRITE(15,*) "LAMPH:", LAMPH

      RES1=THREEJ(INT_NC, LAMPH, INT_NR)

      RES2= DSQRT((2.0D0*DBLE(INT_NC)+1.0D0)*(2.0D0*DBLE(INT_NR)+1.0D0))

      IF(INT_NPHC == (INT_NPHR-1)) THEN
      
      RES3=-DBLE(POLZN)*THRJ(DBLE(INT_NC), DBLE(LAMPH), DBLE(INT_NR),
     &                    -DBLE(INT_MNC), DBLE(POLZN), DBLE(INT_MNR))

      ELSE              
      RES3=0.0D0        
      ENDIF        

      IF(INT_NPHC == (INT_NPHR+1)) THEN

      RES3=DBLE(POLZN)*THRJ(DBLE(INT_NC), DBLE(LAMPH), DBLE(INT_NR),
     &                  -DBLE(INT_MNC), -DBLE(POLZN), DBLE(INT_MNR))

      ELSE
      RES3=0.0D0
      ENDIF        

      RES4=RES1*RES2*RES3

      IF(MOD(INT_MNC,2) /= 0) RES4=-RES4

      RES4=-RES4                ! From the definition


      EAC_D_IN_MHZ=OMEGA/RES4

      EAC_D_IN_INV_CM=FOMGA/RES4



      WRITE(15,*)
      WRITE(15,*) "======================================"
      WRITE(15,*)
      WRITE(15,*) 'Input Rabi frequency in MHz:', OMEGA
      WRITE(15,*)
      WRITE(15,*) 'Eac*d in MHz:',  EAC_D_IN_MHZ
      WRITE(15,*)
      WRITE(15,*) 'Evaluating Rabi frequncy back:',
     &                 EAC_D_IN_MHZ*RES4
      WRITE(15,*)
      WRITE(15,*) "*********************************"
      WRITE(15,*)
      WRITE(15,*) 'Input Rabi frequency in Cm^-1:', FOMGA
      WRITE(15,*)
      WRITE(15,*) 'Eac*d in Cm^-1:',  EAC_D_IN_INV_CM
      WRITE(15,*)
      WRITE(15,*) 'Evaluating Rabi frequncy in Cm^-1 back:',
     &                 EAC_D_IN_INV_CM*RES4
      WRITE(15,*)
      WRITE(15,*) "======================================"


      END SUBROUTINE        
!=========================================================================
      SUBROUTINE MONOMER_MAPPING

      USE ROTORS_MW_MOD, ONLY: NMAX, NMONBASIS, NBLMN, IQN, JQN, NQNS,
     & NMAXSELC, NMONBASISSELC, NDIM, MNIDX, NNIDX, NMNIDXMONBASIS,
     & IDX_SELC_BIG,NTHR,MNTHR

      IMPLICIT NONE

      INTEGER :: IFUNC, KK, NN, II, NTOT

      INTEGER :: ISELC, JSELC, IISELC

      INTEGER, ALLOCATABLE :: MNN(:)

      ! Total number of monomer blocks of given mn.
      NBLMN = 2*NMAX+1

      IF(.NOT. ALLOCATED(MNN)) ALLOCATE(MNN(NBLMN))
      MNN = -999

      ! Total number of monomer n quantum nos.
      NTOT = NMAX+1
      MNN(1) = 0
      KK = 1
      DO II = 1, NMAX
         MNN(KK+1) = -II
         MNN(KK+2) = II
         KK = KK+2
      ENDDO

      IF(.NOT. ALLOCATED(IQN)) ALLOCATE(IQN(NMONBASIS,NQNS-1))
      IQN = -999
      IF(.NOT. ALLOCATED(JQN)) ALLOCATE(JQN(2*NTOT,NBLMN))
      JQN = -999
      IF(.NOT. ALLOCATED(NDIM)) ALLOCATE(NDIM(NBLMN))
      NDIM = -999
      IF(.NOT. ALLOCATED(NNIDX)) ALLOCATE(NNIDX(0:NMAX,NBLMN))      
      NNIDX = -999
      IF(.NOT. ALLOCATED(MNIDX)) ALLOCATE(MNIDX(-NMAX:NMAX))      
      MNIDX = -999
      IF(.NOT. ALLOCATED(NMNIDXMONBASIS)) THEN
      ALLOCATE(NMNIDXMONBASIS(0:NMAX,-NMAX:NMAX))      
      NMNIDXMONBASIS = -999
      ENDIF

      WRITE(15,*) 
      WRITE(15,*)
      WRITE(15,*) "=================================================="
      WRITE(15,*) 
      WRITE(15,*) "Testing new array NMNIDXMONBASIS"
      WRITE(15,*) 


      WRITE(15,*) "****************************************************"
      WRITE(15,*) 
      WRITE(15,*) "Sl No., MN, N, N index in block, Sl No. using N, MN" 
      WRITE(15,*) 
      WRITE(15,*) "****************************************************"
      WRITE(15,*) 
      WRITE(15,*) 
      ! Store the quantum no. n,mn in arrays IQN and JQN to be used
      ! later.
      ! Store dimension of each block in array NDIM.
      IFUNC = 0
      DO II = 1, NBLMN
         KK = 0
         DO NN = 0, NMAX
            IF (NN .LT. ABS(MNN(II))) CYCLE
            IFUNC = IFUNC+1
            IQN(IFUNC,1) = NN
            IQN(IFUNC,2) = MNN(II)
            KK = KK+1
            JQN(KK,II) = NN
            JQN(KK+NTOT,II) = MNN(II)
            NNIDX(NN,II) = KK
            NMNIDXMONBASIS(NN,MNN(II))=IFUNC ! Added by Joy Dutta
            ! The following write statement is added by Joy Dutta
            WRITE(15,*) IFUNC, MNN(II), NN, NNIDX(NN, II), 
     &      NMNIDXMONBASIS(NN, MNN(II)) 
            !
         ENDDO
            WRITE(15,*) "++++++++++++++++++++++++++++++++++++++++++++++"
         NDIM(II) = KK
         MNIDX(MNN(II)) = II
      ENDDO      
      ! The following write statements are added by Joy Dutta
      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) "=================================================="
      WRITE(15,*) 
      WRITE(15,*) 'NMAXSELC:', NMAXSELC 
      WRITE(15,*) 
      WRITE(15,*) "=================================================="
      !

      ! The following section is added by Joy Dutta
      NMONBASISSELC=0 
!     DO ISELC=0,NMAXSELC
      DO ISELC=1,NMAXSELC  ! For N=1<===>N=2 Transition: Joy Dutta
        DO JSELC=-NMAXSELC,NMAXSELC
        IF(ABS(JSELC) > ISELC) CYCLE
        NMONBASISSELC=NMONBASISSELC+1
        ENDDO
      ENDDO

      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) 'NMONBASISSELC:', NMONBASISSELC 
      WRITE(15,*) 
      WRITE(15,*) "=================================================="


      IF(.NOT. ALLOCATED(IDX_SELC_BIG))  THEN 
      ALLOCATE(IDX_SELC_BIG(NMONBASISSELC))      
      IDX_SELC_BIG = -999
      ENDIF

      IF(.NOT. ALLOCATED(MNTHR)) THEN
      ALLOCATE(MNTHR(NMONBASISSELC))
      MNTHR=-999
      ENDIF

      IF(.NOT. ALLOCATED(NTHR)) THEN
      ALLOCATE(NTHR(NMONBASISSELC))
      NTHR=-999
      ENDIF

      WRITE(15,*) 
      WRITE(15,*) "=================================================="
      WRITE(15,*) 
      WRITE(15,*) "Before defining IDX_SELC_BIG array" 
      WRITE(15,*) 
      WRITE(15,*) IDX_SELC_BIG 
      WRITE(15,*) 
      WRITE(15,*) "=================================================="
      WRITE(15,*) 
      WRITE(15,*) "After defining IDX_SELC_BIG array" 
      WRITE(15,*)  
      WRITE(15,*) "Mapping of selected N,MN to Big array serial Number" 
      WRITE(15,*)  
      WRITE(15,*)  

! We map newly selected monomer numbers to old monomer numbers.
! Old monomer numbers are already mapped to old monomer quantum numbers
! Therefore, we could map newly selected monomer numbers to old 
! quantum numbers. 
! This mapping is necessary as info about old big monomer details 
! are required when external electric field will be present  

      IISELC=0 
!     DO ISELC=0,NMAXSELC
      DO ISELC=1,NMAXSELC  ! For N=1<===>N=2 Transition: Joy Dutta
        DO JSELC=-NMAXSELC,NMAXSELC
        IF(ABS(JSELC) > ISELC) CYCLE
        IISELC=IISELC+1
        IDX_SELC_BIG(IISELC)=NMNIDXMONBASIS(ISELC,JSELC)
        MNTHR(IISELC)=JSELC       ! This is stored for THRSH9
        NTHR(IISELC)=ISELC        ! This is stored for THRSH9
        WRITE(15,*) ISELC, JSELC, IDX_SELC_BIG(IISELC) 
        WRITE(25,*) IISELC,NTHR(IISELC),MNTHR(IISELC),
     &              IDX_SELC_BIG(IISELC)
        ENDDO
      ENDDO

      WRITE(15,*) 
      WRITE(15,*) 
      WRITE(15,*) "=================================================="
      !---- The above encapsulated section is added by Joy Dutta
      

      IF(ALLOCATED(MNN)) DEALLOCATE(MNN)

      RETURN
      END SUBROUTINE MONOMER_MAPPING
!=========================================================================
      SUBROUTINE MONOMER_ENERGIES
      USE ROTORS_MW_MOD, ONLY: NPLVL, JPLVL, IQN, POLZN, EVALREF,
     &                ETRNS, FDLTA, NMNIDXMONBASIS, IDX_SELC_BIG,
     &                NMAXSELC,NMONBASISSELC,NTHR,MNTHR,NPHMX,
     &                QMTTL,QNBTL,QMNBTL,QNPHBTL,MTTHR,
     &                EVALDRSDREF,EVECDRSDREF,NMONPHBASIS

      USE PHYSICAL_CONSTANTS, ONLY: HZ_IN_INV_CM

      INTEGER :: ICOL, IROW, NN, IFUNC, IIA, JJA   ! Local loop variables 

      INTEGER :: IIC,IIC1,IBLMT,ICPHT,ICNT,ICMONT,NDBL     ! Local variables

      INTEGER :: IELBLMT ! Local variables

      INTEGER :: MTMIN,MTMAX

      INTEGER, ALLOCATABLE :: IPHT(:)     ! Array of photon Quantum number 
      
!     INTEGER, ALLOCATABLE :: MTTHR(:)    ! Array of proj of rot Quantum number 

      INTEGER :: NBLMT                    ! Total number of MT blocks

      INTEGER :: NTL, MNTL                        

      INTEGER :: NTOTSELC


      INTEGER, ALLOCATABLE :: NDBLMT(:)   ! Array of dimensions of each MT block

      INTEGER, ALLOCATABLE :: MTIDX(:)    ! MT values connect to MT block indices

      ! These are the quantum numbers stored in function of ph-dressed 
      ! monomer basis
      INTEGER, ALLOCATABLE :: QNTL(:), QMNTL(:), QNPHTL(:) 

!     INTEGER, ALLOCATABLE :: QNBTL(:,:),QMNBTL(:,:),QNPHBTL(:,:)

      INTEGER, ALLOCATABLE :: QNMNIDX(:,:,:) 

      DOUBLE PRECISION, ALLOCATABLE :: HMONTOT(:,:) 

      DOUBLE PRECISION, ALLOCATABLE :: EVALEMW(:),EVECEMW(:,:)
      ! This encapsulated part is written by Joy Dutta                                      
!     WRITE(19,*) "MONQN:"                                              
!     WRITE(19,*)                                                      
!     WRITE(19,*) MONQN                                               
!     WRITE(19,*)                                                    
!     WRITE(19,*)                                                   
!     WRITE(19,*) "Sl no. of monomers in MONQN"                    
!     WRITE(19,*) "i.e., |0,0>|0,0> in big basis:"               
!     WRITE(19,*)                                               
!                                                              
      WRITE(25,*)                                             
      WRITE(25,*)                                         
      WRITE(25,*) "TEST FOR STORED N and MN for ONLY"        
      WRITE(25,*) "MICROWAVE FIELD CALCULATIONS"            
      WRITE(25,*)                                          
                                                         
      DO ICOL=1,NMONBASISSELC                           
      WRITE(25,*) ICOL,NTHR(ICOL),MNTHR(ICOL),IDX_SELC_BIG(ICOL)   
      ENDDO                                                       
                                                                 
      IF(.NOT. ALLOCATED(IPHT)) THEN                            
      ALLOCATE(IPHT(NPHMX))                                    
      IPHT=-999                                               
      ENDIF                                                 
                                                           
      WRITE(25,*)                                         
      WRITE(25,*) "========================="            
      WRITE(25,*)    "Photon quantum no. "              
      WRITE(25,*) "========================="          
      WRITE(25,*)                                    
      IIC=0                                         
      DO ICOL=NPHMX,1,-1                           
!     DO ICOL=-NPHMX,0,+1                           
      IIC=IIC+1                                   
      IPHT(IIC)=ICOL-NPHMX                                  
!     IPHT(IIC)=ICOL!-NPHMX                                  
      WRITE(25,*) IIC,IPHT(IIC)                            
      ENDDO                                                
                                                           
      MTMIN=MINVAL(MNTHR)+MINVAL(POLZN*IPHT)               
      MTMAX=MAXVAL(MNTHR)+MAXVAL(POLZN*IPHT)               
                                                          
                                                         
      WRITE(25,*)                                       
      WRITE(25,*) "========================="          
      WRITE(25,*) "MIN & MAX VAL of MT ARRAY"         
      WRITE(25,*) "========================="        
      WRITE(25,*)                                   
      WRITE(25,*) MTMIN, MTMAX                     


      NBLMT=0                                                          
      DO IIC1=MTMIN,MTMAX                                             
      NBLMT=NBLMT+1                                                  
      ENDDO                                                         
                                                                   
                                                                  
      WRITE(25,*)                                               
      WRITE(25,*) "========================="                 
      WRITE(25,*) "TOTAL NUM OF MT BLOCK"                   
      WRITE(25,*) "========================="                
      WRITE(25,*)                                                      
      WRITE(25,*) NBLMT                                               
                                                                      
                                                                      
      IF(.NOT. ALLOCATED(MTTHR)) THEN                                 
      ALLOCATE(MTTHR(NBLMT))                                         
      MTTHR=-999                                                     
      ENDIF                                                          
                                                                     
      IIC=0                                                          
      DO IIC1=MTMIN,MTMAX                                            
      IIC=IIC+1                                                      
      MTTHR(IIC)=IIC1                                                
      ENDDO                                                          
                                                                     
      WRITE(25,*)                                                    
      WRITE(25,*) "========================="                        
      WRITE(25,*) " MT VALUE FOR EACH BLOCK"                         
      WRITE(25,*) "========================="                        
      WRITE(25,*)                                                    
      WRITE(25,*) (MTTHR(IIC),IIC=1,NBLMT)                          

      ! NOW MAPPING OF MONOMER STARTS
      ! WE WILL HAVE EVERY INFO OF VARIOUS MONOMER STATES 
      ! WE WILL ALSO HAVE MAPPING INDEX

      NTOTSELC=(NMAXSELC+1)*(NMAXSELC+1)

      NMONPHBASIS=NTOTSELC*(NPHMX)  ! Total no of Ph-dressed monomer basis

      IF(.NOT. ALLOCATED(NDBLMT)) THEN                                 
      ALLOCATE(NDBLMT(NBLMT))                                         
      NDBLMT=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(MTIDX)) THEN                                 
      ALLOCATE(MTIDX(MTMIN:MTMAX))                                     
      MTIDX=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(QMTTL)) THEN                                 
      ALLOCATE(QMTTL(NMONPHBASIS))                                     
      QMTTL=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(QNTL)) THEN                                 
      ALLOCATE(QNTL(NMONPHBASIS))                                     
      QNTL=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(QMNTL)) THEN                                 
      ALLOCATE(QMNTL(NMONPHBASIS))                                     
      QMNTL=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(QNPHTL)) THEN                                 
      ALLOCATE(QNPHTL(NMONPHBASIS))                                     
      QNPHTL=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(QNBTL)) THEN                                 
      ALLOCATE(QNBTL(NMONPHBASIS,NBLMT)) 
      QNBTL=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(QMNBTL)) THEN                                 
      ALLOCATE(QMNBTL(NMONPHBASIS,NBLMT)) 
      QMNBTL=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(QNPHBTL)) THEN                                 
      ALLOCATE(QNPHBTL(NMONPHBASIS,NBLMT)) 
      QNPHBTL=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(QNMNIDX)) THEN                                 
      ALLOCATE(QNMNIDX(0:NMAXSELC,(1-NPHMX):0,NBLMT)) 
      QNPHBTL=-999                                                     
      ENDIF                                                          

      IF(.NOT. ALLOCATED(HMONTOT)) THEN                                 
      ALLOCATE(HMONTOT(NMONPHBASIS,NMONPHBASIS)) 
      HMONTOT=0.0D0
      ENDIF                                                          

      IF(.NOT. ALLOCATED(EVALEMW)) THEN                                 
      ALLOCATE(EVALEMW(NMONPHBASIS)) 
      EVALEMW=0.0D0
      ENDIF                                                          

      IF(.NOT. ALLOCATED(EVECEMW)) THEN                                 
      ALLOCATE(EVECEMW(NMONPHBASIS,NMONPHBASIS)) 
      EVECEMW=0.0D0
      ENDIF                                                          

      IF(.NOT. ALLOCATED(EVECDRSDREF)) THEN
      ALLOCATE(EVECDRSDREF(NMONPHBASIS,NMONPHBASIS,NBLMT)) 
      EVECDRSDREF=0.0D0
      ENDIF                                                          

      IF(.NOT. ALLOCATED(EVALDRSDREF)) THEN
      ALLOCATE(EVALDRSDREF(NMONPHBASIS)) 
      EVALDRSDREF=0.0D0
      ENDIF                                                          

      WRITE(25,*)                                                    
      WRITE(25,*) "========================="                        
      WRITE(25,*) "MN and IPH quantum no. for each MT block"
      WRITE(25,*) "========================="                        
      WRITE(25,*)                                                    
      WRITE(25,*) "mon sl. no., block indx, mt, n, mn, nph" 
      WRITE(25,*) 
      WRITE(25,*) 

      ICMONT=0                   ! a counter index
      DO IBLMT=1,NBLMT           ! loop over mt block index
       IELBLMT=0                 ! counter for number of elems in each mt block
       DO ICPHT=1,NPHMX          ! loop over Photon quantum no.
       MNTL=MTTHR(IBLMT) - (POLZN*IPHT(ICPHT))
        DO ICNT=0,NMAXSELC       ! loop over rot. quantum no.
        NTL=ICNT
        IF(ABS(MNTL) > NTL) CYCLE ! cond. of Quantum Mechanics
        ICMONT=ICMONT+1          ! counter for total number of mon basis
        IELBLMT=IELBLMT+1        ! counter for number of elems in each mt block

        QMTTL(ICMONT)=MTTHR(IBLMT)  !  Mon index ==> MT quantum no.
        QNTL(ICMONT)=NTL            !  Mon index ==> N quantum no. 
        QMNTL(ICMONT)=MNTL          !  Mon index ==> MN quantum no. 
        QNPHTL(ICMONT)=IPHT(ICPHT)  !  Mon index ==> PH quantum no. 

        QNBTL(IELBLMT,IBLMT)=NTL    ! Block index ==> N quantum no.
        QMNBTL(IELBLMT,IBLMT)=MNTL  ! Block index ==> MN quantum no.      
        QNPHBTL(IELBLMT,IBLMT)=IPHT(ICPHT)! Block index ==> NPH quan no.       

        QNMNIDX(QNTL(ICMONT),QNPHTL(ICMONT),IBLMT)=IELBLMT

        WRITE(25,'(6I4)') ICMONT,IBLMT,MTTHR(IBLMT),NTL,MNTL,IPHT(ICPHT)

        ENDDO                    ! end loop over rot. quantum no.
       ENDDO                     ! end of loop over Photon quantum no.
       NDBLMT(IBLMT)=IELBLMT     ! dimension of each MT block
       MTIDX(MTTHR(IBLMT))=IBLMT ! Index relate MT value to index of MT block
      ENDDO                  ! end of loop over mt block index 

      WRITE(25,*)
      WRITE(25,*) "---------------------------"
      WRITE(25,*) "Dimension of each mt block"
      WRITE(25,*) "---------------------------"
      WRITE(25,*)
      WRITE(25,*) "Index of each MT block"
      WRITE(25,*)
      WRITE(25,*) (IBLMT, IBLMT=1,NBLMT)
      WRITE(25,*)
      WRITE(25,*) "Corresponsing dimensions:"
      WRITE(25,*)
      WRITE(25,*) (NDBLMT(IBLMT),IBLMT=1,NBLMT)
      WRITE(25,*)
      WRITE(25,*) "MT value of each MT block"
      WRITE(25,*)
      WRITE(25,*) (MTTHR(IBLMT),IBLMT=1,NBLMT)
      WRITE(25,*)
      WRITE(25,*) "Corresponsing indices:"
      WRITE(25,*)
      WRITE(25,*) (MTIDX(MTTHR(IBLMT)),IBLMT=1,NBLMT)
      WRITE(25,*)
      WRITE(25,*)
      WRITE(25,*)


      WRITE(25,*)"===================================================="
      WRITE(25,*)"No. Ph-dressed basis, index of block, MT, N, MN, NPH"
      WRITE(25,*)"===================================================="
      WRITE(25,*)

      DO IIC=1,NMONPHBASIS 
      WRITE(25,'(6I5)') IIC,MTIDX(QMTTL(IIC)),QMTTL(IIC),QNTL(IIC),
     &            QMNTL(IIC),QNPHTL(IIC) 
      ENDDO


      WRITE(25,*)"===================================================="
      WRITE(25,*)"Accessing monomer details by looping over"
      WRITE(25,*)"total no. of MT blocks"
      WRITE(25,*)"===================================================="
      WRITE(25,*)

      DO IBLMT=1,NBLMT
       DO IIC1=1,NDBLMT(IBLMT)
       WRITE(25,*) IBLMT, IIC1, QNBTL(IIC1,IBLMT), 
     &  QMNBTL(IIC1,IBLMT), QNPHBTL(IIC1,IBLMT),
     &  QNMNIDX(QNBTL(IIC1,IBLMT),QNPHBTL(IIC1,IBLMT),IBLMT)
     
       ENDDO
      ENDDO

      IFUNC=1
      DO IIC=1,NBLMT
      NN=NDBLMT(IIC)
      IIA=IFUNC
      JJA=IFUNC+NN-1
      CALL HMONBLK(IIC,NN,HMONTOT)
      CALL DIAGVC(HMONTOT(:NN,:NN),NN,NN,EVALEMW(:NN),EVECEMW(:NN,:NN))
      EVECDRSDREF(:NN,:NN,IIC) = EVECEMW(:NN,:NN)
      EVALDRSDREF(IIA:JJA)     = EVALEMW(:NN)
      IFUNC=IFUNC+NN
      ENDDO

      DO IIC=1, NMONPHBASIS
      WRITE(33,*) IIC, (EVALDRSDREF(IIC)/(HZ_IN_INV_CM*1.0D6))
      ENDDO

      RETURN
      END SUBROUTINE MONOMER_ENERGIES
!=========================================================================
      SUBROUTINE HMONBLK(IBLK,NDIMBLMT,HMONTOT)
      
      USE ROTORS_MW_MOD, ONLY: QMTTL,QNBTL,QMNBTL,QNPHBTL, 
     &  NMONPHBASIS,FROT,FDLTA,POLZN,MTTHR,EAC_D_IN_INV_CM 
     
      USE PHYSICAL_CONSTANTS, ONLY: HZ_IN_INV_CM

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IBLK, NDIMBLMT

      INTEGER :: IIC1, ICOL, IROW, NCOL, NROW

      INTEGER :: NC, NR, MNC, MNR, NPHC, NPHR 

      DOUBLE PRECISION, INTENT(OUT) :: HMONTOT(NDIMBLMT,NDIMBLMT)

      LOGICAL :: MW_COND=.TRUE., EF_COND=.FALSE.

      DOUBLE PRECISION, ALLOCATABLE :: HDIAGMW(:,:)

      DOUBLE PRECISION, ALLOCATABLE :: HDIAGROT(:,:)

      DOUBLE PRECISION, ALLOCATABLE :: HDIAG(:,:)

      DOUBLE PRECISION, ALLOCATABLE :: HINTMW(:,:)

 
      DOUBLE PRECISION :: TERMP, TERMM

      IF(.NOT. ALLOCATED(HDIAGMW)) THEN                                 
      ALLOCATE(HDIAGMW(NDIMBLMT,NDIMBLMT)) 
      HDIAGMW=0.0D0
      ENDIF                                                          

      IF(.NOT. ALLOCATED(HDIAGROT)) THEN    
      ALLOCATE(HDIAGROT(NDIMBLMT,NDIMBLMT)) 
      HDIAGROT=0.0D0
      ENDIF                                                          

      IF(.NOT. ALLOCATED(HDIAG)) THEN    
      ALLOCATE(HDIAG(NDIMBLMT,NDIMBLMT)) 
      HDIAG=0.0D0
      ENDIF                                                          

      IF(.NOT. ALLOCATED(HINTMW)) THEN    
      ALLOCATE(HINTMW(NDIMBLMT,NDIMBLMT)) 
      HINTMW=0.0D0
      ENDIF                                                          

!     DO IIC1=1,NDIMBLMT
!     WRITE(26,*) IBLK,IIC1,QNBTL(IIC1,IBLK),QMNBTL(IIC1,IBLK),
!    &           QNPHBTL(IIC1,IBLK)
!     ENDDO

      NCOL=0
      DO ICOL=1,NDIMBLMT
      NCOL=NCOL+1
 
      NC=QNBTL(NCOL,IBLK)
      MNC=QMNBTL(NCOL,IBLK)
      NPHC=QNPHBTL(NCOL,IBLK)

!     WRITE(26,*) IBLK, ICOL, NC, MNC, NPHC

       NROW=0
       DO IROW=1,ICOL
       NROW=NROW+1

       NR=QNBTL(NROW,IBLK)
       MNR=QMNBTL(NROW,IBLK)
       NPHR=QNPHBTL(NROW,IBLK)
 
!      WRITE(27,*) IBLK, IROW, NR, MNR, NPHR

       IF(IROW==ICOL) THEN
       HDIAGROT(IROW,ICOL)=FROT*DBLE(NC*(NC+1))
       IF(MW_COND) THEN
       HDIAGMW(IROW,ICOL)=(2.0D0*FROT+FDLTA)*DBLE(NPHC)
       ELSE
       HDIAGMW(IROW,ICOL)=0.0D0
       ENDIF
       HDIAG(IROW,ICOL)=HDIAGROT(IROW,ICOL)+HDIAGMW(IROW,ICOL)
       WRITE(28,'(3I3, 3f12.6)') NR, NC, NPHC,
     & (HDIAGROT(IROW,ICOL)/(HZ_IN_INV_CM*1.0D6)), 
     & (HDIAGMW(IROW,ICOL)/(HZ_IN_INV_CM*1.0D6)), 
     & (HDIAG(IROW,ICOL)/(HZ_IN_INV_CM*1.0D6))
       ENDIF

       IF(MW_COND) THEN
       IF(IROW /= ICOL) THEN

       IF(MTTHR(IBLK) /= MTTHR(IBLK)) CYCLE

       IF (ABS(NC-NR) /= 1) CYCLE
       IF (ABS(NPHR-NPHC) /= 1) CYCLE

       WRITE(29,'(8I4)') NC, MNC, NPHC, MTTHR(IBLK), 
     & NR, MNR, NPHR, MTTHR(IBLK)

       WRITE(30,'(8I4)') NC, MNC, NPHC, MNC+POLZN*NPHC, 
     & NR, MNR, NPHR, MNR+POLZN*NPHR
!!!!!!!!!!!

        IF (NPHR-NPHC.EQ.1) THEN
           CALL MATCOUPLE(NC,MNC,NR,MNR,POLZN,TERMP)
        IF (POLZN.EQ.1) TERMP = -TERMP
        ELSE                !
           TERMP=0.0D0     ! By J.D
        ENDIF
        IF (NPHC-NPHR.EQ.1) THEN
           CALL MATCOUPLE(NC,MNC,NR,MNR,-POLZN,TERMM)
        IF (POLZN.EQ.-1) TERMM = -TERMM
        ELSE                !
           TERMM=0.0D0     ! By J.D
        ENDIF
           HINTMW(IROW,ICOL)=0.5D0*(-EAC_D_IN_INV_CM*(TERMP+TERMM)) ! Proper val
!          HINTMW(IROW,ICOL)=(-EAC_D_IN_INV_CM*(TERMP+TERMM)) ! Hutson and Karman
!      HINTMW(IROW,ICOL)=(1.0D0/SQRT(2.0D0))*
!    & (-EAC_D_IN_INV_CM*(TERMP+TERMM)) ! Proper val

       WRITE(31,'(6I4,f12.6)') NC, MNC, NPHC, 
     & NR, MNR, NPHR, (HINTMW(IROW,ICOL)/(HZ_IN_INV_CM*1.0D6)) 
!!!!!
       ENDIF ! ENDIF cond for IROW /= ICOL
       ENDIF ! ENDIF cond for MW_COND 

!!!!!  HINT for Electric Field Dressing is not inserted  
           HINTMW(ICOL,IROW)=HINTMW(IROW,ICOL)

       ENDDO
      ENDDO
      
      HMONTOT=HDIAG+HINTMW
     
      NCOL=0
      DO ICOL=1,NDIMBLMT
      NCOL=NCOL+1
 
      NC=QNBTL(NCOL,IBLK)
      MNC=QMNBTL(NCOL,IBLK)
      NPHC=QNPHBTL(NCOL,IBLK)

!     WRITE(26,*) IBLK, ICOL, NC, MNC, NPHC

      NROW=0
      DO IROW=1,NDIMBLMT
      NROW=NROW+1

      NR=QNBTL(NROW,IBLK)
      MNR=QMNBTL(NROW,IBLK)
      NPHR=QNPHBTL(NROW,IBLK)
 
      WRITE(32,'(7I4,f12.6)') IBLK, NC, MNC, NPHC, 
     & NR, MNR, NPHR, (HMONTOT(IROW,ICOL)/(HZ_IN_INV_CM*1.0D6)) 
      ENDDO
      ENDDO

      WRITE(32,*)
      WRITE(32,*)
      WRITE(32,*)

      RETURN 
      END SUBROUTINE HMONBLK
!=========================================================================
      SUBROUTINE BASIS_BUILDER

      ! Constructs the pair energies and eigenvectors for dressed basis
!     USE ROTORS_MW_MOD, ONLY: NMONBASIS, FLDREF, NBLMN, IQN, NDIM,
!    &                         EVECREF, EVALREF, LRSTRCT, MPAIR, IMLST, 
!    &                         NMAX, IPAIR, SLOPE, NPAIR, ITRNS, ETRNS,
!    &                         NMAXSELC,IDX_SELC_BIG,NMNIDXMONBASIS,
!    &                         FROT,NMONBASISSELC,MNIDX,NNIDX 
      USE ROTORS_MW_MOD, ONLY: NMONBASIS, FLDREF, NBLMN, IQN, NDIM,
     &                         EVECREF, EVALREF, LRSTRCT, MPAIR, IMLST, 
     &                         NMAX, IPAIR, NPAIR, ITRNS, ETRNS,
     &                         NMAXSELC,IDX_SELC_BIG,NMNIDXMONBASIS,
     &                         FROT,NMONBASISSELC,MNIDX,NNIDX 
      USE PHYSICAL_CONSTANTS, ONLY: HZ_IN_INV_CM

      IMPLICIT NONE

      DOUBLE PRECISION, PARAMETER :: DFLD = 0.1D0 ! in kV/cm

      DOUBLE PRECISION :: EFLD

      DOUBLE PRECISION, ALLOCATABLE :: EVALNBR(:), DIPSPFXD(:), EVAL(:),
     &                                 HMON(:,:),HSTARK(:,:),EVEC(:,:),
     &                                 DIP(:),EPAIR_REF(:), EPAIR_NBR(:)

      INTEGER :: IBL, II, JJ, NN, IFUNC, NSYM, IA1, IA2, IM1, IM2

      INTEGER :: IISELC, JJSELC, ISELC, JSELC, IMON1, IMON2

      INTEGER :: KK, IDX

      INTEGER, ALLOCATABLE :: IPAIR_IDX(:)

      ! Allocate relevant arrays
      IF (.NOT. ALLOCATED(EVALREF)) ALLOCATE(EVALREF(NMONBASIS))
      EVALREF = 0.D0
      IF (.NOT. ALLOCATED(EVALNBR)) ALLOCATE(EVALNBR(NMONBASIS))
      EVALNBR = 0.D0
      IF (.NOT. ALLOCATED(DIP)) ALLOCATE(DIP(NMONBASIS))
      DIP = 0.D0
      IF (.NOT. ALLOCATED(DIPSPFXD)) ALLOCATE(DIPSPFXD(NMONBASIS))
      DIPSPFXD = 0.D0
      IF (.NOT. ALLOCATED(HMON)) ALLOCATE(HMON(NMAX+1,NMAX+1))
      HMON = 0.D0
      IF (.NOT. ALLOCATED(EVECREF)) 
     &    ALLOCATE(EVECREF(NMAX+1,NMAX+1,NBLMN))
      EVECREF = 0.D0
      IF (.NOT. ALLOCATED(EVEC)) ALLOCATE(EVEC(NMAX+1,NMAX+1))
      EVEC = 0.D0
      IF (.NOT. ALLOCATED(EVAL)) ALLOCATE(EVAL(NMAX+1))
      EVAL = 0.D0
      IF (.NOT. ALLOCATED(HSTARK)) ALLOCATE(HSTARK(NMAX+1,NMAX+1))
      HSTARK = 0.D0

      ! First calculate the monomer eigenvalues and eigenvectors
      ! --------------------------------------------------------

      ! Reference field calculations
      EFLD  = FLDREF
      IFUNC = 1
      DO IBL = 1, NBLMN
        ! Size of this mn block of IBL 
        NN = NDIM(IBL)
        II = IFUNC
        JJ = IFUNC+NN-1
        IF ((IBL.EQ.1) .OR. (MOD(IBL,2).EQ.0)) THEN
           CALL HAMSTARK(IBL,NN,EFLD,HMON(:NN,:NN),HSTARK(:NN,:NN))
           CALL DIAGVC(HMON(:NN,:NN),NN,NN,EVAL(:NN),EVEC(:NN,:NN))
           CALL DIPOLE_SPFXD(NN,EVEC(:NN,:NN),HSTARK(:NN,:NN),DIP(:NN))
        ENDIF
        EVECREF(:NN,:NN,IBL) = EVEC(:NN,:NN)
        EVALREF(II:JJ)       = EVAL(:NN)
        DIPSPFXD(II:JJ)      = DIP(:NN)
        IFUNC = IFUNC+NN
      ENDDO

      ! Neighbouring field calculations
      EFLD  = FLDREF+DFLD
      IFUNC = 1
      DO IBL = 1, NBLMN
        NN = NDIM(IBL)
        II = IFUNC
        JJ = IFUNC+NN-1
        IF ((IBL.EQ.1) .OR. (MOD(IBL,2).EQ.0)) THEN
           CALL HAMSTARK(IBL,NN,EFLD,HMON(:NN,:NN),HSTARK(:NN,:NN))
           CALL DIAGVC(HMON(:NN,:NN),NN,NN,EVAL(:NN),EVEC(:NN,:NN))
        ENDIF
        EVALNBR(II:JJ) = EVAL(:NN)
        IFUNC = IFUNC+NN
      ENDDO

      WRITE(6,107)
 107  FORMAT(/"  EIGENVALUES OF MONOMER HAMILTONIAN AT"  ,
     &         " REFERENCE FIELD"                                  ,
     &       /"  -------------------------------------------------",
     &       /"        a      n      mn      ENERGY(cm-1)         ",
     &       /"  -------------------------------------------------")
      DO II = 1, NMONBASIS
         WRITE(6,'(2X,3I7,F18.8)')II-1,IQN(II,1),IQN(II,2),EVALREF(II)
      ENDDO

      WRITE(6,1017)
 1017  FORMAT(/"  EIGENVALUES OF MONOMER HAMILTONIAN AT"  ,
     &         " REFERENCE FIELD"                                  ,
     &       /"  -------------------------------------------------",
     &       /"        a      n      mn      ENERGY(MHz)         ",
     &       /"  -------------------------------------------------")
      DO II = 1, NMONBASIS
         WRITE(6,'(2X,3I7,F18.8)')II-1,IQN(II,1),IQN(II,2),
     &   EVALREF(II)/(HZ_IN_INV_CM*1.0D6)
      ENDDO

      ETRNS = EVALREF(ITRNS(2)+1)-EVALREF(ITRNS(1)+1)

      WRITE(6,108)
 108  FORMAT(/"  SPACE FIXED DIPOLES AT REFERENCE FIELD "
     &       /"  -------------------------------------------------",
     &       /"        n      mn        D(debye)                  ",
     &       /"  -------------------------------------------------")
      DO II = 1, NMONBASIS
         WRITE(6,'(2X,2I7,F18.8)')IQN(II,1),IQN(II,2),DIPSPFXD(II)
      ENDDO


      WRITE(6,*)
      WRITE(6,*) "====================================================="
      WRITE(6,*) "IDENTIFYING INTERACTING PAIR OF MONOMER STATES:"
      WRITE(6,*) "TARGET STATES: |1,1> and |2,2>"
      WRITE(6,*) "====================================================="
      WRITE(6,*) 
      WRITE(6,*) "-----------------------------------------------------" 
      WRITE(6,*) "                      n       mn      ENERGY(MHz)" 
      WRITE(6,*) "-----------------------------------------------------" 
      WRITE(6,*) 
      WRITE(6,*) ITRNS(2)+1, IQN(ITRNS(2)+1,1), IQN(ITRNS(2)+1,2), 
     &           EVALREF(ITRNS(2)+1)/(HZ_IN_INV_CM*1.0D6)
      WRITE(6,*)
      WRITE(6,*) ITRNS(1)+1, IQN(ITRNS(1)+1,1), IQN(ITRNS(1)+1,2), 
     &           EVALREF(ITRNS(1)+1)/(HZ_IN_INV_CM*1.0D6)
      WRITE(6,*) 
      WRITE(6,*) "====================================================="
      WRITE(6,*) "ENERGY DIFFERENCE (In MHz) BETWEEN: |2,2> and |1,1>"
      WRITE(6,*)
      WRITE(6,*) "ETRNS in MHz:", ETRNS/(HZ_IN_INV_CM*1.0D6)
      WRITE(6,*)
      WRITE(6,*) "NOTE: ENERGY DIFFERENCE = 4 B_rot"

!     The following part of the code is written by Joy Dutta 
!     A Check of new parameter NMAXSELC and IDX_SELC_BIG   
      WRITE(15,*) "=================================================="
      WRITE(15,*) "A check of NMAXSELC and IDX_SELC_BIG in Basis build"
      WRITE(15,*) "=================================================="
      IISELC=0                                                          
!     DO ISELC=0,NMAXSELC                                               
      DO ISELC=1,NMAXSELC        ! For N=1<==>N=2 Transition Case: JOY DUTTA                                   
        DO JSELC=-NMAXSELC,NMAXSELC                                     
        IF(ABS(JSELC) > ISELC) CYCLE                                    
        IISELC=IISELC+1                                                 
        JJSELC=NMNIDXMONBASIS(ISELC,JSELC)
      WRITE(15,'(4I3,2F12.6,6I3 )') ISELC, JSELC, IDX_SELC_BIG(IISELC),
     & JJSELC,EVALREF(JJSELC), 
     & FROT*DBLE(ISELC*(ISELC+1))/(HZ_IN_INV_CM*1.0D6), 
     & IQN(JJSELC,1), IQN(JJSELC,2), MNIDX(JSELC), NDIM(MNIDX(JSELC)),
     & NNIDX(IQN(JJSELC,1), MNIDX(JSELC))
        ENDDO                                                           
      ENDDO                                                             

                                                                        
      WRITE(15,*)                                                       
      WRITE(15,*)                                                       
      WRITE(15,*) "=================================================="
      WRITE(15,*)                                                       
      WRITE(15,*) "NMONBASISSELC:", NMONBASISSELC
      WRITE(15,*)                                                       
!     The above encapsulated part of the code is written by Joy Dutta 

     


      ! Calculation of monomer eigenvalues and eigenvectors finishes here
      ! ----------------------------------------------------------------

      ! ----------------------------------------------------------------
      ! Assigns monomer quantum numbers and energies to pair states
      ! ----------------------------------------------------------------
      ! Pair states selection can be done by two mechanisms:
      ! 1. All pair states are included (usually a large basis set will
      !    be created) - LRSTRCT must be false.
      ! 2. Subset of pair states are to be selected using NPAIR
      !    mechanism. NPAIR should be a positive integer specifying
      !    total number of required pair states. IPAIR then should
      !    contain sets of a1, a2 quantum numbers for those pair states.

      IF (.NOT. LRSTRCT) THEN    ! Building full pair basis functions
        
         WRITE(6,109)
 109     FORMAT(/"  PAIR HAMILTONIAN BASIS SET WILL BE CONSTRUCTED ",
     &           "USING FULL DIRECT PRODUCT OF MONOMER BASIS ",
     &           "FUNCTIONS.")
      ! This was the previous case, when full space of NMAX and 
      ! NMONBASIS were used. But here we select a particular space from
      ! the E-field dressed monomer states. For this reason, Joy Dutta
      ! blocked the following line
      !   NSYM  = NMONBASIS*(NMONBASIS+1)/2

      ! Instead the is modified as the following by Joy Dutta
         NSYM  = NMONBASISSELC*(NMONBASISSELC+1)/2
        
         IF (ALLOCATED(IPAIR)) DEALLOCATE(IPAIR)
         ALLOCATE(IPAIR(NSYM*2))
!        ALLOCATE(IPAIR(NSYM))
         IPAIR = -999
         IF (.NOT. ALLOCATED(EPAIR_REF)) ALLOCATE(EPAIR_REF(NSYM))
         EPAIR_REF = 0.D0
         IF (.NOT. ALLOCATED(EPAIR_NBR)) ALLOCATE(EPAIR_NBR(NSYM))
         EPAIR_NBR = 0.D0
!        IF (.NOT. ALLOCATED(SLOPE)) ALLOCATE(SLOPE(NSYM))
!        SLOPE = 0.D0

      ! By Joy Dutta
       WRITE(15,*) "Sl. no. of pair state, Sl. no. of mon1 in big set",
     & "Sl. no. of mon2 in big set, mon1 N, mon1 MN, mon2 N, mon2 MN"

         II = 0
         JJ = 0
      ! DO IA1 = 1, NMONBASIS
        DO IA1 = 1, NMONBASISSELC
            IMON1=IDX_SELC_BIG(IA1)  ! Added by Joy Dutta
      !     DO IA2 = IA1, NMONBASIS
            DO IA2 = IA1, NMONBASISSELC
            IMON2=IDX_SELC_BIG(IA2)  ! Added by Joy Dutta

               II  = II+1

      !        EPAIR_REF(II) = EVALREF(IA1)+EVALREF(IA2)
               EPAIR_REF(II) = EVALREF(IMON1)+EVALREF(IMON2) ! By J.D
      !        EPAIR_NBR(II) = EVALNBR(IA1)+EVALNBR(IA2)
               EPAIR_NBR(II) = EVALNBR(IMON1)+EVALNBR(IMON2) ! By J.D
      !        IPAIR(JJ+1) = IA1-1
               IPAIR(JJ+1) = IMON1-1     ! By J.D
      !        IPAIR(JJ+2) = IA2-1
               IPAIR(JJ+2) = IMON2-1     ! By J,D
       WRITE(15,'(7I3)') II, IMON1, IPAIR(JJ+1), IMON2, IPAIR(JJ+2) 
               JJ = JJ+2
    
            ENDDO
            WRITE(15,*) 
         ENDDO
!        SLOPE = (EPAIR_NBR-EPAIR_REF)/DFLD
         NPAIR = II
            WRITE(15,*) 
            WRITE(15,*) 
            WRITE(15,'(A,I5)') 'Total pair state:', NPAIR

      ! In the section following ELSE, I (J.D) have not 
      ! modified NMONBASIS with NMONBASISSELC
      ELSE    ! Building restricted basis set using NPAIR mechanism

         WRITE(6,112)
 112     FORMAT(/"  PAIR HAMILTONIAN BASIS SET WILL BE CONSTRUCTED ",
     &           "USING USER INPUT NPAIR AND IPAIR IN &BASIS9"
     &          /"  THROUGH NPAIR MECHANISM.")

         IF (.NOT. ALLOCATED(IPAIR_IDX)) ALLOCATE(IPAIR_IDX(NPAIR))
         IPAIR_IDX = -999
         IF (.NOT. ALLOCATED(EPAIR_REF)) ALLOCATE(EPAIR_REF(NPAIR))
         EPAIR_REF = 0.D0
         IF (.NOT. ALLOCATED(EPAIR_NBR)) ALLOCATE(EPAIR_NBR(NPAIR))
         EPAIR_NBR = 0.D0
!        IF (.NOT. ALLOCATED(SLOPE)) ALLOCATE(SLOPE(NPAIR))
!        SLOPE = 0.D0

         KK = 0
         DO II = 1, NPAIR
            IA1 = MIN(IPAIR(KK+1),IPAIR(KK+2))
            IA2 = MAX(IPAIR(KK+1),IPAIR(KK+2))
            IPAIR(KK+1) = IA1
            IPAIR(KK+2) = IA2
            IPAIR_IDX(II) = (IA1+1)+IA2*(IA2+1)/2
            EPAIR_REF(II) = EVALREF(IA1+1)+EVALREF(IA2+1)
            EPAIR_NBR(II) = EVALNBR(IA1+1)+EVALNBR(IA2+1)
            KK = KK+2
         ENDDO
!        SLOPE = (EPAIR_NBR-EPAIR_REF)/DFLD

         ! Number of functions excluded from the actual basis
         ! functions 
         MPAIR = NMONBASIS*(NMONBASIS+1)/2 - NPAIR
         IF (.NOT. ALLOCATED(IMLST)) ALLOCATE(IMLST(MPAIR,2))
         IMLST = -999

         ! Now store the information of the functions excluded from
         ! the restricted basis set
         KK = 0
         DO IM1 = 1, NMONBASIS
            DO IM2 = IM1, NMONBASIS
            IDX = IM1+IM2*(IM2-1)/2
            ! Inclusion of this function if tag is not present in
            ! the above list ('ipair_idx' array)
            IF (ANY(IPAIR_IDX-IDX .EQ. 0)) CYCLE
            KK = KK+1
            IMLST(KK,1) = IM1
            IMLST(KK,2) = IM2
            ENDDO
         ENDDO

      ENDIF ! IF (.NOT. LRSTRCT) block

      ! Purposes of the following arrays finished so deallocate
      ! Purposes of the following arrays finished so deallocate
      IF (ALLOCATED(EVALNBR)) DEALLOCATE(EVALNBR)
      IF (ALLOCATED(DIP)) DEALLOCATE(DIP)
      IF (ALLOCATED(DIPSPFXD)) DEALLOCATE(DIPSPFXD)
      IF (ALLOCATED(HMON)) DEALLOCATE(HMON)
      IF (ALLOCATED(HSTARK)) DEALLOCATE(HSTARK)
      IF (ALLOCATED(EVEC)) DEALLOCATE(EVEC)
      IF (ALLOCATED(EVAL)) DEALLOCATE(EVAL)
      IF (ALLOCATED(EPAIR_REF)) DEALLOCATE(EPAIR_REF)
      IF (ALLOCATED(EPAIR_NBR)) DEALLOCATE(EPAIR_NBR)
      IF (ALLOCATED(IPAIR_IDX)) DEALLOCATE(IPAIR_IDX)

      RETURN
      END SUBROUTINE BASIS_BUILDER
!=========================================================================
      SUBROUTINE HAMSTARK(IB, NHAM, EFLD, HAM, HSTARK)

      USE ROTORS_MW_MOD, ONLY: FROT, FSTARK, JQN, NMAX

      IMPLICIT NONE

      INTEGER, INTENT(IN)           :: IB, NHAM

      DOUBLE PRECISION, INTENT(IN)  :: EFLD

      DOUBLE PRECISION, INTENT(OUT) :: HAM(NHAM,NHAM), HSTARK(NHAM,NHAM)

      DOUBLE PRECISION              :: HROT(NHAM,NHAM)

      INTEGER                       :: IROW, ICOL, NC, MNC, NR, MNR, KC,
     &                                 KR

      ! Functions
      DOUBLE PRECISION              :: STARK
      
      HROT = 0.D0
      HSTARK = 0.D0

      KC = 0
      DO ICOL = 1, NHAM
         KC = KC+1
         NC  = JQN(KC,IB)
         MNC = JQN(KC+NMAX+1,IB)
         KR  =  0
         DO IROW = 1, ICOL
            KR  = KR+1
            NR  = JQN(KR,IB)
            MNR = JQN(KR+NMAX+1,IB)
            IF (IROW.EQ.ICOL) HROT(IROW,ICOL) = DBLE((NC+1)*NC)
            IF (MNC.EQ.MNR .AND. ABS(NC-NR).LE.1 .AND. NC.NE.NR) 
     &      HSTARK(IROW,ICOL) = STARK(NC,NR,MNC)
         ENDDO
      ENDDO

      DO ICOL = 1, NHAM
         DO IROW = ICOL+1, NHAM
            HSTARK(IROW,ICOL) = HSTARK(ICOL,IROW)
         ENDDO
      ENDDO

      ! In cm-1 now:
      HROT   = FROT * HROT
      HSTARK = EFLD * FSTARK * HSTARK
      HAM    = HROT + HSTARK

      RETURN
      END SUBROUTINE HAMSTARK
!=========================================================================
      DOUBLE PRECISION FUNCTION STARK(N,NP,MN)

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: N, NP, MN

      DOUBLE PRECISION    :: THRJ, THREEJ

      STARK = SQRT(DBLE((2*N+1)*(2*NP+1)))
     &      * THREEJ(N,1,NP)
     &      * THRJ(  DBLE(N),   1.D0,   DBLE(NP),
     &             -DBLE(MN),   0.D0,   DBLE(MN))

      IF (MOD(MN,2).NE.0) STARK = -STARK

      RETURN
      END FUNCTION STARK
!=========================================================================
      SUBROUTINE DIPOLE_SPFXD(N,EVEC,HAM,DIPIND)

      USE ROTORS_MW_MOD, ONLY: FLDREF, FSTARK, DIPOLE, ZERO
      IMPLICIT NONE

      INTEGER, INTENT(IN)           :: N

      DOUBLE PRECISION, INTENT(IN)  :: EVEC(N,N)

      DOUBLE PRECISION, INTENT(IN)  :: HAM(N,N)

      DOUBLE PRECISION, INTENT(OUT) :: DIPIND(N)

      INTEGER                       :: I

      DOUBLE PRECISION              :: FLD, HH(N,N), DIPMAT(N,N)

      FLD = FLDREF
      IF (ABS(FLDREF) .LT. ZERO) FLD = ZERO

      ! In units of debye of dipole
      HH = DIPOLE*HAM/FLD/FSTARK

      DIPMAT = MATMUL(TRANSPOSE(EVEC),MATMUL(HH,EVEC))

      FORALL(I=1:N) DIPIND(I) = DIPMAT(I,I)

      RETURN
      END SUBROUTINE DIPOLE_SPFXD
!=========================================================================
!****************************** FACTORIAL ******************************
      FUNCTION FACTORIAL(N)
      IMPLICIT NONE
!-----------------------------------------------------------------------
!  DATE (LAST UPDATE): 19/06/12       STATUS: FINISHED, TESTED         |
!  AUTHOR: MAYKEL LEONARDO GONZALEZ MARTINEZ                           |
!-----------------------------------------------------------------------
!             THIS FUNCTION EVALUATES THE FACTORIAL OF N               |
!    (IT USES A VERY SIMPLE ALGORITHM, IT SHOULD BE OK FOR SMALL N)    |
!-----------------------------------------------------------------------
      INTEGER N,FCT10(0:10),I
      INTEGER*8 FACTORIAL
      DATA FCT10/1,1,2,6,24,120,720,5040,40320,362880,3628800/

      SELECT CASE (N)
      CASE (:-1)
         STOP " ERROR IN FACTORIAL: N < 0"
      CASE (0:10)
         FACTORIAL=FCT10(N)
      CASE DEFAULT
         FACTORIAL=FCT10(10)
         DO I=11,N
            FACTORIAL=FACTORIAL*I
         ENDDO
      END SELECT

      RETURN
      END FUNCTION FACTORIAL
!=========================================================================
      SUBROUTINE VMFDR(NDC,MNC,NPHC,NDR,MNR,NPHR,V)

      USE ROTORS_MW_MOD, ONLY: MNIDX, NDIM, NNIDX, JQN, EVECREF, ZERO,
     &                         FOMGA, POLZN,EAC_D_IN_INV_CM
     
      USE PHYSICAL_CONSTANTS, ONLY: HZ_IN_INV_CM

      ! EAC_D_IN_IV_CM is calculated and borrowed by J.D

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NDC, MNC, NPHC, NDR, MNR, NPHR

      DOUBLE PRECISION, INTENT(OUT) :: V

      INTEGER :: IBC, NDIMC, IBR, NDIMR, INC, INR, IC, IR, NC, NR

      DOUBLE PRECISION :: TERMP, TERMM, TTERMP, TTERMM, TERMLOC

      V = 0.D0
      TERMLOC=0.0D0

      ! Pointers to the block for MNC and MNR
      IBC = MNIDX(MNC)
      IBR = MNIDX(MNR)
      ! Sizes of the blocks for MNC and MNR
      NDIMC = NDIM(IBC)
      NDIMR = NDIM(IBR)
      ! Pointers to the indices for NDC and NDR
      INC = NNIDX(NDC,IBC)
      INR = NNIDX(NDR,IBR)


!     WRITE(16,*)                                     
!     WRITE(16,*)                                     
!     WRITE(16,*) "*********************************"  
!     WRITE(16,*)                                     
!     WRITE(16,*) 'Input Rabi frequency in Cm^-1:', FOMGA
!     WRITE(16,*)                                       
!     WRITE(16,*) 'Eac*d in Cm^-1:',  EAC_D_IN_INV_CM  
!     WRITE(16,*) "*********************************"  
!     WRITE(16,*)                                     

!     STOP 123

      TTERMP=0.0D0
      TTERMM=0.0D0

       TERMP=0.0D0
       TERMM=0.0D0


      DO IC = 1, NDIMC
         IF (ABS(EVECREF(IC,INC,IBC)).LT.ZERO) CYCLE
         NC  = JQN(IC,IBC)
         DO IR = 1, NDIMR
            IF (ABS(EVECREF(IR,INR,IBR)).LT.ZERO) CYCLE
            NR = JQN(IR,IBR)
            IF (ABS(NC-NR).NE.1) CYCLE
            IF (NPHR-NPHC.EQ.1) THEN
               CALL MATCOUPLE(NC,MNC,NR,MNR,POLZN,TERMP)
               TTERMP = EVECREF(IC,INC,IBC)*EVECREF(IR,INR,IBR)
     &                * TERMP          
               IF (POLZN.EQ.1) TERMP = -TERMP
               TTERMP = EVECREF(IC,INC,IBC)*EVECREF(IR,INR,IBR)
     &                * TERMP          
            ELSE                !
               TERMP=0.0D0     ! By J.D
               TTERMP = EVECREF(IC,INC,IBC)*EVECREF(IR,INR,IBR)
     &                * TERMP          
            ENDIF
            IF (NPHC-NPHR.EQ.1) THEN
               CALL MATCOUPLE(NC,MNC,NR,MNR,-POLZN,TERMM)
               TTERMM = EVECREF(IC,INC,IBC)*EVECREF(IR,INR,IBR)
     &                * TERMM          
               IF (POLZN.EQ.-1) TERMM = -TERMM
               TTERMM = EVECREF(IC,INC,IBC)*EVECREF(IR,INR,IBR)
     &                * TERMM          
            ELSE                !
               TERMM=0.0D0     ! By J.D
               TTERMM = EVECREF(IC,INC,IBC)*EVECREF(IR,INR,IBR)
     &                * TERMM          
            ENDIF
            TERMLOC=0.5D0*(-EAC_D_IN_INV_CM*(TTERMP+TTERMM))
!           V = V-0.5D0*FOMGA*(TTERMP+TTERMM)
!           V = V+0.5D0*ABS(-EAC_D_IN_INV_CM*(TTERMP+TTERMM)) ! Abs val
            V = V+TERMLOC ! Proper val
!           V = V+(-EAC_D_IN_INV_CM*(TTERMP+TTERMM)) ! Tijs and Hutson; 
!                                                   with +/- values too
!           V = V+0.5D0*FOMGA
!           V = V+FOMGA  ! Only Positive and Hutson-Karman formula
            WRITE(16,*)"==============================================" 
            WRITE(16,*) 
            WRITE(16,*) NC, MNC, NPHC, NR, MNR, NPHR
            WRITE(16,*) 
            WRITE(16,*) FOMGA, FOMGA/(HZ_IN_INV_CM*1.D6),
     &      (-EAC_D_IN_INV_CM)/(HZ_IN_INV_CM*1.D6), TERMP, TERMM, 
     &      (-EAC_D_IN_INV_CM*(TERMP+TERMM))/(HZ_IN_INV_CM*1.D6),
     &      (V/(HZ_IN_INV_CM*1.D6))
            WRITE(16,*) 
         ENDDO
      ENDDO

      RETURN
      END SUBROUTINE VMFDR
!=========================================================================
      SUBROUTINE MATCOUPLE(NC,MNC,NR,MNR,POLZN,V)

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NC, MNC, NR, MNR, POLZN

      DOUBLE PRECISION, INTENT(OUT) :: V

      DOUBLE PRECISION    :: THRJ, THREEJ

      V = THREEJ(NC,1,NR)
     &  * THRJ  (DBLE(NC),1.D0,DBLE(NR),
     &          -DBLE(MNC),DBLE(POLZN),DBLE(MNR)) 

      V = V*SQRT((2.D0*DBLE(NC)+1.D0)*(2.D0*DBLE(NR)+1.D0))

      IF (MOD(MNC,2).NE.0) V = -V

      RETURN
      END SUBROUTINE MATCOUPLE
!=========================================================================
      SUBROUTINE VMMDR(N1C, MN1C, N2C, MN2C, N1R, MN1R, N2R, MN2R,
     &                 LCOLDIF, ANORMC, IETAC, LROWDIF, ANORMR, IETAR,
     &                 LC, LR, LM1, LM2, LMSUM, MTOT, V) 
      ! Multipole-multipole interaction term

      USE ROTORS_MW_MOD, ONLY: MNIDX, NDIM, NNIDX, JQN, NMAX,
     &                         EVECREF, ZERO 

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: N1C, MN1C, N2C, MN2C

      INTEGER, INTENT(IN) :: N1R, MN1R, N2R, MN2R

      INTEGER, INTENT(IN) :: LC, LR, IETAC, IETAR

      INTEGER, INTENT(IN) :: LM1, LM2, LMSUM, MTOT

      LOGICAL, INTENT(IN) :: LCOLDIF, LROWDIF

      DOUBLE PRECISION, INTENT(IN) :: ANORMC, ANORMR

      DOUBLE PRECISION, INTENT(OUT) :: V 

      INTEGER :: IB1C, ND1C, IN1C, IB2C, ND2C, IN2C

      INTEGER :: IB1R, IN1R, IB2R, IN2R

      INTEGER :: I1C, I2C, I1R, I2R

      INTEGER :: N1, MN1, N2, MN2, ML

      INTEGER :: N1P, MN1P, N2P, MN2P, MLP

      INTEGER :: MN1DIF, MN2DIF

      DOUBLE PRECISION :: P1, P2, P3, P4, P5

      DOUBLE PRECISION :: TERM, ELEM 

      ! External functions
      DOUBLE PRECISION :: THREEJ, THRJ

      V = 0.D0

      ! Pointer to the block for this MN1C
      IB1C = MNIDX(MN1C)
      ! Size of the block for this MN1C
      ND1C = NDIM(IB1C)
      ! Pointer to the index for this N1C within block IB1C
      IN1C = NNIDX(N1C,IB1C)

      IB2C = MNIDX(MN2C)
      ND2C = NDIM(IB2C)
      IN2C = NNIDX(N2C,IB2C)

      IB1R = MNIDX(MN1R)
      IN1R = NNIDX(N1R,IB1R)

      IB2R = MNIDX(MN2R)
      IN2R = NNIDX(N2R,IB2R)

      ! Now loops over primitive basis set
      ! Just loop over the relevant block for MN1C
      DO I1C = 1, ND1C
         ! Get N1 from JQN array
         N1  = JQN(I1C,IB1C)
         MN1 = MN1C
         IF (ABS(EVECREF(I1C,IN1C,IB1C)).LT.ZERO) CYCLE
         P1  = SQRT(DBLE(2*N1+1))

      DO I2C = 1, ND2C
         N2  = JQN(I2C,IB2C)
         MN2 = MN2C
         IF (ABS(EVECREF(I2C,IN2C,IB2C)).LT.ZERO) CYCLE
         P2  = P1*SQRT(DBLE(2*N2+1))
         ML  = MTOT-MN1-MN2
         IF (MOD(ML+LM1+MN1+MN2,2).NE.0) P2 = -P2
  
         ! Loop over primed quantum numbers according to 3-j symbol
         ! restrictions
         DO N1P = ABS(N1-LM1), N1+LM1
         
            ! Considering <AC|BC|V|AR|BR> term
  
            IF(N1P .GT. NMAX) CYCLE
            MN1P = MN1R
            IF(N1P .LT. ABS(MN1P)) CYCLE
            ! Pointer to the index for this N1P within block IB1R
            I1R  = NNIDX(N1P,IB1R)
            IF (ABS(EVECREF(I1R,IN1R,IB1R)).LT.ZERO) CYCLE
            ! 3j restriction
            IF (N1+N1P .EQ. 0) CYCLE
            P3 = P2*THREEJ(N1,LM1,N1P)*SQRT(DBLE(2*N1P+1))
            IF (ABS(P3) .LT. ZERO) CYCLE
            MN1DIF = MN1-MN1P
            P4 = P3*THRJ( DBLE(N1),     DBLE(LM1), DBLE(N1P),
     &                  -DBLE(MN1),  DBLE(MN1DIF),DBLE(MN1P))
            IF (ABS(P4) .LT. ZERO) CYCLE
            
         DO N2P = ABS(N2-LM2), N2+LM2
            IF(N2P .GT. NMAX) CYCLE
            MN2P = MN2R
            IF(N2P .LT. ABS(MN2P)) CYCLE
            I2R  = NNIDX(N2P,IB2R)
            IF (ABS(EVECREF(I2R,IN2R,IB2R)).LT.ZERO) CYCLE
            MLP  = MTOT-MN1P-MN2P
            ! 3j restriction
            IF (N2+N2P .EQ. 0) CYCLE
            P5 = P4*THREEJ(N2,LM2,N2P)*SQRT(DBLE(2*N2P+1))
            IF (ABS(P5) .LT. ZERO) CYCLE
            MN2DIF = MN2-MN2P
  
            ELEM = P5*THRJ( DBLE(N2),     DBLE(LM2), DBLE(N2P),
     &               -DBLE(MN2),DBLE(MN2DIF),DBLE(MN2P))
     &      *THRJ( DBLE(LM1),   DBLE(LM2),         DBLE(LMSUM),
     &      DBLE(MN1DIF),DBLE(MN2DIF),-DBLE(MN1DIF+MN2DIF))
     &      *THRJ(DBLE(LC),DBLE(LMSUM),DBLE(LR),
     &           -DBLE(ML),-DBLE(MN1DIF+MN2DIF),DBLE(MLP))
  
            IF (ABS(ELEM) .LT. ZERO) CYCLE
  
            TERM = EVECREF(I1C,IN1C,IB1C)*EVECREF(I2C,IN2C,IB2C)*
     &             EVECREF(I1R,IN1R,IB1R)*EVECREF(I2R,IN2R,IB2R)*
     &             ELEM*ANORMC*ANORMR
  
            V = V+TERM
            
         ENDDO
         ENDDO ! Bra primitive basis set loop
  
         ! Next to consider exchange symmetry if both column and
         ! row pair states are non-identical
         ! If one of column or row pair state is identical then
         ! skip the rest
         IF ((.NOT.LCOLDIF) .OR. (.NOT.LROWDIF)) CYCLE
  
         ! Recalculate the primed primitive pair states by
         ! exchanging them against the original explicit row basis
         ! set
         DO N1P = ABS(N1-LM1), N1+LM1

            ! Considering <AC|BC|V|BR|AR> term
  
            IF(N1P .GT. NMAX) CYCLE
            MN1P = MN2R
            IF(N1P .LT. ABS(MN1P)) CYCLE
            I1R  = NNIDX(N1P,IB2R)
            IF (ABS(EVECREF(I1R,IN2R,IB2R)).LT.ZERO) CYCLE
            ! 3j restriction
            IF (N1+N1P .EQ. 0) CYCLE
            P3 = P2*THREEJ(N1,LM1,N1P)*SQRT(DBLE(2*N1P+1))
            IF (ABS(P3) .LT. ZERO) CYCLE
            MN1DIF = MN1-MN1P
            P4 = P3*THRJ( DBLE(N1),     DBLE(LM1), DBLE(N1P),
     &                  -DBLE(MN1),  DBLE(MN1DIF),DBLE(MN1P))
            IF (ABS(P4) .LT. ZERO) CYCLE
            
         DO N2P = ABS(N2-LM2), N2+LM2
            IF(N2P .GT. NMAX) CYCLE
            MN2P = MN1R
            IF(N2P .LT. ABS(MN2P)) CYCLE
            I2R  = NNIDX(N2P,IB1R)
            IF (ABS(EVECREF(I2R,IN1R,IB1R)).LT.ZERO) CYCLE
            MLP  = MTOT-MN1P-MN2P
            ! 3j restriction
            IF (N2+N2P .EQ. 0) CYCLE
            P5 = P4*THREEJ(N2,LM2,N2P)*SQRT(DBLE(2*N2P+1))
            IF (ABS(P5) .LT. ZERO) CYCLE
            MN2DIF = MN2-MN2P
  
            ELEM = P5*THRJ( DBLE(N2),     DBLE(LM2), DBLE(N2P),
     &                    -DBLE(MN2),  DBLE(MN2DIF),DBLE(MN2P))
     &      *THRJ( DBLE(LM1),   DBLE(LM2),         DBLE(LMSUM),
     &      DBLE(MN1DIF),DBLE(MN2DIF),-DBLE(MN1DIF+MN2DIF))
     &      *THRJ(DBLE(LC),DBLE(LMSUM),DBLE(LR),
     &           -DBLE(ML),-DBLE(MN1DIF+MN2DIF),DBLE(MLP))
  
            IF (ABS(ELEM) .LT. ZERO) CYCLE
  
            TERM = EVECREF(I1C,IN1C,IB1C)*EVECREF(I2C,IN2C,IB2C)*
     &             EVECREF(I1R,IN2R,IB2R)*EVECREF(I2R,IN1R,IB1R)*
     &             ELEM*ANORMC*ANORMR*IETAR
  
            V = V+TERM
            
         ENDDO
         ENDDO ! Bra primitive basis set loop
  
      ENDDO
      ENDDO ! Primitive basis set loop

      ! Due to the symmetry of the like multipole-multipole
      ! potential form, exchanging the column or row or 
      ! both pair states will have same value
      IF (LM1.EQ.LM2 .AND. (LCOLDIF.OR.LROWDIF)) V = 2.D0*V

      ! Next consider exchange symmetry due to unlike
      ! multipole-multipole terms
      ! Skip if this is like multipole-multipole term
      IF (LM1 .EQ. LM2) RETURN 

      !!! CAUTION: This part of the code is not tested !!!

      ! Skip if both column and row pair states are identical
      IF ((.NOT.LCOLDIF) .AND. (.NOT.LROWDIF)) RETURN

      ! Exchange column pair state
      DO I1C = 1, ND2C
         N1  = JQN(I1C,IB2C)
         MN1 = MN2C
         IF (ABS(EVECREF(I1C,IN2C,IB2C)).LT.ZERO) CYCLE
         P1  = SQRT(DBLE(2*N1+1))

      DO I2C = 1, ND1C
         N2  = JQN(I2C,IB1C)
         MN2 = MN1C
         IF (ABS(EVECREF(I2C,IN1C,IB1C)).LT.ZERO) CYCLE
         P2  = P1*SQRT(DBLE(2*N2+1))
         ML  = MTOT-MN1-MN2
         IF (MOD(ML+LM1+MN1+MN2,2).NE.0) P2 = -P2

         ! Exchange row pair state
         DO N1P = ABS(N1-LM1), N1+LM1

            ! Considering <BC|AC|V|BR|AR> term

            IF(N1P .GT. NMAX) CYCLE
            MN1P = MN2R
            IF(N1P .LT. ABS(MN1P)) CYCLE              
            I1R  = NNIDX(N1P,IB2R)
            IF (ABS(EVECREF(I1R,IN2R,IB2R)).LT.ZERO) CYCLE
            ! 3j restriction
            IF (N1+N1P .EQ. 0) CYCLE
            P3 = P2*THREEJ(N1,LM1,N1P)*SQRT(DBLE(2*N1P+1))
            IF (ABS(P3) .LT. ZERO) CYCLE
            MN1DIF = MN1-MN1P
            P4 = P3*THRJ( DBLE(N1),     DBLE(LM1), DBLE(N1P),
     &                  -DBLE(MN1),DBLE(MN1DIF),DBLE(MN1P))
            IF (ABS(P4) .LT. ZERO) CYCLE
            
         DO N2P = ABS(N2-LM2), N2+LM2
            IF(N2P .GT. NMAX) CYCLE
            MN2P = MN1R
            IF(N2P .LT. ABS(MN2P)) CYCLE              
            I2R  = NNIDX(N2P,IB1R)
            IF (ABS(EVECREF(I2R,IN1R,IB1R)).LT.ZERO) CYCLE
            MLP  = MTOT-MN1P-MN2P
            ! 3j restriction
            IF (N2+N2P .EQ. 0) CYCLE
            P5 = P4*THREEJ(N2,LM2,N2P)*SQRT(DBLE(2*N2P+1))
            IF (ABS(P5) .LT. ZERO) CYCLE
            MN2DIF = MN2-MN2P
        
            ELEM = P5*THRJ( DBLE(N2),     DBLE(LM2), DBLE(N2P),
     &                  -DBLE(MN2),DBLE(MN2DIF),DBLE(MN2P))
     &      *THRJ( DBLE(LM1),   DBLE(LM2),         DBLE(LMSUM),
     &      DBLE(MN1DIF),DBLE(MN2DIF),-DBLE(MN1DIF+MN2DIF))
     &      *THRJ(DBLE(LC),DBLE(LMSUM),DBLE(LR),
     &           -DBLE(ML),-DBLE(MN1DIF+MN2DIF),DBLE(MLP))
        
            IF (ABS(ELEM) .LT. ZERO) CYCLE
        
            TERM = EVECREF(I1C,IN2C,IB2C)*EVECREF(I2C,IN1C,IB1C)*
     &             EVECREF(I1R,IN2R,IB2R)*EVECREF(I2R,IN1R,IB1R)*
     &             ELEM*ANORMC*ANORMR*IETAC*IETAR
        
            V = V+TERM
            
         ENDDO
         ENDDO ! Bra primitive basis set loop
        
         IF (.NOT. LROWDIF) CYCLE
        
         ! No exchange for row pair state
         DO N1P = ABS(N1-LM1), N1+LM1

            ! Considering <BC|AC|V|AR|BR> term
        
            IF(N1P .GT. NMAX) CYCLE              
            MN1P = MN1R
            IF(N1P .LT. ABS(MN1P)) CYCLE              
            I1R  = NNIDX(N1P,IB1R)
            IF (ABS(EVECREF(I1R,IN1R,IB1R)).LT.ZERO) CYCLE
            ! 3j restriction
            IF (N1+N1P .EQ. 0) CYCLE
            P3 = P2*THREEJ(N1,LM1,N1P)*SQRT(DBLE(2*N1P+1))
            IF (ABS(P3) .LT. ZERO) CYCLE
            MN1DIF = MN1-MN1P
            P4 = P3*THRJ( DBLE(N1),     DBLE(LM1), DBLE(N1P),
     &                  -DBLE(MN1),DBLE(MN1DIF),DBLE(MN1P))
            IF (ABS(P4) .LT. ZERO) CYCLE
            
         DO N2P = ABS(N2-LM2), N2+LM2
            IF(N2P .GT. NMAX) CYCLE              
            MN2P = MN2R
            IF(N2P .LT. ABS(MN2P)) CYCLE              
            I2R  = NNIDX(N2P,IB2R)
            IF (ABS(EVECREF(I2R,IN2R,IB2R)).LT.ZERO) CYCLE
            MLP  = MTOT-MN1P-MN2P
            ! 3j restriction
            IF (N2+N2P .EQ. 0) CYCLE
            P5 = P4*THREEJ(N2,LM2,N2P)*SQRT(DBLE(2*N2P+1))
            IF (ABS(P5) .LT. ZERO) CYCLE
            MN2DIF = MN2-MN2P
        
            ELEM = P5*THRJ( DBLE(N2),     DBLE(LM2), DBLE(N2P),
     &                  -DBLE(MN2),DBLE(MN2DIF),DBLE(MN2P))
     &      *THRJ( DBLE(LM1),   DBLE(LM2),         DBLE(LMSUM),
     &      DBLE(MN1DIF),DBLE(MN2DIF),-DBLE(MN1DIF+MN2DIF))
     &      *THRJ(DBLE(LC),DBLE(LMSUM),DBLE(LR),
     &           -DBLE(ML),-DBLE(MN1DIF+MN2DIF),DBLE(MLP))
        
            IF (ABS(ELEM) .LT. ZERO) CYCLE
        
            TERM = EVECREF(I1C,IN2C,IB2C)*EVECREF(I2C,IN1C,IB1C)*
     &             EVECREF(I1R,IN1R,IB1R)*EVECREF(I2R,IN2R,IB2R)*
     &             ELEM*ANORMC*ANORMR*IETAC
        
            V = V+TERM
            
         ENDDO
         ENDDO ! Bra primitive basis set loop

      ENDDO
      ENDDO ! Primitive basis set loop

      RETURN
      END SUBROUTINE VMMDR
!=========================================================================
      SUBROUTINE BUBBLE_SORT(N, VEC, IDX)
      ! Sort an array of floats into increasing order using bubble
      ! algorithm. Works efficiently for usually small arrays with 
      ! nearly sorted elements.

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: N

      DOUBLE PRECISION, INTENT(INOUT) :: VEC(N)

      INTEGER, INTENT(INOUT) :: IDX(N)

      INTEGER :: BUBBLE, J, LSUP

      DOUBLE PRECISION :: TEMP

      LSUP = N

      DO WHILE (LSUP .GT. 1)
        BUBBLE = 0 !bubble is the greatest element out of order
        DO J = 1, LSUP-1
          IF (VEC(J) .GT. VEC(J+1)) THEN
            TEMP = VEC(J)
            VEC(J) = VEC(J+1)
            VEC(J+1) = TEMP
            BUBBLE = J
            TEMP = IDX(J)
            IDX(J) = IDX(J+1)
            IDX(J+1) = TEMP
          ENDIF 
        ENDDO
        LSUP = BUBBLE   
      ENDDO
 
      END SUBROUTINE BUBBLE_SORT              
!=========================================================================
