      PROGRAM CASTEP
       INCLUDE 'param.inc'
C========================================================================
C   
C                                CASTEP 
C
C                 CAmbridge Serial Total Energy Package
C                             Cambridge 1991
C========================================================================
C
C     THIS IS A MODIFIED VERSION OF MC PAYNE'S CONJUGATE GRADIENT PROGRAM
C     EFFORTS OF FOLLOWING COLLABORATORS WHO TOOK PART IN THE 
C     DEVELOPMENT OF THIS PACKAGE SHOULD BE MENTIONED :
C     LATE XUDONG WENG, VICTOR MILMAN, GRAHAM FRANCIS, UWE BERTRAM,
C     BJORK HAMMER ET AL.
C
C     CURRENT MAINTENANCE AND SUBSEQUENT CHANGES ARE DUE TO VICTOR MILMAN
C     E-MAIL   VYM10@UK.AC.CAM.PHX     OR
C              MILMAN@UK.AC.CAMBRIDGE.WEST.PHYSICS
C
C========================================================================
C
C                               CGION 8
C
C========================================================================
C
C     THE MAJOR CHANGES ARE:
C     1) EITHER LOCAL OR NONLOCAL PSEUDOPOTENTIALS COULD BE USED, WITH 
C        THE LATER IN THE KLEIMANN-BYLANDER FORM.
C     2) GAUSSIAN BROADENING OF BANDS ARE INCLUDED TO ACCOUNT POSSIBLE
C        PARTIAL FILLING OF BANDS AROUND THE FERMI ENERGY.
C     3) THIS VERSION ALLOWS TO PERFORM EFFICIENT STATIC IONIC RELAXATION
C     4) SYMMETRIZATION OF THE CHARGE DENSITY IS INTRODUCED THUS
C        IT IS POSSIBLE TO USE ONLY SPECIAL K-POINTS WHICH ARE
C        INEQUIALENT WITH RESPECT TO SYMMETRY OPERATIONS
C     5) CONSTRAINED IONIC RELAXATION IS ALLOWED WITH ARBITRARY
C        CONSTRAINTS SUPPLIED BY USER (IN CASE OF SYMMETRIZED 
C        VERSION SYMMETRY IS PRESERVED AUTOMATICALLY)
C      
C-------------------------------------------------------------------------
C     THE CHANGES ARE LISTED HERE FOR YOUR REFERENCE.  
C
C     16-MAR-90 COPYIED FROM GEGBSC.FOR (WENG'S LATEST VERSION OF PAYNE'S
C               LOCAL PSEUDOPTENTIAL CONJUGATE GRADIENT PROGRAM)
C     16-MAR-90 THE FOLLOWING CHANGES ARE MADE
C		-- DIMENSION STATEMENTS ADDED FOR NON-LOCAL CALCULATION
C	        -- NEW READ IN STATEMENTS, FROM UNIT 11
C		-- DNLG, DNLKG ADDED WHEN CALLING GENSP1
C		-- CALL FNLFOR ADDED (FROM 'DO 3131' TO '3130 CONTINUE')
C               -- CALL NLCV ADDED TO SET UP CVL(G,Rn)=V(G)EXP(IG.Rn)
C               -- CALL VNL ADDED TO CALCULATE <G'|Vnl|G>
C
C     31-MAR-90 TESTED ON SILICON (OK).   
C
C     14-APR-90 VNL ADDED, REPLACING THE AUX IN CONGRA
C     20-APR-90 VNL INITIALISED, NECESSARY FOR NKPTS>1
C               DO 8364 ... 
C
C     23-MAY-90 Several modifations for re-starting
C               CELEN   ( DO 35 ....) 
C               GENSP1 - > GENSP
C
C     15-JUN-90 IVPTYP = 0 using local pseudopotential     -- tested and
C	               = 1 using non-local pseudopotential -- worked
C
C     17-JUN-90 CV0, CV1, AND CV2 ARE REDUCED TO ONE K-POINT A TIME
C               TO REDUCE THE DEMAND ON MEMORY SPACE.
C               ( SEE CONGRA FOR COMMENTS )
C
C     18-JUN-90 GAUSSIAN BROADENING OF BANDS ADDED, ACTIVATED BY SETTING
C               IOCCUP=1. TEST ON AL FCC 4 ATOMUNIT CELL AGREED WITH 
C               B. HAMMER'S INDEPENDENT RESULTS.
C
C     26-JUN-90 RESTRUCTURED THE NON-LOCAL PART. THE USAGE OF CV0,CV1 AND
C               CV2 IS ABANDONED. WE NOW USE
C 	        1) CPHSGR(NRPLWV,NIONS,NSPEC) = EXP(iq.Rn)
C	        2) VGNL(RPLWV,0:2,NNSPEC) = Vnl(l,q)
C	        IN ORDER TO REDUCE THE REQUIRED SPACE. TESTED OK
C
C     03-JUN-90 NON-ORTHOGONAL UNIT CELL PART OF HAMMER'S PROGRAM IS
C		INCORPORATED.
C
C     03-JUN-90 CPU TIME COUNTER ADDED. TO ACTIVATE, SET ICLOCK=1
C
C     03-JUN-90 NEW VERSION OF NON-LOCAL POTENTIAL (STILL FOR CUBIC 
C		UNIT CELL ONLY). 
C
C     22-JUL-90 RE-START WITH DIFFERENT CUT-OFF ENERGY IS SUCESSFUL.
C               HOWEVER, MUST WITH NGX,NGY,NGZ,NRPLWV UNCHANGED.
C
C
C     14-OCT-90 PARAMETER STATEMENT IN SUBROUTINE EWALD.FOR IS MOVED
C	        INTO THE MAIN PROGRAM. 
C
C     11-JAN-91 SEVERAL CHANGES ARE MADE, MAINLY RELATED TO THE ADDITION
C               OF NON-LOCAL CONTRIBUTION OF THE STRESS 
C               AND A LITTLE RE-STRUCTURING.
C
C     13-FEB-91 IONIC RELAXATION INCLUDED. THIS VERSION DOES NOT NEED
C               IMDTR ROUTINE. NEW PARAMETER NIONCG SHOULD BE SPECIFIED
C               IN A COMMAND FILE. TESTED FOR SILICON.
C
C     26-APR-91 NONLOCAL POTENTIAL PART IS RESTRUCTURED, AND
C               ORDER OF INDICES IS CHANGED FOR CORRESPONDING ARRAYS.
C               AS A RESULT THE NUMBER OF FORMAL PARAMETERS FOR THE
C               CALL OF VNLWAV HAS BEEN CHANGED.
C
C               NITFIX VARIABLE INTRODUCED TO PERFORM NITFIX ITERATIONS
C               OF CONGRA BEFORE MOVING IONS.
C
C      1-MAY-91 RANDOM GENERATION OF INITIAL WAVE FUNCTION IS INCLUDED
C               FOR INRAND=1 (NEW PARAMETER) AND IS PROVIDED BY
C               WFINRN SUBROUTINE
C
C      1-JUL-91 UNIT NUMBERS FOR FILES WITH WAVEFUNCTIONS ARE CHANGED
C
C     16-JUL-91 NEW SCHEME FOR STATIC IONIC RELAXATION WITH THE ENERGY
C               MINIMISATION IN BOTH IONIC AND ELECTRONIC COORDINATES 
C               SPACE IS INTRODUCED. TESTED FOR SEMICONDUCTORS.
C
C     15-SEP-91 ADDITIONAL CALLS TO SUBROT AND EFERMI ARE INSERTED IN
C               THE IONIC RELAXATION PART TO ACCELERATE CONVERGENCE
C               FOR METALS (IOCCUP=1). TESTED FOR Al
C
C     10-OCT-91 MIXED LOCAL-NONLOCAL POTENTIALS ARE ENABLED. READING
C               FROM FILE 11 HAS BEEN CHANGED FOR LOCAL POTENTIAL AND
C               ARRAY DPSP HAS BEEN REMOVED.
C
C      9-NOV-91 BUG IN THE MAIN CODE IS REMOVED WHICH WAS AFFECTING
C               SIMULATIONS WITH ODD NUMBER OF ELECTRONS (LOOP 5210).
C
C               VARIABLE NPRINT HAS NOW MEANING OF THE FREQUENCY OF
C               FORCE AND STRESS CALCULATION WHICH IS PERFORMED EACH
C               NPRINT'TH ITERATION BEFORE NDELAY. AFTER NDELAY
C               IT HAS NO MEANING. DURING RELAXATION, STRESS IS
C               CALCULATED ONLY IF UNIT CELL IS SUPPOSED TO RELAX,
C               OTHERWISE IT IS CALCULATED AT THE LAST ITERATION ONLY.
C
C     10-NOV-91 CONSTRAINED IONIC RELAXATION IS ENABLED. USER DEFINED
C               SUBROUTINE "CONSTR" IS REQUIRED WHICH SPECIFIES MATRIX
C               OF DERIVATIVES OF OLD COORDINATES WITH RESPECT TO NEW
C               ONES. SAMPLE "CONSTR" CORRESPONDS TO ELIMINATION OF
C               FIXED ATOMS AND/OR TO PHONON DISPLACEMENTS.
C               
C     13-JAN-92 SYMMETRIZED SCHEME IS ENABLED. ONLY SYMMETRY
C               UNRELATED SPECIAL K-POINTS NEED TO BE USED, AND
C               PROPER CHARGE DENSITY, IONIC FORCES AND STRESSES
C               ARE CALCULATED USING SYMMETRY OPERATIONS.
C               EXTRA DATA FILE (UNIT 3) IS REQUIRED WITH
C               CRYSTALLOGRAPHIC INFORMATION. THIS FILE SHOULD
C               BE CREATED BY K290 PROGRAM (MODIFIED PROGRAM
C               BY R.J.NEEDS,K.KUNC,O.H.NIELSEN,R.M.MARTIN)
C
C***********************************************************************
C***********************************************************************
C
C  GEGBS
C
C***********************************************************************
C***********************************************************************
C=======================================================================
C THIS IS THE MAIN PROGRAM OF A SET OF PROGRAMS WRITTEN TO PERFORM TOTAL
C ENERGY CALCULATIONS USING THE MOLECULAR DYNAMICS METHOD INTRODUCED BY
C R. CAR AND M. PARRINELLO (PHYSICAL REVIEW LETTERS 55, 2471 (1985)).
C=======================================================================
C THIS VERSION OF THE PROGRAMS WAS WRITTEN FOR USE ON THE FPS-164
C COMPUTER, CALLS TO FPS SUBROUTINES WILL BE IDENTIFIED WHEREVER THEY
C APPEAR (WHICH IS ONLY IN THE MAIN PROGRAM AND THE ORTHOGONALISATION
C SUBROUTINE). THE FPS-164 IS A 64 BIT MACHINE SO THESE SINGLE PRECISION
C PROGRAMS ACTUALLY RUN AT THE ACCURACY OF FORTRAN DOUBLE PRECISION.
C FOR USE ON NON-64 BIT MACHINES THE VARIABLES SHOULD BE DECLARED DOUBLE
C PRECISION OR THE PROGRAMS COMPILED USING AN AUTO-DOUBLE COMPILER.
C=======================================================================
C THIS IS THE SIMPLEST FORM OF THE PROGRAMS. THE UNIT CELL MUST BE
C ORTHORHOMBIC AND THE WAVEFUNCTIONS MUST ALL FIT INTO CORE MEMORY. THE
C METHOD IS MEMORY INTENSIVE AND EXPERIENCE HAS SHOWN THAT CALCULATING
C THE REAL SPACE WAVEFUNCTIONS TWICE EACH TIMESTEP IS FASTER THAN
C WRITING AND READING THEM FROM DISK. IF THE COMPUTER HAS A SOLID STATE
C DISK THIS WILL PROBABLY NO LONGER BE TRUE. THE PROGRAM USES A LOCAL
C PSEUDOPOTENTIAL TO EXPLOIT THE INCREASE IN SPEED OBTAINED BY USING
C FAST FOURIER TRANSFORMS. THE PROGRAM IS WRITTEN FOR A SINGLE SPECIES
C OF ION. PROGRAMS FOR NON-ORTHORHOMBIC UNIT CELLS, FOR MORE THAN ONE
C SPECIES OF ION AND FOR HANDLING MORE WAVEFUNCTIONS THAN FIT IN CORE
C MEMORY HAVE BEEN WRITTEN AND ARE SCATTERED BETWEEN OTHER MEMBERS OF
C PROFESSOR JOHN JOANNOPOULOS'S RESEARCH GROUP AT THE MASSACHUSETTS
C INSTITUTE OF TECHNOLOGY
C=======================================================================
C THIS PROGRAM IS SET UP TO RELAX THE ELECTRONS TO THEIR GROUNDSTATE
C WHILE KEEPING THE IONS AND THE SIZE OF THE  UNIT CELL FIXED. THE
C ELECTRONIC SYSTEM MUST BE RELAXED TO ITS GROUNDSTATE BEFORE THE FORCES
C EXERTED BY THE ELECTRONS ON THE IONS AND THE UNIT CELL ARE CORRECTLY
C GIVEN BY THE HELLMANN-FEYNMAN THEOREM. THE PROGRAM FUFGMSP IS SET UP
C TO READ IN THE RELAXED WAVEFUNCTIONS AND PROCEED TO PERFORM MOLECULAR
C DYNAMICS OR SIMULATED ANNEALING ON THE IONIC SYSTEM OR TO RELAX THE
C SIZE OF THE UNIT CELL. THESE OPERATIONS MAY BE PERFORMED USING THE
C PRESENT PROGRAM BY SIMPLY REMOVING THE FORTRAN COMMENTS FROM THE
C APPROPRIATE SUBROUTINE CALLS
C=======================================================================
C THE UNITS USED IN THE PROGRAMS ARE ELECTRON-VOLTS AND ANGSTROMS
C=======================================================================
C THE STRUCTURE OF THE PROGRAM IS ALMOST TOTALLY DETERMINED BY THE USE
C OF FAST FOURIER TRANSFORMS. THE INDEXING OF THE RECIPROCAL LATTICE IS
C IMPOSED BY THE FAST FOURIER TRANSFORM. THE FFT INTRODUCES A FICTITIOUS
C PERIODICITY IN RECIPROCAL SPACE, IN THE USUAL LABELLING THE RECIPROCAL
C LATTICE POINTS WITH NEGATIVE RECIPROCAL LATTICE VECTOR ARE SHIFTED BY
C THE WAVEVECTOR OF THE ARTIFICIAL PERIODICITY TO BECOME THE COMPONENTS
C LABELLED BY LARGE POSITIVE RECIPROCAL LATTICE VECTORS. THIS EXPLAINS
C THE NEED FOR THE LOOP COUNTERS LPCTX,Y,Z WHICH KEEP TRACK OF THE
C TRUE RECIPROCAL LATTICE VECTORS REPRESENTED BY THE RECIPROCAL LATTICE
C POINTS. THE CALL TO THE FAST FOURIER TRANSFORM IS THROUGH FFT3D
C WHICH IS A SUBROUTINE THAT REORDERS THE DATA IN THE THREE DIMENSIONAL
C FOURIER TRANSFORM GRID TO MAKE OPTIMUM USE OF CFFT99 WHICH IS A
C VECTORISED FAST FOURIER TRANSFORM ROUTINE. PARAMETER STATMENTS IN
C MUST BE CHANGED BEFORE RUNNING THE PROGRAMS
C
C IN FFT3D: ISIGN = -1 : FORWARD : REAL --> RECIP : EXP(+IQR)
C           ISIGN = +1 : BACKWARD: RECIP --> REAL : EXP(-IQR)
C
C=======================================================================
C
C **************TERMS IN THE PARAMETER STATEMENT************************
C
C NGX,NGY,NGZ = THE NUMBER OF RECIPROCAL LATTICE POINTS (AND DUE TO THE
C      PROPERTIES OF THE FAST FOURIER TRANSFORM) THE NUMBER OF POINTS ON
C      THE REAL SPACE GRID IN THE X,Y AND Z DIRECTIONS, RESPECTIVELY.
C NIONS = THE NUMBER OF IONS IN THE UNIT CELL.
C NPSPTS = THE NUMBER OF VALUES IN THE ARRAY USED TO INTERPOLATE THE
C             PSEUDOPOTENTIAL.
C NBANDS = THE NUMBER OF VALENCE BANDS
C NKPTS = THE NUMBER OF SPECIAL K POINTS USED FOR AVERAGING OVER THE
C             BRILLOUIN ZONE.
C NEWPTS = THE NUMBER OF VALUES IN THE ARRAY USED TO INTERPOLATE
C             THE EWALD INTEGRALS
C NEDDPT = IS THE NUMBER OF VALUES USED TO INTERPOLATE THE DATA FOR THE
C             EQUATIONS OF MOTION OF THE ELECTRONIC STATES.
C NPLWV = THE TOTAL NUMBER OF POINTS IN THE REAL AND RECIPROCAL SPACE
C             GRIDS
C NRPLWV = THE MAXIMUM NUMBER OF PLANE WAVE COMPONENTS THAT CAN BE
C      INCLUDED AT ANY K-POINT (THIS WILL LIMIT THE CUT-OFF ENERGY USED)
C=======================================================================
C      IMPLICIT COMPLEX*8 (C)
C=======================================================================
C           Dimensional parameters and control parameters
C                 SI 8 atom unit cell
C
C      PARAMETER(NGX=8,NGY=8,NGZ=8,
C     &   NPLWV=NGX*NGY*NGZ,MPLWV=NGX*NGY*(NGZ+1),
C     &   NPSPTS=501,
C     &   NRPLWV=100,
C     &   NEWPTS=1000, NEDDPT=2000,
C     &   NIONS=8,    NSPEC=1,
C     &   NBANDS=16,   NBANOC=16,    NKPTS=1)
C      DATA NAT /8/
C      DATA ISBROT/0/
C      DATA IOCCUP,DELMIN,DELMAX,NDEL/0,0.2,0.8,3/
C      DATA IPRINT/2/
C      DATA ICLOCK/0/
C      DATA IION,IBOX/1,0/
C      DATA ISTART,NITMAX,NDELAY/0,2,15/
C      DATA NITFIX/1/
C      DATA INRAND/0/
C      DATA NMDCG/1/
C      DATA NIONCG/2/
C=======================================================================
C  THE FOLLOWING ARE PARAMETERS ORIGINALLY INSIDE SUBROUTINE EWLALTR.F
C=======================================================================
C      PARAMETER(MAXCX=3, MAXCY=3, MAXCZ=3,
C     &          MAXGPX=3,MAXGPY=3,MAXGPZ=3,
C     &   NCELLX=(2*MAXCX)+1,NCELLY=(2*MAXCY)+1,NCELLZ=(2*MAXCZ)+1,
C     &   NCELLS=NCELLX*NCELLY*NCELLZ,
C     &   NGPTX=(2*MAXGPX)+1,NGPTY=(2*MAXGPY)+1,NGPTZ=(2*MAXGPZ)+1,
C     &   NGPTS=NGPTX*NGPTY*NGPTZ)
C=====================================================================
C=======================================================================
C FOR AL 4 ATOM CELL
C=======================================================================
C      PARAMETER(NGX=16,NGY=16,NGZ=16,NIONS=4,NSPEC=1,NBANOC=12,
C     &   NPSPTS=600,NBANDS=12,NKPTS=1,NEWPTS=1000,NEDDPT=2000,
C     &   NPLWV=NGX*NGY*NGZ,MPLWV=NGX*NGY*(NGZ+1),NRPLWV=250)
C=======================================================================
C Parameters are added:
C=======================================================================
C  NITMAX: ITERATION NUMBER IN CONGRA
C  NDELAY: ITERATION NUMBER DELAYED, BEFORE ANY DYNAMICS STARTS.
C  ISTART: 0 FOR STARTING FROM BEGIN
C          1 FOR RE-START, NEED FOR013 AND FOR016
C  IVPTYP: 0 FOR LOCAL PSEUDOPOTENTIAL
C          1 FOR NON-LOCAL PSEUDOPOTENTIAL
C  ISBROT: SUB-SPACE ROTATION, 0=NO, 1=YES  (SHOULD BE 1, IF IOCCUP=1)
C  IOCCUP: USING OCCUPATION NUMBER 0=NO, 1=YES
C  IPRINT: THE EXTENT OF PRINT-OUT, 1 FOR GENERAL USE 
C=======================================================================
C      DATA ISTART,NITMAX,NDELAY/0,3,10/
C      DATA IVPTYP/0/
C      DATA ISBROT/1/
C      DATA IOCCUP,DELMIN,DELMAX,NDEL/1,0.026,0.2,4/
C      DATA IPRINT/1/
C      DATA ICLOCK/0/
C=======================================================================
C  THE FOLLOWING ARE PARAMETERS ORIGINALLY INSIDE SUBROUTINE EWLALTR.F
C							14-OCT-90
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,
C     &          MAXGPX=4,MAXGPY=4,MAXGPZ=6,
C     &   NCELLX=(2*MAXCX)+1,NCELLY=(2*MAXCY)+1,NCELLZ=(2*MAXCZ)+1,
C     &   NCELLS=NCELLX*NCELLY*NCELLZ,
C     &   NGPTX=(2*MAXGPX)+1,NGPTY=(2*MAXGPY)+1,NGPTZ=(2*MAXGPZ)+1,
C     &   NGPTS=NGPTX*NGPTY*NGPTZ)
C=======================================================================
C
C  * VARIABLES AND DIMENSION STATEMENTS FOR ARRAYS USED IN SUBROUTINES *
C
C SIZEX,SIZEY,SIZEZ = THE SIZE OF THE UNIT CELL IN X,Y,Z DIRECTIONS
C SIZELX,SIZELY,SIZELZ = THE SIZE OF THE UNIT CELL AT THE LAST TIMESTEP
C SIZEIX,SIZEIY,SIZEIZ = THE INITIAL SIZE OF THE UNIT CELL, THESE MUST
C BE USED IN THE CALL TO GENSP TO ENSURE THAT THE INDICES FOR THE
C COMPONENTS OF THE WAVEFUNCTIONS ARE CORRECT. (IF A CONSTANT CUT-OFF
C ENERGY FOR THE PLANE WAVE BASIS STATES IS REQUIRED THE WAVEFUNCTIONS
C CAN BE READ IN USING THE INDEXING SCHEME BASED ON SIZEIX,Y,Z AND THEN
C           PUT INTO A NEW ARRAY BASED ON THE PRESENT SIZE OF THE BOX).
C POSION(3,NIONS) = THE COORDINATES OF THE IONS IN THE UNIT CELL
C POSIOL(3,NIONS) = THE COORDINATES OF THE IONS AT THE LAST TIMESTEP
C CSTRF(NPLWV) = THE STRUCTURE FACTOR AT THE GRID OF RECIPROCAL LATTICE
C                POINTS
C PSP(NPSPTS) = THE ARRAY OF DATA VALUES USED TO INTERPOLATE THE
C                    PSEUDOPOTENTIAL
C PSCORE = THE G=0 COMPONENT OF THE PSEUDOPOTENTIAL REMAINING AFTER
C          SUBTRACTION OF THE COULOMB POTENTIAL. THIS CONTRIBUTES TO
C                    THE AVERAGE POTENTIAL IN THE PSEUDO-SOLID
C PSGMAX = THE G VECTOR CORRESPONDING TO THE LAST DATA VALUE IN THE
C                    PSEUDOPOTENTIAL ARRAYS
C SIGAL(6) = THE STRESS ON THE UNIT CELL DUE TO THE PSEUDOPTENTIAL CORE
C PSCENC = THE CONTRIBUTION TO THE TOTAL ENERGY DUE TO THE
C                    PSEUDOPOTENTIAL CORE
C VPS(NPLWV) = THE PSEUDOPOTENTIAL ON THE GRID OF RECIPROCAL LATTICE
C                    VECTORS
C DVPS(NPLWV) = THE DERIVATIVE OF THE PSEUDOPOTENTIAL ON THE GRID OF
C                    RECIPROCAL LATTICE VECTORS
C CV(NPLWV) = THE TOTAL POTENTIAL ON THE GRID OF RECIPROCAL LATTICE
C             VECTORS (IE PSEUDOPOTENTIAL*STRUCTURE FACTOR + HARTREE +
C                    EXCHANGE-CORRELATION POTENTIAL)
C CVD(NPLWV) = THE HARTREE POTENTIAL ON THE GRID OF RECIPROCAL LATTICE
C                    VECTORS
C CHDENR(NPLWV) = THE CHARGE DENSITY ON THE GRID OF REAL LATTICE VECTORS
C CHDENG(NPLWV) = THE CHARGE DENSITY ON THE GRID OF RECIPROCAL LATTICE
C                    VECTORS
C DSIF(3) = THE FORCE ON THE UNIT CELL DUE TO THE HARTREE ENERGY
C DENC = THE CORRECTION TO THE TOTAL ENERGY DUE TO OVERCOUNTING THE
C        HARTREE ENERGY ON SUMMING THE ENERGY EIGENVALUES OF THE
C                     ELECTRONIC STATES
C XCSIF(3) = THE FORCE ON THE BOX DUE TO THE EXCHANGE CORRELATION ENERGY
C XCENC = THE CORRECTION TO THE TOTAL ENERGY DUE TO OVERCOUNTING THE
C         EXCHANGE CORRELATION  ENERGY ON SUMMING THE ENERGY EIGENVALUES
C                     OF THE ELECTRONIC STATES
C VKPT(3,NKPTS) = THE COORDINATES OF THE SPECIAL K POINTS USED FOR
C                     BRILLOUIN ZONE AVERAGES
C WTKPT(NKPTS) = THE WEIGHT TO BE ATTACHED TO EACH SPECIAL K POINT
C CPTWFP(NRPLWV,NBANDS,NKPTS) = THE COMPONENTS OF THE ELECTRONIC
C                     WAVEFUNCTIONS AT THE PRESENT TIMESTEP
C CPTWFL(NRPLWV,NBANDS,NKPTS) = THE COMPONENTS OF THE ELECTRONIC
C                     WAVEFUNCTIONS AT THE LAST TIMESTEP
C SIKEF(3) = THE FORCE ON THE UNIT CELL DUE TO THE KINETIC ENERGY OF THE
C                     ELECTRONIC STATES
C TIME = THE TIMESTEP USED IN INTEGRATING THE EQUATIONS OF MOTION FOR
C                     THE ELECTRONIC STATES
C ELDAMP = THE DAMPING USED WHEN INTEGRATING THE EQUATIONS OF MOTION FOR
C                     THE ELECTRONIC STATES
C SIDAMP = THE DAMPIMG FACTOR USED IN INTEGRATING THE EQUATIONS OF
C                     MOTION FOR THE SIZE OF THE UNIT CELL
C SIMASS = THE MASS ASSOCIATED WITH THE SIZE OF THE UNIT CELL
C PODAMP = THE DAMPING FACTOR USED IN INTEGRATING THE EQUATIONS OF
C                     MOTION FOR THE POSITIONS OF THE IONS
C POMASS = THE MASS ASSOCIATED WITH THE IONS
C CELEN(NBANDS,NKPTS) = THE ENERGY EIGENVALUES OF THE ELECTRONIC STATES
C EIFOR(3,NIONS) = THE FORCE EXERTED ON THE IONS BY THE ELECTRONS
C EISIF(3) = THE FORCE EXERTED ON THE UNIT CELL DUE THE ELECTRON-ION
C                     ENERGY
C EWRLEN(NEWPTS),EWRCSS(NEWPTS),EWRLSS(NEWPTS),EWRCSI(NEWPTS) = THE DATA
C           ARRAYS USED TO INTERPOLATE THE VALUES OF THE EWALD INTEGRALS
C TEWEN = THE EWALD ENERGY
C EWIFOR(3,NIONS) = THE FORCE ON THE IONS DUE TO THE COULOMB FORCES
C                     FROM THE OTHER IONS
C EWSIF(3) = THE FORCE ON THE UNIT CELL DUE TO THE EWALD ENERGY
C LPCTX(NGX),LPCTY(NGY),LPCTZ(NGZ) = ARRAYS USED TO LABEL THE RECIPROCAL
C           LATTICE VECTORS. DUE TO THE PERIODICITY IMPOSED BY THE FAST
C           FOURIER TRANSFORMS THE NEGATIVE RECIPROCAL LATTICE VECTORS
C           GX=-1,-2..,N,..(2*PI/SIZEX) ARE CONTAINED IN THE NGX+1-N
C           COMPONENTS OF THE ARRAY TO BE TRANSFORMED. THESE ARRAYS GIVE
C           THE CORRECT LABELLING OF THE RECIPROCAL LATTICE VECTORS
C LPCTFX(NGX),LPCTFY(NGY),LPCTFZ(NGZ) = AS ABOVE BUT WITH THE (NGX/2)+1
C           COMPONENTS SET EQUAL TO 0 SO THAT THE ABSENCE OF A
C           -((NGX/2)+1) WAVEVECTOR COMPONENT IN THE CHARGE DENSITY
C           DOESN'T LEAD TO A SPURIOUS FORCE ON THE IONS FROM THE
C           ELECTRONS
C NGPTAR(3) = THE NUMBER OF RECIPROCAL LATTICE POINTS IN THE X,Y,Z
C             DIRECTIONS, USED IN THE CALL TO THE FAST FOURIER TRANSFORM
C NDIM = THE DIMENSION OF THE FAST FOURIER TRANSFORM
C POTIM = A VARIABLE TIMEPERIOD USED IN INTEGRATING THE EQUATIONS OF
C         MOTION OF THE IONS
C PODISP = THE LARGEST DISTANCE (AS A FRACTION OF THE UNIT CELL SIZE)
C          THAT AN ION IS ALLOWED TO MOVE DURING ONE TIMESTEP. THIS
C          MUST BE LIMITED TO A VALUE THAT ALLOWS THE ELECTRONS TO
C          FOLLOW THE IONIC MOTIONS WITHOUT GETTING TOO FAR FROM THEIR
C          GROUNDSTATE OR THE HELLMANN-FEYNMAN FORCES WILL BE INCORRECT
C ENMAX = THE KINETIC ENERGY CUTOFF FOR THE PLANE WAVE BASIS SET
C NINDPW(NRPLWV,NKPTS) = THE NUMBER INDICES IN THE RECIPROCAL LATTICE
C           GRID OF THE PLANE WAVE BASIS STATES. THIS ARRAY IS USED TO
C           PLACE THE "SPHERE" OF PLANE WAVE BASIS STATES AT EACH K
C           POINT INTO THE BOX USED FOR THE FOURIER TRANSFORMS
C NPLWKP(NKPTS) = THE NUMBER OF PLANE WAVE BASIS STATES AT EACH K POINT
C DATAKE(7,NRPLWV,NKPTS) = THE KINETIC OF THE PLANE WAVE BASIS STATES
C           AND THE X,Y AND Z COMPONENTS OF THE KINETIC ENERGY AT EACH
C           K POINT
C RMOVE(NIONS*NSPEC) = 1.0 IF ION IS ALLOWED TO MOVE OTHERWISE = 0.0
C STATZD = DISTANCE FROM THE EDGE OF THE UNIT CELL TO THE FIXED
C           LAYERS OF ATOMS
C=======================================================================
      PARAMETER (NINSP=NIONS*NSPEC)  
      DIMENSION POSION(3,NIONS,NSPEC)
      DIMENSION POSIOL(3,NIONS,NSPEC)
      DIMENSION CSTRF(NPLWV,NSPEC)
      DIMENSION PSP(NPSPTS,NSPEC)
      DIMENSION VPS(NPLWV,NSPEC), DVPS(NPLWV,NSPEC)
      DIMENSION CV(NPLWV)
      DIMENSION CVD(NPLWV)
      DIMENSION CHDENR(NPLWV), CHDENG(NPLWV)
      DIMENSION CVION(NPLWV)
      DIMENSION VKPT(3,NKPTS), WTKPT(NKPTS)
      DIMENSION CPTWFP(NRPLWV,NBANDS)
      DIMENSION CPTWFL(NRPLWV)
      DIMENSION SIKEF(3)
      DIMENSION CELEN(NBANDS,NKPTS)
      DIMENSION EIFOR(3,NIONS,NSPEC)
      DIMENSION EWRLEN(NEWPTS),EWRCSS(NEWPTS),EWRLSS(NEWPTS),
     &  EWRCSI(NEWPTS)
      DIMENSION EWIFOR(3,NIONS,NSPEC)
      DIMENSION LPCTX(NGX),LPCTY(NGY),LPCTZ(NGZ)
      DIMENSION LPCTFX(NGX),LPCTFY(NGY),LPCTFZ(NGZ)
      DIMENSION NGPTAR(3)
      DIMENSION NINDPW(NRPLWV,NKPTS)
      DIMENSION NPLWKP(NKPTS)
      DIMENSION DATAKE(7,NRPLWV,NKPTS)
      DIMENSION PSCORE(NSPEC)
      DIMENSION PSGMAX(NSPEC)
      DIMENSION ICHARG(NSPEC)
      DIMENSION NIONSP(NSPEC)
      DIMENSION POMASS(NSPEC)
      DIMENSION POSIC(3,NINSP)
      DIMENSION POSICL(3,NINSP)
      DIMENSION EWIFC(3,NINSP)
      DIMENSION EIFC(3,NINSP)
      DIMENSION POMASC(NINSP)
      DIMENSION ICHARC(NINSP)
      DIMENSION RMOVE(NINSP), RMOVEO(NINSP)
      DIMENSION IVPTYN(NSPEC),NPSPTN(NSPEC)
      DIMENSION EXCDAT(2000),XCFDAT(2000),XCPDAT(2000)
      DIMENSION DIRDAT(NPLWV)
C
C
C********************* DIMENSION STATEMENTS ****************************
C
C TIFOR(3,NIONS) = THE TOTAL FORCE ON THE IONS
C
      DIMENSION TIFOR(3,NIONS,NSPEC)
      DIMENSION TIFORC(3,NINSP)
C
C
C     DIMENSION STATEMENTS FOR WORK ARRAYS USED IN SUBROUTINES
C
C     (the dimension of NWORK1 and WORK1 is tripled for random
C     generating of initial wave functions, X.WENG) 
C
      DIMENSION NWORK1(3*NBANDS)
      DIMENSION WORK1(3*NBANDS),WORK3(3,NKPTS),WORK7(3,NIONS,NSPEC),
     &    WORK8(3,NIONS,NSPEC),WORK9(3,NIONS,NSPEC)
      DIMENSION CWORK(MPLWV),CWORK1(MPLWV),CWORK2(MPLWV)
      DIMENSION CWORK3(NGX),CWORK4(NGY),CWORK5(NGZ)
      DIMENSION CWORK6(NRPLWV),CWORK7(NRPLWV),CWORK8(NRPLWV),
     &    CWORK9(NRPLWV),CWOR10(NRPLWV),CWOR11(NRPLWV),WORK12(NPLWV)
C=======================================================================
C     SPACE USED IN REASSN                                WENG 23-JUL-90
C=======================================================================
      DIMENSION IASSGN(NRPLWV), RECO(3,3)
C=======================================================================
C     SPACE USED IN CONGRA (in SUBROT)                    WENG 27-FEB-90
C=======================================================================
      DIMENSION HR(NBANDS, NBANDS), HI(NBANDS, NBANDS), AUX(NBANDS),
     & FV1(NBANDS), FV2(NBANDS), FV3(NBANDS), CH0(NBANDS, NBANDS)
C=======================================================================
C     SPACE USED FOR OCCUPATION NUMBER                    WENG 12-MAR-90
C=======================================================================
      REAL EIGVAL(NBANDS,NKPTS),OCC(NBANDS,NKPTS),SORT(NBANDS,NKPTS)
C=======================================================================
C     DIMENSION STATEMENTS FOR THE NON-LOCAL CALCULATION
C=======================================================================
      DIMENSION PSCALE(0:2,NSPEC),PSPNL(NPSPTS,0:2,NSPEC)
      DIMENSION DNLG(NRPLWV,3,NKPTS),DNLKG(NRPLWV,0:3,NKPTS)
      DIMENSION LOCAL(NSPEC)
C========================================================================
C ARRAYS NEEDED FOR REAL SPACE PROJECTION OF THE NON-LOCAL POTENTIALS
C========================================================================
C
C VRLNL(NRLPTS,MXRLNL,NSPEC) - THE NON-LOCAL PSEUDOPOTENTIALS ON A
C                              RADIAL GRID
C NRLNL(NSPEC) - THE NUMBER OF NON-LOCAL POTENTIALS FOR EACH SPECIES
C IRLNL(MXRLNL,NSPEC) - THE L NUMBER OF EACH NON-LOCAL POTENTIAL
C RLCORE(NSPEC) - THE CORE RADIUS FOR EACH SPECIES
C RMAX(NSPEC) - THE LARGEST RADIUS IN THE NON-LOCAL POTENTIAL ARRAYS
C VRLGRD(NRGRPT,MXRLSH,NIONST) - THE NON-LOCAL POTENTIALS ON THE RADIAL
C              SPACE GRID FOR EACH SPHERICAL HARMONIC FOR EACH POTENTIAL
C              FOR EACH SPECIES
C CPHGRD(NRGRPT,NIONST) - THE PHASE FACTORS FOR MULTIPLYING THE
C              WAVEFUNCTIONS ON THE REAL SPACE GRID
C DVRLGR(NRGRPT,3,MXRLSH,NIONST) - THE X.Y.Z. DERIVATIVES OF THE NON-
C              LOCAL POTENTIALS ON THE REAL SPACE GRID
C NSHSP(NSPEC) - THE TOTAL NUMBER OF NON-LOCAL PROJECTIONS (INCLUDING
C              'M' VALUES FOR EACH SPECIES
C NADGRD(NRGRPT,NIONST) - THE ADDRESS ARRAY THAT INDEXES POINTS ON THE
C              REAL SPACE GRID WITHIN THE CORE RADIUS OF EACH ATOM
C NRLPPI(NIONST) - THE TOTAL NUMEBR OF REAL SPACE POINTS WITHIN THE
C              CORE RADIUS AROUND EACH ATOM
C FNLRL(3,NIONST) - THE FORCES ON THE IONS
C=======================================================================
      DIMENSION  VRLNL(NRLPTS,MXRLNL,NSPEC)                        
      DIMENSION  NRLNL(NSPEC),PRLSCA(MXRLNL,NSPEC)                 
      DIMENSION  IRLNL(MXRLNL,NSPEC),RLCORE(NSPEC),RMAX(NSPEC)     
      DIMENSION  VRLGRD(NRGRPT,MXRLSH,NIONST),CPHGRD(NRGRPT,NIONST)
      DIMENSION  DVRLGR(NRGRPT,3,MXRLSH,NIONST),NSHSP(NSPEC)       
      DIMENSION  NADGRD(NRGRPT,NIONST),NRLPPI(NIONST)              
      DIMENSION  FNLRL(3,NIONST)                                   
C======================================================================
      DIMENSION CPHSGR(NRPLWV,NIONS,NSPEC),VGNL(NRPLWV,0:2,NSPEC),
     &                                     DVGNL(NRPLWV,0:2,NSPEC)
      DIMENSION CELFRC(NRPLWV)
      DIMENSION FNLEIF(3,NIONS,NSPEC),FORWK(3,NIONS)
      DIMENSION CWRK20(NIONS),CWRK21(3,NIONS)
      DIMENSION CWRK22(NIONS),CWRK23(3,3,NIONS)
      DIMENSION VNL(NBANDS,NKPTS)
C========================================================================
C WORK ARRAYS FOR REAL SPACE NON-LOCAL (PARALLEL)
C=======================================================================
      DIMENSION CESAVE(NIONST,36)
C=======================================================================
C     THE WORK ARRAYS ADDED FOR THE NON-LOCAL STRESS, BY U. BERTRAM
C=======================================================================
      DIMENSION CWRK30(NBANDS,NIONS),CWRK31(NBANDS,NIONS,6)
      DIMENSION CWRK32(NBANDS,NIONS,3),CWRK33(NBANDS,NIONS,6,3)
      DIMENSION CWRK34(NBANDS,NIONS,6,3),CWRK35(NBANDS,NIONS,3,3)
      DIMENSION CWRK36(NBANDS,NIONS,6,3,3),CWRK37(NBANDS,NIONS)
      DIMENSION CWRK38(NBANDS,NIONS,6),CWRK39(NBANDS,NIONS,6,3,3)
      DIMENSION SIGWK(6)
C========================================================================
C          ** DIMENSIONS ADDED FOR TRICLINIC CELL  **
C========================================================================
C
C DIRC(3,3) = CURRENT BASIS VECTORS IN REAL SPACE FOR THE SUPER-CELL
C DIRI(3,3) = INITIAL BASIS VECTORS IN REAL SPACE FOR THE SUPER-CELL
C DIRL(3,3) = LAST BASIS VECTORS IN REAL SPACE FOR THE SUPER-CELL
C RECC(3,3) = CURRENT BASIS VECTORS IN RECIPROCAL SPACE
C RECI(3,3) = INITIAL BASIS VECTORS IN RECIPROCAL SPACE
C VOLC = CURRENT VOLUME OF THE SUPER-CELL
C VOLI = INITIAL VOLUME OF THE SUPER-CELL
C        THE INITIAL-VALUES WILL BE USED IN THE CALL TO GENSP TO ENSURE
C        THAT THE INDICES FOR THE COMPONENTS OF THE WAVEFUNCTIONS ARE
C        CORRECT. (IF A CONSTANT CUT-OFF ENERGY FOR THE PLANE WAVE BASIS
C        STATES IS REQUIRED THE WAVEFUNCTIONS CAN BE READ IN USING THE
C        INDEXING SCHEME BASED ON DIRI() AND THEN PUT INTO A NEW ARRAY
C        BASED ON THE PRESENT BASUS OF THE BOX).
C SIGKE(6) = STRESS ON UNIT CELL DUE TO THE KIN.ENERGY OF THE ELEC.STATES
C SIGHA(6) = STRESS ON UNIT CELL DUE TO THE HARTREE ENERGY
C SIGXC(6) = STRESS ON UNIT CELL DUE TO THE EXCHANGE CORRELATION ENERGY
C SIGEI(6) = STRESS ON UNIT CELL DUE THE ELECTRON-ION ENERGY
C SIGAL(6) = STRESS ON UNIT CELL DUE TO THE PSEUDOPTENTIAL CORE
C SIGEW(6) = STRESS ON UNIT CELL DUE TO THE EWALD ENERGY
C========================================================================
      DIMENSION DIRC(3,3),RECC(3,3),DIRL(3,3)
      DIMENSION DIRI(3,3),RECI(3,3)
      DIMENSION SIGKE(6),SIGEI(6),SIGHA(6),SIGNL(6)
      DIMENSION SIGXC(6),SIGAL(6),SIGEW(6),SIGTO(6)
C========================================================================
C    FOR THE USE OF SUBROUTINE "EWALD.FOR" ONLY
C========================================================================
      DIMENSION RFORCE(3,NINSP,NINSP)
      DIMENSION CFORCE(3,NINSP,NINSP)
      DIMENSION FORCEG(3,NGPTS),FORSIG(6,NGPTS),ENERG(NGPTS)
      DIMENSION CPHFX(NGPTX), CPHFY(NGPTY), CPHFZ(NGPTZ)
C========================================================================
C NEW ARRAYS NEEDED FOR CONJUGATE GRADIENTS MINIMISATION OF IONIC
C ENERGIES SUBJECT TO FIXED ELECTRONIC CONFIGURATION
C
C GRION(3,NIONS*NSPEC) ESSENTIALLY THE FORCES FOR THE PRESENT IONIC POSITIONS
C      THESE BEING THE STEEPEST DESCENTS DIRECTIONS FOR MOVING THE IONS
C GRIONP(3,NIONS*NSPEC) THE SAME BUT FOR THE PREVIOUS IONIC POSITIONS
C DIRION(3,NIONS*NSPEC) THE CONJUGATE DIRECTION FOR MOVING THE IONS
C DIPION(3,NIONS*NSPEC) THE PREVIOUS CONJUGATE DIRECTION FOR MOVING THE IONS
C POSIOT(3,NIONS,NSPEC) THE TRIAL POSITIONS FOR DETERMINING STEP LENGTH
C POSITC(3,NIONS*NSPEC) SAME INFO IN COMPACT FORM
C
C  ORIGINAL VERSION OF GRAHAM FRANCIS, IMPLEMENTED BY VICTOR MILMAN
C========================================================================
      DIMENSION CVTRUE(NPLWV)
      DIMENSION CVOLD(NPLWV)
      DIMENSION CVIONT(NPLWV)
      DIMENSION GRION(3,NINSP)
      DIMENSION GRIONP(3,NINSP)
      DIMENSION DIRION(3,NINSP)
      DIMENSION DIPION(3,NINSP)
      DIMENSION GRNEW(3*NINSP)
      DIMENSION GRNEWP(3*NINSP)
      DIMENSION DIRNEW(3*NINSP)
      DIMENSION DIPNEW(3*NINSP)
      DIMENSION POSIOT(3,NIONS,NSPEC)
      DIMENSION POSITC(3,NINSP)
      DIMENSION RCNSTR(3*NINSP,3*NINSP)
      DIMENSION VNLT(NBANDS,NKPTS)
C========================================================================
C  FOR THE USE OF SUBROUTINE "ROSYM.FOR" ONLY  (A.QTEISH)
C========================================================================
C      PARAMETER (NDIM9=10,NDIM10=10)
      DIMENSION LPCTXI(-NGX/2:NGX/2)
      DIMENSION LPCTYI(-NGY/2:NGY/2)
      DIMENSION LPCTZI(-NGZ/2:NGZ/2)
      INTEGER F0(48,NDIM9),LWGHT(NDIM10),IB(48),ITEXT(80)
      DIMENSION WVLIST(3,NDIM10),WVKO(3)
      DIMENSION V(3,48),RB(48,3,3),R(49,3,3)
      DIMENSION A1(3),A2(3),A3(3),B1(3),B2(3),B3(3)
      DIMENSION FORCEWK1(3,NINSP)
      DIMENSION FORCEWK2(3,NINSP)
C  USE MACHINE NAMES IN ALL PROGRAMS:
      COMMON /MACH1/ MACTYP             
      COMMON /MACH2/ MACHIN(0:10)       
      CHARACTER*7  MACHIN               
      DATA TWOPI /6.2831853072/
C
      MACHIN(0) = 'UNKNOWN'             
      MACHIN(1) = 'IBM'                 
      MACHIN(2) = 'CRAY'                
      MACHIN(3) = 'VAX'                 
      MACHIN(4) = 'CDC'                 
      MACHIN(5) = 'UNIVAC'              
      MACHIN(6) = 'VAX+FPS'             
      MACHIN(7) = 'CRAY-2'              
      MACHIN(8) = 'FUJITSU'             
      MACHIN(9) = 'VAX+MAP'             
C
C
      WRITE (*,*)'   ------------------------------------------------'
      WRITE (*,*)'   =             Welcome to CASTEP                ='
      WRITE (*,*)'   =     CAmbridge Serial Total Energy Package    ='
      WRITE (*,*)'   =     =====================================    ='
      WRITE (*,*)'   =     Ab initio Molecular Dynamics Program     ='
      WRITE (*,*)'   =  with Conjugate Gradient Minimisation (CG)   ='
      WRITE (*,*)'   =          for the Electronic Energy           ='
      WRITE (*,*)'   =                                              ='
      WRITE (*,*)'   =     Version of: 8 July 1992                  ='
      WRITE (*,*)'   ------------------------------------------------'
      WRITE (*,*)'  '
      WRITE (*,1859) NGX,NGY,NGZ
 1859 FORMAT(1X,' FOURIER TRANSFORM GRID: ',I3,'x',I3,'x',I3)
C
C
C
C************ THE FOLLOWING PARAMETERS ARE USER VARIABLE ***************
C
C***** THE VARIABLES ARE DEFINED IN THE COMMON BLOCK SECTION ABOVE *****
C
C UNIT 3  = CRYSTALLOGRAPHIC INFORMATION -NEEDED FOR SYMMETRIZATION ONLY
C UNIT 6  = OUTPUT FILE
C UNIT 10 = DATA FOR INTERPOLATING THE EWALD INTEGRAL
C UNIT 11 = DATA FOR INTERPOLATING THE PSEUDOPOTENTIAL
C UNIT 12 = USED TO STORE WAVEFUNCTIONS AT THE END OF THE PROGRAM AND TO
C           READ IN WAVEFUNCTIONS IF THEY ALREADY EXIST
C UNIT 13 = USED TO STORE IONIC POSITIONS, COORDINATES OF THE SPECIAL
C           K POINTS, THE SIZE OF THE UNIT CELL, THE ELECTRONIC
C           EIGENVALUES AT THE END OF PROGRAM EXECUTION, THE COORDINATES
C           OF THE SPECIAL K POINTS AND THEIR WEIGHTS  AND TO READ IN
C           THESE VALUES IF WRITTEN BY A PREVIOUSLY EXECUTED PROGRAM
C UNIT 14 = FORMATTED FILE USED TO READ THE NUMBER OF ITERATIONS, THE
C           LENGTH OF THE TIMESTEPS, DAMPING FACTORS, ETC., ETC.
C UNIT 15 = FORMATTED FILE USED TO READ IN THE SIZE OF THE UNIT CELL,
C           THE POSITIONS OF THE IONS, THE COORDINATES OF THE SPECIAL
C           K POINTS AND THE WEIGHT OF EACH K POINT
C UNIT 16 = USED TO STORE THE ELECTRONIC CHARGE DENSITY AT THE END OF
C           A SET OF ITERATIONS
C UNIT 17 = USED TO STORE WAVEFUNCTIONS AT THE END OF THE PROGRAM AND TO
C           READ IN WAVEFUNCTIONS IF THEY ALREADY EXIST
C
C
C****************** READ IN DATA FROM EXTERNAL UNITS *******************
C
C        
C=========================================================================
C     READ IN EWALD DATASET FROM UNIT 10
C=========================================================================
      OPEN (10,FORM='FORMATTED')
      REWIND 10
      READ (10,*) (EWRLEN(I),I=1,NEWPTS)
      READ (10,*) (EWRCSS(I),I=1,NEWPTS)
      READ (10,*) (EWRLSS(I),I=1,NEWPTS)
      READ (10,*) (EWRCSI(I),I=1,NEWPTS)
      CLOSE (10)
C=========================================================================
C     READ IN FROM UNIT 14
C=========================================================================
      OPEN (14,FORM='FORMATTED')
      REWIND 14
      READ (14,*) NITER
      READ (14,*) NPRINT
      READ (14,*) ENMAX
      READ (14,*) NLPOT
      READ (14,*) NITMAX
      READ (14,*) NDELAY
      READ (14,*) ISTART
      READ (14,*) ISBROT
      READ (14,*) IOCCUP
      READ (14,*) DELMIN
      READ (14,*) DELMAX
      READ (14,*) NDEL
      WRITE (*,*)' ISBROT: ',ISBROT,' IOCCUP: ',IOCCUP
      IF (IOCCUP.NE.0) THEN
        WRITE (*,*)' DELMIN, DELMAX:',DELMIN,DELMAX
        WRITE (*,*)' NDEL (FREQUENCY OF DEL HALVING): ',NDEL
      ENDIF
      READ (14,*) IION
      READ (14,*) IBOX
      WRITE (*,*)' RELAXATION: IONS= ',IION,'  UNIT CELL= ',IBOX
      READ (14,*) IPRINT
      READ (14,*) NIONCG
      WRITE (*,*)' NIONCG: ',NIONCG
      READ (14,*) NITFIX
      WRITE (*,*)' NITFIX: ',NITFIX
      READ (14,*) INRAND
      WRITE (*,*)' INRAND: ',INRAND
      READ (14,*) ICLOCK
      WRITE (*,*)' ICLOCK: ',ICLOCK,' IPRINT: ',IPRINT
      IF (ICLOCK.EQ.1) CALL SCLOCK
      READ (14,*) SITIM
      READ (14,*) SIDAMP
      READ (14,*) SIMASS
      READ (14,*) SIDISP
      READ (14,*) POTIM
      READ (14,*) PODISP
      WRITE (*,*) '  ENMAX =', ENMAX
      WRITE (*,*) '  SITIM =', SITIM
      WRITE (*,*) '  POTIM =', POTIM
      DO 452 NSP = 1 , NSPEC
        READ (14,*) IDUMMY
        ICHARG(NSP) = IDUMMY
 452  CONTINUE
      DO 454 NSP = 1 , NSPEC
        READ (14,*) DUMMY
        POMASS(NSP) = DUMMY
 454  CONTINUE
      DO 456 NSP = 1 , NSPEC
        READ (14,*) IDUMMY
        NIONSP(NSP) = IDUMMY
        IF (NIONSP(NSP).GT. NIONS) STOP
     &    'NIONSP(NSP)>NIONS! INCREASE NIONS IN THE PARAMETER STATEMENT'
 456  CONTINUE
C========================================================================
C  IVPTYP = 0 IF ALL POTENTIALS ARE LOCAL
C           1    THERE IS AT LEAST ONE NONLOCAL PSEUDOPOTENTIAL
C========================================================================
      IVPTYP=0
      DO 457 NSP=1,NSPEC
        READ(14,*) IDUMMY
        IVPTYN(NSP) = IDUMMY
        IVPTYP = IVPTYP + IDUMMY * IDUMMY
        WRITE(*,*)' IVPTYN ',IVPTYN(NSP), ' NIONSP ',NIONSP(NSP)
 457  CONTINUE
      IF (IVPTYP.GT.0) IVPTYP = 1
      IF (IVPTYP.EQ.0) NLPOT = 0
C========================================================================
C  READ DIMENSIONS OF PSEUDOPOTENTIAL ARRAYS
C========================================================================
      DO 458 NSP=1,NSPEC
        READ(14,*) IDUMMY
        NPSPTN(NSP)=IDUMMY
        WRITE(*,*)' NPSPTN',NPSPTN(NSP)
        IF (IDUMMY.GT.NPSPTS) STOP'NPSPTN>NPSPTS - INCREASE NPSPTS!'
 458  CONTINUE
      READ (14,*) SCAL
      READ (14,*) ISYMM
      READ (14,*) IBANS
      IF (IBANS.EQ.1) THEN
        ISTART = 1
        WRITE (*,*) ' IBANS=1, SO WE READ CHARGE DENSITY FROM FILE 13'
      END IF
      IF (IION.NE.0) THEN
        WRITE (*,*) ' ICNSTR = 0 (NO CONSTRAINTS) OR 1 (CONSTRAINTS)'
        ICNSTR = 0
        READ (14,*,END=459,ERR=459) IDUMMY
        ICNSTR = IDUMMY
        WRITE (*,*) ICNSTR
        GO TO 455
 459    CONTINUE
        WRITE (*,*) ' UNCONSTRAINED MINIMISATION (DEFAULT VALUE)'
 455    CONTINUE
      END IF
C=========================================================================
C  NEXTWR   -  FREQUENCY OF WRITING TO EXTERNAL FILES (13,16,17,28 ETC.)
C              DATA FOR THE MD RESTART ARE WRITTEN TO 13 FILE
C=========================================================================
      READ (14,*,ERR=8330,END=8330) NEXTWR
      GO TO 8331
 8330 NEXTWR = 5
 8331 IF (NEXTWR.EQ.0) NEXTWR = 1
      WRITE (*,*) ' NEXTWR=',NEXTWR
      CLOSE (14)
C=========================================================================
C CALCULATE THE TOTAL NUMBER OF ELECTRONS, NELECT, THE TOTAL CHARGE ON
C THE IONS, NIONCH, AND THE SUM OF THE SQUARES OF THE CHARGES ON THE
C IONS, NICHSQ, AND THE TOTAL NUMBER OF IONS, NIONST.
C=========================================================================
      NELECT = 0
      NIONCH = 0
      NICHSQ = 0
      NIONST1 = 0
      DO 3105 NSP = 1 , NSPEC
        NELECT = NELECT + ICHARG(NSP) * NIONSP(NSP)
        NIONCH = NIONCH + ICHARG(NSP) * NIONSP(NSP)
        NICHSQ = NICHSQ + ICHARG(NSP) * ICHARG(NSP) * NIONSP(NSP)
        NIONST1 = NIONST1 + NIONSP(NSP)
 3105 CONTINUE
      IF (NIONST1.NE.NIONST)STOP'CHECK NIONST!'
      IF (NBANDS*2.LT.NELECT)STOP'NBANDS*2 < NELECT, INCREASE NBANDS' 
C=========================================================================
C  READ DATA FOR SYMMETRIZATION
C=========================================================================
      IF (ISYMM.NE.0) THEN
        WRITE (*,*) ' READ INFORMATION NEEDED FOR SYMMETRIZATION'
        CALL SPPTRD(ITEXT,IHG,IHC,ISY,LI,NC,IB,V,F0,R,NIONST,NDIM9,
     &              NDIM10,NUMKPT,IQ1,IQ2,IQ3,WVKO,LWGHT,WVLIST,
     &              A1,A2,A3,B1,B2,B3)
        WRITE (*,*) ' NAT=',NIONST,' NDIM9=',NDIM9,' NDIM10=',NDIM10
        WRITE (*,*) ' IQ 1-3: ',IQ1,IQ2,IQ3
        WRITE (*,*) ' ISY,NC: ',ISY,NC
        WRITE (*,*) ' FILE CONTAINS ',NUMKPT,' SPECIAL POINTS:'
        ISUM = 0
        DO 4020 I = 1 , NUMKPT
          ISUM = ISUM + LWGHT(I)
          DO 4021 M = 1 , 3
            WORK3(M,1) = 0.0
 4021     CONTINUE
          DO 4022 M = 1 , 3
            WORK3(1,1) = WORK3(1,1) + WVLIST(M,I) * A1(M)
            WORK3(2,1) = WORK3(2,1) + WVLIST(M,I) * A2(M)
            WORK3(3,1) = WORK3(3,1) + WVLIST(M,I) * A3(M)
 4022     CONTINUE
          DO 4023 M = 1 , 3
            WVLIST(M,I) = WORK3(M,1)
 4023     CONTINUE
 4020   CONTINUE
        SUMWEI = 1.0 / FLOAT(ISUM)
        DO 4024 I = 1 , NUMKPT
          WEIGHT = FLOAT(LWGHT(I)) * SUMWEI
          WRITE (*,4025) I,WEIGHT,(WVLIST(M,I),M=1,3)
 4024   CONTINUE
        SUMWEI = 2.0 * SUMWEI
        WRITE (*,*) ' SPPTRD FINISHED'
        WRITE (*,*) ' NATURAL WEIGHT: ',SUMWEI
      END IF
 4025 FORMAT(1X,'NKP:',I3,' WEIGHT=',F10.7,' K=(',3F10.7,')')
C=========================================================================
C     READ IN PSEUDOPOTENTIAL DATASET FROM UNIT 11  (X. WENG 14-JUN-90)
C=========================================================================
      OPEN (11,FORM='FORMATTED')
      REWIND 11
      DO 10 NSP=1,NSPEC
        IF(IVPTYN(NSP).EQ.0) THEN
C=========================================================================
C     1) READ IN LOCAL PSEUDOPOTENTIAL                
C     PSCORE MUST BE STORED AS THE SECOND ELEMENT OF THE PSEUDOPOTENTIAL,
C     AND PSGMAX AS THE FIRST ELEMENT 
C=========================================================================
          READ(11,*) PSGMAX(NSP)
          READ(11,*) (PSP(I,NSP), I=1,NPSPTN(NSP))
          PSCORE(NSP)=PSP(1,NSP)
          WRITE(*,*)' PSGMAX(NSP) =',PSGMAX(NSP)
          WRITE(*,*)' PSCORE(NSP) =',PSCORE(NSP)
          PSP(1,NSP)=0.0
          PSP(2,NSP)=0.0
        ELSE
C=========================================================================
C     2) READ IN NON-LOCAL PSEUDOPOTENTIAL
C     PSCORE MUST BE STORED AS THE FIRST ELEMENT OF THE PSEUDOPOTENTIAL,
C=========================================================================
          LOCAL(NSP)=-1
          READ(11,*) PSGMAX(NSP)
          READ(11,*) (PSP(J,NSP),J=1,NPSPTN(NSP))
          DO 40 L=0,2
            READ(11,*) PSCALE(L,NSP)
            IF (ABS(PSCALE(L,NSP)).LT.0.00001) LOCAL(NSP)=L
            READ(11,*) (PSPNL(J,L,NSP),J=1,NPSPTN(NSP))
   40     CONTINUE
          PSCORE(NSP)=PSP(1,NSP)
          PSP(1,NSP)=0.0
          PSP(2,NSP)=0.0
          WRITE(*,*)' PSCORE(NSP)  =', PSCORE(NSP)
          WRITE(*,*)' PGMAX(NSP)   =', PSGMAX(NSP)
          WRITE(*,*)' PSCALE(0,NSP)=', PSCALE(0,NSP)
          WRITE(*,*)' PSCALE(1,NSP)=', PSCALE(1,NSP)
          WRITE(*,*)' PSCALE(2,NSP)=', PSCALE(2,NSP)
          WRITE(*,*)'THE LOCAL POTENTIAL IS THE ONE WITH L=',LOCAL(NSP)
C=======================================================================
C  IF NLPOT=1 WE USE REAL-SPACE IMPLEMENTATION OF NON-LOCAL POTENTIAL
C=======================================================================
          IF (NLPOT.EQ.1) THEN
            READ(18,*) RLCORE(NSP)
            READ(18,*) RMAX(NSP)
C========================================================================
C THE ROUTINE CAN HANDLE ARBITRARY NUMBERS OF NON-LOCAL POTENTIALS
C FOR L=0,1,2 SO ENHANCED KLEINMAN-BYLANDER PROJECTIONS OR VANDERBILT
C POTENTIALS SHOULD CAUSE FEW PROBLEMS
C========================================================================
C THE TOTAL NUMBER OF PROJECTIONS (INCLUDING M VALUES) IS RECORDED BY
C NSPHAR - START BY SETTING IT TO 0
C========================================================================
            NSPHAR=0
C========================================================================
C NOW START READING THE NON-LOCAL POTENTIALS, THE FIRST VALUE READ IS
C THE L VALUE, A VALUE OF 1000 SHOULD BE USED TO MARK THE END OF THE
C POTENTIAL FILE AT WHICH POINT THE PROGRAM MOVES TO THE NEXT SPECIES
C OF ION AFTER RECORDING THE VALUE OF NSPHAR AS NSHSP(NSP)
C========================================================================
            DO 31 NRLIND=1,MXRLNL
              READ(18,*) ISP
              IF(ISP.EQ.1000) GO TO 32
              IRLNL(NRLIND,NSP)=ISP
              NSPHAR=NSPHAR+(1+2*ISP)
              READ(18,*) PRLSCA(NRLIND,NSP)
              READ(18,*) (VRLNL(I,NRLIND,NSP),I=1,NRLPTS)
              WRITE (*,*) '  PRLSCA(NRLIND,NSP) =',PRLSCA(NRLIND,NSP)
              WRITE (*,*) '  IRLNL(NRLIND,NSP)  =',IRLNL(NRLIND,NSP)
              NRLNL(NSP)=NRLIND
 31         CONTINUE
 32         CONTINUE
C========================================================================
C CHECK THAT ARRAYS ARE LARGE ENOUGH TO HOLD ALL SPHERICAL HARMONICS
C========================================================================
            WRITE (*,*) '  NRLNL(NSP)  =',NRLNL(NSP)
            IF(NSPHAR.GT.MXRLSH) THEN
              WRITE(*,*) ' MXRLSH NOT LARGE ENOUGH  '
            END IF
            NSHSP(NSP)=NSPHAR
          END IF
        END IF
C========================================================================
C END OF LOOP FOR ION NSP
C========================================================================
   10 CONTINUE
      CLOSE (11)
C========================================================================
      POTIM0 = POTIM
      TIMEMD = 0.0
      EK = 0.0
C========================================================================
C INSTALL NUMBER OF CONJUGATE GRADIENTS ITERATIONS 
C========================================================================
      IF (NIONCG.EQ.0) THEN
        NIOCG1 = 3
        IF (IOCCUP.EQ.1) NIOCG1 = 2
      ELSE
        NIOCG1 = NIONCG
      END IF
      NIOCG0 = NIOCG1
C========================================================================
C     READ IN THE UNIT CELL VECTORS 
C     1.THE CURRENT ONES (DIRC) AND THE INITIAL ONES (DIRI)
C========================================================================
      OPEN (15,FORM='FORMATTED')
      REWIND 15
      READ (15,*) DIRC(1,1), DIRC(2,1), DIRC(3,1)
      READ (15,*) DIRC(1,2), DIRC(2,2), DIRC(3,2)
      READ (15,*) DIRC(1,3), DIRC(2,3), DIRC(3,3)
      CALL BASTR(DIRC,RECC,VOLC)
      READ (15,*) DIRI(1,1), DIRI(2,1), DIRI(3,1)
      READ (15,*) DIRI(1,2), DIRI(2,2), DIRI(3,2)
      READ (15,*) DIRI(1,3), DIRI(2,3), DIRI(3,3)
      CALL BASTR(DIRI,RECI,VOLI)
C
      WRITE (*,*) 'THE BASIS-SET USED WHEN PLACING ATOMS ( DIRC(I,J) ):'
      WRITE (*,101) ( (DIRC(I,J), I=1,3), J=1,3)
 101  FORMAT(1X,3F10.6)
      WRITE (*,*) 'THE RECIPROCAL SPACE BASIS ( RECC(I,J) ): '
      WRITE (*,101) ( (RECC(I,J), I=1,3), J=1,3)
      WRITE (*,*) 'THE BASIS-SET USED WHEN CUTTING OFF PLANE-WAVES ',
     &        '( DIRI(I,J) ): '
      WRITE (*,101) ( (DIRI(I,J), I=1,3), J=1,3)
 7007 FORMAT(1X,3F10.6)
      WRITE (*,*) ' VOLUME OF CURRENT UNIT CELL:', VOLC
      WRITE (*,*) ' VOLUME OF INITIAL UNIT CELL:', VOLI
C========================================================================
C     2.READ IN ATOMIC POSITIONS (SAME AS PAYNE'S ORIGINAL ONE)
C========================================================================
      NG = 1
      WRITE (*,*) ' '
      DO 7006 NSP = 1 , NSPEC
        DO 7008 NI = 1 , NIONSP(NSP)
          READ (15,*)    (POSION(I,NI,NSP), I=1,3), RMOVE(NG)
          WRITE (*,7004) NI,NSP,(POSION(I,NI,NSP), I=1,3), RMOVE(NG)
          NG = NG + 1
 7008   CONTINUE
        DO 7010 NI = 1 , NIONSP(NSP)
          READ (15,*)    (POSIOL(I,NI,NSP), I=1,3)
 7010   CONTINUE
 7006 CONTINUE
      DO 7011 NKP = 1 , NKPTS
        READ (15,*) VKPT(1,NKP),VKPT(2,NKP),VKPT(3,NKP)
 7011 CONTINUE
      WRITE (*,*) ' '
      DO 7012 NKP = 1 , NKPTS
        READ (15,*) WTKPT(NKP)
        WRITE (*,7014) VKPT(1,NKP),VKPT(2,NKP),VKPT(3,NKP), WTKPT(NKP)  
 7012 CONTINUE
      WRITE (*,*) ' '
      CLOSE (15)
C
 7004 FORMAT(2X,'ION ',I2,' TYPE ',I1,' AT (',3F9.6,')   RMOVE =',F5.2)
 7009 FORMAT(1X,4F10.6)
 7014 FORMAT('  SPECIAL K-POINT:',3F10.5,'   WEIGHT =', F10.5)
 7019 FORMAT(1X,3F10.6)
C=========================================================================
C     FOR RE-START, READ IN FROM FOR013 AND FOR016      (X.WENG 05-MAR-90)
C=========================================================================
      IF (ISTART.NE.1 .AND. IBANS.EQ.0) GO TO 6006
      WRITE (*,*)' ----------------------------------------------------'
      WRITE (*,*)'     SINCE ISTART=1, READ IN FROM FOR013,16,17,28....'
      WRITE (*,*)' ----------------------------------------------------'
      REWIND 13
      READ (13) ( (DIRC(I,J), I=1,3), J=1,3 )
      READ (13) VOLC
      READ (13) ( (RECC(I,J), I=1,3), J=1,3 )
      READ (13) ( (DIRI(I,J), I=1,3), J=1,3 )
      READ (13) ( ( (POSIOL(I,J,K), I=1,3), J=1,NIONS), K=1,NSPEC )
      READ (13) ( ( (POSION(I,J,K), I=1,3), J=1,NIONS), K=1,NSPEC )
      READ (13) (RMOVEO(J), J=1,NIONS*NSPEC)
      IF (IBANS.NE.0) THEN
        NGXO = NGX
        NGYO = NGY
        NGZO = NGZ
        NRPLWVO = NRPLWV
        GO TO 6005
      END IF
      READ (13) ( (CELEN(I,J),I=1,NBANDS), J=1,NKPTS)
      READ (13) ((VKPT(I,J),I=1,3),J=1,NKPTS)
      READ (13) (WTKPT(I),I=1,NKPTS)
      READ (13) ENMAXO
      READ (13) NGXO,NGYO,NGZO,NRPLWO
      READ (13) ( (RECO(I,J), I=1,3), J=1,3)
C
      WRITE (*,162) ( (DIRC(I,J), I=1,3), J=1,3)
      WRITE (*,162) ( (DIRI(I,J), I=1,3), J=1,3)
  162 FORMAT(2X,3F12.5)
63    FORMAT(1X,' KINETIC ENERGY IS ',2F10.6)
62    FORMAT(1X,6F10.6)
      WRITE (*,*) ' VOLC =',VOLC
      NG = 1
      DO 6000 K = 1 , NSPEC
        DO 6002 J = 1 , NIONSP(K)
          WRITE (*,62) (POSION(I,J,K), I=1,3), RMOVEO(NG)
          NG = NG + 1
 6002   CONTINUE  
        DO 6004 J = 1 , NIONSP(K)
          WRITE (*,62) (POSIOL(I,J,K),I=1,3)
 6004   CONTINUE
 6000 CONTINUE
      IF (IPRINT.GE.2)
     &           WRITE (*,63) ( (CELEN(I,J),I=1,NBANDS), J=1,NKPTS)
      WRITE (*,101) ( (VKPT(I,J),I=1,3), J=1,NKPTS)
      WRITE (*,102) (WTKPT(I),I=1,NKPTS)
 102  FORMAT(1X,F10.6)
C
 6005 CONTINUE
      REWIND 16
      NPLWVO = NGXO * NGYO * NGZO
      READ (16) (CHDENR(I),I=1,NPLWVO), (CHDENG(I),I=1,NPLWVO)
 6006 CONTINUE
C=========================================================================
C SET THE SIZE OF THE BOX AT THE LAST TIMESTEP TO THE VALUE AT THE
C PRESENT TIMESTEP AND SET THE ELECTRONIC ENERGY EIGENVALUES TO ZERO TO
C COMMENCE THE ELECTRONIC RELAXATION
C=========================================================================
      DO 446 J = 1 , 3
        DO 445 I = 1 , 3
          DIRL(I,J) = DIRC(I,J)
 445    CONTINUE
 446  CONTINUE
      IF (ISTART.EQ.0 .OR. IBANS.NE.0) THEN
        DO 35 I = 1 , NKPTS
          DO 36 II = 1 , NBANDS
            CELEN(II,I) = (0.0,0.0)
            EIGVAL(II,I) = 0.0
 36       CONTINUE
 35     CONTINUE
      END IF
C=========================================================================
C FORM PACKED ARRAYS CONTAINING THE MASSES OF THE IONS, POMASC, AND THE
C CHARGES ON THE IONS, ICHARC.
C=========================================================================
      NINDX = 1
      DO 3106 NSP = 1 , NSPEC
        DO 3107 NI = 1 , NIONSP(NSP)
          DO 3108 M = 1 , 3
            POSIC(M,NINDX) = POSION(M,NI,NSP)
 3108     CONTINUE
          POMASC(NINDX) = POMASS(NSP)
          ICHARC(NINDX) = ICHARG(NSP)
          NINDX = NINDX + 1
 3107   CONTINUE
 3106 CONTINUE
C=========================================================================
C     INITIALISE (ISTART=0) OR READ IN (ISTART=1) THE OCCUPATION NUMBER 
C=========================================================================
      IF (NELECT/2*2.NE.NELECT .AND. IOCCUP.EQ.0) WRITE (*,9988)
 9988 FORMAT(' WARNING - YOU HAVE ODD NUMBER OF ELECTRONS.'/
     &       ' THE LAST BAND WILL BE HALF-FILLED AUTOMATICALLY.'/
     &       ' CHECK YOUR RESULTS FOR PHYSICAL MEANINGNESS!!'/
     &       ' (PERHAPS IOCCUP=1 WOULD BE A BETTER CHOICE)')
C
      DO 5210 K = 1 , NKPTS
        DO 5211 I = 1 , NBANDS
          OCC(I,K) = 0.0
 5211   CONTINUE
        DO 5212 I = 1 , NELECT/2
          OCC(I,K) = 1.0                                   
 5212   CONTINUE
        IF (NELECT/2*2.NE.NELECT) OCC(NELECT/2+1,K) = 0.5      
 5210 CONTINUE 
      IF (ISTART*IOCCUP.NE.0) THEN
        REWIND 17
        DO 5240 NKP = 1 , NKPTS
          DO 5220 NBD = 1 , NBANDS
            READ (17,*) I,J,DUM
            IF (I.NE.NKP.OR.J.NE.NBD)
     &                   STOP'MISMATCH WHEN READING OCC FROM UNIT 17'
            OCC(NBD,NKP) = DUM
 5220     CONTINUE
 5240   CONTINUE
      END IF
C=========================================================================
C INITIALIZE THE LOOP COUNTERS LPCTX,LPCTY,LPCTZ AND LPCTFX,ETC THAT
C LABEL THE NUMBER OF THE RECIPROCAL LATTICE VECTORS IN THE X,Y,Z
C DIRECTIONS, RESPECTIVELY. FOR THE X DIRECTION THE RECIPROCAL LATTICE
C VECTORS CORRESPONDING TO THE FIRST,SECOND,...,NGXTH ELEMENTS IN ALL
C OF THE RECIPROCAL LATTICE ARRAYS ARE 0,1,..,(NGX/2),-((NGX/2-1),..,-1
C TIMES THE X RECIPROCAL LATTICE VECTOR RECC(1,*)
C=========================================================================
      DO 7020 NX = 1 , (NGX/2)+1
        LPCTX(NX)  = NX - 1
        LPCTFX(NX) = NX - 1
 7020 CONTINUE
      DO 7021 NX = (NGX/2)+2 , NGX
        LPCTX(NX)  = NX - 1 - NGX
        LPCTFX(NX) = NX - 1 - NGX
 7021 CONTINUE
      DO 7022 NY = 1 , (NGY/2)+1
        LPCTY(NY)  = NY - 1
        LPCTFY(NY) = NY - 1
 7022 CONTINUE
      DO 7023 NY = (NGY/2)+2 , NGY
        LPCTY(NY)  = NY - 1 - NGY
        LPCTFY(NY) = NY - 1 - NGY
 7023 CONTINUE
      DO 7024 NZ = 1 , (NGZ/2)+1
        LPCTZ(NZ)  = NZ - 1
        LPCTFZ(NZ) = NZ - 1  
 7024 CONTINUE
      DO 7025 NZ = (NGZ/2)+2 , NGZ
        LPCTZ(NZ)  = NZ - 1 - NGZ
        LPCTFZ(NZ) = NZ - 1 - NGZ
 7025 CONTINUE
C=========================================================================
C SET THE UNBALANCED RECIPROCAL LATTICE VECTOR COMPONENTS TO ZERO IN THE
C LPCTFX,ETC ARRAYS TO REMOVE THE UNPHYSICAL CONTRIBUTION TO THE
C HELLMANN-FEYNMAN FORCE FROM THE NGX/2 RECIPROCAL LATTICE VECTOR WHICH
C IS MATHEMATICALLY EQUALLY THE -NGX/2 LATTICE VECTOR DUE TO THE
C PERIODICITY INTRODUCED BY THE FAST FOURIER TRANSFORM
C=========================================================================
      LPCTFX((NGX/2)+1) = 0
      LPCTFY((NGY/2)+1) = 0
      LPCTFZ((NGZ/2)+1) = 0
C=======================================================================
C      CHECK THE GRID SIZE
C 	TPI = (NOT SURPRISINGLY) 2*PI
C 	HSQDTM = (PLANCK'S CONSTANT/(2*PI))**2/(2*ELECTRON MASS)
C              = 3.810033 EV/ANGSTROM**2
C       ENERGY=HSQDTM*(GY**2)
C   Changed by Carsten Kruse and is valid for arbitrary cells
C
C======================================================================
      HSQDTM = 3.810033
      NGX1=2*SQRT(ENMAX/HSQDTM*(DIRC(1,1)**2+DIRC(1,2)**2+DIRC(1,3)**2))
     *     /TWOPI
      NGY1=2*SQRT(ENMAX/HSQDTM*(DIRC(2,1)**2+DIRC(2,2)**2+DIRC(2,3)**2))
     *     /TWOPI
      NGZ1=2*SQRT(ENMAX/HSQDTM*(DIRC(3,1)**2+DIRC(3,2)**2+DIRC(3,3)**2))
     *     /TWOPI
C     NGX1=2*SQRT(ENMAX/HSQDTM/(RECC(1,1)**2+RECC(1,2)**2+RECC(1,3)**2))
C     NGY1=2*SQRT(ENMAX/HSQDTM/(RECC(2,1)**2+RECC(2,2)**2+RECC(2,3)**2))
C     NGZ1=2*SQRT(ENMAX/HSQDTM/(RECC(3,1)**2+RECC(3,2)**2+RECC(3,3)**2))
      NGX1 = 2 * NGX1
      NGY1 = 2 * NGY1
      NGZ1 = 2 * NGZ1
C======================================================================
C  Skip recommended values of NG which has simple factors other
C  than 2,3, and 5. Fast Fourier Transform can't work in its present
C  form with bigger factors.
C======================================================================
      IF (NGX1.EQ.14 .OR. NGX1.EQ.22 .OR. NGX1.EQ.28 .OR. 
     &    NGX1.EQ.34 .OR. NGX1.EQ.38 .OR. NGX1.EQ.46 .OR.
     &    NGX1.EQ.52. OR. NGX1.EQ.58 .OR. NGX1.EQ.62) NGX1 = NGX1 + 2
      IF (NGX1.EQ.26. OR. NGX1.EQ.44 .OR. NGX1.EQ.56) NGX1 = NGX1 + 4
      IF (NGX1.EQ.42) NGX1 = NGX1 + 6
C
      IF (NGY1.EQ.14 .OR. NGY1.EQ.22 .OR. NGY1.EQ.28 .OR. 
     &    NGY1.EQ.34 .OR. NGY1.EQ.38 .OR. NGY1.EQ.46 .OR.
     &    NGY1.EQ.52. OR. NGY1.EQ.58 .OR. NGY1.EQ.62) NGY1 = NGY1 + 2
      IF (NGY1.EQ.26. OR. NGY1.EQ.44 .OR. NGY1.EQ.56) NGY1 = NGY1 + 4
      IF (NGY1.EQ.42) NGY1 = NGY1 + 6
C
      IF (NGZ1.EQ.14 .OR. NGZ1.EQ.22 .OR. NGZ1.EQ.28 .OR. 
     &    NGZ1.EQ.34 .OR. NGZ1.EQ.38 .OR. NGZ1.EQ.46 .OR.
     &    NGZ1.EQ.52. OR. NGZ1.EQ.58 .OR. NGZ1.EQ.62) NGZ1 = NGZ1 + 2
      IF (NGZ1.EQ.26. OR. NGZ1.EQ.44 .OR. NGZ1.EQ.56) NGZ1 = NGZ1 + 4
      IF (NGZ1.EQ.42) NGZ1 = NGZ1 + 6
C
      WRITE (*,3810) NGX1, NGY1, NGZ1
3810  FORMAT(/,' THE MINIMUM GRID SHOULD BE',I4,' x',I4,' x',I4)
      IF (NGX.GE.NGX1 .AND. NGY.GE.NGY1 .AND. NGZ.GE.NGZ1) THEN
        WRITE (*,*)' THE CURRENT GRID IS LARGE ENOUGH FOR THIS ENMAX.'
      ELSE
        WRITE (*,*)'**> THE CURRENT GRID IS TOO SMALL FOR THIS ENMAX<**'
        WRITE (*,*)'                          ^^^^^^^^^                '
      END IF
C=========================================================================
C  INITIALIZE THE DATA USED FOR THE FAST FOURIER TRANSFORMS
C=========================================================================
      NGPTAR(1) = NGX
      NGPTAR(2) = NGY
      NGPTAR(3) = NGZ
C=========================================================================
C
C ]]]]]]]]]]]]]]]]]]]]]] SUBROUTINE GENSP [[[[[[[[[[[[[[[[[[[[[[[[[[[[[[
C
C SUBROUTINE GENSP PERFORMS A NUMBER OF TASKS. THE DATA ARRAYS COSC AND
C DECAY FOR THE ELECTRON DYNAMICS ARE INITIALISED. THE INDEXING SYSTEM
C FOR PADDING THE SPHERES OF PLANE WAVES AT EACH K POINT INTO THE BOX
C USED FOR THE FAST FOURIER TRANSFORMS IS ALSO COMPUTED AS ARE THE
C KINETIC ENERGIES OF THE PLANE WAVE BASIS STATES AT EACH K POINT
C=========================================================================
      CALL GENBTR(NRPLWV,NGX,NGY,NGZ,NKPTS,ENMAX,NINDPW,NPLWKP,VKPT,
     &LPCTX,LPCTY,LPCTZ,DATAKE,RECC,RECI,IPRINT,DNLG,DNLKG)
      WRITE (*,*) 'NPLWKP:', NPLWKP
C=========================================================================
C     IF THE CUT-OFF ENERGY OF THE PREVIOUS CALCULATION IS DIFFERENT
C     FROM THE CURRENT ONE, WE HAVE TO RE-ASSIGN THE WAVEFUNCTION 
C     ACCORDING TO THE NEW INDEXING SEQUENCE.
C     N.B. ASSUMING THE FFT GRID (I.E. NGX, NGY, NGZ) AND THE UNIT CELL
C           ARE KEPT THE SAME
C
C     ADDED BY X. WENG 22-JUL-90
C
C     NEWOLD=0 MEANS NO REASSIGNMENT NEEDED WHEN USING O WAVEFUNCTION
C           =1 MEANS    REASSIGNMENT  "      "     "    "     "
C
C=========================================================================
      NEWOLD = 0
      IF (ISTART.NE.0 .AND. IBANS.EQ.0) THEN
        WRITE (*,*)' '
        IF (          (ENMAX.EQ.ENMAXO)
     &      .AND.      (NGX.EQ.NGXO)
     &      .AND.      (NGY.EQ.NGYO)
     &      .AND.      (NGZ.EQ.NGZO)
     &      .AND.   (NRPLWV.EQ.NRPLWO) ) THEN
          NEWOLD = 0
        ELSE 
          NEWOLD = 1
        END IF
C
        IF (NEWOLD.EQ.0) THEN
          WRITE (*,5402)
        ELSE 
          WRITE (*,5404) ENMAX,ENMAXO,NGX,NGXO,NGY,NGYO,NGZ,NGZO,
     &                  NRPLWV,NRPLWO
          CALL REASSN(NKPTS,VKPT,
     &                NRPLWV,  NGX,     NGY,     NGZ,    ENMAX,
     &                NRPLWO,  NGXO,    NGYO,    NGZO,   ENMAXO,
     &                LPCTX,   LPCTY,   LPCTZ,   RECI,   RECO,
     &                IASSGN,CWORK6,CPTWFP,NBANDS,IPRINT)
        END IF       
      END IF
C
 5402 FORMAT(/,' NO REASSIGNMENT NEEDED WHEN USING OLD WAVEFUNCTION',/)
 5404 FORMAT(/,' REASSIGNMENT NEEDED WHEN USING OLD WAVEFUNCTION',   /,
     &         '   ENMAX =', F12.5,5X,'ENMAXO =',F12.5,              /,
     &         '   NGX   =', I6,  11X,'NGXO   =',I6,                 /,
     &         '   NGY   =', I6,  11X,'NGYO   =',I6,                 /,
     &         '   NGZ   =', I6,  11X,'NGZO   =',I6,                 /,
     &         '   NRPLWV=', I6,  11X,'NRPLWO =',I6,                 /,
     &         ' WE CALL "REASSN" TO RE-ASSIGN THE WAVE FUNCTIONS',  /)
C=========================================================================
C SET THE REAL (CHDENR) AND RECIPROCAL (CHDENG) SPACE CHARGE DENSITIES 
C TO ZERO PRIOR TO INITIALISING THE ELECTRONIC WAVEFUNCTIONS AND 
C CALCULATING THE INITIAL CHARGE DENSITIES
C=========================================================================
      IF ( ISTART.NE.1 .OR. NEWOLD.NE.0 ) THEN
        DO 5200 M = 1 , NPLWV  
          CHDENG(M) = (0.0,0.0)
          CHDENR(M) = (0.0,0.0)
 5200   CONTINUE
      END IF
      IF (NLPOT.EQ.1) THEN
C=======================================================================
C CALL SUBROUTINE TO SET UP DATA ARRAYS FOR REAL SPACE NON-LOCAL
C PSEUDOPOTENTIALS
C=======================================================================
        CALL 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)
      END IF                                                             
C=========================================================================
C START THE LOOP OVER K POINTS FOR WAVEFUNCTION INITIALISATION AND
C CHARGE DENSITY GENERATION
C=========================================================================
      DO 4100 NKP = 1 , NKPTS
        IF (ISTART.EQ.1 .AND. IBANS.EQ.1) GO TO 4113
        IF (ISTART.EQ.1 .AND. (NEWOLD.EQ.0) )  GO TO 4112
        IF (ISTART.EQ.1 .AND. (NEWOLD.EQ.1) )  THEN
          NPUNIT = 19 + NKP
          REWIND NPUNIT
          READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
          GO TO 4115
        END IF
C=========================================================================
C INITIALISE THE WAVEFUNCTIONS TO ZERO
C=========================================================================
 4113   CONTINUE
        DO 4110 M = 1 , NBANDS
          DO 4111 MM = 1 , NRPLWV
            CPTWFP(MM,M) = (0.0,0.0)
 4111     CONTINUE
 4110   CONTINUE
C=========================================================================
C
C """"""""""""""""""""""""" SUBROUTINE WFINIT """"""""""""""""""""""""""
C
C THIS SUBROUTINES INITIALISES THE WAVEFUNCTIONS FOR THE MOLECULAR
C DYNAMICS BY FILLING THE BANDS AT EACH K POINT WITH THE LOWEST ENERGY
C PLANE WAVE BASIS STATES. IF THE SYSTEM HAS HIGH SYMMETRY THIS METHOD
C MAY NOT BE APPROPRIATE BECAUSE THE LOWEST ENERGY PLANE WAVES MAY NOT
C SPAN THE LOWEST BANDS.
C (NOW WORKING IN RECIPROCAL SPACE)
C=========================================================================
        IF (INRAND.NE.1 .AND. IBANS.EQ.0) THEN
          CALL WFINTR(NRPLWV,NBANDS,NKPTS,DATAKE(1,1,NKP),NBANOC,CPTWFP,
     &         NPLWKP(NKP),WORK1,NWORK1,SCAL,CELEN(1,NKP),NGX, NGY, NGZ,
     &         LPCTX,LPCTY,LPCTZ,NINDPW(1,NKP), RECC )
        ELSE
          CALL WFINRN(NRPLWV,NBANDS,NKPTS,DATAKE(1,1,NKP),NBANOC,CPTWFP,
     &         NPLWKP(NKP),WORK1,NWORK1,CELEN(1,NKP),NGX, NGY, NGZ,
     &         LPCTX,LPCTY,LPCTZ,NINDPW(1,NKP), RECC,CWOR10,CWOR11)
        END IF
C=========================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE ORSP }}}}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE ORTHOGONALISES THE WAVEFUNCTIONS JUST INITIALISED
C (NOW WORKING IN RECIPROCAL SPACE)
C=========================================================================
C
C        CALL ORSP(NBANDS,NKPTS,NRPLWV,NPLWKP(NKP),CPTWFP,
C     &                                        CWORK1,CWORK2)
C=========================================================================
C SUBROUTINE CHSP CONSTRUCTS THE ELECTRONIC CHARGE DENSITY AND COMPUTES
C THE REAL SPACE WAVEFUNCTIONS WHICH ARE USED AGAIN IN THE ELECTRON
C DYNAMICS SUBROUTINE (CHGBO.FOR)
C 
C (Now: CHDENG(r) CONTAINS THE CONTRIBUTION TO THE real-SPACE CHARGE-
C                 DENSITY FROM THIS K-POINT, (only once here)
C       CHDENR(r) CONTAINS THE SUM OF CONTRIBUTIONS TO THE REAL-SPACE 
C                 CHARGE-DENSITY FROM K-POINTS TREATED UP TILL NOW)
C=========================================================================
4115    CONTINUE
        IF (IBANS.EQ.0)
     &  CALL CHSP(NBANOC,NPLWV,NRPLWV,MPLWV,WTKPT(NKP)
     &      ,CPTWFP,CHDENG,CHDENR,NGPTAR,NINDPW(1,NKP),NPLWKP(NKP),
     &       CWORK1,CWORK2,OCC(1,NKP) )
C=========================================================================
C     WRITE WAVEFUNCTION TO EXTERNAL MEMORY
C=========================================================================
        IF (NKPTS.GT.1) THEN
          NPUNIT = 19 + NKP
          REWIND NPUNIT
          WRITE (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
        END IF
        WRITE (*,*) ' NPLWKP=',NPLWKP(NKP)
 4112   CONTINUE
C=========================================================================
C     IF RE-START, READ IN EXISTING WAVEFUNCTIONS
C=========================================================================
        IF (ISTART.EQ.1 .AND. IBANS.EQ.0) THEN 
          NPUNIT = 19 + NKP
          REWIND NPUNIT
          READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
        END IF
        IF (IVPTYP.EQ.0) GOTO  4100
        IF (NLPOT.EQ.0) THEN
C=========================================================================
C  FOR RECIPROCAL SPACE NONLOCAL POTENTIAL
C
C  INITIALISE ENVNL: NON-LOCAL CONTRIBUTION TO TOTAL ENERGY (XW 30-APR-90)
C  1) PHASGR SETS UP  CPHSGR(NRPLWV,NIONS,NSPEC) = EXP(iq.Rn)
C  2) SETVG  SETS UP  VGNL(NRPLWV,0:2,NSPEC) = Vnl(l,q)
C=========================================================================
          IF (ICLOCK.EQ.1) CALL PCLOCK(1)
          CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &              DNLG(1,1,NKP),POSION,DIRC, CPHSGR,IVPTYN)
          CALL SETVG(NRPLWV,NSPEC,NPSPTS,NIONSP,NPLWKP(NKP),
     &            DNLKG(1,0,NKP),PSPNL,PSCALE, PSGMAX, VGNL,DVGNL,
     &            NPSPTN,IVPTYN)
          DO 4121 NB = 1 , NBANDS
            CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP(NKP),
     &         DNLKG(1,0,NKP),VOLC,PSCALE,CPTWFP,CWRK20,CWRK21,
     &         CWRK22,CWRK23,CPHSGR,VGNL,CELFRC,CWOR10,CWOR11,NB,
     &         IVPTYN)
            CVNL = (0.0,0.0)
            DO 4122 M = 1 , NPLWKP(NKP)
             CVNL = CVNL + CONJG ( CPTWFP(M,NB) ) * CELFRC(M)
 4122       CONTINUE
            VNL(NB,NKP) = REAL (CVNL)
 4121     CONTINUE          
          IF (ICLOCK.EQ.1) CALL PCLOCK(2)
        ELSE
C===================================================================
C CALL SUBROUTINE TO SET UP PHASE FACTOR ARRAY FOR THIS K POINT     
C===================================================================
          CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,POSION,       
     &          NRLPPI,NIONSP,RLCORE,CPHGRD,VKPT(1,NKP),NGX,NGY,NGZ)
C                                                                   
          DO 4120 NB=1,NBANDS                                           
C===================================================================
C PREPARE ARRAYS FOR THE FOURIER TRANSFORM                          
C===================================================================
            DO 5345 M = 1 , MPLWV  
              CWORK1(M) = (0.0,0.0)
              CWORK2(M) = (0.0,0.0)
 5345       CONTINUE               
            NINDW = NRPLWV*(NB-1)                                           
            DO 222 M=1,NPLWKP(NKP)                                        
              CWORK1(NINDPW(M,NKP))=CPTWFP(M,NB)                        
 222        CONTINUE                                                      
C===================================================================
C TRANSFORM THE WAVEFUNCTION INTO REAL SPACE                          
C=====================================================================
            CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)                              
            CALL ENRLNL(VOLC,NGX,NGY,NGZ,VNL,NB,NKP,MPLWV,NRPLWV,           
     &         NRGRPT,NIONSP,NRLPPI,NBANDS,NKPTS,NSPEC,CWORK2,CWORK1,      
     &         CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,NADGRD,
     &         MXRLSH,CESAVE)              
 4120     CONTINUE                                                       
C                                                                    
        END IF
        IF (IPRINT.GE.2) WRITE (*,4125) (VNL(NB,NKP),NB=1,NBANDS)
C
 4100 CONTINUE
 4125 FORMAT(' Vnl =',6F12.5)
C=========================================================================
C
C                  END OF K-POINT-LOOP
C
C NOW:  CHDENR(r)IS THE REAL SPACE CHARGE DENSITY SUMMED OVER ALL K-POINTS.
C
C=========================================================================
C      IF (IPRINT.LT.5) GOTO 4304 
      WRITE (*,4300) 
      WRITE (90,4302) CHDENR
 4304 CONTINUE
 4300 FORMAT(/,1X,20(1H-),' CHDENR (REAL SPACE CHARGE DENSITY)',20(1H-))
 4302 FORMAT(10F8.3)
C=========================================================================
C TRANSFORM THE CHARGE DENSITY TO RECIPROCAL SPACE
C=========================================================================
      DO 5348 M = 1 , MPLWV
        CWORK1(M) = (0.0,0.0)
        CWORK2(M) = (0.0,0.0)
 5348 CONTINUE
      RINPLW = 1.0 / NPLWV
      DO 5350 M = 1 , NPLWV
        CWORK1(M) = CHDENR(M) * RINPLW
 5350 CONTINUE
      CALL FFT3D(CWORK1,CWORK2,NGPTAR,-1)
      DO 5360 M = 1 , NPLWV
        CHDENG(M) = CWORK1(M)
 5360 CONTINUE
      IF (ISYMM.EQ.1) THEN
C======================================================================
C     >>>>>>>>>>>>>>>>>>> SYMMETRIZE CHARGE DENSITY  <<<<<<<<<<<<<<<<<<
C======================================================================
        CALL ROSYM4(A1,A2,A3,B1,B2,B3,ISY,NC,IB,V,R,RB,NPLWV,
     +              NGX,NGY,NGZ,CHDENG,CWORK1,LPCTX,LPCTY,
     +              LPCTZ,LPCTXI,LPCTYI,LPCTZI,IPRINT)
C======================================================================
C    TRANSFORM THE CHARGE DENSITY TO REAL SPACE                        
C======================================================================
        DO 5949 M=1,MPLWV                                     
          CWORK1(M)=(0.0,0.0)                                   
          CWORK2(M)=(0.0,0.0)                                   
 5949   CONTINUE                                              
        DO 5951 M=1,NPLWV                                     
          CWORK1(M)=CHDENG(M)                                   
 5951   CONTINUE                                              
        CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)                    
        DO 5961 M=1,NPLWV                                     
          CHDENR(M)=CWORK1(M)                                   
 5961   CONTINUE
      END IF
C
      IF (IPRINT.GE.2)
     &  WRITE (*,*)' CHDENG(G=0)=', CHDENG(1),'   NELECT =',NELECT
C=========================================================================
C
C (AGAIN:  CHDENR(r)IS CHARGE-DENSITY IN REAL-SPACE
C          CHDENG(G) IS      -//-      IN RECIPROCAL SPACE
C
C  - READY FOR CALCULATION OF THE KOHN-SHAM-POTENTIAL
C
C=========================================================================
C     SETUP LOOK-UP TABLES FOR EX-CORR ENERGY AND FORCE
C=========================================================================
      CALL XCDAT(EXCDAT, XCFDAT, XCPDAT)
C
      IF (IION.NE.0) THEN
C=========================================================================
C     SETUP CONSTRAINTS WHICH ARE SUPPOSED TO BE LINEAR WITH
C     RESPECT TO COORDINATES. THEN INSTEAD OF THE PHYSICAL COORDINATES
C     Ri,i=1,3*N, WE INTRODUCE A SET OF AUXILIARY COORDINATES, 
C     Pi,i=1,M, M < 3*N. SUBROUTINE CONSTR SETS UP THE MATRIX OF
C     DERIVATIVES RCONSTR(i,j)=dRi/dPj.     
C
C     NCNSTR - NUMBER OF CONSTRAINTS 
C     NEWCOO = M = 3 * N - NCNSTR
C=========================================================================
        IF (ICNSTR.NE.0) THEN
          CALL CONSTR(NSPEC,NIONS,NIONST,RCNSTR,DIRC,POSIC,
     &                RMOVE,NCNSTR,NEWCOO,WORK7)
          IF (NCNSTR.EQ.0) ICNSTR = 0
        ELSE
          NEWCOO = 3 * NIONST
          NCNSTR = 0
        END IF
      END IF
C=========================================================================
C
C ++++++++++++ DO 1000 N=1,REQUIRED NUMBER OF TIMESTEPS ++++++++++++++++
C
C THIS IS THE MAIN LOOP OF THE PROGRAM DURING WHICH ONE COMPLETE STEP OF
C THE ELECTRON DYNAMICS IS PERFORMED. AT THE END OF THE REQUIRED NUMBER
C OF STEPS THE ELECTRONIC WAVEFUNCTIONS, POSITIONS OF THE IONS ETC ARE
C WRITTEN TO AN EXTERNAL FILE
C=========================================================================
C
C AFTER NDELAY ITERATIONS THERE IS AN OPTION OF PERFORMING ELECTRONIC
C MINIMISATION MORE OFTEN THAN IONIC (IMDTR) OR UNIT CELL (BMDTR) 
C RELAXATION. THIS CAN BE DONE BY SETTING PARAMETER NITFIX>1 WHICH
C WILL MEAN THAT YOU PERFORM NITFIX ELECTRONIC ITERATIONS AT FIXED
C GEOMETRY. IF NITFIX IS NOT SPECIFIED, IT WILL BE SET EQUAL TO 1,
C AS IN PREVIOUS VERSION OF CGMAIN.
C=========================================================================
C
      NITFXX = NITFIX
      IF (NITFXX.EQ.0 .OR. IBANS.NE.0) NITFXX = 1
      N = 0
 1000 CONTINUE
      N = N + 1
C==================================================================
C  THIS SERVES TO CONTROL THE RUN FROM OUTSIDE                     
C==================================================================
      OPEN(UNIT=1,FILE='runinf',STATUS='OLD',ERR=1101)             
      READ(1,*,ERR=1101) NITER                                     
      READ(1,*,ERR=1101) IION                                      
      READ(1,*,ERR=1101) POTIM                                     
      READ(1,*,ERR=1101) NDELAY                                    
1101  CONTINUE                                                     
      CLOSE(1)                                                     
C
        WRITE (*,90) N
        NITFI = 0
C
        IF (IBANS.NE.0 .AND. N.GT.1) GO TO 4219
C
        IF (ICLOCK.EQ.1) CALL PCLOCK(1000)
        IF (IPRINT.EQ.5) WRITE (*,91) CHDENG
  90  FORMAT(//,1X,20(1H>),' ITERATION =',I5,1X,20(1H<),/)
  91  FORMAT(1X,9F8.4)
C========================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE FEWALD }}}}}}}}}}}}}}}}}}}}}}}}}}}}
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 FORCE ON THE UNIT CELL THAT IS RELATED TO THE CHANGE IN THE EWALD
C ENERGY ON CHANGING THE SIZE OF THE CELL ARE ALSO CALCULATED.
C========================================================================
        NINDX = 1
        DO 3190 NSP = 1 , NSPEC
          DO 3191 NI = 1 , NIONSP(NSP)
            DO 3192 M = 1 , 3
              POSIC(M,NINDX) = POSION(M,NI,NSP)
 3192       CONTINUE
          NINDX = NINDX + 1
 3191     CONTINUE
 3190   CONTINUE
        IF ( N.EQ.1 .OR. 
     &     ( N.GT.NDELAY .AND. IBOX.EQ.1 .AND. IION.EQ.0 ) )
     &    CALL EWALTR(POSIC,EWIFC,EWRLEN,EWRCSS,EWRLSS,EWRCSI,DIRC,RECC,
     &         VOLC,SIGEW,TEWEN,NIONST,ICHARC,NIONCH,NICHSQ,CPHFX,CPHFY,
     &         CPHFZ,IPRINT,NIONST,NEWPTS,MAXCX,MAXCY,MAXCZ,MAXGPX,
     &         MAXGPY,MAXGPZ,NGPTS,RFORCE,CFORCE,FORCEG,FORSIG,ENERG)
        NINDX = 1
        DO 3193 NSP = 1 , NSPEC
          DO 3194 NI = 1 , NIONSP(NSP)
            DO 3195 M = 1 , 3
              EWIFOR(M,NI,NSP) = EWIFC(M,NINDX)
 3195       CONTINUE
            NINDX = NINDX + 1
            IF (IPRINT.GE.2) WRITE (*,49) (EWIFOR(M,NI,NSP),M=1,3)
 3194     CONTINUE
 3193   CONTINUE
   49 FORMAT('  EWALD FORCE ON ION :',3F15.6)
C=========================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE FDIR16 }}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE CALCULATES THE HARTREE POTENTIAL FROM THE ELECTRONIC
C CHARGE DENSITY CALCULATED IN SUBROUTINE CHSP. THE CORRECTION TO THE
C TOTAL ENERGY DUE TO OVERCOUNTING THE HARTREE ENERGY ON SUMMING THE
C ELECTRONIC EIGENVALUES IS ALSO COMPUTED (HARTREE CONTRIBUTION TO THE
C TOTAL ENERGY = 0.5*SUM (VH(G)*RHO(-G)) WHEREAS SUM OF EIGENVALUES
C GIVES SUM (VH(G)*RHO(-G)) WHERE VH(G) IS THE HARTREE POTENTIAL AT
C WAVEVECTOR G AND RHO(G) IS THE CHARGE DENSITY AT WAVEVECTOR G
C=========================================================================
 4801   CONTINUE
        NITFI = NITFI + 1
        CALL FDIRTR(NGX,NGY,NGZ,NPLWV,RECC,VOLC,CVD,CHDENG,
     &              SIGHA,DENC,LPCTX,LPCTY,LPCTZ,DIRDAT)
C=========================================================================
C ( CVD(G) IS NOW THE HARTREE-POTENTIAL IN RECIPROCAL SPACE)
C ( DIRDAT(G) IS THE 1/G**2 FACTORS FOR RECALC. OF CVD(G) FOR
C             THE SAME UNITCELL)
C=========================================================================
        IF (IPRINT.GE.2) WRITE (*,*)'HARTREE POT. IN RECIPROCAL SPACE'
        IF (IPRINT.GE.5) WRITE (*,47) CVD
  47  FORMAT(1X,10F8.4)
C=========================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE FEXCTR }}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE CALCULATES THE EXCHANGE CORRELATION POTENTIAL FROM THE
C REAL SPACE CHARGE DENSITY. THE CORRECTION TO THE TOTAL ENERGY DUE TO 
C OVERCOUNTING THE EXCHANGE CORRELATION ENERGY ON SUMMING THE ELECTRONIC
C EIGENVALUES IS ALSO COMPUTED.
C
C  FEXCTR: AT INPUT CV(r) IS THE REAL-SPACE CHARGE DENSITY 
C          AT OUPUT CV(r) IS THE REAL-SPACE XC-POTENTIAL
C
C=========================================================================
        DO 9111 M = 1 , NPLWV
          CV(M) = CHDENR(M)
 9111   CONTINUE
        CALL FEXCTR(NPLWV,CV,VOLC,SIGXC,XCENC,XCENER,EXCDAT,
     &              XCFDAT,XCPDAT)
        IF (IPRINT.GE.2) WRITE (*,*)' EX-CORR POTENTIAL IN REAL SPACE'
        IF (IPRINT.GE.5) WRITE (*,47) CV
C=========================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE FSTFTR }}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE CALCULATES THE STRUCTURE FACTOR ON THE GRID OF
C RECIPROCAL LATTICE VECTORS
C=========================================================================
        DO 3110 NSP = 1 , NSPEC
          CALL FSTFTR(NGX,NGY,NGZ,NIONSP(NSP),NPLWV,POSION(1,1,NSP),
     &                CSTRF(1,NSP),LPCTX,LPCTY,LPCTZ)
 3110   CONTINUE
        IF (IPRINT.GE.2) WRITE (*,*) ' STRUCTURE FACTOR'
        IF (IPRINT.GE.5) WRITE (*,47) CSTRF  
C=========================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE FVPTR }}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE SETS UP THE PSEUDOPOTENTIAL ON THE GRID OF RECIPROCAL
C LATTICE VECTORS. THIS SUBROUTINE MUST BE CALLED WHENEVER RECC(i,j),
C I.E. THE UNIT CELL, HAS BEEN CHANGED
C
C     calc.   V  (G)
C              alpha
C
C (NOW WORKING IN RECIPROCAL SPACE)
C
C=========================================================================
        IF (N.EQ.1 .OR. (N.GE.NDELAY.AND.IBOX.EQ.1)) THEN
          PSCENC = 0.0
          DO 3120 NSP = 1 , NSPEC
            CALL FVPTR(NGX,NGY,NGZ,NPSPTN(NSP),RECC,VOLC,PSP(1,NSP),
     &                 PSGMAX(NSP),VPS(1,NSP),DVPS(1,NSP),LPCTX,LPCTY,
     &                 LPCTZ,ICHARG(NSP),IVPTYN(NSP))
 3120     CONTINUE
        END IF
C=========================================================================
C
C CALCULATE THE CONTRIBUTION TO THE TOTAL ENERGY FROM THE NON-COULOMB
C PART OF THE G=0 COMPONENT OF THE PSEUDOPOTENTIALS AND THE FORCE ON THE
C UNIT CELL DUE TO THE CHANGE IN THIS ENERGY AS THE SIZE OF THE CELL
C CHANGES
C
C               N
C                elec   --
C     E       = -----   >   N      * v
C      core      vol    --   alpha    alpha,core
C                      alpha
C
C=========================================================================
        PSCENC = 0.0
        DO 3125 NSP = 1 , NSPEC
          PSCENC = PSCENC + PSCORE(NSP) * NIONSP(NSP)
 3125   CONTINUE
        PSCENC = PSCENC * NELECT / VOLC
        IF (IPRINT.GE.2) WRITE (*,*)' PSCENC=',PSCENC
C=========================================================================
C    STRESSES ON THE UNIT CELL DUE TO THE ALPHA*Z TERM
C
C                                --  
C     Sigma      = - delta       >    N      * v
C          ij             ij     --    alpha    alpha,core
C                               alpha
C
C     alpha loops over all species.  i,j = x,y,z
C=========================================================================
        SIGAL(1) = - PSCENC / VOLC
        SIGAL(2) = - PSCENC / VOLC
        SIGAL(3) = - PSCENC / VOLC
        SIGAL(4) = 0.0
        SIGAL(5) = 0.0
        SIGAL(6) = 0.0
C=========================================================================
C CALCULATE THE TOTAL POTENTIAL WHICH IS THE SUM OF THE HARTREE
C POTENTIAL, THE EXCHANGE-CORRELATION POTENTIAL AND THE STRUCTURE
C FACTOR*THE PSEUDOPOTENTIAL
C
C STEP 1.  CALCULATE THE TOTAL IONIC POTENTIAL IN RECIPROCAL SPACE
C
C                --
C     V   (G) =  >   S  (G) * v  (G)
C      ion       --   alpha    alpha
C               alpha
C
C (WORKING IN RECIPROCAL SPACE)
C
C=========================================================================
        DO 3126 NP = 1 , MPLWV
          CWORK1(NP) = (0.0,0.0)
          CWORK2(NP) = (0.0,0.0)
 3126   CONTINUE
        DO 3127 NSP = 1 , NSPEC
          DO 3128 NP = 1 , NPLWV
            CWORK1(NP) = CWORK1(NP) + VPS(NP,NSP) * CSTRF(NP,NSP)
 3128     CONTINUE
 3127   CONTINUE
C=========================================================================
C STEP 2. ADD THE HARTREE POTENTIAL AND TRANSFORM TO REAL SPACE
C=========================================================================
        DO 3129 NP = 1 , NPLWV
          CVD(NP) = CWORK1(NP) + CVD(NP)
 3129   CONTINUE
C=========================================================================
C CVD(g) = CONTAINS THE RECIPROCAL SPACE IONIC AND HARTREE POT)
C CVION(g) = CWORK1 (g) = SUM [VPS*CSTRF]
C
C=========================================================================
C BOTH OF THEM ARE TRANSFORMED IN TO REAL SPACE
C=========================================================================
        CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)
        DO 3124 NP = 1 , NPLWV
          CVION(NP) = CWORK1(NP)
 3124   CONTINUE
        DO 5135 M = 1 , NPLWV
          CWORK1(M) = CVD(M)
 5135   CONTINUE
        CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)
C=========================================================================
C STEP 4. ADD THE EXCHANGE-CORRELATION POTENTIAL IN REAL SPACE
C         (REMEMBER CV(r) IS XC-POTENTIAL)
C=========================================================================
        DO 5136 NP = 1 , NPLWV
          CV(NP) = CWORK1(NP) + CV(NP)
 5136   CONTINUE
C=========================================================================
C
C    CV(r) IS (HARTREE + XC + ION) POTENTIAL FOR THE Kohn-Sham-EQUATION 
C CVION(r) IS THE IONIC POTENTIAL ALONE
C
C=========================================================================
C INIT THE KINETIC-ENERGY CONTRIBUTION TO THE STRESS ON THE UNIT CELL
C=========================================================================
        DO 4210 M = 1 , 6
          SIGKE(M) = 0.0
 4210   CONTINUE
C
        IF((N.GT.NDELAY).AND.(IBOX.EQ.1 .OR. IION.EQ.1)                     
     &                  .AND.(IVPTYP.EQ.1 .AND. NLPOT.EQ.1)) THEN
C=======================================================================
C CALL SUBROUTINE TO SET UP DATA ARRAYS FOR REAL SPACE NON-LOCAL
C PSEUDOPOTENTIALS
C=======================================================================
          CALL 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)
        ENDIF                                                             
C=========================================================================
C K-POINT LOOP STARTS
C=========================================================================
 4219   CONTINUE
        DO 4220 NKP = 1 , NKPTS
C=========================================================================
C     READ WAVEFUNCTION FROM EXTERNAL MEMORY
C=========================================================================
          NPUNIT = 19 + NKP
          NUNIT  = 39 + NKP
          IF (NKPTS.GT.1) THEN
            REWIND NPUNIT
            READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
          END IF
C=========================================================================
C  CPHSGR AND VGNL ARE STORED, SO WE HAVE TO RE-CALCULATE THEM EVERY TIME.
C=========================================================================
          IF (IVPTYP.EQ.0) GO TO 4222
          IF (NLPOT.EQ.0) THEN
            IF (NKPTS.GT.1 .OR. N.GE.NDELAY) 
     &      CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &                  DNLG(1,1,NKP),POSION,DIRC, CPHSGR,IVPTYN)
            IF (NKPTS.EQ.1 .AND. IBOX.EQ.0) GO TO 4222
            CALL SETVG(NRPLWV,NSPEC,NPSPTS,NIONSP,NPLWKP(NKP),
     &           DNLKG(1,0,NKP),PSPNL,PSCALE, PSGMAX, VGNL,DVGNL,
     &           NPSPTN,IVPTYN)
          ELSE
c===================================================================
C CALL SUBROUTINE TO SET UP PHASE FACTOR ARRAY FOR THIS K POINT     
C===================================================================
            IF (NKPTS.GT.1 .OR. N.GE.NDELAY)
     &       CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,POSION,       
     &          NRLPPI,NIONSP,RLCORE,CPHGRD,VKPT(1,NKP),NGX,NGY,NGZ)
          END IF
 4222     CONTINUE
C
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE CONGRA }}}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE PERFORMS THE ENERGY MINIMISATION OF THE ELECTRONIC 
C SYSTEM USING THE CONJUGATE GRADIENT METHOD.
C THE SUBROUTINE ALSO CALCULATES THE ELECTRONIC EIGENVALUES AND
C THE STRESS ON THE UNIT CELL RESULTING FROM THE CHANGE IN THE KINETIC
C ENERGY OF THE PLANE WAVE BASIS STATES AS THE SIZE OF THE CELL CHANGES
C
C=========================================================================
         IF (IBANS.EQ.0) THEN
           IF (N.GT.4) THEN   
             SUMW = WTKPT(NKP)
           ELSE               
             SUMW = SUMWEI    
           END IF             
           CALL CONGRA(NGX,NGY,NGZ,NBANOC,NKPTS,NPLWV,MPLWV,
     &       NRPLWV,ENMAX,NINDPW(1,NKP),NPLWKP(NKP),WTKPT,
     &       CV,CPTWFP,CPTWFL,DIRC,RECC,VOLC,CELEN,NGPTAR,
     &       SUMW,DATAKE(1,1,NKP),CWORK1,
     &       CWORK2,CWORK,NBANOC, CWORK6,
     &       CWORK7,CWORK8,CWORK9,CWOR10, WORK12,CVION,CVD,CHDENR,
     &       DENC,XCENC,TEWEN,PSCENC,NKP,CHDENG,
     &       XCFDAT,XCPDAT,DIRDAT,CWOR11,
     &       HR, HI, AUX, FV1, FV2, FV3, CH0, NITMAX,TOTEN1,
     &       NSPEC, NIONS,NIONSP,PSCALE,
     &       DNLKG, CPHSGR,VGNL, CELFRC, CWRK20,
     &       CWRK21, CWRK22, CWRK23,VNL, IVPTYP, IPRINT,ISBROT,
     &       IOCCUP,OCC,EIGVAL, ISYMM,NLPOT,
     &       SIGKE, SIGXC,EXCDAT,ICLOCK,IVPTYN,NRGRPT,NRLPPI,
     &       CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,
     &       NADGRD,MXRLSH,CESAVE)
         ELSE
           CALL CONGRABS(NBANOC,NKPTS,NPLWV,MPLWV,
     &       NRPLWV,NINDPW(1,NKP),NPLWKP(NKP),WTKPT,
     &       CV,CPTWFP,CPTWFL,DIRC,RECC,VOLC,CELEN,NGPTAR,
     &       DATAKE(1,1,NKP),CWORK1,CWORK2, CWORK, CWORK6,CWORK7,
     &       CWORK8,CWORK9,CWOR10, WORK12,CVD,NKP,CWOR11,
     &       HR, HI, AUX, FV1, FV2, FV3, CH0, NITMAX,
     &       NSPEC, NIONS, NIONSP,PSCALE,
     &       DNLKG, CPHSGR,VGNL, CELFRC, CWRK20,
     &       CWRK21, CWRK22, CWRK23,VNL, IVPTYP, IPRINT,ISBROT,
     &       OCC,EIGVAL,NLPOT,
     &       ICLOCK,IVPTYN,NITER,N,NGX,NGY,NGZ,NRGRPT,NRLPPI, 
     &       CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,  
     &       NADGRD,MXRLSH,CESAVE)
         END IF
C=========================================================================
C     WRITE WAVEFUNCTION TO EXTERNAL MEMORY
C=========================================================================
          IF (NKPTS.GT.1) THEN
            REWIND NPUNIT
            WRITE (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
          END IF
C======================================================================
C     WAVEWR PICKS OUT THE NL(=6) MOST SIGNIFICANT PLANE WAVE
C     COMPONENTS OF EACH BAND. 
C======================================================================
C      IF (N.EQ.NITER.OR.N.EQ.NDELAY)
C     & CALL WAVEWR(NBANDS, NRPLWV, NPLWKP, NGX, NGY, NGZ,
C     & SIZEX,SIZEY,SIZEZ,CPTWFP,LPCTX,LPCTY,LPCTZ,NINDPW(1,NKP),NKP)
C
 4220   CONTINUE
C======================================================================
C     >>>>>>>>>>>>>>>>>>> K-POINT LOOP ENDS HERE <<<<<<<<<<<<<<<<<<
C======================================================================
        IF (IBANS.EQ.1) GO TO 4218
C
C======================================================================
C     UPDATE OCCUPATION NUMBER OCC ACCORDING TO THE EIGEN-ENERGY OF 
C     EACH BAND. (EF IS THE FERMI ENERGY DETERMINED)
C======================================================================
        IF (IOCCUP.EQ.0) THEN
          IF ((NELECT/2)*2.EQ.NELECT) THEN 
            GO TO 4800
          ELSE
            GO TO 4802
          END IF
        END IF
C
        IDEL = N / NDEL
        DEL = AMAX1 ( DELMAX / (2**IDEL) , DELMIN ) 
        WRITE (*,*) ' '
        WRITE (*,*) 'THE GAUSSIAN BROADENING ENERGY DEL (eV) = ', DEL
C======================================================================
C     TAKE OUT THE SPIN-DEGENERACY, FOR IT IS CONSIDERED EXPLICITLY
C     IN THIS PROGRAM. BUT NOT IN EFERMI.
C======================================================================
C
C      IN THE MAIN PROGRAM AND ALL OTHER HAMMER'S SUBROUTINE,
C      OCC INCLUDES SPIN-DEGENERACY (2.0), EXCEPT IN C2.F
C
        IF (NBANDS*2.GT.NELECT) THEN 
          CALL EFERMI(NELECT,NBANDS,DEL,NKPTS,NBANDS,NKPTS,WTKPT,
     &                OCC,EF,EIGVAL,SORT) 
          DO 4500 JKP = 1 , NKPTS
            DO 4500 JB = 1 , NBANDS
              OCC(JB,JKP) = OCC(JB,JKP) / 2.0
 4500     CONTINUE
        END IF
C======================================================================
C     CALCULATE THE NEW CHARGE DENSITY ACCORDING TO THE NEW SET OF
C     OCCUPATION NUMBER 
C     (NO NEED TO UPDATE CELEN AND VNL, BECAUSE THEY ARE NOT RELATED TO 
C      CHARGE DENSITY, THEY ARE RELATED TO WAVE FUNCTIONS ONLY.)
C======================================================================
 4802   CONTINUE
        DO 4600 M = 1 , NPLWV
          CHDENR(M) = (0.0,0.0)
          CHDENG(M) = (0.0,0.0)
 4600   CONTINUE
        DO 4602 NKP = 1 , NKPTS
          IF (NKPTS.GT.1) THEN 
            NPUNIT = 19 + NKP
            REWIND NPUNIT
            READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
          END IF
          CALL CHSP(NBANOC,NPLWV,NRPLWV,MPLWV,
     &              WTKPT(NKP),CPTWFP,CHDENG,CHDENR,NGPTAR,
     &              NINDPW(1,NKP),NPLWKP(NKP),CWORK1,CWORK2,
     &              OCC(1,NKP) )
 4602   CONTINUE
C=======================================================================
C FFT (CHDENR) TO GET THE RECIPROCAL SPACE CHARGE DENSITY: CHDENG(g)
C=======================================================================
        DO 4700 M = 1 , MPLWV
          CWORK1(M) = (0.0,0.0)
 4700   CONTINUE
        RINPLW = 1.0 / NPLWV
        DO 4702 M = 1 , NPLWV
          CWORK1(M) = CHDENR(M) * RINPLW
 4702   CONTINUE
        CALL FFT3D(CWORK1,CWORK2,NGPTAR,-1)
        DO 4704 M = 1 , NPLWV
          CHDENG(M) = CWORK1(M)
 4704   CONTINUE 
C=========================================================================
C NOW: CHDENR(r) CONTAINS THE NEW CHARGE DENSITY IN   R E A L  SPACE
C      CHDENG(g) CONTAINS THE NEW CHARGE DENSITY IN RECIPROCAL SPACE
C    -----  Change over for IOCCUP=1
C=========================================================================
C
 4800   CONTINUE
C
C======================================================================
C     >>>>>>>>>>>>>>>>>>> SYMMETRIZE CHARGE DENSITY  <<<<<<<<<<<<<<<<<<
C======================================================================
C
        IF (ISYMM.EQ.1 .AND. N.GE.1) THEN
          IF (IOCCUP.EQ.0) THEN
C=======================================================================
C FFT (CHDENR) TO GET THE RECIPROCAL SPACE CHARGE DENSITY: CHDENG(g)    
C SKIPPED FOR PARTIAL OCCUPANCIES - DONE ALREADY
C=======================================================================
        DO 4670 M = 1 , NPLWV
          CHDENR(M) = (0.0,0.0)
          CHDENG(M) = (0.0,0.0)
 4670   CONTINUE
        DO 4672 NKP = 1 , NKPTS
          IF (NKPTS.GT.1) THEN 
            NPUNIT = 19 + NKP
            REWIND NPUNIT
            READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
          END IF
          CALL CHSP(NBANOC,NPLWV,NRPLWV,MPLWV,
     &              WTKPT(NKP),CPTWFP,CHDENG,CHDENR,NGPTAR,
     &              NINDPW(1,NKP),NPLWKP(NKP),CWORK1,CWORK2,
     &              OCC(1,NKP) )
 4672   CONTINUE
c
            DO 4901 M=1,MPLWV                                      
              CWORK1(M)=(0.0,0.0)                                    
 4901       CONTINUE                                               
            RINPLW=1.0/NPLWV                                       
            DO 4902 M=1,NPLWV                                      
              CWORK1(M)=CHDENR(M)*RINPLW                             
 4902       CONTINUE                                               
            CALL FFT3D(CWORK1,CWORK2,NGPTAR,-1)                    
            DO 4904 M=1,NPLWV                                      
              CHDENG(M)=CWORK1(M)
 4904       CONTINUE          
          END IF
C
          CALL ROSYM4(A1,A2,A3,B1,B2,B3,ISY,NC,IB,V,R,RB,NPLWV,
     +                NGX,NGY,NGZ,CHDENG,CWORK1,LPCTX,LPCTY,
     +                LPCTZ,LPCTXI,LPCTYI,LPCTZI,IPRINT)
C======================================================================
C    TRANSFORM THE CHARGE DENSITY TO REAL SPACE                        
C======================================================================
          DO 5948 M=1,MPLWV                                     
            CWORK1(M)=(0.0,0.0)                                   
            CWORK2(M)=(0.0,0.0)                                   
 5948     CONTINUE                                              
          DO 5950 M=1,NPLWV                                     
            CWORK1(M)=CHDENG(M)                                   
 5950     CONTINUE                                              
          CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)                    
          DO 5960 M=1,NPLWV                                     
            CHDENR(M)=CWORK1(M)                                   
 5960     CONTINUE
        END IF
C
        IF (IPRINT.GE.2) WRITE (*,*)' CHDENG(G=0) =', CHDENG(1)
        IF (NITFI.LT.NITFXX .AND. N.GT.NDELAY) GO TO 4801
C=========================================================================
C RE-CALCULATE THE TOTAL ENERGY WHEN REQUIRED. 
C=========================================================================
C       IF (ISBROT.EQ.0) GOTO 4900
C      IF (ISBROT.EQ.1 .AND. NBANDS*2.EQ.NELECT) GOTO 4900
C======================================================================
C 1. CALC. HARTREE POTENTIAL AND ENERGY CORRECTION IN RECIPROCAL SPACE
C======================================================================
        DENC = 0.0
        DO 4810 M = 1 , NPLWV
          CVD(M) = CHDENG(M) * DIRDAT(M)
          DENC = DENC + REAL ( CVD(M) * CONJG ( CHDENG(M) ) )
 4810   CONTINUE
        DENC = - DENC / 2.0
C======================================================================
C TRANSFORM HARTREE POTENTIAL INTO REAL SPACE
C======================================================================
        DO 4820 M = 1 , MPLWV
          CWORK1(M) = (0.0,0.0)
 4820   CONTINUE
        DO 4822 M = 1 , NPLWV
          CWORK1(M) = CVD(M)
 4822   CONTINUE
        CALL FFT3D(CWORK1,CWORK,NGPTAR,1)
        DO 4824 M = 1 , NPLWV
          CVD(M) = CWORK1(M)
 4824   CONTINUE
C======================================================================
C 2. CALCULATE EXCHANGE CORRELATION ENERGY CORRECTION USING CHDENR(r)
C======================================================================
        DO 4830 M = 1 , NPLWV
          CV(M) = CHDENR(M)
 4830   CONTINUE
        CALL FEXCTR(NPLWV,CV,VOLC,SIGXC,XCENC,XCENER,EXCDAT,XCFDAT,
     &              XCPDAT)
        DO 4835 M = 1 , NPLWV
          CV(M) = CV(M) + CVION(M)
          CV(M) = CV(M) + CVD(M)
 4835   CONTINUE
C======================================================================
C CALCULATE THE TOTAL ENERGY (IN REAL SPACE)
C======================================================================
        ENPOT = 0.0
        DO 4840 M = 1 , NPLWV
          ENPOT = ENPOT + REAL ( CHDENR(M) * CV(M) )
 4840   CONTINUE
        ENPOT = ENPOT / NPLWV
C======================================================================
C THE NON-LOCAL PSEUDOPOTENTIAL ENERGY FOR EACH BAND HAS BEEN UPDATED 
C AT THE END OF SUBROUTINE CONGRA
C======================================================================
        ENKE = 0.0
        ENVNL = 0.0
        DO 4850 NK = 1 , NKPTS
          DO 4850 NB = 1 , NBANDS
            ENKE = ENKE + 2.0 * WTKPT(NK) * 
     &             REAL ( CELEN(NB,NK) ) * OCC(NB,NK)
            IF (IVPTYP.EQ.1) 
     &         ENVNL = ENVNL + 2.0 * WTKPT(NK) * VNL(NB,NK) * OCC(NB,NK)
 4850   CONTINUE
 4900   CONTINUE
C=========================================================================
C WRITE OUT THE ENERGY EIGENVALUES OF THE ELECTRONIC STATES
C=========================================================================
        IF (ISBROT.EQ.0) GO TO 304 
        IF (IPRINT.LT.1 .AND. N.NE.NITER) GO TO 304
 4218   CONTINUE
        DO 302 NK = 1 , NKPTS
          WRITE (*,4217) NK,(VKPT(K,NK),K=1,3),WTKPT(NK)
          DO 301 NB = 1 , NBANDS
            WRITE (*,300) NB,EIGVAL(NB,NK),REAL(CELEN(NB,NK)),OCC(NB,NK)
 301      CONTINUE
 302    CONTINUE
        IF (IBANS.EQ.1) GO TO 1001
 304    CONTINUE
 300  FORMAT('  BAND(',I3,'):  E=',F11.6,'  KE=',F10.6,'  OCC=',F10.6)
 4217 FORMAT(/' K-POINT (',I3,')   (',3F10.6,')  WEIGHT=',F10.6)
C=========================================================================
C      SUM THE CONTRIBUTIONS TO THE TOTAL ENERGY OF THE SYSTEM
C=========================================================================
        TOTEN = ENKE + DENC + XCENC + TEWEN + PSCENC + ENPOT + ENVNL
        TOTEN0 = TOTEN
C=========================================================================
C
C             SELF EXPLANATORY WRITE STATEMENTS
C
C=========================================================================
        WRITE (*,202) N
        WRITE (*,2600) 'TOTAL KINETIC ENERGY        ',  ENKE
        WRITE (*,2600) 'LOCAL POTENTIAL ENERGY      ',  ENPOT
        IF (IVPTYP.NE.0)
     &  WRITE (*,2600) 'NONLOCAL POTENTIAL  ENERGY  ',  ENVNL
        WRITE (*,2600) 'HARTREE ENERGY CORRECTION   ',  DENC
        WRITE (*,2600) 'EX-CORR ENERGY CORRECTION   ',  XCENC
        WRITE (*,2600) 'EWALD ENERGY                ',  TEWEN
        WRITE (*,2600) 'PSEUDOPOTENTIAL CORE ENERGY ',  PSCENC
        WRITE (*,2600) 'TOTAL ENERGY IS             ',  TOTEN
        IF (ICLOCK.EQ.1) CALL PCLOCK(0)
2600  FORMAT(1X,A28,F15.7)
 202  FORMAT(/,1X,'ITERATION',I4)
C=====================================================================
C
C   NOW WE CALCULATE FORCES ON IONS AND STRESSES ON THE UNIT CELL
C
C THE FOLLOWING IF STATEMENT SKIPS THE SECTION OF THE PROGRAM WHICH
C CALLS THE SUBROUTINES THAT PERFORM THE DYNAMICS FOR THE IONIC SYSTEM
C AND THE RELAXATION OF THE SIZE OF THE UNIT CELL. THE HELLMANN-FEYNMAN
C FORCES ARE NOT CORRECT UNTIL THE ELECTRONIC SYSTEM HAS RELAXED TO ITS
C GROUNDSTATE SO THESE SUBROUTINE CALLS SHOULD BE SKIPPED UNTIL THE
C ELECTRONIC SYSTEM IS IN ITS GROUNDSTATE
C=========================================================================
        IF (N.LE.NDELAY .AND. N.NE.(N/NPRINT)*NPRINT) GO TO 8000
C
        DO 3130 NSP = 1 , NSPEC
          CALL FELITR(NGX,NGY,NGZ,NIONSP(NSP),POSION(1,1,NSP),
     &                CHDENG,RECC,VPS(1,NSP),EIFOR(1,1,NSP),LPCTX,LPCTY,
     &                LPCTZ,LPCTFX,LPCTFY,LPCTFZ,CWORK3,CWORK4,CWORK5)
C=========================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE FEICTR }}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE CALCULATES THE FORCE ON THE UNIT CELL WHICH RESULTS
C FROM THE CHANGE IN THE ELECTRON-ION ENERGY ON CHANGING THE SIZE OF
C THE CELL
C=========================================================================
          IF (IBOX.NE.0 .OR. N.EQ.NITER)
     &    CALL FEICTR(NGX,NGY,NGZ,NPLWV,CSTRF,RECC,VOLC,DVPS,VPS,SIGEI,
     &       CHDENG,LPCTX,LPCTY,LPCTZ,NSPEC,CWORK1,NIONSP,POSION,NIONS)
          DO 3135 MU = 1 , NIONS
            DO 3134 M = 1 , 3
              FNLEIF(M,MU,NSP) = 0.0
              FNLRL(M,MU) = 0.0
 3134       CONTINUE
 3135     CONTINUE
 3130   CONTINUE
C=========================================================================
C NEXT, WE CALCULATE FORCES AND STRESSES FROM THE NON-LOCAL PART OF 
C PSEUDOPOTENTIAL, WHICH IS GIVEN IN THE KLEINMAN-BYLANDER FORM.
C=========================================================================
        IF (IVPTYP.NE.0) THEN
C=========================================================================
C INITIALISE THE NON-LOCAL FORCE AND STRESS ARRAYS
C=========================================================================
          DO 3136 M = 1 , 6
            SIGNL(M) = 0.0
 3136     CONTINUE
C
          DO 3131 NKP = 1 , NKPTS
C
            IF (NKPTS.GT.1) THEN 
              NPUNIT = 19 + NKP
              REWIND NPUNIT
              READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
            END IF   
            IF (NLPOT.EQ.0) THEN
              CALL SETVG(NRPLWV,NSPEC,NPSPTS,NIONSP,NPLWKP(NKP),
     &          DNLKG(1,0,NKP),PSPNL,PSCALE, PSGMAX, VGNL,DVGNL,
     &          NPSPTN,IVPTYN)
              CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &                  DNLG(1,1,NKP),POSION,DIRC, CPHSGR,IVPTYN)
C=========================================================================
C FNLFOR CALCULATES THE FORCES
C=========================================================================
              DO 3137 NSP = 1 , NSPEC
                IF (IVPTYN(NSP).EQ.0) GO TO 3137
                CALL FNLFOR(NBANDS,NRPLWV,NIONSP(NSP),NPLWKP(NKP),NIONS,
     &                  VOLC,PSCALE(0,NSP),CPHSGR(1,1,NSP),
     &                  VGNL(1,0,NSP),DNLG(1,1,NKP),DNLKG(1,0,NKP),
     &                  CPTWFP,FORWK,OCC(1,NKP),CWOR10,CWOR11 )
                DO 3133 MU = 1 , NIONSP(NSP)
                  DO 3132 M = 1 , 3
                    FNLEIF(M,MU,NSP) = FNLEIF(M,MU,NSP) + 
     &                                 FORWK(M,MU) * WTKPT(NKP)
 3132             CONTINUE
 3133           CONTINUE
 3137         CONTINUE
C=========================================================================
C FSIGNL CALCULATES THE STRESSES
C=========================================================================
C     IF (ISTRSS.NE.1) GOTO 3131
              IF (IBOX.EQ.0 .AND. N.NE.NITER) GO TO 3131
              CALL FSIGNL(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP(NKP),
     &                  DNLKG(1,0,NKP),VOLC,PSCALE,CPTWFP,CPHSGR,
     &                  VGNL,DVGNL,NBANDS,VNL(1,NKP),OCC(1,NKP),SIGWK,
     &                  CWRK30,CWRK31,CWRK32,CWRK33,CWRK34,CWRK35,
     &                  CWRK36,CWRK37,CWRK38,CWRK39,IVPTYN)
              DO 3138 M = 1 , 6
                SIGNL(M) = SIGNL(M) + SIGWK(M) * WTKPT(NKP)
 3138         CONTINUE
            ELSE
C========================================================================
C CALCULATE THE FORCES ON THE IONS FROM THE NON-LOCAL POTENTIAL USING
C THE REAL SPACE METHOD
C========================================================================
              IF (NKPTS.GT.1) 
     &          CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,POSION,
     &             NRLPPI,NIONSP,RLCORE,CPHGRD,VKPT(1,NKP),NGX,NGY,NGZ)
              CALL FORLNL(VOLC,NGX,NGY,NGZ,CPTWFP,NKP,MPLWV,
     &             NRPLWV,NRGRPT,NIONSP,NRLPPI,NBANDS,NKPTS,NSPEC,
     &             CWORK1,CWORK2,CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,
     &             PRLSCA,VRLGRD,NADGRD,MXRLSH,
     &             FNLRL,DVRLGR,NINDPW,NPLWKP,OCC,IOCCUP,NGPTAR,WTKPT,
     &             CESAVE)
C========================================================================
C FSIGNL CALCULATES THE STRESSES
C========================================================================
            END IF
C
 3131     CONTINUE
          IF (NLPOT.EQ.1) THEN
            WRITE(*,*)' NON-LOCAL STRESS CONTRIBUTION NOT PROGRAMMED'
              NG = 1                              
              DO 3237 NSP = 1 , NSPEC             
                DO 3233 MU = 1 , NIONSP(NSP)      
                  DO 3232 M = 1 , 3               
                    FNLEIF(M,MU,NSP) = FNLRL(M,NG)
 3232             CONTINUE                        
                  NG = NG + 1                     
 3233           CONTINUE                          
 3237         CONTINUE                            
          END IF
        END IF
C=====================================================================
C SYMMETRIZE THE STRESS-TENSOR, I.E. SIGKE,SIGXC,SIGHA,SIGEI,SIGNL    
C=====================================================================
       IF (ISYMM.EQ.1 .AND. (IBOX.EQ.1 .OR. N.EQ.NITER)) THEN     
         CALL STRSYM(NC,IB,R,SIGKE)                                 
         CALL STRSYM(NC,IB,R,SIGXC)                                 
         CALL STRSYM(NC,IB,R,SIGHA)                                 
         CALL STRSYM(NC,IB,R,SIGEI)                                 
         CALL STRSYM(NC,IB,R,SIGNL)                                 
       ENDIF                                                          
C=====================================================================
C SUM THE TOTAL STRESS-TENSOR AND PRINT THEM OUT
C=====================================================================
        DO 2750 M = 1 , 6
          SIGTO(M) = SIGKE(M) + SIGHA(M) + SIGXC(M) + SIGEI(M) + 
     &               SIGAL(M) + SIGEW(M) + SIGNL(M)
2750    CONTINUE
        IF (N.NE.(N/NPRINT)*NPRINT) GO TO 2901 
        IF (N.NE.NITER .AND. IBOX.EQ.0) GO TO 2901
        WRITE (*,2900)
        WRITE (*,2700) 'STRESS', 'xx','yy','zz','yz','zx','xy'
        WRITE (*,2900)
        WRITE (*,2800) 'SIGKE ',SIGKE
        WRITE (*,2800) 'SIGHA ',SIGHA
        WRITE (*,2800) 'SIGXC ',SIGXC
        WRITE (*,2800) 'SIGAL ',SIGAL
        WRITE (*,2800) 'SIGEI ',SIGEI
        WRITE (*,2800) 'SIGEW ',SIGEW
        IF (IVPTYP.NE.0) WRITE (*,2800) 'SIGNL ',SIGNL
        WRITE (*,2900)
        WRITE (*,2800) 'SIGTO ',SIGTO
        WRITE (*,2900)
2901    CONTINUE
2700  FORMAT(4X,A6,4X,6(A2,9X))
2800  FORMAT(1X,A7,6F11.6)
2900  FORMAT(1X,73(1H-))
 115  FORMAT(/,'  Atom      Ewald     Elec-Ion   (Nonlocal)    Total')
C=====================================================================
C SYMMETRIZE THE FORCE,I.E. EIFOR AND FNLEIF                          
C=====================================================================
        IF (ISYMM.EQ.1) THEN                                            
          IF (IPRINT.GE.3) THEN
            WRITE(*,*)'EWIFOR:'                                           
            DO 8820 NSP=1,NSPEC                                           
              DO 8820 NI=1,NIONSP(NSP)                                    
                WRITE(*,8815) NSP, NI, (EWIFOR(M,NI,NSP), M=1,3)          
8820        CONTINUE                                                    
          END IF
          CALL FORSYM (NIONST,NC,IB,F0,R,NDIM9,FORCEWK1,FORCEWK2,
     $         NIONS,NSPEC,NIONSP,EWIFOR)                  
          IF (IPRINT.GE.2) THEN
            WRITE(*,*)'EIFOR:'                                            
            DO 8810 NSP=1,NSPEC                                           
              DO 8810 NI=1,NIONSP(NSP)                                    
                WRITE(*,8815) NSP, NI, (EIFOR(M,NI,NSP), M=1,3)           
8810        CONTINUE                                                    
          END IF
          CALL FORSYM (NIONST,NC,IB,F0,R,NDIM9,FORCEWK1,FORCEWK2,         
     *                   NIONS,NSPEC,NIONSP,EIFOR)                    
          IF (IPRINT.GE.2) THEN
            WRITE(*,*)'EIFOR AFTER SYMMETRIZATION:'                       
            DO 8819 NSP=1,NSPEC                                           
              DO 8819 NI=1,NIONSP(NSP)                                    
                WRITE(*,8815) NSP, NI, (EIFOR(M,NI,NSP), M=1,3)  
8819        CONTINUE                                           
            WRITE(*,*)'FNLEIF:'                                  
            DO 8817 NSP=1,NSPEC                                  
              DO 8817 NI=1,NIONSP(NSP)                           
                WRITE(*,8815) NSP, NI, (FNLEIF(M,NI,NSP), M=1,3) 
8817        CONTINUE                                           
          END IF
          CALL FORSYM (NIONST,NC,IB,F0,R,NDIM9,FORCEWK1,FORCEWK2,
     *                   NIONS,NSPEC,NIONSP,FNLEIF)          
          IF (IPRINT.GE.2) THEN
            WRITE(*,*)'FNLEIF AFTER SYMMETRIZATION:'             
            DO 8818 NSP=1,NSPEC                                  
              DO 8818 NI=1,NIONSP(NSP)                           
                WRITE(*,8815) NSP, NI, (FNLEIF(M,NI,NSP), M=1,3) 
8818        CONTINUE                                           
          END IF
        END IF                                                  
C=====================================================================
C SUM THE TOTAL FORCES ON THE IONS AND PRINT OUT
C=====================================================================
        NG = 1
        IF (IPRINT.GE.2) WRITE (*,115)
        DO 111 NSP = 1 , NSPEC
          DO 112 NI = 1 , NIONSP(NSP)
            DO 113 M = 1 , 3
              EIFOR(M,NI,NSP) = EIFOR(M,NI,NSP) + FNLEIF(M,NI,NSP)
              TIFOR(M,NI,NSP) = EIFOR(M,NI,NSP) + EWIFOR(M,NI,NSP)
 113        CONTINUE
            IF (IPRINT.LT.2) GOTO 112
C
            WRITE (*,116) NG, EWIFOR(1,NI,NSP),EIFOR(1,NI,NSP),
     &                   FNLEIF (1,NI,NSP),TIFOR(1,NI,NSP), 'Fx'
            WRITE (*,116) NG, EWIFOR(2,NI,NSP),EIFOR(2,NI,NSP),
     &                   FNLEIF (2,NI,NSP),TIFOR(2,NI,NSP), 'Fy'
            WRITE (*,116) NG, EWIFOR(3,NI,NSP),EIFOR(3,NI,NSP),
     &                   FNLEIF (3,NI,NSP),TIFOR(3,NI,NSP), 'Fz'
C 
            NG = NG + 1
 112      CONTINUE
 111    CONTINUE
      IF (ICLOCK.EQ.1) CALL PCLOCK(0)
 116  FORMAT(1X,I3,4F13.6,3X,A2)
C=====================================================================
C
C              *    DYNAMICS OF UNIT CELL    *
C
C RELAX THE UNIT CELL ALONG THE DIRECTIONS OF STRESSES (IF IBOX.EQ.1)
C=====================================================================
        IF (IBOX.EQ.0 .OR. N.LE.NDELAY) GOTO 7161
C=====================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE BMDTR }}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE INTEGRATES THE EQUATIONS OF MOTION FOR THE SIZE OF THE
C UNIT CELL
C=====================================================================
        CALL BMDTR(DIRL,DIRC,RECC,VOLC,SIGTO,SITIM,SIDAMP,SIMASS,SIDISP)
        CALL BASTR(DIRC,RECC,VOLC)
C=====================================================================
C WARNING: 
C IF YOU CALL BMDTR YOU MUST CALL GENSP AGAIN TO CALCULATE THE KINETIC
C ENERGIES OF THE BASIS STATES FOR THE NEW BOX SIZE
C=====================================================================
        CALL GENBTR(NRPLWV,NGX,NGY,NGZ,NKPTS,ENMAX,NINDPW,NPLWKP,VKPT,
     &           LPCTX,LPCTY,LPCTZ,DATAKE,RECC,RECI,IPRINT,DNLG,DNLKG)
C=====================================================================
C ONE MUST ALSO CALL FVPTR AGAIN BECAUSE THE PSEUDOPOTENTIAL HAS
C NOW TO BE KNOWN ON THE CHANGED GRID IN RECIPROCAL SPACE 
C=====================================================================
        DO 2650 NSP = 1 , NSPEC
          CALL FVPTR(NGX,NGY,NGZ,NPSPTN(NSP),RECC,VOLC,PSP(1,NSP),
     &               PSGMAX(NSP),VPS(1,NSP),DVPS(1,NSP),LPCTX,LPCTY,
     &               LPCTZ,ICHARG(NSP),IVPTYN(NSP))
 2650   CONTINUE
        WRITE (*,101) ((DIRC(I,J), I=1,3), J=1,3)
7161    CONTINUE
        IF (ICLOCK.EQ.1) CALL PCLOCK(0)
C=====================================================================
C
C           *   *       DYNAMICS OF IONS       *   *
C
C IION   = 1 IF ION-POSITIONS RELAXATION IS INCLUDED
C=====================================================================
        IF (IION.EQ.0 .OR. N.LE.NDELAY) GOTO 8000
        NINDX = 1
        DO 3930 NSP = 1 , NSPEC
          DO 3931 NI = 1 , NIONSP(NSP)
            DO 3932 M = 1 , 3
              POSIC(M,NINDX) = POSION(M,NI,NSP)
              POSICL(M,NINDX) = POSIOL(M,NI,NSP)
              EWIFC(M,NINDX) = EWIFOR(M,NI,NSP)
              EIFC(M,NINDX) = EIFOR(M,NI,NSP)
 3932       CONTINUE
            NINDX = NINDX + 1
 3931     CONTINUE
 3930   CONTINUE
C=======================================================================
C
C DO A PRESET NUMBER OF CONJUGATE GRADIENTS ITERATIONS MOVING IONS FIRST
C SINCE THE CHANGE IN UNIT CELL SIZE CAN BE HANDLED IN A SINGLE STEP
C
C=======================================================================
C FIRST CALCULATE THE ENERGY IN THE LOCAL IONIC POTENTIAL
C=======================================================================
        ENION = 0.0
        ENPOT = 0.0
        DO 3929 NPW = 1 , NPLWV
          ENION = ENION + REAL ( CVION(NPW) ) * REAL ( CHDENR(NPW) )
          ENPOT = ENPOT + REAL ( CV   (NPW) ) * REAL ( CHDENR(NPW) )
          CVOLD(NPW) = CV(NPW)
 3929   CONTINUE
        ENION = ENION * RINPLW
        ENPOT = ENPOT * RINPLW
C========================================================================
C NOW START THE CONJUGATE GRADIENTS LOOP
C=======================================================================
        ENVNL1 = 0.0
        ENVNL2 = 0.0
        ENVNLT = 0.0
        PODISI = PODISP
        DO 3000 NCGI = 1 , NIOCG1
          IF (NCGI.EQ.1) GO TO 9138 
          DO 9130 NSP = 1 , NSPEC
            CALL FELITR(NGX,NGY,NGZ,NIONSP(NSP),POSION(1,1,NSP),
     &                  CHDENG,RECC,VPS(1,NSP),EIFOR(1,1,NSP),LPCTX,
     &          LPCTY,LPCTZ,LPCTFX,LPCTFY,LPCTFZ,CWORK3,CWORK4,CWORK5)
 9130     CONTINUE
C=======================================================================
C  INITIALISE NON-LOCAL FORCE ARRAY
C=======================================================================
          DO 9136 NSP = 1 , NSPEC
            DO 9135 MU = 1 , NIONS
              DO 9134 M = 1 , 3
                FNLEIF (M,MU,NSP) = 0.0
 9134         CONTINUE
 9135       CONTINUE
 9136     CONTINUE
C=======================================================================
          IF (IVPTYP.NE.0) THEN
            IF (NLPOT.EQ.1) THEN
              CALL 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)                                    
              DO 9150 NI = 1 , NIONST
                DO 9151 M = 1 , 3    
                  FNLRL(M,NI) = 0.0  
 9151           CONTINUE             
 9150         CONTINUE               
            END IF
C
            DO 9131 NKP = 1 , NKPTS
              IF (NKPTS.GT.1) THEN
                NPUNIT = 19 + NKP
                REWIND NPUNIT
                READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV) , J=1,NBANDS)
              END IF
              IF (NLPOT.EQ.0) THEN
                CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &                    DNLG(1,1,NKP),POSION,DIRC, CPHSGR,IVPTYN)
                CALL SETVG(NRPLWV,NSPEC,NPSPTS,NIONSP,NPLWKP(NKP),
     &                 DNLKG(1,0,NKP),PSPNL,PSCALE, PSGMAX, VGNL,
     &                 DVGNL,NPSPTN,IVPTYN)
C
                DO 9137 NSP = 1 , NSPEC
                  IF (IVPTYN(NSP).EQ.0) GO TO 9137
                  CALL FNLFOR(NBANDS,NRPLWV,NIONSP(NSP),NPLWKP(NKP),
     &                    NIONS,VOLC,PSCALE(0,NSP),CPHSGR(1,1,NSP),
     &                    VGNL(1,0,NSP),DNLG(1,1,NKP),DNLKG(1,0,NKP),
     &                    CPTWFP,FORWK,OCC(1,NKP),CWOR10,CWOR11)
                  DO 9133 MU = 1 , NIONSP(NSP)
                    DO 9132 M = 1 , 3
                      FNLEIF(M,MU,NSP) = FNLEIF(M,MU,NSP) +
     &                                  WTKPT(NKP) * FORWK(M,MU)
 9132               CONTINUE
 9133             CONTINUE
 9137           CONTINUE
              ELSE
                CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,POSION,
     &             NRLPPI,NIONSP,RLCORE,CPHGRD,VKPT(1,NKP),NGX,NGY,NGZ)
                CALL FORLNL(VOLC,NGX,NGY,NGZ,CPTWFP,NKP,MPLWV,           
     &             NRPLWV,NRGRPT,NIONSP,NRLPPI,NBANDS,NKPTS,NSPEC,     
     &             CWORK1,CWORK2,CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,     
     &             PRLSCA,VRLGRD,NADGRD,MXRLSH,                        
     &             FNLRL,DVRLGR,NINDPW,NPLWKP,OCC,IOCCUP,NGPTAR,WTKPT, 
     &             CESAVE)                                              
                NG = 1                              
                DO 9140 NSP = 1 , NSPEC             
                  DO 9141 MU = 1 , NIONSP(NSP)      
                    DO 9142 M = 1 , 3               
                      FNLEIF(M,MU,NSP) = FNLRL(M,NG)
 9142               CONTINUE                        
                    NG = NG + 1                     
 9141             CONTINUE                          
 9140           CONTINUE                            
              END IF
C
 9131       CONTINUE
C
          END IF
C
C=====================================================================
C SYMMETRIZE THE FORCE,I.E. EIFOR AND FNLEIF                          
C=====================================================================
          IF (ISYMM.EQ.1) THEN                                            
            IF (IPRINT.GE.3) THEN
              WRITE(*,*)'EWIFOR:'                                           
              DO 8816 NSP=1,NSPEC                                           
                DO 8816 NI=1,NIONSP(NSP)                                    
                  WRITE(*,8815) NSP, NI, (EWIFOR(M,NI,NSP), M=1,3)          
8816          CONTINUE                                                    
              CALL FORSYM (NIONST,NC,IB,F0,R,NDIM9,FORCEWK1,
     $             FORCEWK2,NIONS,NSPEC,NIONSP,EWIFOR)
              WRITE(*,*)'EIFOR:'                                            
              DO 8811 NSP=1,NSPEC                                           
                DO 8811 NI=1,NIONSP(NSP)                                    
                  WRITE(*,8815) NSP, NI, (EIFOR(M,NI,NSP), M=1,3)           
8811          CONTINUE                                                    
            END IF
            CALL FORSYM (NIONST,NC,IB,F0,R,NDIM9,FORCEWK1,FORCEWK2,         
     *                   NIONS,NSPEC,NIONSP,EIFOR)                    
            IF (IPRINT.GE.3) THEN
              WRITE(*,*)'EIFOR AFTER SYMMETRIZATION:'                       
              DO 8813 NSP=1,NSPEC                                           
                DO 8813 NI=1,NIONSP(NSP)                                    
                  WRITE(*,8815) NSP, NI, (EIFOR(M,NI,NSP), M=1,3)  
8813          CONTINUE                                           
              WRITE(*,*)'FNLEIF:'                                  
              DO 8812 NSP=1,NSPEC                                  
                DO 8812 NI=1,NIONSP(NSP)                           
                  WRITE(*,8815) NSP, NI, (FNLEIF(M,NI,NSP), M=1,3) 
8812          CONTINUE                                           
            END IF
            CALL FORSYM (NIONST,NC,IB,F0,R,NDIM9,FORCEWK1,FORCEWK2,
     *                   NIONS,NSPEC,NIONSP,FNLEIF)          
            IF (IPRINT.GE.3) THEN
              WRITE(*,*)'FNLEIF AFTER SYMMETRIZATION:'             
              DO 8814 NSP=1,NSPEC                                  
                DO 8814 NI=1,NIONSP(NSP)                           
                  WRITE(*,8815) NSP, NI, (FNLEIF(M,NI,NSP), M=1,3) 
8814          CONTINUE                                           
            END IF
          END IF                                                  
C=====================================================================
          NINDX = 1
          DO 9190 NSP = 1 , NSPEC
            DO 9191 NI = 1 , NIONSP(NSP)
              DO 9192 M = 1 , 3
                EIFC(M,NINDX)  = EIFOR(M,NI,NSP) + FNLEIF (M,NI,NSP)
                POSIC(M,NINDX) = POSION(M,NI,NSP)
 9192         CONTINUE
              NINDX = NINDX + 1
 9191       CONTINUE
 9190     CONTINUE
C=======================================================================
 9138     CONTINUE
C=======================================================================
          DO 3011 NI = 1 , NIONST
            DO 3010 M = 1 , 3
              GRION(M,NI) = RMOVE(NI) * ( EIFC(M,NI) + EWIFC(M,NI) )
              TIFORC(M,NI) = EIFC(M,NI) + EWIFC(M,NI)
 3010       CONTINUE
 3011     CONTINUE
C
          IF (ICNSTR.NE.0) THEN
            DO 3012 NI = 1 , NEWCOO
              GRNEW(NI) = 0.0
              DO 3013 NJ = 1 , NIONST
                DO 3014 K = 1 , 3
                  GRNEW(NI) = GRNEW(NI) + 
     &                        RCNSTR(K+3*(NJ-1),NI) * TIFORC(K,NJ)
 3014           CONTINUE
 3013         CONTINUE
              IF (IPRINT.GT.0) WRITE (*,*) 'NEW FORCE:',NI,GRNEW(NI)
 3012       CONTINUE
          END IF
C=======================================================================
C  CG STEP PERFORMED FOR UNCONSTRAINED SEARCH
C=======================================================================
          IF (ICNSTR.EQ.0) THEN
            IF (NCGI.EQ.1) THEN
              DO 9020 NI = 1 , NIONST
                DO 9021 M = 1 , 3
                  DIRION(M,NI) = GRION(M,NI)
 9021           CONTINUE
 9020         CONTINUE
            ELSE
              DOT1 = 0.0
              DOT2 = 0.0
              DO 9030 NI = 1 , NIONST
                DO 9031 M = 1 , 3
                  DOT1 = DOT1 + GRION(M,NI)  * GRION(M,NI)
                  DOT2 = DOT2 + GRIONP(M,NI) * GRIONP(M,NI)
 9031           CONTINUE
 9030         CONTINUE
              GAMMA = DOT1 / DOT2
              DO 9035 NI = 1 , NIONST
                DO 9036 M = 1 , 3
                  DIRION(M,NI) = GRION(M,NI) + GAMMA * DIPION(M,NI)
 9036           CONTINUE
 9035         CONTINUE
            ENDIF
C
            DIRMAX = 0.0
            FORDDI = 0.0
            DO 9037 NI = 1 , NIONST
              DO 9038 M = 1 , 3
                DIPION(M,NI) = DIRION(M,NI)
                GRIONP(M,NI) = GRION(M,NI)
                IF (ABS(DIRION(M,NI)).GT.DIRMAX) 
     &              DIRMAX = ABS (DIRION(M,NI))
 9038         CONTINUE
 9037       CONTINUE
            SCALE = PODISP / DIRMAX
C
            DO 9045 NI = 1 , NIONST
              DO 9046 M = 1 , 3
                DIRION(M,NI) = SCALE * DIRION(M,NI)
                FORDDI = FORDDI + GRION(M,NI) * DIRION(M,NI)
 9046         CONTINUE
 9045       CONTINUE
C=======================================================================
C  CG STEP PERFORMED FOR CONSTRAINED SEARCH
C=======================================================================
          ELSE
            IF (NCGI.EQ.1) THEN
              DO 6020 NI = 1 , NEWCOO
                DIRNEW(NI) = GRNEW(NI)
 6020         CONTINUE
            ELSE
              DOT1 = 0.0
              DOT2 = 0.0
              DO 6030 NI = 1 , NEWCOO
                DOT1 = DOT1 + GRNEW(NI)  * GRNEW(NI)
                DOT2 = DOT2 + GRNEWP(NI) * GRNEWP(NI)
 6030         CONTINUE
              GAMMA = DOT1 / DOT2
              DO 6035 NI = 1 , NEWCOO
                DIRNEW(NI) = GRNEW(NI) + GAMMA * DIPNEW(NI)
 6035         CONTINUE
            ENDIF
C
            DIRMAX = 0.0
            FORDDI = 0.0
            DO 6037 NI = 1 , NEWCOO
              DIPNEW(NI) = DIRNEW(NI)
              GRNEWP(NI) = GRNEW(NI)
              IF (ABS (DIRNEW(NI)).GT.DIRMAX) DIRMAX = ABS(DIRNEW(NI))
 6037       CONTINUE
            SCALE = PODISP / DIRMAX
C
            DO 6045 NI = 1 , NEWCOO
              DIRNEW(NI) = SCALE * DIRNEW(NI)
              FORDDI = FORDDI + GRNEW(NI) * DIRNEW(NI)
 6045       CONTINUE
C
C RECONSTRUCT DIRION FOR COMPATIBILITY WITH THE PREVIOUS VERSION
C
            DO 6046 NI = 1 , NIONST
              DO 6047 M = 1 , 3
                DIRION(M,NI) = 0.0
                DO 6048 NJ = 1 , NEWCOO
                  DIRION(M,NI) = DIRION(M,NI) +
     &                           DIRNEW(NJ) * RCNSTR(M+3*(NI-1),NJ)
 6048           CONTINUE
 6047         CONTINUE
 6046       CONTINUE
          END IF
C=======================================================================
C  BELOW WE TRANSFORM DISPLACEMENT (DIRION) FROM CARTESIAN COORDINATES
C  TO TRICLINIC SYSTEM DEFINED WITH DIRC-RECC VECTORS.
C  WE USE Ai x Bj = 2 x Pi x Delta(ij) RELATION, WERE Ai ARE DIRC,
C  Bj ARE RECC VECTORS. DIRION IS IN CARTESIAN AS FORCES ARE.
C=======================================================================
          DO 9050 NI = 1 , NIONST
            DO 9051 M = 1 , 3
              DUMMY = 0.0
              DO 9052 K = 1 , 3
                DUMMY = DUMMY + RECC(M,K) * DIRION(K,NI)
 9052         CONTINUE
              POSITC(M,NI) = POSIC(M,NI) + DUMMY / TWOPI
 9051       CONTINUE
 9050     CONTINUE
          NINDX = 1
          DO 9055 NSP = 1 , NSPEC
            DO 9056 NI = 1 , NIONSP(NSP)
              DO 9057 M = 1 , 3
                POSIOT(M,NI,NSP) = POSITC(M,NINDX)
 9057         CONTINUE
            NINDX = NINDX + 1
 9056       CONTINUE
 9055     CONTINUE
C=======================================================================
C COMPUTE NEW STRUCTURE FACTORS
C=======================================================================
          DO 4109 NSP = 1 , NSPEC
            CALL FSTFTR(NGX,NGY,NGZ,NIONSP(NSP),NPLWV,POSIOT(1,1,NSP),
     &                  CSTRF(1,NSP),LPCTX,LPCTY,LPCTZ)
 4109     CONTINUE
C=======================================================================
C CALCULATE THE IONIC POTENTIAL AT TRIAL CONFIGURATION
C=======================================================================
          DO 4126 NP = 1 , MPLWV
            CWORK1(NP) = (0.0,0.0)
            CWORK2(NP) = (0.0,0.0)
 4126     CONTINUE
          DO 4127 NSP = 1 , NSPEC
            DO 4128 NP = 1 , NPLWV
              CWORK1(NP) = CWORK1(NP) + VPS(NP,NSP) * CSTRF(NP,NSP)
 4128       CONTINUE
 4127     CONTINUE
          CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)
          DO 4124 NP = 1 , NPLWV
            CVIONT(NP) = CWORK1(NP)
 4124     CONTINUE
C=======================================================================
C CALCULATE NEW ELECTRON-ION ENERGY
C=======================================================================
          ENIONT = 0.0
          DO 4130 NPW = 1 , NPLWV
           ENIONT = ENIONT + REAL ( CVIONT(NPW) ) * REAL ( CHDENR(NPW) )
 4130     CONTINUE
          ENIONT = ENIONT / NPLWV
          IF (IVPTYP.EQ.0) GO TO 9203
          IF (NLPOT.EQ.1) CALL SERLNL
     &        (NIONST,MXRLNL,MXRLSH,NIONS,NRLPTS,NRGRPT,NSPEC,
     &         DIRC,RECC,RLCORE,RMAX,NRLNL,PRLSCA,POSIOT,VRLNL,
     &         NIONSP,IRLNL,NADGRD,NRLPPI,DVRLGR,VRLGRD,NGX,NGY,NGZ)
C=======================================================================
C K-POINT LOOP STARTS
C=======================================================================
          DO 9200 NKP = 1 , NKPTS
C=======================================================================
C     READ WAVEFUNCTION FROM EXTERNAL MEMORY
C=======================================================================
            IF (NKPTS.GT.1) THEN
              NPUNIT = 19 + NKP 
              REWIND NPUNIT
              READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV) , J=1,NBANDS)
            END IF
C=======================================================================
C  WE DO NOT STORE CPHSGR AND VGNL, SO WE HAVE TO RE-CALCULATE THEM EVER
C  TIME.
C=======================================================================
            IF (NLPOT.EQ.0) THEN
              CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &                  DNLG(1,1,NKP),POSIOT,DIRC, CPHSGR,IVPTYN)
              IF (NKPTS.NE.1) 
     &         CALL SETVG(NRPLWV,NSPEC,NPSPTS,NIONSP,NPLWKP(NKP),
     &                 DNLKG(1,0,NKP),PSPNL,PSCALE, PSGMAX, VGNL,
     &                 DVGNL,NPSPTN,IVPTYN)
              DO 9201 NB = 1 , NBANDS
                CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP(NKP),
     &          DNLKG(1,0,NKP),VOLC,PSCALE,CPTWFP,CWRK20,CWRK21,
     &          CWRK22,CWRK23,CPHSGR,VGNL,CELFRC,CWOR10,CWOR11,NB,
     &          IVPTYN)
                VNL(NB,NKP) = 0.0
                DO 9202 M = 1 , NPLWKP(NKP)
                  VNL(NB,NKP) = VNL(NB,NKP) +
     &                        REAL ( CONJG (CPTWFP(M,NB)) * CELFRC(M) )
 9202           CONTINUE
 9201         CONTINUE
            ELSE
              CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,POSIOT,     
     &          NRLPPI,NIONSP,RLCORE,CPHGRD,VKPT(1,NKP),NGX,NGY,NGZ)  
              DO 4123 NB=1,NBANDS                                         
                DO 5346 M = 1 , MPLWV                                     
                  CWORK1(M) = (0.0,0.0)                                   
                  CWORK2(M) = (0.0,0.0)                                   
 5346           CONTINUE                                                  
                NINDW = NRPLWV*(NB-1)                                     
                DO 223 M=1,NPLWKP(NKP)                                    
                  CWORK1(NINDPW(M,NKP))=CPTWFP(M,NB)                      
 223            CONTINUE                                                  
                CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)                        
                CALL ENRLNL(VOLC,NGX,NGY,NGZ,VNL,NB,NKP,MPLWV,NRPLWV,     
     &         NRGRPT,NIONSP,NRLPPI,NBANDS,NKPTS,NSPEC,CWORK2,CWORK1, 
     &         CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,NADGRD, 
     &         MXRLSH,CESAVE)                                         
 4123         CONTINUE
            END IF
 9200     CONTINUE
C=======================================================================
C ENVNLT - NONLOCAL ENERGY FOR NEW IONIC POSITIONS AND OLD WAVEFUNCTIONS
C=======================================================================
          ENVNLT = 0.0
          DO 9210 NK = 1 , NKPTS
            DO 9220 NB = 1 , NBANDS
              ENVNLT = ENVNLT + 2.0 * WTKPT(NK) * VNL(NB,NK) 
     &                                          * OCC(NB,NK)
 9220       CONTINUE
 9210     CONTINUE
C
 9203     CONTINUE
C=======================================================================
C
C CALCULATE NEW EWALD ENERGY
C
C=======================================================================
         CALL EWALTR(POSITC,EWIFC,EWRLEN,EWRCSS,EWRLSS,EWRCSI,DIRC,RECC,
     &       VOLC,SIGEW,TEWENT,NIONST,ICHARC,NIONCH,NICHSQ,CPHFX,CPHFY,
     &       CPHFZ,IPRINT,NIONST,NEWPTS,MAXCX,MAXCY,MAXCZ,MAXGPX,MAXGPY,
     &       MAXGPZ,NGPTS,RFORCE,CFORCE,FORCEG,FORSIG,ENERG)
C=======================================================================
C  UPDATE TOTAL POTENTIAL
C=======================================================================
          DO 7128 M = 1 , NPLWV
            CVTRUE(M) = CV(M) + CVIONT(M) - CVION(M)
 7128     CONTINUE   
C=======================================================================
C  UPDATE WITH A SCREENED POTENTIAL IN RECIPROCAL SPACE
C=======================================================================
          DO 7250 M = 1 , NPLWV
            CWORK(M) = CV(M)
 7250     CONTINUE
          CALL FFT3D(CWORK,CWORK1,NGPTAR,1)
          DO 7251 M = 1 , NPLWV
            CV(M) = CWORK(M) * RINPLW
 7251     CONTINUE
C
          DO 7252 M = 1 , NPLWV
            CWORK(M) = CVION(M)
 7252     CONTINUE
          CALL FFT3D(CWORK,CWORK1,NGPTAR,1)
          DO 7253 M = 1 , NPLWV
            CVION(M) = CWORK(M) * RINPLW
 7253     CONTINUE
C
          DO 7254 M = 1 , NPLWV
            CWORK(M) = CVIONT(M)
 7254     CONTINUE
          CALL FFT3D(CWORK,CWORK1,NGPTAR,1)
          DO 7255 M = 1 , NPLWV
            CVIONT(M) = CWORK(M) * RINPLW
 7255     CONTINUE
C
          DO 7256 M = 1 , NPLWV
            CWORK(M) = CVOLD(M)
 7256     CONTINUE
          CALL FFT3D(CWORK,CWORK1,NGPTAR,1)
          DO 7257 M = 1 , NPLWV
            CVOLD(M) = CWORK(M) * RINPLW
 7257     CONTINUE
C
          DO 7258 M = 1 , NPLWV
            CV(M) = CV(M) + ( CVIONT(M) - CVION(M) ) * SQRT ( MIN 
     &              ( 1.0 , ABS ( CVOLD(M) / ( CVION(M) + 1.0E-6 ) ) ) )
     &                
 7258     CONTINUE
C=======================================================================
C   BACK TO THE REAL SPACE
C=======================================================================
          DO 7350 M = 1 , NPLWV
            CWORK(M) = CV(M)
 7350     CONTINUE
          CALL FFT3D(CWORK,CWORK1,NGPTAR,-1)
          DO 7351 M = 1 , NPLWV
            CV(M) = CWORK(M)
 7351     CONTINUE
C
          DO 7352 M = 1 , NPLWV
            CWORK(M) = CVION(M)
 7352     CONTINUE
          CALL FFT3D(CWORK,CWORK1,NGPTAR,-1)
          DO 7353 M = 1 , NPLWV
            CVION(M) = CWORK(M)
 7353     CONTINUE
C
          DO 7354 M = 1 , NPLWV
            CWORK(M) = CVIONT(M)
 7354     CONTINUE
          CALL FFT3D(CWORK,CWORK1,NGPTAR,-1)
          DO 7355 M = 1 , NPLWV
            CVIONT(M) = CWORK(M)
 7355     CONTINUE
C
          DO 7356 M = 1 , NPLWV
            CWORK(M) = CVOLD(M)
 7356     CONTINUE
          CALL FFT3D(CWORK,CWORK1,NGPTAR,-1)
          DO 7357 M = 1 , NPLWV
            CVOLD(M) = CWORK(M)
 7357     CONTINUE
C=======================================================================
C  INITIALISE ELECTRONIC GRADIENT TO ZERO AND SET ALPHA - TRIAL DISTANCE
C  FOR ELECTRONS
C=======================================================================
          ENGRSI = 0.0
          ALPHA = 0.01
          RCOSTH = COS (ALPHA)
          RSINTH = SIN (ALPHA)
C=======================================================================
C  START THE K-POINT LOOP
C=======================================================================
          IF (IVPTYP.EQ.1 .AND. NLPOT.EQ.1)
     &          CALL SERLNL (NIONST,MXRLNL,MXRLSH,NIONS,NRLPTS,       
     &                       NRGRPT,NSPEC,DIRC,RECC,RLCORE,RMAX,NRLNL,
     &                       PRLSCA,POSIOT,VRLNL,NIONSP,IRLNL,NADGRD, 
     &                       NRLPPI,DVRLGR,VRLGRD,NGX,NGY,NGZ)        
          DO 7100 NKP = 1 , NKPTS
            NPUNIT = 19 + NKP 
            NUNIT  = 39 + NKP
            REWIND NUNIT
            IF (NKPTS.GT.1) THEN
              REWIND NPUNIT
              READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV) , J=1,NBANDS)
            END IF
C=======================================================================
C  WE DO NOT STORE CPHSGR AND VGNL, SO WE HAVE TO RE-CALCULATE THEM EVER
C  TIME.
C=======================================================================
            IF (IVPTYP.NE.0) THEN
              IF (NLPOT.EQ.0) THEN
                CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &                    DNLG(1,1,NKP),POSIOT,DIRC, CPHSGR,IVPTYN)
                IF (NKPTS.GT.1)
     &           CALL SETVG(NRPLWV,NSPEC,NPSPTS,NIONSP,NPLWKP(NKP),
     &                 DNLKG(1,0,NKP),PSPNL,PSCALE, PSGMAX, VGNL,
     &                 DVGNL,NPSPTN,IVPTYN)
              ELSE IF (NKPTS.GT.1) THEN
                CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,
     &                      POSIOT,NRLPPI,NIONSP,RLCORE,CPHGRD,
     &                      VKPT(1,NKP),NGX,NGY,NGZ)
              END IF
            END IF          
C=======================================================================
C
C            SUBROUTINE  WFTRI
C
C  SUBROUTINE WFTRI CALCULATES THE KINETIC ENERGY AND CHARGE DENSITY FOR
C  TOTAL WAVEFUNCTIONS AT A FIXED DISTANCE ALONG THE GRADIENT DIRECTION
C  FOR THE NEW IONIC COORDINATES
C
C=======================================================================
            CALL WFTRI(NGX,NGY,NGZ,NBANOC,NKPTS,NPLWV,
     &               MPLWV,NRPLWV,NINDPW(1,NKP),
     &               NPLWKP(NKP),WTKPT,CV,CPTWFP,CPTWFL,VOLC,CELEN,
     &               NGPTAR,DATAKE(1,1,NKP),CWORK1,CWORK2,CWORK,
     &               NUNIT,NBANOC, CWORK6,CWORK7,CWORK8,CWORK9,CWOR10,
     &               WORK12,CHDENR,NKP,CWOR11,
     &               NSPEC,NIONS,NIONSP,PSCALE,
     &               DNLKG,CPHSGR,VGNL,CELFRC,
     &               CWRK20,CWRK21,CWRK22,CWRK23,VNL,IVPTYP,OCC,
     &               ENGRSI,ALPHA,CVTRUE,IVPTYN,NLPOT,NRGRPT,NRLPPI,
     &               CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,
     &               NADGRD,MXRLSH,CESAVE)
C
            IF (IVPTYP.NE.0) THEN
C=======================================================================
C NONLOCAL ENERGY FOR OLD IONIC POSITIONS AND NEW WAVE FUNCTIONS IS
C IN ENVNL2, AND IT IS CALCULATED FROM VNLT(NB,NKP) ARRAY.
C=======================================================================
C NONLOCAL ENERGY FOR NEW IONIC POSITIONS AND NEW WAVE FUNCTIONS IS
C IN ENVNL1, AND IT IS CALCULATED FROM VNL(NB,NKP) ARRAY. VNL ARE NOT
C UPDATED AS THEY WERE RECALCULATED IN WFTRI.
C=======================================================================
              IF (NLPOT.EQ.0)
     &          CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &                    DNLG(1,1,NKP),POSION,DIRC, CPHSGR,IVPTYN)
              REWIND NUNIT
              DO 9205 NB = 1 , NBANDS
                READ (NUNIT) (CPTWFL(I),I=1,NRPLWV)
                DO 9206 M = 1 , NRPLWV
                  CPTWFL(M) = CPTWFL(M) * RSINTH + CPTWFP(M,NB) * RCOSTH
 9206           CONTINUE
                IF (NLPOT.EQ.0) THEN
                  CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP(NKP),
     &                      DNLKG(1,0,NKP),VOLC,PSCALE,CPTWFL,
     &                      CWRK20,CWRK21,CWRK22,CWRK23,CPHSGR,VGNL,
     &                      CELFRC,CWOR10,CWOR11,1,IVPTYN)
                  VNLT(NB,NKP) = 0.0
                  DO 9106 M = 1 , NPLWKP(NKP)
                    VNLT(NB,NKP) = VNLT(NB,NKP) +
     &                           REAL ( CONJG (CPTWFL(M)) * CELFRC(M) )
 9106             CONTINUE
                END IF
 9205         CONTINUE
            END IF
 7100     CONTINUE
C
          IF (IVPTYP.EQ.1 .AND. NLPOT.EQ.1) THEN
            CALL 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)        
            DO 7101 NKP = 1 , NKPTS
              NUNIT = NKP + 39
              REWIND NUNIT
              IF (NKPTS.GT.1) THEN                                     
                NPUNIT = 19 + NKP
                REWIND NPUNIT                                          
                READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV) , J=1,NBANDS)
              END IF                                                   
              CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,POSION,
     &            NRLPPI,NIONSP,RLCORE,CPHGRD,VKPT(1,NKP),NGX,NGY,NGZ) 
              DO 7102 NB = 1 , NBANDS
                READ (NUNIT) (CPTWFL(I),I=1,NRPLWV)
                DO 7103 M = 1 , NRPLWV
                  CPTWFL(M) = CPTWFL(M) * RSINTH + CPTWFP(M,NB) * RCOSTH
 7103           CONTINUE
                DO 7104 M = 1 , MPLWV
                  CWORK1(M) = (0.0,0.0)
                  CWORK2(M) = (0.0,0.0)
 7104           CONTINUE
                DO 224 M=1,NPLWKP(NKP)           
                  CWORK1(NINDPW(M,NKP))=CPTWFL(M)
 224            CONTINUE                         
                CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)             
                CALL ENRLNL(VOLC,NGX,NGY,NGZ,VNLT,NB,NKP,MPLWV,
     &                      NRPLWV,NRGRPT,NIONSP,NRLPPI,NBANDS,
     &                      NKPTS,NSPEC,CWORK2,CWORK1,CPHGRD,  
     &                      NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,  
     &                      VRLGRD,NADGRD,MXRLSH,CESAVE)       
 7102         CONTINUE
 7101       CONTINUE
          END IF
C
          IF (IVPTYP.EQ.1) THEN
            ENVNL1 = 0.0
            ENVNL2 = 0.0
            DO 9103 NK = 1 , NKPTS
              DO 9104 NB = 1 , NBANDS
                ENVNL1 = ENVNL1 + 2.0 * WTKPT(NK) * VNL(NB,NK) 
     &                                            * OCC(NB,NK)
                ENVNL2 = ENVNL2 + 2.0 * WTKPT(NK) * VNLT(NB,NK) 
     &                                            * OCC(NB,NK)
 9104         CONTINUE
 9103       CONTINUE
          END IF
C=======================================================================
C  CALCULATE THE TOTAL KINETIC ENERGY FOR THE TRIAL WAVEFUNCTIONS
C=======================================================================
          ENKET = 0.0
          DO 7110 NK = 1 , NKPTS
            DO 7111 NB = 1 , NBANDS
              ENKET = ENKET + 2.0 * REAL ( CELEN(NB,NK) ) * WTKPT(NK)
     &                            * OCC(NB,NK)
 7111       CONTINUE
 7110     CONTINUE
C=======================================================================
C  TRANSFORM CHARGE DENSITY TO RECIPROCAL SPACE AND CALCULATE HARTREE
C  AND EXCHANGE-CORRELATION ENERGY AND TIDY UP POTENTIALS
C=======================================================================
C
          DO 7120 M = 1 , NPLWV
            CWORK1(M) = CHDENR(M) * RINPLW
 7120     CONTINUE
          CALL FFT3D(CWORK1,CWORK,NGPTAR,-1)
C
          IF (ISYMM.EQ.1) THEN
C
C======================================================================
C     >>>>>>>>>>>>>>>>>>> SYMMETRIZE CHARGE DENSITY  <<<<<<<<<<<<<<<<<<
C======================================================================
C
            CALL ROSYM4(A1,A2,A3,B1,B2,B3,ISY,NC,IB,V,R,RB,NPLWV,
     +                  NGX,NGY,NGZ,CWORK1,CWORK,LPCTX,LPCTY,
     +                  LPCTZ,LPCTXI,LPCTYI,LPCTZI,IPRINT)
C======================================================================
C    TRANSFORM THE CHARGE DENSITY TO REAL SPACE                        
C======================================================================
            DO 4910 M = 1 , NPLWV
              CWORK(M) = CWORK1(M)
 4910       CONTINUE
            CALL FFT3D(CWORK,CWORK2,NGPTAR,1)                    
            DO 4911 M = 1 , NPLWV
              CHDENR(M) = CWORK(M)
 4911       CONTINUE
          END IF  
C======================================================================
C  END OF SYMMETRIZATION ROUTINE
C  CHARGE DENSITY IS IN CHDENR (REAL) AND CWORK1(RECIPROCAL)
C======================================================================
          DENCT = 0.0
          DO 7122 M = 1 , NPLWV
            CVD(M) = CWORK1(M) * DIRDAT(M)
            DENCT = DENCT + REAL ( CVD(M) * CONJG ( CWORK1(M) ) )
 7122     CONTINUE
          DENCT = - DENCT / 2.0
C
          DO 7130 M = 1 , NPLWV
            CWORK1(M) = CVD(M)
 7130     CONTINUE
          CALL FFT3D(CWORK1,CWORK,NGPTAR,1)
          DO 7135 M = 1 , NPLWV
            CVD(M) = CWORK1(M)
 7135     CONTINUE
          DO 7140 M = 1 , NPLWV
            CV(M) = CHDENR(M)
 7140     CONTINUE
          CALL FEXCTR(NPLWV,CV,VOLC,SIGXC,XCENCT,XCENER,EXCDAT,XCFDAT,
     &              XCPDAT)
          DO 7150 M = 1 , NPLWV
            CV(M) = CV(M) + CVION(M)
            CV(M) = CV(M) + CVD(M)
 7150     CONTINUE
C=======================================================================
C  CALCULATE NEW TOTAL POTENTIAL ENERGY AND IONIC ENERGIES FOR THE NEW
C  CHARGE DENSITY
C=======================================================================
          ENPOTT = 0.0
          ENITEL = 0.0
          ENIELT = 0.0
          DO 7155 M = 1 , NPLWV
            ENPOTT = ENPOTT + REAL ( CV    (M) ) * REAL ( CHDENR(M) )
            ENITEL = ENITEL + REAL ( CVIONT(M) ) * REAL ( CHDENR(M) )
            ENIELT = ENIELT + REAL ( CVION (M) ) * REAL ( CHDENR(M) )
 7155     CONTINUE
          ENPOTT = ENPOTT * RINPLW
          ENITEL = ENITEL * RINPLW
          ENIELT = ENIELT * RINPLW
C
          EN1 = ENPOT + DENC + XCENC + TEWEN + ENKE + ENVNL
          EN2 = EN1 + ENIONT - ENION + TEWENT - TEWEN + ENVNLT - ENVNL
          EN3 = ENPOTT + DENCT + XCENCT + TEWEN + ENKET + ENVNL2
          EN4 = EN3 + TEWENT - TEWEN + ENITEL - ENIELT + ENVNL1 - ENVNL2
          WRITE (*,*) 'EN 1-4:',EN1,EN2,EN3,EN4
C=======================================================================
C USE RESULTS TO FIND MINIMUM IN THE TWO-DIMENSIONAL ENERGY SPACE
C=======================================================================
          E = EN4 - EN3 - EN2 + EN1
          A = 2.0 * ALPHA * ENGRSI - E
          B = EN3 - EN1 - A
          F = - FORDDI
          D = EN2 - ( EN1 + F )
C=======================================================================
C CALCULATE OPTIMUM STEPLENGTH
C=======================================================================
          STEPSI = ( 2.0 * A * D - E * F ) / ( E * E - 4.0 * B * D )
          IF (ABS(D).GT.1.E-6) THEN
            STEPIO = - ( F + E * STEPSI ) / ( 2.0 * D )
          ELSE 
            IF (ABS(E).GT.1.E-6) THEN
              STEPIO = - ( 2.0 * B * STEPSI + A ) / E
            ELSE
              STEPIO = 0.0
            END IF
          END IF
          IF (IPRINT.GE.1) THEN
            IF (IVPTYP.NE.0) THEN
              WRITE (*,*) ' ENVNL,ENVNLT  ',ENVNL,ENVNLT
              WRITE (*,*) ' ENVNL1,ENVNL2 ',ENVNL1,ENVNL2
            END IF
            WRITE (*,*) ' ENPOT,ENPOTT  ',ENPOT,ENPOTT
            WRITE (*,*) ' ENION,ENIONT  ',ENION,ENIONT
            WRITE (*,*) ' ENITEL,ENIELT ',ENITEL,ENIELT
            WRITE (*,*) ' DENC,DENCT    ',DENC,DENCT
            WRITE (*,*) ' XCENC,XCENCT  ',XCENC,XCENCT
            WRITE (*,*) ' TEWEN,TEWENT  ',TEWEN,TEWENT
            WRITE (*,*) ' ENKE,ENKET    ',ENKE,ENKET
            WRITE (*,*) ' ENGRSI,FORDDI ',ENGRSI,FORDDI
            WRITE (*,*) ' A,B,F ', A,B,F
            WRITE (*,*) ' D,E   ', D,E
            WRITE (*,*) ' STEPSI,STEPIO ',STEPSI,STEPIO
          END IF
C=======================================================================
C  SET REAL-SPACE CHARGE DENSITY TO ZERO BEFORE ELECTRONIC UPDATING
C=======================================================================
          DO 7160 M = 1 , NPLWV
            CHDENR(M) = (0.0,0.0)
            CHDENG(M) = (0.0,0.0)
 7160     CONTINUE
C=======================================================================
C  START K-POINT LOOP
C=======================================================================
          DO 7170 NKP = 1 , NKPTS
            NPUNIT = 19 + NKP
            NUNIT  = 39 + NKP
            IF (NKPTS.GT.1) THEN
              REWIND NPUNIT
              READ (NPUNIT) ((CPTWFP(I,J),I=1,NRPLWV) , J=1,NBANDS)
            END IF
            REWIND NUNIT
C=======================================================================
C  UPDATE WAVEFUNCTIONS
C=======================================================================
            DO 7171 NB = 1 , NBANDS
              READ (NUNIT) (CPTWFL(I), I = 1,NRPLWV)
              THETA = ALPHA * STEPSI
              RCOSTH = COS (THETA)
              RSINTH = SIN (THETA)
              DO 7172 M = 1 , NRPLWV
                CPTWFP(M,NB) = RCOSTH * CPTWFP(M,NB) + 
     &                         RSINTH * CPTWFL(M)
 7172         CONTINUE
 7171       CONTINUE
C=======================================================================
C  ORTHOGONALISE WAVEFUNCTIONS AND CALCULATE CHARGE DENSITY
C=======================================================================
            CALL ORSP(NBANDS,NKPTS,NRPLWV,NPLWKP(NKP),CPTWFP,
     &                                        CWORK1,CWORK2)
            CALL CHSP(NBANOC,NPLWV,NRPLWV,MPLWV,
     &                WTKPT(NKP),CPTWFP,CHDENG,CHDENR,NGPTAR,
     &                NINDPW(1,NKP),NPLWKP(NKP),CWORK1,CWORK2,
     &                OCC(1,NKP) )
C=======================================================================
C  CALCULATE NEW KINETIC ENERGIES OF WAVEFUNCTIONS
C=======================================================================
            DO 7180 NB = 1 , NBANDS
              CELEN(NB,NKP) = 0.0
              DO 7181 M = 1 , NPLWKP(NKP)
                CELEN(NB,NKP) = CELEN(NB,NKP) + DATAKE(1,M,NKP) * 
     &                          CPTWFP(M,NB) * CONJG ( CPTWFP(M,NB) )
 7181         CONTINUE
 7180       CONTINUE
C=======================================================================
C  STORE WAVEFUNCTIONS
C=======================================================================
            IF (NKPTS.GT.1) THEN
              REWIND NPUNIT
              WRITE (NPUNIT) ((CPTWFP(I,J),I=1,NRPLWV) , J=1,NBANDS)
            END IF
C=======================================================================
C  END OF K-POINT LOOP
C=======================================================================
 7170     CONTINUE
C=======================================================================
C SET NEW IONIC POSITIONS AND CALCULATE NEW STARTING ENERGIES ETC
C=======================================================================
          DO 4140 NI = 1 , NIONST
            DO 4141 M = 1 , 3
              POSIC(M,NI) = POSIC(M,NI) + STEPIO * 
     &           ( RECC(M,1) * DIRION(1,NI) + RECC(M,2) * DIRION(2,NI)
     &           + RECC(M,3) * DIRION(3,NI) ) / TWOPI
 4141       CONTINUE
 4140     CONTINUE
          PODISP = STEPIO * PODISP
          WRITE (*,*) ' PODISP,PODISI ',PODISP,PODISI
          DO 1650 NI = 1 , NIONST
            DO 1651 M = 1 , 3
              IF (POSIC(M,NI).LT.0.0) POSIC(M,NI) = POSIC(M,NI) + 1.0 -
     &                                INT ( POSIC(M,NI) )
              IF (POSIC(M,NI).GT.1.0) POSIC(M,NI) = POSIC(M,NI) -
     &                                INT ( POSIC(M,NI) )
 1651       CONTINUE
 1650     CONTINUE
          NINDX = 1
          DO 4150 NSP = 1 , NSPEC
            DO 4151 NI = 1 , NIONSP(NSP)
              DO 4152 M = 1 , 3
                POSION(M,NI,NSP) = POSIC(M,NINDX)
 4152         CONTINUE
              NINDX = NINDX + 1
 4151       CONTINUE
 4150     CONTINUE
C=======================================================================
C COMPUTE STRUCTURE FACTOR AT FINAL ION POSITIONS
C=======================================================================
          DO 5110 NSP = 1 , NSPEC
            CALL FSTFTR(NGX,NGY,NGZ,NIONSP(NSP),NPLWV,POSION(1,1,NSP),
     &                  CSTRF(1,NSP),LPCTX,LPCTY,LPCTZ)
 5110     CONTINUE
C=======================================================================
C CALCULATE THE IONIC POTENTIAL AT FINAL CONFIGURATION
C=======================================================================
          DO 5126 NP = 1 , MPLWV
            CWORK1(NP) = (0.0,0.0)
            CWORK2(NP) = (0.0,0.0)
 5126     CONTINUE
          DO 5127 NSP = 1 , NSPEC
            DO 5128 NP = 1 , NPLWV
              CWORK1(NP) = CWORK1(NP) + VPS(NP,NSP) * CSTRF(NP,NSP)
 5128       CONTINUE
 5127     CONTINUE
          CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)
          DO 5124 NP = 1 , NPLWV
            CVION(NP) = CWORK1(NP)
 5124     CONTINUE
C=======================================================================
C  TRANSFORM CHARGE DENSITY TO RECIPROCAL SPACE AND CALCULATE HARTREE
C  AND EXCHANGE-CORRELATION ENERGY AND TIDY UP POTENTIALS
C=======================================================================
          DO 5129 NP = 1 , MPLWV
            CWORK1(NP) = (0.0,0.0)
            CWORK(NP)  = (0.0,0.0)
 5129     CONTINUE
          DO 7220 M = 1 , NPLWV
            CWORK1(M) = CHDENR(M) * RINPLW
 7220     CONTINUE
          CALL FFT3D(CWORK1,CWORK,NGPTAR,-1)
C
          IF (ISYMM.EQ.1) THEN
C======================================================================
C     >>>>>>>>>>>>>>>>>>> SYMMETRIZE CHARGE DENSITY  <<<<<<<<<<<<<<<<<<
C======================================================================
C
            CALL ROSYM4(A1,A2,A3,B1,B2,B3,ISY,NC,IB,V,R,RB,NPLWV,
     +                  NGX,NGY,NGZ,CWORK1,CWORK,LPCTX,LPCTY,
     +                  LPCTZ,LPCTXI,LPCTYI,LPCTZI,IPRINT)
C======================================================================
C    TRANSFORM THE CHARGE DENSITY TO REAL SPACE                        
C======================================================================
            DO 4920 M = 1 , NPLWV
              CWORK(M) = CWORK1(M)
 4920       CONTINUE
            CALL FFT3D(CWORK,CWORK2,NGPTAR,1)                    
            DO 4921 M = 1 , NPLWV
              CHDENR(M) = CWORK(M)
 4921       CONTINUE
C======================================================================
C  END OF SYMMETRIZATION ROUTINE
C  CHARGE DENSITY IS IN CHDENR (REAL) AND CWORK1(RECIPROCAL)
C======================================================================
          END IF  
C
          DENC = 0.0
          DO 7222 M = 1 , NPLWV
            CHDENG(M) = CWORK1(M)
            CVD(M) = CWORK1(M) * DIRDAT(M)
            DENC = DENC + REAL ( CVD(M) * CONJG ( CWORK1(M) ) )
 7222     CONTINUE
          DENC = - DENC / 2.0
C
          DO 7231 M = 1 , NPLWV
            CWORK1(M) = CVD(M)
 7231     CONTINUE
          CALL FFT3D(CWORK1,CWORK,NGPTAR,1)
          DO 7236 M = 1 , NPLWV
            CVD(M) = CWORK1(M)
 7236     CONTINUE
          DO 7241 M = 1 , NPLWV
            CV(M) = CHDENR(M)
 7241     CONTINUE
          CALL FEXCTR(NPLWV,CV,VOLC,SIGXC,XCENC,XCENER,EXCDAT,XCFDAT,
     &              XCPDAT)
          DO 7243 M = 1 , NPLWV
            CV(M) = CV(M) + CVD(M) + CVION(M)
 7243     CONTINUE
C=======================================================================
C  CALCULATE EIGENVALUES IN CASE OF SUBSPACE ROTATION SWITCHED ON
C=======================================================================
          IF (IVPTYP.EQ.1 .AND. NLPOT.EQ.1)  CALL 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)                                    
          IF (ISBROT.EQ.1 .OR. IOCCUP.EQ.1) THEN
            DO 7190 NKP = 1 , NKPTS
              IF (NKPTS.GT.1) THEN
                NPUNIT = 19 + NKP
                REWIND NPUNIT
                READ (NPUNIT) ((CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
              END IF
              IF (IVPTYP.NE.0) THEN
                IF (NLPOT.EQ.0) THEN
                  CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &                      DNLG(1,1,NKP),POSION,DIRC, CPHSGR,IVPTYN)
                  IF (NKPTS.NE.1) 
     &              CALL SETVG(NRPLWV,NSPEC,NPSPTS,NIONSP,
     &                  NPLWKP(NKP),DNLKG(1,0,NKP),PSPNL,PSCALE, PSGMAX, 
     &                  VGNL,DVGNL,NPSPTN,IVPTYN)
                ELSE
                 CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,POSION,   
     &           NRLPPI,NIONSP,RLCORE,CPHGRD,VKPT(1,NKP),NGX,NGY,NGZ)
                END IF
              END IF
              CALL SUBROT(NBANDS,NKPTS,NPLWV,MPLWV,NRPLWV,NINDPW(1,NKP),
     &                  NPLWKP(NKP),CV,CPTWFP,CPTWFL,VOLC,CELEN,
     &                  VNL,NGPTAR,DATAKE(1,1,NKP),CWORK1,CWORK2,
     &                  CWORK,CWORK6,NKP,
     &                  HR, HI, AUX, FV1, FV2, FV3, CH0, NSPEC, NIONS, 
     &                  NIONSP,PSCALE,
     &                  DNLKG,CPHSGR,VGNL,CELFRC,CWRK20,
     &                  CWRK21, CWRK22, CWRK23,IVPTYP,IVPTYN,IPRINT,
     &                  NLPOT,NGX,NGY,NGZ,NRGRPT,NRLPPI,CPHGRD,NRLNL,
     &                  NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,NADGRD,
     &                  MXRLSH,CESAVE)
              DO 4502 NB = 1 , NBANDS
                EIGVAL(NB,NKP) = AUX(NB)
 4502         CONTINUE
C  STORE WAVEFUNCTIONS
C=======================================================================
              IF (NKPTS.GT.1) THEN
                REWIND NPUNIT
                WRITE (NPUNIT) ((CPTWFP(I,J),I=1,NRPLWV) , J=1,NBANDS)
              END IF
C=======================================================================
C  END OF K-POINT LOOP
C=======================================================================
 7190       CONTINUE
C======================================================================
C  RECALCULATE OCCUPATION NUMBERS AND CHARGE DENSITY
C======================================================================
            IF (NBANDS*2.GT.NELECT .AND. IOCCUP.EQ.1) THEN
              CALL EFERMI(NELECT,NBANDS,DEL,NKPTS,
     &                    NBANDS,NKPTS,WTKPT,OCC,EF,EIGVAL,SORT) 
C======================================================================
C     TAKE OUT THE SPIN-DEGENERACY, FOR IT IS CONSIDERED EXPLICITLY
C     IN THIS PROGRAM. BUT NOT IN EFERMI.
C======================================================================
              DO 4501 JKP = 1 , NKPTS
                DO 4501 JB = 1 , NBANDS
                  OCC(JB,JKP) = OCC(JB,JKP) / 2.0
 4501         CONTINUE
            END IF
C======================================================================
C     CALCULATE THE NEW CHARGE DENSITY ACCORDING TO THE NEW SET OF
C     OCCUPATION NUMBER 
C     (NO NEED TO UPDATE CELEN AND VNL, BECAUSE THEY ARE NOT RELATED TO 
C      CHARGE DENSITY, THEY ARE RELATED TO WAVE FUNCTIONS ONLY.)
C======================================================================
            DO 4601 M = 1 , NPLWV
              CHDENR(M) = (0.0,0.0)
              CHDENG(M) = (0.0,0.0)
 4601       CONTINUE
            DO 4603 NKP = 1 , NKPTS
              IF (NKPTS.GT.1) THEN 
                NPUNIT = 19 + NKP
                REWIND NPUNIT
                READ (NPUNIT) ( (CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
              END IF
              CALL CHSP(NBANOC,NPLWV,NRPLWV,MPLWV,
     &                  WTKPT(NKP),CPTWFP,CHDENG,CHDENR,NGPTAR,
     &                  NINDPW(1,NKP),NPLWKP(NKP),CWORK1,CWORK2,
     &                  OCC(1,NKP) )
 4603       CONTINUE
C=========================================================================
C NOW: CHDENR(r) CONTAINS THE NEW CHARGE DENSITY IN   R E A L  SPACE
C    -----  Change over for IOCCUP=1
C=========================================================================
          END IF
C=======================================================================
C  CALCULATE THE TOTAL KINETIC ENERGY
C=======================================================================
          ENKE = 0.0
          DO 7210 NK = 1 , NKPTS
            DO 7211 NB = 1 , NBANDS
              ENKE = ENKE + 2.0 * REAL ( CELEN(NB,NK) ) * WTKPT(NK)
     &                          * OCC(NB,NK)
 7211       CONTINUE
 7210     CONTINUE
C          IF ( (ISBROT.EQ.1 .OR. IOCCUP.EQ.1) .AND. 
C     &         NCGI.EQ.NIOCG1 ) THEN
          IF (ISBROT.EQ.1 .OR. IOCCUP.EQ.1) THEN
C=======================================================================
C  TRANSFORM CHARGE DENSITY TO RECIPROCAL SPACE AND CALCULATE HARTREE
C  AND EXCHANGE-CORRELATION ENERGY AND TIDY UP POTENTIALS
C=======================================================================
          DO 7224 M = 1 , MPLWV
            CWORK1(M) = (0.0,0.0)
            CWORK2(M) = (0.0,0.0)
 7224     CONTINUE
          DO 7221 M = 1 , NPLWV
            CWORK1(M) = CHDENR(M) * RINPLW
 7221     CONTINUE
          CALL FFT3D(CWORK1,CWORK,NGPTAR,-1)
          DENC = 0.0
          DO 7223 M = 1 , NPLWV
            CHDENG(M) = CWORK1(M)
            CVD(M) = CWORK1(M) * DIRDAT(M)
            DENC = DENC + REAL ( CVD(M) * CONJG ( CWORK1(M) ) )
 7223     CONTINUE
          DENC = - DENC / 2.0
C
          DO 7230 M = 1 , NPLWV
            CWORK1(M) = CVD(M)
 7230     CONTINUE
          CALL FFT3D(CWORK1,CWORK,NGPTAR,1)
          DO 7235 M = 1 , NPLWV
            CVD(M) = CWORK1(M)
 7235     CONTINUE
          DO 7240 M = 1 , NPLWV
            CV(M) = CHDENR(M)
 7240     CONTINUE
          CALL FEXCTR(NPLWV,CV,VOLC,SIGXC,XCENC,XCENER,EXCDAT,XCFDAT,
     &              XCPDAT)
          DO 7242 M = 1 , NPLWV
            CV(M) = CV(M) + CVD(M) + CVION(M)
 7242     CONTINUE
          END IF
C=======================================================================
C CALCULATE NEW POTENTIAL ENERGY
C=======================================================================
          ENPOT = 0.0
          ENION = 0.0
          DO 5130 NPW = 1 , NPLWV
            ENION = ENION + REAL ( CVION(NPW) ) * REAL ( CHDENR(NPW) )
            ENPOT = ENPOT + REAL ( CV   (NPW) ) * REAL ( CHDENR(NPW) )
 5130     CONTINUE
          ENION = ENION * RINPLW
          ENPOT = ENPOT * RINPLW
C=======================================================================
C  NO NEED TO UPDATE VNL IF SUBROT WAS CALLED BEFORE
C=======================================================================
          IF (IVPTYP.EQ.0 .OR. ISBROT.EQ.1 .OR. IOCCUP.EQ.1) 
     &        GO TO 9303
C=======================================================================
C K-POINT LOOP STARTS
C=======================================================================
          DO 9300 NKP = 1 , NKPTS
C=======================================================================
C     READ WAVEFUNCTION FROM EXTERNAL MEMORY
C=======================================================================
            IF (NKPTS.GT.1) THEN
              NPUNIT = 19 + NKP
              REWIND NPUNIT
              READ (NPUNIT) ((CPTWFP(I,J), I=1,NRPLWV), J=1,NBANDS)
            END IF
C=======================================================================
C  WE DO NOT STORE CPHSGR AND VGNL,SO WE HAVE TO RE-CALCULATE THEM EVERY
C  TIME.
C=======================================================================
            IF (NLPOT.EQ.0) THEN
              CALL PHASGR(NRPLWV,NSPEC,NIONSP,NIONS,NPLWKP(NKP),
     &                  DNLG(1,1,NKP),POSION,DIRC, CPHSGR,IVPTYN)
              IF (NKPTS.NE.1) 
     &          CALL SETVG(NRPLWV,NSPEC,NPSPTS,NIONSP,NPLWKP(NKP),
     &                 DNLKG(1,0,NKP),PSPNL,PSCALE, PSGMAX, VGNL,
     &                 DVGNL,NPSPTN,IVPTYN)
            ELSE IF (NKPTS.NE.1) THEN
              CALL SERLPH(NIONST,NIONS,NRGRPT,NSPEC,DIRC,RECC,
     &                    POSION,NRLPPI,NIONSP,RLCORE,CPHGRD,
     &                    VKPT(1,NKP),NGX,NGY,NGZ)
            END IF
            DO 9301 NB = 1 , NBANDS
              IF (NLPOT.EQ.0) THEN
                CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP(NKP),
     &           DNLKG(1,0,NKP),VOLC,PSCALE,CPTWFP,CWRK20,CWRK21,
     &           CWRK22,CWRK23,CPHSGR,VGNL,CELFRC,CWOR10,CWOR11,NB,
     &           IVPTYN)
                VNL(NB,NKP) = 0.0
                DO 9302 M = 1 , NPLWKP(NKP)
                  VNL(NB,NKP) = VNL(NB,NKP) +
     &                        REAL ( CONJG (CPTWFP(M,NB)) * CELFRC(M) )
 9302           CONTINUE
              ELSE
                DO 5349 M = 1 , MPLWV                          
                  CWORK1(M) = (0.0,0.0)                        
                  CWORK2(M) = (0.0,0.0)                        
 5349           CONTINUE                                       
                DO 225 M=1,NPLWKP(NKP)                         
                  CWORK1(NINDPW(M,NKP))=CPTWFP(M,NB)              
 225            CONTINUE                                       
                CALL FFT3D(CWORK1,CWORK2,NGPTAR,1)             
                CALL ENRLNL(VOLC,NGX,NGY,NGZ,VNL,NB,NKP,MPLWV,
     &                      NRPLWV,NRGRPT,NIONSP,NRLPPI,NBANDS,
     &                      NKPTS,NSPEC,CWORK2,CWORK1,CPHGRD,  
     &                      NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,  
     &                      VRLGRD,NADGRD,MXRLSH,CESAVE)       
              END IF      
 9301       CONTINUE
 9300     CONTINUE
C
 9303     CONTINUE
C
          ENVNL = 0.0
          DO 9310 NK = 1 , NKPTS
            DO 9320 NB = 1 , NBANDS
              ENVNL = ENVNL + 2.0 * WTKPT(NK) * VNL(NB,NK) 
     &                                        * OCC(NB,NK)
 9320       CONTINUE
 9310     CONTINUE
C
C=======================================================================
C CALCULATE NEW EWALD ENERGY
C=======================================================================
          CALL EWALTR(POSIC,EWIFC,EWRLEN,EWRCSS,EWRLSS,EWRCSI,DIRC,RECC,
     &        VOLC,SIGEW,TEWEN,NIONST,ICHARC,NIONCH,NICHSQ,CPHFX,CPHFY,
     &        CPHFZ,IPRINT,NIONST,NEWPTS,MAXCX,MAXCY,MAXCZ,MAXGPX,
     &        MAXGPY,MAXGPZ,NGPTS,RFORCE,CFORCE,FORCEG,FORSIG,ENERG)
C=======================================================================
C  CALCULATE TOTAL ENERGY
C=======================================================================
          TOTEN = ENPOT + DENC + XCENC + TEWEN + ENKE + ENVNL + PSCENC
C=======================================================================
C   END OF IONIC CONJUGATE GRADIENTS STEP
C=======================================================================
          WRITE (*,*)'At end of CG loop'
          WRITE (*,2600) 'TOTAL KINETIC ENERGY        ',  ENKE
          WRITE (*,2600) 'LOCAL POTENTIAL ENERGY      ',  ENPOT
          IF (IVPTYP.NE.0)
     &    WRITE (*,2600) 'NONLOCAL POTENTIAL  ENERGY  ',  ENVNL
          WRITE (*,2600) 'HARTREE ENERGY CORRECTION   ',  DENC
          WRITE (*,2600) 'EX-CORR ENERGY CORRECTION   ',  XCENC
          WRITE (*,2600) 'EWALD ENERGY                ',  TEWEN
          WRITE (*,2600) 'TOTAL ENERGY IS             ',  TOTEN
C========================================================================
C  END OF CONJUGATE GRADIENTS STEP
C========================================================================
 3000   CONTINUE
        PODISP = MIN (ABS(PODISP), PODISI)
        PODISP = MAX (PODISP,1.0E-3)
        RCONV = ABS ( (TOTEN - TOTEN0) / NIONS )
        IF (RCONV .LE. 1.E-5) NIOCG1 = MIN (NIOCG0 , 2) 
        IF (RCONV .LE. 1.E-6) NIOCG1 = 1 
C=======================================================================
C  RECALCULATED OCCUPATION NUMBERS AND CHARGE DENSITY
C=======================================================================
 8000   CONTINUE
C========================================================================
C
C             MORE SELF EXPLANATORY WRITE STATEMENTS
C      
C========================================================================
C PRINT OUT NEW UNIT CELL AND NEW IONIC POSITIONS
C========================================================================
        IF (N.LT.NITER .AND. NPRINT*(N/NPRINT).NE.N) GOTO 8020
        IF (IBOX.NE.0) WRITE (*,101) ((DIRC(I,J), I=1,3), J=1,3)
C
        WRITE (*,2900) 
        WRITE (*,8005) 'a1', 'a2', 'a3', 'Fx','Fy', 'Fz'
        WRITE (*,2900) 
        DO 8010 NSP = 1 , NSPEC
          DO 8010 NI = 1 , NIONSP(NSP)
            WRITE (*,8015) NSP, NI, (POSION(M,NI,NSP), M=1,3),
     &                       (TIFOR(M,NI,NSP), M=1,3)
 8010   CONTINUE
        WRITE (*,2900) 
 8020   CONTINUE
        IF (NEXTWR*(N/NEXTWR).EQ.N .OR. N.EQ.NITER) THEN
          REWIND 13
          WRITE (13) ((DIRC(I,J), I=1,3), J=1,3)
          WRITE (13) VOLC
          WRITE (13) ((RECC(I,J), I=1,3), J=1,3)
          WRITE (13) ((DIRI(I,J), I=1,3), J=1,3)
          WRITE (13) (((POSIOL(I,J,K), I=1,3), J=1,NIONS), K=1,NSPEC)
          WRITE (13) (((POSION(I,J,K), I=1,3), J=1,NIONS), K=1,NSPEC)
          WRITE (13) (RMOVE(J),J=1,NIONS*NSPEC)
          WRITE (13) ((CELEN(I,J), I=1,NBANDS), J=1,NKPTS)
          WRITE (13) ((VKPT(I,J), I=1,3), J=1,NKPTS)
          WRITE (13) (WTKPT(I),I=1,NKPTS)
          WRITE (13) ENMAX
          WRITE (13) NGX,NGY,NGZ,NRPLWV
          WRITE (13) ((RECI(I,J), I=1,3), J=1,3)
          REWIND 16
          WRITE (16) (CHDENR(I),I=1,NPLWV),(CHDENG(I),I=1,NPLWV)
          IF (NKPTS.EQ.1) THEN
            NPUNIT = 19 + NKPTS
            REWIND NPUNIT
            WRITE (NPUNIT) ((CPTWFP(I,J),I=1,NRPLWV), J=1,NBANDS)
          END IF
C========================================================================
C WRITE OUT THE OCCUPANCIES FOR EACH BAND AND K-POINT
C========================================================================
          REWIND 17
          DO 1080 NKP = 1 , NKPTS
            DO 1070 NBD = 1 , NBANDS
              WRITE (17,*) NKP,NBD,OCC(NBD,NKP)
1070        CONTINUE
1080      CONTINUE
        END IF
C========================================================================
C
C THIS IS THE END OF THE LOOP FOR ONE TIMESTEP OF THE ELECTRON DYNAMICS
C
C========================================================================
 1001 CONTINUE
      IF (N.LT.NITER) GOTO 1000
      IF (ICLOCK.EQ.1) CALL PCLOCK(0)
C========================================================================
C AT THE END OF THE REQUIRED NUMBER OF TIMESTEPS WRITE THE WAVEFUNCTIONS
C THE POSITIONS OF THE IONS, THE SIZE OF THE UNIT CELL, THE KINETIC
C ENERGY OF THE IONS, THE ELECTRONIC EIGENVALUES, THE COORDINATES OF THE
C SPECIAL K POINTS AND THEIR WEIGHTS AND THE STRING OF RANDOM NUMBERS
C TO EXTERNAL FILES.
C========================================================================
 109  FORMAT(1X,'TIME IS ', F12.6,' PICOSECOND')
 8005 FORMAT(1X,'NSP',' ATOM ',4X,6(A2,9X))
 8015 FORMAT(1X,I2,I4,1X,6F11.6)
 8815 FORMAT(1X,I2,I4,1X,3F11.6)
1090  FORMAT(1X,I3,I3,F8.4)
C========================================================================
C  FILE 18 IS USED FOR BAND STRUCTURE CALCULATIONS
C========================================================================
C      IF (IBANS.NE.1) THEN
C        NINDX = 1
C        DO 9106 NSP = 1 , NSPEC
C          DO 9107 NI = 1 , NIONSP(NSP)
C            DO 9108 M = 1 , 3
C              POSIC(M,NINDX) = POSION(M,NI,NSP)
C 9108       CONTINUE
C            NINDX = NINDX + 1
C 9107     CONTINUE
C 9106   CONTINUE
C        REWIND 18
C        WRITE (18) (CHDENR(I),I=1,NPLWV)
C        WRITE (18) (CV(I),I=1,NPLWV)
C        WRITE (18) NIONST
C        WRITE (18) NGX,NGY,NGZ
C        WRITE (18) A1,A2,A3,B1,B2,B3
C        WRITE (18) ((POSIC(I,J),I=1,3),J=1,NIONST)
C      END IF
C
C
C
C           BREATHE A SIGH OF RELIEF - YOU HAVE FINISHED
C
C
C
C
      STOP
      END
