C 22 FEB 06 - LVS - IMPLEMENT A STABLE OPTIMIZATION ALGORITHM
C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA MMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 13 FEB 05 - MWS - PAD COMMON BLOCK HERMIT, WERMIT, NSHEL, ROOT
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 26 MAR 03 - MWS - WRITE SCREEN PARAMETERS TO -IP-, NOT -IPIRC-
C 14 JAN 03 - MWS - FIX BUG READING DMA FILES
C 12 DEC 02 - MWS - ALLOW IFTTYP/IFTFIX INPUT SELECTION
C  7 AUG 02 - MWS - NEW MODULE FOR SCREENING TERM FOR CHARGE PENETRATION
C
C
C*MODULE CHGPEN  *DECK CGPX
      SUBROUTINE CGPX
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SCREEN
C
      PARAMETER (MXATM=2000,MXEFC=200,NMAX=50)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /PDCPAR/ CENTER(3),DPOLE(3),QPOLE(6),RMAX,DELR,CONSTR,
     *                PTSEL,VDWSCL,PTDENS,VDWINC,NFREQ,LAYER,NPDC
      COMMON /STNBUF/ STNPNT(4,2*MXATM),BIGEXP,NPTSTN,NBUFFM
C
      DIMENSION ALPINI(NPTSTN),COEFINI(NPTSTN),NLAYER(LAYER+1)
C
      DATA CHECK/8HCHECK   /
      DATA GEODES /8HGEODESIC/
C
C              ---- CHARGE PENETRATION SCREENING ----
C          M.A.FREITAG, M.S.GORDON, J.H.JENSEN, W.H.STEVENS
C                 J.CHEM.PHYS. 112, 7300-7306(2000)
C     THE CLASSICAL ELECTROSTATIC POTENTIAL DUE TO MULTIPOLES IS
C     SCREENED BY A DAMPING TERM SUCH THAT IT PRODUCES A BETTER
C     FIT TO THE QUANTUM MECHANICAL ELECTROSTATIC POTENTIAL.
C
      NEFC = NPTSTN
C
      IPIRC=4
C
      IF(GOPARR) THEN
         IF(MASWRK) WRITE(IW,*) 'SORRY, SCREENING FIT IS SERIAL ONLY'
         RETURN
      END IF
C
C         --- READ INPUT CONTROLLING GENERATION OF THIS FIT ---
C     IF THERE IS NO SUCH INPUT THE USER MUST WANT TO SKIP THIS TERM.
C     ESTIMATE THE NUMBER OF MULTIPOLE EXPANSION POINTS AS BEING
C     AT EVERY ATOM AND AT EVERY ATOM-PAIR CENTER
C
      CALL VALFM(LOADFM)
      LEFNAM = LOADFM + 1
      LEFPOS = LEFNAM + NEFC
      LEFC   = LEFPOS + NEFC*3
      LEFD   = LEFC   + NEFC
      LEFQ   = LEFD   + NEFC*3
      LEFO   = LEFQ   + NEFC*6
      LIEQ   = LEFO   + NEFC*10
      LAST = LIEQ   + NEFC
      NEEDV  = LAST - LOADFM - 1
      CALL GETFM(NEEDV)
C
      CALL CGPINP(SCREEN,IPIRC,X(LEFNAM),X(LEFPOS),X(LEFC),X(LEFD),
     *            X(LEFQ),X(LEFO),X(LIEQ))
      IF(.NOT.SCREEN) GO TO 888
C
C         SAVE ORIGINAL VALUES OF ALP & COEF
      DO I=1,NEFC
         ALPINI(I)=ALP(I)
         COEFINI(I)=COEF(I)
      ENDDO
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
C
C          PURE STONE ANALYSIS REQUIRES SAMPLING POINTS CLOSER TO
C          THE NUCLEI THAN IF NUMERICAL STONE ANALYSIS WAS USED.
C
      IF(VDWSCL.EQ.0.0D+00) THEN
         IF(BIGEXP.EQ.0.0D+00) THEN
            VDWSCL=0.5D+00
         ELSE
            VDWSCL=0.8D+00
         END IF
      END IF
C
C         DO CYCLE OVER DIFFERENT VALUES OF RMIN HERE
C
      DO 666 M=1,NRUN
         WRITE(IW,499) PTSEL
         WRITE(IW,501) VDWSCL
C         
C          GRID WILL CONTAIN NDPC POINTS, NPDC<MAXPNT
C
      MAXPNT=NPDC
      CALL VALFM(LOADFM)
      IXYZ    = LOADFM + 1
      LRAD    = IXYZ   + 3*MAXPNT
      LCICOS  = LRAD   + NEFC
      LASTPDC = LCICOS + 3*2562
      NEEDPDC = LASTPDC - LOADFM - 1
      WRITE(IW,*)' MEMORY FOR PDCPTS: ',NEEDPDC
      CALL GETFM(NEEDPDC)
C      
      IF(EXETYP.EQ.CHECK) GO TO 654
      CALL PDCPTS(X(IXYZ),X(LRAD),X(LCICOS),MAXPNT,NEFC,NLAYER)
      KOUNT = NPDC
      NPDC   = MAXPNT
 654  CONTINUE  
      CALL RETFM(NEEDPDC)
C
C        NOW GRID CONTAINS KOUNT POINTS
C
      NEW_KOUNT = 0
      IF(PTSEL.EQ.GEODES) THEN
         DO ILAYER=1,LAYER
            SCALED = (VDWSCL+VDWINC*(ILAYER-1))/VDWSCL
            KOUNT_L = INT(SCALED*SCALED*
     *                   (NLAYER(ILAYER+1)-NLAYER(ILAYER)))
            NEW_KOUNT = NEW_KOUNT + KOUNT_L
         END DO
         WRITE(IW,502) NEW_KOUNT
      ELSE
         NEW_KOUNT = KOUNT
      END IF
C
      NVX = NV   
C
C        X,Y,Z = GRID COORDINATES
C        QUANT = AB INITIO ELECTROSTATIC POTENTIAL ON GRID
C        CLASS = UNSCREENED MULTIPOLAR ELECTROSTATIC POTENTIAL ON GRID
C        DIFF  = QUANT - CLASS
C        FITTD = SCREENED MULTIPOLAR ELECTROSTATIC POTENTIAL
C        RAW   = UNSCREENED MULTIPOLAR POTENTIAL (VESTIGE OF ISSCHG=1)
C        NEED THESE ARRAYS FOR FITTING PROCEDURE
C
      CALL VALFM(LOADFM)
      IXYZ   = LOADFM +1
      LX     = IXYZ   + 3*KOUNT
      LY     = LX     + KOUNT
      LZ     = LY     + KOUNT
      LQUANT = LZ     + KOUNT
      LCLASS = LQUANT + KOUNT
      LDIFF  = LCLASS + KOUNT
      LFITTD = LDIFF  + KOUNT
      LRAW   = LFITTD + KOUNT
      LP     = LRAW   + KOUNT
      LXI    = LP     + NVX
      LPT    = LXI    + NVX*NVX
      LPPT   = LPT    + NVX
      LXIT   = LPPT   + NVX
      LTMP   = LXIT   + NVX
      LASTPW = LTMP   + NVX
      NEEDPW   = LASTPW - LOADFM - 1
      IF(MASWRK) WRITE(IW,9000) NEEDPW
      CALL GETFM(NEEDPW)
C
      DO I=1,KOUNT
         X(LX+I-1)= X(IXYZ+3*(I-1))
         X(LY+I-1)= X(IXYZ+3*(I-1)+1)
         X(LZ+I-1)= X(IXYZ+3*(I-1)+2)
      END DO
C
      CALL VALFM(LOADFM)
      LDEN   = LOADFM + 1
      LBETA  = LDEN   + L2
      LG     = LBETA  + L2
      LAST1  = LG     + 225*1000
      NEED1  = LAST1 - LOADFM - 1
C
      IF(MASWRK) WRITE(IW,9001) NEED1
      CALL GETFM(NEED1)
      IF(EXETYP.EQ.CHECK) GO TO 700
C
C        LOAD TOTAL DENSITY MATRIX FIRST
C
      CALL DENDD1(X(LDEN),X(LBETA),L2)
C
      IF(MASWRK) WRITE(IW,9010) KOUNT
      CALL FLSHBF(IW)
C
C     ----- COMPUTE QUANTUM ELECTROSTATIC POTENTIAL ON THE GRID -----
C
      CALL TSECND(TIM0)
      CALL CGPQUA(X(LQUANT),X(LX),X(LY),X(LZ),KOUNT,X(LDEN),L2,X(LG))
      CALL TSECND(TIM1)
      TQUANT = TIM1-TIM0
      TIM0 = TIM1
      IF(MASWRK) WRITE(IW,9020)
      CALL FLSHBF(IW)
C
C       MEMERY CHECK
 700  CONTINUE
      CALL RETFM(NEED1)
C
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     ----- COMPUTE CLASSICAL ELECTROSTATIC POTENTIAL ON THE GRID -----
C
      CALL CGPCLS(X(LCLASS),X(LX),X(LY),X(LZ),KOUNT,
     *            X(LEFPOS),X(LEFC),X(LEFD),X(LEFQ),X(LEFO),NEFC)
      CALL TSECND(TIM1)
      TCLASS = TIM1-TIM0
      TIM0 = TIM1
      IF(MASWRK) WRITE(IW,9030) TQUANT,TCLASS
      CALL FLSHBF(IW)
C
C         FORM DIFFERENCE POTENTIAL = QUANTUM - CLASSICAL
C
      CALL VSUB(X(LCLASS),1,X(LQUANT),1,X(LDIFF),1,KOUNT)
C
C        SELECT EXPONENTIAL/GAUSSIAN
C        TO GET THE EXPONENT VALUES SETTLED DOWN, THE FIRST FIT IS
C        ALWAYS BE REQUIRED TO HAVE LINEAR COEFFICIENTS.
C      
      NFIT  = 0
 100  CONTINUE
      NFIT = NFIT+1
C
      IERRF = IFTTYP(NFIT)
      ICFIX = IFTFIX(NFIT)
      IF(IERRF.NE.0  .AND.  IERRF.NE.2) GO TO 820
C
      IF(MASWRK) THEN
         IF(IERRF.EQ.0) WRITE(IW,75)
         IF(IERRF.EQ.2) WRITE(IW,74)
         IF(IFTFIX(NFIT).NE.1) WRITE(IW,76)
      END IF
C
C         RESTORE ORIGINAL ALPHA VALUES FOR MULTIPLE RUNS  
C         GO BACK TO ORIGINAL GUESS VALUES, FOR BOTH FITS,
C         IN CASE THE FIRST FIT BECAME CRAQZY ALONG THE WAY.
C
C--   IF(NFIT.EQ.1) THEN
         DO I=1,NEFC
             ALP(I)  = ALPINI(I)
             COEF(I) = COEFINI(I)
         END DO
C--   END IF   
C      
C     ---- OPTIMIZE SCREENING EXPONENTS -----
C
      IF(NV.NE.0) THEN
         IF(MASWRK) WRITE(IW,9040)
C
C        NUMBER OF FUNCTION EVALUATIONS   
      NF=0
C      
C        CALCULATE INITIAL VALUES OF OPTIMIZED PARAMETERS 
      CALL CALC_ALP(X(LP),X(LIEQ))
C
C        CALCULATE INITIAL SET OF DIRECTIONS XI
      CALL VCLR(X(LXI),1,NVX*NVX)    
      DO I=1,NVX
         X(LXI+(I-1)*NVX+I-1) = 0.1D+00
      END DO
C
C        CONVERGENCE CRITERIUM
      FTOL = 1.0D-06
C
      CALL POWELL(X(LP),X(LXI),NVX,FTOL,ITER,FRET,
     *            X(LPT),X(LPPT),X(LXIT),
     *            X(LX),X(LY),X(LZ),X(LQUANT),X(LFITTD),X(LDIFF),KOUNT,
     *            X(LTMP),NLAYER,X(LEFPOS),X(LEFC),X(LEFD),X(LEFQ),
     *            X(LEFO))
C
         CALL TSECND(TIM1)
         TFIT = TIM1-TIM0
         TIM0 = TIM1
         IF(MASWRK) WRITE(IW,9050) TFIT
         CALL FLSHBF(IW)
      END IF
C
C     --- GENERATE THE FINAL FITTED POTENTIAL, PRINT/PUNCH SUMMARY ---
C
      CALL CGPFIT(X(LX),X(LY),X(LZ),X(LFITTD),KOUNT,
     *            X(LEFPOS),X(LEFC),X(LEFD),X(LEFQ),X(LEFO))
C
      IF(IERRF.EQ.0) WRITE(IP,11) VDWSCL
      IF(IERRF.EQ.2) WRITE(IP,12) VDWSCL
      CALL CGPPRT(X(LX),X(LY),X(LZ),X(LQUANT),X(LFITTD),X(LCLASS),
     *            X(LEFNAM),X(LEFPOS),KOUNT,NEW_KOUNT)
      WRITE(IP,31)
C
C        NOW GO BACK AND DO THE SECOND FIT
C
      IF(NFIT.LT.2) GO TO 100
C
  820 CONTINUE
      IF(MASWRK) THEN
         WRITE(IW,9200)
         CALL TIMIT(1)
      END IF
C
  800 CONTINUE
C
      CALL RETFM(NEEDPW)
C
C           THIS ENDS THE CYCLE OF VDWSCL
        IF(NRUN .GT. 1) VDWSCL=VDWSCL+0.1D+00
  666 CONTINUE
C
  888 CONTINUE
      CALL RETFM(NEEDV)
      RETURN
C
   74 FORMAT(/' *** SINGLE EXPONENTIAL USED AS DAMPING FUNCTION ***')
   75 FORMAT(/' *** SINGLE GAUSSIAN USED AS DAMPING FUNCTION ***')
   76 FORMAT(/1X,'INITIAL FIT WILL RESTRICT LINEAR COEFS TO UNITY,'/
     *        1X,'THIS RESTRICTION WILL BE REMOVED IN THE SECOND FIT.')
   11 FORMAT('SCREEN       (FROM VDWSCL=',F8.3,')')
   12 FORMAT('SCREEN2      (FROM VDWSCL=',F8.3,')')
   31 FORMAT('STOP')
  499 FORMAT(/1X,'QUANTUM ELECTROSTATIC POTENTIAL GRID IS',
     *              ' GENERATED BY $PDC PTSEL=',A8)
  501 FORMAT(1X,'POINT SELECTION IS USING VDWSCL= ',F6.3)
  502 FORMAT(/1X,' TO SIMULATE EVEN POINT DENSITY IN GEODES GRID'/
     *        1X,' TOTAL NUMBER OF POINTS WAS ARTIFICIALLY INCREASED'/
     *        1X,' TO ',I10)   
 9000 FORMAT(1X,'DAMPING COMPUTATION NEEDS',I12,' WORDS')
 9001 FORMAT(1X,'PLUS ',I12,' WORDS')
 9010 FORMAT(/1X,'----- THE ELECTROSTATIC POTENTIAL GRID CONTAINS',I9,
     *          ' POINTS -----')
 9020 FORMAT(1X,'----- THE QUANTUM POTENTIAL HAS BEEN GENERATED',
     *          ' ON THE GRID -----')
 9030 FORMAT(1X,'----- THE UNSCREENED CLASSICAL POTENTIAL GRID HAS',
     *          ' BEEN GENERATED -----'/
     *       1X,'ELECTROSTATIC POTENTIAL CPU TIME: QUANTUM=',F8.1,
     *          ', CLASSICAL=',F8.2)
 9040 FORMAT(/1X,'FITTING SCREENED MULTIPOLE CLASSICAL POTENTIAL',
     *          ' TO THE QUANTUM POTENTIAL...')
 9050 FORMAT(1X,'POWELL PROCEDURE REQUIRED',F8.1,' SECONDS')
 9200 FORMAT(1X,'..... DONE WITH CHARGE PENETRATION SCREENING .....')
      END
C
C*MODULE CHGPEN  *DECK CGPINP
      SUBROUTINE CGPINP(SCREEN,IPIRC,EFNAM,EFPOS,EFC,EFD,EFQ,EFO,IEQ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SCREEN
C
      CHARACTER*8 STR8,WORD2,WORD3
      CHARACTER*8  COEFS(0:1)
      CHARACTER*11 TYPES(0:3)
C
      DIMENSION EFNAM(NEFC),EFPOS(3,NEFC),EFC(NEFC),
     *          EFD(3,NEFC),EFQ(6,NEFC),EFO(10,NEFC),
     *          IEQ(NEFC),MAIN_ALP(NEFC)
C
      PARAMETER (NMAX=50, MXEFC=200)
      PARAMETER (MXATM=2000)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     DATA VALUES IN /CGPPAR/ ARE AS FOLLOWS:
C         ISCCHG = 1 IF ONLY ATOMS BUT NOT BONDS ARE TO BE SCREENED.
C         THRSH  = CONTROLS PRINTING OF LARGE ERRORS BETWEEN AB INITIO
C                  AND SCREENED MULTIPOLES ON THE GRID
C         MAXITS = NON-ZERO CONTROLS THE NUMBER OF MAXIMUM POWELL
C                  STEPS IN THE "STEPIT" ROUTINE.
C         ICFIX  = 1 TO FREEZE THE LINEAR COEFFICIENTS AT 1.0.
C                  IF SCREEN=[ 1.0 - A*EXP(-B*R**2) ], THEN SCREEN
C                  DOESN'T = 0.0 AT THE ORIGIN UNLESS A=1.0.  THIS
C                  IS RARELY USED FOR EXPONENTIALS, BUT IS THE
C                  TRADITIONAL CHOICE FOR GAUSSIANS.  IT MAY BE
C                  VERY SENSIBLE TO FIX GAUSSIAN COEFS TO UNITY,
C                  HOWEVER, ESPECIALLY IF ANY OF THE LINEAR COEFS
C                  GO TO NEGATIVE NUMBERS.
C                  THIS ISN'T OPERATIVE FOR ERROR FUNCTION SCREENING.
C         IERRF  = 0 FOR GAUSSIAN FUNCTION SCREENING
C                = 1 FOR ERROR FUNCTION SCREENING
C                = 2 FOR EXPONENTIAL SCREENING
C                = 3 SEE CODE FOR THE FORMULA
C         IFTTYP AND IFTFIX STORE THE IERRF/ICFIX DATA FOR THE TWO FITS.
C         NEFC   = NUMBER OF SCREENED MULTIPOLE CENTERS
C
C     FOR THE READING OF $DAMP NAMELIST
C
      PARAMETER (NNAM=8)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA DAMP/8HDAMP    /
      DATA QNAM/8HISCCHG  ,8HTHRSH   ,8HMAXIT   ,8HIFTTYP  ,
     *          8HIFTFIX  ,8HNMAIN   ,8HNRUN    ,8HSCREEN  /
      DATA KQNAM/1,3,1,21,  21,1,1,0/
C
      DATA BLANK/8H        /
C
      DATA ZERO,ONE,TWO,THREE,FIVE/0.0D+00,1.0D+00,2.0D+00,3.0D+00,
     *                            5.0D+00/
      DATA TYPES/'   GAUSSIAN', 'ERROR FUNCT', 'EXPONENTIAL',
     *           '  SEE CODE'/
      DATA COEFS/'VARIABLE','   FIXED'/
C
C       MAIN DRIVER TO FIT THE CHARGE PENETRATION SCREENING
C
C SET DEFAULT VALUES
C
      ISCCHG = 0
      THRSH  = 1.0D+02
      MAXITS = 30
      NMAIN  = NAT
      SCREEN = .TRUE.
C
C        FIRST PASS IS AN EXPONENTIAL FIT, FOR EFP2 SCREENING.
C        SECOND PASS IS A GAUSSIAN FIT, FOR THE ORIGINAL EFP1 SCREENING.
C
C        THE GAUSSIAN FIT WILL BE USED FOR AB INITIO-FRAGMENT SCREENING,
C        SEE EQUATION (4) OF J.PHYS.CHEM. 105, 293-307(2001), WHILE
C        THE EXPONENTIAL FIT WILL BE USED FOR FRAGMENT-FRAGMENT CHARGE
C        PENETRATION SCREENING, SEE EQUATION (28).  NOTE THAT THESE
C        HAVE A SINGLE TERM, AND (28) LACKS A LINEAR COEFFICIENT.
C
      IFTTYP(1)=2
      IFTTYP(2)=0
      IFTFIX(1)=1
      IFTFIX(2)=1
C
C             -NRUN- IS AN UNDOCUMENTED DEVELOPMENT PARAMETER.
C         IT DEFINES THE NUMBER OF RUNS TO BE MADE, STARTING FROM
C         DIFFERENT -VDWSCL- VALUES (GIVEN IN $PDC).  EACH RUN WILL
C         INCREMENT THE MINIMUM RADIUS (VDWSCL) BY 0.1 ANGSTROMS.
C
      NRUN = 1
C
      DO 40 J=1,NEFC
         EFNAM(J)=BLANK
         EFC(J)=ZERO
         DO I=1,3
            EFPOS(I,J)=ZERO
            EFD(I,J)=ZERO
         ENDDO
         DO I=1,10
            EFO(I,J)=ZERO
         ENDDO
         DO I=1,6
            EFQ(I,J)=ZERO
         ENDDO
40    CONTINUE
C
      CALL ERRT
C
      JRET = 0
      CALL NAMEIO(IR,JRET,DAMP,NNAM,QNAM,KQNAM,
     *            ISCCHG,THRSH,MAXITS,IFTTYP,IFTFIX,NMAIN,NRUN,SCREEN,
     *            0,          0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *            0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *            0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
      IF(JRET.EQ.2) THEN
         IF(MASWRK) WRITE(IW,900)
         CALL ABRT
      END IF
C
      IF(MASWRK) THEN
         WRITE(IW,70)
         WRITE(IW,77) SCREEN,ISCCHG,THRSH,MAXITS,NMAIN
         WRITE(IW,80) TYPES(IFTTYP(1)),COEFS(IFTFIX(1))
         IF(IFTTYP(2).GE.0)
     *   WRITE(IW,81) TYPES(IFTTYP(2)),COEFS(IFTFIX(2))
         IF(ISCCHG.EQ.1) WRITE(IW,85)
      END IF
C
      IF(.NOT.SCREEN) RETURN
C
C     ----- READ IN DMA MULTIPOLE INFORMATION PREVIOUSLY COMPUTED -----
C
      IEOF  = 0
      IERR  = 0
      STR8  = '        '
      KSIZE = -8
      CALL OPNCRD(IPIRC,-IW)
      CALL SEQREW(IPIRC)
C
   86 CONTINUE
      CALL RDCARD('CGP CTR1',IEOF)
      IF(IEOF.NE.0) GO TO 1000
      CALL GSTRNG(STR8,KSIZE)
C
      IF(STR8.NE.'RUNTYP  ') GO TO 86
C
C---- SKIP THREE LINES
C
      CALL RDCARD('SKIPPING',IEOF)
      CALL RDCARD('SKIPPING',IEOF)
      CALL RDCARD('SKIPPING',IEOF)
C
C   USE DENSITY BASED MULTIPOLAR EXPANSION INSTEAD 
C
C      IF(STR8.NE.'DENSITY ') GO TO 86
C      CALL RDCARD('SKIPPING',IEOF)
C
C---- GET DMA EXPANSION CENTERS
C
      IF(MASWRK) WRITE(IW,105)
      DO IEFC=1,NEFC
         CALL RDCARD('CGP CTR2',IEOF)
         CALL GSTRNG(STR8,KSIZE)
         READ(UNIT=STR8,FMT='(A8)') EFNAM(IEFC)
         EFPOS(1,IEFC) = RFIND('X CENTER',IERR)
         EFPOS(2,IEFC) = RFIND('Y CENTER',IERR)
         EFPOS(3,IEFC) = RFIND('Z CENTER',IERR)
         IF(MASWRK) WRITE(IW,370) EFNAM(IEFC),(EFPOS(J,IEFC),J=1,3)
      END DO
C
C---- GET CHARGES (ASSUME THEY ARE IN THE INPUT FILE)
C
      CALL RDCARD('SKIPPING',IEOF)
      CALL RDCARD('SKIPPING',IEOF)
      IF(MASWRK) WRITE(IW,145)
      DO IEFC=1,NEFC
         CALL RDCARD('CGP CHRG',IEOF)
         CALL GSTRNG(STR8,KSIZE)
         EFC(IEFC) = RFIND('AT. CHG.',IERR)
         IF(MASWRK) WRITE(IW,370) STR8,EFC(IEFC)
      ENDDO
C
C---- GET DIPOLES (DON'T ASSUME THEY ARE IN THE INPUT FILE)
C
      CALL SEQREW(IPIRC)
  160 CONTINUE
      CALL RDCARD('CGP DIP1',IEOF)
      IF(IEOF.NE.0) THEN
         IF(MASWRK) WRITE(IW,170)
         GO TO 200
      END IF
      CALL GSTRNG(STR8,KSIZE)
      IF(STR8.NE.'DIPOLES') GO TO 160
C
      IF(MASWRK) WRITE(IW,180)
      DO 190 IEFC=1,NEFC
         CALL RDCARD('CGP DIP2',IEOF)
         CALL GSTRNG(STR8,KSIZE)
         DO I=1,3
            EFD(I,IEFC) = RFIND('DIPOLE  ',IERR)
         ENDDO
         IF(MASWRK) WRITE(IW,370) STR8,(EFD(J,IEFC),J=1,3)
  190 CONTINUE
C
C---- GET QUADRUPOLES (DON'T ASSUME THEY ARE IN THE INPUT FILE)
C
  200 CONTINUE
      CALL SEQREW(IPIRC)
  210 CONTINUE
      CALL RDCARD('CGP QUA1',IEOF)
      IF(IEOF.NE.0) THEN
         IF(MASWRK) WRITE(IW,220)
         GO TO 300
      END IF
      CALL GSTRNG(STR8,KSIZE)
      IF(STR8.NE.'QUADRUPO') GO TO 210
C
      IF(MASWRK) WRITE(IW,230)
      DO 290 IEFC=1,NEFC
         CALL RDCARD('CGP QUA2',IEOF)
         CALL GSTRNG(STR8,KSIZE)
         DO I=1,4
            EFQ(I,IEFC) = RFIND('QUADRPLE',IERR)
         ENDDO
         CALL RDCARD('CGP QUA3',IEOF)
         DO I=5,6
            EFQ(I,IEFC) = RFIND('QUADRPLE',IERR)
         ENDDO
         IF(MASWRK) WRITE(IW,370) STR8,(EFQ(J,IEFC),J=1,6)
290   CONTINUE
C
C---- CONVERT QUADRUPOLES
C
      DO 260 NEFX=1,NEFC
         XX=EFQ(1,NEFX)
         YY=EFQ(2,NEFX)
         ZZ=EFQ(3,NEFX)
         XY=EFQ(4,NEFX)
         XZ=EFQ(5,NEFX)
         YZ=EFQ(6,NEFX)
         DUM=XX+YY+ZZ
         EFQ(1,NEFX)=THREE*XX-DUM
         EFQ(2,NEFX)=THREE*YY-DUM
         EFQ(3,NEFX)=THREE*ZZ-DUM
         EFQ(4,NEFX)=THREE*XY
         EFQ(5,NEFX)=THREE*XZ
         EFQ(6,NEFX)=THREE*YZ
  260 CONTINUE
C
C---- GET OCTUPOLES (DON'T ASSUME THEY ARE IN THE INPUT FILE)
C
  300 CONTINUE
      CALL SEQREW(IPIRC)
  310 CONTINUE
      CALL RDCARD('CGP OCT1',IEOF)
      IF(IEOF.NE.0) THEN
         IF(MASWRK) WRITE(IW,320)
         GO TO 392
      END IF
      CALL GSTRNG(STR8,KSIZE)
      IF(STR8.NE.'OCTUPOLE') GO TO 310
C
      IF(MASWRK) WRITE(IW,330)
      DO 350 IEFC=1,NEFC
         CALL RDCARD('CGP OCT2',IEOF)
         CALL GSTRNG(STR8,KSIZE)
         DO K=1,4
            EFO(K,IEFC) = RFIND('OCTUPOLE',IERR)
         ENDDO
         CALL RDCARD('CGP OCT3',IEOF)
         DO K=5,8
            EFO(K,IEFC) = RFIND('OCTUPOLE',IERR)
         ENDDO
         DO K=9,10
            EFO(K,IEFC) = RFIND('OCTUPOLE',IERR)
         ENDDO
         IF(MASWRK) WRITE(IW,380) STR8,(EFO(K,IEFC),K=1,10)
  350 CONTINUE
C
C---- CONVERT OCTUPOLES
C
      DO 390 NEFX=1,NEFC
         XXX=EFO(1,NEFX)
         YYY=EFO(2,NEFX)
         ZZZ=EFO(3,NEFX)
         XXY=EFO(4,NEFX)
         XXZ=EFO(5,NEFX)
         XYY=EFO(6,NEFX)
         YYZ=EFO(7,NEFX)
         XZZ=EFO(8,NEFX)
         YZZ=EFO(9,NEFX)
         XYZ=EFO(10,NEFX)
         DUMX=XXX+XYY+XZZ
         DUMY=XXY+YYY+YZZ
         DUMZ=XXZ+YYZ+ZZZ
         EFO(1,NEFX)=FIVE*XXX-THREE*DUMX
         EFO(2,NEFX)=FIVE*YYY-THREE*DUMY
         EFO(3,NEFX)=FIVE*ZZZ-THREE*DUMZ
         EFO(4,NEFX)=FIVE*XXY-DUMY
         EFO(5,NEFX)=FIVE*XXZ-DUMZ
         EFO(6,NEFX)=FIVE*XYY-DUMX
         EFO(7,NEFX)=FIVE*YYZ-DUMZ
         EFO(8,NEFX)=FIVE*XZZ-DUMX
         EFO(9,NEFX)=FIVE*YZZ-DUMY
         EFO(10,NEFX)=FIVE*XYZ
  390 CONTINUE
C
  392 CONTINUE
C
C         WE ARE NOW DONE READING MULTIPOLE DATA,
C         RETURN INPUT SCANNER TO NORMAL CARD INPUT FILE.
C
C  400 CONTINUE
      CALL SEQREW(IPIRC)
      CALL OPNCRD(IR,-IW)
C
C------------------------------------------------
C READ IN INITIAL GUESS SCREENING DATA IN $DAMPGS
C------------------------------------------------
C
      DO J=1,NEFC
         ISETA(J)=J
         COEF(J)=ONE
         IEQ(J)=J
         IF(J.LE.NMAIN) THEN
            MAIN_ALP(J) = 1
         ELSE
            MAIN_ALP(J) = 0
            IF(ISCCHG.EQ.1) COEF(J) = ZERO
         END IF
      END DO
C         
C          DEFAULT VALUES FOR ALPHAS 
C
      ALPNEW = TWO
      DO I=1,NEFC
         IF (MAIN_ALP(I).EQ.1) THEN 
            ALP(I) = ALPNEW
         ELSE   
            IF(ISCCHG.EQ.1) THEN
               ALP(I) = 10.0D+00
            ELSE
               ALP(I) = 2.0D+00*ALPNEW
            END IF
         ENDIF
      END DO
C
      IF(ISCCHG.EQ.1) THEN 
         NEFCC = NMAIN
      ELSE 
         NEFCC = NEFC
      END IF
C         NUMBER OF EQUIVALENT PARAMETERS      
      NEQUIV=0
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $DAMPGS',IEOF)
      IF(IEOF.NE.0) THEN
         IF(MASWRK) WRITE(IW,505)
         GOTO 16
      END IF
C
 12   CONTINUE
      CALL RDCARD('DAMPGS 1',IEOF)
      IF(IEOF.NE.0) CALL ABRT
C
      WORD2='        '
      LEN=-8
      CALL GSTRNG(WORD2,LEN)
      IF(WORD2.EQ.'$END    ') GOTO 16
C
      READ(UNIT=WORD2,FMT='(A8)') WORD
      DO 470 IEFC=1,NEFCC
         IF(WORD.EQ.EFNAM(IEFC)) THEN
            WORD3='        '
            LGSTR=-8
            CALL GSTRNG(WORD3,LGSTR)
            IF(WORD3.EQ.'        ') THEN
               CALL REREAD
               IEQ(IEFC) = IEFC
C                 READ THE TERMS DEFINING THE SCREENING FUCNTION
C                 COEF = THE LINEAR COEFFICIENT
C                 ALP  = THE EXPONENT
               CALL RDCARD('DAMPGS 2',IEOF)
               COEF(IEFC) = RFIND('COEF    ',IERR)
               ALP(IEFC)  = RFIND('ALP     ',IERR)
            ELSE
               READ(UNIT=WORD3,FMT='(A8)') WORD1
               IEQ(IEFC)=-1
               DO J=1,NEFCC
                  IF(WORD1.EQ.EFNAM(J)) IEQ(IEFC)=J
               ENDDO
               IF(IEQ(IEFC).EQ.-1) THEN
                  IF(MASWRK) WRITE(IW,9002) WORD3,IEFC
                  CALL ABRT
               END IF
               NEQUIV=NEQUIV+1
            END IF
C
            GO TO 12
         END IF
 470  CONTINUE
C
C        IF WE ARE HERE, SOMETHING IS WRONG...
      IF(MASWRK) WRITE(IW,9003) WORD2 
      CALL ABRT
C
C        ALL INPUT IS READ
  16  CONTINUE        
C
      NV = NEFCC-NEQUIV
      WRITE(IW,430) NV
C
      IF(MASWRK) WRITE(IW,435)
C
C  NALP COUNTS THE NUMBER OF VARIABLE EXPONENTS
      NALP=0
C
C ISETA WAS DEFINED SO THAT ONLY THE "ACTIVE" EXPONENTS ARE PASSED
C TO POWELL.  ISETA MAPS THE CURRENT EXPONENT TO THE NEXT ACTIVE
C EXPONENT.  THERE ARE "NALP" ACTIVE EXPONENTS.
C
      NALP = 0
      DO I=1,NEFCC
         IF (IEQ(I) .EQ. I) THEN
            NALP=NALP+1
            ISETA(I)=NALP
            IF(MASWRK) WRITE(IW,440) EFNAM(I),COEF(I),ALP(I)
         ELSE
            ISETA(I)=ISETA(IEQ(I))
            ALP(I)  = ALP(IEQ(I))
            COEF(I) = COEF(IEQ(I))
            IF(MASWRK) WRITE(IW,455) EFNAM(I),EFNAM(IEQ(I))
         END IF
      ENDDO
C
      RETURN
C
 1000 CONTINUE
      IF(MASWRK) WRITE(IW,1010)
      CALL ABRT
      RETURN
C
  900 FORMAT(1X,'ERROR IN $DAMP NAMELIST INPUT - STOP')
   70 FORMAT(/'     ---------------------------------------'/
     *        '     EFFECTIVE POTENTIAL SCREENING PARAMETER'/
     *        '              OPTIMIZATION PROGRAM'/
     *        '     ---------------------------------------'/
     *        '     AUTHORS: WALT STEVENS, MARK FREITAG,',
     *         ' LYUDMILA SLIPCHENKO')
   77 FORMAT(/1X,'THE $DAMP INPUT DATA FOR THIS RUN IS:'/
     *        1X,'SCREEN = ',L5/
     *        1X,'ISCCHG = ',I5/
     *        1X,'THRSH  = ',F5.1/
     *        1X,'MAXIT  = ',I5/
     *        1X,'NMAIN  = ',I5)
   80 FORMAT(1X,'INITIAL FIT WILL USE ',A11,' FORM WITH ',A8,
     *          ' LINEAR COEFICIENTS')
   81 FORMAT(1X,' SECOND FIT WILL USE ',A11,' FORM WITH ',A8,
     *          ' LINEAR COEFICIENTS')
   85 FORMAT(/1X,'*******************************************',
     *       /1X,'NOTE: ONLY THE ATOMIC CENTERS ARE SCREENED!',
     *       /1X,'*******************************************')
  105 FORMAT(/1X,'THE CLASSICAL MULTIPOLE EXPANSION TO BE SCREENED IS'//
     *        1X,'CENTERS OF THE DISTRIBUTED MULTIPOLE EXPANSION,',
     *           ' X,Y,Z')
  145 FORMAT(/1X,'CHARGES             Q')
  170 FORMAT(/1X,'NO DIPOLES FOUND')
  180 FORMAT(/1X,'DIPOLES             X         Y         Z')
  220 FORMAT(/1X,'NO QUADRUPOLES FOUND')
  230 FORMAT(/1X,'QUADRUPOLES        XX        YY        ZZ        XY',
     *           '        XZ        YZ')
  320 FORMAT(/1X,'NO OCTUPOLES FOUND')
  330 FORMAT(/1X,'OCTUPOLES')
  370 FORMAT(5X,A8,6F10.5)
  380 FORMAT(5X,A8,5F10.5/13X,5F10.5)
  430 FORMAT(/1X,'NUMBER OF OPTIMIZED PARAMETERS IS ',I3)
  435 FORMAT(/1X,'INITIAL VALUES OF SCREENING PARAMETERS ARE:'/
     *        1X,'NAME        COEF       EXP')
  440 FORMAT(1X,A8,2F10.5,I6)
  455 FORMAT(1X,A8,' = ',A8)
  505 FORMAT(//25X,'*****   WARNING   *****'/
     *        1X,'NO $DAMPGS CHARGE PENETRATION SCREENING',
     *           ' INPUT WAS FOUND, CONTINUING...'/
     *        3X,'THIS INPUT GROUP LETS YOU ASSIGN EXACTLY EQUAL',
     *           ' SCREENING EXPONENTS'/
     *        3X,'TO ANY SYMMETRY EQUIVALENT ATOMS AND BOND CENTERS,',
     *           ' AND IT IS USUALLY'/
     *        3X,'A VERY GOOD IDEA TO ENTER THIS DATA.'/)
 1010 FORMAT(///'NO GAMESS DMA INPUT FOUND')
 9002 FORMAT(1X,'*** ERROR ***'/1X,'EQUIVALENT POINT NAME ',A8,
     *          ' FOR NEXT SCREENING POINT',I5/
     *          1X,'WAS NOT DEFINED EARLIER IN $DAMPGS')
 9003 FORMAT(1X,'*** ERROR ***'/1X,'POINT NAME ',A8,
     *          ' IS NOT RECOGNIZED')
      END
C
C
C*MODULE CHGPEN  *DECK CGPQUA
      SUBROUTINE CGPQUA(POT,CCX,CCY,CCZ,KOUNT,DAO,L2,G)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL NORM,DOUBLE,IANDJ
C
      DIMENSION POT(KOUNT),CCX(KOUNT),CCY(KOUNT),CCZ(KOUNT),
     *          DAO(L2),G(225,1000)
      DIMENSION XIN(250),YIN(250),ZIN(250),
     *          DIJ(225),FIJ(225),IJX(225),IJY(225),IJZ(225),
     *          IX(35),IY(35),IZ(35),JX(35),JY(35),JZ(35)
C
      PARAMETER (MXGTOT=20000, MXSH=5000, MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /XYZSTV/ XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,
     *                CX,CY,CZ
C
      DATA ZERO,ONE,TWO/0.0D+00,1.0D+00,2.0D+00/
      DATA PI212 /1.1283791670955D+00/
      DATA SQRT3 /1.73205080756888D+00/
      DATA SQRT5 /2.23606797749979D+00/
      DATA SQRT7 /2.64575131106459D+00/
C
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     1          3, 0, 0, 2, 2, 1, 0, 1, 0, 1,
     2          4, 0, 0, 3, 3, 1, 0, 1, 0, 2,
     3          2, 0, 2, 1, 1/
      DATA IX / 1, 6, 1, 1,11, 1, 1, 6, 6, 1,
     1         16, 1, 1,11,11, 6, 1, 6, 1, 6,
     2         21, 1, 1,16,16, 6, 1, 6, 1,11,
     3         11, 1,11, 6, 6/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     1          0, 3, 0, 1, 0, 2, 2, 0, 1, 1,
     2          0, 4, 0, 1, 0, 3, 3, 0, 1, 2,
     3          0, 2, 1, 2, 1/
      DATA IY / 1, 1, 6, 1, 1,11, 1, 6, 1, 6,
     1          1,16, 1, 6, 1,11,11, 1, 6, 6,
     2          1,21, 1, 6, 1,16,16, 1, 6,11,
     3          1,11, 6,11, 6/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     1          0, 0, 3, 0, 1, 0, 1, 2, 2, 1,
     2          0, 0, 4, 0, 1, 0, 1, 3, 3, 0,
     3          2, 2, 1, 1, 2/
      DATA IZ / 1, 1, 1, 6, 1, 1,11, 1, 6, 6,
     1          1, 1,16, 1, 6, 1, 6,11,11, 6,
     2          1, 1,21, 1, 6, 1, 6,16,16, 1,
     3         11,11, 6, 6,11/
C
      ICMIN=1
      ICMAX=1000
      IF(ICMAX.GT.KOUNT) ICMAX=KOUNT
      DO IC=1,KOUNT
         POT(IC)=ZERO
      ENDDO
C
C---- PROGRAM LOOPS BACK TO HERE IN ORDER TO CALCULATE ELECTROSTATIC
C---- POTENTIAL ON GRID POINTS IN BLOCKS OF 1000 POINTS.
C
    1 CONTINUE
      NUMIC=ICMAX-ICMIN+1
      NORM=.TRUE.
C
C     ----- I SHELL
C
      DO 9000 II=1,NSHELL
      I=KATOM(II)
      XI=C(1,I)
      YI=C(2,I)
      ZI=C(3,I)
      I1=KSTART(II)
      I2=I1+KNG(II)-1
      LIT=KTYPE(II)
      MINI=KMIN(II)
      MAXI=KMAX(II)
      LOCI=KLOC(II)-MINI
C
C     ----- J SHELL
C
      DO 8000 JJ=1,II
      J=KATOM(JJ)
      XJ=C(1,J)
      YJ=C(2,J)
      ZJ=C(3,J)
      J1=KSTART(JJ)
      J2=J1+KNG(JJ)-1
      LJT=KTYPE(JJ)
      MINJ=KMIN(JJ)
      MAXJ=KMAX(JJ)
      LOCJ=KLOC(JJ)-MINJ
      NROOTS=(LIT+LJT-2)/2+1
      RR=(XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
      IANDJ=II.EQ.JJ
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
      IJ=0
      MAX=MAXJ
      DO 50 I=MINI,MAXI
      NX=IX(I)
      NY=IY(I)
      NZ=IZ(I)
      IF(IANDJ) MAX=I
      DO 50 J=MINJ,MAX
      IJ=IJ+1
      IJX(IJ)=NX+JX(J)
      IJY(IJ)=NY+JY(J)
      IJZ(IJ)=NZ+JZ(J)
   50 CONTINUE
      DO 65 IC=1,NUMIC
      DO 60 I=1,IJ
      G(I,IC)=ZERO
   60 CONTINUE
   65 CONTINUE
C
C     ----- I PRIMITIVE
C
      JGMAX=J2
      DO 7000 IG=I1,I2
      AI=EX(IG)
      ARRI=AI*RR
      AXI=AI*XI
      AYI=AI*YI
      AZI=AI*ZI
      CSI=CS(IG)
      CPI=CP(IG)
      CDI=CD(IG)
      CFI=CF(IG)
      CGI=CG(IG)
C
C     ----- J PRIMITIVE
C
      IF(IANDJ) JGMAX=IG
      DO 6000 JG=J1,JGMAX
      AJ=EX(JG)
      AA=AI+AJ
      AA1=ONE/AA
      DUM=AJ*ARRI*AA1
C     IF(DUM.GT.TOL) GO TO 6000
      FAC=EXP(-DUM)
      CSJ=CS(JG)
      CPJ=CP(JG)
      CDJ=CD(JG)
      CFJ=CF(JG)
      CGJ=CG(JG)
      AX=(AXI+AJ*XJ)*AA1
      AY=(AYI+AJ*YJ)*AA1
      AZ=(AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
      DOUBLE=IANDJ.AND.IG.NE.JG
      MAX=MAXJ
      NN=0
      DO 310 I=MINI,MAXI
      GO TO ( 70, 80,180,180, 90,180,180,100,180,180,
     1       110,180,180,120,180,180,180,180,180,130,
     2       140,180,180,150,180,180,180,180,180,160,
     3       180,180,170,180,180),I
   70 DUM1=CSI*FAC
      GO TO 180
   80 DUM1=CPI*FAC
      GO TO 180
   90 DUM1=CDI*FAC
      GO TO 180
  100 IF(NORM) DUM1=DUM1*SQRT3
      GO TO 180
  110 DUM1=CFI*FAC
      GO TO 180
  120 IF(NORM) DUM1=DUM1*SQRT5
      GO TO 180
  130 IF(NORM) DUM1=DUM1*SQRT3
      GO TO 180
  140 DUM1=CGI*FAC
      GO TO 180
  150 IF(NORM) DUM1=DUM1*SQRT7
      GO TO 180
  160 IF(NORM) DUM1=DUM1*SQRT5/SQRT3
      GO TO 180
  170 IF(NORM) DUM1=DUM1*SQRT3
  180 IF(IANDJ) MAX=I
      DO 310 J=MINJ,MAX
      GO TO (190,200,300,300,210,300,300,220,300,300,
     1       230,300,300,240,300,300,300,300,300,250,
     2       260,300,300,270,300,300,300,300,300,280,
     3       300,300,290,300,300),J
  190 DUM2=DUM1*CSJ
      IF(.NOT.DOUBLE) GO TO 300
      IF(I.GT.1) GO TO 195
      DUM2=DUM2+DUM2
      GO TO 300
  195 DUM2=DUM2+CSI*CPJ*FAC
      GO TO 300
  200 DUM2=DUM1*CPJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 300
  210 DUM2=DUM1*CDJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 300
  220 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 300
  230 DUM2=DUM1*CFJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 300
  240 IF(NORM) DUM2=DUM2*SQRT5
      GO TO 300
  250 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 300
  260 DUM2=DUM1*CGJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 300
  270 IF(NORM) DUM2=DUM2*SQRT7
      GO TO 300
  280 IF(NORM) DUM2=DUM2*SQRT5/SQRT3
      GO TO 300
  290 IF(NORM) DUM2=DUM2*SQRT3
  300 NN=NN+1
  310 DIJ(NN)=DUM2
C
C     ----- NUCLEAR ATTRACTION
C
      DUM=PI212*AA1
      DO 400 I=1,IJ
  400 FIJ(I)=DIJ(I)*DUM
      AAX=AA*AX
      AAY=AA*AY
      AAZ=AA*AZ
      ICX=0
      DO 450 IC=ICMIN,ICMAX
      ICX=ICX+1
      CX=CCX(IC)
      CY=CCY(IC)
      CZ=CCZ(IC)
      XX=AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
C
      IF (NROOTS.LE.3) CALL RT123
      IF (NROOTS.EQ.4) CALL ROOT4
      IF (NROOTS.EQ.5) CALL ROOT5
C
      MM=0
      DO 420 K=1,NROOTS
      UU=AA*U(K)
C     WW=W(K)*ZNUC
      WW=W(K)
      TT=ONE/(AA+UU)
      T= SQRT(TT)
      X0=(AAX+UU*CX)*TT
      Y0=(AAY+UU*CY)*TT
      Z0=(AAZ+UU*CZ)*TT
      IN=-5+MM
      DO 410 I=1,LIT
      IN=IN+5
      NI=I
      DO 410 J=1,LJT
      JN=IN+J
      NJ=J
      CALL STVXYZ
      XIN(JN)=XINT
      YIN(JN)=YINT
      ZIN(JN)=ZINT*WW
  410 CONTINUE
  420 MM=MM+25
      DO 440 I=1,IJ
      NX=IJX(I)
      NY=IJY(I)
      NZ=IJZ(I)
      DUM=ZERO
      MM=0
      DO 430 K=1,NROOTS
      NXM=NX+MM
      NYM=NY+MM
      NZM=NZ+MM
      DUM=DUM+XIN(NXM)*YIN(NYM)*ZIN(NZM)
  430 MM=MM+25
  440 G(I,ICX)=G(I,ICX)+DUM*FIJ(I)
  450 CONTINUE
C
C           END OF PRIMITIVE LOOPS
C
6000  CONTINUE
7000  CONTINUE
      ICX=0
      DO 7600 IC=ICMIN,ICMAX
         ICX=ICX+1
         MAX=MAXJ
         NN=0
         DO 7500 I=MINI,MAXI
            LI=LOCI+I
            IN=(LI*(LI-1))/2
            IF(IANDJ) MAX=I
            FACT=TWO
            DO 7400 J=MINJ,MAX
               LJ=LOCJ+J
               JN=LJ+IN
               NN=NN+1
               IF(IANDJ.AND.I.EQ.J) FACT=ONE
               POT(IC)=POT(IC)-G(NN,ICX)*DAO(JN)*FACT
               FACT=TWO
 7400       CONTINUE
 7500    CONTINUE
 7600 CONTINUE
C         END OF SHELL LOOPS
 8000 CONTINUE
 9000 CONTINUE
      IF(ICMAX.EQ.KOUNT) RETURN
      ICMIN=ICMIN+1000
      ICMAX=ICMAX+1000
      IF(ICMAX.GT.KOUNT) ICMAX=KOUNT
      GO TO 1
      END
C
C*MODULE CHGPEN  *ECK STVXYZ
      SUBROUTINE STVXYZ
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /XYZSTV/ XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,
     *                CX,CY,CZ
      COMMON /HERMIT/ H1,H2(2),H3(3),H4(4),H5(5),H6(6),H7(7),
     *                H8(8),H9(9),H10(10)
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
C
      DIMENSION H(28),W(28),MIN(7),MAX(7)
C
      EQUIVALENCE (H(1),H1),(W(1),W1)
C
      DATA MIN /1,2,4, 7,11,16,22/
      DATA MAX /1,3,6,10,15,21,28/
      DATA ZERO /0.0D+00/
C
C     ----- GAUSS-HERMITE QUADRATURE USING MINIMUM POINT FORMULA -----
C
      XINT=ZERO
      YINT=ZERO
      ZINT=ZERO
      NPTS=(NI+NJ-2)/2+1
      IMIN=MIN(NPTS)
      IMAX=MAX(NPTS)
      DO 13 I=IMIN,IMAX
      DUM=W(I)
      PX=DUM
      PY=DUM
      PZ=DUM
      DUM=H(I)*T
      PTX=DUM+X0
      PTY=DUM+Y0
      PTZ=DUM+Z0
      AX=PTX-XI
      AY=PTY-YI
      AZ=PTZ-ZI
      BX=PTX-XJ
      BY=PTY-YJ
      BZ=PTZ-ZJ
      GO TO (5,4,3,2,1),NI
    1 PX=PX*AX
      PY=PY*AY
      PZ=PZ*AZ
    2 PX=PX*AX
      PY=PY*AY
      PZ=PZ*AZ
    3 PX=PX*AX
      PY=PY*AY
      PZ=PZ*AZ
    4 PX=PX*AX
      PY=PY*AY
      PZ=PZ*AZ
    5 GO TO (12,11,10,9,8,7,6),NJ
    6 PX=PX*BX
      PY=PY*BY
      PZ=PZ*BZ
    7 PX=PX*BX
      PY=PY*BY
      PZ=PZ*BZ
    8 PX=PX*BX
      PY=PY*BY
      PZ=PZ*BZ
    9 PX=PX*BX
      PY=PY*BY
      PZ=PZ*BZ
   10 PX=PX*BX
      PY=PY*BY
      PZ=PZ*BZ
   11 PX=PX*BX
      PY=PY*BY
      PZ=PZ*BZ
   12 CONTINUE
      XINT=XINT+PX
      YINT=YINT+PY
      ZINT=ZINT+PZ
   13 CONTINUE
      RETURN
      END
C*MODULE CHGPEN  *DECK CGPCLS
      SUBROUTINE CGPCLS(CLASS,X,Y,Z,KOUNT,EFPOS,EFC,EFD,EFQ,EFO,NEFC)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION CLASS(KOUNT),X(KOUNT),Y(KOUNT),Z(KOUNT),EFPOS(3,NEFC),
     *          EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),EFO(10,NEFC)
C
      PARAMETER (ZERO=0.0D+00)
C
C     THE CLASSICAL POTENTIAL IS CALCUALATED USING THE MULTIPOLES
C
      DO 350 K=1,KOUNT
         VAL=ZERO
         DO IEFC=1,NEFC
            VAL=VAL
     *         +CGPMUL(K,IEFC,X,Y,Z,KOUNT,EFPOS,EFC,EFD,EFQ,EFO,NEFC)
         ENDDO
         CLASS(K)=VAL
  350 CONTINUE
      RETURN
      END
C*MODULE CHGPEN  *DECK CGPMUL
      DOUBLE PRECISION FUNCTION CGPMUL(K,IPT,X,Y,Z,KOUNT,
     *                                 EFPOS,EFC,EFD,EFQ,EFO,NEFC)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),EFPOS(3,NEFC),
     *          EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),EFO(10,NEFC)
C
      DATA PT5/0.5D+00/
      DATA ZERO,TWO/0.0D+00,2.0D+00/
      DATA THREE,FIVE,SIX/3.0D+00,5.0D+00,6.0D+00/
C
C       GENERATE ELECTROSTATIC POTENTIAL FOR THE DMA EXPANSION
C
      VAL=ZERO
      XI=X(K)-EFPOS(1,IPT)
      YI=Y(K)-EFPOS(2,IPT)
      ZI=Z(K)-EFPOS(3,IPT)
      RR=XI*XI+YI*YI+ZI*ZI
      R1=SQRT(RR)
      R3=R1*RR
      R5=RR*RR*R1
      R7=RR*R5
C
      VALO=EFC(IPT)/R1
      VAL=VAL+VALO
C
      TERM=ZERO
      TERM=TERM+EFD(1,IPT)*XI
      TERM=TERM+EFD(2,IPT)*YI
      TERM=TERM+EFD(3,IPT)*ZI
      VALO=TERM/R3
      VAL=VAL+VALO
C
      XX=EFQ(1,IPT)
      YY=EFQ(2,IPT)
      ZZ=EFQ(3,IPT)
      XY=EFQ(4,IPT)
      XZ=EFQ(5,IPT)
      YZ=EFQ(6,IPT)
      XXI=XI*XI
      YYI=YI*YI
      ZZI=ZI*ZI
      XYI=XI*YI
      XZI=XI*ZI
      YZI=YI*ZI
      TERM=ZERO
      TERM=TERM+XX*(THREE*XXI-RR)
      TERM=TERM+YY*(THREE*YYI-RR)
      TERM=TERM+ZZ*(THREE*ZZI-RR)
      TERM=TERM+XY*THREE*XYI*TWO
      TERM=TERM+XZ*THREE*XZI*TWO
      TERM=TERM+YZ*THREE*YZI*TWO
      VALO=PT5*TERM/(THREE*R5)
      VAL=VAL+VALO
C
      XXX=EFO(1,IPT)
      YYY=EFO(2,IPT)
      ZZZ=EFO(3,IPT)
      XXY=EFO(4,IPT)
      XXZ=EFO(5,IPT)
      XYY=EFO(6,IPT)
      YYZ=EFO(7,IPT)
      XZZ=EFO(8,IPT)
      YZZ=EFO(9,IPT)
      XYZ=EFO(10,IPT)
      XXXI=XI*XI*XI
      YYYI=YI*YI*YI
      ZZZI=ZI*ZI*ZI
      XXYI=XI*XI*YI
      XXZI=XI*XI*ZI
      XYYI=XI*YI*YI
      YYZI=YI*YI*ZI
      XZZI=XI*ZI*ZI
      YZZI=YI*ZI*ZI
      XYZI=XI*YI*ZI
      RRXI=RR*XI
      RRYI=RR*YI
      RRZI=RR*ZI
      TERM=ZERO
      TERM=TERM+XXX*(FIVE*XXXI-THREE*RRXI)
      TERM=TERM+YYY*(FIVE*YYYI-THREE*RRYI)
      TERM=TERM+ZZZ*(FIVE*ZZZI-THREE*RRZI)
      TERM=TERM+XXY*(FIVE*XXYI-RRYI)*THREE
      TERM=TERM+XXZ*(FIVE*XXZI-RRZI)*THREE
      TERM=TERM+XYY*(FIVE*XYYI-RRXI)*THREE
      TERM=TERM+YYZ*(FIVE*YYZI-RRZI)*THREE
      TERM=TERM+XZZ*(FIVE*XZZI-RRXI)*THREE
      TERM=TERM+YZZ*(FIVE*YZZI-RRYI)*THREE
      TERM=TERM+XYZ*(FIVE*XYZI)*SIX
      VALO=PT5*TERM/(FIVE*R7)
      VAL=VAL+VALO
      CGPMUL=VAL
      RETURN
      END
C
C*MODULE CHGPEN  *DECK CGPFIT
      SUBROUTINE CGPFIT(X,Y,Z,FITTED,KOUNT,
     *                  EFPOS,EFC,EFD,EFQ,EFO)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     DIMENSION RAW(KOUNT)
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),FITTED(KOUNT),
     *          EFPOS(3,NEFC),EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),
     *          EFO(10,NEFC)
C
      PARAMETER (NMAX=50,MXEFC=200)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
C
      DATA PT5/0.5D+00/
      DATA ZERO,ONE,TWO/0.0D+00,1.0D+00,2.0D+00/
      DATA THREE,FIVE,SIX/3.0D+00,5.0D+00,6.0D+00/
C
C LOOP OVER ALL GRID POINTS
C
      DO 420 K=1,KOUNT
      VAL=ZERO
C      RAWVAL=ZERO
C
C LOOP OVER ALL MULTIPOLE POINTS
C
      DO 410 IPT=1,NEFC
      XI=X(K)-EFPOS(1,IPT)
      YI=Y(K)-EFPOS(2,IPT)
      ZI=Z(K)-EFPOS(3,IPT)
      RR=XI*XI+YI*YI+ZI*ZI
      R1=SQRT(RR)
      R3=R1*RR
      R5=RR*RR*R1
      R7=RR*R5
C
      XDAMP=ZERO
C LOOP OVER ALL TERMS IN DAMPING FUNCTION
      IF(ALP(IPT).EQ.ZERO) GO TO 405
C
C COMPUTE PROPER DAMPING FUNCTION
C
      ALPRR=ALP(IPT)*RR
      IF(IERRF.EQ.0) XDAMP=XDAMP+COEF(IPT)*EXP(-ALPRR)
      IF(IERRF.EQ.1) XDAMP=XDAMP+COEF(IPT)*ERRF(SQRT(ALPRR))
      IF(IERRF.EQ.2) XDAMP=XDAMP+COEF(IPT)*EXP(-ALP(IPT)*SQRT(RR))
      IF(IERRF.EQ.3) XDAMP=XDAMP+COEF(IPT)*EXP(-TWO*ALP(IPT)
     *                       *SQRT(RR))*(ONE+ALP(IPT)*SQRT(RR))
405   CONTINUE
      IF(IERRF.EQ.0) XDAMP=ONE-XDAMP
      IF(IERRF.EQ.2) XDAMP=ONE-XDAMP
      IF(IERRF.EQ.3) XDAMP=ONE-XDAMP
C
C CALCULATE THE MULTIPOLAR ELECTROSTATIC POTENTIAL
C
      VALO=EFC(IPT)/R1
      VAL=VAL+VALO*XDAMP
C      RAWVAL=RAWVAL+VALO
C
      TERM=ZERO
      TERM=TERM+EFD(1,IPT)*XI
      TERM=TERM+EFD(2,IPT)*YI
      TERM=TERM+EFD(3,IPT)*ZI
      VALO=TERM/R3
      VAL=VAL+VALO*XDAMP
C      RAWVAL=RAWVAL+VALO
C
      XX=EFQ(1,IPT)
      YY=EFQ(2,IPT)
      ZZ=EFQ(3,IPT)
      XY=EFQ(4,IPT)
      XZ=EFQ(5,IPT)
      YZ=EFQ(6,IPT)
      XXI=XI*XI
      YYI=YI*YI
      ZZI=ZI*ZI
      XYI=XI*YI
      XZI=XI*ZI
      YZI=YI*ZI
      TERM=ZERO
      TERM=TERM+XX*(THREE*XXI-RR)
      TERM=TERM+YY*(THREE*YYI-RR)
      TERM=TERM+ZZ*(THREE*ZZI-RR)
      TERM=TERM+XY*THREE*XYI*TWO
      TERM=TERM+XZ*THREE*XZI*TWO
      TERM=TERM+YZ*THREE*YZI*TWO
      VALO=PT5*TERM/(THREE*R5)
      VAL=VAL+VALO*XDAMP
C      RAWVAL=RAWVAL+VALO
C
      XXX=EFO(1,IPT)
      YYY=EFO(2,IPT)
      ZZZ=EFO(3,IPT)
      XXY=EFO(4,IPT)
      XXZ=EFO(5,IPT)
      XYY=EFO(6,IPT)
      YYZ=EFO(7,IPT)
      XZZ=EFO(8,IPT)
      YZZ=EFO(9,IPT)
      XYZ=EFO(10,IPT)
      XXXI=XI*XI*XI
      YYYI=YI*YI*YI
      ZZZI=ZI*ZI*ZI
      XXYI=XI*XI*YI
      XXZI=XI*XI*ZI
      XYYI=XI*YI*YI
      YYZI=YI*YI*ZI
      XZZI=XI*ZI*ZI
      YZZI=YI*ZI*ZI
      XYZI=XI*YI*ZI
      RRXI=RR*XI
      RRYI=RR*YI
      RRZI=RR*ZI
      TERM=ZERO
      TERM=TERM+XXX*(FIVE*XXXI-THREE*RRXI)
      TERM=TERM+YYY*(FIVE*YYYI-THREE*RRYI)
      TERM=TERM+ZZZ*(FIVE*ZZZI-THREE*RRZI)
      TERM=TERM+XXY*(FIVE*XXYI-RRYI)*THREE
      TERM=TERM+XXZ*(FIVE*XXZI-RRZI)*THREE
      TERM=TERM+XYY*(FIVE*XYYI-RRXI)*THREE
      TERM=TERM+YYZ*(FIVE*YYZI-RRZI)*THREE
      TERM=TERM+XZZ*(FIVE*XZZI-RRXI)*THREE
      TERM=TERM+YZZ*(FIVE*YZZI-RRYI)*THREE
      TERM=TERM+XYZ*(FIVE*XYZI)*SIX
      VALO=PT5*TERM/(FIVE*R7)
C APPLY DAMPING
      VAL=VAL+VALO*XDAMP
C      RAWVAL=RAWVAL+VALO
410   CONTINUE
      FITTED(K)=VAL
C      RAW(K)=RAWVAL
C--      IF(MOD(K,1000).EQ.0) WRITE(6,415) K,FITTED(K)
C--  415 FORMAT(1X,'CGPFIT:',I8,F12.5)
420   CONTINUE
      RETURN
      END
C
C*MODULE CHGPEN  *DECK CGPINV
      SUBROUTINE CGPINV(A,N,B,NMAX,M,DET,INDEX1,INDEX2,IPIVOT,PIVOT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION INDEX1(NMAX),INDEX2(NMAX),IPIVOT(NMAX),A(NMAX*NMAX),
     *          B(NMAX)
C
      EQUIVALENCE (IROW,JROW),(ICOLUM,JCOLUM),(AMAX,T,SWAP),(NL1,NROW)
C
C   SOLUTION OF MATRIX EQUATION AX=B  OBTAIN A(-1) & X
C     A     FIRST ADDRESS OF A  -  A(-1) UPON RETURN
C               A IS DESTROYED DURING SOLUTION
C     N     ORDER OF A AT THIS ENTRANCE
C     B     VECTORS B   -  LOCATION OF X UPON RETURN
C               B IS DESTROYED DURING SOLUTION
C     NMAX  EXTERNAL DIMENSION OF A
C     M     NUMBER OF COLUMN VECTORS IN B
C     DET   COMPUTED VALUE OF THE DETERMINANT OF A
C     INDEX1, INDEX2, & IPIVOT   INTEGER WORK ARRAYS NMAX IN LENGTH
C     PIVOT ON WORD WORK SPACE
C
C     INITIALIZATION
C
      DET = 1.0D+00
      DO 20 J=1,N
   20 IPIVOT(J)=0
      DO 550 I=1,N
C
C     SEARCH FOR PIVOT ELEMENT
C
      AMAX = 0.0D+00
      DO 105 J=1,N
      IF (IPIVOT(J)-1) 60,105,60
   60 DO 100 K=1,N
      IF (IPIVOT(K)-1) 71,100,740
   71 LL=(K-1)*NMAX+J
      IF ( ABS(AMAX)- ABS(A(LL))) 85,100,100
   85 IROW=J
      ICOLUM=K
      AMAX=A (LL)
  100 CONTINUE
  105 CONTINUE
      IF (AMAX) 110,111,110
  111 DET=0.0D+00
      RETURN
  110 IPIVOT(ICOLUM)=IPIVOT(ICOLUM)+1
C
C     INTERCHANGE ROWS TO PUT ELEMENT ON DIAGONAL
C
      IF (IROW-ICOLUM) 140,260,140
  140 DET=-DET
      DO 200 L=1,N
      JCON=(L-1)*NMAX
      NROW=JCON+IROW
      NCOL=JCON+ICOLUM
      SWAP=A(NROW)
      A(NROW)=A(NCOL)
  200 A(NCOL)=SWAP
      IF (M) 260,260,210
  210 DO 250 L=1,M
      JCON=(L-1)*NMAX
      NROW=JCON+IROW
      NCOL=JCON+ICOLUM
      SWAP=B(NROW)
      B(NROW)=B(NCOL)
  250 B(NCOL)=SWAP
  260 INDEX1(I)=IROW
      INDEX2(I)=ICOLUM
      LL=(ICOLUM-1)*NMAX+ICOLUM
      PIVOT=A(LL)
      DET=DET*PIVOT
C
C     DIVIDE PIVOT ROW BY PIVOT ELEMENT
C
      A(LL)=1.0D+00
      DO 350 L=1,N
      NCOL=(L-1)*NMAX+ICOLUM
  350 A(NCOL)=A(NCOL)/PIVOT
      IF(M) 371,371,360
  360 DO 370 L=1,M
      NCOL=(L-1)*NMAX+ICOLUM
  370 B(NCOL)=B(NCOL)/PIVOT
C
C     REDUCE NON-PIVOT ROWS
C
  371 JCON=(ICOLUM-1)*NMAX
      DO 550 L1=1,N
      IF (L1-ICOLUM) 391,550,391
  391 LL=JCON+L1
      T=A(LL)
      A(LL)=0.0D+00
      DO 450 L=1,N
      NCON=(L-1)*NMAX
      NL1=NCON+L1
      NCOL=NCON+ICOLUM
  450 A(NL1)=A(NL1)-A(NCOL)*T
      IF (M) 550,550,460
  460 DO 500 L=1,M
      NCON=(L-1)*NMAX
      NL1=NCON+L1
      NCOL=NCON+ICOLUM
  500 B(NL1)=B(NL1)-B(NCOL)*T
  550 CONTINUE
C
C     INTERCHANGE COLUMNS
C
      DO 710 I=1,N
      L=N+1-I
      IF (INDEX1(L)-INDEX2(L)) 630,710,630
  630 JROW=INDEX1(L)
      JCOLUM=INDEX2(L)
      JCON=(JROW-1)*NMAX
      NCON=(JCOLUM-1)*NMAX
      DO 705 K=1,N
      NROW=JCON+K
      NCOL=NCON+K
      SWAP=A(NROW)
      A(NROW)=A(NCOL)
      A(NCOL)=SWAP
  705 CONTINUE
  710 CONTINUE
  740 RETURN
      END
C
C*MODULE CHGPEN  *DECK CALC_CHISQ
      SUBROUTINE CALC_CHISQ(X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XVAL,NLAYER,
     *                  EFPOS,EFC,EFD,EFQ,EFO)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (NMAX=50,MXEFC=200)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
      COMMON /PDCPAR/ CENTER(3),DPOLE(3),QPOLE(6),RMAX,DELR,CONSTR,
     *                PTSEL,VDWSCL,PTDENS,VDWINC,NFREQ,LAYER,NPDC
C
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),
     *          QUANT(KOUNT),FITTED(KOUNT),DIFF(KOUNT),
     *          XVAL(NV),EFPOS(3,NEFC),
     *          EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),EFO(10,NEFC),
     *          NLAYER(LAYER+1)
      DIMENSION INDEX1(NV),INDEX2(NV),IPIVOT(NV),XQ(NV,NV),
     *          BQ(NV),VQ(NV)
      DATA GEODES /8HGEODESIC/
C
C      APPLIES CURRENT EXPONENTS
C
      DATA ZERO,ONE,TWO/0.0D+00,1.0D+00,2.0D+00/
C
C---- STORE CURRENT DAMPING EXPONENTS
C LOOP OVER MULTIPOLAR CENTERS
C LOOP OVER TERMS IN DAMPING FUNCTION
C IF THIS TERM HAS ACTIVE EXPONENT, GET CURRENT VALUE FROM -XVAL-
C ISETA MAPS ACTIVE EXPONENTS ONTO PROPER DAMPING CENTER
C
      IF(ISCCHG.EQ.1) THEN
         NEFCC=NMAIN
      ELSE 
         NEFCC=NEFC
      END IF   
C 
      DO 10 IEFC=1,NEFCC
            IND=ISETA(IEFC)
            IF(IND.EQ.0) GO TO 10
            ALP(IEFC)=XVAL(IND)
C            WRITE(6,5) IEFC,ALP(J,IEFC)
C    5       FORMAT(1X,'CALC_CHISQ: IEFC,ALP = ',I6,F12.5)
   10 CONTINUE
C
C ----- SET COEFFICIENTS TO 1.0 IF ICFIX IS TURNED ON
C
      IF (ICFIX.EQ.1) THEN
         DO IEFC=1,NEFCC
               COEF(IEFC)=ONE
         ENDDO
      ELSE
C
C       IF ICFIX IS NOT ON, DO MATRIX INVERSION LEAST SQUARES
C       ZERO OUT ARRAYS USED IN LINEAR COEFFICIENT DETERMINATION -----
C
      DO 100 I=1,NV
         VQ(I)=ZERO
         BQ(I)=ZERO
         DO 100 J=1,NV
            XQ(J,I)=ZERO
100   CONTINUE
C
C       ACCUMULATE SUMS FOR LINEAR LEAST SQUARES -----
C
      DO 500 K=1,KOUNT
         IF(IERRF.EQ.0) EVAL=DIFF(K)
         IF(IERRF.EQ.1) EVAL=QUANT(K)
         IF(IERRF.EQ.2) EVAL=DIFF(K)
         IF(IERRF.EQ.3) EVAL=DIFF(K)
         DO 101 I=1,NV
            VQ(I)=ZERO
101      CONTINUE
         DO 450 IEFC=1,NEFC
            XX=(X(K)-EFPOS(1,IEFC))**2
            YY=(Y(K)-EFPOS(2,IEFC))**2
            ZZ=(Z(K)-EFPOS(3,IEFC))**2
            RRI=XX+YY+ZZ
C
            VALI=CGPMUL(K,IEFC,X,Y,Z,KOUNT,EFPOS,EFC,EFD,EFQ,EFO,NEFC)
C
            IND=ISETA(IEFC)
C
C    WE ONLY GET COEFFICIENTS FOR ACTIVE (NON-FROZEN) TERMS
            IF(IND.EQ.0)GO TO 450
C
C    THE FOLLOWING CHANGE IS FOR ERROR FUNCTIONS
            ALPRRI=ALP(IEFC)*RRI
            IF(IERRF.EQ.0) DAMPI=EXP(-ALPRRI)
            IF(IERRF.EQ.1) DAMPI=-ERRF(SQRT(ALPRRI))
            IF(IERRF.EQ.2) DAMPI=EXP(-ALP(IEFC)*SQRT(RRI))
            IF(IERRF.EQ.3) DAMPI=EXP(-TWO*ALP(IEFC)*SQRT(RRI))*
     *                     (ONE+ALP(IEFC)*SQRT(RRI))
            DVALI=VALI*DAMPI
C
            VQ(IND)=VQ(IND)+DVALI
            BQ(IND)=BQ(IND)-DVALI*EVAL
450      CONTINUE
C
         DO 350 I=1,NV
            DO 350 J=1,NV
                XQ(I,J)=XQ(I,J)+VQ(I)*VQ(J)
350      CONTINUE
500   CONTINUE
C
C
C       DO MATRIX INVERSION AND STORE COEFFICIENTS -----
C
      CALL CGPINV(XQ,NV,BQ,NV,1,DET,INDEX1,INDEX2,IPIVOT,PIVOT)
C
      DO 600 IEFC=1,NEFC
         IND=ISETA(IEFC)
         IF(IND.EQ.0)GO TO 600
         COEF(IEFC)=BQ(IND)
C         WRITE(6,510) IEFC,COEF(IEFC)
C  510    FORMAT('AFTER CGPINV, IEFC,COEF = ',I6,F12.5)
600   CONTINUE
C
      ENDIF
C
C ----- CALCULATE CHISQ -----
C
C      CGPFIT ACTUALLY CALCULATES THE ELECTROSTATIC POTENTIAL DUE TO
C      THE DAMPED MULTIPOLES
C
      CALL CGPFIT(X,Y,Z,FITTED,KOUNT,
     *            EFPOS,EFC,EFD,EFQ,EFO)
C
C      WRITE(6,*)'IN CALC_CHISQ: NLAYER',NLAYER
      CHISQ=ZERO
      NEW_KOUNT = 0
      IF(PTSEL.EQ.GEODES) THEN
         DO ILAYER=1,LAYER
            SCALED = (VDWSCL+VDWINC*(ILAYER-1))/VDWSCL
            CHISQ_L = ZERO
            DO N=NLAYER(ILAYER)+1,NLAYER(ILAYER+1)
               XDIFF=(QUANT(N)-FITTED(N))*627.52D+00
               CHISQ_L=CHISQ_L+XDIFF**2
             END DO
             KOUNT_L = INT(SCALED*SCALED*
     *                    (NLAYER(ILAYER+1)-NLAYER(ILAYER)))
             NEW_KOUNT = NEW_KOUNT + KOUNT_L
             CHISQ_L = CHISQ_L * SCALED*SCALED
             CHISQ = CHISQ + CHISQ_L
          END DO
C
      ELSE
      DO 700 K=1,KOUNT
C         XDIFF=(QUANT(K)-FITTED(K))
         XDIFF=(QUANT(K)-FITTED(K))*627.52D+00
         CHISQ=CHISQ+XDIFF**2
700   CONTINUE
      NEW_KOUNT = KOUNT
      END IF
C
      CHISQ=SQRT(CHISQ/(NEW_KOUNT-1))
C
C       THE NEXT IS A KEY DEBUGGING CLAUSE, TO CHECK EVERY STEP WE TAKE.
C
C     WRITE(6,900) CHISQ
C     WRITE(6,901) (ALP(J),J=1,NEFC)
C 900 FORMAT(1X,'CHISQ=',F16.8)
C 901 FORMAT(1X,9F8.3)
C
      RETURN
      END
C
C*MODULE CHGPEN  *DECK CGPPRT
      SUBROUTINE CGPPRT(X,Y,Z,QUANT,FITTED,RAW,
     *                  EFNAM,EFPOS,KOUNT,NEW_KOUNT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),QUANT(KOUNT),
     *          FITTED(KOUNT),RAW(KOUNT),EFNAM(NEFC),EFPOS(3,NEFC)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      PARAMETER (NMAX=50,MXEFC=200)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA ZERO/0.0D+00/
C
C         PRINT OUT FINAL RESULTS OF CHARGE SCREENING
C
      IF(.NOT.MASWRK) RETURN
C
      WRITE(IW,10)
      DO 30 K=1,NEFC
         WRITE(IW,20) EFNAM(K),(EFPOS(I,K),I=1,3),
     *                COEF(K),ALP(K)
         WRITE(IP,21) EFNAM(K),COEF(K),ALP(K)
30    CONTINUE
C
      SUM=ZERO
      RSUM=ZERO
      RSQ=ZERO
      RSQ1=ZERO
      CHISQ=ZERO
      CHISQ1=ZERO
      NUM1=0
      NUMR1=0
      RSQ2=ZERO
      CHISQ2=ZERO
      NUM2=0
      NUMR2=0
      RSQ3=0
      CHISQ3=ZERO
      NUM3=0
      NUMR3=0
      RSQ4=ZERO
      CHISQ4=ZERO
      NUMR4=0
      NUM4=0
      DO 100 K=1,KOUNT
      XDIFF=(QUANT(K)-FITTED(K))*627.52D+00
      YDIFF=(QUANT(K)-RAW(K))*627.52D+00
      SUM=SUM+ABS(XDIFF)
      RSUM=RSUM+ABS(YDIFF)
      CHISQ=CHISQ+XDIFF**2
      RSQ=RSQ+YDIFF**2
      IF(ABS(XDIFF).GT.THRSH/8.0D+00) THEN
         CHISQ1=CHISQ1+XDIFF**2
         NUM1=NUM1+1
      END IF
      IF(ABS(XDIFF).GT.THRSH/4.0D+00) THEN
         CHISQ2=CHISQ2+XDIFF**2
         NUM2=NUM2+1
      END IF
      IF(ABS(XDIFF).GT.THRSH/2.0D+00) THEN
         CHISQ3=CHISQ3+XDIFF**2
         NUM3=NUM3+1
      END IF
      IF(ABS(XDIFF).GT.THRSH) THEN
         CHISQ4=CHISQ4+XDIFF**2
         NUM4=NUM4+1
      END IF
      IF(ABS(YDIFF).GT.THRSH/8.0D+00) THEN
         RSQ1=RSQ1+YDIFF**2
         NUMR1=NUMR1+1
      END IF
      IF(ABS(YDIFF).GT.THRSH/4.0D+00) THEN
         RSQ2=RSQ2+YDIFF**2
         NUMR2=NUMR2+1
      END IF
      IF(ABS(YDIFF).GT.THRSH/2.0D+00) THEN
         RSQ3=RSQ3+YDIFF**2
         NUMR3=NUMR3+1
      END IF
      IF(ABS(YDIFF).GT.THRSH) THEN
         RSQ4=RSQ4+YDIFF**2
         NUMR4=NUMR4+1
      END IF
100   CONTINUE
C
      WRITE(IW,105)
      RMS1 = ZERO
      RMS2 = ZERO
      RMS3 = ZERO
      RMS4 = ZERO
                     RMSNEW =SQRT(RSQ /(NEW_KOUNT-1))
                     RMS    =SQRT(RSQ /(KOUNT-1))
      IF(NUMR1.GT.1) RMS1   =SQRT(RSQ1/(NUMR1-1))
      IF(NUMR2.GT.1) RMS2   =SQRT(RSQ2/(NUMR2-1))
      IF(NUMR3.GT.1) RMS3   =SQRT(RSQ3/(NUMR3-1))
      IF(NUMR4.GT.1) RMS4   =SQRT(RSQ4/(NUMR4-1))
      AVG=RSUM/KOUNT
      AVGNEW=RSUM/NEW_KOUNT
      WRITE(IW,110) KOUNT,AVG,KOUNT,RMS,NUMR1,RMS1,NUMR2,RMS2,NUMR3,
     *              RMS3,NUMR4,RMS4,NEW_KOUNT,AVGNEW,NEW_KOUNT,RMSNEW
C
      WRITE(IW,106)
      RMS1 = ZERO
      RMS2 = ZERO
      RMS3 = ZERO
      RMS4 = ZERO
                    RMS    =SQRT(CHISQ /(KOUNT-1))
                    RMSNEW =SQRT(CHISQ /(NEW_KOUNT-1))
      IF(NUM1.GT.1) RMS1   =SQRT(CHISQ1/(NUM1 -1))
      IF(NUM2.GT.1) RMS2   =SQRT(CHISQ2/(NUM2 -1))
      IF(NUM3.GT.1) RMS3   =SQRT(CHISQ3/(NUM3 -1))
      IF(NUM4.GT.1) RMS4   =SQRT(CHISQ4/(NUM4 -1))
      AVG=SUM/KOUNT
      AVGNEW=SUM/NEW_KOUNT
      WRITE(IW,110) KOUNT,AVG,KOUNT,RMS,NUM1,RMS1,NUM2,RMS2,NUM3,
     *              RMS3,NUM4,RMS4,NEW_KOUNT,AVGNEW,NEW_KOUNT,RMSNEW
C
      WRITE(IW,120) THRSH
      DO 130 K=1,KOUNT
         SDIFF=(QUANT(K)-FITTED(K))*627.52D+00
         RDIFF=(QUANT(K)-RAW(K))*627.52D+00
C
         NN = -1
         TDIST = 1.0D+23
         DO 35 N=1,NEFC
             ADIST = SQRT((X(K)-EFPOS(1,N))**2+(Y(K)-EFPOS(2,N))**2+
     *              (Z(K)-EFPOS(3,N))**2)
             IF(ADIST .LT. TDIST) THEN
                TDIST = ADIST
                NN = N
             END IF
 35      CONTINUE       
C
         IF((ABS(SDIFF).LT.THRSH) .AND. 
     *       (ABS(RDIFF).LT.THRSH*2.0D+00)) GO TO 130
         WRITE(IW,125) K,EFNAM(NN),TDIST,X(K),Y(K),Z(K),
     *                 QUANT(K)*627.52D+00,
     *                 SDIFF,RDIFF
  130 CONTINUE
C      
      RETURN
C
   10 FORMAT(/15X,'----- OPTIMIZED SCREENING PARAMETERS ----- '/
     *       1X,'NAME',10X,'X',8X,'Y',8X,'Z',15X,'COEF',17X,'ALPHA')
   20 FORMAT(1X,A8,3(1X,F8.4),2(1X,F20.10),:,(/36X,2F20.10))
  105 FORMAT(/1X,'STATISTICS FOR UNSCREENED MULTIPOLES'/
     1        1X,'------------------------------------')
  106 FORMAT(/1X,'STATISTICS FOR SCREENED MULTIPOLES'/
     1        1X,'----------------------------------')
  110 FORMAT(' AVERAGE UNSIGNED ERROR        (',I10,' PTS) = ',1P,E12.6/
     * ' RMS DEVIATION                 (',I10,' PTS) = ',1P,E12.6/
     * ' RMS DEVIATION (ERR > THRSH/8) (',I10,' PTS) = ',1P,E12.6/
     * ' RMS DEVIATION (ERR > THRSH/4) (',I10,' PTS) = ',1P,E12.6/
     * ' RMS DEVIATION (ERR > THRSH/2) (',I10,' PTS) = ',1P,E12.6/
     * ' RMS DEVIATION (ERR > THRSH)   (',I10,' PTS) = ',1P,E12.6/
     *       ' '/
     *  ' EFFECTIVE AV. UNSIGNED ERROR  (',I10,' PTS) = ',1P,E12.6/
     *  ' EFFECTIVE RMS DEVIATION       (',I10,' PTS) = ',1P,E12.6/)
  120 FORMAT(/1X,'GRID POINTS WITH ERRORS (QUANTUM - EITHER FIT) ',
     *        'GREATER THAN',F8.2,' KCAL/MOL'/
     *        1X,'  POINT  CENTER  DISTANCE   X     Y       Z',
     *           '    QUANTUM      DIFF      DIFF   '/
     *        1X,'                 FROM ORIGIN              ',
     *           '                SCREEN  NOSCREEN'/)
  125 FORMAT(1X,I8,2X,A6,F7.3,3F7.3,3F10.2)
C
   21 FORMAT(1X,A8,2F14.9)
      END
C
C*MODULE CHGPEN  *DECK POWELL
      SUBROUTINE POWELL(P,XI,N,FTOL,ITER,FRET,PT,PTT,XIT,
     *                  X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *                  EFPOS,EFC,EFD,EFQ,EFO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C      EXTERNAL CALC_CHISQ
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION P(N),XI(N,N),IFROZN(N)
      DIMENSION PT(N),PTT(N),XIT(N)
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),QUANT(KOUNT),FITTED(KOUNT),
     *          DIFF(KOUNT), 
     *          XT(N),EFPOS(3,NEFC),EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),
     *          EFO(10,NEFC),NLAYER(1000)
C
      PARAMETER (NMAX=50,MXEFC=200)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     MINIMIZATION OF A FUNCTION OF N VARIABLES. INPUT CONSISTS 
C     OF AN INITIAL STARTING POINT P(1:N), AN INITIAL MATRIX XI(1:N,1:N)
C     WITH PHYSICAL DIMENTIONS NP BY NP, AND WHOSE COLUMNS CONTAIN
C     THE INITIAL SET OF DIRECTIONS (USUALLY THE N UNIT VECTORS); AND 
C     FTOL, THE FRACTIONAL TOLERANCE IN THE FUNCTION VALUE SUCH THAT 
C     FAILURE TO DECREASE BY MORE THAN THIS AMOUNT ON ONE ITERATION 
C     SIGNALS DONENESS. 
C     ON OUTPUT, P IS SET TO THE BEST POINT, XI IS THE THEN-CURRENT 
C     DIRECTION SET, FRET IS THE RETURNED FUNCTION VALUE AT P, AND 
C     ITER IS THE NUMBER OF ITERATIONS TAKEN. 
C     THE ROUTINE LINMIM IS USED.
C     ADOPTED FROM 'NUMERICAL RECIPIES'
C     LVS 2006
C
C
      DATA ZERO,TWO/0.0D+00,2.0D+00/
C
C
C      WRITE(IW,*)'START POWELL'
C       WRITE(IW,*)'NEFC=',NEFC
C       WRITE(IW,*) 'ALPHAS',(ALP(I),I=1,NEFC)
C       WRITE(IW,*) 'COEF',(COEF(I),I=1,NEFC)
C       WRITE(IW,*) 'ISETA',(ISETA(I),I=1,NEFC)
C       WRITE(IW,*) 'INITIAL SET OF DIRECTION',(XI(I,I),I=1,N)
C
C       CALCULATE INITIAL CHISQ
C
      CALL CALC_CHISQ(X,Y,Z,QUANT,FITTED,DIFF,KOUNT,P,NLAYER,
     *                EFPOS,EFC,EFD,EFQ,EFO)
      FRET = CHISQ
C      FRET = FUNC(P)
      WRITE(IW,100) FRET
C
C       SAVE THE INITIAL POINT
      DO 11 J=1,N
         PT(J)=P(J)
         IFROZN(J) = 0
  11  ENDDO 
      ITER = 0
  1   ITER = ITER + 1
      FP = FRET
      IBIG = 0
C       WILL BE THE BIGGEST FUNCTION DECREASE
      DEL = ZERO
C      
C         IN EACH ITERATION, LOOP OVER ALL DIRECTIONS IN THE SET
C         COPY THE DIRECTION, MINIMIZE ALONG IT, 
C         AND RECORD IT IF IT IS THE LARGEST DECREASE SO FAR
C
      DO 13 I=1,N
         IF(IFROZN(I).EQ.0) THEN
         DO 12 J=1,N
            XIT(J)= XI(J,I)
 12      ENDDO 
         FPTT = FRET
         CALL CGPLINMIN(IFFAIL,P,XIT,N,FRET,X,Y,Z,QUANT,FITTED,DIFF,
     *                  KOUNT,XT,NLAYER,EFPOS,EFC,EFD,EFQ,EFO)
C         IF(IFFAIL.EQ.0) THEN
            IF(ABS(FPTT-FRET).GT.DEL) THEN 
               DEL = ABS(FPTT-FRET)
               IBIG = I
            END IF
         END IF
 13   ENDDO
C
C         TERMINATION CRITERION
      IF(TWO*ABS(FP-FRET).LE.FTOL*(ABS(FP)+ABS(FRET))) THEN
         IF(MASWRK) WRITE(IW,9000) NF
C             NORNAL CASE
         IF(IFFAIL.EQ.0) THEN 
            IF(MASWRK) WRITE(IW,9001) ITER,FRET,(P(I),I=1,N)
C             ALMOST CONVERGED WHEN PARAMETER GOES TO +INFINITY
         ELSE         
            ITER = ITER -1
C            DO I=1,N
C               P(I)=PT(I)
C            END DO
C            CALL CALC_CHISQ(X,Y,Z,QUANT,FITTED,DIFF,KOUNT,P,NLAYER,
C     *                EFPOS,EFC,EFD,EFQ,EFO)
C            FRET = CHISQ
            IF(MASWRK) WRITE(IW,9002) ITER,FRET,(P(I),I=1,N)    
         END IF
         RETURN
      END IF
      IF(ITER.EQ.MAXITS) THEN
         WRITE(IW,*)'POWELL EXCEEDING MAXIMUM ITERATIONS'
         RETURN
      END IF
      IF(MASWRK) WRITE(IW,9010) ITER,FRET
C
C        CONSTRUCT THE EXTRAPOLATED POINT AND THE AVERAGE 
C        DIRECTION MOVED. SAVE THE OLD STARTING POINT.    
C
      DO J=1,N
         IF(IFROZN(J).EQ.0) THEN
            PTT(J)=TWO*P(J)-PT(J)
            XIT(J)=    P(J)-PT(J)
             PT(J)=    P(J)
         END IF
      ENDDO      
C
C        REJECT NEGATIVE EXPONENTS OUT OF HAND
C
      DO J=1,N
        IF(PTT(J).LE.ZERO) GO TO 1
      ENDDO
C
C       FUNCTION VALUE AT EXTRAPOLATED POINT
C
      CALL CALC_CHISQ(X,Y,Z,QUANT,FITTED,DIFF,KOUNT,PTT,NLAYER,
     *                EFPOS,EFC,EFD,EFQ,EFO)
      FPTT = CHISQ
C
C        ONE REASON NOT TO USE NEW DIRECTION
      IF(FPTT.GE.FP) GOTO 1
C        OTHER REASON NOT TO USE NEW DIRECTION
      T=TWO*(FP-TWO*FRET+FPTT)*(FP-FRET-DEL)**2-DEL*(FP-FPTT)**2
      IF(T.GE.ZERO) GOTO 1
C
C        MOVE TO THE MINIMUM OF THE NEW DIRECTION,
C        AND SAVE THE NEW DIRECTION
C---  WRITE(IW,*)'  POWELL: NEW DIRECTION IS ACCEPTED'
C
      CALL CGPLINMIN(IFFAIL,P,XIT,N,FRET,X,Y,Z,QUANT,FITTED,DIFF,
     *               KOUNT,XT,NLAYER,EFPOS,EFC,EFD,EFQ,EFO)
      DO 15 J=1,N
C         XI(J,IBIG)=XI(J,N)
C         XI(J,N)=XIT(J)
         XI(J,IBIG)=XIT(J)
 15   ENDDO 
C
C        BACK TO ANOTHER ITERATION
      GOTO 1
C
 100  FORMAT(1X,'INITIAL CHI^2 = ',F16.8)
 9000 FORMAT(1X,'FUNCTION WAS EVALUATED ',I6,' TIMES')
 9001 FORMAT(1X,'POWELL MINIMIZATION PROCEDURE CONVERGED',
     *      ' IN',I6,' ITERATIONS',/
     *      1X,'FINAL VALUE OF CHI^2 IS',F16.8/
     *      1X,'OPTIMIZED NONLINEAR PARAMETERS ARE',/
     *      1X,8(1X,F8.4))
 9002 FORMAT(1X,'POWELL MINIMIZATION PROCEDURE ALMOST CONVERGED',
     *      ' IN',I6,' ITERATIONS',/
     *      1X,'FINAL VALUE OF CHI^2 IS',F16.8/
     *      1X,'OPTIMIZED NONLINEAR PARAMETERS ARE',/
     *      1X,8(1X,F8.4))
 9010 FORMAT(1X,'POWELL PROCEDURE, ITERATION ',I6,' CHI^2 = ',F16.8)
      END
C
C
C*MODULE CHGPEN  *DECK CGPLINMIN
      SUBROUTINE CGPLINMIN(IFFAIL,P,XI,N,FRET,X,Y,Z,QUANT,FITTED,DIFF,
     *                     KOUNT,XT,NLAYER,EFPOS,EFC,EFD,EFQ,EFO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C          MAXIMUM ANTICIPATED N, AND TOL PASSED TO BRENT
      PARAMETER(TOL=1.0D-04)
      PARAMETER (NMAX=50,MXEFC=200)
C
      DIMENSION P(N),XI(N)
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),QUANT(KOUNT),FITTED(KOUNT),
     *          DIFF(KOUNT),
     *          XT(NV),EFPOS(3,NEFC),EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),
     *          EFO(10,NEFC),NLAYER(1000)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
C
C     GIVEN AN N-DIMENSIONAL POINT P(1:N) AND AN N-DIMENSIONAL DIRECTION 
C     XI(1:N), MOVES AND RESETS P TO WHERE THE FUNCTION FUNC(P) TAKES ON 
C     A MINIMUM ALONG THE DIRECTION XI FROM P, AND REPLACES XI BY THE 
C     ACTUAL VECTOR DISPLACEMENT THAT P WAS MOVED. ALSO RETURNS AS FRET
C     THE VALUE OF FUNC AT THE RETURNED LOCATION P. THIS IS ACTUALLY ALL
C     ACCOMPLISHED BY CALLING THE ROUTINES MNBRAK AND BRENT.
C     ADOPTED FROM 'NUMERICAL RECIPIES'
C     LVS 2006
C     
      DATA ZERO/0.0D+00/
C
C       SET UP COMMON BLOCK
      N=NV
      DO 11 J =1,N
         PCOM(J)=P(J)
         XICOM(J)=XI(J)
  11  ENDDO 
C
      IFFAIL = 0
C        INITIAL GUESS FOR BRACKETS
      AX = 0.0D+00
      XX = 0.5D+00
      CALL MNBRAK(IFFAIL,AX,XX,BX,FA,FX,FB,X,Y,Z,QUANT,FITTED,
     *           DIFF, KOUNT,XT,NLAYER,EFPOS,EFC,EFD,EFQ,EFO)
C
C      WRITE(IW,*)'  AX=',AX,'XX=',XX,'BX=',BX
      IF(IFFAIL.EQ.0) THEN
         XMIN = ZERO
         CALL BRENT(FRET,AX,XX,BX,TOL,XMIN,X,Y,Z,QUANT,FITTED,DIFF,
     *              KOUNT,XT,NLAYER,EFPOS,EFC,EFD,EFQ,EFO)
C
C        CONSTRUCT THE VECTOR RESULTS TO RETURN
      DO 12 J =1,N
         XI(J)=XMIN*XI(J)
         P(J)=P(J)+XI(J)
 12   ENDDO 
      ELSE
         XMIN = XX
         FRET = FX
C        CONSTRUCT THE VECTOR RESULTS TO RETURN
      DO 13 J =1,N
         XI(J)=XICOM(J)
         P(J)=PCOM(J)
 13   ENDDO 
      END IF
C
C      WRITE(IW,*)'XI=',XI,'P=',P
C
      RETURN
      END
C
C*MODULE CHGPEN  *DECK F1DIM
      SUBROUTINE F1DIM(RETVAL,XX,X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,
     *                 NLAYER,EFPOS,EFC,EFD,EFQ,EFO,IFFAIL,ALPMAX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (NMAX=50, MXEFC=200)
C
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),QUANT(KOUNT),FITTED(KOUNT),
     *          DIFF(KOUNT),
     *          XT(NV),EFPOS(3,NEFC),EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),
     *          EFO(10,NEFC),NLAYER(1000)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
C
C     USED BY CGPLINMIN AS THE FUNCTION PASSED TO MNBRAK AND BRENT
C     
      DATA ONE /1.0D+00/
C
C          ALPMIN DEPENDS ON TYPE OF FIT: EXP. OR GAUSSIAN
C
      IF(IERRF .EQ. 2) THEN
         ALPMIN = 0.5D+00
         GAUSS_SCALE = ONE
      ELSE 
         ALPMIN = 0.5D+00*0.5D+00*2.709498091D-01
         GAUSS_SCALE = 2.709498091D-01
      END IF
C      
      IFFAIL = 0
      DO 11 J =1,NV
         XT(J)=PCOM(J)+XX*XICOM(J)
         IF (XT(J).GT.ALPMAX) THEN 
            IFFAIL = 1
C            WRITE(IW,*)'   F1DIM: ALP(',J,')=',XT(J)
            XT(J) = ALPMAX - ONE*GAUSS_SCALE
C            JBAD =J
            XICOM(J)=-0.1D+00*XICOM(J)
C            XICOM(J)=0.0D+00
            PCOM(J) = XT(J)
         END IF
         IF (XT(J).LT.ALPMIN) THEN
            IFFAIL =1
C            WRITE(IW,*)'   F1DIM: ALP(',J,')=',XT(J)
            XT(J) = ALPMIN+ONE*GAUSS_SCALE
            XICOM(J)=-0.1D+00*XICOM(J)
            PCOM(J) = XT(J)
         END IF
 11   ENDDO
C
      CALL CALC_CHISQ(X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *                EFPOS,EFC,EFD,EFQ,EFO)
      RETVAL = CHISQ
      NF=NF+1
C
      RETURN
      END
C
C
C*MODULE CHGPEN  *DECK MNBRAK
      SUBROUTINE MNBRAK(IFFAIL,AX,BX,CX,FA,FB,FC,
     *                  X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *                  EFPOS,EFC,EFD,EFQ,EFO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER(GOLD=1.618034D+00)
      PARAMETER (NMAX=50,MXEFC=200)
C
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),QUANT(KOUNT),FITTED(KOUNT),
     *          DIFF(KOUNT),
     *          XT(NV),EFPOS(3,NEFC),EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),
     *          EFO(10,NEFC),NLAYER(1000)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     GIVEN A FUNCTION FUNC, AND GIVEN DISTINCT INITIAL POINTS AX 
C     AND BX, THIS ROUTINE SEARCHES IN THE DOWNHILL DIRECTION (DEFINED
C     BY THE FUNCTION AS EVALUATED AT THE INITIAL POINTS) AND RETURNS
C     NEW POINTS AX, BX, CX THAT BRACKET A MINIMUM OF THE FUNCTION.
C     ALSO RETURNED ARE THE FUNCTION VALUES AT THE THREE POINTS, FA,FB, 
C     AND FC.
C     PARAMETERS: GOLD IS THE DEFAULT RATIO BY WHICH SUCCESSIVE
C     INTERVALS ARE MAGNIFIED; GLIMIT IS THE MAXIMUM MAGNIFICATION
C     ALLOWED FOR A PARABOLIC-FIT STEP.
C     ADOPTED FROM 'NUMERICAL RECIPIES'
C     LVS 2006
C     
C        VALUE OF MAXIMUM ALLOWED VALUE OF ALPHA'S (IN CALC_CHISQ)
C        IT DEPENDS ON THE RADIUS OF FIRST SPHERE (VDWSCL)
C--    XMAX = 6.0D+00/VDWSCL
C        IT MIGHT BE MADE TO DEPEND ON TYPE OF FIT (EXP. OR GAUSSIAN),
C        BUT FOR NOW IT IS KEPT THE SAME FOR BOTH
C--   IF(IERRF .EQ. 2) THEN
C--      XMAX = 11.0D+00
C--   ELSE 
C--      XMAX=11.0D+00*11.0D+00*2.709498091D-01
C--   END IF
C                 NOTE, THIS BLOCK REALLY SETS XMAX TO 10.0 (SEE F1DIM)
      IF(IERRF.EQ.2) THEN
         GAUSS_SCALE = 1.0D+00
      ELSE 
         GAUSS_SCALE = 2.709498091D-01
      END IF
      XMAX = 10.0D+00 + GAUSS_SCALE
C
C        IF MNBRAK FAILED (IFFAIL=1) TO BRACKET A MINIMUM
      IFFAIL = 0
C
C---  WRITE(IW,*) 'IN MNBRAK...'
C
      CALL F1DIM(FA,AX,X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *          EFPOS,EFC,EFD,EFQ,EFO,IFFAIL,XMAX)
      CALL F1DIM(FB,BX,X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *          EFPOS,EFC,EFD,EFQ,EFO,IFFAIL,XMAX)
C
      IF(IFFAIL.EQ.1) IFFAIL=0
C
C        SWITCH ROLES OF A AND B SO THAT WE CAN GO DOWNHILL IN THE 
C        DIRECTION FROM A TO B.
      IF(FB.GT.FA) THEN
         DUM = AX
         AX = BX
         BX = DUM
         DUM = FB
         FB = FA
         FA = DUM
      ENDIF
C
C        FIRST GUESS FOR C
      CX = BX +GOLD*(BX-AX)
      CALL F1DIM(FC,CX,X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *          EFPOS,EFC,EFD,EFQ,EFO, IFFAIL,XMAX)
C
C        "DO WHILE": KEEP RETURNING HERE UNTIL WE BRACKET
      ITER = 0
  1   IF(FB.GE.FC) THEN
C         
C               REJECT PARABOLIC U, USE DEFAULT MAGNIFICATION
            U=CX+GOLD*(CX-BX)
            CALL F1DIM(FU,U,X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *                 EFPOS,EFC,EFD,EFQ,EFO,  IFFAIL,XMAX)
C
C             ELIMINATE OLDEST POINT AND CONTINUE
         AX=BX
         BX=CX
         CX=U
         FA=FB
         FB=FC
         FC=FU
         ITER = ITER + 1
C
C     THIS NEXT ONE HAPPENS FAIRLY OFTEN SO THE PRINTING IS TURNED OFF
C
            IF(IFFAIL.EQ.1) THEN
               BX=U
               FB=FU
C--               WRITE(IW,*) 'MNBRAK: FAILED, ALP VALUE TOO BIG'
               RETURN
            END IF 
C
         IF (ITER.GT.20) THEN
            WRITE(IW,*) 'MNBRAK: FAILED, MORE THAN 20 ITERATIONS'
            IFFAIL = 1
C--           BX = 6.0D+00
            RETURN
         END IF
         GOTO 1
      ENDIF
C
      IFFAIL=0
      RETURN
      END
C
C
C*MODULE CHGPEN  *DECK BRENT
      SUBROUTINE BRENT(RETVAL,AX,BX,CX,TOL,XMIN,X,Y,Z,QUANT,FITTED,
     *                 DIFF,KOUNT,XT,NLAYER,EFPOS,EFC,EFD,EFQ,EFO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER(ITMAX=100,CGOLD=0.3819660D+00,ZEPS=1.0D-05)
      PARAMETER (NMAX=50,MXEFC=200)
C
      DIMENSION X(KOUNT),Y(KOUNT),Z(KOUNT),QUANT(KOUNT),FITTED(KOUNT),
     *          DIFF(KOUNT),
     *          XT(NV),EFPOS(3,NEFC),EFC(NEFC),EFD(3,NEFC),EFQ(6,NEFC),
     *          EFO(10,NEFC),NLAYER(1000)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     GIVEN A FUNCTION F, AND GIVEN A BRACKETING TRIPLET OF ABSCISSAS 
C     AX,BX,CX, THIS ROUTINE ISOLATES THE MINIMUM TO A FRACTIONSL 
C     PRECISION OF ABOUT TOL USING BRENT'S METHOD. THE ABSCISSA OF THE 
C     MINIMUM IS RETURNED AS XMIN, AND THE MINIMUM FUNCTION VALUE IS 
C     RETURNED AS BRENT, THE RETURNED FUNCTION VALUE.
C     PARAMETERS: MAX. ALLOWED NUMBER OF ITERATIONS, GOLDEN RATIO, AND 
C     A SMALL NUMBER THAT PROTECTS AGAINST TRYING TO ACHIEVE FRACTIONAL 
C     ACCURACY FOR A MINIMUM THAT HAPPENS TO BE EXACTLY ZERO. 
C     ADOPTED FROM 'NUMERICAL RECIPIES'
C     LVS 2006
C     
      DATA ZERO,HALF,TWO/0.0D+00,0.5D+00,2.0D+00/
C
C---  WRITE(IW,*) 'IN BRENT...'
C
      IFFAIL =0
      XMAX = 100.0D+00
C        A AND B MUST BE IN ASCENDING ORDER, THOUGH THE INPUT 
C        ABSCISSAS NEED NOT BE.
      A = MIN(AX,CX)
      B = MAX(AX,CX)
C
C        INITIALIZATION...
      V=BX
      W=V
      XXX=V
      D=ZERO
C
C         THIS WILL BE THE DISTANCE MOVED ON THE STEP BEFORE LAST
      E=ZERO
      CALL F1DIM(FX,XXX,X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *           EFPOS,EFC,EFD,EFQ,EFO,IFFAIL,XMAX)
      FV=FX
      FW=FX
C
C         MAIN LOOP
      DO 11 ITER=1,ITMAX
         XM=HALF*(A+B)
         TOL1=TOL*ABS(XXX)+ZEPS
         TOL2=TWO*TOL1
C
C            TEST FOR DONE HERE
         IF(ABS(XXX-XM).LE.(TOL2-HALF*(B-A))) GOTO 3
C
C            CONSTRUCT A TRIAL PARABOLIC FIT
         IF(ABS(E).GT.TOL1) THEN
            R=(XXX-W)*(FX-FV)
            Q=(XXX-V)*(FX-FW)
            P=(XXX-V)*Q-(XXX-W)*R
            Q=TWO*(Q-R)
            IF(Q.GT.ZERO) P=-P
            Q=ABS(Q)
            ETEMP = E
            E=D
            IF(ABS(P).GE.ABS(HALF*Q*ETEMP).OR.P.LE.Q*(A-XXX).OR.
     *         P.GE.Q*(B-XXX)) GOTO 1
C
C             THE ABOVE CONDITIONS DETERMINE THEACCEPTABILITY OF THE 
C             PARABOLIC FIT. HERE IT IS O.K.:
C             TAKE THE PARABOLIC STEP.
            D=P/Q
            U=XXX+D
            IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-XXX)
C
C              SKIP OVER THE GOLDEN SECTION STEP
            GOTO 2
         ENDIF
C
C           GOLDEN SECTION STEP
   1      IF(XXX.GE.XM) THEN
            E=A-XXX
         ELSE
            E=B-XXX
         ENDIF
         D=CGOLD*E
C
C            ARRIVE HERE WITH D COMPUTED EITHER FROM PARABOLIC FIT,
C            OR ELSE FROM GOLDEN SECTION
  2      IF(ABS(D).GE.TOL1) THEN
            U=XXX+D
         ELSE
            U=XXX+SIGN(TOL1,D)
         ENDIF
C
C           FUNCTION EVALUATION AND HOUSEKEEPING
         CALL F1DIM(FU,U,X,Y,Z,QUANT,FITTED,DIFF,KOUNT,XT,NLAYER,
     *              EFPOS,EFC,EFD,EFQ,EFO, IFFAIL,XMAX)
         IF(FU.LE.FX) THEN
            IF(U.GE.XXX) THEN
               A=XXX
            ELSE 
               B=XXX
            ENDIF
            V=W
            FV=FW
            W=XXX
            FW=FX
            XXX=U
            FX=FU
         ELSE
            IF(U.LT.XXX) THEN
               A=U
            ELSE
               B=U
            ENDIF
            IF(FU.LE.FW .OR. W.EQ.XXX) THEN
               V=W
               FV=FW
               W=U
               FW=FU
            ELSE IF(FU.LE.FV .OR. V.EQ.XXX .OR. V.EQ.W) THEN
               V=U
               FV=FU
            ENDIF
         ENDIF
 11   ENDDO 
C      
      WRITE(IW,*) 'BRENT EXCEED MAXIMUM ITERATIONS'
      CALL ABRT
C
C        EXIT WITH BEST VALUES
  3   XMIN=XXX
      RETVAL=FX
      RETURN
      END
C      
C*MODULE CHGPEN  *DECK CALC_ALP
      SUBROUTINE CALC_ALP(ALP_INI,IEQ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION ALP_INI(NV),IEQ(NEFC)
C
      PARAMETER (NMAX=50,MXEFC=200)
C
      COMMON /CGPPAR/ CHISQ,THRSH,ALP(MXEFC),COEF(MXEFC),
     *                PCOM(NMAX),XICOM(NMAX),NEFC,NMAIN,
     *                ISCCHG,MAXITS,ICFIX,IERRF,NV,NF,NFIT,NRUN,
     *                IFTTYP(2),IFTFIX(2),ISETA(MXEFC)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C        CALCULATES INITIAL VALUES OF PARAMETERS GIVEN TO POWELL
C        BASED ON USER DEFINED ALPHAS
C
      NALPX=0
      IF(ISCCHG.EQ.1) THEN
         NEFCC=NMAIN
      ELSE 
         NEFCC=NEFC
      END IF   
      DO 38 IEFC=1,NEFCC
         IF(IEQ(IEFC).NE.IEFC) GO TO 38
               NALPX = NALPX+1
               IF(NFIT.EQ.1) THEN
                  ALP_INI(NALPX)=ALP(IEFC)
               ELSE
                  IF(IFTTYP(1).EQ.2) THEN
                     ALP_INI(NALPX)=ALP(IEFC)*
     *                              ALP(IEFC)*2.709498091D-01
                  ELSE
                     ALP_INI(NALPX)=SQRT(ALP(IEFC)/2.709498091D-01)
                  END IF
               END IF
   38 CONTINUE
C
      WRITE(IW,101) (ALP_INI(I),I=1,NALPX)
 101  FORMAT(1X,' INITIAL ALPHAS ARE ',15F8.3)
      RETURN
      END
C
