      SUBROUTINE VINIT(I,RRM,EEPSIL)
      USE physical_constants, ONLY: bohr_to_Angstrom, hartree_in_inv_cm
      USE potential, ONLY: LAMBDA
      USE efvs, ONLY: SCALAM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE
C-----------------------------------------------------------------------
C  DATE (LAST UPDATE): 18/10/19                        STATUS: WORKING |
C  AUTHOR: M D Frye                                                    |
C  BASED ON MGLM'S VSTAR FOR LI+YB INTERACTION AS EDITED BY DGG        |
C-----------------------------------------------------------------------
C     THIS SUBROUTINE INCLUDES TO BOUND & MOLSCAT, VIA THE VSTAR       |
C     MECHANISM,                                                       |
C       POTENTIALS FOR 1S + L ATOM. POTENTIALS CAN BE DEFINED FOR      |
C       DIFFERENT |LAMBDA|, THE PROJECTION OF THE ATOMIC L ONTO THE    |
C       INTERATOMIC AXIS. THESE ARE CONVERTED INTO ANISOTRPIC          |
C       COEFFICIENTS (SEE EG EQ. (5) OF PRA 88, 053420 (2013).         |
C       ALTERNATIVELY, THE ANISOTROPIC COEFFICIENTS CAN BE DEFINED     |
C       DIRECTLY.                                                      |
C       CURRENTLY ONLY A LENNARD-JONES FUNCTIONAL FORM IS SUPPORTED    |
C       BUT THIS COULD EASILY BE GENERALISED.                          |
C       DESIGNED TO WORK WITH ROUTINES IN                              |
C       base9-Sat_Las_wf_whf_cpld_MDF.f                                |
C       NB EARLIER VERSIONS OF A SIMILAR BASE9 MAY NOT WORK DUE TO     |
C       CHANGES IN THE MAPPING OF THE POTENTIAL.                       |
C                                                                      |
C       (AT INPUT, X IN ANGSTROMS.  AT OUTPUT, ENERGIES IN CM-1)       |
C-----------------------------------------------------------------------
C USES V3Sg_INIT, V3Pg_INIT,
C      V3Sg_EVAL, V3Pg_EVAL,
C PARAMETERS
      PARAMETER (MXL=100,MXNPOT=10)
C
C  20-09-2016: UPDATED TO USE MODULE (physical_constants) THAT CONTAINS
C              CONSISTENT AND UP-TO-DATE VALUES
      PARAMETER (A02AA=bohr_to_Angstrom)
      LOGICAL FLG_ANISO
C 'EXTRA' VARIABLES:
      DIMENSION LAM(1,MXLAM), NPTS(NDIM), XPT(MXPT,NDIM),
     1          XWT(MXPT,NDIM), XX(MX)
C NAMELISTS
      DIMENSION C6(MXNPOT), DE(MXNPOT), C12(MXNPOT)
      NAMELIST /POTL9/ ITYPE, C6, DE, NPOT, C6_REF, DE_REF, flg_aniso
C
      RM=RRM
      EPSIL=EEPSIL
      IF (I.NE.1) RETURN
      WRITE(*,*) " IN VSTAR FOR 1S+L atom."
      IF(FLG_ANISO)THEN
        WRITE(*,*) " FLG_ANISO=T, SO ANISOTROPY REPRESENTATION."
        WRITE(*,*) " CONTRIBUTIONS FOR EVEN K UP TO 2*L"
        WRITE(*,*) " EACH OF THESE WILL BE REPRESENTED BY A "
     1           "LENNARD-JONES POTENTIAL"
        WRITE(*,*) "  V(R)=C6/R**6+C12/R**12"
        WRITE(*,*) " INPUT POTENTIAL PARAMETERS ARE:"
        DO II=1,NPOT
          WRITE(*,*) "POTENTIAL TERM",II," WITH K=",2*(II-1)
          WRITE(*,*) "C6 = ",C6(II),"(ATOMIC UNITS); DEPTH DE=",
     1             DE(II),"CM-1"
          C6(II)=C6(II)*A02AA**6*HARTREE_IN_INV_CM
          IF (DE(II).EQ.0.0) THEN
            C12(II)=0.0
          ELSE
            C12(II)=C6(II)**2/(4.0*DE(II))
          ENDIF
          C6(II)=C6(II)/(EPSIL*RM**6)
          C12(II)=C12(II)/(EPSIL*RM**12)
          WRITE(*,*) "CONVERTED TO RM AND EPSIL UNITS:"
          WRITE(*,*) "C6 = ",C6(II),"; C12=",C12(II)
          WRITE(*,*)
        ENDDO
      else
        WRITE(*,*) " FLG_ANISO=F, SO LAMBDA POTENTIAL REPRESENTATION."
        WRITE(*,*) " THERE ARE NPOT INTERACTION POTENTIALS, EACH WITH A"
     1           " DIFFERENT MAGNITUDE OF THE PROJECTION OF L_B ONTO "
     2           "THE INTERNUCLEAR AXIS, CALLED LAMBDA."
        WRITE(*,*) " EACH OF THESE POTENTIALS WILL BE REPRESENTED BY A "
     1           "LENNARD-JONES POTENTIAL"
        WRITE(*,*) "  V(R)=C6/R**6+C12/R**12"
        WRITE(*,*) " INPUT POTENTIAL PARAMETERS ARE:"
        DO II=1,NPOT
          WRITE(*,*) "POTENTIAL",II," WITH |LAMBDA|=",II-1
          WRITE(*,*) "C6 = ",C6(II),"(ATOMIC UNITS); DEPTH DE=",
     1             DE(II),"CM-1"
          C6(II)=C6(II)*A02AA**6*HARTREE_IN_INV_CM
          C12(II)=C6(II)**2/(4.0*DE(II))
          C6(II)=C6(II)/(EPSIL*RM**6)
          C12(II)=C12(II)/(EPSIL*RM**12)
          WRITE(*,*) "CONVERTED TO RM AND EPSIL UNITS:"
          WRITE(*,*) "C6 = ",C6(II),"; C12=",C12(II)
          WRITE(*,*)
        ENDDO
      endif
      RETURN
C-----------------------------------------------------------------------
      ENTRY VSTAR(I,X,SUM)
C-----------------------------------------------------------------------
      XA0=X*RM/A02AA ! CONVERTS R/X FROM RM UNITS TO BOHRS
      LMBD=LAMBDA(I)
C FOR SEVERAL SPIN MULTIPLICITIES, LAMBDA IS MAPPED AS:
C  LAMBDA = (I2S1TMP-1)*10 + K/2, WITH
C  (1) I2S1TMP => ~ (I2S1+1)TH SPIN MULTIPLICITY (MAX. 10 DIFFERENT I2S1)
C  (2) K => 0, 2, ..., 2*L2 TERM IN THE POTENTIAL EXPANSION, GIVEN I2S1
CCC NB This is a DIFFERENT mapping from other similar vstar routines
CCC (Yb+Yb* or Li+Yb*) and it MUST be run with the appropriately
CCC modified base9 (currently named base9-Sat_Lat_wf_whf_cpld_MDF.f)
      I2S1TMP=INT(LMBD/10)+1
      K=(LMBD-10*(I2S1TMP-1))*2
      l2=npot-1 ! NB not doubled
      xk=dble(k)
      xl2=dble(l2)
      SUM=0.0
      IF (I2S1TMP.EQ.1) THEN
        IF (FLG_ANISO) THEN
          II=K/2+1
          SUM=C6(II)/X**6+SCALAM*C12(II)/X**12
        ELSE
          PREFACTOR=(2*K+1)/((2*L2+1)*THREEJ(L2,K,L2))
          DO LL=-L2,L2
            XLL=DBLE(LL)
            FAC=THRJ(XL2,XK,XL2,XLL,0.D0,-XLL)*(-1)**(LL)
            LLL=ABS(LL)+1
            POTL=C6(LLL)/X**6+SCALAM*C12(LLL)/X**12
            SUM=SUM+PREFACTOR*FAC*POTL
          ENDDO
        ENDIF
      ELSE
        PRINT*, '* * * ERROR * * * SUBROUTINE VSTAR'
        STOP    ' ONLY ONE MULTIPLICITY INCLUDED IN THIS VERSION'
      ENDIF
      SUM=SUM/SCALAM
      RETURN
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, MXLAM, NPTS, NDIM, XPT, XWT,
     1                  MXPT, IVMIN, IVMAX, L1MAX, L2MAX,
     2                  MXLMB, XX, MX, IXFAC)
C-----------------------------------------------------------------------
      ITYPE=1
      NPOT=0
      DO II=1,MXNPOT
        C6(II)=0.D0
        C12(II)=0.D0
        DE(II)=0.D0
      ENDDO
      C6_REF=0.D0
      DE_REF=0.D0
      FLG_ANISO=.FALSE.
      READ(5,POTL9)
      IF (C6_REF.NE.0.D0) THEN
        DO II=1,MXNPOT
        C6(II)=C6(II)*C6_REF
        ENDDO
      ENDIF
      IF (DE_REF.NE.0.D0) THEN
        DO II=1,MXNPOT
        DE(II)=DE(II)*DE_REF
        ENDDO
      ENDIF

      ITYPP=ITYPE
      END SUBROUTINE VINIT
