C----------------------------------------------------------------------C   
C File : serlnl.f                                                      C   
C Created: 08/05/91 IS                                                 C   
C History: See below                                                   C   
C----------------------------------------------------------------------C   
C                                                                       
C   02-JUN-91  PERIODIC BOUNDARY CONDITIONS APPLIED TO IONIC POSITIONS  
C              AND TO RXDIS,RYDIS,RZDIS.   (IS)                         
C
C   10-DEC-91  PERIODIC BOUNDARY CONDITIONS ELIMINATED FROM THE SEARCH
C              OVER REAL-SPACE GRID POINTS. THE OLD VERSION WAS
C              INCORRECT IF RLCORE WAS LARGER THAN HALF THE BOX SIDE.
C              (IS)
C---------------------------------------------------------------------- 
      SUBROUTINE SERLNL(NIONST,MXRLNL,MXRLSH,NIONS,NRLPTS,NRGRPT,NSPEC,
     &       DIRC,RECC,RLCORE,RMAX,NRLNL,PRLSCA,POSION,VRLNL,
     &       NIONSP,IRLNL,NADGRD,NRLPPI,DVRLGR,VRLGRD,NGX,NGY,NGZ)
C
      IMPLICIT REAL*8 (A,B,D-H,O-Z)  
      IMPLICIT COMPLEX*16 (C)        
      LOGICAL WRDIM
C                            
C
C========================================================================
C THIS SUBROUTINE SETS UP THE POTENTIAL AND DERIVATIVES ON THE REAL SPACE
C GRID AROUND EACH ATOM IN PREPARATION FOR CALCULATING H*PSI AND THE
C FORCES ON THE ATOMS - IN FUTURE TO ALLOW MAXIMUM FELXIBILTY THESE TWO
C PROCESSES COULD BE PERFORMED IN SEPARATE SUBROUTINES TO ALLOW THE CHOICE
C OF REAL SPACE PROJECTIONS FOR H*PSI AND RECIPROCAL SPACE FOR THE FORCES
C THE PHASE FACTORS ARE CALCULATED IN A SEPARATE BUT SIMILAR SUBROUTINE
C SINCE ONLY THE PHASE FACTORS CHANGE BETWEEN K-POINTS
C========================================================================
      DIMENSION DIRC(3,3),RECC(3,3)
      DIMENSION RLCORE(NSPEC),RMAX(NSPEC),NRLNL(NSPEC)
      DIMENSION PRLSCA(MXRLNL,NSPEC)
      DIMENSION POSION(3,NIONS,NSPEC),VRLNL(NRLPTS,MXRLNL,NSPEC)
      DIMENSION IRLNL(MXRLNL,NSPEC),NADGRD(NRGRPT,NIONST),NRLPPI(NIONST)
      DIMENSION DVRLGR(NRGRPT,3,MXRLSH,NIONST)
      DIMENSION NIONSP(NSPEC),VRLGRD(NRGRPT,MXRLSH,NIONST)
      DIMENSION POSTMP(3)
C======================================================================== 
C                                                                         
      DATA TPI,HALF,SIXTH / 6.283185307179,0.5,0.166666666667 /
      DATA RTH,RFIV,RFIF / 1.732050807568,2.236067977499,3.872983346207/
      DATA WRDIM /.FALSE./
C
C========================================================================
C INITIALISE THE INDEX FOR THE IONS, NPION, AND CONSTANTS, AND
C START THE LOOP OVER SPECIES OF IONS
C========================================================================
      RLENV1=DSQRT(DIRC(1,1)**2+DIRC(1,2)**2+DIRC(1,3)**2)
      RLENV2=DSQRT(DIRC(2,1)**2+DIRC(2,2)**2+DIRC(2,3)**2)
      RLENV3=DSQRT(DIRC(3,1)**2+DIRC(3,2)**2+DIRC(3,3)**2)
      RLENG1=DSQRT(RECC(1,1)**2+RECC(1,2)**2+RECC(1,3)**2)
      RLENG2=DSQRT(RECC(2,1)**2+RECC(2,2)**2+RECC(2,3)**2)
      RLENG3=DSQRT(RECC(3,1)**2+RECC(3,2)**2+RECC(3,3)**2)
      NPION=1
      DO 1000 NSP=1,NSPEC
C========================================================================
C START BY DETERMINING THE SIZE OF THE REAL SPACE GRID TO BE SEARCHED
C FOR THIS SPECIES OF ION THE ENSURE THAT ALL THE REAL SPACE POINTS
C WITHIN RLCORE(NSP) ARE SAMPLED - THE TRICK TO THIS PROCESS IS TO GO
C ALONG EACH LATTICE VECTOR UNTIL IT HITS THE PLANE DEFINED BY THE OTHER
C LATTICE VECTORS THAT IS TANGENTIAL TO A SPHERE OF RADIUS RLCORE AND
C THEN ADD HALF A LATTICE VECTOR TO ALLOW FOR THE CENTRE OF THE ION
C BEING UP TO HALF A LATTICE VECTOR FROM THE NEAREST GRID POINT
C
C TPI/RLENV1*RLENG1 IS THE ANGLE BETWEEN THE FIRST LATTICE VECTOR AND THE
C NORMAL TO THE PLANE DEFINED BY THE OTHER LATTICE VECTORS
C USE THIS TO DETERMINE SIZE OF GRID TO BE SEARCHED
C========================================================================
      NUMXPT=INT(5.D-1+RLCORE(NSP)*NGX*RLENV1*RLENG1/(TPI*RLENV1))
      NUMYPT=INT(5.D-1+RLCORE(NSP)*NGY*RLENV2*RLENG2/(TPI*RLENV2))
      NUMZPT=INT(5.D-1+RLCORE(NSP)*NGZ*RLENV3*RLENG3/(TPI*RLENV3))
C========================================================================
C CALCULATE THE SCALING FACTOR FOR INTERPOLATING POTENTIALS
C========================================================================
      ARGSC=DFLOAT(NRLPTS-1)/RMAX(NSP)
C========================================================================
C START LOOP OVER THE IONS FOR THIS SPECIES
C========================================================================
      DO 1100 NI=1,NIONSP(NSP)
C========================================================================
      DO 1101 IM=1,3
      POSTMP(IM)=MOD(POSION(IM,NI,NSP),1.D0)  
      POSTMP(IM)=POSTMP(IM)-INT(2.D0*POSTMP(IM))
 1101 CONTINUE
C========================================================================
C FIRST IDENTIFY THE NEAREST GRID POINT TO THE CURRENT ION
C========================================================================
      NGRPX=1+MOD(INT(5.D-1+NGX*POSTMP(1)),NGX)
      NGRPY=1+MOD(INT(5.D-1+NGY*POSTMP(2)),NGY)
      NGRPZ=1+MOD(INT(5.D-1+NGZ*POSTMP(3)),NGZ)
C========================================================================
C START SEARCH OVER GRID TO FIND POINTS WITHIN RLCORE OF THE ION
C INDAD INDICATES THE NUMBER OF GRID POINTS THAT HAVE BEEN FOUND
C WITHIN RLCORE AND IS USED AS THE ADDRESS INDEX FOR PUTTING THE
C POTENTIALS ETC INTO THE DATA ARRAYS
C========================================================================
      INDAD=1
      DO 1200 NZ=-NUMZPT,NUMZPT
      RZDIS=((1.D0*(NGRPZ+NZ-1))/(1.D0*NGZ)) - POSTMP(3)
      RZDIS=MOD(RZDIS,1.D0)
      RZDIS=RZDIS-INT(2.D0*RZDIS)
      DO 1201 NY=-NUMYPT,NUMYPT
      RYDIS=((1.D0*(NGRPY+NY-1))/(1.D0*NGY)) - POSTMP(2)
      RYDIS=MOD(RYDIS,1.D0)
      RYDIS=RYDIS-INT(2.D0*RYDIS)
      DO 1202 NX=-NUMXPT,NUMXPT
      RXDIS=((1.D0*(NGRPX+NX-1))/(1.D0*NGX)) - POSTMP(1)
      RXDIS=MOD(RXDIS,1.D0)
      RXDIS=RXDIS-INT(2.D0*RXDIS)
      XCOR=DIRC(1,1)*RXDIS+DIRC(2,1)*RYDIS+DIRC(3,1)*RZDIS
      YCOR=DIRC(1,2)*RXDIS+DIRC(2,2)*RYDIS+DIRC(3,2)*RZDIS
      ZCOR=DIRC(1,3)*RXDIS+DIRC(2,3)*RYDIS+DIRC(3,3)*RZDIS
      RAD=DSQRT(XCOR**2+YCOR**2+ZCOR**2)
C========================================================================
C CHECK TO SEE WHETHER GRID POINT IS WITHIN RANGE OF NON-LOCAL POTENTIAL
C========================================================================
      IF(RAD.GT.RLCORE(NSP)) GO TO 1202
C========================================================================
C FIND AND STORE ADDRESS OF THIS POINT IN THE REAL SPACE WAVEFUNCTION
C GRID
C========================================================================
      INDDX=1+MOD(NGRPX+NGX+NX-1,NGX)    
      INDDY=1+MOD(NGRPY+NGY+NY-1,NGY)
      INDDZ=1+MOD(NGRPZ+NGZ+NZ-1,NGZ)
      NADGRD(INDAD,NPION)=(INDDZ-1)*NGY*NGX+(INDDY-1)*NGX+INDDX
C========================================================================
C START LOOP OVER NON-LOCAL POTENTIAL COMPONENTS
C========================================================================
      NIND=1
      DO 1210 NNL=1,NRLNL(NSP)
C========================================================================
C INTERPOLATE POTENTIALS AND CALCULATE SPHERICAL HARMONICS ETC
C FOR THIS REAL SPACE DATA POINT
C========================================================================
C INTERPOLATE POTENTIAL AND DERIVATIVE
C========================================================================
      ADDR=1 + RAD*ARGSC
      NADDR=INT(ADDR)
      REM=ADDR-NADDR
      IF(NADDR.EQ.1) THEN
      V1=VRLNL(1,NNL,NSP)
      V2=VRLNL(2,NNL,NSP)
      V3=VRLNL(3,NNL,NSP)
      T0=V1
      T1=(4.D0*V2-3.D0*V1-V3)*HALF
      T2=(V3+V1-2.D0*V2)*HALF
      VRLPR=T0+REM*T1+REM*REM*T2
      DVRLPR=(T1+2.D0*REM*T2)*ARGSC
      ELSE
      V1=VRLNL(NADDR-1,NNL,NSP)
      V2=VRLNL(NADDR,NNL,NSP)
      V3=VRLNL(NADDR+1,NNL,NSP)
      V4=VRLNL(NADDR+2,NNL,NSP)
      T0=V2
      T1=((6.D0*V3)-(2.D0*V1)-(3.D0*V2)-V4)*SIXTH
      T2=(V1+V3-(2.D0*V2))*HALF
      T3=(V4-V1+(3.D0*(V2-V3)))*SIXTH
      VRLPR=T0+REM*(T1+REM*(T2+REM*T3))
      DVRLPR=ARGSC*(T1+REM*((2.D0*T2)+(3.D0*REM*T3)))
      ENDIF
C========================================================================
C 'S' POTENTIAL
C========================================================================
      IF(IRLNL(NNL,NSP).EQ.0) THEN
      VRLGRD(INDAD,NIND,NPION)=VRLPR
      IF(RAD.LE.1.D-6) THEN
      DVRLGR(INDAD,1,NIND,NPION)=0.D0
      DVRLGR(INDAD,2,NIND,NPION)=0.D0
      DVRLGR(INDAD,3,NIND,NPION)=0.D0
      NIND=NIND+1
      GO TO 1210
      ENDIF
      DVRLGR(INDAD,1,NIND,NPION)=-DVRLPR*XCOR/RAD
      DVRLGR(INDAD,2,NIND,NPION)=-DVRLPR*YCOR/RAD
      DVRLGR(INDAD,3,NIND,NPION)=-DVRLPR*ZCOR/RAD
      NIND=NIND+1
      GO TO 1210
      ENDIF
C========================================================================
C 'P' POTENTIAL
C========================================================================
      IF(IRLNL(NNL,NSP).EQ.1) THEN
      IF(RAD.LE.1.D-6) THEN
      VRLGRD(INDAD,NIND,NPION)=0.D0
      VRLGRD(INDAD,NIND+1,NPION)=0.D0
      VRLGRD(INDAD,NIND+2,NPION)=0.D0
      DVRLGR(INDAD,1,NIND,NPION)=-RTH*DVRLPR
      DVRLGR(INDAD,2,NIND,NPION)=0.D0
      DVRLGR(INDAD,3,NIND,NPION)=0.D0
      DVRLGR(INDAD,1,NIND+1,NPION)=0.D0
      DVRLGR(INDAD,2,NIND+1,NPION)=-RTH*DVRLPR
      DVRLGR(INDAD,3,NIND+1,NPION)=0.D0
      DVRLGR(INDAD,1,NIND+2,NPION)=0.D0
      DVRLGR(INDAD,2,NIND+2,NPION)=0.D0
      DVRLGR(INDAD,3,NIND+2,NPION)=-RTH*DVRLPR
      NIND=NIND+3
      GO TO 1210
      ENDIF
      VRLGRD(INDAD,NIND,NPION)=VRLPR*RTH*XCOR/RAD
      VRLGRD(INDAD,NIND+1,NPION)=VRLPR*RTH*YCOR/RAD
      VRLGRD(INDAD,NIND+2,NPION)=VRLPR*RTH*ZCOR/RAD
      DVRLGR(INDAD,1,NIND,NPION)=-RTH*(DVRLPR*(XCOR/RAD)**2 +
     &          (1.D0 - (XCOR/RAD)**2)*VRLPR/RAD)
      DVRLGR(INDAD,2,NIND,NPION)=-(RTH*XCOR*YCOR/RAD**2) *
     &        (DVRLPR - VRLPR/RAD)
      DVRLGR(INDAD,3,NIND,NPION)=-(RTH*XCOR*ZCOR/RAD**2) *
     &        (DVRLPR - VRLPR/RAD)
      DVRLGR(INDAD,1,NIND+1,NPION)=-(RTH*YCOR*XCOR/RAD**2) *
     &        (DVRLPR - VRLPR/RAD)
      DVRLGR(INDAD,2,NIND+1,NPION)=-RTH*(DVRLPR*(YCOR/RAD)**2 +
     &          (1.D0 - (YCOR/RAD)**2)*VRLPR/RAD)
      DVRLGR(INDAD,3,NIND+1,NPION)=-(RTH*YCOR*ZCOR/RAD**2) *
     &        (DVRLPR - VRLPR/RAD)
      DVRLGR(INDAD,1,NIND+2,NPION)=-(RTH*ZCOR*XCOR/RAD**2) *
     &        (DVRLPR - VRLPR/RAD)
      DVRLGR(INDAD,2,NIND+2,NPION)=-(RTH*ZCOR*YCOR/RAD**2) *
     &        (DVRLPR - VRLPR/RAD)
      DVRLGR(INDAD,3,NIND+2,NPION)=-RTH*(DVRLPR*(ZCOR/RAD)**2 +
     &          (1.D0 - (ZCOR/RAD)**2)*VRLPR/RAD)
      NIND=NIND+3
      GO TO 1210
      ENDIF
C========================================================================
C 'D' POTENTIAL
C========================================================================
      IF(IRLNL(NNL,NSP).EQ.2) THEN
      IF(RAD.LE.1.D-6) THEN
      VRLGRD(INDAD,NIND,NPION)=0.D0
      VRLGRD(INDAD,NIND+1,NPION)=0.D0
      VRLGRD(INDAD,NIND+2,NPION)=0.D0
      VRLGRD(INDAD,NIND+3,NPION)=0.D0
      VRLGRD(INDAD,NIND+4,NPION)=0.D0
      DVRLGR(INDAD,1,NIND,NPION)=0.D0
      DVRLGR(INDAD,2,NIND,NPION)=0.D0
      DVRLGR(INDAD,3,NIND,NPION)=0.D0
      DVRLGR(INDAD,1,NIND+1,NPION)=0.D0
      DVRLGR(INDAD,2,NIND+1,NPION)=0.D0
      DVRLGR(INDAD,3,NIND+1,NPION)=0.D0
      DVRLGR(INDAD,1,NIND+2,NPION)=0.D0
      DVRLGR(INDAD,2,NIND+2,NPION)=0.D0
      DVRLGR(INDAD,3,NIND+2,NPION)=0.D0
      DVRLGR(INDAD,1,NIND+3,NPION)=0.D0
      DVRLGR(INDAD,2,NIND+3,NPION)=0.D0
      DVRLGR(INDAD,3,NIND+3,NPION)=0.D0
      DVRLGR(INDAD,1,NIND+4,NPION)=0.D0
      DVRLGR(INDAD,2,NIND+4,NPION)=0.D0
      DVRLGR(INDAD,3,NIND+4,NPION)=0.D0
      NIND=NIND+5
      GO TO 1210
      ENDIF
      VRLGRD(INDAD,NIND,NPION)=(RFIF/2.D0)*VRLPR*
     &       (XCOR**2-YCOR**2)/RAD**2
      VRLGRD(INDAD,NIND+1,NPION)=(RFIV/2.D0)*VRLPR*
     &       (2.D0*ZCOR**2-XCOR**2-YCOR**2)/RAD**2
      VRLGRD(INDAD,NIND+2,NPION)=RFIF*VRLPR*XCOR*YCOR/RAD**2
      VRLGRD(INDAD,NIND+3,NPION)=RFIF*VRLPR*XCOR*ZCOR/RAD**2
      VRLGRD(INDAD,NIND+4,NPION)=RFIF*VRLPR*YCOR*ZCOR/RAD**2
      XYRSQ=(XCOR**2 - YCOR**2)/RAD**2
      DVRLGR(INDAD,1,NIND,NPION)=
     &          -(RFIF/2.D0)*(DVRLPR*XCOR*XYRSQ/RAD +
     &          VRLPR*2.D0*XCOR*(1.D0-XYRSQ)/RAD**2)
      DVRLGR(INDAD,2,NIND,NPION)=
     &          -(RFIF/2.D0)*(DVRLPR*YCOR*XYRSQ/RAD -
     &          VRLPR*2.D0*YCOR*(XYRSQ+1.D0)/RAD**2)
      DVRLGR(INDAD,3,NIND,NPION)=
     &         -(RFIF/2.D0)*(DVRLPR*ZCOR*XYRSQ/RAD -
     &         VRLPR*2.D0*ZCOR*XYRSQ/RAD**2)
      ZXYRSQ=(2.D0*ZCOR**2 - XCOR**2 - YCOR**2)/RAD**2
      DVRLGR(INDAD,1,NIND+1,NPION)=
     &         -(RFIV/2.D0)*(DVRLPR*XCOR*ZXYRSQ/RAD -
     &         VRLPR*2.D0*XCOR*(1.D0+ZXYRSQ)/RAD**2)
      DVRLGR(INDAD,2,NIND+1,NPION)=
     &         -(RFIV/2.D0)*(DVRLPR*YCOR*ZXYRSQ/RAD -
     &         VRLPR*2.D0*YCOR*(1.D0+ZXYRSQ)/RAD**2)
      DVRLGR(INDAD,3,NIND+1,NPION)=
     &         -(RFIV/2.D0)*(DVRLPR*ZCOR*ZXYRSQ/RAD +
     &         VRLPR*2.D0*ZCOR*(2.D0-ZXYRSQ)/RAD**2)
      XTYDRS=XCOR*YCOR/RAD**2
      DVRLGR(INDAD,1,NIND+2,NPION)=
     &         -RFIF*(DVRLPR*XCOR*XTYDRS/RAD +
     &         VRLPR*(YCOR-2.D0*XCOR*XTYDRS)/RAD**2)
      DVRLGR(INDAD,2,NIND+2,NPION)=
     &         -RFIF*(DVRLPR*YCOR*XTYDRS/RAD +
     &         VRLPR*(XCOR-2.D0*YCOR*XTYDRS)/RAD**2)
      DVRLGR(INDAD,3,NIND+2,NPION)=
     &         -RFIF*(DVRLPR*ZCOR*XTYDRS/RAD -
     &         VRLPR*2.D0*ZCOR*XTYDRS/RAD**2)
      XTZDRS=XCOR*ZCOR/RAD**2
      DVRLGR(INDAD,1,NIND+3,NPION)=
     &         -RFIF*(DVRLPR*XCOR*XTZDRS/RAD +
     &         VRLPR*(ZCOR-2.D0*XCOR*XTZDRS)/RAD**2)
      DVRLGR(INDAD,2,NIND+3,NPION)=
     &         -RFIF*(DVRLPR*YCOR*XTZDRS/RAD -
     &         VRLPR*2.D0*YCOR*XTZDRS/RAD**2)
      DVRLGR(INDAD,3,NIND+3,NPION)=
     &         -RFIF*(DVRLPR*ZCOR*XTZDRS/RAD +
     &         VRLPR*(XCOR-2.D0*ZCOR*XTZDRS)/RAD**2)
      YTZDRS=YCOR*ZCOR/RAD**2
      DVRLGR(INDAD,1,NIND+4,NPION)=
     &         -RFIF*(DVRLPR*XCOR*YTZDRS/RAD -
     &         VRLPR*2.D0*XCOR*YTZDRS/RAD**2)
      DVRLGR(INDAD,2,NIND+4,NPION)=
     &         -RFIF*(DVRLPR*YCOR*YTZDRS/RAD +
     &         VRLPR*(ZCOR-2.D0*YCOR*YTZDRS)/RAD**2)
      DVRLGR(INDAD,3,NIND+4,NPION)=
     &         -RFIF*(DVRLPR*ZCOR*YTZDRS/RAD +
     &         VRLPR*(YCOR-2.D0*ZCOR*YTZDRS)/RAD**2)
      NIND=NIND+5
      ENDIF
C========================================================================
C MOVE ONTO NEXT NON-LOCAL COMPONENT
C========================================================================
 1210 CONTINUE
C========================================================================
C INCREMENT ADDRESS COUNTER
C========================================================================
      INDAD=INDAD+1
C========================================================================
C END LOOPS OVER X,Y,Z RESPECTIVELY
C========================================================================
 1202 CONTINUE
 1201 CONTINUE
 1200 CONTINUE
C========================================================================
C CHECK THAT NRGRPT IS LARGE ENOUGH TO STORE DATA
C========================================================================
      IF(INDAD-1.GT.NRGRPT) THEN
      WRITE(*,1500) INDAD-1
      WRDIM=.TRUE.
 1500 FORMAT(' NRGRPT NOT BIG ENOUGH SHOULD BE:',I5)
      ENDIF
C========================================================================
C RECORD NUMBER OF REAL SPACE GRID POINTS TO BE SAMPLED FOR PRESENT ION
C========================================================================
      NRLPPI(NPION)=INDAD-1
C========================================================================
C INCREMENT ION COUNTER BY ONE
C========================================================================
      NPION=NPION+1
C========================================================================
C END ION LOOP
C========================================================================
 1100 CONTINUE
C========================================================================
C END SPECIES LOOP
C========================================================================
 1000 CONTINUE
C========================================================================
C CHECK THAT TOTAL NUMBER OF IONS HAS BEEN CORRECTLY SET
C========================================================================
      IF(NPION-1.GT.NIONST) THEN
      WRITE(*,1600) NPION-1
      WRDIM=.TRUE.
 1600 FORMAT(' WRONG VALUE FOR NIONST, SHOULD BE:',I3)
      ENDIF
      IF(WRDIM) STOP
      RETURN
      END

