      SUBROUTINE VINIT(I,RM,EPSIL)
      USE physical_constants, ONLY: bohr_to_Angstrom
      USE COMMONS, ONLY: RSOFLG, A1, A2
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: I
      DOUBLE PRECISION, INTENT(INOUT) :: RM, EPSIL
      !-----------------------------------------------------------------------
      !  DATE (LAST UPDATE): 21/10/21       STATUS: FINISHED                 |
      !  AUTHOR: BIJIT MUKHERJEE                                             |
      !  BASED ON MGLM'S VSTAR FOR LI+YB INTERACTION AS EDITED BY DGG        |
      !-----------------------------------------------------------------------
      !      THIS SUBROUTINE INCLUDES TO BOUND & MOLSCAT, VIA THE VSTAR      |
      !      MECHANISM, 4 MODEL 'DIABATIC' POTENTIAL CURVES (3{2Sigma+},     |
      !       2{2Pi}, 1{4Pi}, 1{4Sigma+}), AND MODEL R-DEPENDENT SPIN        |
      !      ORBIT COUPLING FUNCTION FOR Rb(2S)+Yb(3P) ARE CALCULATED        |
      !       (AT INPUT, X IN ANGSTROMS. AT OUTPUT, ENERGIES IN CM-1)        |
      !    DESIGNED TO WORK WITH ROUTINES IN base9-Sat_Lat_wf_whf_cpld_BM.f  |
      !-----------------------------------------------------------------------
      INTEGER, PARAMETER :: MXL = 100, IXMX = 200, IEXMX = 200,
     &                      MXDIM = 3
      INTEGER :: CFLAG, ICNSYM, ICNSY2, IHOMO, IHOMO2, IVMIN, IVMAX,
     &           LAMBDA(MXL), LMAX, L1MAX, L2MAX, MMAX, MXLAM, MXSYM,
     &           NHAM, NPOWER(IXMX), NPTS(MXDIM), NPS, NPT, NTERM(MXL)
      DOUBLE PRECISION :: A(IXMX), E(IEXMX)
      LOGICAL :: LVRTP
      CHARACTER(10) :: RMNAME, EPNAME
      INTEGER :: LMBD
      ! VARIABLES FOR ENTRY POINT VSTAR
      DOUBLE PRECISION :: X, SUM, FSO
      DOUBLE PRECISION :: XA0
      INTEGER :: I2S1, K
      DOUBLE PRECISION :: V14SP, V14P, V32SP, V22P

      SAVE
! NAMELISTS (NOTE THAT THIS LIST SHOULD MATCH THAT IN potenl.f)
      NAMELIST /POTL/ A,      CFLAG,  E,     EPNAME, EPSIL,
     &               ICNSYM, ICNSY2, IHOMO, IHOMO2, IVMIN,
     &               IVMAX,  LAMBDA, LMAX,   L1MAX, L2MAX,
     &               LVRTP,  MMAX,   MXLAM,  MXSYM, NHAM,
     &               NPOWER, NPS,    NPTS,   NPT,   NTERM,
     &               RM,     RMNAME

      IF (I.NE.1) RETURN
      REWIND(5)

      !DEFAULT VALUES:
      A = 0.D0
      E = 0.D0
      LAMBDA = 0
      NPOWER = 0
      NPTS = 0
      NTERM = 0

      READ(5,POTL)

      CALL V_POT_INIT

      !IF R-DEP SO CALCULATIONS NEED TO BE DONE
      IF (RSOFLG) THEN
         CALL VSO_INIT
      ENDIF

      RETURN
! C-----------------------------------------------------------------------
      ENTRY VSTAR(I,X,SUM)
! C-----------------------------------------------------------------------
      LMBD = LAMBDA(I)
      IF (LMBD .EQ. -4 .OR. LMBD .EQ. -2) THEN
         IF (.NOT. RSOFLG) THEN
            SUM = 0.D0
            RETURN
         ELSE
            CALL VSO_EVAL(X,FSO)
            IF (LMBD .EQ. -4) THEN
                  SUM = A1 * FSO
                  RETURN
            ELSE
                  SUM = -A2 * FSO
                  RETURN
            ENDIF
         ENDIF
      ELSEIF (LMBD .EQ. -1) THEN
         XA0 = X / bohr_to_Angstrom ! CONVERTS R/X FROM ANGSTROMS TO BOHRS
! C CP1 ELECTRON SPIN-CP2 ELECTRON SPIN R-DEPENDENCE
         SUM = XA0 ** (-3) ! X HAS BEEN CONVERTED FROM ANGSTROMS TO BOHRS
         RETURN
      ELSEIF (LMBD .GE. 0) THEN
! C FOR SEVERAL SPIN MULTIPLICITIES, LAMBDA IS MAPPED AS:
! C  LAMBDA = (I2S1-1)*10 + K, WITH
! C  (1) I2S1 => ~ (I2S1+1)TH SPIN MULTIPLICITY (MAX. 10 DIFFERENT I2S1)
! C  (2) K => 0, 2, ..., 2*L2 TERM IN THE POTENTIAL EXPANSION, GIVEN I2S1
         I2S1 = INT(LMBD/10) + 1
         K = LMBD - 10 * (I2S1 - 1)
         IF (I2S1 .EQ. 1) THEN             ! QUARTET STATES
            CALL V14SP_EVAL(X,V14SP)
            CALL V14P_EVAL(X,V14P)
            IF (K .EQ. 0) THEN
               SUM = (V14SP + 2.D0 * V14P) / 3.D0
            ELSEIF (K .EQ. 2) THEN
               SUM = 5.D0 * (V14SP - V14P) / 3.D0
            ELSE
               WRITE(6,'(/"ONE OF THE NON-NEGATIVE LAMBDA VALUES IS",
     &                    " INCORRECT")')
               STOP
            ENDIF
         ELSEIF (I2S1 .EQ. 2) THEN         ! DOUBLET STATES
            CALL V32SP_EVAL(X,V32SP)
            CALL V22P_EVAL(X,V22P)
            IF (K .EQ. 0) THEN
               SUM = (V32SP + 2.D0 * V22P) / 3.D0
            ELSEIF (K .EQ. 2) THEN
               SUM = 5.D0 * (V32SP - V22P) / 3.D0
            ELSE
               WRITE(6,'(/"ONE OF THE NON-NEGATIVE LAMBDA VALUES IS",
     &                    " INCORRECT")')
               STOP
            ENDIF
         ELSE
            PRINT*, '* * * ERROR * * * SUBROUTINE VSTAR'
            STOP    ' ONLY 2S + 1 = 2, 4 INCLUDED IN THIS VERSION'
         ENDIF
         RETURN
      ENDIF
! C DUMMY ROUTINES
      ENTRY VSTAR1(I,X,SUM)
      ENTRY VSTAR2(I,X,SUM)
      PRINT*, '* * * ERROR * * * SUBROUTINE VSTAR'
      STOP    ' VSTAR: DERIVATIVES NOT IMPLEMENTED'
! C
      RETURN
! ! C-----------------------------------------------------------------------
!       ENTRY POTIN9(ITYPP, LAM, MXLAM9, NPTS9, NDIM, XPT, XWT,
!      &             MXPT, IVMIN9, IVMAX9, L1MAX9, L2MAX9,
!      &             MXLMB, XX, MX, IXFAC)
! ! C-----------------------------------------------------------------------
!       READ(5,POTL9)
      END SUBROUTINE VINIT
! C****************************** V_POT_INIT *******************************
      SUBROUTINE V_POT_INIT
      USE physical_constants, ONLY: bohr_to_Angstrom, hartree_in_inv_cm
      IMPLICIT NONE
      SAVE
! C-----------------------------------------------------------------------
! C  DATE (LAST UPDATE): 21/10/21       STATUS: FINISHED                 |
! C  AUTHOR: BIJIT MUKHERJEE                                             |
! C-----------------------------------------------------------------------
! C    THIS SUBROUTINE EVALUATES THE 4 MODEL R-DEPENDENT SPIN ORBIT      |
! C                        COUPLING FUNCTIONS:                           |
! C-----------------------------------------------------------------------
      CHARACTER*100 :: FILE_IN
      DOUBLE PRECISION :: DE_2S, DE_2P, DE_4P, DE_4S,
     &                    RE_2S, RE_2P, RE_4P, RE_4S,
     &                    BETA_2S, BETA_2P, BETA_4P, BETA_4S,
     &                    AP_2S, AP_2P, AP_4P, AP_4S,
     &                    AQ_2S, AQ_2P, AQ_4P, AQ_4S,
     &                    ALPHA_2S, ALPHA_2P, ALPHA_4P, ALPHA_4S,
     &                    C6_SIGMA, C6_PI, C8_SIGMA, C8_PI
!     VARIABLES FOR ENTRY V_EVAL
      DOUBLE PRECISION :: X, V32SP, V22P, V14P, V14SP,
     &                    P6(6), DP6(6), P8(8), DP8(8)
      DOUBLE PRECISION :: HULBERT_HIRSCHFELDER
      
      NAMELIST /V_PARAM/ DE_2S, DE_2P, DE_4P, DE_4S,
     &                   RE_2S, RE_2P, RE_4P, RE_4S,
     &                   BETA_2S, BETA_2P, BETA_4P, BETA_4S,
     &                   AP_2S, AP_2P, AP_4P, AP_4S,
     &                   AQ_2S, AQ_2P, AQ_4P, AQ_4S,
     &                   ALPHA_2S, ALPHA_2P, ALPHA_4P, ALPHA_4S,
     &                   C6_SIGMA, C6_PI, C8_SIGMA, C8_PI

      FILE_IN = "PEC_param_RbYb.ini"
      OPEN(9, FILE = TRIM(FILE_IN), STATUS = 'OLD', ERR = 2010)

      WRITE(6,'(/"  PARAMETERS FOR ELECTROSATIC POTENTIALS",
     &           " READ FROM FILE: ",A20)') FILE_IN

      READ(9,V_PARAM)
      CLOSE(9)

      C6_SIGMA = C6_SIGMA * hartree_in_inv_cm * bohr_to_Angstrom**6
      C8_SIGMA = 80.D0 * C6_SIGMA * bohr_to_Angstrom**2
      C6_PI = C6_PI * hartree_in_inv_cm * bohr_to_Angstrom**6
      C8_PI = 80.D0 * C6_PI * bohr_to_Angstrom**2
! C
      RETURN
! C------------------------------ V32SP_EVAL -----------------------------
      ENTRY V32SP_EVAL(X,V32SP)
! C-----------------------------------------------------------------------

      V32SP = HULBERT_HIRSCHFELDER(X, DE_2S, RE_2S, BETA_2S, AP_2S,
     &                             AQ_2S)
      CALL DAMP(1, 6, ALPHA_2S, X, P6, DP6)
      CALL DAMP(1, 8, ALPHA_2S, X, P8, DP8)
      V32SP = V32SP - P6(6) * C6_SIGMA / X**6
     &              - P8(8) * C8_SIGMA / X**8
! C
      RETURN
! C------------------------------ V22P_EVAL -----------------------------
      ENTRY V22P_EVAL(X,V22P)
! C-----------------------------------------------------------------------

      V22P = HULBERT_HIRSCHFELDER(X, DE_2P, RE_2P, BETA_2P, AP_2P,
     &                            AQ_2P)
      CALL DAMP(1, 6, ALPHA_2P, X, P6, DP6)
      CALL DAMP(1, 8, ALPHA_2P, X, P8, DP8)
      V22P = V22P - P6(6) * C6_PI / X**6
     &            - P8(8) * C8_PI / X**8
! C
      RETURN
! C------------------------------ V14P_EVAL -----------------------------
      ENTRY V14P_EVAL(X,V14P)
! C-----------------------------------------------------------------------

      V14P = HULBERT_HIRSCHFELDER(X, DE_4P, RE_4P, BETA_4P, AP_4P,
     &                            AQ_4P)
      CALL DAMP(1, 6, ALPHA_4P, X, P6, DP6)
      CALL DAMP(1, 8, ALPHA_4P, X, P8, DP8)
      V14P = V14P - P6(6) * C6_PI / X**6
     &            - P8(8) * C8_PI / X**8
! C
      RETURN
! C------------------------------ V14SP_EVAL -----------------------------
      ENTRY V14SP_EVAL(X,V14SP)
! C-----------------------------------------------------------------------

      V14SP = HULBERT_HIRSCHFELDER(X, DE_4S, RE_4S, BETA_4S, AP_4S,
     &                            AQ_4S)
      CALL DAMP(1, 6, ALPHA_4S, X, P6, DP6)
      CALL DAMP(1, 8, ALPHA_4S, X, P8, DP8)
      V14SP = V14SP - P6(6) * C6_SIGMA / X**6
     &              - P8(8) * C8_SIGMA / X**8
! C
      RETURN
 2010 PRINT*, " ERROR IN V_POT_INIT: COULD NOT OPEN FILE ",FILE_IN
      STOP
      END SUBROUTINE V_POT_INIT
! C****************************** VSO_INIT **********************************
      SUBROUTINE VSO_INIT
      ! C-----------------------------------------------------------------------
      ! C  DATE (LAST UPDATE): 21/10/21       STATUS: FINISHED                 |
      ! C  AUTHOR: BIJIT MUKHERJEE                                             |
      ! C-----------------------------------------------------------------------
      ! C-----------------------------------------------------------------------
      ! C  THIS SUBROUTINE EVALUATES THE R-DEPENDENT SPIN-ORBIT COEFFICIENTS   |
      ! C----------------------------------------------------------------------|
      USE COMMONS, ONLY: A1
      IMPLICIT NONE
      SAVE
      DOUBLE PRECISION :: EPS, SIG, R0
!     VARIABLES FOR ENTRY VSO_EVAL
      DOUBLE PRECISION :: X, FSO
      CHARACTER*100 :: FILE_IN
      
      NAMELIST /SO_PARAM/ A1, EPS, SIG, R0

      FILE_IN = "SO_param_RbYb.ini"

      OPEN(9, FILE = TRIM(FILE_IN), STATUS = 'OLD', ERR = 2014)

      WRITE(6,'(/"  PARAMETERS FOR R-DEPENDENT SPIN-ORBIT",
     &           " COEFFICIENTS READ FROM FILE: ",A20)') FILE_IN

      READ(9,SO_PARAM)
      CLOSE(9)

      RETURN
! C------------------------------ VSO_EVAL -----------------------------
      ENTRY VSO_EVAL(X,FSO)
! C-------------------------------------------------------------------------
      FSO = EPS * (1.D0 - TANH(SIG * (X - R0)))
      RETURN
 2014 PRINT*, " ERROR IN VSO_INIT: COULD NOT OPEN FILE ",FILE_IN
      END SUBROUTINE VSO_INIT

! C************************** HULBERT_HIRSCHFELDER ***********************
      DOUBLE PRECISION FUNCTION HULBERT_HIRSCHFELDER(R,DE,RE,BETA,P,Q)
! C-----------------------------------------------------------------------
! C  DATE (LAST UPDATE): 21/10/21        STATUS: FINISHED                |
! C  AUTHOR: BIJIT MUKHERJEE                                             |
! C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(IN) :: R, DE, RE, BETA, P, Q
      DOUBLE PRECISION :: XX, V

      XX = BETA * (R - RE)
      V = DE * (EXP(-2.D0 * XX) - 2.D0 * EXP(-XX) + P * XX**3
     &  * EXP(-2.D0 * XX) * (1.D0 + Q * XX))
      HULBERT_HIRSCHFELDER = V

      RETURN
      END FUNCTION HULBERT_HIRSCHFELDER
