c**PROGRAMA DE RESERVA, VERSION ORIGINAL.

C     ABHTATOMOPM.  A PROGRAM TO COMPUTE VARIATIONALLY OPTIMIZED EFFECTIVE    
C     ATOMIC POTENTIALS.  J.D. TALMAN.                                    
C     REF. IN COMP. PHYS. COMMUN. 54 (1989) 85                          ABHT0000
C                                                                       ABHT0017
C  THIS PROGRAM COMPUTES A VARIATIONALLY OPTIMIZED CENTRAL              ABHT0018
C  POTENTIAL FOR ATOMS OR IONS.  TWO INPUT FILES, OPMDATA AND           ABHT0019
C  INPOT ARE REQUIRED AND TWO OUTPUT FILES, OPMOUT AND POTOUT           ABHT0020
C  ARE WRITTEN.  OPMDATA SPECIFIES THE ATOM OR ION, AND INPOT           ABHT0021
C  IS AN INPUT POTENTIAL, WHICH MAY BE TAKEN TO BE ZERO.                ABHT0022
C  OPMOUT IS THE OUTPUT LISTING, AND POTOUT CONTAINS THE R              ABHT0023
C  VALUES FOR WHICH THE POTENTIAL IS CALCULATED, THE OPTIMIZED          ABHT0024
C  POTENTIAL, THE EXCHANGE PART OF THE POTENTIAL, THE SINGLE--          ABHT0025
C  PARTICLE ENERGIES, AND THE RESULTING REDUCED RADIAL WAVE             ABHT0026
C  FUNCTIONS.                                                           ABHT0027
C                                                                       ABHT0028
C  THIS PROGRAM IS BASED ON THE THEORY DEVELOPED IN THE ARTICLE         ABHT0029
C  'AN OPTIMIZED EFFECTIVE ATOMIC CENTRAL POTENTIAL' BY J.D TALMAN      ABHT0030
C  AND W.F. SHADWICK, PHYSICAL REVIEW A14, 36(1976).                    ABHT0031
C                                                                       ABHT0032
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               ABHT0035
      DIMENSION IMAXA(20),Y1(800),Y2(800),V(800),VB(800),DNUM(800),VDIR(ABHT0037
     *   800)                                                           ABHT0038
      DIMENSION NP(20),EG(20),LA(20),IO(20),O(20),FAC(0:20)             ABHT0039
      DIMENSION AKERN(5050),QA(800)                                     ABHT0040
      DIMENSION PHI2(200),PSI(200),S2(800),PHI(800,20),HA(800),T(800),S1ABHT0041
     *   (800),RHO(800),AR(800),WR(800),SQRWR(800)                      ABHT0042
      CHARACTER TITLE*50,CH1*2,CH4*8,CH5*6,CH2*1,CH3*1,CH6*4,CH7*1,CH8*7ABHT0043
      LOGICAL ENDING                                                    ABHT0044
      DATA CH4,CH5,CH6,CH7 / '12345678','SPDFGH','PHI(',')'/            ABHT0045
C                                                                       ABHT0046
C     OPEN INPUT FILE "OPMDATA", OUTPUT FILE "OPMOUT"                   ABHT0047

c..... Dario
C.....OPEN INPUT FILE "talman.inp"
c.....OUTPUT FILE "talman.out"
C                                                                       
      OPEN (UNIT=11,FILE='talman.inp',status='old')    !datos entrada                  
      OPEN (UNIT=8,FILE='talman.out',status='unknown')                                
      OPEN (UNIT=7,FILE='potout.dat')                                      
      OPEN(UNIT=10,FILE='pot.dat')  !potencial final
      open(unit=9,file='fun.dat')  !funcion
                                                                        ABHT0052
C  ND1,ND2: DIMENSIONS OF PHI ARRAY                                     ABHT0053
C                                                                       ABHT0054
      ND1 = 800                                                         ABHT0055
      ND2 = 20                                                          ABHT0056
C                                                                       ABHT0057
C  READ PROGRAM INFORMATION INPUT                                       ABHT0058
C                                                                       ABHT0059
      READ (11,'(A50)') TITLE                                           ABHT0060
      READ (11,*,END=460) NM                                            ABHT0061
      IF (NM.EQ.0) STOP                                                 ABHT0062
      READ (11,*,END=460) (LA(I),I=1,NM)                                ABHT0063
      READ (11,*,END=460) (NP(I),I=1,NM)                                ABHT0064
      READ (11,*,END=460) Z                                             ABHT0065
      READ (11,*,END=460) (IO(I),I=1,NM)                                ABHT0066
      READ (11,*,END=460) NDIM                                          ABHT0067
      IF (NDIM.GT.ND1/4) STOP 'STOP-INSUFFICIENT ARRAY DIMENSIONS'      ABHT0068
      READ (11,*,END=460) ROMIN,H,ALPHA                                 ABHT0069
      READ (11,*,END=460) POTACC,NITER                                  ABHT0070
C                                                                       ABHT0071
C  WRITE DETAILS OF THE CALCULATION                                     ABHT0072
C                                                                       ABHT0073
      WRITE(8,*) '1  OPTIMIZED POTENTIAL MODEL CALCULATION'            
      WRITE(8,*) '     '//TITLE                                       
      WRITE(8,6000)                                                    
      DO 100 I = 1, NM                                                  ABHT0077
         NPI = NP(I)                                                    ABHT0078
         LI = LA(I)+1                                                   ABHT0079
         CH2 = CH4(NPI:NPI)                                             ABHT0080
         CH3 = CH5(LI:LI)                                               
         CH1 = CH2//CH3                                                 
         WRITE(8,6100) CH1,IO(I)                                       
         EG(I) = 0.0                                                    ABHT0084
         O(I) = IO(I)                                                   ABHT0085
         NP(I) = NP(I)-LA(I)                                            ABHT0086
  100 CONTINUE                                                          ABHT0087
      WRITE(8,6200) INT(Z)                                              ABHT0088
      WRITE(8,6300) NDIM                                                ABHT0089
      WRITE(8,6400) ROMIN                                               ABHT0090
      WRITE(8,6500) H                                                   ABHT0091
      WRITE(8,6600) ALPHA                                               ABHT0092
      WRITE(8,6700) POTACC                                              ABHT0093
      WRITE(8,6800) NITER                                               ABHT0094
c      WRITE(8,*) INT(Z)                                                  
c      WRITE(8,*) NDIM                                                    
c      WRITE(8,*) ROMIN                                                   
c      WRITE(8,*) H                                                      
c      WRITE(8,*) ALPHA                                             
c      WRITE(8,*) POTACC                                             
c      WRITE(8,*) NITER                                              

C                                                                       ABHT0095
C  CALCULATE NUMBER OF ELECTRONS (CHA) AND FACTORIAL TABLE              ABHT0096
C                                                                       ABHT0097
      CHA = 0.0                                                         ABHT0098
      DO 110 I = 1, NM                                                  ABHT0099
         CHA = CHA+O(I)                                                 ABHT0100
  110 CONTINUE                                                          ABHT0101
      FAC(0) = 1.0                                                      ABHT0102
      DO 120 I = 1, 20                                                  ABHT0103
         FAC(I) = I*FAC(I-1)                                            ABHT0104
  120 CONTINUE                                                          ABHT0105
C                                                                       ABHT0106
C  SET UP RADIAL MESH AR(I) AND ARRAYS WR(I), SQRWR(I)                  ABHT0107
C                                                                       ABHT0108
      CALL RADMSH (ROMIN,H,ALPHA,NDIM,AR,WR,SQRWR)                      ABHT0109
      WRITE(8,8300) (AR(I),I=1,NDIM)                                   
C                                                                       ABHT0111
C  READ INITIAL POTENTIAL                                               ABHT0112
C                                                                       ABHT0113
      CALL POTRD (V,AR,NDIM)                                            ABHT0114
C                                                                       ABHT0115
      ENDING = .FALSE.                                                  ABHT0116
      WGHT = 0.5                                                        ABHT0117
      IRS = 1                                                           ABHT0118

C                                                                       ABHT0119
C  ENTRY POINT TO BEGIN A NEW ITERATION                                 ABHT0120
C                                                                       ABHT0121
  130 CONTINUE                                                          ABHT0122
      IF (ENDING) GO TO 140                                             ABHT0123
      WRITE(8,6900) IRS                                                
      WRITE(8,7000) (V(I),I=1,NDIM)
  140 WRITE(8,7100)                                                    
C                                                                       ABHT0127
C  SUBROUTINE RADEQS CALCULATES THE SINGLE-PARTICLE ENERGIES            ABHT0128
C  AND WAVE FUNCTIONS AND STORES THEM IN THE ARRAYS EG AND PHI          ABHT0129
C  IMAXA(N) GIVES THE MAXIMUM R VALUE AT WHICH A PARTICULAR             ABHT0130
C  WAVE FUNCTION IS CALCULATED.                                         ABHT0131
C                                                                       ABHT0132
      CALL RADEQS (PHI,NM,NDIM,ND1,ND2,V,AR,WR,SQRWR,EG,LA,NP,H,ALPHA,Z,ABHT0133
     *   T,VB,DNUM,Y1,Y2,S1,S2,IMAXA,ENDING)                            ABHT0134
C  CALCULATE NPOT, THE NUMBER OF MESH POINTS IN THE POTENTIAL           ABHT0136
C  CALCULATION.                                                         ABHT0137
C                                                                       ABHT0138
      NPOT = 0                                                          ABHT0139
      DO 150 N = 1, NM                                                  ABHT0140
         IF (O(N).EQ.0.0) GO TO 150                                     ABHT0141
         NPOT = MAX(NPOT,IMAXA(N))                                      ABHT0142
  150 CONTINUE                                                          ABHT0143
C                                                                       ABHT0144
C  INITIALIZE KERNEL AND RHS TO ZERO                                    ABHT0145
C                                                                       ABHT0146
      VEX = 0.0                                                         ABHT0147
      IF (ENDING) GO TO 180                                             ABHT0148
      IND = 0                                                           ABHT0149
      DO 160 I = 1, NDIM                                                ABHT0150
         DO 160 J = I, NDIM                                             ABHT0151
            IND = IND+1                                                 ABHT0152
            AKERN(IND) = 0.0                                            ABHT0153
  160 CONTINUE                                                          ABHT0154
      DO 170 I = 1, NDIM                                                ABHT0155
         QA(I) = 0.0                                                    ABHT0156
  170 CONTINUE                                                          ABHT0157
  180 CONTINUE                                                          ABHT0158

C                                                                       ABHT0159
C  BEGIN THE SUM OVER THE ORBITALS TO OBTAIN THE EXCHANGE POTENTIAL     ABHT0160
C  ENERGY (VEX), THE KERNEL (AKERN) AND THE RHS (QA) OF THE INTEGRAL    ABHT0161
C  EQUATION                                                             ABHT0162
C                                                                       ABHT0163
      DO 250 N = 1, NM                                                  ABHT0164
         IF (O(N).EQ.0.0) GO TO 250                                     ABHT0165
         E = EG(N)                                                      ABHT0166
         L = LA(N)                                                      ABHT0167
C                                                                       ABHT0168
C  DETERMINE RANGE OF ORBITAL N AND STORE PHI(I,N) IN T(I)              ABHT0169
C                                                                       ABHT0170
         DO 190 I = 1, NPOT                                             ABHT0171
            T(I) = PHI(I,N)                                             ABHT0172
            IF (T(I).NE.0.0) NPHI = I                                   ABHT0173
  190    CONTINUE                                                       ABHT0174
         IF (ENDING) GO TO 210                                          ABHT0175
C                                                                       ABHT0176
C  CALL SR SECSOL TO CALCULATE THE SECOND SOLUTION PHI2 AND THE         ABHT0177
C  AUXILIARY FUNCTION PSI                                               ABHT0178
C                                                                       ABHT0179
         CALL SECSOL (PHI2,PSI,NDIM,E,Z,L,NPHI,V,ND1,AR,H,ALPHA,Y1,Y2,S1ABHT0180
     *      ,S2,T,VB,WR,SQRWR,DNUM)                                     ABHT0181
C                                                                       ABHT0182
C  SET UP KERNEL OF INTEGRAL EQUATION (SEE EQS. (11) AND (13))          ABHT0183
C                                                                       ABHT0184
         IND = 0                                                        ABHT0185
         DO 200 I = 1, NPOT                                             ABHT0186
            AA = PHI(I,N)                                               ABHT0187
            DO 200 J = I, NPOT                                          ABHT0188
               IND = IND+1                                              ABHT0189
               IF (J.GT.NPHI) GO TO 200                                 ABHT0190
               CC = PHI(J,N)*(PHI2(I)-PSI(I))-PSI(J)*PHI(I,N)           ABHT0191
               AKERN(IND) = AKERN(IND)+O(N)*AA*PHI(J,N)*CC              ABHT0192
  200    CONTINUE                                                       ABHT0193
  210    CONTINUE                                                       ABHT0194
C                                                                       ABHT0195
C  CALCULATE THE HARTREE-FOCK EXCHANGE POTENTIAL FOR THE ORBITAL N      ABHT0196
C                                                                       ABHT0197
         CALL EXPOT (HA,N,LA,O,NDIM,NPHI,NM,PHI,ND1,ND2,FAC,H,AR,WR,Y1, ABHT0198
     *      Y2,S1,S2,ENDING)                                            ABHT0199
C                                                                       ABHT0200
C  ACCUMULATE THE CONTRIBUTION TO THE EXCHANGE POTENTIAL ENERGY (VEX)   ABHT0201
C                                                                       ABHT0202
         DO 220 I = 1, NPHI                                             ABHT0203
            Y2(I) = PHI(I,N)*HA(I)                                      ABHT0204
  220    CONTINUE                                                       ABHT0205
         CALL XINT (Y2,S2,H,NPHI,AR,WR,DEFI,ENDING)                     ABHT0206
         CA = DEFI                                                      ABHT0207
         VEX = VEX-0.5*CA                                               ABHT0208
         IF (ENDING) GO TO 250                                          ABHT0209
C                                                                       ABHT0210
C  ACCUMULATE THE CONTRIBUTION TO THE RIGHT HAND SIDE OF THE INTEGRAL   ABHT0211
C  EQUATION Q(R) (SEE EQ.(12))                                          ABHT0212
C                                                                       ABHT0213
         DO 230 I = 1, NPHI                                             ABHT0214
            T(I) = PSI(I)*HA(I)                                         ABHT0215
            Y1(I) = PHI2(I)*HA(I)                                       ABHT0216
  230    CONTINUE                                                       ABHT0217
         CALL XINT (T,VB,H,NPHI,AR,WR,DEFI,ENDING)                      ABHT0218
         CB = DEFI                                                      ABHT0219
         CALL XINT (Y1,S1,H,NPHI,AR,WR,DEFI,ENDING)                     ABHT0220
         DO 240 I = 1, NPHI                                             ABHT0221
            XY = PHI(I,N)*S1(I)+(CA-S2(I))*PHI2(I)                      ABHT0222
            XY = XY-CA*PSI(I)-CB*PHI(I,N)                               ABHT0223
            QA(I) = QA(I)-PHI(I,N)*XY                                   ABHT0224
  240    CONTINUE                                                       ABHT0225
  250 CONTINUE                                                          ABHT0226
C                                                                       ABHT0227
C  TOTAL ENERGY CALCULATION (SEE EQ.(2)).  ESP IS THE SINGLE PARTICLE   ABHT0228
C  ENERGY, THE KINETIC ENERGY + THE NUCLEAR ATTRACTION ENERGY.  THIS    ABHT0229
C  IS CALCULATED USING THE RADIAL SCHROEDINGER EQUATION.  EVD IS THE    ABHT0230
C  DIRECT TWO-PARTICLE ENERGY.  RHO(R) IS THE ELECTRONIC CHARGE DENSITY ABHT0231
C  (MULTIPLIED BY 4*PI*R**2) AND VDIR(R) IS THE ELECTROSTATIC POTENTIAL ABHT0232
C  PRODUCED BY RHO.                                                     ABHT0233
C                                                                       ABHT0234
      DO 270 I = 1, NPOT                                                ABHT0235
         S = 0.0                                                        ABHT0236
         DO 260 N = 1, NM                                               ABHT0237
            S = S+O(N)*PHI(I,N)**2                                      ABHT0238
  260    CONTINUE                                                       ABHT0239
         RHO(I) = S                                                     ABHT0240
  270 CONTINUE                                                          ABHT0241
      DO 280 I = 1, NPOT                                                ABHT0242
         Y1(I) = RHO(I)                                                 ABHT0243
         Y2(I) = RHO(I)/AR(I)                                           ABHT0244
  280 CONTINUE                                                          ABHT0245
      CALL XINT (Y1,S1,H,NPOT,AR,WR,DEFI,ENDING)                        ABHT0246
      CALL XINT (Y2,S2,H,NPOT,AR,WR,DEFI,ENDING)                        ABHT0247
      CA = DEFI                                                         ABHT0248
      DO 290 I = 1, NPOT                                                ABHT0249
         VDIR(I) = 2.0*(S1(I)/AR(I)+CA-S2(I))                           ABHT0250
  290 CONTINUE                                                          ABHT0251
      ESP = 0.0                                                         ABHT0252
      DO 300 N = 1, NM                                                  ABHT0253
         ESP = ESP+O(N)*EG(N)                                           ABHT0254
  300 CONTINUE                                                          ABHT0255
      DO 310 I = 1, NPOT                                                ABHT0256
         Y1(I) = RHO(I)*V(I)                                            ABHT0257
  310 CONTINUE                                                          ABHT0258
      CALL XINT (Y1,S1,H,NPOT,AR,WR,DEFI,ENDING)                        ABHT0259
      ESP = ESP-DEFI                                                    ABHT0260
      DO 320 I = 1, NPOT                                                ABHT0261
         Y1(I) = RHO(I)*VDIR(I)                                         ABHT0262
  320 CONTINUE                                                          ABHT0263
      CALL XINT (Y1,S1,H,NPOT,AR,WR,DEFI,ENDING)                        ABHT0264
      EVD = 0.5*DEFI                                                    ABHT0265
      DO 330 I = 1, NPOT                                                ABHT0266
         Y1(I) = 2.0*Z*RHO(I)/AR(I)                                     ABHT0267
  330 CONTINUE                                                          ABHT0268
      CALL XINT (Y1,S1,H,NPOT,AR,WR,DEFI,ENDING)                        ABHT0269
      EKIN = ESP+DEFI                                                   ABHT0270
      ETOT = EVD+ESP+VEX                                                ABHT0271
      WRITE(8,7200) ETOT                                               
      WRITE(8,7300) EKIN                                               
      WRITE(8,7400) EVD                                                
      WRITE(8,7500) VEX
	IF (ENDING) STOP                                                  


c      if (ending)goto 1435
C                                                                       ABHT0277
C  SOLVE THE INTEGRAL EQUATION FOR THE NEW EXCHANGE POTENTIAL           ABHT0278
C                                                                       ABHT0279
      CALL SOLVE (AKERN,QA,NPOT)                                        ABHT0280
      DO 340 I = 1, NPOT                                                ABHT0281
         QA(I) = QA(I)/(H*WR(I))                                        ABHT0282
  340 CONTINUE                                                          ABHT0283
C                                                                       ABHT0284
C  CORRECT THE SOLUTION AT THE ENDS                                     ABHT0285
C                                                                       ABHT0286
      QA(1) = QA(2)                                                     ABHT0287
      NPOT = NPOT-1                                                     ABHT0288
      CCC = -2.0/AR(NPOT)-QA(NPOT)                                      ABHT0289
      DO 350 I = 1, NPOT                                                ABHT0290
         QA(I) = QA(I)+CCC                                              ABHT0291
  350 CONTINUE                                                          
      WRITE(8,7600) (QA(I),I=1,NPOT)                                   
C                                                                       ABHT0294
C  CALCULATE THE POTENTIAL FOR THE NEXT ITERATION.  AMXDEV IS THE       ABHT0295
C  LARGEST DIFFERENCE BETWEEN SUCCESSIVE ITERATIONS FOR THE POTENTIAL   ABHT0296
C                                                                       ABHT0297
      AMXDEV = 0.0                                                      ABHT0298
      DO 360 I = 1, NPOT                                                ABHT0299
         CCC = VDIR(I)+QA(I)                                            ABHT0300
         VSAVE = V(I)                                                   ABHT0301
         V(I) = (1.0-WGHT)*V(I)+WGHT*CCC                                ABHT0302
         Y1(I) = V(I)-VSAVE                                             ABHT0303
         AMXDEV = MAX(AMXDEV,ABS(Y1(I)))                                ABHT0304
  360 CONTINUE                                                          ABHT0305
C                                                                       ABHT0306
C  EXTEND THE POTENTIAL TO LARGER R VALUES                              ABHT0307
C                                                                       ABHT0308
      DO 370 I = NPOT, NDIM                                             ABHT0309
         V(I) = 2.0*(CHA-1.0)/AR(I)                                     ABHT0310
         QA(I) = -2.0/AR(I)                                             ABHT0311
         Y1(I) = 0.0                                                    ABHT0312
  370 CONTINUE                                                          ABHT0313
C                                                                       ABHT0314
C  ESTIMATE SINGLE PARTICLE ENERGIES FOR NEW POTENTIAL USING FIRST      ABHT0315
C  ORDER PERTURBATION THEORY                                            ABHT0316
C                                                                       ABHT0317
      DO 390 N = 1, NM                                                  ABHT0318
         DO 380 I = 1, NDIM                                             ABHT0319
            Y2(I) = Y1(I)*PHI(I,N)**2                                   ABHT0320
  380    CONTINUE                                                       ABHT0321
         CALL XINT (Y2,S2,H,NDIM,AR,WR,DEFI,ENDING)                     ABHT0322
         EG(N) = EG(N)+DEFI                                             ABHT0323
  390 CONTINUE                                                          ABHT0324
      IRS = IRS+1                                                       ABHT0325
      IF (AMXDEV.LT.POTACC) GO TO 400                                   ABHT0326
      IF (IRS.LE.NITER) GO TO 130                                       ABHT0327
C                                                                       ABHT0328
C   COMPLETION OF POTENTIAL ITERATIONS                                  ABHT0329
C                                                                       ABHT0330
  400 CONTINUE                                                          ABHT0331

      WRITE(8,7700)                                                    
C                                                                       ABHT0333
C  PRINT THE REDUCED RADIAL WAVE FUNCTIONS ON THE RHO MESH              ABHT0334
C                                                                       ABHT0335
      DO 410 N = 1, NM                                                  ABHT0336
         NPI = NP(N)+LA(N)                                              ABHT0337
         LI = LA(N)+1                                                   ABHT0338
         CH2 = CH4(NPI:NPI)                                             ABHT0339
         CH3 = CH5(LI:LI)                                               ABHT0340
         CH8 = CH6//CH2//CH3//CH7                                       ABHT0341
c         WRITE(8,7800) CH8,(PHI(I,N),I=1,NDIM)
c      WRITE(8,*) (PHI(I,N),I=1,NDIM)                         
  410 CONTINUE                                                          ABHT0343
C                                                                       ABHT0344
C  WRITE RADIAL COORDINATES, EFFECTIVE POTENTIAL, EXCHANGE POTENTIAL,   ABHT0345
C  SINGLE PARTICLE ENERGIES, RADIAL WAVE FUNCTIONS  TO FILE 'POTOUT'    ABHT0346
C                                                                       ABHT0347
      WRITE (7,7900) NDIM,NM                                            ABHT0348
      WRITE (7,8000) (AR(I),I=1,NDIM)                                   ABHT0349
      WRITE (7,8000) (V(I),I=1,NDIM)

c      do i=1,ndim
c	write(6,*)v(i)
c	end do
	WRITE (7,8000) (QA(I),I=1,NDIM)                                   
      DO 420 N = 1, NM                                                  ABHT0352
         WRITE (7,8100) LA(N),NP(N)+LA(N),EG(N)                         ABHT0353
         WRITE (7,8000) (PHI(I,N),I=1,NDIM)                             ABHT0354
  420 CONTINUE                                                          ABHT0355

      do i=1,ndim
c	write(10,8000)ar(i),V(i,1),V(i,2),V(i,3),V(i,4)
	write(9,8435)ar(i),phi(i,1),phi(i,2),phi(i,3),phi(i,4)
	enddo
C                                                                       ABHT0356
C  CARRY OUT THE ENERGY CALCULATION OF THE REFINED MESH                 ABHT0357
C                                                                       ABHT0358

	ending = .true.
	do i=1,ndim
c		y2(i) = (1.0d0/AR(i) ) *PHI(i,2)**2
		y2(i) = (AR(i) ) *PHI(i,2)**2
	enddo
          CALL XINT (Y2,S2,H,NDIM,AR,WR,DEFI,ENDING)
	print*,'norma = ',defi

 

      ENDING = .TRUE.                                                   ABHT0359
      WRITE(8,8200)                                                    

C                                                                       ABHT0361
C  INTERPOLATE V(R) ONTO THE REFINED MESH                               ABHT0362
C                                                                       ABHT0363
      DO 430 I = 1, NDIM                                                ABHT0364
         T(I) = V(I)                                                    ABHT0365
  430 CONTINUE                                                          ABHT0366
      V(2) = 0.75*T(1)+0.25*T(2)                                        ABHT0367
      V(3) = 0.5*(T(1)+T(2))                                            ABHT0368
      V(4) = 0.25*T(1)+0.75*T(2)                                        ABHT0369
      NDIMM = NDIM-3                                                    ABHT0370
      DO 440 I = 1, NDIMM                                               ABHT0371
         A = T(I)                                                       ABHT0372
         B = T(I+1)                                                     ABHT0373
         C = T(I+2)                                                     ABHT0374
         D = T(I+3)                                                     ABHT0375
         V(4*I+1) = B                                                   ABHT0376
         V(4*I+2) = (-7.0*A+105.0*B+35.0*C-5.0*D)/128.0                 ABHT0377
         V(4*I+3) = (-A+9.0*(B+C)-D)/16.0                               ABHT0378
         V(4*I+4) = (-5.0*A+35.0*B+105.0*C-7.0*D)/128.0                 ABHT0379
  440 CONTINUE                                                          ABHT0380
C                                                                       ABHT0381
C  SET UP THE REFINED MESH POINTS AND WR AND SQRWR                      ABHT0382
C                                                                       ABHT0383
      H = H/4.0                                                         ABHT0384
      NDIM = 4*NDIM                                                     ABHT0385
      CALL RADMSH (ROMIN,H,ALPHA,NDIM,AR,WR,SQRWR)                      ABHT0386
C                                                                       ABHT0387
C  EXTRAPOLATE THE POTENTIAL FOR LARGE R VALUES                         ABHT0388
C                                                                       ABHT0389
      NDIMM = NDIM-8                                                    ABHT0390
      DO 450 I = NDIMM, NDIM                                            ABHT0391
         V(I) = 2.0*(CHA-1.0)/AR(I)                                     ABHT0392
  450 CONTINUE                                                          ABHT0393
     	DO I=1,ndim
	WRITE(10,*)ar(i),V(I)
	END DO

     	DO I=1,ndim
	WRITE(105,*)ar(i),V(I)*ar(i)/2.0d0 - 11.0d0
	END DO
	
	
      GO TO 130                                                         ABHT0394

C                                                                       ABHT0395
C  END-OF-FILE REACHED ON INPUT DATA FILE                               ABHT0396
C                                                                       ABHT0397
  460 WRITE(8,8400)                                                    
   	

	STOP                                                              
C                                                                       ABHT0400
C                                                                       ABHT0401
 6000 FORMAT (/,'        ORBITAL  OCCUPATION')                          ABHT0402
 6100 FORMAT ('          ',A2,'        ',I2)                            ABHT0403
 6200 FORMAT (/,'        NUCLEAR CHARGE: ',I2)                          ABHT0404
 6300 FORMAT (/,'        NUMBER OF MESH POINTS:',I3)                    ABHT0405
 6400 FORMAT (/,'        MINIMUM RHO VALUE:',F10.5)                     ABHT0406
 6500 FORMAT (/,'        DELTA RHO VALUE:',F10.5)                       ABHT0407
 6600 FORMAT (/,'        ALPHA VALUE:',F10.5)                           ABHT0408
 6700 FORMAT (/,'        POTENTIAL ACCURACY PARAMETER:',F10.5)          ABHT0409
 6800 FORMAT (/,'        MAXIMUM NUMBER OF ITERATIONS:',I2)             ABHT0410
 6900 FORMAT (/,'  ITERATION NUMBER ',I2)                               ABHT0411
 7000 FORMAT (/,'    POTENTIAL VALUES',//(' ',10F10.5))                 ABHT0412
 7100 FORMAT (/,'    SINGLE PARTICLE ENERGIES',/)                       ABHT0413
 7200 FORMAT (/,'    TOTAL ENERGY:',F12.5)                              ABHT0414
 7300 FORMAT (/,'    KINETIC ENERGY:',F12.5)                            ABHT0415
 7400 FORMAT (/,'    DIRECT TWO-PARTICLE ENERGY:',F12.5)                ABHT0416
 7500 FORMAT (/,'    EXCHANGE ENERGY:',F12.5)                           ABHT0417
 7600 FORMAT (/,'    NEW EXCHANGE POTENTIAL',//(' ',10F10.5))           ABHT0418
 7700 FORMAT (/,'    COMPLETION OF ITERATIONS FOR EFFECTIVE POTENTIAL') ABHT0419
 7800 FORMAT (/,'    ',A7,//(' ',10F10.5))                              ABHT0420
 7900 FORMAT (2I5)                                                      ABHT0421
 8000 FORMAT (5G14.6)                                                   ABHT0422
 8100 FORMAT (2I5,G14.6)                                                ABHT0423
 8200 FORMAT (/,' ENERGY CALCULATION ON THE REFINED MESH')              ABHT0424
 8300 FORMAT (/,'    RADIAL MESH COORDINATES',//(' ',10F10.5))          ABHT0425
 8400 FORMAT (/,' STOP--END-OF-FILE REACHED ON INPUT DATA FILE')        ABHT0426
 8435 FORMAT(G14.6,G14.6,G14.6,G14.6,G14.6) 
C                                                                       ABHT0427


	END                                                               
      SUBROUTINE SOLVE (A,R,ND)                                         ABHT0429
C                                                                       ABHT0430
C  SOLVE A SYSTEM OF ND EQUATIONS IN ND UNKNOWNS WITH THE COEFFICIENTS  ABHT0431
C  GIVEN IN SYMMETRIC STORAGE MODE IN ARRAY A AND THE RIGHT HAND SIDE   ABHT0432
C  IN R.  THE SOLUTION IS RETURNED IN R.  THIS SUBROUTINE CAN BE        ABHT0433
C  REPLACED BY A STANDARD PACKAGE.                                      ABHT0434
C                                                                       ABHT0435
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               ABHT0436
C                                                                       ABHT0437
      DIMENSION A(*),R(*)                                               ABHT0438
      IS = 1                                                            ABHT0439
      NA = ND                                                           ABHT0440
      DO 120 IA = 1, ND                                                 ABHT0441
         NAM = NA-1                                                     ABHT0442
         IC = IS                                                        ABHT0443
         B = A(IS)                                                      ABHT0444
         IF (NAM.EQ.0) GO TO 130                                        ABHT0445
         DO 100 I = 1, NAM                                              ABHT0446
            IC = IC+1                                                   ABHT0447
            A(IC) = A(IC)/B                                             ABHT0448
  100    CONTINUE                                                       ABHT0449
         DO 110 KA = 1, NAM                                             ABHT0450
            C = B*A(IS+KA)                                              ABHT0451
            DO 110 KB = KA, NAM                                         ABHT0452
               IC = IC+1                                                ABHT0453
               A(IC) = A(IC)-C*A(IS+KB)                                 ABHT0454
  110    CONTINUE                                                       ABHT0455
         IS = IS+NA                                                     ABHT0456
         NA = NA-1                                                      ABHT0457
  120 CONTINUE                                                          ABHT0458
  130 IC = 0                                                            ABHT0459
      DO 150 IA = 1, ND                                                 ABHT0460
         C = R(IA)                                                      ABHT0461
         IC = IC+1                                                      ABHT0462
         R(IA) = R(IA)/A(IC)                                            ABHT0463
         IF (IA.EQ.ND) GO TO 150                                        ABHT0464
         IAP = IA+1                                                     ABHT0465
         DO 140 K = IAP, ND                                             ABHT0466
            IC = IC+1                                                   ABHT0467
            R(K) = R(K)-C*A(IC)                                         ABHT0468
  140    CONTINUE                                                       ABHT0469
  150 CONTINUE                                                          ABHT0470
      IC = ND*(ND+1)/2+1                                                ABHT0471
      DO 170 IB = 1, ND                                                 ABHT0472
         IA = ND-IB                                                     ABHT0473
         IC = IC-1                                                      ABHT0474
         IF (IA.EQ.0) GO TO 180                                         ABHT0475
         K = ND                                                         ABHT0476
         DO 160 I = 1, IB                                               ABHT0477
            IC = IC-1                                                   ABHT0478
            R(IA) = R(IA)-A(IC)*R(K)                                    ABHT0479
            K = K-1                                                     ABHT0480
  160    CONTINUE                                                       ABHT0481
  170 CONTINUE                                                          ABHT0482
  180 RETURN                                                            ABHT0483
      END                                                               ABHT0484
      SUBROUTINE RADEQS (PHI,NM,NDIM,ND1,ND2,V,AR,WR,SQRWR,EG,LA,NP,H,  ABHT0485
     *   ALPHA,Z,T,VB,DNUM,Y1,Y2,S1,S2,IMAXA,ENDING)                    ABHT0486
C                                                                       ABHT0487
C  THIS SUBROUTINE CALCULATES THE ENERGY EIGENVALUES AND EIGEN-         ABHT0488
C  FUNCTIONS FOR THE POTENTIAL GIVEN IN ARRAY V.  THE RESULTS           ABHT0489
C  ARE RETURNED IN ARRAYS PHI AND EG.  THE ARRAY IMAXA GIVES THE        ABHT0490
C  MESH POINT AT WHICH ORBITAL N BECOMES LESS THAN 0.0001.              ABHT0491
C                                                                       ABHT0492
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               ABHT0493
C                                                                       ABHT0494
      DIMENSION PHI(ND1,ND2),EG(ND2),LA(ND2),NP(ND2),V(ND1),T(ND1),     ABHT0495
     *   VB(ND1),DNUM(ND1),Y1(ND1),Y2(ND1),S1(ND1),S2(ND1),SQRWR(ND1),  ABHT0496
     *   IMAXA(ND2),AR(ND1),WR(ND1)                                     ABHT0497
      LOGICAL ENDING                                                    ABHT0498
      CHARACTER LBL*6,CH1*2,CH2*1,CH3*1,CH4*2,CH5*8,CH6*6               ABHT0499
      DATA CH1,CH4,CH5,CH6 / 'E(','):','12345678','SPDFGH'/             ABHT0500
      DO 330 N = 1, NM                                                  ABHT0501
         E = EG(N)                                                      ABHT0502
         L = LA(N)                                                      ABHT0503
         NPS = NP(N)                                                    ABHT0504
         EUB = 10.0                                                     ABHT0505
         ELB = -Z*Z/(L+NPS)**2                                          ABHT0506
         ICOU = 0                                                       ABHT0507
         AL = L                                                         ABHT0508
C                                                                       ABHT0509
C  ENTRY POINT FOR NEW ITERATION IN EIGENVALUE CALCULATION              ABHT0510
C                                                                       ABHT0511
  100    ICOU = ICOU+1                                                  ABHT0512
C                                                                       ABHT0513
         IF (ICOU.GT.30) STOP 'CONVERGENCE FAILURE IN SR RADEQS'        ABHT0514
         DO 110 I = 1, NDIM                                             ABHT0515
            PHI(I,N) = 0.0                                              ABHT0516
            R = AR(I)                                                   ABHT0517
            VVV = V(I)-2.0*Z/R-E+AL*(AL+1.0)/R**2                       ABHT0518
            T(I) = VVV                                                  ABHT0519
            VVV = VVV+(0.25+ALPHA*R)/(R*(1.0+ALPHA*R))**2               ABHT0520
            VB(I) = (H*WR(I))**2*VVV                                    ABHT0521
            VB(I) = MIN(VB(I),4.0)                                      ABHT0522
            DNUM(I) = 1.0-VB(I)/12.0                                    ABHT0523
            VB(I) = VB(I)/DNUM(I)                                       ABHT0524
  110    CONTINUE                                                       ABHT0525
C                                                                       ABHT0526
C  DETERMINE NDIMR WHICH IS THE LAST MESH POINT FOR THE CALCULATION     ABHT0527
C  OF THE CURRENT ORBITAL                                               ABHT0528
C                                                                       ABHT0529
         DO 120 I = NDIM, 1, -1                                         ABHT0530
            ITURN = I                                                   ABHT0531
            IF (T(I).LT.0.0) GO TO 130                                  ABHT0532
  120    CONTINUE                                                       ABHT0533
  130    ITURN = ITURN+1                                                ABHT0534
         NDIMR = NDIM                                                   ABHT0535
         IF (ITURN.GE.NDIM) GO TO 150                                   ABHT0536
         SS = 0.0                                                       ABHT0537
         DO 140 I = ITURN, NDIM                                         ABHT0538
            SS = SS+H*SQRT(T(I))*WR(I)                                  ABHT0539
            NDIMR = I                                                   ABHT0540
            IF (SS.GT.10.0) GO TO 150                                   ABHT0541
  140    CONTINUE                                                       ABHT0542
  150    CONTINUE                                                       ABHT0543
C                                                                       ABHT0544
C  INTEGRATE THE RADIAL EQUATION OUTWARD FROM THE ORIGIN AND            ABHT0545
C  COUNT THE NODES IN THE SOLUTION                                      ABHT0546
C                                                                       ABHT0547
         B1 = -Z/(L+1)                                                  ABHT0548
         R = AR(1)                                                      ABHT0549
         Y1(1) = (B1*R+1.0)*R**(L+1)*DNUM(1)/SQRWR(1)                   ABHT0550
         R = AR(2)                                                      ABHT0551
         NNODES = 0                                                     ABHT0552
         Y1(2) = (B1*R+1.0)*R**(L+1)*DNUM(2)/SQRWR(2)                   ABHT0553
         DO 160 I = 3, NDIMR                                            ABHT0554
            Y1(I) = (2.0+VB(I-1))*Y1(I-1)-Y1(I-2)                       ABHT0555
            IF (Y1(I)*Y1(I-1).LT.0.0) NNODES = NNODES+1                 ABHT0556
  160    CONTINUE                                                       ABHT0557
         IF (NNODES.GE.NPS) EUB = MIN(EUB,E)                            ABHT0558
         IF (NNODES.LT.NPS) ELB = MAX(ELB,E)                            ABHT0559
         IF (NNODES.LE.NPS) GO TO 170                                   ABHT0560
         E = (EUB+ELB)/2.0                                              ABHT0561
         GO TO 100                                                      ABHT0562
  170    S1(1) = (Y1(1)*WR(1))**2                                       ABHT0563
         DO 180 I = 2, NDIMR                                            ABHT0564
            S1(I) = S1(I-1)+(Y1(I)*WR(I))**2                            ABHT0565
  180    CONTINUE                                                       ABHT0566
C                                                                       ABHT0567
C  INTEGRATE THE RADIAL EQUATION INWARD                                 ABHT0568
C                                                                       ABHT0569
         AK = SQRT(-E)                                                  ABHT0570
         Y2(NDIMR) = EXP(-AK*AR(NDIMR))                                 ABHT0571
         Y2(NDIMR-1) = EXP(-AK*AR(NDIMR-1))                             ABHT0572
         DO 190 I = (NDIMR-2), 1, -1                                    ABHT0573
            Y2(I) = (2.0+VB(I+1))*Y2(I+1)-Y2(I+2)                       ABHT0574
  190    CONTINUE                                                       ABHT0575
         S2(NDIMR) = 0.0                                                ABHT0576
         DO 200 I = (NDIMR-1), 1, -1                                    ABHT0577
            S2(I) = S2(I+1)+(WR(I+1)*Y2(I+1))**2                        ABHT0578
  200    CONTINUE                                                       ABHT0579
C                                                                       ABHT0580
C  CALCULATE A NEW ESTIMATE FOR THE ENERGY EIGENVALUE BY                ABHT0581
C  FINDING THE OPTIMUM POINT TO MATCH THE INWARD AND OUTWARD            ABHT0582
C  SOLUTIONS.  EXPH IS THE EXPECTATION VALUE OF THE SINGLE              ABHT0583
C  PARTICLE HAMILTONIAN FOR THE WAVE FUNCTION FORMED BY MATCHING        ABHT0584
C  THE INWARD AND OUTWARD SOLUTIONS.                                    ABHT0585
C                                                                       ABHT0586
         ENEW = 10.0                                                    ABHT0587
         IMIN = 5                                                       ABHT0588
         IMAX = NDIMR-5                                                 ABHT0589
         IMID = NDIMR/2                                                 ABHT0590
         W = (Y1(IMID+1)*Y2(IMID)-Y2(IMID+1)*Y1(IMID))/H                ABHT0591
         DO 230 I = IMIN, IMAX                                          ABHT0592
            XX = H*(S1(I)*Y2(I)**2+S2(I)*Y1(I)**2)                      ABHT0593
            EXPH = E+W*Y1(I)*Y2(I)/XX                                   ABHT0594
C                                                                       ABHT0595
C  IF THE ORBITAL IS NOT THE LOWEST FOR THE PARTICULAR L, PROJECT       ABHT0596
C  OUT LOWER STATES OF THE SAME L VALUE                                 ABHT0597
C                                                                       ABHT0598
            SN = 0.0                                                    ABHT0599
            SD = 0.0                                                    ABHT0600
            IF (NPS.EQ.1) GO TO 220                                     ABHT0601
            NQ = N-1                                                    ABHT0602
            DO 210 N1 = 1, NQ                                           ABHT0603
               IF (LA(N1).NE.L) GO TO 210                               ABHT0604
               EPR = EG(N1)                                             ABHT0605
               OLAPSQ = (PHI(I,N1)*W/(EPR-E))**2/(XX*WR(I))             ABHT0606
               SN = SN+EPR*OLAPSQ                                       ABHT0607
               SD = SD+OLAPSQ                                           ABHT0608
  210       CONTINUE                                                    ABHT0609
            IF (SD.GE.1.0) GO TO 230                                    ABHT0610
  220       EXPH = (EXPH-SN)/(1.0-SD)                                   ABHT0611
            IF (EXPH.GT.ENEW) GO TO 230                                 ABHT0612
            ENEW = EXPH                                                 ABHT0613
            MPT = I                                                     ABHT0614
  230    CONTINUE                                                       ABHT0615
         IF (ENEW.GT.EUB) ENEW = 0.5*(EUB+ELB)                          ABHT0616
C                                                                       ABHT0617
C  TEST FOR COMPLETION OF ENERGY CALCULATION.  TE GIVES CHANGE          ABHT0618
C  IN ENERGY IN THE CURRENT ITERATION                                   ABHT0619
C                                                                       ABHT0620
         TE = ABS(E-ENEW)/MAX(1.0,ABS(E))                               ABHT0621
         E = ENEW                                                       ABHT0622
         IF (TE.GT.0.00001) GO TO 100                                   ABHT0623
C                                                                       ABHT0624
C  CALCULATE NORMALIZED ORBITAL                                         ABHT0625
C                                                                       ABHT0626
         C1 = Y2(MPT)                                                   ABHT0627
         C2 = Y1(MPT)                                                   ABHT0628
         DO 240 I = 1, MPT                                              ABHT0629
            T(I) = C1*Y1(I)*SQRWR(I)/DNUM(I)                            ABHT0630
  240    CONTINUE                                                       ABHT0631
         DO 250 I = MPT, NDIMR                                          ABHT0632
            T(I) = C2*Y2(I)*SQRWR(I)/DNUM(I)                            ABHT0633
  250    CONTINUE                                                       ABHT0634
         S = 0.0                                                        ABHT0635
         DO 260 I = 1, NDIMR                                            ABHT0636
            Y1(I) = T(I)**2                                             ABHT0637
  260    CONTINUE                                                       ABHT0638
         CALL XINT (Y1,S1,H,NDIMR,AR,WR,DEFI,ENDING)                    ABHT0639
         CC = SQRT(DEFI)                                                ABHT0640
         DO 270 I = 1, NDIMR                                            ABHT0641
            T(I) = T(I)/CC                                              ABHT0642
  270    CONTINUE                                                       ABHT0643
C                                                                       ABHT0644
C  CHECK THAT THE ORBITAL HAS THE CORRECT NUMBER OF NODES               ABHT0645
C                                                                       ABHT0646
         NNODES = 0                                                     ABHT0647
         DO 280 I = 2, NDIMR                                            ABHT0648
            IF (T(I)*T(I-1).LT.0.0) NNODES = NNODES+1                   ABHT0649
  280    CONTINUE                                                       ABHT0650
         IF (NNODES.EQ.NPS-1) GO TO 290                                 ABHT0651
         IF (NNODES.GE.NPS) E = 0.5*(E+ELB)                             ABHT0652
         IF (NNODES.LT.NPS-1) E = 0.5*(E+EUB)                           ABHT0653
         GO TO 100                                                      ABHT0654
  290    CONTINUE                                                       ABHT0655
C                                                                       ABHT0656
C  FIND RANGE OF THE ORBITAL - IS IT LESS THAN 0.0001?                  ABHT0657
C                                                                       ABHT0658
         DO 300 I = NDIMR, 1, -1                                        ABHT0659
            IMAX = I                                                    ABHT0660
            IF (ABS(T(I)).GT.0.0001) GO TO 310                          ABHT0661
  300    CONTINUE                                                       ABHT0662
  310    CONTINUE                                                       ABHT0663
         IMAXA(N) = IMAX                                                ABHT0664
         DO 320 I = 1, NDIMR                                            ABHT0665
            PHI(I,N) = T(I)                                             ABHT0666
  320    CONTINUE                                                       ABHT0667
         EG(N) = E                                                      ABHT0668
         CH2 = CH5(NPS+L:NPS+L)                                         ABHT0669
         CH3 = CH6(L+1:L+1)                                             ABHT0670
         LBL = CH1//CH2//CH3//CH4                                       ABHT0671
         WRITE(8,6000) LBL,E,NDIMR,ICOU                                
  330 CONTINUE                                                          ABHT0673
      RETURN                                                            ABHT0674
C                                                                       ABHT0675
C                                                                       ABHT0676
C                                                                       ABHT0677
 6000 FORMAT ('    ',A6,F12.5,'   NMAX:',I4,I6,' ITERATIONS')           ABHT0678
C                                                                       ABHT0679
      END                                                               ABHT0680
      SUBROUTINE XINT (X,Y,H,N,AR,WR,DEFI,ENDING)                       ABHT0681
C                                                                       ABHT0682
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               ABHT0683
C                                                                       ABHT0684
      DIMENSION X(*),Y(*)                                               ABHT0685
C                                                                       ABHT0686
C  SUBROUTINE TO CARRY OUT INDEFINITE INTEGRAL NUMERICALLY.             ABHT0687
C  THE FUNCTION TO BE INTEGRATED IS IN ARRAY X, AND THE RESULT          ABHT0688
C  IS IN Y.  H IS DELTA RHO AND N IS THE NUMBER OF MESH POINTS.         ABHT0689
C  IF ENDING=.FALSE. THE TRAPEZOIDAL RULE IS USED.  IF ENDING=          ABHT0690
C  .TRUE. A CUBIC PROCEDURE IS USED.                                    ABHT0691
C  A CORRECTION FOR THE INTEGRAL RMAX TO INFINITY IS ESTIMATED          ABHT0692
C  ASSUMING THE INTEGRAND DECREASES EXPONENTIALLY LIKE EXP(-AK          ABHT0693
C  *R) FOR LARGE R.  THE VALUE OF THE INTEGRAL OVER ALL R IS            ABHT0694
C  RETURNED IN DEFI.                                                    ABHT0695
C                                                                       ABHT0696
      DIMENSION AR(N),WR(N)                                             ABHT0697
      LOGICAL ENDING                                                    ABHT0698
      AK = 1.0                                                          ABHT0699
      IF (X(N)*X(N-1).EQ.0.0) GO TO 100                                 ABT0700
      AK = LOG(ABS(X(N)))-LOG(ABS(X(N-1)))                              ABHT0701
      AK = AK/(AR(N-1)-AR(N))                                           ABHT0702
  100 IF (ENDING) GO TO 120                                             ABHT0703
C                                                                       ABHT0704
C  USE TRAPEZOIDAL RULE                                                 ABHT0705
C                                                                       ABHT0706
      CA = WR(1)*X(1)                                                   ABHT0707

      QQ = 0.5*H                                                        ABHT0708
      Y(1) = 0.5*CA                                                     ABHT0709
      DO 110 I = 2, N                                                   ABHT0710
         CB = WR(I)*X(I)                                                ABHT0711
         Y(I) = Y(I-1)+QQ*(CA+CB)                                       ABHT0712
         CA = CB                                                        ABHT0713
  110 CONTINUE                                                          ABHT0714
      DEFI = Y(N)+X(N)/AK                                               ABHT0715
      RETURN                                                            ABHT0716
C                                                                       ABHT0717
C  USE A CUBIC RULE                                                     ABHT0718
C                                                                       ABHT0719
  120 CA = WR(1)*X(1)                                                   ABHT0720
      CB = WR(2)*X(2)                                                   ABHT0721
      CC = WR(3)*X(3)                                                   ABHT0722
      Y(1) = CA/2.0                                                     ABHT0723
      Y(2) = Y(1)+H*(CA+CB)/2.0                                         ABHT0724
      WA = -H/24.0                                                      ABHT0725
      WB = 13.0*H/24.0                                                  ABHT0726
      NMI = N-1                                                         ABHT0727
      DO 130 I = 3, NMI                                                 ABHT0728
         CD = WR(I+1)*X(I+1)                                            ABHT0729
         Y(I) = Y(I-1)+WA*(CA+CD)+WB*(CB+CC)                            ABHT0730
         CA = CB                                                        ABHT0731
         CB = CC                                                        ABHT0732
         CC = CD                                                        ABHT0733
  130 CONTINUE                                                          ABHT0734
      Y(N) = Y(N-1)+H*(CB+CC)/2.0                                       ABHT0735
      DEFI = Y(N)+X(N)/AK                                               ABHT0736
      RETURN                                                            ABHT0737
      END                                                               ABHT0738
      SUBROUTINE POTRD (V,AR,NDIM)                                      ABHT0739
C                                                                       ABHT0740
C  SUBROUTINE TO READ AN INPUT POTENTIAL GIVEN AT R VALUES IN           ABHT0741
C  ARRAY ARIN AND INTERPOLATE IT TO A POTENTIAL AT R VALUES IN          ABHT0742
C  ARRAY AR.                                                            ABHT0743
C                                                                       ABHT0744
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               ABHT0745
C                                                                       ABHT0746
      DIMENSION V(NDIM),AR(NDIM),VIN(200),ARIN(200)                     ABHT0747
      OPEN (UNIT=1,FILE='INPOT.DAT')                                        
      NIN1 = 200                                                        ABHT0749
      READ (1,6000,END=150) NIN                                         ABHT0750
      IF (NIN.GT.NIN1) GO TO 160                                        ABHT0751
      READ (1,6200,END=150) (ARIN(I),I=1,NIN)                           ABHT0752
      READ (1,6200,END=150) (VIN(I),I=1,NIN)                            ABHT0753
      DO 140 I = 1, NDIM                                                ABHT0754
         RR = AR(I)                                                     ABHT0755
         IF (RR.GT.ARIN(1)) GO TO 100                                   ABHT0756
         V(I) = VIN(1)                                                  ABHT0757
         GO TO 140                                                      ABHT0758
  100    IF (RR.LT.ARIN(NIN)) GO TO 110                                 ABHT0759
         V(I) = VIN(NIN)*ARIN(NIN)/RR                                   ABHT0760
         GO TO 140                                                      ABHT0761
  110    DO 120 KK = 1, NIN                                             ABHT0762
            K = KK                                                      ABHT0763
            IF (ARIN(K).GT.RR) GO TO 130                                ABHT0764
  120    CONTINUE                                                       ABHT0765
  130    IF (K.EQ.2) K = 3                                              ABHT0766
         IF (K.EQ.NIN) K = NIN-1                                        ABHT0767
         A = ARIN(K-2)                                                  ABHT0768
         B = ARIN(K-1)                                                  ABHT0769
         C = ARIN(K)                                                    ABHT0770
         D = ARIN(K+1)                                                  ABHT0771
         VA = VIN(K-2)                                                  ABHT0772
         VB = VIN(K-1)                                                  ABHT0773
         VC = VIN(K)                                                    ABHT0774
         VD = VIN(K+1)                                                  ABHT0775
         VINT = VA*(RR-B)*(RR-C)*(RR-D)/((A-B)*(A-C)*(A-D))+VB*(RR-A)*  ABHT0776
     *      (RR-C)*(RR-D)/((B-A)*(B-C)*(B-D))+VC*(RR-A)*(RR-B)*(RR-D)/((ABHT0777
     *      C-A)*(C-B)*(C-D))+VD*(RR-A)*(RR-B)*(RR-C)/((D-A)*(D-B)*(D-C)ABHT0778
     *      )                                                           ABHT0779
         V(I) = VINT                                                    ABHT0780
  140 CONTINUE                                                          ABHT0781
      RETURN                                                            ABHT0782
C                                                                       ABHT0783
C  END-OF-FILE REACHED ON POTENTIAL INPUT FILE                          ABHT0784
C                                                                       ABHT0785
  150 WRITE(8,6300)                                                    
      STOP                                                              ABHT0787
C                                                                       ABHT0788
  160 WRITE(8,6100)                                                        
      STOP                                                              ABHT0790
C                                                                       ABHT0791
 6000 FORMAT (I5)                                                       ABHT0792
 6100 FORMAT (/,'  STOP - THE DIMENSIONS OF ARRAYS VIN AND ARIN         ABHT0793
     * IN SR POTRD ARE TOO SMALL')                                      ABHT0794
 6200 FORMAT (5G14.6)                                                   ABHT0795
 6300 FORMAT (/,'  STOP--END-OF-FILE REACHED ON POTENTIAL INPUT FILE')  ABHT0796
C                                                                       ABHT0797
      END                                                               ABHT0798
      SUBROUTINE RADMSH (ROMIN,H,ALPHA,NDIM,AR,WR,SQRWR)                ABHT0799
C                                                                       ABHT0800
C  THIS SUBROUTINE SETS UP THE RADIAL MESH POINTS AND THE ARRAYS        ABHT0801
C  WR AND SQRWR.  THE VALUES OF R ARE CALCULATED FROM THE RHO VALUES    ABHT0802
C  USING THE NEWTON-RAPHSON METHOD.                                     ABHT0803
C                                                                       ABHT0804
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)  
	DIMENSION AR(NDIM),WR(NDIM),SQRWR(NDIM)                           
      DO 110 I = 1, NDIM                                                ABHT0806
         RHX = ROMIN+(I-1)*H                                            ABHT0807
         IF (I.EQ.1) R = EXP(ROMIN)                                     ABHT0808
         IF (I.NE.1) R = AR(I-1)                                        ABHT0809
         DO 100 K = 1, 5                                                ABHT0810
            RHOP = LOG(R)+ALPHA*R                                       ABHT0811
C                                                                       ABHT0812
C           RHOP = DLOG(R)+ALPHA*R                                      ABHT0813
C                                                                       ABHT0814
            R = R+(RHX-RHOP)*R/(1.0+ALPHA*R)                            ABHT0815
  100    CONTINUE                                                       ABHT0816
         AR(I) = R                                                      ABHT0817
         WR(I) = R/(1.0+ALPHA*R)                                        ABHT0818
         SQRWR(I) = SQRT(WR(I))                                         ABHT0819
  110 CONTINUE                                                          ABHT0820
      RETURN                                                            ABHT0821
      END                                                               ABHT0822
      SUBROUTINE SECSOL (PHI2,PSI,NDIM,E,Z,L,NPHI,V,ND1,AR,H,ALPHA,Y1,Y2ABHT0823
     *   ,S1,S2,T,VB,WR,SQRWR,DNUM)                                     ABHT0824
C                                                                       ABHT0825
C  CALCULATE A SECOND SOLUTION FOR THE REDUCED RADIAL SCHROEDINGER      ABHT0826
C  EQUATION INWARD AND OUTWARD FROM THE OUTER MAXIMUM OF THE FIRST      ABHT0827
C  SOLUTION AND FIX THE WRONSKIAN OF THE TWO SOLUTIONS TO BE 1          ABHT0828
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                                                  ABHT0829
      DIMENSION PHI2(NDIM),PSI(NDIM),T(ND1),AR(ND1),WR(ND1),SQRWR(ND1), ABHT0830
     *   DNUM(ND1),Y1(ND1),Y2(ND1),V(ND1),VB(ND1),S1(ND1),S2(ND1)       ABHT0831
      LOGICAL ENDING                                                    ABHT0832
      AL = L                                                            ABHT0833
      DO 100 I = 1, NPHI                                                ABHT0834
         R = AR(I)                                                      ABHT0835
         VVV = V(I)-2.0*Z/R-E+AL*(AL+1.0)/R**2                          ABHT0836
         VVV = VVV+(0.25+ALPHA*R)/(R*(1.0+ALPHA*R))**2                  ABHT0837
         VB(I) = (H*WR(I))**2*VVV                                       ABHT0838
         DNUM(I) = 1.0-VB(I)/12.0                                       ABHT0839
         VB(I) = VB(I)/DNUM(I)                                          ABHT0840
  100 CONTINUE                                                          ABHT0841
      AA = 0.0                                                          ABHT0842
      DO 110 I = NPHI, 1, -1                                            ABHT0843
         IK = I                                                         ABHT0844
         CC = ABS(T(I))                                                 ABHT0845
         IF (CC.LT.AA) GO TO 120                                        ABHT0846
         AA = CC                                                        ABHT0847
  110 CONTINUE                                                          ABHT0848
  120 Y1(IK+1) = 1.0                                                    ABHT0849
      Y1(IK) = 0.0                                                      ABHT0850
      DO 130 I = (IK+2), NPHI                                           ABHT0851
         Y1(I) = (2.0+VB(I-1))*Y1(I-1)-Y1(I-2)                          ABHT0852
  130 CONTINUE                                                          ABHT0853
      DO 140 I = (IK-1), 1, -1                                          ABHT0854
         Y1(I) = (2.0+VB(I+1))*Y1(I+1)-Y1(I+2)                          ABHT0855
  140 CONTINUE                                                          ABHT0856
      DO 150 I = 1, NPHI                                                ABHT0857
         Y1(I) = SQRWR(I)*Y1(I)/DNUM(I)                                 ABHT0858
  150 CONTINUE                                                          ABHT0859
      D1 = T(IK+1)-T(IK-1)                                              ABHT0860
      D2 = Y1(IK+1)-Y1(IK-1)                                            ABHT0861
      AA = (T(IK)*D2-Y1(IK)*D1)/(2.0*H*WR(IK))                          ABHT0862
      DO 160 I = 1, NPHI                                                ABHT0863
         PHI2(I) = Y1(I)/AA                                             ABHT0864
  160 CONTINUE                                                          ABHT0865
C                                                                       ABHT0866
C  CALCULATE AUXILIARY FUNCTION PSI (SEE EQ.(14))                       ABHT0867
C                                                                       ABHT0868
      DO 170 I = 1, NPHI                                                ABHT0869
         Y1(I) = T(I)*PHI2(I)                                           ABHT0870
         Y2(I) = T(I)**2                                                ABHT0871
  170 CONTINUE                                                          ABHT0872
      ENDING = .FALSE.                                                  ABHT0873
      CALL XINT (Y1,S1,H,NPHI,AR,WR,DEFI,ENDING)                        ABHT0874
      CALL XINT (Y2,S2,H,NPHI,AR,WR,DEFI,ENDING)                        ABHT0875
      AA = DEFI                                                         ABHT0876
      DO 180 I = 1, NPHI                                                ABHT0877
         Y1(I) = T(I)*S1(I)+PHI2(I)*(AA-S2(I))                          ABHT0878
         Y2(I) = Y1(I)*T(I)                                             ABHT0879
  180 CONTINUE                                                          ABHT0880
      CALL XINT (Y2,S2,H,NPHI,AR,WR,DEFI,ENDING)                        ABHT0881
      CA = DEFI/2.0                                                     ABHT0882
      DO 190 I = 1, NPHI                                                ABHT0883
         PSI(I) = Y1(I)-CA*T(I)                                         ABHT0884
  190 CONTINUE                                                          ABHT0885
      RETURN                                                            ABHT0886
      END                                                               ABHT0887
      SUBROUTINE EXPOT (HA,N,LA,O,NDIM,NPHI,NM,PHI,ND1,ND2,FAC,H,AR,WR, ABHT0888
     *   Y1,Y2,S1,S2,ENDING)                                            ABHT0889
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                                                 ABHT0890
C  SET UP THE ARRAY HA - ESSENTIALLY THE EXCHANGE POTENTIAL OF          ABHT0891
C  THE HARTREE-FOCK EQUATION (SEE EQ.(7))                               ABHT0892
C                                                                            
                                                            
      DIMENSION HA(ND1),AR(ND1),WR(ND1),Y1(ND1),Y2(ND1),S1(ND1),S2(ND1),
     *   LA(20),O(20),FAC(0:20),PHI(ND1,ND2)                            
      LOGICAL ENDING                                                    
      L = LA(N)                                                         
      DO 100 I = 1, NDIM                                                
         HA(I) = 0.0                                                    ABHT0899
  100 CONTINUE                                                          ABHT0900
      DO 130 NS = 1, NM                                                 ABHT0901
         LP = LA(NS)                                                    ABHT0902
         LCMAX = L+LP                                                   ABHT0903
         LCMIN = ABS(L-LP)                                              ABHT0904
         DO 130 LC = LCMIN, LCMAX, 2                                    ABHT0905
            LD = (LC+L+LP)/2                                            ABHT0906

            TA = FAC(LD-L)*FAC(LD-LP)*FAC(LD-LC)/FAC(LD)                ABHT0907
            LD = LC+L+LP                                                ABHT0908
            TB = FAC(LD-2*L)*FAC(LD-2*LP)*FAC(LD-2*LC)/FAC(LD+1)        ABHT0909
            THREEJ = TB/TA**2                                           ABHT0910
            VV = 2*(2*L+1)                                              ABHT0911
            IF (NS.NE.N) WAZ = 0.5*O(N)*O(NS)*THREEJ                    ABHT0912
            IF ((NS.EQ.N).AND.(LC.NE.0)) WAZ = 0.5*THREEJ*O(N)*(O(N)-1.0ABHT0913
     *         )*VV/(VV-1.0)                                            ABHT0914
            IF ((NS.EQ.N).AND.(LC.EQ.0)) WAZ = O(N)                     ABHT0915
            DO 110 I = 1, NPHI                                          ABHT0916
               AA = PHI(I,N)*PHI(I,NS)                                  ABHT0917
               AB = AR(I)**(LC+1)                                       ABHT0918
               Y1(I) = AA*AB/AR(I)                                      ABHT0919
               Y2(I) = AA/AB                                            ABHT0920
  110       CONTINUE                                                    ABHT0921
            CALL XINT (Y1,S1,H,NPHI,AR,WR,DEFI,ENDING)                  ABHT0922
            CALL XINT (Y2,S2,H,NPHI,AR,WR,DEFI,ENDING)                  ABHT0923
            CA = DEFI                                                   ABHT0924
            DO 120 I = 1, NPHI                                          ABHT0925
               AB = AR(I)**(LC+1)                                       ABHT0926
               XY = S1(I)/AB+(CA-S2(I))*AB/AR(I)                        ABHT0927
               XY = 2.0*WAZ*XY                                          ABHT0928
               HA(I) = HA(I)+XY*PHI(I,NS)                               ABHT0929
  120       CONTINUE                                                    ABHT0930
  130 CONTINUE                                                          ABHT0931
      RETURN                                                            ABHT0932
      END                                                               ABHT0933
