C=======================================================================
C FILE: EWALTR.FOR
C
C HISTORY: 23/2/90: COPIED FROM FEWALD.FOR
C                   FEWALD() CHANGED TO EWALTR()
C                   EWSIF CHANGED TO SIGEW
C                   STATZD DELETED TOGETHER WITH ITS ASSOCIATED CODE
C                   NIONS CHANGED TO 4
C                   TRICLINIC GEOMETRY
C          25/2/90: OFF-DIAGONAL ELEMENTS IN STRESS TENSOR
C           8/3/90: SIGN OF STRESS-CONTRIB. CHANGED
C          24/4/90: SCALE IS MULTIPLIED BY 1.5 CHOSEN FOR BEST
C                     CONVERGENCE
C           6/5/90: CONDITIONED PRINTING OUT OF PARAMETERS
C
C          19/7/90: ELECT is changed to a more accuarate value
C
C         14/10/90: PARAMETER STATEMENT IS REMOVED. MATRICES ARE 
C                   FIRST DECLARED IN THE MAIN PROGRAM.
C=======================================================================
      SUBROUTINE EWALTR(POSION,EWIFC,EWRLEN,EWRCSS,EWRLSS,EWRCSI,
     &    DIRC,RECC,VOLC,SIGEW,TEWEN,NIONST,ICHARC,NIONCH,NICHSQ,
     &    CPHFX,CPHFY,CPHFZ,IPRINT,
     &    NIONS,NEWPTS,MAXCX,MAXCY,MAXCZ,MAXGPX,MAXGPY,MAXGPZ,NGPTS,
     &    RFORCE,CFORCE,FORCEG,FORSIG,ENERG)
C=======================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE EWALTR }}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE CALCULATES THE EWALD ENERGY DUE TO COULOMB ENERGY
C BETWEEN THE IONS AND THE NEUTRALISING BACKGROUND. THE COULOMB FORCES
C ON THE IONS DUE TO COULOMB INTERACTION WITH THE OTHER IONS AND THE
C STRESS ON THE UNIT CELL THAT IS RELATED TO THE CHANGE IN THE EWALD
C ENERGY ON CHANGING THE SIZE OF THE CELL
C          
C                         -G**2
C                     (-  -----)              2G  G
C           PI   --  e    4 eta [ --     iGR][  i  j    G**2         ]
C Sigma  = ----  >  ----------  [ >  Z  e   ][ -----  ( ---- ) - &   ]
C      ij  2eta  --   G**2      [ --  I     ][ G**2    4 eta      ij ]
C               G<>0( ----- )     I
C                     4 eta
C
C                                             D  D !
C                1    1/2 --            1/2    i  j!
C         + Vol ---eta    >   Z Z H'(eta   D)------!
C                2        --   I J            D**2 !
C                        I,J,l                     !D=R -R +l <>0
C                                                      J  I
C
C             PI  [ --    ]2
C         + ----- [ >  Z  ]  &
C           2 eta [ --  I ]   ij
C                   I
C=======================================================================
C THE EWALD METHOD OF PERFORMING COULOMB SUMS OVER A LATTICE IS
C EXPLAINED IN MANY TEXTS, ONE IS 'DYNAMICAL THEORY OF CRYSTAL LATTICES'
C BY BORN AND HUANG. THE ONLY THING TO BE CAREFUL ABOUT IS THAT IN A
C TOTAL ENERGY CALCULATION THE EWALD ENERGY INCLUDES THE ENERGY OF
C INTERACTION BETWEEN THE IONS AND A NEUTRALISING BACKGROUND.
C THE EWALD METHOD WORKS BY WRITING THE COULOMB ENERGY AS A SUM OF A
C REAL SPACE PART (WHICH IS SUMMED OVER ALL SPACE) AND A RECIPROCAL
C SPACE PART (WHICH IS SUMMED OVER ALL THE RECIPROCAL LATTICE POINTS).
C THIS PART OF THE METHOD IS EXACT AND DOESN'T REDUCE THE AMOUNT OF
C WORK NEEDED TO PERFORM THE SUMMATION. HOWEVER, THE REAL SPACE
C CONTRIBUTIONS FALL OFF RAPIDLY WITH INCREASING DISTANCE SO THAT ONLY
C A FINITE NUMBER OF UNIT CELLS NEED TO BE INCLUDED IN THE SUM AND THE
C RECIPROCAL SPACE CONTRIBUTIONS FALL OFF RAPIDLY WITH THE MAGNITUDE OF
C THE RECIPROCAL LATTICE VECTOR SO THAT ONLY A FINITE NUMBER OF
C RECIPROCAL LATTICE POINTS NEED TO BE INCLUDED IN THE RECIPROCAL SPACE
C SUM
C A VARIABLE (SCALE IN THIS PROGRAM) DETERMINES HOW RAPIDLY THE EWALD
C SUMS FALL OFF WITH INCREASING DISTANCE IN REAL AND RECIPROCAL SPACE.
C THE NUMBER OF REAL SPACE CELLS SUMMED OVER IS DETERMINED BY MAXCX,Y,Z
C AND THE NUMBER OF RECIPROCAL SPACE CELLS BY MAXGPX,Y,Z. THE VALUE OF
C SCALE SHOULD BE CHOSEN TO MINIMISE THE TOTAL NUMBER OF OPERATIONS
C NEEDED TO COMPUTE THE COULOMB LATTICE SUMS.
C=======================================================================
C RATHER THAN COMPUTE THE VALUES OF THE EWALD INTEGRALS THIS PROGRAM
C INTERPOLATES THE VALUES OF THE INTEGRALS FROM A SET OF DATA ARRAYS
C EWRLEN,EWRCSS,EWRLSS AND EWRCSI
C=======================================================================
      IMPLICIT COMPLEX (C)
C=======================================================================
C
C             TERMS IN THE PARAMETER STATEMENT
C
C MAXCX,MAXCY,MAXCZ = THE NUMBER OF CELLS IN EACH DIRECTION FROM THE
C CENTRAL CELL INCLUDED IN THE EVALUATION OF THE REAL SPACE EWALD SUMS
C THE VALUES OF THESE PARAMETERS SHOULD BE CHOSEN SO THAT
C MAXCX(Y,Z)*SQRT(PI)*SIZEX(Y,Z)/MIN(SIZEX,SIZEY,SIZEZ) ARE ALL GREATER
C THAN 3, (THE LARGER THE VALUE THE CLOSER TO THE CORRECT COULOMB SUM)
C MAXGPX,MAXGPY,MAXGPZ = THE NUMBER OF RECIPROCAL LATTICE POINTS IN EACH
C DIRECTION FROM G=0 INCLUDED IN THE EVALUATION OF THE RECIPROCAL SPACE
C EWALD SUMS. THE VALUES OF THESE PARAMETERS SHOULD BE CHOSEN SO THAT
C NGPTX(Y,Z)*MIN(SIZEX,SIZEY,SIZEZ)/(SQRT(PI)*SIZEX(Y,X) ARE ALL GREATER
C THAN 3, (THE LARGER THE VALUE THE CLOSER TO THE CORRECT COULOMB SUM)
C ARGMRE IS THE ARGUMENT CORRESPONDING TO THE LAST ELEMENT IN THE REAL
C SPACE EWALD INTEGRAL DATA ARRAYS
C ARGMRC IS THE ARGUMENT CORRESPONDING TO THE LAST ELEMENT IN THE
C RECIPROCAL SPACE EWALD INTEGRAL DATA ARRAYS
C NEWPTS IS THE NUMBER OF VALUES IN THE EWALD INTEGRAL DATA ARRAYS
C=======================================================================
C      PARAMETER(MAXCX=4,MAXCY=4,MAXCZ=2,NIONS=32,
C     &   NCELLX=(2*MAXCX)+1,NCELLY=(2*MAXCY)+1,NCELLZ=(2*MAXCZ)+1,
C     &   NCELLS=NCELLX*NCELLY*NCELLZ,
C     &   MAXGPX=4,MAXGPY=4,MAXGPZ=6,
C     &   NGPTX=(2*MAXGPX)+1,NGPTY=(2*MAXGPY)+1,NGPTZ=(2*MAXGPZ)+1,
C     &   NGPTS=NGPTX*NGPTY*NGPTZ,
C     &   NEWPTS=1000,ARGMRE=4.0,ARGMRC=4.0)
C
C
      DIMENSION POSION(3,*)
      DIMENSION EWIFC(3,*)
      DIMENSION ICHARC(*)
      DIMENSION EWRLEN(*),EWRCSS(*),EWRLSS(*),EWRCSI(*)
      DIMENSION DIRC(3,*)
      DIMENSION RECC(3,*)
      DIMENSION SIGEW(*)
      DIMENSION CPHFX(*),CPHFY(*),CPHFZ(*)
C=======================================================================
C
C                  DIMENSION STATEMENTS
C
C RFORCE(L,A,B) = THE CONTRIBUTION TO THE TOTAL COULOMB FORCE IN THE 'L'
C        DIRECTION ON ION A FROM THE PERIODIC ARRAY OF IONS B
C        FROM THE PART OF THE EWALD SUMS PERFORMED IN REAL SPACE
C CFORCE(L,A,B) = THE CONTRIBUTION TO THE TOTAL COULOMB FORCE IN THE 'L'
C        DIRECTION ON ION A FROM THE PERIODIC ARRAY OF IONS B
C        FROM THE PART OF THE EWALD SUMS PERFORMED IN RECIPROCAL SPACE
C FORCEG = THE RECIPROCAL SPACE LATTICE VECTOR MULTIPLIED BY AN EWALD
C        INTEGRAL (USED FOR CALCULATING THE FORCES ON THE IONS)
C FORSIG = THE COMBINATION OF EWALD INTEGRALS AND FUNCTIONS OF THE
C        RECIPROCAL LATTICE VECTOR USED FOR CALCULATING THE FORCES ON
C        THE UNIT CELL
C RSIFOR = THE FORCE ON THE UNIT CELL FROM THE REAL SPACE PARTS OF THE
C        EWALD SUMMATIONS
C CSIFOR = THE FORCE ON THE UNIT CELL FROM THE RECIPROCAL SPACE PARTS OF
C        THE EWALD SUMMATIONS
C ENERG = THE EWALD INTEGRAL USED TO CALCULATE THE CONTRIBUTION TO THE
C        TOTAL ENERGY FROM THE RECIPROCAL SPACE PART OF THE EWALD SUM
C=======================================================================
      DIMENSION RFORCE(3,NIONS,NIONS),CFORCE(3,NIONS,NIONS)
      DIMENSION FORCEG(3,NGPTS),FORSIG(6,NGPTS),RSIFOR(6),CSIFOR(6)
      DIMENSION ENERG(NGPTS)
      DATA   ARGMRE/4.0/,ARGMRC/4.0/
C=======================================================================
C ELECT = (THE CHARGE ON THE ION * THE ELECTRONIC CHARGE)/(4*PI * THE
C          PERMITTIVITY OF FREE SPACE)
C=======================================================================
      DATA ELECT/14.39976868/
      DATA HALF,SIXTH / 0.5,0.166666666667 /
      DATA PI,TPI,CITPI/3.1415926536,6.2831853072,(0.0,6.2831853072)/
C=======================================================================
C=======================================================================
      NCELLX=(2*MAXCX)+1
      NCELLY=(2*MAXCY)+1
      NCELLZ=(2*MAXCZ)+1
      NCELLS=NCELLX*NCELLY*NCELLZ
      NGPTX=(2*MAXGPX)+1
      NGPTY=(2*MAXGPY)+1
      NGPTZ=(2*MAXGPZ)+1
      IF(IPRINT.GT.1) THEN
        WRITE(*,1856) MAXCX,MAXCX,MAXCY,MAXCY,MAXCZ,MAXCZ
        WRITE(*,1857) MAXGPX,MAXGPX,MAXGPY,MAXGPY,MAXGPZ,MAXGPZ
      END IF
 1856 FORMAT(/,1X,'EWALD  R E A L   SPACE CELLS: ',
     &      '(-',I2,',',I2,')x(-',I2,',',I2,')x(-',I2,',',I2,')')
 1857 FORMAT(1X,'EWALD RECIPROCAL SPACE CELLS: ',
     &      '(-',I2,',',I2,')x(-',I2,',',I2,')x(-',I2,',',I2,')')
C
       IF(IPRINT.GE.2)WRITE(*,*)' VOLC=',VOLC
C=======================================================================
C INITIALISE THE VALUES OF THE DATA ARRAYS
C=======================================================================
      DO 5001 N=1,NIONST
        DO 5002 NN=1,NIONST
          DO 5003 M=1,3
            CFORCE(M,NN,N)=(0.0,0.0)
            RFORCE(M,NN,N)=0.0
 5003     CONTINUE
 5002   CONTINUE
 5001 CONTINUE
      DO 5004 N=1,NIONST
        DO 5005 M=1,3
          EWIFC(M,N)=0.0
 5005   CONTINUE
 5004 CONTINUE
      DO 5006 M=1,6
        RSIFOR(M)=0.0
        CSIFOR(M)=(0.0,0.0)
        SIGEW(M)=0.0
 5006 CONTINUE
      TEWEN=0.0
      REWEN=0.0
      CEWEN=(0.0,0.0)
C=======================================================================
C CHOOSE THE VARIABLE 'SCALE' THAT DETERMINES THE RATE OF DECAY OF THE
C REAL AND RECIPROCAL PARTS OF THE EWALD SUMS. THE CHOICE OF SCALE USED
C BELOW IS GOOD FOR SYSTEMS WHERE THE LENGTHS OF THE LATTICE VECTORS
C ARE SIMILAR, FOR ODDER SHAPE SYSTEMS A DIFFERENT CHOICE FOR THE VALUE
C OF SCALE MAY REDUCE THE AMOUNT OF WORK INVOLVED IN CALCULATING THE
C COULOMB LATTICE SUMS. DONT FORGET THAT WHATEVER VALUE IS USED FOR
C SCALE THE VALUES OF MAXCX,Y,Z AND MAXGPX,Y,Z MUST BE SUFFICIENTLY
C LARGE TO OBTAIN THE LEVEL OF CONVERGENCE DESIRED
C=======================================================================
C was as follows
C
C      SIZMIN=SQRT(MIN( (DIRC(1,1)**2+DIRC(1,2)**2+DIRC(1,3)**2),
C    &                  (DIRC(2,1)**2+DIRC(2,2)**2+DIRC(2,3)**2),
C    &                  (DIRC(3,1)**2+DIRC(3,2)**2+DIRC(3,3)**2) ))
C
C      SIZMIN=SQRT( (DIRC(1,1)**2+DIRC(1,2)**2+DIRC(1,3)**2)
C     &            +(DIRC(2,1)**2+DIRC(2,2)**2+DIRC(2,3)**2)
C     &            +(DIRC(3,1)**2+DIRC(3,2)**2+DIRC(3,3)**2) )/3.0
C
      SIZMIN=    ( (DIRC(1,1)**2+DIRC(1,2)**2+DIRC(1,3)**2)
     &            *(DIRC(2,1)**2+DIRC(2,2)**2+DIRC(2,3)**2)
     &            *(DIRC(3,1)**2+DIRC(3,2)**2+DIRC(3,3)**2) )**(1./6.)
      SCALE=SQRT(PI)/SIZMIN*1.0
      PIDSCA=PI/SCALE
      IF(IPRINT.GT.1) WRITE(*,*)'EWALD ETA: ', SCALE
C=======================================================================
C     CHECK TO SEE WHETHER THE PARAMETERS MAXCX,..., ARE BIG ENOUGH
C     SUGGESTION WILLBE MADE.		             X.WENG 24-OCT-1990
C=======================================================================
      CALL MCHECK(DIRC,RECC,ARGMRE,ARGMRC,MAXCX,MAXCY,MAXCZ,
     &            MAXGPX,MAXGPY,MAXGPZ,SCALE,PIDSCA,IPRINT)
C=======================================================================
C CALCULATE THE CONTRIBUTION TO THE TOTAL ENERGY, THE FORCES ON THE IONS
C AND THE FORCES ON THE UNIT CELL DUE TO THE REAL SPACE PART OF THE
C EWALD SUMS
C=======================================================================
      DO 5010 NI=1,NIONST
      DO 5011 NNI=NI,NIONST
      ICHSQ=ICHARC(NI)*ICHARC(NNI)
C=======================================================================
C THE ENERGY OF INTERACTION BETWEEN EACH PAIR OF IONS MUST BE DIVIDED
C EQUALLY BETWEEN THE TWO IONS. SINCE THE LOOPS OVER THE IONS GIVE THE
C ENERGY OF INTERATION OF ION NI WITH THE PERIODIC ARRAY OF IONS NNI
C WHERE THE INDEX NNI IS GREATER THAN OR EQUAL TO NI WE INCLUDE HALF
C THE ENERGY IF NI IS EQUAL TO NNI BUT ALL OF OF IT IF NNI IS GREATER
C THAN NI TO ACCOUNT FOR ENERGY DUE TO ION NNI INTERACTING WITH THE
C ARRAY OF IONS NI
C=======================================================================
      DCOUNT=ICHSQ
      IF(NNI.EQ.NI) DCOUNT=DCOUNT/2.0
C=======================================================================
C INITIALISE THE DISTANCES IN THE X,Y,Z DIRECTIONS BETWEEN ION NI AND IO
C NNI IN THE CORNER CELL OF THE BLOCK OF UNIT CELLS INCLUDED IN THE REAL
C SPACE PART OF THE EWALD SUM
C=======================================================================
      A1DIFF=POSION(1,NI)-POSION(1,NNI)+MAXCX
      A2DIFF=POSION(2,NI)-POSION(2,NNI)+MAXCY
      A3DIFF=POSION(3,NI)-POSION(3,NNI)+MAXCZ
      XDIFF=DIRC(1,1)*A1DIFF+DIRC(2,1)*A2DIFF+DIRC(3,1)*A3DIFF
      YDIFF=DIRC(1,2)*A1DIFF+DIRC(2,2)*A2DIFF+DIRC(3,2)*A3DIFF
      ZDIFF=DIRC(1,3)*A1DIFF+DIRC(2,3)*A2DIFF+DIRC(3,3)*A3DIFF
C=======================================================================
C CALCULATE THE REAL SPACE CONTRIBUTION TO THE TOTAL ENERGY ETC BY
C SUMMING THE CONTRIBUTIONS DUE TO ION NI INTERACTING WITH ALL THE
C IONS NNI IN THE BLOCK OF UNIT CELLS RETAINED IN THE REAL SPACE SUM
C=======================================================================
      DO 5100 N=1,NCELLZ
      DO 5101 NN=1,NCELLY
      DO 5102 NNN=1,NCELLX
C=======================================================================
C CONVERT THE DISTANCE BETWEEN THE IONS INTO A DIMENSIONLESS FORM BY
C MULTIPLYING BY 'SCALE'
C=======================================================================
      ARG=SQRT(XDIFF**2+YDIFF**2+ZDIFF**2)*SCALE
C=======================================================================
C CONVERT THE DIMENSIONLESS DISTANCE TO AN ADDRESS IN THE EWALD INTEGRAL
C ARRAYS TO INTERPOLATE THE VALUES OF THE EWALD INTEGRALS, IF THE
C DISTANCE IS OUTSIDE THE RANGE OF VALUES IN THE EWALD INTEGRAL DATA
C ARRAYS MOVE ONTO ION NNI IN THE NEXT UNIT CELL
C=======================================================================
      IF(ARG.GT.3.98) GO TO 5103
      ARG=(ARG*NEWPTS/ARGMRE)+2.0
      NADDR=INT(ARG)
      REM=AMOD(ARG,1.0)
C=======================================================================
C INTERPOLATE THE VALUE OF THE EWALD INTEGRAL USED IN THE REAL SPACE
C PART OF THE EWALD SUMMATION FOR THE FORCES ON THE IONS AND THE UNIT
C CELL
C=======================================================================
      V1=EWRLSS(NADDR-1)
      V2=EWRLSS(NADDR)
      V3=EWRLSS(NADDR+1)
      V4=EWRLSS(NADDR+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
      EWRL=(T0+(T1*REM)+(T2*(REM**2))+(T3*(REM**3)))
C=======================================================================
C ADD THE CONTRIBUTION TO THE TOTAL FORCE ON ION NI DUE TO THE PERIODIC
C ARRAY OF IONS NNI FROM ION NNI IN THE PRESENT REAL SPACE CELL
C=======================================================================
      RFORCE(1,NNI,NI)=RFORCE(1,NNI,NI)+XDIFF*EWRL
      RFORCE(2,NNI,NI)=RFORCE(2,NNI,NI)+YDIFF*EWRL
      RFORCE(3,NNI,NI)=RFORCE(3,NNI,NI)+ZDIFF*EWRL
C=======================================================================
C ADD THE CONTRIBUTION TO THE TOTAL FORCE ON THE UNIT CELL DUE TO THE
C COULOMB INTERACTION BETWEEN ION NI AND ION NNI IN THE PRESENT REAL
C SPACE CELL
C=======================================================================
      D1X=DIRC(1,1)*A1DIFF
      D1Y=DIRC(1,2)*A1DIFF
      D1Z=DIRC(1,3)*A1DIFF
      D2X=DIRC(2,1)*A2DIFF
      D2Y=DIRC(2,2)*A2DIFF
      D2Z=DIRC(2,3)*A2DIFF
      D3X=DIRC(3,1)*A3DIFF
      D3Y=DIRC(3,2)*A3DIFF
      D3Z=DIRC(3,3)*A3DIFF
C
      D1=D1X+D2X+D3X
      D2=D1Y+D2Y+D3Y
      D3=D1Z+D2Z+D3Z
C
      RSIFOR(1)=RSIFOR(1)-EWRL*(D1**2)*DCOUNT
      RSIFOR(2)=RSIFOR(2)-EWRL*(D2**2)*DCOUNT
      RSIFOR(3)=RSIFOR(3)-EWRL*(D3**2)*DCOUNT
      RSIFOR(4)=RSIFOR(4)-EWRL*(D2*D3)*DCOUNT
      RSIFOR(5)=RSIFOR(5)-EWRL*(D3*D1)*DCOUNT
      RSIFOR(6)=RSIFOR(6)-EWRL*(D1*D2)*DCOUNT
C=======================================================================
C INTERPOLATE THE VALUE OF THE EWALD INTEGRAL USED IN THE REAL SPACE
C PART OF THE EWALD SUMMATION FOR THE FORCES ON THE IONS AND THE UNIT
C CELL
C=======================================================================
      V1=EWRLEN(NADDR-1)
      V2=EWRLEN(NADDR)
      V3=EWRLEN(NADDR+1)
      V4=EWRLEN(NADDR+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
      VALUE=(T0+(T1*REM)+(T2*(REM**2))+(T3*(REM**3)))
C=======================================================================
C ADD THE CONTRIBUTION TO THE TOTAL ENERGY DUE TO THE COULOMB
C INTERACTION BETWEEN ION NI AND ION NNI IN THE PRESENT REAL SPACE CELL
C=======================================================================
      REWEN=REWEN+VALUE*DCOUNT
 5103 CONTINUE
C=======================================================================
C MOVE ONTO ION NNI IN THE NEXT REAL SPACE CELL IN THE A1 DIRECTION
C=======================================================================
      XDIFF=XDIFF-DIRC(1,1)
      YDIFF=YDIFF-DIRC(1,2)
      ZDIFF=ZDIFF-DIRC(1,3)
      A1DIFF=A1DIFF-1.0
 5102 CONTINUE
C=======================================================================
C MOVE ONTO THE NEXT ROW OF CELLS IN THE A2 DIRECTION AND BACK TO THE
C FIRST CELL IN THE ROW
C=======================================================================
      XDIFF=XDIFF+DIRC(1,1)*NCELLX-DIRC(2,1)
      YDIFF=YDIFF+DIRC(1,2)*NCELLX-DIRC(2,2)
      ZDIFF=ZDIFF+DIRC(1,3)*NCELLX-DIRC(2,3)
      A1DIFF=A1DIFF+NCELLX
      A2DIFF=A2DIFF-1.0
 5101 CONTINUE
C=======================================================================
C MOVE ONTO THE NEXT PLANE OF CELLS IN THE A3 DIRECTION AND BACK TO THE
C FIRST CELL IN THE PLANE
C=======================================================================
      XDIFF=XDIFF+DIRC(2,1)*NCELLY-DIRC(3,1)
      YDIFF=YDIFF+DIRC(2,2)*NCELLY-DIRC(3,2)
      ZDIFF=ZDIFF+DIRC(2,3)*NCELLY-DIRC(3,3)
      A2DIFF=A2DIFF+NCELLY
      A3DIFF=A3DIFF-1.0
 5100 CONTINUE
C=======================================================================
C END OF THE LOOP OVER ION NNI, MOVE ONTO THE NEXT ION
C=======================================================================
 5011 CONTINUE
C=======================================================================
C END OF THE LOOP OVER ION NI, MOVE ONTO THE NEXT ION
C=======================================================================
 5010 CONTINUE
C=======================================================================
C CALCULATE THE SCALING FACTORS THAT WILL BE USED TO CONVERT THE VALUES
C OF THE FORCES AND THE ENERGY CALCULATED ABOVE TO REAL FORCES AND
C ENERGIES. THE FACTORS TAKE ACCOUNT OF THE NUMERICAL CONSTANTS
C INTRODUCED BY THE EWALD SEPARATION OF THE COULOMB LATTICE SUMS, THE
C RESCALING OF ALL OF THE EWALD INTEGRALS ONTO THE RANGE ZERO TO ONE
C AS WELL PUTTING IN THE ELECTROSTATIC CONSTANTS VIA 'ELECT'
C=======================================================================
      RSCALE=(SCALE**3)*4.0*ELECT/SQRT(PI)
      RENSCA=SCALE*2.0*ELECT/SQRT(PI)
C=======================================================================
C
C                      ***********
C NOW MOVE ONTO THE RECIPROCAL SPACE PART OF THE EWALD SUMMATIONS
C                      ***********
C
C=======================================================================
C INITIALISE THE X,Y,Z RECIPROCAL LATTICE VECTORS TO THE RECIPROCAL
C LATTICE POINT AT THE CORNER OF THE GRID USED FOR THE RECIPROCAL SPACE
C CONTRIBUTIONS TO THE EWALD SUMS
C=======================================================================
      G1X=-(RECC(1,1)*MAXGPX)/TPI
      G1Y=-(RECC(1,2)*MAXGPX)/TPI
      G1Z=-(RECC(1,3)*MAXGPX)/TPI
      G2X=-(RECC(2,1)*MAXGPY)/TPI
      G2Y=-(RECC(2,2)*MAXGPY)/TPI
      G2Z=-(RECC(2,3)*MAXGPY)/TPI
      G3X=-(RECC(3,1)*MAXGPZ)/TPI
      G3Y=-(RECC(3,2)*MAXGPZ)/TPI
      G3Z=-(RECC(3,3)*MAXGPZ)/TPI
      GX=G1X+G2X+G3X
      GY=G1Y+G2Y+G3Y
      GZ=G1Z+G2Z+G3Z
C=======================================================================
C CALCULATE THE SCALING FACTORS THAT WILL BE USED TO CALCULATE THE
C FORCES ON THE UNIT CELL. THE FACTORS TAKE ACCOUNT OF THE NUMERICAL
C CONSTANTS INTRODUCED BY THE EWALD SEPARATION OF THE COULOMB LATTICE
C SUMS, THE RESCALING OF ALL OF THE EWALD INTEGRALS ONTO THE RANGE ZERO
C TO ONE AS WELL AS PUTTING IN THE ELECTROSTATIC CONSTANTS VIA 'ELECT'
C=======================================================================
      GSISC1=-4.0*(PI**3)*ELECT/(VOLC*SCALE**4)
      GSISC2=2.0*PI*ELECT/(VOLC*SCALE**2)
C=======================================================================
C NG IS AN INDEX USED TO LABEL THE RECIPROCAL LATTICE VECTORS
C INITIALISE TO ONE
C=======================================================================
      NG=1
C=======================================================================
C INTERPOLATE THE VALUES OF THE RECIPROCAL SPACE EWALD INTEGRALS AT THE
C GRID OF RECIPROCAL LATTICE POINTS AND CALCULATE FUNCTIONS OF THE EWALD
C INTEGRALS AND RECIPROCAL LATTICE VECTORS TO BE USED FOR CALULATING THE
C FORCES ON THE IONS AND THE UNIT CELL AND THE TOTAL ENERGY
C=======================================================================
      DO 5200 N=1,NGPTZ
      DO 5201 NN=1,NGPTY
      DO 5202 NNN=1,NGPTX
C=======================================================================
C CALCULATE THE DIMENSIONLESS LENGTH OF THE RECIPROCAL LATTICE VECTOR
C=======================================================================
      GX=G1X+G2X+G3X
      GY=G1Y+G2Y+G3Y
      GZ=G1Z+G2Z+G3Z
      GARG=SQRT(GX**2+GY**2+GZ**2)*PIDSCA
C=======================================================================
C CONVERT THE DIMENSIONLESS LENGTH TO AN ADDRESS IN THE EWALD INTEGRAL
C ARRAYS TO INTERPOLATE THE VALUES OF THE EWALD INTEGRALS, IF THE
C DISTANCE IS OUTSIDE THE RANGE OF VALUES IN THE EWALD INTEGRAL DATA
C ARRAYS SET THE VALUES OF THE EWALD INTEGRALS TO ZERO
C=======================================================================
      IF(GARG.GT.3.98) THEN
      EWRC=0.0
      EWSI=0.0
      ELSE
      GARG=(GARG*NEWPTS/ARGMRC)+2.0
      NADDR=INT(GARG)
      REM=AMOD(GARG,1.0)
C=======================================================================
C INTERPOLATE THE VALUES OF THE EWALD INTEGRALS USED IN THE RECIPROCAL
C SPACE PART OF THE EWALD SUMMATION
C=======================================================================
      V1=EWRCSS(NADDR-1)
      V2=EWRCSS(NADDR)
      V3=EWRCSS(NADDR+1)
      V4=EWRCSS(NADDR+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
      EWRC=(T0+(T1*REM)+(T2*(REM**2))+(T3*(REM**3)))
C
      V1=EWRCSI(NADDR-1)
      V2=EWRCSI(NADDR)
      V3=EWRCSI(NADDR+1)
      V4=EWRCSI(NADDR+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
      EWSI=(T0+(T1*REM)+(T2*(REM**2))+(T3*(REM**3)))
      ENDIF
C=======================================================================
C USING THE VALUES OF THE EWALD INTEGRALS FOR THIS RECIPROCAL LATTICE
C VECTOR CALCULATE THE FUNCTIONS OF THE RECIPROCAL LATTICE VECTOR AND
C THE EWALD INTEGRALS USED TO CALCULATE THE FORCES ETC
C=======================================================================
      FORCEG(1,NG)=GX*EWRC
      FORCEG(2,NG)=GY*EWRC
      FORCEG(3,NG)=GZ*EWRC
      G1=GX
      G2=GY
      G3=GZ
      FORSIG(1,NG)=-EWSI*(G1**2)*GSISC1-(EWRC*GSISC2)
      FORSIG(2,NG)=-EWSI*(G2**2)*GSISC1-(EWRC*GSISC2)
      FORSIG(3,NG)=-EWSI*(G3**2)*GSISC1-(EWRC*GSISC2)
      FORSIG(4,NG)=-EWSI*(G2*G3)*GSISC1
      FORSIG(5,NG)=-EWSI*(G3*G1)*GSISC1
      FORSIG(6,NG)=-EWSI*(G1*G2)*GSISC1
      ENERG(NG)=EWRC
C=======================================================================
C INCREASE THE INDEX NG BY ONE FOR THE NEXT RECIPROCAL LATTICE VECTOR
C=======================================================================
      NG=NG+1
C=======================================================================
C MOVE ONTO THE NEXT RECIPROCAL LATTICE POINT IN THE X DIRECTION
C=======================================================================
      G1X=G1X+RECC(1,1)/TPI
      G1Y=G1Y+RECC(1,2)/TPI
      G1Z=G1Z+RECC(1,3)/TPI
 5202 CONTINUE
C=======================================================================
C MOVE ONTO THE NEXT ROW OF RECIPROCAL LATTICE VECTORS IN THE Y
C DIRECTION AND BACK TO THE FIRST POINT IN THE ROW
C=======================================================================
      G1X=-(RECC(1,1)*MAXGPX)/TPI
      G1Y=-(RECC(1,2)*MAXGPX)/TPI
      G1Z=-(RECC(1,3)*MAXGPX)/TPI
      G2X=G2X+RECC(2,1)/TPI
      G2Y=G2Y+RECC(2,2)/TPI
      G2Z=G2Z+RECC(2,3)/TPI
 5201 CONTINUE
C=======================================================================
C MOVE ONTO THE NEXT PLANE OF RECIPROCAL LATTICE VECTORS IN THE Z
C DIRECTION AND BACK TO THE POINT IN THE CORNER OF THE PLANE
C=======================================================================
      G2X=-(RECC(2,1)*MAXGPY)/TPI
      G2Y=-(RECC(2,2)*MAXGPY)/TPI
      G2Z=-(RECC(2,3)*MAXGPY)/TPI
      G3X=G3X+RECC(3,1)/TPI
      G3Y=G3Y+RECC(3,2)/TPI
      G3Z=G3Z+RECC(3,3)/TPI
 5200 CONTINUE
C=======================================================================
C USING THE FUNCTIONS COMPUTED ABOVE CALCULATE THE FORCES ON THE IONS
C AND ON THE UNIT CELL AND THE CONTRIBUTION TO THE TOTAL ENERGY WHICH
C ARISES FROM THE RECIPROCAL SPACE PART OF THE EWALD SUMMATION
C=======================================================================
C======================================================================
C THE FOLLOWING ARE THE LOOPS TO CALCULATE THE FORCES BETWEEN ALL PAIRS
C OF ATOMS IN THE UNIT CELL AND THE ENERGY
C=======================================================================
      DO 5300 NI=1,NIONST
      DO 5301 NNI=NI,NIONST
      ICHSQ=ICHARC(NI)*ICHARC(NNI)
C=======================================================================
C AS BEFORE, DCOUNT ACCOUNTS FOR THE ENERGY OF INTERACTION BETWEEN PAIRS
C OF IONS FOR THE RESTRICTED SUM NNI GREATER OR EQUAL TO NI
C=======================================================================
      DCOUNT=ICHSQ
      IF(NI.EQ.NNI) DCOUNT=DCOUNT/2.0
C=======================================================================
C CX,Y,ZPHF ARE THE CHANGES IN PHASE FACTORS ON MOVING ONE RECIPROCAL
C LATTICE VECTOR IN THE X,Y,Z DIRECTIONS, RESPECTIVELY
C=======================================================================
      CXPHF=(POSION(1,NI)-POSION(1,NNI))*CITPI
      CYPHF=(POSION(2,NI)-POSION(2,NNI))*CITPI
      CZPHF=(POSION(3,NI)-POSION(3,NNI))*CITPI
C=======================================================================
C INITIALISE THE TOTAL PHASE FACTOR TO THE VALUE FOR THE RECIPROCAL
C LATTICE POINT AT THE CORNER OF THE RECIPROCAL LATTICE GRID AND
C INITIALISE THE RECIPROCAL LATTICE VECTOR INDEX TO ONE
C=======================================================================
      CPHF=-CXPHF*MAXGPX
      DO 5330 N=1,NGPTX
      CPHFX(N)=CEXP(CPHF)
      CPHF=CPHF+CXPHF
 5330 CONTINUE
      CPHF=-CYPHF*MAXGPY
      DO 5331 N=1,NGPTY
      CPHFY(N)=CEXP(CPHF)
      CPHF=CPHF+CYPHF
 5331 CONTINUE
      CPHF=-CZPHF*MAXGPZ
      DO 5332 N=1,NGPTZ
      CPHFZ(N)=CEXP(CPHF)
      CPHF=CPHF+CZPHF
 5332 CONTINUE
      CFPHFZ=-CITPI*MAXGPZ
      NG=1
C=======================================================================
C FOR EACH PAIR OF ATOMS SUM THE FORCES AND THE CONTRIBUTION TO THE
C TOTAL ENERGY OVER THE GRID OF RECIPROCAL LATTICE POINTS
C=======================================================================
      DO 5310 N=1,NGPTZ
      DO 5311 NN=1,NGPTY
      CPHFZY=CPHFY(NN)*CPHFZ(N)
      DO 5312 NNN=1,NGPTX
C=======================================================================
C CALCULATE THE EXPONENTIAL OF THE PHASE FACTOR
C=======================================================================
      CEXPHF=CPHFX(NNN)*CPHFZY
C=======================================================================
C CALCULATE THE CONTRIBUTION TO THE FORCE BETWEEN ION NI AND IONS NNI,
C AND THE FORCE ON THE UNIT CELL DUE TO THE INTERACTION BETWEEN THESE
C IONS FROM THE RECIPROCAL LATTICE PART OF THE EWALD SUM AT THE PRESENT
C K POINT
C=======================================================================
       DO 5320 M=1,3
       CFORCE(M,NNI,NI)=CFORCE(M,NNI,NI)+FORCEG(M,NG)*CEXPHF
       CSIFOR(M)=CSIFOR(M)+(FORSIG(M,NG)*CEXPHF)*DCOUNT
 5320  CONTINUE
       DO 5325 M=4,6
       CSIFOR(M)=CSIFOR(M)+(FORSIG(M,NG)*CEXPHF)*DCOUNT
 5325  CONTINUE
C=======================================================================
C CALCULATE THE CONTRIBUTION TO THE TOTAL ENERGY DUE TO THE INTERACTION
C BETWEEN ION NI AND IONS NNI FROM THE RECIPROCAL SPACE PART OF THE
C EWALD SUM AT THE PRESENT K POINT
C=======================================================================
      CEWEN=CEWEN+(ENERG(NG)*CEXPHF)*DCOUNT
C=======================================================================
C INCREASE THE RECIPROCAL LATTICE POINT INDEX BY ONE AND CHANGE THE
C PHASE FACTOR BY THE PHASE CORRESPONDING TO MOVING ONE RECIPROCAL
C LATTICE VECTOR IN THE X DIRECTION, IE TO THE NEXT RECIPROCAL LATTICE
C POINT IN THE X DIRECTION
C=======================================================================
      NG=NG+1
 5312 CONTINUE
C=======================================================================
C MOVE TO THE NEXT ROW OF RECIPROCAL LATTICE POINTS IN THE Y DIRECTION
C AND THE FIRST RECIPROCAL LATTICE POINT IN THE ROW
C=======================================================================
 5311 CONTINUE
C=======================================================================
C MOVE TO THE NEXT PLANE OF RECIPROCAL LATTICE POINTS IN THE Z DIRECTION
C AND THE RECIPROCAL LATTICE POINT AT THE CORNER OF THE PLANE
C=======================================================================
 5310 CONTINUE
C=======================================================================
C MOVE ONTO THE NEXT ION NNI
C=======================================================================
 5301 CONTINUE
C=======================================================================
C MOVE ONTO THE NEXT ION NI
C=======================================================================
 5300 CONTINUE
C=======================================================================
C USE INVERSION THROUGH A POINT MIDWAY BETWEEN ION NI AND NNI TO OBTAIN
C THE FORCES ON ION NNI DUE TO INTERACTION WITH THE PERIODIC ARRAY OF
C IONS NNI. ITS BASICALLY NEWTON'S THIRD LAW BUT SLIGHTLY SUBTLER
C=======================================================================
C     DO 5400 NI=1,NIONST      << was
C
      DO 5400 NI=2,NIONST
      DO 5401 NNI=1,NI-1
      DO 5402 M=1,3
      CFORCE(M,NNI,NI)=-CFORCE(M,NI,NNI)
      RFORCE(M,NNI,NI)=-RFORCE(M,NI,NNI)
 5402 CONTINUE
 5401 CONTINUE
 5400 CONTINUE
C=======================================================================
C CALCULATE THE SCALING FACTORS THAT WILL BE USED TO CONVERT THE VALUES
C OF THE FORCES AND THE ENERGY CALCULATED ABOVE TO REAL FORCES AND
C ENERGIES. THE FACTORS TAKE ACCOUNT OF THE NUMERICAL CONSTANTS
C INTRODUCED BY THE EWALD SEPARATION OF THE COULOMB LATTICE SUMS, THE
C RESCALING OF ALL OF THE EWALD INTEGRALS ONTO THE RANGE ZERO TO ONE
C AS WELL PUTTING IN THE ELECTROSTATIC CONSTANTS VIA 'ELECT'
C=======================================================================
      IF(IPRINT.GE.5)WRITE(*,*)'VOLC,SCALE:', VOLC,SCALE
      GENSCA=2.0*PI*ELECT/(VOLC*SCALE**2)
      GSCALE=((2.0*PI)**2)*ELECT/(VOLC*SCALE**2)
C=======================================================================
C CALCULATE THE TOTAL FORCES ON THE IONS
C=======================================================================
      DO 5410 NI=1,NIONST
      DO 5411 NNI=1,NIONST
      ICHSQ=ICHARC(NI)*ICHARC(NNI)
      DO 5412 M=1,3
      EWIFC(M,NI)=EWIFC(M,NI)
     &+( RFORCE(M,NNI,NI)*RSCALE+AIMAG(CFORCE(M,NNI,NI))*GSCALE )*ICHSQ
 5412 CONTINUE
 5411 CONTINUE
 5410 CONTINUE
C=======================================================================
C CALCULATE THE TOTAL STRESS ON THE UNIT CELL
C=======================================================================
      DUM=NIONCH**2*GENSCA/4.0
      SIGEW(1)=(RSIFOR(1)*RSCALE+REAL(CSIFOR(1))+DUM)/VOLC
      SIGEW(2)=(RSIFOR(2)*RSCALE+REAL(CSIFOR(2))+DUM)/VOLC
      SIGEW(3)=(RSIFOR(3)*RSCALE+REAL(CSIFOR(3))+DUM)/VOLC
      SIGEW(4)=(RSIFOR(4)*RSCALE+REAL(CSIFOR(4))    )/VOLC
      SIGEW(5)=(RSIFOR(5)*RSCALE+REAL(CSIFOR(5))    )/VOLC
      SIGEW(6)=(RSIFOR(6)*RSCALE+REAL(CSIFOR(6))    )/VOLC
C=======================================================================
C CALCULATE THE TOTAL COULOMB ENERGY PER CELL INCLUDING THE ENERGY DUE
C THE INTERACTION BETWEEN THE IONS AND THE UNIFORM NEUTRALISING
C BACKGROUND (GIVEN BY THE FINAL TWO TERMS IN THE SUM BELOW)
C=======================================================================
      TEWEN=RENSCA*REWEN+REAL(CEWEN)*GENSCA-(NIONCH**2*GENSCA/4.0)
     &   -(SCALE*ELECT*NICHSQ/SQRT(PI))
      WRITE(*,*)' '
      RETURN
      END
C
C
C
C
C
C=======================================================================
C     CHECK TO SEE WHETHER THE PARAMETERS MAXCX,..., ARE BIG ENOUGH
C     SUGGESTION WILLBE MADE.		             X.WENG 24-OCT-1990
C=======================================================================
      SUBROUTINE MCHECK(DIRC,RECC,ARGMRE,ARGMRC,MAXCX,MAXCY,MAXCZ,
     &            MAXGPX,MAXGPY,MAXGPZ,SCALE,PIDSCA,IPRINT)
      DIMENSION DIRC(3,*), RECC(3,*)
C
      A1=SQRT(DIRC(1,1)**2+DIRC(1,2)**2+DIRC(1,3)**2)
      A2=SQRT(DIRC(2,1)**2+DIRC(2,2)**2+DIRC(2,3)**2)
      A3=SQRT(DIRC(3,1)**2+DIRC(3,2)**2+DIRC(3,3)**2)
C
      B1=SQRT(RECC(1,1)**2+RECC(1,2)**2+RECC(1,3)**2)
      B2=SQRT(RECC(2,1)**2+RECC(2,2)**2+RECC(2,3)**2)
      B3=SQRT(RECC(3,1)**2+RECC(3,2)**2+RECC(3,3)**2)
C
      IF(IPRINT.GT.1) THEN
      WRITE(*,5)'THE LENGTH OF UNIT  CELL VECTORS (A):', A1,A2,A3
      WRITE(*,5)'THE LENGTH OF RECIPORCAL VECTORS (A):', B1,B2,B3
    5 FORMAT(1X,A37,1X,3F10.5)
      END IF
C=======================================================================
C     CHECK TO SEE IF MAXCX,MAXCY, MAXCZ ARE LARGE ENOUGH  23-OCT-1990
C=======================================================================
      A=SQRT( ( ARGMRE/SCALE)**2 
     &      + ( ABS(DIRC(1,1))+ABS(DIRC(2,1))+ABS(DIRC(3,1)) )**2
     &      + ( ABS(DIRC(1,2))+ABS(DIRC(2,2))+ABS(DIRC(3,2)) )**2
     &      + ( ABS(DIRC(1,3))+ABS(DIRC(2,3))+ABS(DIRC(3,3)) )**2 )
      MRX=A/A1+1
      MRY=A/A2+1
      MRZ=A/A3+1
      IF( MAXCX.GE.MRX .AND. MAXCY.GE.MRY .AND. MAXCZ.GE.MRZ) THEN
      IF(IPRINT.GT.1) THEN
      WRITE(*,*)'THE CURRENT REAL SPACE SUMMATION IS LARGE ENOUGH'
      WRITE(*,10) MRX, MRY, MRZ
      END IF
      ELSE
      WRITE(*,*)'THE CURRENT REAL SPACE SUMMATION IS TOO SMALL <****'
      WRITE(*,10) MRX, MRY, MRZ
      END IF
   10 FORMAT(' THE MINIMUM EWALD  R E A L   SPACE SUMMATION:',3I4)
C=======================================================================
C     DO THE SAME FOR MAXGPX,MAXGPY, MAXGPZ 
C=======================================================================
      B=SQRT( ( ARGMRC/PIDSCA*(2.0*3.1415926545897)          )**2
     &      + ( ABS(RECC(1,1))+ABS(RECC(2,1))+ABS(RECC(3,1)) )**2
     &      + ( ABS(RECC(1,2))+ABS(RECC(2,2))+ABS(RECC(3,2)) )**2
     &      + ( ABS(RECC(1,3))+ABS(RECC(2,3))+ABS(RECC(3,3)) )**2 )
      MGX=B/B1+1
      MGY=B/B2+1
      MGZ=B/B3+1
      IF( MAXGPX.GE.MGX .AND. MAXGPY.GE.MGY .AND. MAXGPZ.GE.MGZ) THEN
      IF(IPRINT.GT.1) THEN
      WRITE(*,*)'THE CURRENT RECIPROCAL SUMMATION IS LARGE ENOUGH'
      WRITE(*,20) MGX, MGY, MGZ
      END IF
      ELSE
      WRITE(*,*)'THE CURRENT RECIPROCAL SUMMATION IS TOO SMALL <****'
      WRITE(*,20) MGX, MGY, MGZ
      END IF
   20 FORMAT(' THE MINIMUM EWALD RECIPROCAL SPACE SUMMATION:',3I4)
      RETURN
      END
