C----------------------------------------------------------------------C   
C File : serlph.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   15-JUL-91  NODE INDEX REMOVED FROM ALL LOOKUP TABLES; I.E. EACH NODE 
C              STORES ONLY ITS OWN LOOKUP TABLE.   (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 SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,POSION,
     &           NRLPPI,NIONSP,RLCORE,CPHGRD,VKPT,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 PHASE FACTORS ON THE REAL SPACE
C GRID AROUND EACH ATOM IN PREPARATION FOR CALCULATING H*PSI AND THE
C FORCES ON THE ATOMS - THE ABSOLUTE VALUE OF THE PHASE FACTOR IS
C IRRELEVANT - ONLY THE RELATIVE VALUE BETWEEN GRIDPOINTS MATTERS
C
C========================================================================
      DIMENSION DIRC(3,3),RECC(3,3)
      DIMENSION POSION(3,NIONS,NSPEC),NRLPPI(NIONST)
      DIMENSION CPHGRD(NRGRPT,NIONST)
      DIMENSION NIONSP(NSPEC),RLCORE(NSPEC)
      DIMENSION VKPT(3)
      DIMENSION POSTMP(3)
C                                                                         
      DATA CITPI,TPI / (0.0,6.283185307179),6.283185307179 /
      DATA WRDIM /.FALSE./
C========================================================================
C INITIALISE THE INDEX FOR THE IONS, NPION, 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 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 NOW 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 PHASE FACTORS 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)
      IF(RAD.LE.RLCORE(NSP)) THEN
      CZPH=EXP(-CITPI*VKPT(3)*((1.D0*(NGRPZ+NZ-1))/(1.D0*NGZ)))
      CYPH=EXP(-CITPI*VKPT(2)*((1.D0*(NGRPY+NY-1))/(1.D0*NGY)))
      CXPH=EXP(-CITPI*VKPT(1)*((1.D0*(NGRPX+NX-1))/(1.D0*NGX)))
      CPHGRD(INDAD,NPION)=CXPH*CYPH*CZPH
C========================================================================
C INCREMENT ADDRESS COUNTER
C========================================================================
      INDAD=INDAD+1
      ENDIF
C========================================================================
C END LOOPS OVER X,Y,Z RESPECTIVELY
C========================================================================
 1202 CONTINUE
 1201 CONTINUE
 1200 CONTINUE
C========================================================================
C CHECK THAT POTENTIAL AND PHASE ARRAYS AGREE ON NUMBER OF GRID POINTS
C========================================================================
      IF(INDAD-1.NE.NRLPPI(NPION)) THEN
      WRITE(*,1500) NPION,INDAD-1,NRLPPI(NPION)
      WRDIM=.TRUE.
 1500 FORMAT(' # POTENTIAL AND PHASE POINTS DISAGREE:',
     &       ' ATOM:',I3,I5,I5)
      ENDIF
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


