C***********************************************************************
C file: FVPTR.FOR
C
C history: 15/2/90: COPIED FROM FVPGB.FOR, NEW SUBROUTINE-NAME
C          15/2/90: unused parameters removed by BH
C          15/2/90: wrong comments removed by BH
C          22/2/90: TRICLINIC GEOMETRY
C
C          18/7/90: ICHARG (valence charge) added
C                   DACOUL replaced ENCOUL
C          18/7/90: ARGSC=NPSPTS/PSGMAX is changed into 
C                   ARGSC=(NPSPTS-1)/PSGMAX, to be constant with the 
C                   grid of pseudopotential
C***********************************************************************
      SUBROUTINE FVPTR(NGX,NGY,NGZ,NPSPTS,RECC,VOLC,PSP,PSGMAX,VPS,
     &                 DVPS,LPCTX,LPCTY,LPCTZ,ICHARG,IVPTYP)
C=======================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE FVPTR "}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE CALCULATES THE PSEUDOPOTENTIAL ON THE GRID OF
C RECIPROCAL LATTICE VECTORS.
C THIS SUBROUTINE MUST BE CALLED WHENEVER RECC(i,j) HAS CHANGED
C
C=======================================================================
      IMPLICIT COMPLEX (C)
      DIMENSION RECC(3,*)
      DIMENSION PSP(*)
      DIMENSION VPS(*)
      DIMENSION DVPS(*)
      DIMENSION LPCTX(*),LPCTY(*),LPCTZ(*)
C=======================================================================
      DATA TPI,HALF,SIXTH / 6.2831853072,0.5,0.166666666667 /
      DATA DACOUL /180.95283/
C      COMMON /PCONST/PI,EV,EPSLON0
C
C     DACOUL=EV/EPSLON0 = 180.95270
C           
C     1EV=1.6021892E-12 erg 
C     EPSLON0(PERMITTIVITY OF VACUUM) = 8.85418782E-12 F m-1  
C                                     = 8.85418782E-2  F (Angstrom)-1  
C
C     DATA ENCOUL /324299.886/
C=======================================================================
      VOLINV=1.0/VOLC
      ZEFFG2=PSP(NPSPTS)*PSGMAX**2
C=======================================================================
C CALCULATE THE SCALING FACTOR ARGSC THAT CONVERTS THE MAGNITUDE OF A
C RECIPROCAL LATTICE VECTOR TO THE CORREPONDING POSITION IN THE
C PSEUDOPOTENTIAL ARRAYS
C=======================================================================
      IF (IVPTYP.EQ.0) THEN
        ARGSC=NPSPTS/PSGMAX
      ELSE
        ARGSC=(NPSPTS-1)/PSGMAX
      END IF
C
C      ENCOUL=ICHARG*DACOUL/(ARGSC**2)
C
      ENCOUL=ICHARG*DACOUL*(ARGSC**2)
      NPSPT2=NPSPTS-2
      NI=1
      DO 9100 N=1,NGZ
      GZX=RECC(3,1)*LPCTZ(N)
      GZY=RECC(3,2)*LPCTZ(N)
      GZZ=RECC(3,3)*LPCTZ(N)
      DO 9101 NN=1,NGY
      GYX=RECC(2,1)*LPCTY(NN)
      GYY=RECC(2,2)*LPCTY(NN)
      GYZ=RECC(2,3)*LPCTY(NN)
      DO 9102 NNN=1,NGX
      GXSQ=(RECC(1,1)*LPCTX(NNN)+GYX+GZX)**2
      GYSQ=(RECC(1,2)*LPCTX(NNN)+GYY+GZY)**2
      GZSQ=(RECC(1,3)*LPCTX(NNN)+GYZ+GZZ)**2
C=======================================================================
C CALCULATE THE MAGNITUDE OF THE RECIPROCAL LATTICE VECTOR (SINCE WE ARE
C USING A LOCAL PSEUDOPOTENTIAL THE VALUE OF THE POTENTIAL DEPENDS ONLY
C ON THE MAGNITUDE OF G)
C=======================================================================
      G=SQRT(GXSQ+GYSQ+GZSQ)
      IF(G.GT.GMAX) GMAX=G
C=======================================================================
C CONVERT THE MAGNITUDE OF THE RECIPROCAL LATTTICE VECTOR TO A POSITION
C IN THE PSEUDOPOTENTIAL ARRAYS AND INTERPOLATE THE PSEUDOPOTENTIAL AND
C ITS DERIVATIVE
C=======================================================================
      ARG=(G*ARGSC)+1.0
      NADDR=INT(ARG)
      REM=ARG-NADDR
C
      IF(NADDR.LT.3) THEN
      VPS(NI) = 0.0
      DVPS(NI)= 0.0
      GO TO 9103
      END IF
C
C     Vps = Zeff/G**2,  dVps/dG = -2 Zeff/G**3 = -2*Vps/G
C
      IF(NADDR.GE.NPSPT2) THEN
      VPS(NI)=  ZEFFG2/G**2
      DVPS(NI)=-2.0*VPS(NI)/G
      GO TO 9103
      END IF
C
      V1=PSP(NADDR-1)+(ENCOUL/((NADDR-2)**2))
      V2=PSP(NADDR)+(ENCOUL/((NADDR-1)**2))
      V3=PSP(NADDR+1)+(ENCOUL/((NADDR)**2))
      V4=PSP(NADDR+2)+(ENCOUL/((NADDR+1)**2))
      T0=V2
      T1=((6.0*V3)-(2.0*V1)-(3.0*V2)-V4)*SIXTH
      T2=(V1+V3-(2.0*V2))*HALF
      T3=(V4-V1+(3.0*(V2-V3)))*SIXTH
      VPS(NI)=((-ENCOUL/((ARG-1.0)**2))+T0+
     &         REM*(T1+REM*(T2+REM*T3)))*VOLINV
      DVPS(NI)=((2.0*ENCOUL/(((ARG-1.0)**2)*G))+
     &         (ARGSC*(T1+REM*((2.0*T2)+(3.0*REM*T3)))))*VOLINV
      IF(VPS(NI).GT.10.0) WRITE(*,*) VPS(NI),NI,NADDR
 9103 CONTINUE
      NI=NI+1
 9102 CONTINUE
 9101 CONTINUE
 9100 CONTINUE
      IF(GMAX.GT.PSGMAX) WRITE(*,9999) GMAX,PSGMAX
 9999 FORMAT(/,' WARNING IN FVP_TR:',
     &/,' THE LARGEST G VALUE OF THE FFT GRID (GMAX)=',F12.4,
     &/,' EXCEEDS THE UPLIMIT OF THE POTENTIAL DATASET (PSGMAX)=',F12.4,
     &/,' WE USE EXTRAPOLATION ASSUMING COULOMBLISH BEHAVIOR (Z/G**2):',
     &/,'   V(G)    =  V(PSGMAX)*(PSGMAX**2/G**2)   FOR  G > PSGMAX',
     &/,'  dV(G)/dG = -2 V(G)/G',/)
      RETURN
      END
