C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C 13 JUL 05 - MWS - SVPODN: FINAL LOOP SHOULD ONLY RUN TO 10
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 27 JUN 05 - DMC - KEYWORD IOPPRD TO ALLOW SMALL/LARGE MEMORY USE
C  1 JUN 05 - DMC - ADD REWINDS, SMALL FIX IN SVPODN
C 30 APR 05 - DMC - FINISHING TOUCHES FOR THE INITIAL IMPLEMENTATION
C 22 FEB 05 - MWS - FIX MEMORY LEAK
C 13 FEB 05 - DMC - READ INPUT FOR SVP, AND GENERATE CAVITY INFO
C
C     -----      SVP REACTION FIELD CODE MODULE      -----
C     ----- VERSION 01.6 (08/06/01) FOR USE IN HONDO -----
C     -----   CONVERTED FEB 05 FOR USE IN GAMESS     -----
C     -----       DEVELOPED AND MAINTAINED BY        -----
C     -----     D. M. CHIPMAN, CHIPMAN.1@ND.EDU      -----
C
C*MODULE SVPINP  *DECK SVPINP
      SUBROUTINE SVPINP
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- SVP REACTION FIELD INPUT -----
C
      PARAMETER   (MXCAV=5810)
C
      LOGICAL     SOME
      LOGICAL     OUT
      LOGICAL     MORE
      LOGICAL     DBUG
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SVPOPT/ ISVP,NSVP
      COMMON/SVPFLT/RHOISO,TOLCAV,RSTART,RSPACE,RADSPH,RADSLV,TOLCHG,
     1              DIELST,DIELOP,ZIONP,ZIONM,TKELVN,ROTCAV(3),
     2              TRNCAV(3),CSDIAG,CVGLIN
      COMMON/SVPTGR/METHOD,IQORP,NSTAGE,IOPPRD,IRENRM,ISHAPE,ITRCAV,
     1       NDRCAV,ICVICK,IFRQUP,IRDRF,IONSTR,INTCAV,NPTLEB,NPTTHE,
     2       NPTPHI,IROTGR,ITRNGR,ICRD,LINEQ,IPNRF,IPRRF,IPNDN,IPRDN,
     3       IEX
      COMMON/SVPPAR/CAVC(3,MXCAV),CAVW(MXCAV),CAVQ(MXCAV),
     1              CAVN(3,MXCAV),SPHW(MXCAV),CAVP(MXCAV),
     2              CENCAV(3),RNTCAV(3),NCAV,IQP
      COMMON/SVPPR0/CAVC0(3,MXCAV),CAVW0(MXCAV),CAVQ0(MXCAV),
     1              CAVN0(3,MXCAV),SPHW0(MXCAV),CAVP0(MXCAV),
     2              CENCV0(3),RNTCV0(3),NCAV0,IQP0
      COMMON/SVPPRP/SURF,VOL,RHOEXT,CHGBAR,CHGRUL,CHGDEV
      COMMON/SVPDBG/IDBG
      COMMON/SVPCAV/IFPNCH
C
C         NAMELIST SIMULATION
C
      PARAMETER (NNAM=46)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (TOANGS=0.52917724924D+00,TOBOHR=ONE/TOANGS)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA SVPKWD/8HSVP     /
      DATA QNAM/8HMETHOD  ,8HIQORP   ,8HIOPPRD  ,8HNSTAGE  ,8HIRENRM  ,
     *          8HISHAPE  ,8HRHOISO  ,8HTOLCAV  ,8HITRCAV  ,8HNDRCAV  ,
     *          8HICVICK  ,8HRSTART  ,8HRSPACE  ,8HRADSPH  ,8HRADSLV  ,
     *          8HTOLCHG  ,8HIFRQUP  ,8HIRDRF   ,8HDIELST  ,8HDIELOP  ,
     *          8HINTCAV  ,8HNPTLEB  ,8HNPTTHE  ,8HNPTPHI  ,8HIROTGR  ,
     *          8HROTTHE  ,8HROTPHI  ,8HROTCHI  ,8HITRNGR  ,8HTRANX   ,
     *          8HTRANY   ,8HTRANZ   ,8HIONSTR  ,8HZIONP   ,8HZIONM   ,
     *          8HTKELVN  ,8HICRD    ,8HLINEQ   ,8HCVGLIN  ,8HCSDIAG  ,
     *          8HIPNRF   ,8HIPRRF   ,8HIPNDN   ,8HIPRDN   ,8HIEX     ,
     *          8HIFPNCH  /
      DATA KQNAM/1,1,1,1,1,   1,3,3,1,1,  1,3,3,3,3,  3,1,1,3,3,
     *           1,1,1,1,1,   3,3,3,1,3,  3,3,1,3,3,  3,1,1,3,3,
     *           1,1,1,1,1,   1/
C
      DBUG=.FALSE.
      MORE=.FALSE.
      MORE=MORE.OR.DBUG
      OUT =.FALSE.
      OUT =OUT.OR.MORE
      SOME=.FALSE.
      SOME=SOME.OR.OUT
C
      IF(DBUG) THEN
         IDBG=4
      ELSEIF(MORE) THEN
         IDBG=3
      ELSEIF(OUT) THEN
         IDBG=2
      ELSEIF(SOME) THEN
         IDBG=1
      ELSE
         IDBG=0
      ENDIF
C
C     -- SET DEFAULT VALUES OF PARAMETERS IN GROUP $SVP --
C
C     -- METHOD = APPROXIMATION TO SVP SOLVATION METHOD
C     -- IQORP  = TYPE OF SURFACE DISTRIBUTION FOR REPRESENTATION
C                 OF RXN FIELD
C     -- IOPPRD = FORM OF DS OR SD* OPERATOR PRODUCT
C     -- NSTAGE = NUMBER OF STAGES IN SOLUTION OF SYSTEM EQUATION
C     -- IRENRM = OPTION FOR AD HOC RENORMALIZATION OF SURFACE CHARGES
C     -- ISHAPE = SHAPE OF CAVITY SURFACE
C     -- RHOISO = ELECTRONIC ISODENSITY CONTOUR USED TO
C                 SPECIFY CAVITY
C     -- TOLCAV = CONVERGENCE CRITERION ON RHOISO
C     -- ITRCAV = MAX NUMBER OF ITERATIONS ALLOWED TO FIND RHOISO
C     -- NDRCAV = NUMBER OF ANALYTIC DENSITY DERIVATIVES USED
C                 TO FIND RHOISO
C     -- ICVICK = CHECK VALIDITY OF SINGLE CENTER EXPANSION
C     -- RSTART = RADIUS OF INNERMOST SHELL IN ISODENSITY CAVITY CHECK
C     -- RSPACE = SPACING BETWEEN SHELLS IN ISODENSITY CAVITY CHECK
C     -- RADSPH = RADIUS OF SPHERE USED TO CHARACTERIZE CAVITY
C                 SURFACE
C     -- RADSLV = RADIUS OF SHPERE USED TO MODEL ROLLING SOLVENT
C                 MOLECULE
C     -- TOLCHG = CONVERGENCE CRITERION FOR SURFACE CHARGES
C                 OR DIPOLES
C     -- IFRQUP = FREQUENCY OF UPDATING CAVITY SURFACE
C                 DETERMINATION
C     -- IRDRF  = READ REACTION FIELD FROM EXTERNAL FILE
C     -- DIELST = STATIC DIELECTRIC CONSTANT
C     -- DIELOP = OPTICAL DIELECTRIC CONSTANT
C     -- INTCAV = METHOD FOR INTEGRATION OF CAVITY SURFACE
C     -- NPTLEB = NUMBER OF LEBEDEV-STYLE POINTS
C     -- NPTTHE = NUMBER OF THETA POINTS IN (THETA,PHI)
C     -- NPTPHI = NUMBER OF PHI POINTS IN (THETA,PHI)
C     -- IROTGR = SPECIFY ROTATION OF CAVITY GRID
C     -- ROTTHE =
C     -- ROTPHI =
C     -- ROTCHI = EULER ANGLES FOR ROTATION OF CAVITY GRID
C     -- ITRNGR = SPECIFY TRANSLATION OF CAVITY GRID
C     -- TRANX  =
C     -- TRANY  =
C     -- TRANZ  = COORDINATES FOR TRANSLATON OF CAVITY GRID
C     -- IONSTR = METHOD FOR TREATMENT OF IONIC STRENGTH
C     -- ZIONP  = CATION STRENGTH IN MOLES/LITER
C     -- ZIONM  = ANIONIC STRENGTH IN MOLES/LITER
C     -- TKELVN = TEMPERATURE IN DEGREES KELVIN
C     -- ICRD   = METHOD FOR TREATMENT OF
C                 CAVITATION/REPULSION/DISPERSION
C     -- LINEQ  = METHOD FOR SOLVING LINEAR EQUATIONS
C                 FOR CHARGES
C     -- CVGLIN = CONVERGENCE CRITERION FOR SOLUTION OF
C                 LINEAR EQUATIONS
C     -- CSDIAG = FACTOR TO MULTIPLY DIAGONAL ELEMENTS OF
C                 SURF POT MATRIX S
C     -- IPNRF  = PUNCH REACTION FIELD
C     -- IPRRF  = PRINT REACTION FIELD
C     -- IPNDN  = PUNCH DENSITY PROPERTIES
C     -- IPRDN  = PRINT DENSITY PROPERTIES
C     -- IFPNCH = PUNCH CAVITY DATA ONLY (POINT+NORMAL+AREA)
C     -- IEX    = DISPLAY OPTIONS NOT YET IMPLEMENTED OR RESERVED
C                 FOR EXPERTS
C
      METHOD  = 0
      IQORP   = 0
      IOPPRD  = 0
      NSTAGE  = 0
      IRENRM  = 0
      ISHAPE  = 0
      RHOISO  = 0.001D+00
      TOLCAV  = 1.0D-10
      ITRCAV  = 99
      NDRCAV  = 0
      ICVICK  = 2
      RSTART  = 0.2D+00*TOBOHR
      RSPACE  = 0.2D+00*TOBOHR
      RADSPH  = 2.64561733D+00
      RADSLV  = 2.83458999D+00
      TOLCHG  = 1.0D-07
      IFRQUP  = 0
      IRDRF   = 0
      DIELST  = 78.39D+00
      DIELOP  = 1.776D+00
      INTCAV  = 0
      NPTLEB  = 110
      NPTTHE  = 8
      NPTPHI  = 16
      ITRNGR  = 2
      IROTGR  = 2
      ROTTHE  = 0.0D+00
      ROTPHI  = 0.0D+00
      ROTCHI  = 0.0D+00
      TRANX   = 0.0D+00
      TRANY   = 0.0D+00
      TRANZ   = 0.0D+00
      IONSTR  = 0
      ZIONP   = 0.0D+00
      ZIONM   = 0.0D+00
      TKELVN  = 298.0D+00
      ICRD    = 0
      LINEQ   = 1
      CVGLIN  = 1.0D-07
      CSDIAG  = 1.104D+00
      IPNRF   = 0
      IPRRF   = 0
      IPNDN   = 0
      IPRDN   = 0
      IFPNCH  = 0
      IEX     = 0
C
C     ----- READ INPUT GROUP -SVP- -----
C
      ISVP = 0
C
      JRET=0
      CALL NAMEIO(IR,JRET,SVPKWD,NNAM,QNAM,KQNAM,
     *            METHOD,IQORP,IOPPRD,NSTAGE,IRENRM,ISHAPE,
     *            RHOISO,TOLCAV,ITRCAV,NDRCAV,ICVICK,RSTART,RSPACE,
     *            RADSPH,RADSLV,TOLCHG,IFRQUP,IRDRF,DIELST,DIELOP,
     *            INTCAV,NPTLEB,NPTTHE,NPTPHI,
     *            IROTGR,ROTTHE,ROTPHI,ROTCHI,
     *            ITRNGR,TRANX,TRANY,TRANZ,
     *            IONSTR,ZIONP,ZIONM,TKELVN,ICRD,
     *            LINEQ,CVGLIN,CSDIAG,IPNRF,IPRRF,IPNDN,IPRDN,IEX,
     *            IFPNCH,
     *        0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
C
C        DETECTION OF THE $SVP GROUP ENABLES THE SVP COMPUTATION
C
      IF(JRET.EQ.1) THEN
         RETURN
      ELSE
         ISVP = 1
      END IF
C
      IF(JRET.EQ.2) THEN
        IF (MASWRK) WRITE(IW,*) 'ERROR READING $SVP INPUT, CHECK ABOVE'
        CALL ABRT
      END IF
C
      IF(NPROC.GT.1) THEN
        IF(MASWRK) WRITE(IW,*)
     *       'SORRY, SVP MODEL DOES NOT RUN IN PARALLEL YET.'
        CALL ABRT
      END IF
C
C     ----- PRINT -SVP- PARAMETERS -----
C
      IF(MASWRK) THEN
      WRITE(IW,9999)
      IF(IEX.EQ.0) THEN
         WRITE(IW,9997) ISHAPE,ITRCAV,NDRCAV,
     2                  INTCAV,NPTLEB,NPTTHE,NPTPHI,
     3                  ITRNGR,IROTGR,IOPPRD,
     4                  LINEQ ,IRDRF , IPNRF
         WRITE(IW,9996) RHOISO,TOLCAV,RADSPH,TOLCHG,
     1                  DIELST,CVGLIN,CSDIAG,
     2                  ROTTHE,ROTPHI,ROTCHI,
     3                  TRANX ,TRANY ,TRANZ
      ELSEIF(IEX.EQ.1) THEN
         WRITE(IW,9995) METHOD,IQORP ,IOPPRD,NSTAGE,
     2                  IRENRM,ISHAPE,ITRCAV,NDRCAV,
     3                  ICVICK,IFRQUP,IRDRF ,INTCAV,
     4                  NPTLEB,NPTTHE,NPTPHI,IROTGR,
     5                  ITRNGR,IONSTR,ICRD  ,LINEQ ,
     6                  IPNRF ,IPRRF ,IPNDN ,IPRDN ,
     7                  IFPNCH,IEX
         WRITE(IW,9994) RHOISO,TOLCAV,RSTART,RSPACE,
     1                  RADSPH,TOLCHG,RADSLV,DIELST,DIELOP,
     2                  ROTTHE,ROTPHI,ROTCHI,
     3                  TRANX ,TRANY ,TRANZ ,
     4                  CVGLIN,CSDIAG,
     5                  ZIONP ,ZIONM ,TKELVN
      ENDIF
      END IF
C
C     ----- CHECK FOR LEGAL VALUES OF -SVP- PARAMETERS -----
C
      IF(   METHOD.LT.0.OR.METHOD.GT.3  .OR.  IQORP .LT.0.OR.IQORP .GT.2
     1 .OR. IOPPRD.LT.0.OR.IOPPRD.GT.3  .OR.  NSTAGE.LT.0.OR.NSTAGE.GT.2
     2 .OR. IRENRM.LT.0.OR.IRENRM.GT.1  .OR.  ISHAPE.LT.0.OR.ISHAPE.GT.3
     3 .OR. RHOISO.LE.ZERO              .OR.  ITRCAV.LT.0
     4 .OR. NDRCAV.LT.0.OR.NDRCAV.GT.2  .OR.  ICVICK.LT.0.OR.ICVICK.GT.2
     5                                  .OR.  RSPACE.EQ.ZERO
     6 .OR. TOLCAV.LE.ZERO              .OR.  RADSPH.EQ.ZERO
     7 .OR. IFRQUP.LT.0.OR.IFRQUP.GT.4  .OR.  IRDRF .LT.0.OR.IRDRF .GT.2
     8 .OR. DIELST.LT.ONE               .OR.  DIELOP.LT.ONE
     9 .OR. INTCAV.LT.0.OR.INTCAV.GT.1  .OR.  NPTLEB.LE.0
     1 .OR. NPTTHE.LE.0                 .OR.  NPTPHI.LE.0
     2 .OR. IROTGR.LT.0.OR.IROTGR.GT.3  .OR.  ITRNGR.LT.0.OR.ITRNGR.GT.6
     3 .OR. IONSTR.LT.0.OR.IONSTR.GT.2  .OR.  TKELVN.LE.ZERO
     4 .OR. ZIONP .LT.ZERO              .OR.  ZIONM .LT.ZERO
     5 .OR. ICRD  .LT.0.OR.ICRD.  GT.1  .OR.  LINEQ .LT.0.OR.LINEQ .GT.2
     6 .OR. CVGLIN.LE.ZERO              .OR.  CSDIAG.LT.ZERO
     7 .OR. IPNRF .LT.0.OR.IPNRF .GT.2  .OR.  ABS(IPRRF).GT.2
     8 .OR. IPNRF .LT.0.OR.IPNRF .GT.2  .OR.  ABS(IPRDN).GT.2
     9 .OR. IEX   .LT.0.OR.IEX   .GT.1) THEN
         IF(MASWRK) WRITE(IW,*)
     *        ' ILLEGAL VALUE REQUESTED FOR SOME SVP PARAMETER'
         CALL ABRT
      ENDIF
C
C     ----- CHECK FOR INCONSISTENCIES -----
C
      IF(IROTGR.EQ.1.AND.ITRNGR.NE.1) THEN
         IF(MASWRK) WRITE(IW,*)
     *       ' IROTGR.EQ.1 IS INCONSISTENT WITH ITRNGR.NE.1'
         CALL ABRT
      ENDIF
      IF(IROTGR.EQ.2.AND.ITRNGR.NE.2) THEN
         IF(MASWRK) WRITE(IW,*)
     *       ' IROTGR.EQ.2 IS INCONSISTENT WITH ITRNGR.NE.2'
         CALL ABRT
      ENDIF
C
C     ----- SET OPTIONS DEPENDENT ON -$SVP- PARAMETERS -----
C
      IF(RSTART.LT.ZERO) RSTART=-RSTART*TOBOHR
      IF(RSPACE.LT.ZERO) RSPACE=-RSPACE*TOBOHR
      IF(RADSPH.LT.ZERO) RADSPH=-RADSPH*TOBOHR
      IF(RADSLV.LT.ZERO) RADSLV=-RADSLV*TOBOHR
      ROTCAV(1) = ROTTHE
      ROTCAV(2) = ROTPHI
      ROTCAV(3) = ROTCHI
      TRNCAV(1) = TRANX
      TRNCAV(2) = TRANY
      TRNCAV(3) = TRANZ
      IF(INTCAV.EQ.0) NCAV=NPTLEB
      IF(INTCAV.EQ.1) NCAV=NPTTHE*NPTPHI
C
C     ----- INITIALIZE PARTS OF -SVPPRP- ... TO BE SAFE ... -----
C
      RHOEXT=ZERO
      CHGBAR=ZERO
      CHGRUL=ZERO
      CHGDEV=ZERO
C
C     ----- SHAPED CAVITY FROM INTERSECTING SPHERES -----
C
      IF(ISHAPE.EQ.2) THEN
         CALL SVP_CAVINP
      ENDIF
C
C     ----- SHAPED CAVITY FROM INPUT ... -----
C
      IF(ISHAPE.EQ.3) THEN
         CALL SVP_CAV_RD
      ENDIF
C
      IF(NCAV.LE.0.OR.NCAV.GT.MXCAV) THEN
         IF(MASWRK) WRITE(IW,*)
     *       ' ILLEGAL NUMBER OF SURFACE POINTS REQUESTED',NCAV
         CALL ABRT
      ENDIF
C
      NSVP=NCAV
      IQP=IQORP
C
C     ----- INITIALIZE REACTION FIELD TO ZERO -----
C
      DO I=1,NCAV
         CAVQ(I)=ZERO
         CAVP(I)=ZERO
      ENDDO
C
C     -- OBTAIN INITIAL CAVITY SURFACE AND REACTION FIELD EXTERNALLY --
C
      IF(IRDRF.GT.0) THEN
         CALL SVPIRF(IQP0,NCAV0,CAVC0,CAVQ0,CAVW0,CAVN0,CAVP0)
         IF(NCAV0.NE.NCAV) THEN
            IF(MASWRK) WRITE(IW,*)
     *          'NOTE THAT $SVPIRF NCAV0 .NE. $SVP NCAV'
         ENDIF
         IF(IQP0.NE.IQP) THEN
            IF(MASWRK) WRITE(IW,*)
     *          'NOTE THAT $SVPIRF IQP0  .NE. $SVP IQP'
         ENDIF
         SURF=ZERO
         DO I=1,NCAV
            SURF=SURF+CAVW0(I)
         ENDDO
      ENDIF
C
      RETURN
 9999 FORMAT(/,10X,13('-'),
     1       /,10X,'-SVP- SOLVENT',
     2       /,10X,13('-'))
 9997 FORMAT(/,1X,'----- INPUT GROUP $SVP -----',/,
     1 ' ISHAPE=',I4,' ITRCAV=',I4,' NDRCAV=',I4,/
     2 ' INTCAV=',I4,' NPTLEB=',I4,' NPTTHE=',I4,' NPTPHI=',I4,/,
     3 ' ITRNGR=',I4,' IROTGR=',I4,' IOPPRD=',I4,/,
     4 ' LINEQ =',I4,' IRDRF =',I4,' IPNRF =',I4)
 9996 FORMAT(' RHOISO=',F12.6,' TOLCAV=',D12.4,' RADSPH=',F12.6,
     1       ' TOLCHG=',D12.4,/,
     2       ' DIELST=',F12.6,' CVGLIN=',D12.4,' CSDIAG=',F12.6,/,
     3       ' ROTTHE=',F12.6,' ROTPHI=',F12.6,' ROTCHI=',F12.6,/,
     4       ' TRANX =',F12.6,' TRANY =',F12.6,' TRANZ =',F12.6,/)
 9995 FORMAT(/,1X,'----- INPUT GROUP $SVP -----',/,
     1 ' METHOD=',I4,' IQORP =',I4,' IOPPRD=',I4,' NSTAGE=',I4,/,
     2 ' IRENRM=',I4,' ISHAPE=',I4,' ITRCAV=',I4,' NDRCAV=',I4,/,
     3 ' ICVICK=',I4,' IFRQUP=',I4,' IRDRF =',I4,' INTCAV=',I4,/,
     4 ' NPTLEB=',I4,' NPTTHE=',I4,' NPTPHI=',I4,' IROTGR=',I4,/,
     5 ' ITRNGR=',I4,' IONSTR=',I4,' ICRD  =',I4,' LINEQ =',I4,/,
     6 ' IPNRF =',I4,' IPRRF =',I4,' IPNDN =',I4,' IPRDN =',I4,/,
     7 ' IFPNCH=',I4,' IEX   =',I4)
 9994 FORMAT(' RHOISO=',F12.6,' TOLCAV=',D12.4,' RSTART=',F12.6,/,
     1       ' RSPACE=',F12.6,' RADSPH=',F12.6,' TOLCHG=',D12.4,/,
     2       ' RADSLV=',F12.6,' DIELST=',F12.6,' DIELOP=',F12.6,/,
     3       ' ROTTHE=',F12.6,' ROTPHI=',F12.6,' ROTCHI=',F12.6,/,
     4       ' TRANX =',F12.6,' TRANY =',F12.6,' TRANZ =',F12.6,/,
     5       ' CVGLIN=',D12.4,' CSDIAG=',F12.6,/,
     6       ' ZIONP =',F12.6,' ZIONM =',F12.6,' TKELVN=',F12.6)
      END
C
C*MODULE SVPINP *DECK SVP_CAV_RD
      SUBROUTINE SVP_CAV_RD
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER  (MXATM=2000)
      PARAMETER  (MXCAV=5810)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TEN3=1.0D+03)
      PARAMETER (TOANGS=0.52917724924D+00,TOBOHR=ONE/TOANGS)
C
      CHARACTER*8 WORD
      CHARACTER*8 WRDSVP
      CHARACTER*8 WRDEND
C
      LOGICAL     SOME
      LOGICAL     OUT
      LOGICAL     DBUG
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON/SVPPAR/CAVC(3,MXCAV),CAVW(MXCAV),CAVQ(MXCAV),
     1              CAVN(3,MXCAV),SPHW(MXCAV),CAVP(MXCAV),
     2              CENCAV(3),RNTCAV(3),NCAV,IQP
      COMMON/SVPPRP/SURF,VOL,RHOEXT,CHGBAR,CHGRUL,CHGDEV
C
      DATA WRDSVP /' $SVPXYZ'/
      DATA WRDEND /' $END   '/
      DATA IUNIT  /1/
C
      DBUG=.FALSE.  .AND.  MASWRK
      OUT =.FALSE.  .AND.  MASWRK
      OUT =OUT.OR.DBUG
      SOME=.TRUE.
      SOME=(SOME.OR.OUT)  .AND.  MASWRK
      IF(SOME) THEN
         WRITE(IW,9997)
      ENDIF
C
C     ----- PRELIMINARY COUNT ... ----
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,WRDSVP,IEOF)
      IF(IEOF.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR, CANNOT FIND $SVPXYZ GROUP'
         CALL ABRT
      END IF
      GO TO 30
C
   20 IF(MASWRK) WRITE(IW,*) 'ERROR PROCESSING THE $SVPXYZ DATA GROUP'
      CALL ABRT
      STOP
C     ****
C
   30 CONTINUE
      NLINE=0
      IF(MASWRK) READ(IR,9996,END=20,ERR=20) WORD
      IF(WORD.EQ.WRDEND) GO TO 40
      NLINE=NLINE+1
      GO TO 30
C
   40 CONTINUE
      NCAV=NLINE
      IF(MASWRK) WRITE(IW,9995) NCAV
C
C     ----- READ SURFACE POINTS FROM -HND- FILE (IR) ... -----
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,WRDSVP,IEOF)
      IF(IEOF.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR, CANNOT FIND $SVPXYZ GROUP'
         CALL ABRT
      END IF
C
      IF(MASWRK) THEN
      DO ICAV=1,NCAV
         READ(IR,9998) ISVPXYZ,
     1                  CAVC(1,ICAV),CAVC(2,ICAV),CAVC(3,ICAV),
     2                  CAVN(1,ICAV),CAVN(2,ICAV),CAVN(3,ICAV),
     3                  CAVW(  ICAV)
         IF(ISVPXYZ.NE.ICAV) THEN
            IF(MASWRK) WRITE(IW,*) 'LINE NO. MISMATCH IN $SVPXYZ'
            CALL ABRT
         END IF
         SCALE=1.00D+00
         CAVC(1,ICAV)=SCALE   *CAVC(1,ICAV)
         CAVC(2,ICAV)=SCALE   *CAVC(2,ICAV)
         CAVC(3,ICAV)=SCALE   *CAVC(3,ICAV)
         CAVW(  ICAV)=SCALE**2*CAVW(  ICAV)
         IF(IUNIT.GT.0) THEN
            CAVC(1,ICAV)=CAVC(1,ICAV)*TOBOHR
            CAVC(2,ICAV)=CAVC(2,ICAV)*TOBOHR
            CAVC(3,ICAV)=CAVC(3,ICAV)*TOBOHR
            CAVW(  ICAV)=CAVW(  ICAV)*TOBOHR**2
         ENDIF
      ENDDO
      END IF
C        NOTE THAT PARALLELIZATION REQUIRES BROADCASTING THE ABOVE.
C
      IF(OUT) THEN
         WRITE(IW,9999) WRDSVP
         DO ICAV=1,NCAV
            WRITE(IW,9998) ICAV,
     1                   CAVC(1,ICAV),CAVC(2,ICAV),CAVC(3,ICAV),
     2                   CAVN(1,ICAV),CAVN(2,ICAV),CAVN(3,ICAV),
     3                   CAVW(  ICAV)
         ENDDO
         WRITE(IW,9999) WRDEND
      ENDIF
C
      IF(MASWRK) WRITE(IW,9989) SCALE
C
C     ----- MOLECULAR SURFACE ... -----
C
      SRFMOL=ZERO
      DO ICAV=1,NCAV
         SRFMOL=SRFMOL+CAVW(ICAV)
      ENDDO
C
      SURF  =SRFMOL
      VOL   =ZERO
C
      SRFMOL=SRFMOL*TOANGS**2
C
      IF(MASWRK) WRITE(IW,9994) SRFMOL
C
C     ----- CAVITY/DISPERSION FREE ENERGY ---
C           SITKOFF, SHARP, AND HONIG,
C           J.PHYS.CHEM. 98, 1978 (1994)
C
      CAVDSP=0.860D+00+0.005D+00*SRFMOL
      IF(SOME) THEN
         WRITE(IW,9993) CAVDSP
      ENDIF
C
C     ----- CHECK SOLUTE-CAVITY PROXIMITY ... -----
C
      IF(DBUG.AND.SOME) THEN
         WRITE(IW,9991)
      ENDIF
      DO IAT=1,NAT
         DMIN=TEN3
         MIND=NCAV+1
         DO ICAV=1,NCAV
            DIST=SQRT((C(1,IAT)-CAVC(1,ICAV))**2+
     1                (C(2,IAT)-CAVC(2,ICAV))**2+
     2                (C(3,IAT)-CAVC(3,ICAV))**2)*TOANGS
            IF(DIST.LT.DMIN) THEN
               DMIN=DIST
               MIND=ICAV
            ENDIF
         ENDDO
         IF(DBUG.AND.SOME) THEN
            WRITE(IW,9990) IAT,MIND,DMIN
         ENDIF
      ENDDO
C
      IF(MASWRK) WRITE(IW,9992)
      CALL FLSHBF(IW)
C
      RETURN
 9999 FORMAT(A8)
 9998 FORMAT(I5,3F13.8,3F9.6,F9.5)
 9997 FORMAT(10X,13('-'),/,10X,'SHAPED CAVITY',/,10X,13('-'))
 9996 FORMAT(A8)
 9995 FORMAT(' NUMBER OF CAVITY SURFACE POINTS = ',I10)
 9994 FORMAT(' MOLECULAR SURFACE = ',F10.3,' ANGSTROM**2')
 9993 FORMAT(' G(CAV/DISP)       = ',F10.3,' KCAL/MOL')
 9992 FORMAT(/)
 9991 FORMAT(' SOLUTE/CAVITY PROXIMITY',/,1X,23('-'))
 9990 FORMAT(' IAT = ',I5,' IPT_CAV = ',I7,' DIST = ',F8.4)
 9989 FORMAT(' SCALING FACTOR = ',F7.3)
      END
C*MODULE SVP    *DECK SVP_CAVINP
      SUBROUTINE SVP_CAVINP
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (MXATM=2000)
      LOGICAL    SOME
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON/CAVPAR/MINBEM,MAXBEM,IFICOS
      COMMON/CAVRAD/RADIUS(MXATM)
      COMMON/CAVSLV/RSOLV
      COMMON/CAVLAY/RLAYER,NLAYER
      COMMON/CAVPRT/IPRT
      COMMON/SVPDBG/IDBG
      DIMENSION    VDWR(103)
      DIMENSION    VANDER(36)
      DIMENSION    RADII(MXATM)
      EQUIVALENCE (RADII(1),RADIUS(1))
      DATA MXPASS  /6/
      DATA RADTOL  /0.1D+00/
      DATA NPFLG   /0/
      DATA ICALL   /1/
C
      SAVE ICALL
C
C     ----- VAN DER WAALS RADII IN ANGSTROMS -----
C
      DATA DEF  /1.80D+00/
      DATA VDWR /
     1   0.80,0.49,0.00,0.00,0.00,1.65,1.55,1.50,1.50,0.00,
     2   2.30,1.70,2.05,2.10,1.85,1.80,1.80,0.00,2.80,2.75,
     3   0.00,0.00,1.20,0.00,0.00,0.00,2.70,0.00,0.00,0.00,
     4   0.00,0.00,0.00,1.90,1.90,0.00,0.00,0.00,0.00,1.55,
     5   0.00,1.64,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,
     6   0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,
     7   0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,
     8   0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,
     9   0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,
     1   0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,1.65,
     2   0.00,0.00,0.00/
C
C     VDW RADII FROM MERZ/KOLLMAN/SINGH
C
      DATA (VANDER(I),I=1,36)
     1                  /1.20D+00,1.20D+00,1.37D+00,1.45D+00,
     2 1.45D+00,1.50D+00,1.50D+00,1.40D+00,1.35D+00,1.30D+00,
     3 1.57D+00,1.36D+00,1.24D+00,1.17D+00,1.80D+00,1.75D+00,
     4 1.70D+00,19*0.0D+00/
C
C     VDW RADII FROM GAVEZZOTTI (J.AM.CHEM.SOC. 105, 5220 (1983))
C     SUPPLEMENTED BY VALUES FOR BORON AND THIRD AND FOURTH ROWS
C     FROM M.A.SPACKMAN (J.CHEM.PHYS. 85, 6579 (1986))
C
C     DATA (VANDER(I),I=1,36)
C    1                  /1.20D+00,0.00D+00,0.00D+00,0.00D+00,
C    2 1.85D+00,1.50D+00,1.50D+00,1.40D+00,1.35D+00,0.00D+00,
C    3 0.00D+00,0.00D+00,2.07D+00,2.05D+00,1.96D+00,1.89D+00,
C    4 1.80D+00,19*0.0D+00/
C
C     VDW RADII FROM BRENEMAN & WIBERG
C
C     DATA (VANDER(I),I=1,36)
C    1                  /1.45D+00,1.45D+00,1.50D+00,1.50D+00,
C    2 1.50D+00,1.50D+00,1.70D+00,1.70D+00,1.70D+00,1.70D+00,
C    3 2.00D+00,2.00D+00,2.00D+00,2.00D+00,2.00D+00,2.00D+00,
C    4 2.00D+00,19*0.0D+00/
C
C         NAMELIST SIMULATION
C
      PARAMETER (NNAM=9)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA CAVKWD/8HCAVITY  /
      DATA QNAM/8HMINBEM  ,8HMAXBEM  ,8HIFICOS  ,8HRADIUS  ,8HRADII   ,
     *          8HRSOLV   ,8HRLAYER  ,8HNLAYER  ,8HNPFLG   /
      DATA KQNAM/1,1,1,3,3,   3,3,1,1/
C
      SOME=.FALSE.
      SOME=SOME.OR.ICALL.EQ.1
      SOME=SOME.OR.IDBG.GE.1
C
      IFICOS=0
      MINBEM=2
      MAXBEM=3
      RSOLV =ZERO
      RLAYER=ZERO
      NLAYER=1
C
      DO I=1,36
         VDWR(I)=VANDER(I)
      ENDDO
      DO IAT=1,NAT
         NUCZ = INT(ZAN(IAT) + IZCORE(IAT))
         RAD=VDWR(NUCZ)
         IF(RAD.EQ.ZERO) THEN
            RAD=DEF
         ENDIF
         RADIUS(IAT)=RAD
      ENDDO
C
C     ----- READ INPUT GROUP -CAVITY- -----
C
      JRET=0
      CALL NAMEIO(IR,JRET,CAVKWD,NNAM,QNAM,KQNAM,
     *            MINBEM,MAXBEM,IFICOS,RADIUS,RADII,RSOLV,RLAYER,
     *            NLAYER,NPFLG,
     *                    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)
C
      IF(JRET.EQ.2) THEN
        IF (MASWRK) WRITE(IW,*)
     *       'ERROR READING $CAVITY INPUT, CHECK ABOVE'
        CALL ABRT
      END IF
C
      IF(SOME) THEN
         WRITE(IW,9999)
      ENDIF
C
C     ----- TESSELATION -----
C
      IF(MINBEM.LT.1.OR.MINBEM.GT.MXPASS) THEN
         MINBEM=2
      ENDIF
      IF(MAXBEM.LT.1.OR.MAXBEM.GT.MXPASS) THEN
         MAXBEM=4
      ENDIF
      IF(MAXBEM.LT.MINBEM) THEN
         MAXBEM=MINBEM+2
      ENDIF
C
      IFICOS=ABS(IFICOS)
C
      DO IAT=1,NAT
         RADIUS(IAT)=ABS(RADIUS(IAT))
         IF(RADIUS(IAT).LT.RADTOL) THEN
            RADIUS(IAT)=ZERO
         ENDIF
      ENDDO
C
      RSOLV=ABS(RSOLV)
C
      NLAYER=ABS(NLAYER)
      RLAYER=ABS(RLAYER)
      IF(NLAYER.EQ.0) THEN
         NLAYER=1
      ELSEIF(NLAYER.GT.1) THEN
         IF(RLAYER.EQ.ZERO) THEN
            NLAYER=1
         ENDIF
      ENDIF
C
C     ----- PRINT ... -----
C
      IF(SOME) THEN
         WRITE(IW,9991) MINBEM
         WRITE(IW,9992) MAXBEM
         IF(IFICOS.EQ.0) THEN
            WRITE(IW,9990)
         ELSE
            WRITE(IW,9989)
         ENDIF
         WRITE(IW,9994) RSOLV
         WRITE(IW,9988)
         DO IAT=1,NAT
            NUCZ = INT(ZAN(IAT) + IZCORE(IAT))
            WRITE(IW,9987) IAT,NUCZ,RADIUS(IAT)
         ENDDO
         WRITE(IW,9986) NLAYER
         IF(NLAYER.GT.1) THEN
            WRITE(IW,9985) RLAYER
         ENDIF
      ENDIF
C
C     ----- SET PRINT FLAG -----
C
      NPFLG=ABS(NPFLG)
      IF(NPFLG.GT.0) THEN
         IPRT=2
      ELSE
         IF(ICALL.EQ.1) THEN
            IPRT=1
         ENDIF
         IF(IDBG.GT.0) THEN
            IPRT=2
         ENDIF
      ENDIF
C
C     ----- CREATE BOUNDARY ELEMENTS OF CAVITY -----
C
      CALL SVP_CAVMAK
C
C     ----- RESET PRINT FLAG -----
C
      IF(NPFLG.EQ.0) THEN
         IPRT=0
      ENDIF
C
      ICALL=2
C
      RETURN
 9999 FORMAT(/,10X,8('-'),
     1       /,10X,'-CAVITY-',
     2       /,10X,8('-'))
 9994 FORMAT(' SOLVENT RADIUS (ANG.)      = ',F7.3)
 9992 FORMAT(' -BEM- HIGH LEVEL           = ',I3)
 9991 FORMAT(' -BEM- LOW  LEVEL           = ',I3)
 9990 FORMAT(' -BEM- FROM -OCTAHEDRAL-')
 9989 FORMAT(' -BEM- FROM -ICOSAHEDRAL-')
 9988 FORMAT('  ATOMIC RADII  (ANG.)        ',/,
     1       '  --------------------')
 9987 FORMAT(I5,I5,F7.3)
 9986 FORMAT(' NLAYER                     = ',I3)
 9985 FORMAT(' RLAYER (ANG.)              = ',F7.3)
      END
C*MODULE SVP    *DECK SVPIRF
      SUBROUTINE SVPIRF(IQP,NPT,XYZ,CHG,WGT,XYZN,DIP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- READ FORMATTED PSEUDONAMELIST -$SVPIRF- -----
C
C     IQP  = TYPE OF REACTION FIELD FOUND
C     NPT  = NUMBER OF SURFACE GRID POINTS FOUND
C     XYZ  = GRID POINT (X,Y,Z) COORDINATES
C     CHG  = POINT CHARGE ON SURFACE
C     WGT  = SURFACE AREA; SIGN IS THE SIGN OF ZN JUST BELOW
C     XYZN = UNIT VECTORS NORMAL TO SURFACE
C            ONLY XN,YN COMPONENTS ARE READ; |ZN|=SQRT(1-XN**2-YN**2)
C     DIP  = POINT DIPOLE ON SURFACE
C
      PARAMETER (ZERO=0.0D+00)
      CHARACTER*8 RFWRD,ENDWRD,WORD
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DIMENSION XYZ(3,*),CHG(*),WGT(*),XYZN(3,*),DIP(*)
C
      DATA RFWRD     /' $SVPIRF'/
      DATA ENDWRD    /' $END   '/
C
C     ----- FIND INPUT FOR -$SVPIRF- -----
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,RFWRD,IEOF)
      IF(IEOF.NE.0) GO TO 70
C
C     ----- READ ALL -$SVPIRF- CARDS TO GET -NPT- -----
C
      NLINES=0
   20 NLINES=NLINES+1
      READ(IR,9099,END=80,ERR=90) WORD
      IF(WORD.NE.ENDWRD) GOTO 20
      NLINES=NLINES-1
      NPT=NLINES/2
C
C     ----- REREAD ALL -$SVPIRF- CARDS TO GET DATA -----
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,RFWRD,IEOF)
      CHGSUM=ZERO
      DIPSUM=ZERO
      IF(MASWRK) THEN
      DO I=1,NPT
         READ(IR,9096) IPT1,(XYZ(K,I),K=1,3),CHG(I),WGT(I)
         READ(IR,9096) IPT2,(XYZN(K,I),K=1,3),DIP(I)
         IF(IPT1.NE.I.OR.IPT2.NE.I) THEN
            IF(MASWRK) WRITE(IW,*)
     *          ' POINT OUT OF ORDER IN -$SVPIRF- AT I=',I
            CALL ABRT
         ENDIF
         CHGSUM=CHGSUM+CHG(I)
         DIPSUM=DIPSUM+DIP(I)
      ENDDO
      ENDIF
C        NOTE THAT PARALLELIZATION REQUIRES BROADCASTING THE ABOVE.
C
C     ----- DETERMINE TYPE OF REACTION FIELD PRESENT -----
C
      IF(CHGSUM.EQ.ZERO.AND.DIPSUM.EQ.ZERO) IQP=-1
      IF(CHGSUM.NE.ZERO.AND.DIPSUM.EQ.ZERO) IQP=0
      IF(CHGSUM.EQ.ZERO.AND.DIPSUM.NE.ZERO) IQP=1
      IF(CHGSUM.NE.ZERO.AND.DIPSUM.NE.ZERO) IQP=2
C
      IF(MASWRK) WRITE(IW,9999) NPT
      IF(IQP.EQ.-1) THEN
         IF(MASWRK) WRITE(IW,9998)
         IQP=0
      ENDIF
      RETURN
C
   70 IF(MASWRK) WRITE(IW,*) ' NO NAMELIST -$SVPIRF- FOUND. STOP'
      CALL ABRT
   80 IF(MASWRK) WRITE(IW,*) ' EOF READING NAMELIST $SVPIRF'
      CALL ABRT
   90 IF(MASWRK) WRITE(IW,*) ' ERROR READING NAMELIST $SVPIRF'
      CALL ABRT
      STOP
C
 9999 FORMAT(/,' INITIAL GUESS SVP REACTION FIELD READ FROM CARDS',
     1         ' WITH -NPT- =',I6)
 9998 FORMAT(  ' INITIAL GUESS SVP REACTION FIELD IS FOUND TO BE ZERO')
 9099 FORMAT(A8)
 9096 FORMAT(I5,5F15.10)
      END
C*MODULE SVPINP  *DECK SVP_CAV_PUN
      SUBROUTINE SVP_CAV_PUN
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER   (MXATM=2000)
      PARAMETER   (MXCAV=5810)
      PARAMETER (TEN3=1.0D+03)
      PARAMETER (TOANGS=0.52917724924D+00)
      CHARACTER*8  WRDSVP
      CHARACTER*8  WRDEND
      LOGICAL      OUT
      LOGICAL      SOME
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON/SVPCAV/IFPNCH
      COMMON/SVPPAR/CAVC(3,MXCAV),CAVW(MXCAV),CAVQ(MXCAV),
     1              CAVN(3,MXCAV),SPHW(MXCAV),CAVP(MXCAV),
     2              CENCAV(3),RNTCAV(3),NCAV,IQP
      DATA WRDSVP /' $SVPXYZ'/
      DATA WRDEND /' $END   '/
C
      OUT =.FALSE.
      SOME=.FALSE.
      SOME=SOME.OR.OUT
C
C     ----- PRINT SOLUTE-CAVITY PROXIMITY ... -----
C
      IF(SOME) THEN
         WRITE(IW,9991)
      ENDIF
      DO IAT=1,NAT
         DMIN=TEN3
         MIND=NCAV+1
         DO ICAV=1,NCAV
            DIST=SQRT((C(1,IAT)-CAVC(1,ICAV))**2+
     1                (C(2,IAT)-CAVC(2,ICAV))**2+
     2                (C(3,IAT)-CAVC(3,ICAV))**2)*TOANGS
            IF(DIST.LT.DMIN) THEN
               DMIN=DIST
               MIND=ICAV
            ENDIF
         ENDDO
         IF(SOME) THEN
            WRITE(IW,9990) IAT,MIND,DMIN
         ENDIF
      ENDDO
C
      IF(IFPNCH.NE.1) THEN
         RETURN
      ENDIF
C
C     ----- WRITE SURFACE POINTS ON -PUN- FILE ... -----
C
      WRITE(IP,9999) WRDSVP
      DO ICAV=1,NCAV
         WRITE(IP,9998) ICAV,
     1                   CAVC(1,ICAV),CAVC(2,ICAV),CAVC(3,ICAV),
     2                   CAVN(1,ICAV),CAVN(2,ICAV),CAVN(3,ICAV),
     3                   CAVW(  ICAV)
      ENDDO
      WRITE(IP,9999) WRDEND
C
      IF(OUT) THEN
         WRITE(IW,9999) WRDSVP
         DO ICAV=1,NCAV
            WRITE(IW,9998) ICAV,
     1                   CAVC(1,ICAV),CAVC(2,ICAV),CAVC(3,ICAV),
     2                   CAVN(1,ICAV),CAVN(2,ICAV),CAVN(3,ICAV),
     3                   CAVW(  ICAV)
         ENDDO
         WRITE(IW,9999) WRDEND
      ENDIF
C
      RETURN
 9999 FORMAT(A8)
 9998 FORMAT(I5,3F13.8,3F9.6,F9.5)
 9991 FORMAT(' SOLUTE/CAVITY PROXIMITY',/,1X,23('-'))
 9990 FORMAT(' IAT = ',I5,' IPT_CAV = ',I7,' DIST = ',F8.4,' ANGS.')
      END
C*MODULE SVPINP  *DECK SVPPTW
      SUBROUTINE SVPPTW(XYZ,W,NPTTHE,NPTPHI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     --- (THETA,PHI) INTEGRATION POINTS AND WEIGHTS ON UNIT SPHERE ---
C
      PARAMETER (PT5=0.5D+00, ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00)
      PARAMETER (MXTHE=128)
      PARAMETER (MXPHI=256)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
C
      DIMENSION XYZ(3,*),W(*)
      DIMENSION GLROOT(MXTHE),GLWGHT(MXTHE)
      DIMENSION COSPHI(MXPHI),SINPHI(MXPHI)
C
      PI=FOUR*ATAN(ONE)
      FOURPI=FOUR*PI
C
      IF(NPTTHE.GT.MXTHE) THEN
         WRITE(IW,*) 'LOCAL DIMENSION TOO SMALL FOR NPTTHE. STOP'
         CALL ABRT
      ENDIF
      IF(NPTPHI.GT.MXPHI) THEN
         WRITE(IW,*) 'LOCAL DIMENSION TOO SMALL FOR NPTPHI. STOP'
         CALL ABRT
      ENDIF
C
C     ----- GAUSS-LEGENDRE ROOTS AND WEIGHTS FOR THETA -----
C
      CALL GLGQUD(-ONE,+ONE,GLROOT,GLWGHT,NPTTHE)
C
C     ----- ANGLE SINES AND COSINES FOR PHI -----
C
      DPHI=TWO*PI/NPTPHI
      DO JPHI=1,NPTPHI
         PHII=DPHI*(JPHI-PT5)
         COSPHI(JPHI)=COS(PHII)
         SINPHI(JPHI)=SIN(PHII)
      ENDDO
C
C     ----- COMBINE (THETA,PHI) FACTORS -----
C
      IJ=0
      DO ITHE=1,NPTTHE
         COSTHE=GLROOT(ITHE)
         SINTHE=SQRT(ONE-COSTHE**2)
         WTHE  =GLWGHT(ITHE)
         DO JPHI=1,NPTPHI
            IJ=IJ+1
            XYZ(1,IJ)=COSPHI(JPHI)*SINTHE
            XYZ(2,IJ)=SINPHI(JPHI)*SINTHE
            XYZ(3,IJ)=             COSTHE
            WPHI     =DPHI
            W(IJ)    =WTHE*WPHI/FOURPI
         ENDDO
      ENDDO
C
      RETURN
      END
C*MODULE SVPINP  *DECK SVPWRF
      SUBROUTINE SVPWRF(IQP,NPT,XYZ,CHG,WGT,XYZN,DIP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- PRINT REACTION FIELD POINTS AND CHARGES  -----
C
C     IQP   = 0 FOR CHG; 1 FOR DIP; 2 FOR BOTH
C     NPT   = NUMBER OF SURFACE GRID POINTS
C     XYZ   = GRID POINT (X,Y,Z) COORDINATES
C     CHG   = POINT CHARGE ON SURFACE
C     WGT   = SURFACE AREA
C     XYZN  = UNIT VECTORS NORMAL TO SURFACE
C     DIP   = POINT DIPOLE ON SURFACE
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      DIMENSION XYZ(3,*),CHG(*),WGT(*),XYZN(3,*),DIP(*)
C
      IF(IQP.EQ.0) WRITE(IW,9998)
      IF(IQP.GT.0) WRITE(IW,9997)
C
      DO I=1,NPT
         SIG=CHG(I)/WGT(I)
         XMU=DIP(I)/WGT(I)
         IF(IQP.EQ.0) WRITE(IW,9996)
     1      I,(XYZ(K,I),K=1,3),WGT(I),CHG(I),SIG,
     2        (XYZN(K,I),K=1,3)
         IF(IQP.GT.0) WRITE(IW,9996)
     1      I,(XYZ(K,I),K=1,3),WGT(I),CHG(I),SIG,
     2        (XYZN(K,I),K=1,3),DIP(I),XMU
      ENDDO
C
 9998 FORMAT(10X,'X',10X,'Y',10X,'Z',10X,'W',10X,'Q',10X,'SIG',
     1       8X,'XN', 9X,'YN', 9X,'ZN')
 9997 FORMAT(10X,'X',10X,'Y',10X,'Z',10X,'W',10X,'Q',10X,'SIG',
     1       8X,'XN', 9X,'YN', 9X,'ZN', 9X,'P',10X,'MU')
 9996 FORMAT(I5,11F11.6)
      RETURN
      END
C*MODULE SVPINP  *DECK SVPPRF
      SUBROUTINE SVPPRF(IQP,NPT,XYZ,CHG,WGT,XYZN,DIP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- PUNCH REACTION FIELD POINTS AND CHARGES FOR -$SVPIRF- -----
C
C     IQP   = 0 FOR CHG; 1 FOR DIP; 2 FOR BOTH
C     NPT   = NUMBER OF SURFACE GRID POINTS
C     XYZ   = GRID POINT (X,Y,Z) COORDINATES
C     CHG   = POINT CHARGE ON SURFACE
C     WGT   = SURFACE AREA
C     XYZN  = UNIT VECTORS NORMAL TO SURFACE
C     DIP   = POINT DIPOLE ON SURFACE
C
      LOGICAL OUT
      CHARACTER*8 RFWRD,ENDWRD
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      DIMENSION XYZ(3,*),CHG(*),WGT(*),XYZN(3,*),DIP(*)
      DATA RFWRD     /' $SVPIRF'/
      DATA ENDWRD    /' $END   '/
C
      OUT =.FALSE.
      IF(OUT) THEN
         WRITE(IW,*) 'IN -SVPPRF- ... PUNCH $SVPIRF ...'
      ENDIF
C
      WRITE(IP,8998) RFWRD
C
      DO I=1,NPT
         IF(IQP.EQ.0) THEN
            WRITE(IP,8997) I,(XYZ(K,I),K=1,3),CHG(I),WGT(I)
            WRITE(IP,8997) I,(XYZN(K,I),K=1,3)
         ENDIF
         IF(IQP.GT.0) THEN
            WRITE(IP,8997) I,(XYZ(K,I),K=1,3),CHG(I),WGT(I)
            WRITE(IP,8997) I,(XYZN(K,I),K=1,3),DIP(I)
         ENDIF
      ENDDO
C
      WRITE(IP,8998) ENDWRD
C
 8998 FORMAT(A8)
 8997 FORMAT(I5,5F15.10)
      RETURN
      END
C*MODULE SVPINP  *DECK SVPIDN
      SUBROUTINE SVPIDN(NPT,POT,FLD,RHO,DRN,ALP,WGT,DRG,OUT,DDN,DDL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- READ FORMATTED PSEUDONAMELIST -$SVPIDN- -----
C
C     NPT    = NUMBER OF SURFACE GRID POINTS FOUND
C     POT    = ELECTROSTATIC POTENTIAL
C     FLD    = NORMAL ELECTRIC FIELD
C     RHO    = DENSITY
C     DRN    = DERIVATIVE OF DENSITY IN NORMAL DIRECTION
C     ALP    = NEGATIVE OF (RADIAL DENSITY DERIVATIVE)/DENSITY
C     WGT    = WEIGHT OF POINT
C     DRG    = |DEL RHO|
C     OUT    = SOLUTE CHARGE PENETRATING OUTSIDE CAVITY
C     DDN    = 2ND DERIVATIVE OF DENSITY IN NORMAL DIRECTION
C     DDL    = LAPLACIAN OF DENSITY
C
      CHARACTER*8 DNWRD,ENDWRD,WORD
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DIMENSION POT(*),FLD(*)
      DIMENSION RHO(*),DRN(*),ALP(*),WGT(*),DRG(*)
      DIMENSION OUT(*),DDN(*),DDL(*)
      DATA DNWRD  /' $SVPIDN'/
      DATA ENDWRD /' $END   '/
C
C     ----- FIND INPUT FOR -$SVPIDN- -----
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,DNWRD,IEOF)
      IF(IEOF.GT.0) GO TO 70
C
C     ----- READ ALL -$SVPIDN- CARDS TO GET NPT -----
C
      NPT=0
   20 NPT=NPT+1
      IF(MASWRK) READ(IR,9099,END=80,ERR=90) WORD
      IF(WORD.NE.ENDWRD) GOTO 20
      NPT=NPT-1
C
C     ----- REREAD ALL -$SVPIDN- CARDS TO GET DATA -----
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,DNWRD,IEOF)
   30 CONTINUE
      IF(MASWRK) READ(IR,9099,END=70,ERR=70) WORD
      IF(WORD.NE.DNWRD) GO TO 30
      DO I=1,NPT
         READ(IR,9098,END=80,ERR=90)
     1                  IPT1,POT(I),FLD(I),RHO(I),DRN(I),ALP(I)
         READ(IR,9098,END=80,ERR=90)
     1                  IPT2,WGT(I),DRG(I),OUT(I),DDN(I),DDL(I)
         IF(IPT1.NE.I.OR.IPT2.NE.I) THEN
            IF(MASWRK) WRITE(IW,*)
     *         ' POINT FOUND OUT OF ORDER IN -$SVPIDN- AT I=',I
            CALL ABRT
         ENDIF
      ENDDO
C
      IF(MASWRK) WRITE(IW,9999)
      RETURN
C
   70 IF(MASWRK) WRITE(IW,*) ' NO NAMELIST -$SVPIDN- FOUND. STOP'
      CALL ABRT
   80 IF(MASWRK) WRITE(IW,*) ' EOF READING NAMELIST $SVPIDN'
      CALL ABRT
   90 IF(MASWRK) WRITE(IW,*) ' ERROR READING NAMELIST $SVPIDN'
      CALL ABRT
      STOP
 9999 FORMAT(/,' SVP DENSITY-RELATED PROPERTIES READ FROM CARDS')
 9099 FORMAT(A8)
 9098 FORMAT(I5,5E15.10)
      END
C*MODULE SVPINP  *DECK SVPODN
      SUBROUTINE SVPODN(IPRT,IPUN,NCAV,CAVC,CAVW,CAVN,SPHW,CENCAV)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     ----- ADDITIONAL SOLUTE DENSITY PROPERTIES ON CAVITY SURFACE -----
C
C     IPRT = 0 NO PRINT; <0 RESULTS AT EACH POINT; >0 INTEGRATED RESULTS
C     IPUN = 0 NO PUNCH; .NE.0 RESULTS AT EACH POINT
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER   (MXATM=2000)
C
      CHARACTER*8 DNWRD,ENDWRD
C
      LOGICAL     DBUG,OUT,SOME
      LOGICAL     CAVPOT,CAVFLD
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /FMCOM / X(1)
      DIMENSION CAVC(3,*),CAVW(*),CAVN(3,*),SPHW(*),CENCAV(*)
      DIMENSION PRP(10),PRPMIN(10),PRPMAX(10),PRPTOT(10),
     1          PRPSCL(10),PRPAVG(10)
C
      DATA DNWRD  /' $SVPIDN'/
      DATA ENDWRD /' $END   '/
      DATA PT01/1.0D-02/
      DATA SML,BIG    /-1.0D+09,1.0D+09/
C
      DBUG=.FALSE.
      OUT =.FALSE.
      OUT =OUT.OR.DBUG
      SOME=NPRINT.NE.-5
      SOME=SOME.OR.OUT
      IF(.NOT.SOME) RETURN
C
      IF(IPRT.NE.0) WRITE(IW,8999)
      IF(IPRT.NE.0) WRITE(IW,8998)
      IF(IPUN.NE.0) WRITE(IP,7999)
      IF(IPUN.NE.0) WRITE(IP,7998) DNWRD
C
      L2=(NUM*(NUM+1))/2
C
C     ----- SET POINTERS FOR PARTITIONING OF CORE -----
C
C     -DMAT-   AT X(I10)
C     -ELP-    AT X(I20)
C     -ELF-    AT X(I30)
C     -RHO-    AT X(I40)
C     -DRHO-   AT X(I50)
C     -DDRHO-  AT X(I60)
C
C     --- GET FAST MEMORY -----
C
      CALL VALFM(LOADFM)
      I10=LOADFM+1
      I20=I10+L2
      I30=I20+NCAV
      I40=I30+NCAV
      I50=I40+NCAV
      I60=I50+3*NCAV
      LAST=I60+6*NCAV
      LAST=MAX(LAST,I20+L2)
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C     ----- GET TOTAL DENSITY MATRIX OVER AO'S -----
C
      CALL DENDD1(X(I10),X(I20),L2)
C
C     ----- EVALUATE SOLUTE DENSITY PROPERTIES AT CAVITY SURFACE -----
C
C     PRP(1)  ELECTROSTATIC POTENTIAL
C     PRP(2)  NORMAL ELECTRIC FIELD
C     PRP(3)  ELECTRONIC DENSITY
C     PRP(4)  DERIVATIVE OF DENSITY IN NORMAL DIRECTION = NORMAL*(DEL RH
C     PRP(5)  EFFECTIVE EXPONENT OF EXPONENTIAL FIT TO NORMAL DENSITY /
C     PRP(6)  WEIGHT, I.E., SURFACE AREA / 100
C     PRP(7)  |DEL RHO| (REDUNDANT WITH DRHON ON ISODENSITY SURFACE)
C     PRP(8)  SOLUTE CHARGE PENETRATING OUTSIDE CAVITY
C     PRP(9)  2ND DERIVATIVE OF DENSITY IN NORMAL DIRECTION
C     PRP(10) LAPLACIAN OF DENSITY (REDUNDANT WITH DDRHON ON SPHERICAL S
C
      CAVPOT=.TRUE.
      CAVFLD=.TRUE.
      IADDL=I50
      ILAST=I60+6*NCAV
C
C     --- GRAB THE AS YET UNUSED SPACE IN DRHO FOR SCRATCH IN SVPEPF ---
C
      CALL SVPEPF(X(I10),IADDL,ILAST,NCAV,CAVC,CAVN,
     1            CAVPOT,X(I20),CAVFLD,X(I30))
C
      NDER=2
      CALL SVP_ELDEND(NCAV,CAVC,X(I40),X(I50),X(I60),X(I10),NDER)
C
      DO K=1,10
         PRPMIN(K)=BIG
         PRPMAX(K)=SML
         PRPTOT(K)=ZERO
         PRPSCL(K)=ONE
      ENDDO
      WGTTOT=ZERO
      SPHTOT=ZERO
      PRPSCL(5)=PT01
      PRPSCL(6)=PT01
      DO I=1,NCAV
         ELP=X(I+I20-1)
         ELFN=X(I+I30-1)
         RHO=X(I+I40-1)
         DRHOX=X(1+3*(I-1)+I50-1)
         DRHOY=X(2+3*(I-1)+I50-1)
         DRHOZ=X(3+3*(I-1)+I50-1)
         DDRHXX=X(1+6*(I-1)+I60-1)
         DDRHYY=X(2+6*(I-1)+I60-1)
         DDRHZZ=X(3+6*(I-1)+I60-1)
         DDRHXY=X(4+6*(I-1)+I60-1)
         DDRHXZ=X(5+6*(I-1)+I60-1)
         DDRHYZ=X(6+6*(I-1)+I60-1)
         XN=CAVN(1,I)
         YN=CAVN(2,I)
         ZN=CAVN(3,I)
         WGT=CAVW(I)
         DRN=XN*DRHOX+YN*DRHOY+ZN*DRHOZ
         DG2=DRHOX**2+DRHOY**2+DRHOZ**2
C
C        FIT RADIAL DENSITY AND ITS FIRST AND SECOND RADIAL DERIVATIVES
C        AT R0 TO BB*[1+CC*(R-R0)**2]*EXP[-AA*(R-R0)] TO GET CHG PEN
C
         RX=CAVC(1,I)-CENCAV(1)
         RY=CAVC(2,I)-CENCAV(2)
         RZ=CAVC(3,I)-CENCAV(3)
         RR2=RX**2+RY**2+RZ**2
         RR=SQRT(RR2)
         RHO0=RHO
         DRHO0=(RX*DRHOX+RY*DRHOY+RZ*DRHOZ)/RR
         DDRHO0=(RX*DDRHXX*RX+RY*DDRHYY*RY+RZ*DDRHZZ*RZ+
     1          (RX*DDRHXY*RY+RX*DDRHXZ*RZ+RY*DDRHYZ*RZ)*TWO)/RR2
         BB=RHO0
         AA=-DRHO0/BB
         CC=(DDRHO0/BB-AA**2)/TWO
         PRP(1)=ELP
         PRP(2)=ELFN
         PRP(3)=RHO
         PRP(4)=DRN
         PRP(5)=-DRN/RHO
         PRP(6)=WGT
         PRP(7)=SQRT(DG2)
C    FIX A MINOR ERROR (DMC 5/05), HAD BEEN MISSING *RR IN 6.0 TERM
         PRP(8)=-BB*(RR2/AA+TWO*RR/AA**2+TWO/AA**3+
     1           TWO*CC*(RR2/AA**3+6.0D+00*RR/AA**4+12.0D+00/AA**5))
         PRP(9)=XN*DDRHXX*XN+YN*DDRHYY*YN+ZN*DDRHZZ*ZN+
     1      (XN*DDRHXY*YN+XN*DDRHXZ*ZN+YN*DDRHYZ*ZN)*TWO
         PRP(10)=DDRHXX+DDRHYY+DDRHZZ
         DO K=1,10
            IF(PRP(K).LT.PRPMIN(K)) PRPMIN(K)=PRP(K)
            IF(PRP(K).GT.PRPMAX(K)) PRPMAX(K)=PRP(K)
         ENDDO
         DO K=1,7
            PRPTOT(K)=PRPTOT(K)+PRP(K)*WGT
         ENDDO
            PRPTOT(8)=PRPTOT(8)+PRP(8)*SPHW(I)
         DO K=9,10
            PRPTOT(K)=PRPTOT(K)+PRP(K)*WGT
         ENDDO
         WGTTOT=WGTTOT+WGT
         SPHTOT=SPHTOT+SPHW(I)
C
         IF(IPRT.LT.0) THEN
            WRITE(IW,8997) I,(PRP(K)*PRPSCL(K),K=1,10)
         ENDIF
C
C        NOTE THAT PRPSCL FACTORS ARE NOT INCLUDED IN PUNCH OUTPUT
C
         IF(IPUN.NE.0) THEN
            WRITE(IP,7997) I,(PRP(K),K=1,5)
            WRITE(IP,7997) I,(PRP(K),K=6,10)
         ENDIF
      ENDDO
C
      DO K=1,7
         PRPAVG(K)=PRPTOT(K)/WGTTOT
      ENDDO
         PRPAVG(8)=PRPTOT(8)/SPHTOT
      DO K=9,10
         PRPAVG(K)=PRPTOT(K)/WGTTOT
      ENDDO
      PRPTOT(6)=WGTTOT
      IF(IPUN.NE.0) THEN
         WRITE(IP,7998) ENDWRD
      ENDIF
C
      IF(IPRT.NE.0) THEN
         IF(IPRT.LT.0)  WRITE(IW,8998)
         WRITE(IW,8996) (PRPTOT(K)*PRPSCL(K),K=1,10)
         WRITE(IW,8995) (PRPMIN(K)*PRPSCL(K),K=1,10)
         WRITE(IW,8994) (PRPAVG(K)*PRPSCL(K),K=1,10)
         WRITE(IW,8993) (PRPMAX(K)*PRPSCL(K),K=1,10)
      ENDIF
C
C     ----- RESET FAST MEMORY -----
C
      CALL RETFM(NEED)
      CALL FLSHBF(IP)
      CALL FLSHBF(IW)
      RETURN
C
 8999 FORMAT(/,' --- ADDITIONAL SOLUTE DENSITY-RELATED PROPERTIES ON ',
     1         'CAVITY SURFACE ---')
 8998 FORMAT(10X,'POT',9X,'FLD',9X,'RHO',9X,'DRN',9X,'ZET/100',
     1        5X,'WGT/100',5X,'DRG',9X,'CHGPEN',6X,'DDN',9X,'DDL')
 8997 FORMAT(I5,10F12.7)
 8996 FORMAT(' TOT ',10F12.7)
 8995 FORMAT(' MIN ',10F12.7)
 8994 FORMAT(' AVG ',10F12.7)
 8993 FORMAT(' MAX ',10F12.7)
C
 7999 FORMAT(' --- SOLUTE DENSITY-RELATED PROPERTIES FOR -$SVPIDN-',
     1       ' --- ')
 7998 FORMAT(A8)
 7997 FORMAT(I5,5F15.10)
      END
C*MODULE SVPINP  *DECK SVPSET
      SUBROUTINE SVPSET
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- SVP INITIALIZATION FOR CAVITY SURFACE -----
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00, FOUR=4.0D+00)
      PARAMETER   (MXCAV=5810)
      LOGICAL     SOME
      LOGICAL     OUT
      LOGICAL     MORE
      LOGICAL     DBUG
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON/SVPFLT/RHOISO,TOLCAV,RSTART,RSPACE,RADSPH,RADSLV,TOLCHG,
     1              DIELST,DIELOP,ZIONP,ZIONM,TKELVN,ROTCAV(3),
     2              TRNCAV(3),CSDIAG,CVGLIN
      COMMON/SVPTGR/METHOD,IQORP,NSTAGE,IOPPRD,IRENRM,ISHAPE,ITRCAV,
     1       NDRCAV,ICVICK,IFRQUP,IRDRF,IONSTR,INTCAV,NPTLEB,NPTTHE,
     2       NPTPHI,IROTGR,ITRNGR,ICRD,LINEQ,IPNRF,IPRRF,IPNDN,IPRDN,
     3       IEX
      COMMON/SVPPAR/CAVC(3,MXCAV),CAVW(MXCAV),CAVQ(MXCAV),
     1              CAVN(3,MXCAV),SPHW(MXCAV),CAVP(MXCAV),
     2              CENCAV(3),RNTCAV(3),NCAV,IQP
      COMMON/SVPPRP/SURF,VOL,RHOEXT,CHGBAR,CHGRUL,CHGDEV
      COMMON/SVPDBG/IDBG
      DIMENSION ROT(3,3)
C
      DBUG=.FALSE.
      DBUG=DBUG.OR.IDBG.EQ.4
      MORE=.FALSE.
      MORE=MORE.OR.DBUG.OR.IDBG.EQ.3
      OUT =.FALSE.
      OUT =OUT.OR.MORE.OR.IDBG.EQ.2
      SOME=.FALSE.
      SOME=SOME.OR.OUT.OR.IDBG.EQ.1
C
      PI=FOUR*ATAN(ONE)
      FOURPI=FOUR*PI
C
C     --- OBTAIN INITIAL CAVITY SURFACE ---
C         AND REACTION FIELD INTERNALLY
C
C     FOR ISODENSITY CAVITY, GENERATE A SPHERICAL INITIAL GUESS
C     FOR SPHERICAL CAVITY, GENERATE THE FINAL SPHERICAL SURFACE
C     FOR SHAPED CAVITY, GENERATE THE FINAL SHAPED SURFACE
C
      IF(.NOT.(IRDRF.EQ.1.AND.IFRQUP.EQ.4)) THEN
C
         IF(ISHAPE.EQ.0.OR.ISHAPE.EQ.1) THEN
C
C     --- SET UP INTEGRATION POINTS AND WEIGHTS ON A UNIT SPHERE ---
C
            IF(INTCAV.EQ.0) CALL LEBPTW(CAVN,CAVC,SPHW,NPTLEB,MXCAV)
            IF(INTCAV.EQ.1) CALL SVPPTW(CAVC,SPHW,NPTTHE,NPTPHI)
C
C     --- DETERMINE TRANSLATION OF THE UNIT SPHERE ---
C
            CALL SVPTRN(ITRNGR,TRNCAV,CENCAV,OUTRAD)
C
C     --- DETERMINE ROTATION OF THE UNIT SPHERE ---
C
            CALL SVPROT(IROTGR,ROTCAV,RNTCAV,CENCAV,ROT)
C
C     --- DETERMINE FINAL RADIUS OF SPHERE ---
C
            IF(ISHAPE.EQ.0) RADIUS=OUTRAD+TWO
            IF(ISHAPE.EQ.1) RADIUS=RADSPH
            RAD2=RADIUS**2
C
C     --- TRANSLATE, ROTATE, SCALE UNIT SPHERE TO GET FINAL SPHERE ---
C
            DO I=1,NCAV
               CAVXI=CAVC(1,I)
               CAVYI=CAVC(2,I)
               CAVZI=CAVC(3,I)
               CAVN(1,I)=ROT(1,1)*CAVXI+ROT(1,2)*CAVYI+ROT(1,3)*CAVZI
               CAVN(2,I)=ROT(2,1)*CAVXI+ROT(2,2)*CAVYI+ROT(2,3)*CAVZI
               CAVN(3,I)=ROT(3,1)*CAVXI+ROT(3,2)*CAVYI+ROT(3,3)*CAVZI
               CAVC(1,I)=CAVN(1,I)*RADIUS+CENCAV(1)
               CAVC(2,I)=CAVN(2,I)*RADIUS+CENCAV(2)
               CAVC(3,I)=CAVN(3,I)*RADIUS+CENCAV(3)
               SPHW(I)=SPHW(I)*FOURPI
               CAVW(I)=SPHW(I)*RAD2
            ENDDO
            SURF=FOURPI*RAD2
            VOL=FOURPI*RADIUS*RAD2/THREE
C
         ELSEIF(ISHAPE.EQ.2) THEN
C
C     --- SHAPED CAVITY FROM INTERSECTING SPHERES ---
C
            CALL SVP_CAVMAK
C
         ELSEIF(ISHAPE.EQ.3) THEN
C
C     --- SHAPED CAVITY FROM -IR- ---
C
C
         ELSE
C
C     --- ILLEGAL VALUE OF -ISHAPE- HERE ... ---
C
            WRITE(IW,*) 'ILLEGAL -ISHAPE- IN -SVPSET- ... ISHAPE = ',
     1                                                    ISHAPE
            CALL ABRT
         ENDIF
      ENDIF
C
      RETURN
      END
C*MODULE SVPINP  *DECK SVPTRN
      SUBROUTINE SVPTRN(ITRNGR,TRNCAV,CENCAV,OUTRAD)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- LOCATE CENTER OF CAVITY SURFACE INTEGRATION GRID -----
C
      PARAMETER  (MXATM=2000)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00,  TWO=2.0D+00)
      PARAMETER (TOANGS=0.52917724924D+00,TOBOHR=ONE/TOANGS)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DIMENSION C1(3,MXATM),AI(3,3),TR(3)
      DIMENSION OUTMAX(3),OUTMIN(3)
      DIMENSION TRNCAV(3),CENCAV(3)
C
C     ---- DEFAULT TO GRID CENTER AT ORIGIN OF NUCLEAR COORDINATES ----
C
      DO K=1,3
         CENCAV(K) = ZERO
      ENDDO
C
C     ----- GRID CENTER AT CENTER OF NUCLEAR MASS -----
C
      IF(ITRNGR.EQ.1) THEN
         CALL MOLINR(C,C1,NAT,CENCAV,AI,TR)
C
C     ----- GRID CENTER AT CENTER OF NUCLEAR CHARGE -----
C
      ELSEIF(ITRNGR.EQ.2) THEN
         TCHRG=ZERO
         DO IAT=1,NAT
            TCHRG=TCHRG+ZAN(IAT)
            DO K=1,3
               CENCAV(K)=CENCAV(K)+ZAN(IAT)*C(K,IAT)
            ENDDO
          ENDDO
         DO K=1,3
            CENCAV(K)=CENCAV(K)/TCHRG
         ENDDO
C
C     ----- GRID CENTER AT MIDPOINT OF OUTERMOST ATOMS -----
C
      ELSEIF(ITRNGR.EQ.3) THEN
         DO K=1,3
            OUTMIN(K)= 1.0D+10
            OUTMAX(K)=-1.0D+10
            DO IAT=1,NAT
               IF(C(K,IAT).LT.OUTMIN(K)) OUTMIN(K)=C(K,IAT)
               IF(C(K,IAT).GT.OUTMAX(K)) OUTMAX(K)=C(K,IAT)
            ENDDO
            CENCAV(K)=(OUTMIN(K)+OUTMAX(K))/TWO
         ENDDO
C
C     ----- GRID CENTER AT MIDPOINT OF OUTERMOST NONHYDROGEN ATOMS -----
C
      ELSEIF(ITRNGR.EQ.4) THEN
         DO K=1,3
            OUTMIN(K)= 1.0D+10
            OUTMAX(K)=-1.0D+10
            DO IAT=1,NAT
               NUCZ = INT(ZAN(IAT) + IZCORE(IAT))
               IF(NUCZ.NE.1) THEN
                  IF(C(K,IAT).LT.OUTMIN(K)) OUTMIN(K)=C(K,IAT)
                  IF(C(K,IAT).GT.OUTMAX(K)) OUTMAX(K)=C(K,IAT)
               ENDIF
            ENDDO
            CENCAV(K)=(OUTMIN(K)+OUTMAX(K))/TWO
         ENDDO
C
C     ----- GRID CENTER AT USER SPECIFIED POSITION IN BOHR -----
C
      ELSEIF(ITRNGR.EQ.5) THEN
         DO K=1,3
            CENCAV(K)=TRNCAV(K)
         ENDDO
C
C     ----- GRID CENTER AT USER SPECIFIED POSITION IN ANGSTROM -----
C
      ELSEIF(ITRNGR.EQ.6) THEN
         DO K=1,3
            CENCAV(K)=TRNCAV(K)*TOBOHR
         ENDDO
C
      ENDIF
C
C     ----- DETERMINE DISTANCE FROM ORIGIN TO OUTERMOST ATOM -----
C
      OUTRAD=ZERO
      DO IAT=1,NAT
         DIST=SQRT((C(1,IAT)-CENCAV(1))**2+(C(2,IAT)-CENCAV(2))**2
     1            +(C(3,IAT)-CENCAV(3))**2)
         IF(DIST.GT.OUTRAD) OUTRAD=DIST
      ENDDO
C
      RETURN
      END
C*MODULE SVPINP  *DECK SVPROT
      SUBROUTINE SVPROT(IROTGR,ROTCAV,RNTCAV,CENCAV,ROT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     --- DETERMINE ROTATION OF AXES FOR CAVITY INTEGRATION POINTS ---
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, FOUR=4.0D+00)
      PARAMETER  (MXATM=2000)
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      DIMENSION ROTCAV(3),RNTCAV(3),CENCAV(3)
      DIMENSION C1(3,MXATM),TR(3),ROT(3,3)
      DIMENSION AI(3,3),EIG(3)
C
      PI=FOUR*ATAN(ONE)
      RADIAN=PI/180.0D+00
C
C     ----- DEFAULT TO NO ROTATION -----
C
      IF(IROTGR.EQ.0) THEN
         DO K=1,3
            ROT(K,K)=ONE
            DO L=1,3
               IF(L.NE.K) ROT(K,L)=ZERO
            ENDDO
         ENDDO
C
C     ----- FIND ROTATION TO PRINCIPAL AXES OF NUCLEAR INERTIA ---
C
      ELSEIF(IROTGR.EQ.1) THEN
         DO I=1,3
            DO J=1,3
              AI(I,J)=ZERO
            ENDDO
         ENDDO
         CALL MOLINR(C,C1,NAT,CENCAV,AI,TR)
         CALL MOLAXS(AI,ROT,EIG,3,3,3)
C
C     ----- FIND ROTATION TO PRINCIPAL AXES OF NUCLEAR CHARGE ---
C
      ELSEIF(IROTGR.EQ.2) THEN
         DO I=1,3
            DO J=1,3
               AI(I,J)=ZERO
            ENDDO
         ENDDO
         DO IAT=1,NAT
            X =C(1,IAT)-CENCAV(1)
            Y =C(2,IAT)-CENCAV(2)
            Z =C(3,IAT)-CENCAV(3)
            CHG = ZAN(IAT) + IZCORE(IAT)
            AI(1,1)=AI(1,1)+CHG*(Y*Y+Z*Z)
            AI(2,1)=AI(2,1)-CHG* X*Y
            AI(1,2)=AI(2,1)
            AI(3,1)=AI(3,1)-CHG* X*Z
            AI(1,3)=AI(3,1)
            AI(2,2)=AI(2,2)+CHG*(X*X+Z*Z)
            AI(3,2)=AI(3,2)-CHG* Y*Z
            AI(2,3)=AI(3,2)
            AI(3,3)=AI(3,3)+CHG*(X*X+Y*Y)
         ENDDO
         CALL MOLAXS(AI,ROT,EIG,3,3,3)
C
C     ----- FIND ROTATION FOR USER SPECIFIED EULER ANGLES -----
C
      ELSEIF(IROTGR.EQ.3) THEN
         THE=ROTCAV(1)
         PHI=ROTCAV(2)
         CHI=ROTCAV(3)
         COSTHE=COS(THE*RADIAN)
         SINTHE=SIN(THE*RADIAN)
         COSPHI=COS(PHI*RADIAN)
         SINPHI=SIN(PHI*RADIAN)
         COSCHI=COS(CHI*RADIAN)
         SINCHI=SIN(CHI*RADIAN)
         ROT(1,1)= COSTHE*COSPHI*COSCHI-SINPHI*SINCHI
         ROT(1,2)= COSTHE*SINPHI*COSCHI+COSPHI*SINCHI
         ROT(1,3)=-SINTHE       *COSCHI
         ROT(2,1)=-COSTHE*COSPHI*SINCHI-SINPHI*COSCHI
         ROT(2,2)=-COSTHE*SINPHI*SINCHI+COSPHI*COSCHI
         ROT(2,3)= SINTHE       *SINCHI
         ROT(3,1)= SINTHE*COSPHI
         ROT(3,2)= SINTHE*SINPHI
         ROT(3,3)= COSTHE
C
      ENDIF
C
C     ----- EULER ANGLES OF CAVITY ROTATION -----
C
      IF(IROTGR.EQ.1.OR.IROTGR.EQ.2) THEN
         THE=ACOS(ROT(3,3))
         IF(ABS(ROT(3,3)-ONE).GT.1.0D-14) THEN
            PHI=ATAN2(ROT(3,2),ROT(3,1))
            CHI=ATAN2(ROT(2,3),-ROT(1,3))
C        SPECIAL TREATMENT FOR THETA=0
         ELSE
            PHI=ATAN2(ROT(1,2),ROT(1,1))
            CHI=ZERO
         ENDIF
         RNTCAV(1)=THE
         RNTCAV(2)=PHI
         RNTCAV(3)=CHI
      ENDIF
C
      RETURN
      END
C*MODULE SVPINP  *DECK SVP_CAVMAK
      SUBROUTINE SVP_CAVMAK
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C              ----- STARTING FROM -ICOSAHEDRON- -----
C
C     PASS, NAPEX, NFACE, ERROR =   0      12      20      20
C     PASS, NAPEX, NFACE, ERROR =   1      42      80     100    0.4982
C     PASS  NAPEX, NFACE, ERROR =   2     162     320     420    0.1848
C     PASS  NAPEX, NFACE, ERROR =   3     642    1280    1700    0.0523
C     PASS  NAPEX, NFACE, ERROR =   4    2562    5120    6820    0.0135
C     PASS  NAPEX, NFACE, ERROR =   5   10242   20480   27300    0.0034
C
C              ----- STARTING FROM -OCTAHEDRON-  -----
C
C     PASS, NAPEX, NFACE, ERROR =   0       6       8       8
C     PASS, NAPEX, NFACE, ERROR =   1      18      32      40    0.8075
C     PASS  NAPEX, NFACE, ERROR =   2      66     128     168    0.4557
C     PASS  NAPEX, NFACE, ERROR =   3     258     512     680    0.1619
C     PASS  NAPEX, NFACE, ERROR =   4    1026    2048    2728    0.0451
C     PASS  NAPEX, NFACE, ERROR =   5    4098    8192   10920    0.0116
C     PASS  NAPEX, NFACE, ERROR =   6   16386   32768   43688    0.0029
C
      PARAMETER   (MXFACE=10920)
      PARAMETER   (MXAPEX=4098)
      PARAMETER   (MXATM=2000)
      LOGICAL     DBUG
      LOGICAL     OUT
      LOGICAL     SOME
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON/CAVPAR/MINBEM,MAXBEM,IFICOS
      COMMON/CAVPRT/IPRT
      COMMON /FMCOM / X(1)
C
C     XYZATM(3,MXATM)
C       RATM(  MXATM)
C       NSPA(  MXATM)
C       NPPA(  MXATM)
C     IJKFAC(3,MXFACE)
C     XYZSEG(3,MXFACE)
C     IJKSEG(  MXFACE)
C     INSSEG(  MXFACE,MXATM)
C     XYZSPA(3,MXFACE,MXATM)
C     IJKSPA(  MXFACE,MXATM)
C     NUMPPS(  MXFACE,MXATM)
C       APEX(3,MXAPEX)
C
      DBUG=.FALSE.
      OUT =.FALSE.
      OUT =OUT.OR.DBUG.OR.IPRT.EQ.2
      SOME=.FALSE.
      SOME=SOME.OR.OUT.OR.IPRT.EQ.1
C
      IF(SOME) THEN
         WRITE(IW,9999)
         CALL FLSHBF(IW)
      ENDIF
C
      IF(IFICOS.EQ.0.AND.MAXBEM.GT.6) THEN
         WRITE(IW,*) '-MAXBEM- TOO LARGE FOR PARAMETERS IN -CAVMAK-'
         CALL FLSHBF(IW)
         CALL ABRT
      ELSEIF(IFICOS.NE.0.AND.MAXBEM.GT.7) THEN
         WRITE(IW,*) '-MAXBEM- TOO LARGE FOR PARAMETERS IN -CAVMAK-'
         CALL FLSHBF(IW)
         CALL ABRT
      ENDIF
C
      NFAC=MXFACE
      NAPX=MXAPEX
C
C     ----- PARTITION MEMORY -----
C
      CALL VALFM(LOADFM)
      I10 =LOADFM + 1
      I20 =I10 +3*NAT
      I30 =I20 +  NAT
      I40 =I30 +  NAT
      I50 =I40 +  NAT
      I60 =I50 +3*NFAC
      I70 =I60 +3*NFAC
      I80 =I70 +  NFAC
      I90 =I80 +  NFAC*NAT
      I100=I90 +3*NFAC*NAT
      I110=I100+  NFAC*NAT
      I120=I110+  NFAC*NAT
      LAST=I120+3*NAPX
      NEED = LAST - LOADFM - 1
C
      CALL GOTFM(NGOTMX)
      IF(NEED.GT.NGOTMX) THEN
         WRITE(IW,9010) NEED,NGOTMX
         CALL ABRT
      ENDIF
      CALL GETFM(NEED)
C
C     ----- GET -CAVITY- SURFACE -----
C
      CALL SVP_CAVSRF(NAT,C,NAT,NFAC,NAPX,
     1                X(I10),X(I20),X(I30),X(I40),
     2                X(I50),X(I60),X(I70),
     3                X(I80),X(I90),X(I100),X(I110),X(I120),
     4                SOME)
      CALL RETFM(NEED)
      RETURN
C
 9999 FORMAT(/,10X,16('-'),
     1       /,10X,'-CAVITY- SURFACE',
     2       /,10X,16('-'))
 9010 FORMAT(1X,'ERROR: NOT ENOUGH MEMORY IN -CAVMAK-. STOP.'/
     *       1X,'MEMORY NEEDED=',I10,' BUT ONLY HAVE',I10)
      END
C*MODULE SVPINP  *DECK SVP_CAVSRF
      SUBROUTINE SVP_CAVSRF(NAT,C,MXATMX,MXFAC,MXAPX,
     1                      XYZATM,RATM,NSPA,NPPA,
     2                      IJKFAC,XYZSEG,IJKSEG,
     3                      INSSEG,XYZSPA,IJKSPA,NUMPPS,APEX,
     4                      SOME)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (TOANGS=0.52917724924D+00,TOBOHR=ONE/TOANGS)
      PARAMETER (MXATM=2000)
      LOGICAL    SOME
      LOGICAL    OUT
      LOGICAL    MORE
      LOGICAL    DBUG
      LOGICAL    INSSEG
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON/CAVDAT/ADIAG,DSURF,DVOL,SRFMOL,VOLMOL
      COMMON/CAVRAD/RADIUS(MXATM)
      COMMON/CAVSLV/RSOLV
      COMMON/CAVLAY/RLAYER,NLAYER
      DIMENSION      C(3,*    )
      DIMENSION XYZATM(3,MXATMX)
      DIMENSION   RATM(  MXATMX)
      DIMENSION   NSPA(  MXATMX)
      DIMENSION   NPPA(  MXATMX)
      DIMENSION IJKFAC(3,MXFAC)
      DIMENSION XYZSEG(3,MXFAC)
      DIMENSION IJKSEG(  MXFAC)
      DIMENSION INSSEG(  MXFAC,MXATMX)
      DIMENSION XYZSPA(3,MXFAC,MXATMX)
      DIMENSION IJKSPA(  MXFAC,MXATMX)
      DIMENSION NUMPPS(  MXFAC,MXATMX)
      DIMENSION   APEX(3,MXAPX)
C
      DBUG=.FALSE.
      MORE=.FALSE.
      MORE=MORE.OR.DBUG
      OUT =.FALSE.
      OUT =OUT.OR.MORE
      SOME=SOME.OR.OUT
C
C     ----- APPROXIMATE SPHERE WITH SEGMENTS AND POINTS -----
C
      CALL SVP_CAVSPH(NSEG,NFAC,NDIV,
     1                IJKFAC,XYZSEG,IJKSEG,MXFAC,APEX,MXAPX,
     2                DSURF,DVOL,ADIAG)
C
C     ----- DEBUG PRINTING -----
C
      IF(OUT) THEN
         WRITE(IW,9999) NSEG,NFAC,NDIV,DSURF,DVOL
         WRITE(IW,9995) ADIAG
         IF(MORE) THEN
            WRITE(IW,9998)
            DO LSEG=1,NSEG
               WRITE(IW,9997) LSEG,XYZSEG(1,LSEG),XYZSEG(2,LSEG),
     1                             XYZSEG(3,LSEG),IJKSEG(  LSEG)
            ENDDO
         ENDIF
         IF(DBUG) THEN
            WRITE(IW,9996)
            DO LFAC=1,NFAC
               MFAC=LFAC+NSEG
               WRITE(IW,9997) MFAC,XYZSEG(1,MFAC),XYZSEG(2,MFAC),
     1                             XYZSEG(3,MFAC),IJKSEG(  MFAC)
            ENDDO
         ENDIF
         CALL FLSHBF(IW)
      ENDIF
C
C     ----- LOOP OVER LAYERS ... -----
C
      IPTLAY=0
      DO ILAYER=1,NLAYER
C
C     ----- SET MOLECULE -----
C
         DO IAT=1,NAT
            DO I=1,3
               XYZATM(I,IAT)=C(I,IAT)
            ENDDO
         ENDDO
         DO IAT=1,NAT
            IF(RADIUS(IAT).EQ.ZERO) THEN
               RATM(IAT)=ZERO
            ELSE
               RATM(IAT)=(RADIUS(IAT)+RSOLV+RLAYER*(ILAYER-1))*TOBOHR
            ENDIF
         ENDDO
C
C     ----- CREATE -SOLVENT ACCESSIBLE SURFACE- OF THE MOLECULE -----
C
         CALL SVP_CAVSAS(NAT,XYZATM,RATM,MXATMX,
     1                   NSPA,NPPA,XYZSPA,IJKSPA,
     2                   NSEG,NFAC,XYZSEG,IJKSEG,INSSEG,
     3                   NUMPPS,MXFAC,
     4                   ILAYER,IPTLAY,
     5                   SOME)
C
      ENDDO
C
      IF(SOME) THEN
         WRITE(IW,9994)
         CALL FLSHBF(IW)
      ENDIF
C
      RETURN
 9999 FORMAT(' NSEG,NFAC,NDIV=NFAC/NSEG,DSURF,DVOL = ',3I7,2F10.6)
 9998 FORMAT('  SEG  ','      X     ','      Y     ','      Z     ',
     1       ' SEG ',/,1X,47('-'))
 9997 FORMAT(I7,3F12.8,I5,F12.8)
 9996 FORMAT('  FAC  ','      X     ','      Y     ','      Z     ',
     1       ' SEG ',/,1X,47('-'))
 9995 FORMAT(' ADIAG           = ',F12.6)
 9994 FORMAT(/)
      END
C*MODULE SVPINP  *DECK SVP_CAVSAS
      SUBROUTINE SVP_CAVSAS(NAT,XYZATM,RATM,MXATM,
     1                      NSPA,NPPA,XYZSPA,IJKSPA,
     2                      NSEG,NFAC,XYZSEG,IJKSEG,INSSEG,
     3                      NUMPPS,MXFACE,
     4                      ILAYER,IPTLAY,
     5                      SOME)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER   (MXCAV=5810)
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (TOANGS=0.52917724924D+00)
      LOGICAL     LAYER
      LOGICAL     SOME
      LOGICAL     OUT
      LOGICAL     DBUG
      LOGICAL     INSSEG
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON/CAVPRT/IPRT
      COMMON/CAVSLV/RSOLV
      COMMON/CAVDAT/ADIAG,DSURF,DVOL,SRFMOL,VOLMOL
      COMMON/CAVNRG/CAVDSP
      COMMON/CAVBEM/CAVS(MXCAV),IATCAV(MXCAV)
      COMMON/SVPPAR/CAVC(3,MXCAV),CAVW(MXCAV),CAVQ(MXCAV),
     1              CAVN(3,MXCAV),SPHW(MXCAV),CAVP(MXCAV),
     2              CENCAV(3),RNTCAV(3),NCAV,IQP
      COMMON/SVPPRP/SURF,VOL,RHOEXT,CHGBAR,CHGRUL,CHGDEV
      DIMENSION XYZATM(3,MXATM)
      DIMENSION   RATM(  MXATM)
      DIMENSION   NSPA(  MXATM)
      DIMENSION   NPPA(  MXATM)
      DIMENSION XYZSEG(3,MXFACE)
      DIMENSION IJKSEG(  MXFACE)
      DIMENSION INSSEG(  MXFACE,MXATM)
      DIMENSION XYZSPA(3,MXFACE,MXATM)
      DIMENSION IJKSPA(  MXFACE,MXATM)
      DIMENSION NUMPPS(  MXFACE,MXATM)
C
      DIST(XI,YI,ZI,XJ,YJ,ZJ)=SQRT((XJ-XI)**2+(YJ-YI)**2+(ZJ-ZI)**2)
C
      DBUG=.FALSE.
      OUT =.FALSE.
      OUT =OUT.OR.DBUG
      OUT =OUT.OR.IPRT.GT.1
      SOME=SOME.OR.OUT
      SOME=SOME.OR.IPRT.EQ.1
C
      LAYER=.FALSE.
      LAYER=LAYER.OR.ILAYER.EQ.1
C
      IF(SOME.AND.LAYER) THEN
         WRITE(IW,9999)
      ENDIF
C
C     ----- PRINT ATOMIC CENTERS -----
C
      IF(OUT) THEN
         WRITE(IW,9998)
         DO IAT=1,NAT
            WRITE(IW,9997) IAT,XYZATM(1,IAT),XYZATM(2,IAT),
     1                                       XYZATM(3,IAT),
     2                    (RATM(IAT)*TOANGS-RSOLV)
         ENDDO
      ENDIF
C
C     ----- CLEAR ARRAYS ..... -----
C
      DO IAT=1,NAT
         DO I=1,MXFACE
            IJKSPA(I,IAT)=0
            NUMPPS(I,IAT)=0
         ENDDO
      ENDDO
C
C     ----- SIFT THROUGH ATOMIC CENTERS AND DECIDE IF A FACE -----
C           BELONGS TO THE -SAS- OR IS INSIDE THE MOLECULE.
C
      DO IAT=1,NAT
C
         NSPA(IAT)=0
         NPPA(IAT)=0
         IF(RATM(IAT).NE.ZERO) THEN
C
            DO ISEG=1,NSEG
               IJKSPA(       ISEG,IAT)=IJKSEG(       ISEG)
               DO M=1,3
                  XYZSPA(M,     ISEG,IAT)=XYZSEG(M,     ISEG)*RATM(IAT)
     1                                   +XYZATM(M,IAT)
               ENDDO
            ENDDO
            DO IFAC=1,NFAC
               IJKSPA(  IFAC+NSEG,IAT)=IJKSEG(  IFAC+NSEG)
               DO M=1,3
                  XYZSPA(M,IFAC+NSEG,IAT)=XYZSEG(M,IFAC+NSEG)*RATM(IAT)
     1                                   +XYZATM(M,IAT)
               ENDDO
            ENDDO
            IF(DBUG) THEN
               WRITE(IW,9996) IAT
               WRITE(IW,9995) (IJKSPA(IFAC+NSEG,IAT),IFAC=1,NFAC)
            ENDIF
            DO JAT=1,NAT
               DIJ=DIST(XYZATM(1,IAT),
     1                  XYZATM(2,IAT),
     2                  XYZATM(3,IAT),
     3                  XYZATM(1,JAT),
     4                  XYZATM(2,JAT),
     5                  XYZATM(3,JAT))
               IF(JAT.NE.IAT.AND.(DIJ.LT.(RATM(IAT)+RATM(JAT)))) THEN
                  DO IFAC=1,NFAC
                     DUM=DIST(XYZSPA(1,IFAC+NSEG,IAT),
     1                        XYZSPA(2,IFAC+NSEG,IAT),
     2                        XYZSPA(3,IFAC+NSEG,IAT),
     3                        XYZATM(1,JAT),
     4                        XYZATM(2,JAT),
     5                        XYZATM(3,JAT))
                     IF(DUM.LT.RATM(JAT)) THEN
                        IJKSPA(IFAC+NSEG,IAT)=0
                     ENDIF
                  ENDDO
               ENDIF
            ENDDO
            IF(DBUG) THEN
               WRITE(IW,9996) IAT
               WRITE(IW,9995) (IJKSPA(IFAC+NSEG,IAT),IFAC=1,NFAC)
            ENDIF
C
C     ----- CHECK SEGMENTS BELONGING TO -SAS- -----
C
            DO IFAC=1,NSEG+NFAC
               INSSEG(IFAC,IAT)=.TRUE.
            ENDDO
            DO IFAC=1,NFAC
               ISEG=IJKSPA(IFAC+NSEG,IAT)
               IF(ISEG.NE.0) THEN
                  INSSEG(IFAC+NSEG,IAT)=.FALSE.
                  INSSEG(     ISEG,IAT)=.FALSE.
               ENDIF
            ENDDO
            IF(DBUG) THEN
               WRITE(IW,9994) IAT
               WRITE(IW,9993) (INSSEG(IFAC,IAT),IFAC=1,NSEG+NFAC)
            ENDIF
            MSEG=0
            DO ISEG=1,NSEG
               IF(.NOT.INSSEG(     ISEG,IAT)) THEN
                  MSEG=MSEG+1
               ENDIF
            ENDDO
            MFAC=0
            DO IFAC=1,NFAC
               IF(.NOT.INSSEG(IFAC+NSEG,IAT)) THEN
                  MFAC=MFAC+1
               ENDIF
            ENDDO
            NSPA(IAT)=MSEG
            NPPA(IAT)=MFAC
C
C     ----- SURFACE AREA OF SEGMENTS -----
C
            DO ISEG=1,NSEG
               NUMPPS(ISEG,IAT)=0
            ENDDO
            DO IFAC=1,NFAC
               ISEG=IJKSPA(IFAC+NSEG,IAT)
               IF(ISEG.NE.0) THEN
                  NUMPPS(ISEG,IAT)=NUMPPS(ISEG,IAT)+1
               ENDIF
            ENDDO
C
         ENDIF
C
      ENDDO
C
      IF(SOME.AND.LAYER) THEN
         WRITE(IW,9985) NSEG,NFAC
         WRITE(IW,9992)
         DO IAT=1,NAT
            NPP=0
            DO ISEG=1,NSEG
               NPP=NPP+NUMPPS(ISEG,IAT)
            ENDDO
            WRITE(IW,9991) IAT,NSPA(IAT),NPPA(IAT),NPP
         ENDDO
      ENDIF
      IF(DBUG) THEN
         WRITE(IW,9987)
         DO IAT=1,NAT
            DO ISEG=1,NSEG
               WRITE(IW,9986) IAT,ISEG,NUMPPS(ISEG,IAT)
            ENDDO
         ENDDO
      ENDIF
C
C     ----- SAVE COORDINATES OF SURFACE POINTS -----
C           SAVE SEGMENT SURFACES
C           SAVE SEGMENT TO ATOM MAPPING
C
      SRFMOL=ZERO
      VOLMOL=ZERO
      IEF0  =IPTLAY
      IEF   =0
      DO IAT=1,NAT
         IF(RATM(IAT).NE.ZERO) THEN
            DO ISEG=1,NSEG
               IF(.NOT.INSSEG(ISEG,IAT)) THEN
                  IEF=IEF+1
                  IF(IEF+IEF0.GT.MXCAV) THEN
                     WRITE(IW,*) 'TOO MANY -CAVITY- SURFACE POINTS'
                     CALL ABRT
                  ENDIF
C
                  IPP=NUMPPS(ISEG,IAT)
                  CAVS(IEF+IEF0)=         IPP*DSURF*RATM(IAT)**2
                  SRFMOL        =SRFMOL + IPP*DSURF*RATM(IAT)**2
                  VOLMOL        =VOLMOL + IPP*DVOL *RATM(IAT)**3
                  IATCAV(IEF+IEF0)=IAT
C
C     ----- ... INTO THE -SVP- ARRAYS ... -----
C
                  CAVW(IEF+IEF0)=CAVS(IEF+IEF0)
                  DO I=1,3
                     CAVC(I,IEF+IEF0)=XYZSEG(I,ISEG)*RATM(IAT)
     1                               +XYZATM(I,IAT )
                  ENDDO
                  DO I=1,3
                     CAVN(I,IEF+IEF0)=XYZSEG(I,ISEG)
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
      ENDDO
C
      SURF  =SRFMOL
      VOL   =VOLMOL
C
      SRFMOL=SRFMOL*(TOANGS**2)
      VOLMOL=VOLMOL*(TOANGS**3)
C
C
      IPTLAY=IEF+IEF0
      NCAV  =IEF+IEF0
      IF(SOME) THEN
         IF(ILAYER.EQ.1) THEN
            WRITE(IW,9990) NCAV
            WRITE(IW,9984) SRFMOL
            WRITE(IW,9983) VOLMOL
         ELSE
            WRITE(IW,9979) ILAYER,IEF,NCAV
         ENDIF
      ENDIF
C
C     ----- CAVITY/DISPERSION FREE ENERGY ---
C           SITKOFF, SHARP, AND HONIG,
C           J.PHYS.CHEM. 98, 1978 (1994)
C
      CAVDSP=0.860D+00+0.005D+00*SRFMOL
      IF(SOME) THEN
         IF(ILAYER.EQ.1) THEN
            WRITE(IW,9981) CAVDSP
         ENDIF
      ENDIF
C
C     ----- PRINT SEGMENT SURFACES -----
C
      IF(DBUG) THEN
         WRITE(IW,9989)
         DO IEF=1,NCAV
            WRITE(IW,9988) IEF,CAVS(IEF),IATCAV(IEF),
     1                     (CAVC(I,IEF),I=1,3)
         ENDDO
      ENDIF
      IF(OUT) THEN
         WRITE(IW,9982)
         DO IEF=1,NCAV
            WRITE(IW,9988) IEF,CAVS(IEF),IATCAV(IEF),
     1                     (CAVC(I,IEF),I=1,3)
         ENDDO
      ENDIF
C
      IF(OUT) THEN
         WRITE(IW,9980)
      ENDIF
      CALL FLSHBF(IW)
C
      RETURN
 9999 FORMAT(/,1X,'SOLVENT ACCESSIBLE SURFACE',/,1X,26('-'))
 9998 FORMAT(/,1X,'---------- ATOMIC COORDINATES (A.U.) ----------',
     1            '-- VDWR(ANG.) --')
 9997 FORMAT(  1X,I5,3F14.8,F10.3)
 9996 FORMAT(/,1X,'---------- SEGMENTS FOR -IAT- = ',I5)
 9995 FORMAT(16I4)
 9994 FORMAT(/,1X,'-INSSEG- FACES FOR IAT = ',I5)
 9993 FORMAT(16L4)
 9992 FORMAT(  1X,'ATOM',' ( ','  NSPA',',','  NPPA',' )',/,1X,22('-'))
 9991 FORMAT(  1X,   I4 ,' ( ',     I6 ,',',     I6 ,' )',I8)
 9990 FORMAT(  1X,'NUMBER OF -CAVITY- SURFACE POINTS = ',I8)
 9989 FORMAT(  1X,'SEGMENT SURFACES =',/,1X,18('-'))
 9988 FORMAT(I8,F10.5,I5,3F12.6)
 9987 FORMAT(  1X,'NUMBER OF FACES / SEGMENT =',/,1X,27('-'))
 9986 FORMAT(3I5)
 9985 FORMAT(' NUMBER OF SEGMENTS PER ATOM = ',I10,/,
     1       ' NUMBER OF   POINTS PER ATOM = ',I10)
 9984 FORMAT(' MOLECULAR SURFACE = ',F10.3,' ANGSTROM**2')
 9983 FORMAT(' MOLECULAR VOLUME  = ',F10.3,' ANGSTROM**3')
 9982 FORMAT(  1X,'-CAVITY- POINT CHARGES =',/,1X,23('-'))
 9981 FORMAT(' G(CAV/DISP)       = ',F10.3,' KCAL/MOL')
 9980 FORMAT(/)
 9979 FORMAT(  1X,'LAYER = ',I3,' # OF LAYER POINTS = ',I8,
     1                          ' TOTAL # POINTS = ',I8)
      END
C*MODULE SVPINP  *DECK SVP_CAVSPH
      SUBROUTINE SVP_CAVSPH(NSEG,NFAC,NDIV,
     1                      IJKFAC,XYZSEG,IJKSEG,MXFAC,APEX,MXAPX,
     2                      DSURF,DVOL,ADIAG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C              ----- STARTING FROM -ICOSAHEDRON- -----
C
C     PASS, NAPEX, NFACE, ERROR =   0      12      20      20
C     PASS, NAPEX, NFACE, ERROR =   1      42      80     100    0.4982
C     PASS  NAPEX, NFACE, ERROR =   2     162     320     420    0.1848
C     PASS  NAPEX, NFACE, ERROR =   3     642    1280    1700    0.0523
C     PASS  NAPEX, NFACE, ERROR =   4    2562    5120    6820    0.0135
C     PASS  NAPEX, NFACE, ERROR =   5   10242   20480   27300    0.0034
C
C              ----- STARTING FROM -OCTAHEDRON-  -----
C
C     PASS, NAPEX, NFACE, ERROR =   0       6       8       8
C     PASS, NAPEX, NFACE, ERROR =   1      18      32      40    0.8075
C     PASS  NAPEX, NFACE, ERROR =   2      66     128     168    0.4557
C     PASS  NAPEX, NFACE, ERROR =   3     258     512     680    0.1619
C     PASS  NAPEX, NFACE, ERROR =   4    1026    2048    2728    0.0451
C     PASS  NAPEX, NFACE, ERROR =   5    4098    8192   10920    0.0116
C     PASS  NAPEX, NFACE, ERROR =   6   16386   32768   43688    0.0029
C
      PARAMETER (MXPASS=    7)
C
      LOGICAL ICOS
      LOGICAL OCTA
      LOGICAL SOME,OUT,DBUG
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON/CAVPAR/MINBEM,MAXBEM,IFICOS
C
      DIMENSION   APEX(3,MXAPX)
      DIMENSION IJKFAC(3,MXFAC)
      DIMENSION IJKSEG(  MXFAC)
      DIMENSION XYZSEG(3,MXFAC)
      DIMENSION MINFAC(MXPASS)
      DIMENSION MAXFAC(MXPASS)
      DIMENSION MINICO(MXPASS)
      DIMENSION MAXICO(MXPASS)
      DIMENSION MINOCT(MXPASS)
      DIMENSION MAXOCT(MXPASS)
      DIMENSION IJKNEW(3)
      DIMENSION IJKOLD(3)
C
      EQUIVALENCE (IJKOLD(1),IOLD),(IJKOLD(2),JOLD),(IJKOLD(3),KOLD)
      EQUIVALENCE (IJKNEW(1),INEW),(IJKNEW(2),JNEW),(IJKNEW(3),KNEW)
C
      DATA MINICO /    1,   21,  101,  421, 1701, 6821,    0/
      DATA MAXICO /   20,  100,  420, 1700, 6820,27300,    0/
      DATA MINOCT /    1,    9,   41,  169,  681, 2729,10921/
      DATA MAXOCT /    8,   40,  168,  680, 2728,10920,43688/
      DATA ZERO  /0.0D+00/
      DATA ONE   /1.0D+00/
      DATA THREE /3.0D+00/
      DATA FOUR  /4.0D+00/
C
      DIST(XI,YI,ZI,XJ,YJ,ZJ)=SQRT((XJ-XI)**2+(YJ-YI)**2+(ZJ-ZI)**2)
C
      DBUG=.FALSE.
      OUT =.FALSE.
      OUT =OUT.OR.DBUG
      SOME=.FALSE.
      SOME=SOME.OR.OUT
C
      PI=FOUR*ATAN(ONE)
      RAD=ONE
      SRF=FOUR*PI*RAD**2
      VOL=FOUR*PI*RAD**3/THREE
C
      NPASS=MAXBEM
C
C     ----- DEFINE  HEDRON  -----
C           DEFINE -MINFAC-
C           DEFINE -MAXFAC-
C
      ICOS=IFICOS.NE.0
      OCTA=.NOT.ICOS
      IF(ICOS) THEN
         CALL CAV_SPHICO(APEX,NAPEX,IJKFAC,IJKSEG,NFACE)
         DO IPASS=1,MXPASS
            MINFAC(IPASS)=MINICO(IPASS)
            MAXFAC(IPASS)=MAXICO(IPASS)
         ENDDO
      ENDIF
      IF(OCTA) THEN
         CALL CAV_SPHOCT(APEX,NAPEX,IJKFAC,IJKSEG,NFACE)
         DO IPASS=1,MXPASS
            MINFAC(IPASS)=MINOCT(IPASS)
            MAXFAC(IPASS)=MAXOCT(IPASS)
         ENDDO
      ENDIF
      IF(SOME.OR.OUT.OR.DBUG) THEN
         IF(ICOS) THEN
            WRITE(IW,9994)
            CALL FLSHBF(IW)
         ENDIF
         IF(OCTA) THEN
            WRITE(IW,9982)
            CALL FLSHBF(IW)
         ENDIF
         IF(OUT) THEN
            WRITE(IW,*) '-MINBEM- = ',MINBEM
            WRITE(IW,*) '-MAXBEM- = ',MAXBEM
            WRITE(IW,*) '-MINFAC- = ',MINFAC
            WRITE(IW,*) '-MAXFAC- = ',MAXFAC
            WRITE(IW,*) '-NPASS - = ',NPASS
            WRITE(IW,9999)
            DO IAPEX=1,NAPEX
               WRITE(IW,9998) IAPEX,APEX(1,IAPEX),
     1                              APEX(2,IAPEX),
     2                              APEX(3,IAPEX)
            ENDDO
            CALL FLSHBF(IW)
         ENDIF
      ENDIF
C
C     ----- LOOP OVER DIVISIONS TO CREATE SPHERE -----
C
      MXFAC=0
      IPASS=1
  100 IPASS=IPASS+1
         MNFAC=MXFAC+1
         MXFAC=NFACE
         IF(OUT) THEN
            WRITE(IW,9996) IPASS,NAPEX,NFACE,MNFAC,MXFAC
            CALL FLSHBF(IW)
         ENDIF
C
         DMIN =ONE
         MAPEX=NAPEX
         MFACE=NFACE
         DO LFACE=MNFAC,MXFAC
            IOLD=IJKFAC(1,LFACE)
            JOLD=IJKFAC(2,LFACE)
            KOLD=IJKFAC(3,LFACE)
            CALL CAV_SPHAPX(APEX,MAPEX,IJKFAC,IJKSEG,MFACE,LFACE,
     1                  IJKOLD,IJKNEW,DIJK)
            DMIN=MIN(DMIN,DIJK)
         ENDDO
         NAPEX=MAPEX
         NFACE=MFACE
         IF(OUT) THEN
            WRITE(IW,9995) NAPEX,NFACE
            CALL FLSHBF(IW)
         ENDIF
C
C     ----- PRINT OUT NEW APECES -----
C
         IF(DBUG) THEN
            DO IAPEX=1,NAPEX
               WRITE(IW,9998) IAPEX,APEX(1,IAPEX),APEX(2,IAPEX),
     1                              APEX(3,IAPEX)
            ENDDO
            CALL FLSHBF(IW)
         ENDIF
C
C     ----- PRINT APPROXIMATE VOLUME -----
C
         RADAPP=    DMIN
         RADRAT=    DMIN
         RADERR=ONE-RADRAT
         SRFAPP=SRF*DMIN**2
         SRFRAT=    DMIN**2
         SRFERR=ONE-SRFRAT
         VOLAPP=VOL*DMIN**3
         VOLRAT=    DMIN**3
         VOLERR=ONE-VOLRAT
         IF(OUT) THEN
            WRITE(IW,9997) VOL,VOLAPP,VOLRAT,VOLERR
            WRITE(IW,9992) SRF,SRFAPP,SRFRAT,SRFERR
            WRITE(IW,9991) RAD,RADAPP,RADRAT,RADERR
            CALL FLSHBF(IW)
         ENDIF
C
C     ----- ASSIGN EARLY SEGMENT TO EACH FACE -----
C
         IF(IPASS.GT.(MINBEM+1)) THEN
            IF(DBUG) THEN
               WRITE(IW,9981) IPASS
               WRITE(IW,9980) (MINFAC(I),I=1,IPASS)
               WRITE(IW,9979) (MAXFAC(I),I=1,IPASS)
               CALL FLSHBF(IW)
            ENDIF
            MAXSEG=MAXFAC(MINBEM)
            LFACMN=MINFAC(IPASS)
            LFACMX=MAXFAC(IPASS)
            IF(DBUG) THEN
               WRITE(IW,9990) IPASS
               WRITE(IW,9988) (IJKSEG(LFACE),LFACE=LFACMN,LFACMX)
               CALL FLSHBF(IW)
            ENDIF
            DO LFACE=LFACMN,LFACMX
               IJKSEG(LFACE)=IJKSEG(IJKSEG(LFACE))
               IF(IJKSEG(LFACE).GT.MAXSEG) THEN
                  WRITE(IW,9987) LFACE,IJKSEG(LFACE)
               ENDIF
            ENDDO
            IF(DBUG) THEN
               WRITE(IW,9989) IPASS
               WRITE(IW,9988) (IJKSEG(LFACE),LFACE=LFACMN,LFACMX)
               CALL FLSHBF(IW)
            ENDIF
         ENDIF
C
      IF(IPASS.LT.NPASS) GO TO 100
C
C     ----- END OF LOOP OVER TESSALATING PASSES -----
C
      IF(DBUG) THEN
         DO IPASS=1,NPASS
            LFACMN=MINFAC(IPASS)
            LFACMX=MAXFAC(IPASS)
            WRITE(IW,9989) IPASS
            WRITE(IW,*) '-LFACMN- = ',LFACMN
            WRITE(IW,*) '-LFACMX- = ',LFACMX
            WRITE(IW,9988) (IJKSEG(LFACE),LFACE=LFACMN,LFACMX)
            CALL FLSHBF(IW)
         ENDDO
      ENDIF
      IF(SOME.OR.OUT.OR.DBUG) THEN
         WRITE(IW,9993) NPASS,NAPEX,MINFAC(NPASS),MAXFAC(NPASS),
     1                  RADAPP,RADERR,SRFAPP,SRFERR,VOLAPP,VOLERR
         CALL FLSHBF(IW)
      ENDIF
C
C     ----- AT THIS POINT EACH OF THE FACES IS ASSIGNED TO ONE -----
C           SEGMENT. NOW DEFINE CENTERS OF SEGMENTS ...
C
      THIRD =ONE/THREE
      LFACMN= MINFAC(MINBEM)
      LFACMX= MAXFAC(MINBEM)
      DO LFACE=LFACMN,LFACMX
         MFACE=LFACE-LFACMN+1
         IJKSEG(MFACE)=MFACE
         I=IJKFAC(1,LFACE)
         J=IJKFAC(2,LFACE)
         K=IJKFAC(3,LFACE)
         DO M=1,3
            XYZSEG(M,MFACE)=(APEX(M,I)+APEX(M,J)+APEX(M,K))*THIRD
         ENDDO
         DSEG=ONE/DIST(XYZSEG(1,MFACE),XYZSEG(2,MFACE),XYZSEG(3,MFACE),
     1                 ZERO,ZERO,ZERO)
         DO M=1,3
            XYZSEG(M,MFACE)=XYZSEG(M,MFACE)*DSEG
         ENDDO
      ENDDO
      NSEG=(LFACMX-LFACMN+1)
C
      IF(DBUG) THEN
         LFACMN=1
         LFACMX=NSEG
         WRITE(IW,*)    'SEGMENT TO SEGMENT MAPPING = '
         WRITE(IW,9988) (IJKSEG(LFACE),LFACE=LFACMN,LFACMX)
      ENDIF
C
C     ----- NOW THE FACES ... -----
C
      IF(NPASS.GT.MINBEM) THEN
         LFACMN=MINFAC(MINBEM+1)
         LFACMX=MAXFAC(NPASS   )
         DO LFACE=LFACMN,LFACMX
            MFACE=LFACE-LFACMN+1
     1                        +(MAXFAC(MINBEM)-MINFAC(MINBEM)+1)
            IJKSEG(MFACE)=IJKSEG(LFACE)
     1                        -(               MINFAC(MINBEM)-1)
            I=IJKFAC(1,LFACE)
            J=IJKFAC(2,LFACE)
            K=IJKFAC(3,LFACE)
            DO M=1,3
               XYZSEG(M,MFACE)=(APEX(M,I)+APEX(M,J)+APEX(M,K))*THIRD
            ENDDO
            DSEG=ONE/DIST(XYZSEG(1,MFACE),
     1                    XYZSEG(2,MFACE),
     2                    XYZSEG(3,MFACE),ZERO,ZERO,ZERO)
            DO M=1,3
               XYZSEG(M,MFACE)=XYZSEG(M,MFACE)*DSEG
            ENDDO
         ENDDO
         NFAC=(LFACMX-LFACMN+1)
      ELSE
         DO ISEG=1,NSEG
            IFAC=ISEG+NSEG
            IJKSEG(IFAC)=IJKSEG(ISEG)
            DO M=1,3
               XYZSEG(M,IFAC)=XYZSEG(M,ISEG)
            ENDDO
         ENDDO
         NFAC=NSEG
      ENDIF
C
      IF(DBUG) THEN
         LFACMN=NSEG+1
         LFACMX=NSEG+NFAC
         WRITE(IW,*)    ' FACE   TO SEGMENT MAPPING = '
         WRITE(IW,9988) (IJKSEG(LFACE),LFACE=LFACMN,LFACMX)
      ENDIF
C
C     --- CALCULATE -DSURF DVOL- FOR THE CAVITY/DISPERSION ENERGY TERM
C
      NDIV =NFAC/NSEG
      DSURF=SRF/NFAC
      DVOL =VOL/NFAC
      IF(SOME.OR.OUT.OR.DBUG) THEN
         WRITE(IW,9986) NSEG,NFAC,NDIV,DSURF,DVOL
         CALL FLSHBF(IW)
      ENDIF
      IF(OUT) THEN
         WRITE(IW,9985)
         DO I=1,NSEG
            DONE=DIST(XYZSEG(1,I),XYZSEG(2,I),XYZSEG(3,I),
     1                ZERO,ZERO,ZERO)
            WRITE(IW,9984) I,
     1                     XYZSEG(1,I),XYZSEG(2,I),XYZSEG(3,I),
     2                     IJKSEG(I),DONE
         ENDDO
         CALL FLSHBF(IW)
      ENDIF
      IF(DBUG) THEN
         WRITE(IW,9985)
         DO I=NSEG+1,NSEG+NFAC
            DONE=DIST(XYZSEG(1,I),XYZSEG(2,I),XYZSEG(3,I),
     1                ZERO,ZERO,ZERO)
            WRITE(IW,9984) (I-NSEG),
     1                     XYZSEG(1,I),XYZSEG(2,I),XYZSEG(3,I),
     2                     IJKSEG(I),DONE
         ENDDO
      ENDIF
C
C     ----- CALCULATE -ADIAG- OF THE -SOLVATION- THEORY -----
C
      AVGDIA=ZERO
      AVGFAC=ZERO
      DO MSEG=1,NSEG
         SUM=ZERO
         DO LSEG=1,NSEG
            IF(LSEG.NE.MSEG) THEN
               L1=MSEG
               L2=LSEG
         SUM=SUM+RAD/DIST(XYZSEG(1,L2),XYZSEG(2,L2),XYZSEG(3,L2),
     1                    XYZSEG(1,L1),XYZSEG(2,L1),XYZSEG(3,L1))
            ENDIF
         ENDDO
         DMY = NSEG
         FAC=(DMY-SUM)/SQRT(DMY)
         ADIAG=SQRT(FOUR*PI)*FAC
         IF(SOME.OR.OUT.OR.DBUG) THEN
            WRITE(IW,9983) MSEG,ADIAG,FAC,DMY,SUM
            CALL FLSHBF(IW)
         ENDIF
         AVGDIA=AVGDIA+ADIAG
         AVGFAC=AVGFAC+FAC
      ENDDO
      ADIAG=AVGDIA/DMY
      FAC  =AVGFAC/DMY
      IF(SOME.OR.OUT.OR.DBUG) THEN
         WRITE(IW,9978)      ADIAG,FAC
         CALL FLSHBF(IW)
      ENDIF
C--
C--      ADIAG=SQRT(FOUR*PI)*1.07D+00
C--
C
      RETURN
 9999 FORMAT(/,1X,'  APEX',5X,'X',6X,5X,'Y',6X,5X,'Z',6X,/,1X,42('-'))
 9998 FORMAT(1X,I6,F12.8,F12.8,F12.8)
 9997 FORMAT(' VOL, APPROX., RATIO, ERROR = ',2F12.8,2 F8.4)
 9996 FORMAT(' PASS, NAPEX, NFACE, MNFAC, MXFAC = ',I3,4I8)
 9995 FORMAT('       NAPEX, NFACE               = ',3X,2I8)
 9994 FORMAT(1X,'SPHERE FROM -ICOSAHEDRON-',/,1X,25('-'))
 9993 FORMAT(' NPASS = ',I2,' NAPEX = ',I8,
     1       ' MINFAC = ',I8,' MAXFAC = ',I8,/,
     2       ' RAD = ',F10.6,' ERROR = ',F8.4,/,
     3       ' SRF = ',F10.6,' ERROR = ',F8.4,/,
     4       ' VOL = ',F10.6,' ERROR = ',F8.4)
 9992 FORMAT(' SRF, APPROX., RATIO, ERROR = ',2F12.8,2 F8.4)
 9991 FORMAT(' RAD, APPROX., RATIO, ERROR = ',2F12.8,2 F8.4)
 9990 FORMAT(' ABSOLUTE -IJKSEG- , FOR -IPASS- = ',I5)
 9989 FORMAT(' RELATIVE -IJKSEG- , FOR -IPASS- = ',I5)
 9988 FORMAT(12I6)
 9987 FORMAT(' ASSIGNED SEGMENT FOR -LFACE- = ',I7,
     1       ' IS = ',I7,' ( GREATER THAN -MAXSEG- = ',I4,' )')
 9986 FORMAT(' NSEG,NFAC,NDIV=NFAC/NSEG,DSURF,DVOL = ',3I7,2F10.6)
 9985 FORMAT('   PT  ','      X     ','      Y     ','      Z     ',
     1       ' SEG ','    NORM    ',/,1X,59('-'))
 9984 FORMAT(I7,3F12.8,I5,F12.8)
 9983 FORMAT(' MSEG,ADIAG,FAC,M,SUM = ',I7,4F12.6)
 9982 FORMAT(1X,'SPHERE FROM -OCTAHEDRON-',/,1X,24('-'))
 9981 FORMAT(' PASS # = ',I5)
 9980 FORMAT(' MINFAC = ',10I5)
 9979 FORMAT(' MAXFAC = ',10I5)
 9978 FORMAT('      ADIAG,FAC       = ',   2F12.6)
      END
C*MODULE SVPINP  *DECK CAV_SPHOCT
      SUBROUTINE CAV_SPHOCT(APEX,NAPEX,IJKFAC,IJKSEG,NFACE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL OUT
      COMMON /IOFILE/ IR,IW,IP,IJKA,IJKT,IDAF,NAV,IODA(950)
      DIMENSION    XYZ(3,6)
      DIMENSION    IJK(3,8)
      DIMENSION   APEX(3,*)
      DIMENSION IJKFAC(3,*)
      DIMENSION IJKSEG(  *)
      DATA XYZ / 1.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 1.0D+00, 0.0D+00,
     1          -1.0D+00, 0.0D+00, 0.0D+00, 0.0D+00,-1.0D+00, 0.0D+00,
     2           0.0D+00, 0.0D+00, 1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00/
      DATA IJK / 5, 1, 2, 5, 2, 3, 5, 3, 4, 5, 4, 1,
     1           6, 1, 2, 6, 2, 3, 6, 3, 4, 6, 4, 1/
C
      OUT=.FALSE.
      IF(OUT) THEN
         WRITE(IW,9997)
         CALL FLSHBF(IW)
      ENDIF
C
C     ----- SET THE 6 APECES OF AN OCTAHEDRON -----
C
C     1     1.     0.     0.
C     2     0.     1.     0.
C     3    -1.     0.     0.
C     4     0.    -1.     0.
C     5     0.     0.     1.
C     6     0.     0.    -1.
C
      NAPEX=6
      DO IAPEX=1,NAPEX
         APEX(1,IAPEX)=XYZ(1,IAPEX)
         APEX(2,IAPEX)=XYZ(2,IAPEX)
         APEX(3,IAPEX)=XYZ(3,IAPEX)
      ENDDO
      IF(OUT) THEN
         WRITE(IW,9999)
         DO IAPEX=1,NAPEX
            WRITE(IW,9998) IAPEX,APEX(1,IAPEX),APEX(2,IAPEX),
     1                           APEX(3,IAPEX)
         ENDDO
         CALL FLSHBF(IW)
      ENDIF
C
      NFACE=8
      DO LFACE=1,NFACE
         IJKFAC(1,LFACE)=IJK(1,LFACE)
         IJKFAC(2,LFACE)=IJK(2,LFACE)
         IJKFAC(3,LFACE)=IJK(3,LFACE)
         IJKSEG(  LFACE)=      LFACE
      ENDDO
C
      IF(OUT) THEN
         WRITE(IW,*) '...... END OF -CAV_SPHOCT- ......'
      ENDIF
      RETURN
 9999 FORMAT(/,1X,'  APEX',5X,'X',6X,5X,'Y',6X,5X,'Z',6X,/,1X,42('-'))
 9998 FORMAT(1X,I6,F12.8,F12.8,F12.8)
 9997 FORMAT(/,1X,'OCTAHEDRON',/,1X,10('-'))
      END
C*MODULE SVPINP  *DECK CAV_SPHICO
      SUBROUTINE CAV_SPHICO(APEX,NAPEX,IJKFAC,IJKSEG,NFACE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL OUT
      COMMON /IOFILE/ IR,IW,IP,IJKA,IJKT,IDAF,NAV,IODA(950)
      DIMENSION      C(3,12)
      DIMENSION      S(3,12)
      DIMENSION    IJK(3,20)
      DIMENSION   APEX(3,*)
      DIMENSION IJKFAC(3,*)
      DIMENSION IJKSEG(  *)
      DATA C   / 0.0D+00, 1.0D+00, 0.0D+00, 0.0D+00, 1.0D+00, 0.0D+00,
     1           0.0D+00,-1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00, 0.0D+00,
     2           0.0D+00, 0.0D+00, 1.0D+00, 0.0D+00, 0.0D+00, 1.0D+00,
     3           0.0D+00, 0.0D+00,-1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00,
     4           1.0D+00, 0.0D+00, 0.0D+00, 1.0D+00, 0.0D+00, 0.0D+00,
     5          -1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00, 0.0D+00, 0.0D+00/
      DATA S   / 0.0D+00, 0.0D+00, 1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00,
     1           0.0D+00, 0.0D+00, 1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00,
     2           1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00, 0.0D+00, 0.0D+00,
     3           1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00, 0.0D+00, 0.0D+00,
     4           0.0D+00, 1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00, 0.0D+00,
     5           0.0D+00, 1.0D+00, 0.0D+00, 0.0D+00,-1.0D+00, 0.0D+00/
      DATA IJK / 1, 2, 9, 1, 9, 5, 1, 5, 6, 1, 6,11, 1,11, 2,
     1                    2, 9, 7, 2, 7, 8, 2, 8,11,
     2           3, 4,10, 3,10, 5, 3, 5, 6, 3, 6,12, 3,12, 4,
     3                    4,10, 7, 4, 7, 8, 4, 8,12,
     4           9,10, 7, 9, 5,10,11,12, 8,11, 6,12/
      DATA ONE  /1.0D+00/
      DATA TWO  /2.0D+00/
      DATA FIVE /5.0D+00/
C
      OUT=.FALSE.
      IF(OUT) THEN
         WRITE(IW,9997)
      ENDIF
C
C     ----- SET THE 12 APECES OF AN ICOSAHEDRON -----
C
C     1     0.     COSA   SINA
C     2     0.     COSA  -SINA
C     3     0.    -COSA   SINA
C     4     0.    -COSA  -SINA
C     5     SINA   0.     COSA
C     6    -SINA   0.     COSA
C     7     SINA   0.    -COSA
C     8    -SINA   0.    -COSA
C     9     COSA   SINA   0.
C    10     COSA  -SINA   0.
C    11    -COSA   SINA   0.
C    12    -COSA  -SINA   0.
C
      ANG=ACOS(ONE/SQRT(FIVE))/TWO
      COSA=COS(ANG)
      SINA=SIN(ANG)
      NAPEX=12
      DO IAPEX=1,NAPEX
         APEX(1,IAPEX)=COSA*C(1,IAPEX)+SINA*S(1,IAPEX)
         APEX(2,IAPEX)=COSA*C(2,IAPEX)+SINA*S(2,IAPEX)
         APEX(3,IAPEX)=COSA*C(3,IAPEX)+SINA*S(3,IAPEX)
      ENDDO
      IF(OUT) THEN
         WRITE(IW,9999)
         DO IAPEX=1,NAPEX
            WRITE(IW,9998) IAPEX,APEX(1,IAPEX),APEX(2,IAPEX),
     1                           APEX(3,IAPEX)
         ENDDO
      ENDIF
C
      NFACE=20
      DO LFACE=1,NFACE
         IJKFAC(1,LFACE)=IJK(1,LFACE)
         IJKFAC(2,LFACE)=IJK(2,LFACE)
         IJKFAC(3,LFACE)=IJK(3,LFACE)
         IJKSEG(  LFACE)=      LFACE
      ENDDO
C
      IF(OUT) THEN
         WRITE(IW,*) '...... END OF -CAV_SPHICO- ......'
      ENDIF
      RETURN
 9999 FORMAT(/,1X,'  APEX',5X,'X',6X,5X,'Y',6X,5X,'Z',6X,/,1X,42('-'))
 9998 FORMAT(1X,I6,F12.8,F12.8,F12.8)
 9997 FORMAT(/,1X,'ICOSAHEDRON',/,1X,11('-'))
      END
C*MODULE SVPINP  *DECK CAV_SPHAPX
      SUBROUTINE CAV_SPHAPX(APEX,MAPEX,IJKFAC,IJKSEG,MFACE,LFACE,
     1                             IJKOLD,IJKNEW,DMIN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL OUT
      LOGICAL DUPLIC
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      DIMENSION   APEX(3,*)
      DIMENSION IJKFAC(3,*)
      DIMENSION IJKSEG(  *)
      DIMENSION IJKOLD(3)
      DIMENSION IJKNEW(3)
      DIMENSION    XYZ(3,3)
      DIMENSION      D(3)
      DIMENSION XYZIJK(3)
      EQUIVALENCE (XYZ(1,1),XIJ),(XYZ(2,1),YIJ),(XYZ(3,1),ZIJ),
     1            (XYZ(1,2),XJK),(XYZ(2,2),YJK),(XYZ(3,2),ZJK),
     2            (XYZ(1,3),XKI),(XYZ(2,3),YKI),(XYZ(3,3),ZKI)
      DATA ZERO  /0.0D+00/
      DATA PT5   /0.5D+00/
      DATA ONE   /1.0D+00/
      DATA THREE /3.0D+00/
      DATA TOL   /1.0D-04/
C
      DIST(X1,Y1,Z1,X2,Y2,Z2)=SQRT((X2-X1)**2+(Y2-Y1)**2+(Z2-Z1)**2)
C
      OUT=.FALSE.
C
C     ----- CREATE MID-POINT OF THE 3 EDGES -----
C
      IOLD=IJKOLD(1)
      JOLD=IJKOLD(2)
      KOLD=IJKOLD(3)
      DO M=1,3
         XYZ(M,1)=(APEX(M,IOLD)+APEX(M,JOLD))*PT5
         XYZ(M,2)=(APEX(M,JOLD)+APEX(M,KOLD))*PT5
         XYZ(M,3)=(APEX(M,KOLD)+APEX(M,IOLD))*PT5
      ENDDO
C
C     ----- PROJECT ONTO SPHERE -----
C
      D(1)=DIST(XIJ,YIJ,ZIJ,ZERO,ZERO,ZERO)
      D(2)=DIST(XJK,YJK,ZJK,ZERO,ZERO,ZERO)
      D(3)=DIST(XKI,YKI,ZKI,ZERO,ZERO,ZERO)
      D(1)=ONE/D(1)
      D(2)=ONE/D(2)
      D(3)=ONE/D(3)
      DO L=1,3
         DO M=1,3
            XYZ(M,L)=XYZ(M,L)*D(L)
         ENDDO
      ENDDO
C
C     ----- CHECK FOR DUPLICATE APECES -----
C
      NEWAPX=0
      DO IAPX=1,3
         DUPLIC=.FALSE.
         DO LAPEX=1,MAPEX
            DD=DIST(XYZ(1,  IAPX),XYZ(2,  IAPX),XYZ(3,  IAPX),
     1              APEX(1,LAPEX),APEX(2,LAPEX),APEX(3,LAPEX))
            IF(ABS(DD).LT.TOL) THEN
               DUPLIC=.TRUE.
               LDUPLC=LAPEX
            ENDIF
         ENDDO
         IF(DUPLIC) THEN
            IJKNEW(IAPX)=LDUPLC
            IF(OUT) THEN
               WRITE(IW,9999) IAPX,IJKOLD,LDUPLC
            ENDIF
         ELSE
            NEWAPX=NEWAPX+1
            JAPX=MAPEX+NEWAPX
            IJKNEW(IAPX)=JAPX
            DO M=1,3
               APEX(M,JAPX)=XYZ(M,IAPX)
            ENDDO
            IF(OUT) THEN
               WRITE(IW,9998) IAPX,IJKOLD,JAPX,
     1                        APEX(1,JAPX),APEX(2,JAPX),APEX(3,JAPX)
            ENDIF
         ENDIF
      ENDDO
      MAPEX=MAPEX+NEWAPX
C
C     ----- MAKE UP NEW FACES AND THEIR CENTERS -----
C
      THIRD=ONE/THREE
      DMIN =ONE
C
      MFACE=MFACE+1
      IJKSEG(  MFACE)=LFACE
      IJKFAC(1,MFACE)=IJKOLD(1)
      IJKFAC(2,MFACE)=IJKNEW(1)
      IJKFAC(3,MFACE)=IJKNEW(3)
      DO M=1,3
         XYZIJK(M)=(APEX(M,IOLD)+APEX(M,JOLD)+APEX(M,KOLD))*THIRD
      ENDDO
      DIJK=DIST(XYZIJK(1),XYZIJK(2),XYZIJK(3),ZERO,ZERO,ZERO)
      DMIN=MIN(DMIN,DIJK)
C
      MFACE=MFACE+1
      IJKSEG(  MFACE)=LFACE
      IJKFAC(1,MFACE)=IJKOLD(2)
      IJKFAC(2,MFACE)=IJKNEW(1)
      IJKFAC(3,MFACE)=IJKNEW(2)
      DO M=1,3
         XYZIJK(M)=(APEX(M,IOLD)+APEX(M,JOLD)+APEX(M,KOLD))*THIRD
      ENDDO
      DIJK=DIST(XYZIJK(1),XYZIJK(2),XYZIJK(3),ZERO,ZERO,ZERO)
      DMIN=MIN(DMIN,DIJK)
C
      MFACE=MFACE+1
      IJKSEG(  MFACE)=LFACE
      IJKFAC(1,MFACE)=IJKOLD(3)
      IJKFAC(2,MFACE)=IJKNEW(2)
      IJKFAC(3,MFACE)=IJKNEW(3)
      DO M=1,3
         XYZIJK(M)=(APEX(M,IOLD)+APEX(M,JOLD)+APEX(M,KOLD))*THIRD
      ENDDO
      DIJK=DIST(XYZIJK(1),XYZIJK(2),XYZIJK(3),ZERO,ZERO,ZERO)
      DMIN=MIN(DMIN,DIJK)
C
      MFACE=MFACE+1
      IJKSEG(  MFACE)=LFACE
      IJKFAC(1,MFACE)=IJKNEW(1)
      IJKFAC(2,MFACE)=IJKNEW(2)
      IJKFAC(3,MFACE)=IJKNEW(3)
      DO M=1,3
         XYZIJK(M)=(APEX(M,IOLD)+APEX(M,JOLD)+APEX(M,KOLD))*THIRD
      ENDDO
      DIJK=DIST(XYZIJK(1),XYZIJK(2),XYZIJK(3),ZERO,ZERO,ZERO)
      DMIN=MIN(DMIN,DIJK)
C
      IF(OUT) THEN
         WRITE(IW,9997) DMIN,MFACE
      ENDIF
C
      RETURN
 9999 FORMAT(' DUPLICATED APEX =',I2,' FOR FACE ',3I5,'. SAME AS = ',I5)
 9998 FORMAT('    NEW     APEX =',I2,' FOR FACE ',3I5,'.  NEWAPX = ',I5,
     1       /,7X,3F12.8)
 9997 FORMAT(' --- DMIN = ',F12.8,' --- MFACE = ',I10)
      END
C*MODULE SVPINP  *DECK MOLAXS
      SUBROUTINE MOLAXS(A,VEC,EIG,NVEC,N,NDIM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- ROUTINE TO SUBSTITUTE DIAGIV FOR DIAGONALIZATION -----
C           OF SYMMETRIC 3X3 MATRIX A  IN TRIANGULAR FORM
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      DIMENSION A(3,3),VEC(3,NVEC),EIG(3),B(3,2),ACOPY(9),JB(3)
      DATA CONV   /1.0D-10/
C
      IF(N.EQ.3.AND.NDIM.EQ.3) GO TO 10
         WRITE(IW,9999)
         CALL ABRT
   10 CONTINUE
C
      DO 30 I=1,3
         DO 20 J=1,3
            VEC(J,I)=ZERO
   20    CONTINUE
         VEC(I,I)=ONE
   30 CONTINUE
      IJ=0
      DO I=1,3
        DO J=1,I
           IJ=IJ+1
           ACOPY(IJ) = A(I,J)
        ENDDO
      ENDDO
C----      CALL DIAJAC(A,VEC,EIG,N,N,B)
      CALL JACDG(ACOPY,VEC,EIG,JB,B,N,N)
C
C     ----- CHECK FOR RIGHT HANDEDNESS, CORRECT IF NOT -----
C
      TEST =   VEC(1,3)*( VEC(2,1)*VEC(3,2) - VEC(3,1)*VEC(2,2) )
     1       + VEC(2,3)*( VEC(3,1)*VEC(1,2) - VEC(1,1)*VEC(3,2) )
     2       + VEC(3,3)*( VEC(1,1)*VEC(2,2) - VEC(2,1)*VEC(1,2) )
      IF(TEST.GT.ZERO) RETURN
      IF( ABS(EIG(1)-EIG(2)).GT.CONV) GO TO 60
         T = EIG(1)
         EIG(1) = EIG(2)
         EIG(2) = T
         DO 50 I=1,3
            T = VEC(I,1)
            VEC(I,1) = VEC(I,2)
            VEC(I,2) = T
   50    CONTINUE
         RETURN
   60 IF( ABS(EIG(2)-EIG(3)).GT.CONV) GO TO 80
         T = EIG(2)
         EIG(2) = EIG(3)
         EIG(3) = T
         DO 70 I=1,3
            T = VEC(I,2)
            VEC(I,2) = VEC(I,3)
            VEC(I,3) = T
   70    CONTINUE
         RETURN
   80 DO 90 I=1,3
         VEC(I,3) = - VEC(I,3)
   90 CONTINUE
      RETURN
 9999 FORMAT(/,' -DIAAXS- DIAGONALIZATION ONLY SET UP FOR 3X3 MATRIX')
      END
C*MODULE SVPINP  *DECK MOLINR
      SUBROUTINE MOLINR(C,C1,NAT,CI,AI,TR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ----- CENTER AND MOMENTS OF INERTIA -----
C
      PARAMETER (MXATM=2000)
      PARAMETER (ZERO=0.0D+00)
C
      LOGICAL OUT
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /MASSES/ ZMASS(MXATM)
C
      DIMENSION C(3,*),C1(3,*)
      DIMENSION CI(3),AI(3,3)
      DIMENSION TR(3)
C
      OUT=.FALSE.
C
C     -----                GET ATOMIC MASSES                -----
C     ----- CALCULATE CENTER OF MASS AND MOMENTS OF INERTIA -----
C
      DO J=1,3
         CI(J)=ZERO
         DO I=1,3
            AI(J,I)=ZERO
         ENDDO
      ENDDO
      AMASS=ZERO
      DO IAT=1,NAT
         AMASS=AMASS+ZMASS(IAT)
         DO I=1,3
            CI(I)=CI(I)+ZMASS(IAT)*C(I,IAT)
         ENDDO
      ENDDO
      DO I=1,3
         CI(I)=CI(I)/AMASS
         TR(I)=CI(I)
      ENDDO
      DO IAT=1,NAT
         DO I=1,3
            C1(I,IAT)=C(I,IAT)-TR(I)
         ENDDO
      ENDDO
C
      DO IAT=1,NAT
         WT=ZMASS(IAT)
         X =C1(1,IAT)
         Y =C1(2,IAT)
         Z =C1(3,IAT)
         AI(1,1)=AI(1,1)+WT*(Y*Y+Z*Z)
         AI(2,1)=AI(2,1)-WT* X*Y
         AI(1,2)=AI(2,1)
         AI(3,1)=AI(3,1)-WT* X*Z
         AI(1,3)=AI(3,1)
         AI(2,2)=AI(2,2)+WT*(X*X+Z*Z)
         AI(3,2)=AI(3,2)-WT* Y*Z
         AI(2,3)=AI(3,2)
         AI(3,3)=AI(3,3)+WT*(X*X+Y*Y)
      ENDDO
      IF(OUT) THEN
         WRITE(IW,9999)
         WRITE(IW,9998) (CI(I),I=1,3)
         WRITE(IW,9997)
         CALL PRSQ(AI,3,3,3)
      ENDIF
C
      RETURN
 9999 FORMAT(/,1X,'CENTER OF MASS',/,1X,14('-'))
 9998 FORMAT(' X = ',F12.8,' Y = ',F12.8,' Z = ',F12.8)
 9997 FORMAT(/,1X,'MOMENTS OF INERTIA',/,1X,18('-'))
      END
