mopac606_nbo/src/esp.f

   1       SUBROUTINE ESP
   2       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   3       INCLUDE 'SIZES'
   4 C***********************************************************************
   5 C
   6 C     THIS IS A DRIVER ROUTINE FOR ELECTROSTATIC POTENTIAL GENERATION
   7 C     WRITTEN BY K.M.MERZ FEB. 1989 AT UCSF
   8 C
   9 C***********************************************************************
  10       COMMON /KEYWRD/ KEYWRD
  11       CHARACTER*241 KEYWRD
  12 C
  13 C     SET STANDARD PARAMETERS FOR THE SURFACE GENERATION
  14 C
  15       IF(INDEX(KEYWRD,'SCALE=') .NE. 0)THEN
  16          SCALE = READA(KEYWRD,INDEX(KEYWRD,'SCALE='))
  17       ELSE
  18          SCALE = 1.4D0
  19       ENDIF
  20 C
  21       IF(INDEX(KEYWRD,'DEN=') .NE. 0)THEN
  22          DEN = READA(KEYWRD,INDEX(KEYWRD,'DEN='))
  23       ELSE
  24          DEN = 1.0D0
  25       ENDIF
  26 C
  27       IF(INDEX(KEYWRD,'SCINCR=') .NE. 0)THEN
  28          SCINCR = READA(KEYWRD,INDEX(KEYWRD,'SCINCR='))
  29       ELSE
  30          SCINCR = 0.20D0
  31       ENDIF
  32 C
  33       IF(INDEX(KEYWRD,'NSURF=') .NE. 0)THEN
  34          N = READA(KEYWRD,INDEX(KEYWRD,'NSURF='))
  35       ELSE
  36          N = 4
  37       ENDIF
  38 C
  39       TIME1=SECOND()
  40 C
  41 C     NOW CALCULATE THE SURFACE POINTS
  42 C
  43       IF(INDEX(KEYWRD,'WILLIAMS') .NE. 0) THEN
  44          CALL PDGRID
  45       ELSE
  46          DO 10 I = 1,N
  47             CALL SURFAC(SCALE,DEN,I)
  48             SCALE = SCALE + SCINCR
  49    10    CONTINUE
  50       ENDIF
  51 C
  52 C     NEXT CALCULATE THE ESP AT THE POINTS CALCULATED BY SURFAC
  53 C
  54       CALL POTCAL
  55 C
  56 C     END OF CALCULATION
  57 C
  58       TIME1=SECOND()-TIME1
  59       WRITE(6,20) 'TIME TO CALCULATE ESP:',TIME1,' SECONDS'
  60    20 FORMAT(/9X,A,F8.2,A)
  61       RETURN
  62       END
  63       SUBROUTINE PDGRID
  64 C
  65 C     ROUTINE TO CALCULATE WILLIAMS SURFACE
  66 C
  67       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  68       INCLUDE 'SIZES'
  69       DIMENSION IZ(100),XYZ(3,100),VDERW(53),DIST(100)
  70       DIMENSION XMIN(3),XMAX(3),COORD(3,NUMATM)
  71       COMMON /GEOM/   GEO(3,NUMATM)
  72       COMMON /GEOKST/ NATOMS,LABELS(NUMATM), NABC(3*NUMATM)
  73 C
  74       COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
  75       COMMON /WORK1/    POTPT(3,MESP), WORK1D(4*MESP)
  76       COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
  77 C
  78       DATA VDERW/53*0.0D0/
  79       VDERW(1)=2.4D0
  80       VDERW(5)=3.0D0
  81       VDERW(6)=2.9D0
  82       VDERW(7)=2.7D0
  83       VDERW(8)=2.6D0
  84       VDERW(9)=2.55D0
  85       VDERW(15)=3.1D0
  86       VDERW(16)=3.05D0
  87       VDERW(17)=3.0D0
  88       VDERW(35)=3.15D0
  89       VDERW(53)=3.35D0
  90       SHELL=1.2D0
  91       NESP=0
  92       GRID=0.8D0
  93       CLOSER=0.D0
  94 C     CHECK IF VDERW IS DEFINED FOR ALL ATOMS
  95 C
  96 C     CONVERT INTERNAL TO CARTESIAN COORDINATES
  97 C
  98       CALL GMETRY(GEO,COORD)
  99 C
 100 C     STRIP COORDINATES AND ATOM LABEL FOR DUMMIES (I.E. 99)
 101 C
 102       ICNTR = 0
 103       DO 20 I=1,NATOMS
 104          DO 10 J=1,3
 105    10    CO(J,I) = COORD(J,I)
 106          IF(LABELS(I) .EQ. 99) GOTO 20
 107          ICNTR = ICNTR + 1
 108          IAN(ICNTR) = LABELS(I)
 109    20 CONTINUE
 110       NATOM=ICNTR
 111 C
 112       DO 30 I=1,NATOM
 113          J=IAN(I)
 114          IF (VDERW(J).EQ.0.0D0) GO TO 40
 115    30 CONTINUE
 116       GO TO 50
 117    40 CONTINUE
 118       WRITE(6,*) 'VAN DER WAALS'' RADIUS NOT DEFINED FOR ATOM',I
 119       WRITE(6,*) 'IN WILLIAMS SURFACE ROUTINE PDGRID!'
 120       STOP
 121 C     NOW CREATE LIMITS FOR A BOX
 122    50 DO 100 IX = 1,3
 123          XMIN(IX)= 100000.0D0
 124          XMAX(IX)=-100000.0D0
 125          DO 90 IA = 1,NATOM
 126             IF (CO(IX,IA)-XMIN(IX))60,70,70
 127    60       XMIN(IX)=CO(IX,IA)
 128    70       IF (CO(IX,IA)-XMAX(IX))90,90,80
 129    80       XMAX(IX)=CO(IX,IA)
 130    90    CONTINUE
 131   100 CONTINUE
 132 C     ADD (OR SUBTRACT) THE MAXIMUM VDERW PLUS SHELL
 133       VDMAX=0.0D0
 134       DO 110 I=1,53
 135          IF (VDERW(I).GT.VDMAX) VDMAX=VDERW(I)
 136   110 CONTINUE
 137       DO 120 I=1,3
 138          XMIN(I)=XMIN(I)-VDMAX-SHELL
 139   120 XMAX(I)=XMAX(I)+VDMAX+SHELL
 140 C STEP GRID BACK FROM ZERO TO FIND STARTING POINTS
 141       XSTART=0.0D0
 142   130 XSTART=XSTART-GRID
 143       IF (XSTART.GT.XMIN(1)) GO TO 130
 144       YSTART=0.0D0
 145   140 YSTART=YSTART-GRID
 146       IF (YSTART.GT.XMIN(2)) GO TO 140
 147       ZSTART=0.0D0
 148   150 ZSTART=ZSTART-GRID
 149       IF (ZSTART.GT.XMIN(3)) GO TO 150
 150       NPNT=0
 151       ZGRID=ZSTART
 152   160 YGRID=YSTART
 153   170 XGRID=XSTART
 154   180 DO 190 L=1,NATOM
 155          JZ=IAN(L)
 156          DIST(L)=SQRT((CO(1,L)-XGRID)**2+(CO(2,L)-YGRID)**2+
 157      1 (CO(3,L)-ZGRID)**2)
 158 C     REJECT GRID POINT IF ANY ATOM IS TOO CLOSE
 159          IF(DIST(L).LT.(VDERW(JZ)-CLOSER)) GO TO 220
 160   190 CONTINUE
 161 C BUT AT LEAST ONE ATOM MUST BE CLOSE ENOUGH
 162       DO 200 L=1,NATOM
 163          JZ=IAN(L)
 164          IF(DIST(L).GT.(VDERW(JZ)+SHELL)) GO TO 200
 165          GO TO 210
 166   200 CONTINUE
 167       GO TO 220
 168   210 NPNT=NPNT+1
 169       NESP=NESP+1
 170       POTPT(1,NESP)=XGRID
 171       POTPT(2,NESP)=YGRID
 172       POTPT(3,NESP)=ZGRID
 173   220 XGRID=XGRID+GRID
 174       IF (XGRID.LE.XMAX(1)) GO TO 180
 175       YGRID=YGRID+GRID
 176       IF (YGRID.LE.XMAX(2)) GO TO 170
 177       ZGRID=ZGRID+GRID
 178       IF (ZGRID.LE.XMAX(3)) GO TO 160
 179       RETURN
 180       END
 181 C***********************************************************************
 182       SUBROUTINE SURFAC(SCALE,DENS,IPT)
 183       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 184       INCLUDE 'SIZES'
 185 C***********************************************************************
 186 C
 187 C      THIS SUBROUTINE CALCULATES THE MOLECULAR SURFACE OF A MOLECULE
 188 C      GIVEN THE COORDINATES OF ITS ATOMS.  VAN DER WAALS' RADII FOR
 189 C      THE ATOMS AND THE PROBE RADIUS MUST ALSO BE SPECIFIED.
 190 C
 191 C      ON INPUT    SCALE = INITIAL VAN DER WAALS' SCALE FACTOR
 192 C                  DENS  = DENSITY OF POINTS PER UNIT AREA
 193 C
 194 C      THIS SUBROUTINE WAS LIFTED FROM MICHAEL CONNOLLY'S SURFACE
 195 C      PROGRAM FOR UCSF GRAPHICS SYSTEM BY U.CHANDRA SINGH AND
 196 C      P.A.KOLLMAN AND MODIFIED FOR USE IN QUEST. K.M.MERZ
 197 C      ADAPTED AND CLEANED UP THIS PROGRAM FOR USE IN AMPAC/MOPAC
 198 C      IN FEB. 1989 AT UCSF.
 199 C
 200 C***********************************************************************
 201       COMMON /GEOM/   GEO(3,NUMATM)
 202       COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
 203      1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
 204       COMMON /KEYWRD/ KEYWRD
 205 C
 206       COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
 207       COMMON /WORK1/    POTPT(3,MESP), PAD1(2*MESP), RAD(MESP),
 208      1IAS(MESP)
 209       COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
 210 C
 211       CHARACTER*241 KEYWRD
 212 C
 213 C     CARTESIAN COORDINATE AND ATOM LABELS
 214 C
 215       DIMENSION COORD(3,NUMATM),VANDER(100)
 216       DIMENSION CON(3,1000),ROT(3,3)
 217 C
 218 C     NEIGHBOR ARRAYS
 219 C
 220 C     THIS SAME DIMENSION FOR THE MAXIMUM NUMBER OF NEIGHBORS
 221 C     IS USED TO DIMENSION ARRAYS IN THE LOGICAL FUNCTION COLLID
 222 C
 223       DIMENSION INBR(200),CNBR(3,200),RNBR(200)
 224       LOGICAL SNBR(200),MNBR(200)
 225 C
 226 C     ARRAYS FOR ALL ATOMS
 227 C
 228 C     IATOM, JATOM AND KATOM COORDINATES
 229 C
 230       DIMENSION CI(3), IELDAT(56), TEMP0(3)
 231 C
 232 C     GEOMETRIC CONSTRUCTION VECTORS
 233 C
 234       DIMENSION CW(3,2)
 235 C
 236 C     LOGICAL VARIABLES
 237 C
 238       LOGICAL SI
 239 C
 240 C     LOGICAL FUNCTIONS
 241 C
 242       LOGICAL COLLID
 243 C
 244 C     DATA FOR VANDER VALL RADII
 245 C
 246       CHARACTER MARKER*3, MARKSS*3, MYNAM*3, IELDAT*4, NAMATM*4
 247       DATA VANDER/1.20D0,1.20D0,1.37D0,1.45D0,1.45D0,1.50D0,1.50D0,
 248      1            1.40D0,1.35D0,1.30D0,1.57D0,1.36D0,1.24D0,1.17D0,
 249      2            1.80D0,1.75D0,1.70D0,17*0.0D0,2.3D0,65*0.0D0/
 250       DATA MARKER/'A  '/,MARKSS/'SS0'/,MYNAM/'UC '/
 251 C
 252       DATA IELDAT/'  BQ','  H ','  HE','  LI','  BE','  B ',
 253      1            '  C ','  N ','  O ','  F ','  NE','  NA',
 254      2            '  MG','  AL','  SI','  P ','  S ','  CL',
 255      3            '  AR','  K ','  CA','  SC','  TI','  V ',
 256      4            '  CR','  MN','  FE','  CO','  NI','  CU',
 257      5            '  ZN','  GA','  GE','  AS','  SE','  BR',
 258      6            '  KR','  RB','  SR','   Y','  ZR','  NB',
 259      7            '  MO','  TC','  RU','  RH','  PD','  AG',
 260      8            '  CD','  IN','  SN','  SB','  TE','   I',
 261      9            '   X','  CS'/
 262       PI=4.D0*ATAN(1.D0)
 263 C     INSERT VAN DER WAAL RADII FOR ZINC
 264       VANDER(30)=1.00D0
 265 C
 266 C     CONVERT INTERNAL TO CARTESIAN COORDINATES
 267 C
 268       CALL GMETRY(GEO,COORD)
 269 C
 270 C     STRIP COORDINATES AND ATOM LABEL FOR DUMMIES (I.E. 99)
 271 C
 272       ICNTR = 0
 273       DO 20 I=1,NATOMS
 274          DO 10 J=1,3
 275    10    CO(J,I) = COORD(J,I)
 276          IF(LABELS(I) .EQ. 99) GOTO 20
 277          ICNTR = ICNTR + 1
 278          IAN(ICNTR) = LABELS(I)
 279    20 CONTINUE
 280 C
 281 C     ONLY VAN DER WAALS' TYPE SURFACE IS GENERATED
 282 C
 283       IOP = 1
 284       RW =0.0D0
 285       NATOM = ICNTR
 286       DEN = DENS
 287       DO 30 I=1,NATOM
 288          IPOINT = IAN(I)
 289          RAD(I) = VANDER(IPOINT)*SCALE
 290          IF (RAD(I) .LT. 0.01D0) THEN
 291             WRITE(6,'(T2,''VAN DER WAALS'''' RADIUS FOR ATOM '',I3,
 292      1         '' IS ZERO, SUPPLY A VALUE IN SUBROUTINE SURFAC)''
 293      2         )')
 294          ENDIF
 295          IAS(I) = 2
 296    30 CONTINUE
 297 C
 298 C     BIG LOOP FOR EACH ATOM
 299 C
 300       DO 110 IATOM = 1, NATOM
 301          IF (IAS(IATOM) .EQ. 0) GO TO 110
 302 C
 303 C     TRANSFER VALUES FROM LARGE ARRAYS TO IATOM VARIABLES
 304 C
 305          NAMATM =IELDAT(IAN(IATOM)+1)
 306          RI = RAD(IATOM)
 307          SI = IAS(IATOM) .EQ. 2
 308          DO 40 K = 1,3
 309             CI(K) = CO(K,IATOM)
 310    40    CONTINUE
 311 C
 312 C     GATHER THE NEIGHBORING ATOMS OF IATOM
 313 C
 314          NNBR = 0
 315          DO 60 JATOM = 1, NATOM
 316             IF (IATOM .EQ. JATOM .OR. IAS(JATOM) .EQ. 0) GO TO 60
 317             D2 = DIST2(CI,CO(1,JATOM))
 318             IF (D2 .GE. (2*RW+RI+RAD(JATOM)) ** 2) GO TO 60
 319 C
 320 C     WE HAVE A NEW NEIGHBOR
 321 C     TRANSFER ATOM COORDINATES, RADIUS AND SURFACE REQUEST NUMBER
 322 C
 323             NNBR = NNBR + 1
 324             IF (NNBR .GT. 200)THEN
 325                WRITE (6,'(''ERROR'',2X,''TOO MANY NEIGHBORS:'',I5)')NNBR
 326                STOP
 327             ENDIF
 328             INBR(NNBR) = JATOM
 329             DO 50 K = 1,3
 330                CNBR(K,NNBR) = CO(K,JATOM)
 331    50       CONTINUE
 332             RNBR(NNBR) = RAD(JATOM)
 333             SNBR(NNBR) = IAS(JATOM) .EQ. 2
 334    60    CONTINUE
 335 C
 336 C     CONTACT SURFACE
 337 C
 338          IF (.NOT. SI) GO TO 110
 339          NCON = (4 * PI * RI ** 2) * DEN
 340          IF (NCON .GT. 1000) NCON = 1000
 341 C
 342 C     THIS CALL MAY DECREASE NCON SOMEWHAT
 343 C
 344          IF ( NCON .EQ. 0) THEN
 345             WRITE(6,'(T2,''VECTOR LENGTH OF ZERO IN SURFAC'')')
 346             STOP
 347          ENDIF
 348          CALL GENUN(CON,NCON)
 349          AREA = (4 * PI * RI ** 2) / NCON
 350 C
 351 C     CONTACT PROBE PLACEMENT LOOP
 352 C
 353          DO 100 I = 1,NCON
 354             DO 70 K = 1,3
 355                CW(K,1) = CI(K) + (RI + RW) * CON(K,I)
 356    70       CONTINUE
 357 C
 358 C     CHECK FOR COLLISION WITH NEIGHBORING ATOMS
 359 C
 360             IF (COLLID(CW(1,1),RW,CNBR,RNBR,MNBR,NNBR,1,
 361      1      JNBR,KNBR)) GO TO 100
 362             DO 80 KK=1,3
 363                TEMP0(KK) =CI(KK)+RI*CON(KK,I)
 364    80       CONTINUE
 365 C
 366 C     STORE POINT IN POTPT AND INCREMENT NESP
 367 C
 368             NESP = NESP + 1
 369             IF (NESP .GT. MESP) THEN
 370                WRITE(6,90)
 371    90          FORMAT(/'ERROR - TO MANY POINTS GENERATED IN SURFAC')
 372                WRITE(6,'(''    REDUCE NSURF, SCALE, DEN, OR SCINCR'')')
 373                STOP
 374             ENDIF
 375             POTPT(1,NESP) = TEMP0(1)
 376             POTPT(2,NESP) = TEMP0(2)
 377             POTPT(3,NESP) = TEMP0(3)
 378   100    CONTINUE
 379   110 CONTINUE
 380       RETURN
 381       END
 382 C****************************************************************
 383       FUNCTION DIST2(A,B)
 384 C
 385 C     DETERMINE DISTANCES BETWEEN NEIGHBORING ATOMS
 386 C
 387       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 388       DIMENSION A(3)
 389       DIMENSION B(3)
 390       DIST2 = (A(1)-B(1))**2 + (A(2)-B(2))**2 + (A(3)-B(3))**2
 391       RETURN
 392       END
 393 C****************************************************************
 394       LOGICAL FUNCTION COLLID(CW,RW,CNBR,RNBR,MNBR,NNBR,ISHAPE,
 395      1JNBR,KNBR)
 396 C****************************************************************
 397 C
 398 C     COLLISION CHECK OF PROBE WITH NEIGHBORING ATOMS
 399 C     USED BY SURFAC ONLY.
 400 C
 401 C****************************************************************
 402       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 403       DIMENSION CW(3)
 404       DIMENSION CNBR(3,200)
 405       DIMENSION RNBR(200)
 406       LOGICAL MNBR(200)
 407       IF (NNBR .LE. 0) GO TO 20
 408 C
 409 C     CHECK WHETHER PROBE IS TOO CLOSE TO ANY NEIGHBOR
 410 C
 411       DO 10 I = 1, NNBR
 412          IF (ISHAPE .GT. 1 .AND. I .EQ. JNBR) GO TO 10
 413          IF (ISHAPE .EQ. 3 .AND. (I .EQ. KNBR .OR. .NOT. MNBR(I)))
 414      1   GO TO 10
 415          SUMRAD = RW + RNBR(I)
 416          VECT1 = DABS(CW(1) - CNBR(1,I))
 417          IF (VECT1 .GE. SUMRAD) GO TO 10
 418          VECT2 = DABS(CW(2) - CNBR(2,I))
 419          IF (VECT2 .GE. SUMRAD) GO TO 10
 420          VECT3 = DABS(CW(3) - CNBR(3,I))
 421          IF (VECT3 .GE. SUMRAD) GO TO 10
 422          SR2 = SUMRAD ** 2
 423          DD2 = VECT1 ** 2 + VECT2 ** 2 + VECT3 ** 2
 424          IF (DD2 .LT. SR2) GO TO 30
 425    10 CONTINUE
 426    20 CONTINUE
 427       COLLID = .FALSE.
 428       GO TO 40
 429    30 CONTINUE
 430       COLLID = .TRUE.
 431    40 CONTINUE
 432       RETURN
 433       END
 434 C****************************************************************
 435       SUBROUTINE GENUN(U,N)
 436 C****************************************************************
 437 C
 438 C     GENERATE UNIT VECTORS OVER SPHERE. USED BY SURFAC ONLY.
 439 C
 440 C****************************************************************
 441       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 442       DIMENSION U(3,N)
 443       PI=4.D0*ATAN(1.D0)
 444       NEQUAT = SQRT(N * PI)
 445       NVERT = NEQUAT/2
 446       NU = 0
 447       DO 20 I = 1,NVERT+1
 448          FI = (PI * (I-1)) / NVERT
 449          Z = COS(FI)
 450          XY = SIN(FI)
 451          NHOR = NEQUAT * XY
 452          IF (NHOR .LT. 1) NHOR = 1
 453          DO 10 J = 1,NHOR
 454             FJ = (2.D0 * PI * (J-1)) / NHOR
 455             X = DCOS(FJ) * XY
 456             Y = DSIN(FJ) * XY
 457             IF (NU .GE. N) GO TO 30
 458             NU = NU + 1
 459             U(1,NU) = X
 460             U(2,NU) = Y
 461             U(3,NU) = Z
 462    10    CONTINUE
 463    20 CONTINUE
 464    30 CONTINUE
 465       N = NU
 466       RETURN
 467       END
 468 C***********************************************************************
 469       SUBROUTINE POTCAL
 470       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 471       INCLUDE 'SIZES'
 472 C***********************************************************************
 473 C
 474 C     THIS SUBROUTINE CALCULATES THE TOTAL ELECTROSTATIC POTENTIAL
 475 C     THE NUCLEAR CONTRIBUTION IS EVALUATED BY NUCPOT
 476 C     THE ELECTRONIC CONTRIBUTION IS EVALUATED BY ELESP
 477 C     ESPFIT FITS THE QUANTUM POTENTIAL TO A CLASSICAL POINT CHARGE
 478 C     MODEL.
 479 C     THIS SUBROUTINE WAS WRITTEN BY B.H.BESLER AND K.M.MERZ IN FEB.
 480 C     1989 AT UCSF
 481 C
 482 C***********************************************************************
 483       COMMON /KEYWRD/ KEYWRD
 484       COMMON /CORE/ TORE(107)
 485       COMMON /ELEMTS/ ELEMNT(107)
 486       COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
 487       COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
 488       COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
 489       COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
 490       COMMON /DIPSTO/ UX,UY,UZ,CH(NUMATM)
 491       COMMON /ESPF/  AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
 492      1Q(NUMATM+4),QSC(NUMATM+4),CF, ESPFD(MAXORB**2-NUMATM-5)
 493       CHARACTER*241 KEYWRD
 494       CHARACTER *2  ELEMNT
 495       LOGICAL DEBUG,WRTESP,CEQUIV(NUMATM,NUMATM)
 496 C
 497 C     DEBUG PRINTING - RESULTS IN COPIOUS OUTPUT
 498 C
 499       DEBUG = (INDEX(KEYWRD,'DEBUG') .NE. 0)
 500 C
 501 C
 502       CALL ELESP
 503       BOHR = 0.529167D00
 504 C
 505 C     NOW FIT THE ELECTROSTATIC POTENTIAL
 506 C
 507       WRITE(6,'(//12X,''ELECTROSTATIC POTENTIAL CHARGES'',/)')
 508       IZ=0
 509       IF(INDEX(KEYWRD,'CHARGE=') .NE. 0) IZ=READA(KEYWRD,INDEX(KEYWRD,
 510      1'CHARGE='))
 511 C
 512 C     DIPOLAR CONSTRAINTS IF DESIRED
 513 C
 514       IF(INDEX(KEYWRD,'DIPOLE') .NE. 0) THEN
 515          IDIP = 1
 516          IF(IZ .NE. 0)THEN
 517             IDIP = 0
 518             WRITE(6,'(/12X,''  DIPOLE CONSTRAINTS NOT USED'')')
 519             WRITE(6,'(12X,''        CHARGED MOLECULE'',/)')
 520          ENDIF
 521       ELSE
 522          IDIP = 0
 523       ENDIF
 524       IF (IDIP .EQ. 1) THEN
 525          WRITE(6,'(/12X,''DIPOLE CONSTRAINTS WILL BE USED'',/)')
 526       ENDIF
 527 C
 528 C     GET X,Y,Z DIPOLE COMPONENTS IF DESIRED
 529 C
 530       IF(INDEX(KEYWRD,'DIPX=') .NE. 0) THEN
 531          DX = READA(KEYWRD,INDEX(KEYWRD,'DIPX='))
 532       ELSE
 533          DX = UX
 534       ENDIF
 535       IF(INDEX(KEYWRD,'DIPY=') .NE. 0) THEN
 536          DY = READA(KEYWRD,INDEX(KEYWRD,'DIPY='))
 537       ELSE
 538          DY = UY
 539       ENDIF
 540       IF(INDEX(KEYWRD,'DIPZ=') .NE. 0) THEN
 541          DZ = READA(KEYWRD,INDEX(KEYWRD,'DIPZ='))
 542       ELSE
 543          DZ = UZ
 544       ENDIF
 545       CALL ESPFIT(IDIP,NATOM,NESP,IZ,ESP,POTPT,CO,DX,DY,DZ,RMS,RRMS)
 546 C
 547 C     WRITE OUT OUR RESULTS TO CHANNEL 6
 548 C     THE CHARGES ARE SCALED TO REPRODUCE 6-31G* CHARGES FOR MNDO ONLY
 549 C     AM1 AND MINDO/3 CHARGES ARE NOT SCALED DUE TO THE LOW COORELATION
 550 C     COEFFICIENT. SEE BESLER,MERZ,KOLLMAN IN J. COMPUT. CHEM.
 551 C     (IN PRESS)
 552 C
 553       IF((INDEX(KEYWRD,'AM1') .NE. 0) .OR.
 554      1(INDEX(KEYWRD,'MINDO') .NE. 0) .OR.
 555      2(INDEX(KEYWRD,'PDG') .NE. 0) .OR.
 556      3(INDEX(KEYWRD,'PM3') .NE. 0))THEN
 557          WRITE(6,'(15X,''ATOM NO.    TYPE    CHARGE'')')
 558          DO 10 I=1,NATOM
 559             WRITE(6,'(17X,I2,9X,A2,1X,F10.4)')I,ELEMNT(IAN(I)),Q(I)
 560    10    CONTINUE
 561       ELSE
 562 C
 563 C     MNDO CALCULATION-SCALE THE CHARGES. TEST FOR SLOPE KEYWORD
 564 C
 565          IF(INDEX(KEYWRD,'SLOPE=') .NE. 0) THEN
 566             SLOPE = READA(KEYWRD,INDEX(KEYWRD,'SLOPE='))
 567          ELSE
 568             SLOPE = 1.422D0
 569          ENDIF
 570          DO 20 I=1,NATOM
 571             QSC(I) = SLOPE*Q(I)
 572    20    CONTINUE
 573          WRITE(6,'(7X,''ATOM NO.    TYPE    CHARGE   SCALED CHARGE'')')
 574          DO 30 I=1,NATOM
 575             WRITE(6,'(9X,I2,9X,A2,1X,F10.4,2X,F10.4)')I,ELEMNT(IAN(I
 576      1)),   Q(I),QSC(I)
 577    30    CONTINUE
 578       ENDIF
 579       WRITE(6,'(/12X,A,4X,I6)') 'THE NUMBER OF POINTS IS:',NESP
 580       WRITE(6,'(12X,A,4X,F9.4)') 'THE RMS DEVIATION IS:',RMS
 581       WRITE(6,'(12X,A,3X,F9.4)') 'THE RRMS DEVIATION IS:',RRMS
 582 C
 583 C     CALCULATE DIPOLE MOMENT IF NEUTRAL MOLECULE
 584 C
 585       IF (IZ .NE. 0) THEN
 586          GO TO 60
 587       ELSE
 588          WRITE(6,40)
 589    40    FORMAT (//5X,'DIPOLE MOMENT EVALUATED FROM '
 590      1,'THE POINT CHARGES',/)
 591          DO 50 I=1,NATOM
 592             DIPX=DIPX+CO(1,I)*Q(I)/BOHR
 593             DIPY=DIPY+CO(2,I)*Q(I)/BOHR
 594             DIPZ=DIPZ+CO(3,I)*Q(I)/BOHR
 595    50    CONTINUE
 596          DIP=SQRT(DIPX**2+DIPY**2+DIPZ**2)
 597          WRITE(6,'(12X,'' X        Y        Z       TOTAL'')')
 598          WRITE(6,'(8X,4F9.4)')DIPX*CF,DIPY*CF,DIPZ*CF,DIP*CF
 599       ENDIF
 600    60 CONTINUE
 601 C     DETERMINE WHICH CHARGES SHOULD BE EQUIVALENT BY SYMMETRY AND
 602 C     AVERAGE THEM IF DESIRED
 603       IF(INDEX(KEYWRD,'SYMAVG') .NE. 0) THEN
 604          DO 70 I=1,NATOM
 605             DO 70 J=1,NATOM
 606                CEQUIV(I,J)=.FALSE.
 607                IF(ABS(ABS(CH(I))-ABS(CH(J))) .LT. 1.D-5)  CEQUIV(I,J)=.T
 608      1RUE.
 609    70    CONTINUE
 610          DO 90 I=1,NATOM
 611             IEQ=0
 612             QSC(I)=0.D0
 613             DO 80 J=1,NATOM
 614                IF(CEQUIV(I,J)) THEN
 615                   QSC(I)=QSC(I)+ABS(Q(J))
 616                   IEQ=IEQ+1
 617                ENDIF
 618    80       CONTINUE
 619             CH(I)=Q(I)/ABS(Q(I))*QSC(I)/IEQ
 620    90    CONTINUE
 621          WRITE(6,*) ' '
 622          WRITE(6,*)'   ELECTROSTATIC POTENTIAL CHARGES AVERAGED FOR'
 623          WRITE(6,*)'   SYMMETRY EQUIVALENT ATOMS'
 624          WRITE(6,*) ' '
 625          IF((INDEX(KEYWRD,'AM1') .NE. 0) .OR.
 626      1(INDEX(KEYWRD,'MINDO') .NE. 0) .OR.
 627      2(INDEX(KEYWRD,'PDG') .NE. 0) .OR.
 628      3(INDEX(KEYWRD,'PM3') .NE. 0))THEN
 629             WRITE(6,'(7X,''ATOM NO.    TYPE    CHARGE'')')
 630             DO 100 I=1,NATOM
 631                WRITE(6,'(9X,I2,9X,A2,1X,F10.4)')I,ELEMNT(IAN(I)),
 632      1   CH(I)
 633   100       CONTINUE
 634          ELSE
 635             WRITE(6,'(7X,''ATOM NO.    TYPE    CHARGE   SCALED CHARGE'')
 636      1')
 637             DO 110 I=1,NATOM
 638                WRITE(6,'(9X,I2,9X,A2,1X,F10.4,2X,F10.4)')I,ELEMNT(IA
 639      1N(I)),   CH(I),CH(I)*SLOPE
 640   110       CONTINUE
 641          ENDIF
 642       ENDIF
 643       RETURN
 644       END
 645       SUBROUTINE ELESP
 646       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 647 C***********************************************************************
 648 C     ELESP LOADS THE STO-6G BASIS SET ONTO THE ATOMS, PERFOMS THE
 649 C     DEORTHOGONALIZATION OF THE COEFFICIENTS AND EVALUATES THE
 650 C     ELECTRONIC CONTRIBUTION TO THE ESP. IT WAS WRITTEN BY B.H.BESLER
 651 C     AND K.M.MERZ IN FEB. 1989 AT UCSF.
 652 C
 653 C***********************************************************************
 654       CHARACTER*241 KEYWRD
 655       DOUBLE PRECISION NORM,OVL
 656       LOGICAL CALLED,POTWRT,RST,STO3G
 657       INCLUDE 'SIZES'
 658       COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
 659      1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
 660       COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
 661       COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
 662       COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
 663       COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
 664       COMMON /STO6G/  ALLC(6,5,2),ALLZ(6,5,2)
 665       COMMON /VECTOR/ C(MORB2*2+MAXORB*2)
 666       COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
 667      1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
 668      2                NCLOSE,NOPEN,NDUMY,FRACT
 669       COMMON /KEYWRD/ KEYWRD
 670       COMMON /ESPC/  CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
 671      1                EX(MAXPR),ESPI(MAXORB,MAXORB),
 672      2                FV(0:8,821),FAC(0:7),
 673      3                DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR),
 674      4                OVL(MAXORB,MAXORB),FC(MAXPR*6)
 675      6       /CORE  / TORE(107)
 676      7       /EXPONT/ ZS(107),ZP(107),ZD(107)
 677 *
 678 *  END OF MINDO/3 COMMON BLOCKS
 679 *
 680       COMMON /INDX/   INDC(MAXORB)
 681       DIMENSION CESPM2(MAXORB,MAXORB),SLA(10)
 682       DIMENSION CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB)
 683       DATA BOHR/0.529167D0/
 684       PI=4.D0*ATAN(1.D0)
 685 C
 686 C     PUT STO-6G BASIS SET ON ATOM CENTERS
 687 C
 688       DO 10 I=-1,10
 689          DEX(I)=DEX2(I)
 690    10 CONTINUE
 691       DO 20   I=0,7
 692          FAC(I)=1.D0/FAC(I)
 693    20 CONTINUE
 694       DO 30 M=0,8
 695          K=1
 696          FV(M,1)=1.D0/(2.D0*M+1.D0)
 697          DO 30 T=0.05D0,41.D0,0.05D0
 698             K=K+1
 699             CALL FSUB(M,T,FVAL)
 700             FV(M,K)=FVAL
 701    30 CONTINUE
 702 C
 703 C     LOAD BASIS FUNCTIONS INTO ARRAYS
 704 C
 705       STO3G=(INDEX(KEYWRD,'STO3G') .NE. 0)
 706       IF(STO3G) THEN
 707          ICD=3
 708          CALL SETUP3
 709       ELSE
 710          ICD=6
 711          CALL SETUPG
 712       ENDIF
 713       NC=0
 714       NPR=0
 715       DO 80 I=1,NATOM
 716          IF (IAN(I) .LE. 2) THEN
 717             DO 40 J=1,ICD
 718                CC(NPR+J)=ALLC(J,1,1)
 719                EX(NPR+J)=ALLZ(J,1,1)*ZS(1)**2
 720                CEN(NPR+J,1)=CO(1,I)/BOHR
 721                CEN(NPR+J,2)=CO(2,I)/BOHR
 722                CEN(NPR+J,3)=CO(3,I)/BOHR
 723                IAM(NPR+J,1)=0
 724                IAM(NPR+J,2)=0
 725                FC(NPR+J)=I
 726    40       CONTINUE
 727             NC=NC+1
 728             NPR=NPR+ICD
 729          ELSE
 730 C        DETERMINE PRINCIPAL QUANTUM NUMBER(NQN)
 731 C        OF ORBITALS TO BE USED
 732 C
 733             NQN=2
 734             IF(IAN(I) .GT. 10 .AND. IAN(I) .LE. 18) NQN=3
 735             IF(IAN(I) .GT. 18 .AND. IAN(I) .LE. 36) NQN=4
 736             IF(IAN(I) .GT. 36 .AND. IAN(I) .LE. 54) NQN=5
 737 C
 738             DO 50 J=1,ICD
 739                CC(NPR+J)=ALLC(J,NQN,1)
 740                EX(NPR+J)=ALLZ(J,NQN,1)*ZS(IAN(I))**2
 741                CEN(NPR+J,1)=CO(1,I)/BOHR
 742                CEN(NPR+J,2)=CO(2,I)/BOHR
 743                CEN(NPR+J,3)=CO(3,I)/BOHR
 744                IAM(NPR+J,1)=0
 745                IAM(NPR+J,2)=0
 746    50       CONTINUE
 747             NC=NC+1
 748             NPR=NPR+ICD
 749             DO 70 K=1,3
 750                DO 60  J=1,ICD
 751                   CC(NPR+J)=ALLC(J,NQN,2)
 752                   EX(NPR+J)=ALLZ(J,NQN,2)*ZP(IAN(I))**2
 753                   CEN(NPR+J,1)=CO(1,I)/BOHR
 754                   CEN(NPR+J,2)=CO(2,I)/BOHR
 755                   CEN(NPR+J,3)=CO(3,I)/BOHR
 756                   IAM(NPR+J,1)=1
 757                   IAM(NPR+J,2)=K
 758    60          CONTINUE
 759                NC=NC+1
 760                NPR=NPR+ICD
 761    70       CONTINUE
 762          ENDIF
 763    80 CONTINUE
 764 C
 765 C     CALCULATE NORMALIZATION CONSTANTS AND INCLUDE
 766 C     THEM IN THE CONTRACTION COEFFICIENTS
 767 C
 768       DO 90 I=1,NPR
 769          NORM=(2.D0*EX(I)/PI)**0.75D0*(4.D0*EX(I))**(IAM(I,1)/2.D0)/
 770      1   SQRT(DEX(2*IAM(I,1)-1))
 771          CC(I)=CC(I)*NORM
 772    90 CONTINUE
 773       IPR=0
 774 C
 775 C     PERFORM SORT OF PRIMITIVES BY ANGULAR MOMENTUM
 776 C
 777       IS=0
 778       IP=0
 779       IPC=0
 780       ISC=0
 781       J=0
 782       DO 100 I=1,NPR
 783          IF (IAM(I,1) .EQ. 0) THEN
 784             IS=IS+1
 785             IND(IS)=I
 786          ENDIF
 787   100 CONTINUE
 788       IP=IS
 789       DO 110 I=1,NPR
 790          IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 1) THEN
 791             IP=IP+1
 792             IND(IP)=I
 793          ENDIF
 794   110 CONTINUE
 795       DO 120 I=1,NPR
 796          IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 2) THEN
 797             IP=IP+1
 798             IND(IP)=I
 799          ENDIF
 800   120 CONTINUE
 801       DO 130 I=1,NPR
 802          IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 3) THEN
 803             IP=IP+1
 804             IND(IP)=I
 805          ENDIF
 806   130 CONTINUE
 807       DO 140 I=1,NC
 808          IN=I*ICD-ICD+1
 809          IF (IAM(IN,1) .EQ. 0) THEN
 810             ISC=ISC+1
 811             INDC(ISC)=I
 812          ENDIF
 813   140 CONTINUE
 814       IPC=ISC
 815       DO 150 I=1,NC
 816          IN=I*ICD-ICD+1
 817          IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 1) THEN
 818             IPC=IPC+1
 819             INDC(IPC)=I
 820          ENDIF
 821   150 CONTINUE
 822       DO 160 I=1,NC
 823          IN=I*ICD-ICD+1
 824          IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 2) THEN
 825             IPC=IPC+1
 826             INDC(IPC)=I
 827          ENDIF
 828   160 CONTINUE
 829       DO 170 I=1,NC
 830          IN=I*ICD-ICD+1
 831          IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 3) THEN
 832             IPC=IPC+1
 833             INDC(IPC)=I
 834          ENDIF
 835   170 CONTINUE
 836       DO 180 I=1,NPR
 837          TEMP(I)=CC(IND(I))
 838   180 CONTINUE
 839       DO 190 I=1,NPR
 840          CC(I)=TEMP(I)
 841   190 CONTINUE
 842       DO 200 I=1,NPR
 843          TEMP(I)=EX(IND(I))
 844   200 CONTINUE
 845       DO 210 I=1,NPR
 846          EX(I)=TEMP(I)
 847   210 CONTINUE
 848       DO 220 I=1,NPR
 849          TEMP(I)=CEN(IND(I),1)
 850   220 CONTINUE
 851       DO 230 I=1,NPR
 852          CEN(I,1)=TEMP(I)
 853   230 CONTINUE
 854       DO 240 I=1,NPR
 855          TEMP(I)=CEN(IND(I),2)
 856   240 CONTINUE
 857       DO 250 I=1,NPR
 858          CEN(I,2)=TEMP(I)
 859   250 CONTINUE
 860       DO 260 I=1,NPR
 861          TEMP(I)=CEN(IND(I),3)
 862   260 CONTINUE
 863       DO 270 I=1,NPR
 864          CEN(I,3)=TEMP(I)
 865   270 CONTINUE
 866       DO 280 I=1,NPR
 867          ITEMP(I)=IAM(IND(I),1)
 868   280 CONTINUE
 869       DO 290 I=1,NPR
 870          IAM(I,1)=ITEMP(I)
 871   290 CONTINUE
 872       DO 300 I=1,NPR
 873          ITEMP(I)=IAM(IND(I),2)
 874   300 CONTINUE
 875       DO 310 I=1,NPR
 876          IAM(I,2)=ITEMP(I)
 877   310 CONTINUE
 878 C     CALCULATE OVERLAP MATRIX OF STO-6G FUNCTIONS
 879 C
 880       DO 320 J=1,NC
 881          CALL OVLP(J,1,IS,IP,NPR,NC,ICD)
 882   320 CONTINUE
 883 C
 884       DO 330 J=1,NC
 885          DO 330 K=1,NC
 886             CESPM2(INDC(J),INDC(K))=OVL(J,K)
 887   330 CONTINUE
 888       DO 340 J=1,NC
 889          DO 340 K=1,NC
 890             OVL(J,K)=CESPM2(J,K)
 891   340 CONTINUE
 892       L=0
 893       DO 350 I=1,NC
 894          DO 350 J=1,I
 895             L=L+1
 896             CESP(L)=OVL(I,J)
 897   350 CONTINUE
 898 C
 899 C     DEORTHOGONALIZE THE COEFFICIENTS AND REFORM THE DENSITY MATRIX
 900 C
 901       CALL RSP(CESP,NC,1,TEMP,CESPML)
 902       DO 360 I=1,NC
 903          DO 360 J=1,I
 904             SUM=0.D0
 905             DO 360 K=1,NC
 906                SUM=SUM+CESPML(I+(K-1)*NC)/SQRT(TEMP(K))*CESPML(J+(K-1)*N
 907      1C)
 908                CESP(I+(J-1)*NC)=SUM
 909                CESP(J+(I-1)*NC)=SUM
 910   360 CONTINUE
 911       CALL MULT(C,CESP,CESPML,NC)
 912       CALL DENSIT(CESPML,NC,NC,NCLOSE,NOPEN,FRACT,CESP,2)
 913 C
 914 C     NOW CALCULATE THE ELECTRONIC CONTRIBUTION TO THE ELECTROSTATIC POT
 915 C
 916       L=0
 917       DO 370 I=1,NC
 918          DO 370 J=1,I
 919             L=L+1
 920             CESPM(I,J)=CESP(L)
 921             CESPM(J,I)=CESP(L)
 922   370 CONTINUE
 923       IPX=(NPR-IS)/3
 924       IPE=IS+IPX
 925       DO 380 I=1,NESP
 926          ES(I)=0.D0
 927   380 CONTINUE
 928       CALL NAICAS(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
 929       CALL NAICAP(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
 930 C     CALCULATE TOTAL ESP AND FORM ARRAYS FOR ESPFIT
 931       DO 400 I=1,NESP
 932          ESP(I)=0.D0
 933          DO 390 J=1,NATOM
 934             RA=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2+(CO(
 935      13,J)-POTPT(3,I))**2)
 936             ESP(I)=ESP(I)+TORE(IAN(J))/(RA/BOHR)
 937   390    CONTINUE
 938          ESP(I)=ESP(I)-ES(I)
 939          DO 400  J=1,NATOM
 940             RIJ=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2
 941      1+(CO(3,J)-POTPT(3,I))**2)/BOHR
 942             B(J)=B(J)+ESP(I)*1.D0/RIJ
 943   400 CONTINUE
 944 C
 945 C     IF REQUESTED WRITE OUT ELECTRIC POTENTIAL DATA TO
 946 C     UNIT 21
 947 C
 948       POTWRT=(INDEX(KEYWRD,'POTWRT') .NE. 0)
 949       IF(POTWRT) THEN
 950          OPEN(21,FILE='FOR021',STATUS='NEW')
 951          WRITE(21,'(I5)') NESP
 952          DO 410 I=1,NESP
 953   410    WRITE(21,420) ESP(I),POTPT(1,I)/BOHR,POTPT(2,I)/BOHR,
 954      1POTPT(3,I)/BOHR
 955       ENDIF
 956   420 FORMAT(1X,4E16.7)
 957       RETURN
 958       END
 959       DOUBLE PRECISION FUNCTION DEX2(M)
 960       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 961       IF(M .LT. 2) THEN
 962          DEX2=1
 963       ELSE
 964          DEX2=1
 965          DO 10 I=1,M,2
 966    10    DEX2=DEX2*I
 967       ENDIF
 968       RETURN
 969       END
 970       BLOCK DATA ESPBLO
 971       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 972       INCLUDE 'SIZES'
 973       COMMON /ESPC/  CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
 974      1                EX(MAXPR),ESPI(MAXORB,MAXORB),
 975      2                FV(0:8,821),FAC(0:7),
 976      3                DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR),
 977      4                OVL(MAXORB,MAXORB),FC(MAXPR*6)
 978       DATA TF/33.D0,37.D0,41.D0/
 979       DATA FAC/1.D0,1.D0,2.D0,6.D0,24.D0,120.D0,720.D0,5040.D0/
 980       END
 981 C***********************************************************************
 982       SUBROUTINE ESPFIT(IDIP,NATOM,NESP,IZ,ESP,POTPT,CO,
 983      1DX,DY,DZ,RMS,RRMS)
 984       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 985       INCLUDE 'SIZES'
 986 C***********************************************************************
 987 C
 988 C     THIS ROUTINE FITS THE ELECTROSTATIC POTENTIAL TO A MONOPOLE
 989 C     EXPANSION. FITTING TO THE DIPOLE MONENT CAN ALSO BE DONE.
 990 C     THIS ROUTINE WAS WRITTEN BY B.H.BESLER AND K.M.MERZ
 991 C     IN FEB. 1989 AT UCSF.
 992 C
 993 C     ON INPUT:  IDIP = FLAG TO INDICATE IF THE DIPOLE IS FIT
 994 C                NATOM = NUMBER OF ATOMS
 995 C                NESP = NUMBER OF ESP POINTS
 996 C                IZ = MOLECULAR CHARGE
 997 C                ESP = TOTAL ESP AT THE POINTS
 998 C                POTPT = ESP POINTS
 999 C                CO = COORDINATES
1000 C                DX = X COMPONENT OF THE DIPOLE
1001 C                DY = Y COMPONENT OF THE DIPOLE
1002 C                DZ = Z COMPONENT OF THE DIPOLE
1003 C
1004 C     ON OUTPUT: Q = ESP CHARGES
1005 C                RMS = ROOT MEAN SQUARE FIT
1006 C                RRMS = RELATIVE ROOT MEAN SQUARE FIT
1007 C
1008 C     FOR MORE DETAILS SEE: BESLER,MERZ,KOLLMAN J. COMPUT. CHEM.
1009 C     (IN PRESS)
1010 C***********************************************************************
1011       COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
1012      1Q(NUMATM+4),QSC(NUMATM+4),CF, ESPFD(MAXORB**2-NUMATM-5)
1013       DIMENSION CO(3,*),ESP(*),POTPT(3,*)
1014       BOHR = 0.529167D00
1015 C     CONVERSION FACTOR FOR DEBYE TO ATOMIC UNITS
1016       CF=5.2917715D-11*1.601917D-19/3.33564D-30
1017 C
1018 C     THE FOLLOWING SETS UP THE LINEAR EQUATION A*Q=B
1019 C     SET UP THE A(J,K) ARRAY
1020 C
1021       DO 20  K=1,NATOM
1022          DO 10  J=1,NATOM
1023             DO 10  I=1,NESP
1024                RIK=SQRT((CO(1,K)-POTPT(1,I))**2+(CO(2,K)-POTPT(2,I))**2
1025      1      +(CO(3,K)-POTPT(3,I))**2)/BOHR
1026                RIJ=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2
1027      1      +(CO(3,J)-POTPT(3,I))**2)/BOHR
1028                A(J,K)=A(J,K)+1.D0/RIK*1.D0/RIJ
1029    10    CONTINUE
1030          A(NATOM+1,K)=1.D0
1031          A(K,NATOM+1)=1.D0
1032          A(NATOM+1,NATOM+1)=0.D0
1033          IF(IDIP .EQ. 1) THEN
1034             A(NATOM+2,K)=CO(1,K)/BOHR
1035             A(K,NATOM+2)=CO(1,K)/BOHR
1036             A(NATOM+2,NATOM+2)=0.D0
1037             A(NATOM+3,K)=CO(2,K)/BOHR
1038             A(K,NATOM+3)=CO(2,K)/BOHR
1039             A(NATOM+3,NATOM+3)=0.D0
1040             A(NATOM+4,K)=CO(3,K)/BOHR
1041             A(K,NATOM+4)=CO(3,K)/BOHR
1042             A(NATOM+4,NATOM+4)=0.D0
1043          ENDIF
1044    20 CONTINUE
1045       B(NATOM+1)=FLOAT(IZ)
1046       B(NATOM+2)=DX/CF
1047       B(NATOM+3)=DY/CF
1048       B(NATOM+4)=DZ/CF
1049 C
1050 C     INSERT CHARGE AND DIPOLAR (IF DESIRED) CONSTRAINTS
1051 C
1052       IF(IDIP .EQ. 1) THEN
1053          L=0
1054          DO 30 I=1,NATOM+4
1055             DO 30 J=1,NATOM+4
1056                L=L+1
1057    30    AL(L)=A(I,J)
1058       ELSE
1059          L=0
1060          DO 40 I=1,NATOM+1
1061             DO 40 J=1,NATOM+1
1062                L=L+1
1063    40    AL(L)=A(I,J)
1064       ENDIF
1065       IF (IDIP .EQ. 1) THEN
1066          CALL OSINV(AL,NATOM+4,DET)
1067       ELSE
1068          CALL OSINV(AL,NATOM+1,DET)
1069       ENDIF
1070       IF(IDIP .EQ. 1) THEN
1071          L=0
1072          DO 50 I=1,NATOM+4
1073             DO 50 J=1,NATOM+4
1074                L=L+1
1075    50    A(I,J)=AL(L)
1076       ELSE
1077          L=0
1078          DO 60 I=1,NATOM+1
1079             DO 60 J=1,NATOM+1
1080                L=L+1
1081    60    A(I,J)=AL(L)
1082       ENDIF
1083 C
1084 C     SOLVE FOR THE CHARGES
1085 C
1086       IF(IDIP .EQ. 1) THEN
1087          DO 70 I=1,NATOM+4
1088             DO 70 J=1,NATOM+4
1089                Q(I)=Q(I)+A(I,J)*B(J)
1090    70    CONTINUE
1091       ELSE
1092          DO 80 I=1,NATOM+1
1093             DO 80 J=1,NATOM+1
1094                Q(I)=Q(I)+A(I,J)*B(J)
1095    80    CONTINUE
1096       ENDIF
1097 C
1098 C     CALCULATE ROOT MEAN SQUARE FITS AND RELATIVE ROOT MEAN SQUARE FITS
1099 C
1100       CTOT=0.0
1101       DO 100 I=1,NESP
1102          ESPC=0.D0
1103          DO 90 J=1,NATOM
1104             RIJ=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2
1105      1      +(CO(3,J)-POTPT(3,I))**2)/BOHR
1106    90    ESPC=ESPC+Q(J)/RIJ
1107          RMS=RMS+(ESPC-ESP(I))**2
1108   100 RRMS=RRMS+ESP(I)**2
1109       RMS=SQRT(RMS/NESP)
1110       RRMS=RMS/SQRT(RRMS/NESP)
1111       RMS=RMS*627.51D0
1112       RETURN
1113       END
1114 C***********************************************************************
1115       SUBROUTINE FSUB(N,X,FVAL)
1116       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1117 C***********************************************************************
1118 C
1119 C     CALCULATE THE FM(T). KINDLY SUPPLIED BY RUS PITZER AND CLEANED UP
1120 C     BY K.M.MERZ IN FEB. 1989 AT UCSF.
1121 C
1122 C     ON INPUT:  N = INDEX
1123 C                X = EXPONENT
1124 C     ON OUTPUT: FVAL = VALUE OF THE FUNCTION
1125 C
1126 C     FOR MORE DETAILS SEE: OBARA AND SAIKA J. CHEM. PHYS. 1986,84,3963
1127 C***********************************************************************
1128       DIMENSION FF(21),TERM(200),A(10),RT(10)
1129       DATA A0, A1S2, PIE4, A1
1130      1   /0.0D0,0.5D0,0.7853981633974483096156608D0,1.0D0/
1131       DATA XSW /24.0D0/
1132       E=A1S2*EXP(-X)
1133       FAC0=N
1134       FAC0=FAC0+A1S2
1135       IF(X.GT.XSW) GO TO 50
1136 C
1137 C     USE POWER SERIES
1138 C
1139    10 FAC=FAC0
1140       TERM0=E/FAC
1141       SUM=TERM0
1142       KU=(X-FAC0)
1143       IF(KU.LT.1) GO TO 30
1144 C
1145 C     SUM INCREASING TERMS FORWARDS
1146 C
1147       DO 20 K=1,KU
1148          FAC=FAC+A1
1149          TERM0=TERM0*X/FAC
1150          SUM=SUM+TERM0
1151    20 CONTINUE
1152    30 I=1
1153       FAC=FAC+A1
1154       TERM(1)=TERM0*X/FAC
1155       SUMA=SUM+TERM(1)
1156       IF(SUM.EQ.SUMA) GO TO 90
1157    40 I=I+1
1158       FAC=FAC+A1
1159       TERM(I)=TERM(I-1)*X/FAC
1160       SUM1=SUMA
1161       SUMA=SUMA+TERM(I)
1162       IF(SUM1-SUMA) 40,90,40
1163 C
1164 C     USE ASYMPTOTIC SERIES
1165 C
1166    50 SUM=SQRT(PIE4/X)
1167       IF(N.EQ.0) GO TO 70
1168       FAC=-A1S2
1169       DO 60 K=1,N
1170          FAC=FAC+A1
1171          SUM=SUM*FAC/X
1172    60 CONTINUE
1173    70 I=1
1174       TERM(1)=-E/X
1175       SUMA=SUM+TERM(1)
1176       IF(SUM.EQ.SUMA) GO TO 90
1177       FAC=FAC0
1178       KU=(X+FAC0-A1)
1179       DO 80 I=2,KU
1180          FAC=FAC-A1
1181          TERM(I)=TERM(I-1)*FAC/X
1182          SUM1=SUMA
1183          SUMA=SUMA+TERM(I)
1184          IF(SUM1.EQ.SUMA) GO TO 90
1185    80 CONTINUE
1186 C
1187 C     XSW SET TOO LOW. USE POWER SERIES.
1188 C
1189       GO TO 10
1190 C
1191 C     SUM DECREASING TERMS BACKWARDS
1192 C
1193    90 SUM1=A0
1194       DO 100 K=1,I
1195          SUM1=SUM1+TERM(I+1-K)
1196   100 CONTINUE
1197       FF(N+1)=SUM+SUM1
1198 C
1199 C     USE RECURRENCE RELATION
1200 C
1201       IF(N.EQ.0) GOTO 120
1202       DO 110 K=1,N
1203          FAC0=FAC0-A1
1204          FF(N+1-K)=(E+X*FF(N+2-K))/FAC0
1205   110 CONTINUE
1206   120 FVAL=FF(N+1)
1207       RETURN
1208       END
1209       SUBROUTINE SETUP3
1210       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1211       INCLUDE 'SIZES'
1212       COMMON /NATYPE/ NZTYPE(107),MTYPE(30),LTYPE
1213       COMMON /STO6G/ ALLC(6,5,2),ALLZ(6,5,2)
1214 C     SET-UP THE STEWART'S STO-3G EXPANSIONS
1215 C     FROM J. CHEM. PHYS. 52 431.
1216 C                                            1S
1217       ALLZ(1,1,1) =2.227660584D00
1218       ALLZ(2,1,1) =4.057711562D-01
1219       ALLZ(3,1,1) =1.098175104D-01
1220 C
1221       ALLC(1,1,1) =1.543289673D-01
1222       ALLC(2,1,1) =5.353281423D-01
1223       ALLC(3,1,1) =4.446345422D-01
1224 C                                      2S
1225       ALLZ(1,2,1) =2.581578398D00
1226       ALLZ(2,2,1) =1.567622104D-01
1227       ALLZ(3,2,1) =6.018332272D-02
1228 C
1229       ALLC(1,2,1) =-5.994474934D-02
1230       ALLC(2,2,1) =5.960385398D-01
1231       ALLC(3,2,1) =4.581786291D-01
1232 C                                     2P
1233       ALLZ(1,2,2) =9.192379002D-01
1234       ALLZ(2,2,2) =2.359194503D-01
1235       ALLZ(3,2,2) =8.009805746D-02
1236 C
1237       ALLC(1,2,2) =1.623948553D-01
1238       ALLC(2,2,2) =5.661708862D-01
1239       ALLC(3,2,2) =4.223071752D-01
1240 C                                      3S
1241       ALLZ(1,3,1) =5.641487709D-01
1242       ALLZ(2,3,1) =6.924421391D-02
1243       ALLZ(3,3,1) =3.269529097D-02
1244 C
1245       ALLC(1,3,1) =-1.782577972D-01
1246       ALLC(2,3,1) =8.612761663D-01
1247       ALLC(3,3,1) =2.261841969D-01
1248 C                                     3P
1249       ALLZ(1,3,2) =2.692880368D00
1250       ALLZ(2,3,2) =1.489359592D-01
1251       ALLZ(3,3,2) =5.739585040D-02
1252 C
1253       ALLC(1,3,2) =-1.061945788D-02
1254       ALLC(2,3,2) =5.218564264D-01
1255       ALLC(3,3,2) =5.450015143D-01
1256 C                                      4S
1257       ALLZ(1,4,1) =2.267938753D-01
1258       ALLZ(2,4,1) =4.448178019D-02
1259       ALLZ(3,4,1) =2.195294664D-02
1260 C
1261       ALLC(1,4,1) =-3.349048323D-01
1262       ALLC(2,4,1) =1.056744667D00
1263       ALLC(3,4,1) =1.256661680D-01
1264 C                                     4P
1265       ALLZ(1,4,2) =4.859692220D-01
1266       ALLZ(2,4,2) =7.430216918D-02
1267       ALLZ(3,4,2) =3.653340923D-02
1268 C
1269       ALLC(1,4,2) =-6.147823411D-02
1270       ALLC(2,4,2) =6.604172234D-01
1271       ALLC(3,4,2) =3.932639495D-01
1272 C                                      5S
1273       ALLZ(1,5,1) =1.080198458D-01
1274       ALLZ(2,5,1) =4.408119382D-02
1275       ALLZ(3,5,1) =2.610811810D-02
1276 C
1277       ALLC(1,5,1) =-6.617401158D-01
1278       ALLC(2,5,1) =7.467595004D-01
1279       ALLC(3,5,1) =7.146490945D-01
1280 C                                     5P
1281       ALLZ(1,5,2) =2.127482317D-01
1282       ALLZ(2,5,2) =4.729648620D-02
1283       ALLZ(3,5,2) =2.604865324D-02
1284 C
1285       ALLC(1,5,2) =-1.389529695D-01
1286       ALLC(2,5,2) =8.076691064D-01
1287       ALLC(3,5,2) =2.726029342D-01
1288 C
1289       RETURN
1290       END
1291       SUBROUTINE OVLP(IC,IESP,IS,IP,NPR,NC,ICD)
1292       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1293 C***********************************************************************
1294 C
1295 C     OVLP CALCULATES THE OVERLAP INTEGRALS FOR A STO-6G BASIS SET.
1296 C     THE RESULTING INTEGRALS ARE USED IN THE DEORTHOGONALIZATION
1297 C     PROCESS.
1298 C     THE CODE WAS WRITTEN BY B.H.BESLER AND K.M.MERZ IN FEB. 1989
1299 C     AT UCSF.
1300 C
1301 C     ON INPUT:  IC = LOOP INDEX
1302 C                IESP = LOOP INDEX
1303 C                IS = NUMBER OF S ORBITALS
1304 C                IP = NUMBER OF P ORBITALS
1305 C                NPR = NUMBER OF PRIMITIVES
1306 C                NC = NUMBER OF CONTRACTED FUNCTIONS
1307 C
1308 C     ON OUTPUT: OVL IS FILLED WITH THE OVERLAP INTEGRAL VALUE
1309 C
1310 C     FOR FURTHER INFO SEE: OBARA & SAIKA J.CHEM.PHYS. 1986,84,3963
1311 C***********************************************************************
1312       LOGICAL CALLED
1313       DOUBLE PRECISION NAI,NAI1,NAI2
1314       INCLUDE 'SIZES'
1315       COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
1316       COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
1317       COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
1318       COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
1319       COMMON /EXPONT/ ZS(107),ZP(107),ZD(107)
1320       COMMON /STO6G/  ALLC(6,5,2),ALLZ(6,5,2)
1321       COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
1322      1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821),
1323      2FAC(0:7),DEX(-1:96),TF(0:2),
1324      3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),XDMY(MAXPR*6)
1325       COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6),
1326      1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6),
1327      2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6)
1328      3,NAI1(MAXPR,6),NAI2(MAXPR,6)
1329       DATA BOHR/0.529167D0/
1330 C
1331 C     CALCULATE DISTANCE ARRAYS
1332 C
1333       PI=4.D0*ATAN(1.D0)
1334       IPR=IC*ICD-ICD+1
1335       ISTART=IPR
1336       DO 10 I=ISTART,NPR
1337          DX(I)=CEN(IPR,1)-CEN(I,1)
1338          DY(I)=CEN(IPR,2)-CEN(I,2)
1339          DZ(I)=CEN(IPR,3)-CEN(I,3)
1340          TD(I)=DX(I)**2+DY(I)**2+DZ(I)**2
1341    10 CONTINUE
1342 C
1343 C     CALCULATE EXPONENT SUM
1344 C
1345       DO 20 I=ISTART,NPR
1346          DO 20 J=1,ICD
1347             EXS(I,J)=1.D0/(EX(IPR+J-1)+EX(I))
1348             CE(I,J)=EX(IPR+J-1)*EX(I)*EXS(I,J)
1349    20 CONTINUE
1350 C
1351 C     CALCULATE EXPONENT WEIGHTED CENTERS
1352 C
1353       DO 30 I=ISTART,NPR
1354          DO 30 J=1,ICD
1355             EWCX(I,J)=(EX(I)*CEN(I,1)+EX(IPR+J-1)
1356      1*CEN(IPR+J-1,1))*EXS(I,J)
1357             EWCY(I,J)=(EX(I)*CEN(I,2)+EX(IPR+J-1)
1358      1*CEN(IPR+J-1,2))*EXS(I,J)
1359             EWCZ(I,J)=(EX(I)*CEN(I,3)+EX(IPR+J-1)
1360      1*CEN(IPR+J-1,3))*EXS(I,J)
1361    30 CONTINUE
1362       DO 40 I=1,NPR
1363          DO 40 J=1,ICD
1364             EXPN(I,J)=EXP(-CE(I,J)*TD(I))
1365             NAI(I,J)=(PI*EXS(I,J))**1.5D0*EXPN(I,J)
1366             EXPN(I,J)=NAI(I,J)
1367    40 CONTINUE
1368 C
1369 C     CALCULATE (S||P) ESP INTEGRALS
1370 C
1371       IF((IAM(IPR,1) .EQ. 0) .AND. (IS .NE. IP)) THEN
1372          NP=IS+1
1373          DO 80 I=NP,NPR
1374             DO 80 J=1,ICD
1375                GO TO (50,60,70),IAM(I,2)
1376    50          NAI(I,J)=(EWCX(I,J)-CEN(I,1))*EXPN(I,J)
1377                go TO 80
1378    60          NAI(I,J)=(EWCY(I,J)-CEN(I,2))*EXPN(I,J)
1379                GO TO 80
1380    70          NAI(I,J)=(EWCZ(I,J)-CEN(I,3))*EXPN(I,J)
1381    80    CONTINUE
1382       ENDIF
1383 C
1384 C     CALCULATE (P||S) ESP INTEGRALS
1385 C
1386       IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN
1387          NP=IS+1
1388          DO 120 I=      ISTART,NPR
1389             DO 120 J=1,ICD
1390                GO TO (90,100,110),IAM(IPR+J-1,2)
1391    90          NAI(I,J)=(EWCX(I,J)-CEN(IPR+J-1,1))*EXPN(I,J)
1392                GO TO 120
1393   100          NAI(I,J)=(EWCY(I,J)-CEN(IPR+J-1,2))*EXPN(I,J)
1394                GO TO 120
1395   110          NAI(I,J)=(EWCZ(I,J)-CEN(IPR+J-1,3))*EXPN(I,J)
1396   120    CONTINUE
1397       ENDIF
1398 C
1399 C     CALCULATE (P||P) ESP INTEGRALS
1400 C
1401       IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN
1402          DO 160 I=ISTART,NPR
1403             DO 160 J=1,ICD
1404                GO TO (130,140,150),IAM(I,2)
1405   130          NAI(I,J)=(EWCX(I,J)-CEN(I,1))*NAI(I,J)
1406                IF(IAM(IPR+J-1,2) .EQ. IAM(I,2))
1407      1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0
1408      2      *EXPN(I,J)
1409                GO TO 160
1410   140          NAI(I,J)=(EWCY(I,J)-CEN(I,2))*NAI(I,J)
1411                IF(IAM(IPR+J-1,2) .EQ. IAM(I,2))
1412      1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0
1413      2      *EXPN(I,J)
1414                GO TO 160
1415   150          NAI(I,J)=(EWCZ(I,J)-CEN(I,3))*NAI(I,J)
1416                IF(IAM(IPR+J-1,2) .EQ. IAM(I,2))
1417      1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0
1418      2      *EXPN(I,J)
1419   160    CONTINUE
1420       ENDIF
1421       IPS=IC*ICD-ICD+1
1422       DO 180 I=IC,NC
1423          JPS=I*ICD-ICD+1
1424          OVL(IC,I)=0.D0
1425          DO 170 J=JPS,JPS+ICD-1
1426             DO 170 K=IPS,IPS+ICD-1
1427                OVL(IC,I)=OVL(IC,I)+CC(J)*CC(K)*NAI(J,K-IPS+1)
1428   170    CONTINUE
1429          OVL(I,IC)=OVL(IC,I)
1430   180 CONTINUE
1431       RETURN
1432       END
1433       SUBROUTINE NAICAS(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
1434       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1435 C***********************************************************************
1436 C
1437 C     THIS SUBROUTINE EVALUATES (S|S) , (S|P) TYPE NUCLEAR ATTRACTION
1438 C     INTEGRALS FOR A STO-NG BASIS SET
1439 C     WRITTEN BY B.H. BESLER AT FORD SCIENTIFIC RESEARCH LABS IN
1440 C     DECEMBER 1989.
1441 C
1442 C     ON INPUT:  IC = LOOP INDEX OF THE GAUSSIAN
1443 C                IESP = LOOP INDEX OF THE ESP POINT
1444 C                IPE = INDEX OF LAST Px PRIMITIVE
1445 C                IPX = NUMBER OF Px PRIMITIVES
1446 C                IS = NUMBER OS S ORBITALS
1447 C                ISC = NUMBER OF CONTRACTED S ORBITALS
1448 C                IP = NUMBER OF P ORBITALS
1449 C                NPR = NUMBER OF PRIMITIVES
1450 C                NC = NUMBER OF CONTRACTED FUNCTIONS
1451 C
1452 C
1453 C     FOR MORE INFO SEE: OBARA&SAIKA J.CHEM.PHYS. 1986,84,3963.
1454 C***********************************************************************
1455       INCLUDE 'SIZES'
1456       DOUBLE PRECISION NAI,NAI1,NAI2
1457       CHARACTER*241 KEYWRD
1458       COMMON/KEYWRD/ KEYWRD
1459       COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
1460      1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
1461       COMMON /INDX/ INDC(MAXORB)
1462       COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
1463       COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
1464       COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
1465       COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
1466       COMMON /EXPONT/ ZS(107),ZP(107),ZD(107)
1467       COMMON /STO6G/  ALLC(6,5,2),ALLZ(6,5,2)
1468       COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
1469      1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821),
1470      2FAC(0:7),DEX(-1:96),TF(0:2),
1471      3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),EXSR(MAXPR,6)
1472       COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6),
1473      1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6),
1474      2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6)
1475      3,NAI1(MAXPR,6),NAI2(MAXPR,6)
1476       DATA BOHR/0.529167D0/
1477 C
1478 C     CALCULATE DISTANCE ARRAYS
1479 C
1480       WRITE(6,*)
1481       PI=4.D0*ATAN(1.D0)
1482       IPX2=2*IPX
1483 C     IF THIS IS A RESTART RUN, READ IN RESTART INFO
1484       IF(INDEX(KEYWRD,'ESPRST') .NE. 0) THEN
1485          OPEN(UNIT=15,FILE='ESP.DUMP',STATUS='OLD',FORM='UNFORMATTED')
1486          READ(15) JSTART,IESPS
1487          IF(JSTART .EQ. ISC*2) THEN
1488             CLOSE(15)
1489             RETURN
1490          ENDIF
1491          DO 10 I=1,NESP
1492             READ(15) ES(I)
1493    10    CONTINUE
1494          CLOSE(15)
1495 C
1496          JSTART=JSTART+1
1497       ELSE
1498          JSTART=1
1499       ENDIF
1500       NP=IS+1
1501       DO 200 IC=JSTART,ISC
1502          IPR=IC*ICD-ICD+1
1503          ISTART=IPR
1504          DO 20 I=ISTART,IPE
1505             DX(I)=CEN(IPR,1)-CEN(I,1)
1506             DY(I)=CEN(IPR,2)-CEN(I,2)
1507             DZ(I)=CEN(IPR,3)-CEN(I,3)
1508             TD(I)=DX(I)**2+DY(I)**2+DZ(I)**2
1509    20    CONTINUE
1510 C
1511 C     CALCULATE EXPONENT SUM
1512 C
1513          DO 30 I=ISTART,IPE
1514             DO 30 J=1,ICD
1515                EXSR(I,J)=EX(IPR+J-1)+EX(I)
1516                EXS(I,J)=1.D0/EXSR(I,J)
1517                CE(I,J)=EX(IPR+J-1)*EX(I)*EXS(I,J)
1518                EXPN(I,J)=EXP(-CE(I,J)*TD(I))
1519    30    CONTINUE
1520 C
1521 C     CALCULATE EXPONENT WEIGHTED CENTERS
1522 C
1523          DO 40 I=ISTART,IPE
1524             DO 40 J=1,ICD
1525                EWCX(I,J)=(EX(I)*CEN(I,1)+EX(IPR+J-1)
1526      1*CEN(IPR+J-1,1))*EXS(I,J)
1527                EWCY(I,J)=(EX(I)*CEN(I,2)+EX(IPR+J-1)
1528      1*CEN(IPR+J-1,2))*EXS(I,J)
1529                EWCZ(I,J)=(EX(I)*CEN(I,3)+EX(IPR+J-1)
1530      1*CEN(IPR+J-1,3))*EXS(I,J)
1531    40    CONTINUE
1532 C
1533 C     BEGIN LOOP OVER ESP POINTS
1534 C
1535          DO 180 IESP=1,NESP
1536             POTP1=POTPT(1,IESP)/BOHR
1537             POTP2=POTPT(2,IESP)/BOHR
1538             POTP3=POTPT(3,IESP)/BOHR
1539 C
1540 C     BEGIN LOOP OVER COMPONENTS OF CONTRACTED FUNCTION IC
1541 C
1542             DO 150 J=1,ICD
1543 C
1544 C     CALCULATE DISTANCE BETWEEN EXPONENT WEIGHTED AND PROBE POINT
1545 C
1546                DO 50 I=ISTART,IPE
1547                   U(I,J)=((EWCX(I,J)-POTP1)**2+(EWCY(I,J)-POTP2)**2+
1548      1      (EWCZ(I,J)-POTP3)**2)*EXSR(I,J)
1549                   NAI(I,J)=SQRT(PI/U(I,J))
1550    50          CONTINUE
1551 C
1552 C     CALCULATE ESP INTEGRALS
1553 C
1554                DO 70 I=ISTART,IPE
1555                   IF(U(I,J) .LE. TF(0)) THEN
1556                      IREF=DNINT(U(I,J)*20.D0)
1557                      REF=0.05D0*IREF
1558                      RES=U(I,J)-REF
1559                      TERM=1.D0
1560                      F0(I,J)=0.D0
1561                      DO 60 K=0,6
1562                         F=FV(K,IREF+1)
1563                         TS=F*TERM*FAC(K)
1564                         TERM=-TERM*RES
1565                         F0(I,J)=F0(I,J)+TS
1566    60                CONTINUE
1567                   ELSE
1568                      F0(I,J)=NAI(I,J)*0.5D0
1569                   ENDIF
1570    70          CONTINUE
1571                DO 90 I=NP,IPE
1572                   IF(U(I,J) .LE. TF(1)) THEN
1573                      IREF=DNINT(U(I,J)*20.D0)
1574                      REF=0.05D0*IREF
1575                      RES=U(I,J)-REF
1576                      TERM1=1.D0
1577                      F1(I,J)=0.D0
1578                      DO 80 K=0,6
1579                         FI=FV(K+1,IREF+1)
1580                         TS1=FI*TERM1*FAC(K)
1581                         TERM1=-TERM1*RES
1582                         F1(I,J)=F1(I,J)+TS1
1583    80                CONTINUE
1584                   ELSE
1585                      F1(I,J)=NAI(I,J)*0.25D0/U(I,J)
1586                   ENDIF
1587    90          CONTINUE
1588                DO 100 I=ISTART,IS
1589   100          U(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F0(I,J)
1590                NP=IS+1
1591                DO 110 I=NP,IPE
1592                   NAI(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F0(I,J)
1593                   NAI1(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F1(I,J)
1594   110          CONTINUE
1595 C
1596 C     CALCULATE (S||P) ESP INTEGRALS
1597 C
1598                IF((IAM(IPR,1) .EQ. 0) .AND. (IS .NE. IP)) THEN
1599                   DO 120 I=NP,IPE
1600   120             U(I,J)=(EWCX(I,J)-CEN(I,1))*NAI(I,J)
1601      1-(EWCX(I,J)-POTP1)*NAI1(I,J)
1602                   DO 130 I=IPE+1,IPE+1+IPX
1603   130             U(I,J)=(EWCY(I-IPX,J)-CEN(I-IPX,2))*NAI(I-IPX,J)
1604      1-(EWCY(I-IPX,J)-POTP2)*NAI1(I-IPX,J)
1605                   DO 140 I=IPE+1+IPX,NPR
1606   140             U(I,J)=(EWCZ(I-IPX2,J)-CEN(I-IPX2,3))*NAI(I-IPX2,J)
1607      1-(EWCZ(I-IPX2,J)-POTP3)*NAI1(I-IPX2,J)
1608                ENDIF
1609   150       CONTINUE
1610             IPS=IC*ICD-ICD+1
1611             DO 170 I=IC,NC
1612                JPS=I*ICD-ICD+1
1613                ESPI(I,IC)=0.D0
1614                DO 160 J=JPS,JPS+ICD-1
1615                   DO 160 K=IPS,IPS+ICD-1
1616                      ESPI(I,IC)=ESPI(I,IC)+CC(J)*CC(K)*U(J,K-IPS+1)
1617   160          CONTINUE
1618                ES(IESP)=ES(IESP)+2.D0*CESPM(INDC(I),INDC(IC))*ESPI(I,IC)
1619   170       CONTINUE
1620             ES(IESP)=ES(IESP)-CESPM(INDC(IC),INDC(IC))*ESPI(IC,IC)
1621   180    CONTINUE
1622 C     WRITE OUT RESTART INFORMATION
1623          OPEN(UNIT=15,FILE='ESP.DUMP',STATUS='UNKNOWN',FORM='UNFORMATTED
1624      1')
1625          IESPS=0
1626          WRITE(15) IC,IESPS
1627          DO 190 I=1,NESP
1628             WRITE(15) ES(I)
1629   190    CONTINUE
1630          CLOSE(15)
1631 C
1632          WRITE(6,'(A,F6.2,A)')
1633      1'NAICAS DUMPED: ',100.D0/ISC*IC,' PERCENT COMPLETE'
1634   200 CONTINUE
1635       RETURN
1636       END
1637       SUBROUTINE NAICAP(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
1638       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1639 C***********************************************************************
1640 C     THIS ROUTINE EVALUATES (P|P) NUCLEAR ATTRACTION INTEGRALS OVER
1641 C
1642 C     A STO-NG BASIS SET.
1643 C     WRITTEN BY B.H. BESLER AT FORD SCIENTIFIC RESEARCH LABS IN
1644 C     SEPT. 1989
1645 C
1646 C     ON INPUT:  IC = LOOP INDEX OF THE GAUSSIAN
1647 C                ICD = CONTRACTION DEPTH OF BASIS SET
1648 C                IESP = LOOP INDEX OF THE ESP POINT
1649 C                IS = NUMBER OS S PRIMITIVES
1650 C                IPE = INDEX OF LAST PX PRIMITIVE
1651 C                IPX = NUMBER OF PX PRIMITIVES
1652 C                IS = NUMBER OS S PRIMITIVES
1653 C                ISC = NUMBER OF CONTRACTED
1654 C                NPR = NUMBER OF PRIMITIVES
1655 C                NC = NUMBER OF CONTRACTED FUNCTIONS
1656 C
1657 C
1658 C     FOR MORE INFO SEE: OBARA&SAIKA J.CHEM.PHYS. 1986,84,3963.
1659 C***********************************************************************
1660       INCLUDE 'SIZES'
1661       DOUBLE PRECISION NAI,NAI1,NAI2
1662       CHARACTER*241 KEYWRD
1663       COMMON /KEYWRD/ KEYWRD
1664       COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
1665      1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
1666       COMMON /INDX/ INDC(MAXORB)
1667       COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
1668       COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
1669       COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
1670       COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
1671       COMMON /EXPONT/ ZS(107),ZP(107),ZD(107)
1672       COMMON /STO6G/  ALLC(6,5,2),ALLZ(6,5,2)
1673       COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
1674      1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821),
1675      2FAC(0:7),DEX(-1:96),TF(0:2),
1676      3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),EXSR(MAXPR,6)
1677       COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6),
1678      1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6),
1679      2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6)
1680      3,NAI1(MAXPR,6),NAI2(MAXPR,6)
1681       COMMON/FP/ PF0(MAXHES),PF1(MAXHES),PF2(MAXHES),ID(MAXPAR),
1682      1PEXS(MAXHES),PCE(MAXHES),PEXPN(MAXHES),PTD(MAXHES),
1683      2PEWCX(MAXHES),PEWCY(MAXHES),PEWCZ(MAXHES),IRD(MAXHES)
1684       DATA BOHR/0.529167D0/
1685 C     SET NUMBER OF EQUALLY SPACED DUMPS
1686       IDN=10
1687 C
1688       IDC=0
1689       WRITE(6,*)
1690       IPX2=2*IPX
1691       PI=4.D0*ATAN(1.D0)
1692       NP=IS+1
1693 C     SETUP INDEX ARRAY
1694       DO 10 I=NP,IPE
1695          IRD(I)=I-IS
1696          IRD(I+IPX)=I-IS
1697          IRD(I+IPX2)=I-IS
1698    10 CONTINUE
1699 C
1700 C     CALCULATE QUANTITIES INVARIANT WITH ESP POINT FOR
1701 C     (P|P) ESP INTEGRALS
1702 C
1703       IL=L
1704       L=0
1705       DO 30 I=NP,IPE
1706          DO 20 J=I,IPE
1707             L=L+1
1708             PTD(L)=(CEN(I,1)-CEN(J,1))**2+(CEN(I,2)-CEN(J,2))**2+
1709      1(CEN(I,3)-CEN(J,3))**2
1710             PEXS(L)=1.d0/(EX(I)+EX(J))
1711             PCE(L)=EX(I)*EX(J)*PEXS(L)
1712             PEXPN(L)=EXP(-PCE(L)*PTD(L))
1713             PEWCX(L)=(EX(I)*CEN(I,1)+EX(J)*CEN(J,1))*PEXS(L)
1714             PEWCY(L)=(EX(I)*CEN(I,2)+EX(J)*CEN(J,2))*PEXS(L)
1715             PEWCZ(L)=(EX(I)*CEN(I,3)+EX(J)*CEN(J,3))*PEXS(L)
1716    20    CONTINUE
1717 C
1718 C     SET UP OTHER INDEX ARRAY FOR PACKED SYMMETRIC ARRAY
1719 C     STORAGE
1720 C
1721          ID(I-IS)=L-IPX
1722    30 CONTINUE
1723 C
1724 C     READ IN RESTART INFORMATION IF THIS IS A RESTART
1725 C
1726       IF(INDEX(KEYWRD,'ESPRST') .NE. 0) THEN
1727          OPEN(UNIT=15,FILE='ESP.DUMP',STATUS='UNKNOWN',FORM='UNFORMATTED
1728      1')
1729          READ(15) JSTART,IESPS
1730          IF(JSTART .NE. ISC*2) THEN
1731             IESPS=0
1732             CLOSE(15)
1733             GOTO 50
1734          ENDIF
1735          DO 40 I=1,NESP
1736             READ(15) ES(I)
1737    40    CONTINUE
1738          CLOSE(15)
1739          IDC=FLOAT(IESPS)/FLOAT(NESP)*10
1740       ELSE
1741          IESPS=0
1742       ENDIF
1743    50 CONTINUE
1744 C
1745 C     LOOP OVER ESP PROBE POINTS
1746 C
1747       DO 250 IESP=IESPS+1,NESP
1748          POTP1=POTPT(1,IESP)/BOHR
1749          POTP2=POTPT(2,IESP)/BOHR
1750          POTP3=POTPT(3,IESP)/BOHR
1751 C     CALCULATE QUANTITY U
1752 C
1753          L=0
1754          DO 60 I=NP,IPE
1755             DO 60 J=I,IPE
1756                L=L+1
1757                PTD(L)=((PEWCX(L)-POTP1)**2+(PEWCY(L)-POTP2)**2+
1758      1      (PEWCZ(L)-POTP3)**2)/PEXS(L)
1759                PCE(L)=SQRT(PI/PTD(L))
1760    60    CONTINUE
1761 C
1762 C     CALCULATE F0, F1, AND F2(U) USING TAYLOR SERIES
1763 C     OR ASYMPTOTIC EXPANSION
1764 C
1765          IL=L
1766          L=0
1767          DO 100 I=1,IL
1768             IF(PTD(I) .LE. TF(0)) THEN
1769                IREF=DNINT(PTD(I)*20.D0)
1770                REF=0.05D0*IREF
1771                RES=PTD(I)-REF
1772                TERM=1.D0
1773                PF0(I)=0.D0
1774                DO 70 K=0,6
1775                   F=FV(K,IREF+1)
1776                   TS=F*TERM*FAC(K)
1777                   TERM=-TERM*RES
1778                   PF0(I)=PF0(I)+TS
1779    70          CONTINUE
1780             ELSE
1781                PF0(I)=PCE(I)*0.5D0
1782             ENDIF
1783             IF(PTD(I) .LE. TF(1)) THEN
1784                IREF=DNINT(PTD(I)*20.D0)
1785                REF=0.05D0*IREF
1786                RES=PTD(I)-REF
1787                TERM1=1.D0
1788                PF1(I)=0.D0
1789                DO 80 K=0,6
1790                   FI=FV(K+1,IREF+1)
1791                   TS1=FI*TERM1*FAC(K)
1792                   TERM1=-TERM1*RES
1793                   PF1(I)=PF1(I)+TS1
1794    80          CONTINUE
1795             ELSE
1796                PF1(I)=PCE(I)*0.25D0/PTD(I)
1797             ENDIF
1798             IF(PTD(I) .LE. TF(2)) THEN
1799                IREF=DNINT(PTD(I)*20.D0)
1800                REF=0.05D0*IREF
1801                RES=PTD(I)-REF
1802                TERM2=1.D0
1803                PF2(I)=0.D0
1804                DO 90 K=0,6
1805                   FII=FV(K+2,IREF+1)
1806                   TS2=FII*TERM2*FAC(K)
1807                   TERM2=-TERM2*RES
1808                   PF2(I)=PF2(I)+TS2
1809    90          CONTINUE
1810             ELSE
1811                PF2(I)=PCE(I)*0.375D0/(PTD(I)*PTD(I))
1812             ENDIF
1813   100    CONTINUE
1814 C
1815 C     CALCULATE (S||S) TYPE INTEGRALS
1816 C
1817          DO 110 I=1,IL
1818             PF0(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF0(I)
1819             PTD(I)=PF0(I)
1820             PF1(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF1(I)
1821             PF2(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF2(I)
1822   110    CONTINUE
1823 C
1824          DO 230 IC=ISC+1,NC
1825             IPR=IC*ICD-ICD+1
1826             ISTART=IPR
1827             DO 200 J=1,ICD
1828 C
1829 C     CALCULATE (P||S) ESP INTEGRALS
1830 C
1831                IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN
1832                   DO 150 I=ISTART,NPR
1833                      IN=IPR+J-1
1834                      IR=IRD(I)+ID(IRD(IN))
1835                      IR2=ID(IRD(I))+IRD(IN)
1836                      IF(IR2 .LE. IR ) IR=IR2
1837                      GO TO (120,130,140),IAM(IN,2)
1838   120                NAI2(I,J)=(PEWCX(IR)-CEN(IN,1))*PF1(IR)-PF2(IR)*
1839      1      (PEWCX(IR)-POTP1)
1840                      NAI(I,J)=(PEWCX(IR)-CEN(IN,1))*PF0(IR)-PF1(IR)*
1841      1      (PEWCX(IR)-POTP1)
1842                      GO TO 150
1843   130                NAI2(I,J)=(PEWCY(IR)-CEN(IN,2))*PF1(IR)-PF2(IR)*
1844      1      (PEWCY(IR)-POTP2)
1845                      NAI(I,J)=(PEWCY(IR)-CEN(IN,2))*PF0(IR)-PF1(IR)*
1846      1      (PEWCY(IR)-POTP2)
1847                      GO TO 150
1848   140                NAI2(I,J)=(PEWCZ(IR)-CEN(IN,3))*PF1(IR)-PF2(IR)*
1849      1      (PEWCZ(IR)-POTP3)
1850                      NAI(I,J)=(PEWCZ(IR)-CEN(IN,3))*PF0(IR)-PF1(IR)*
1851      1      (PEWCZ(IR)-POTP3)
1852   150             CONTINUE
1853                ENDIF
1854 C
1855 C     CALCULATE (P||P) ESP INTEGRALS
1856 C
1857                IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN
1858                   DO 190 I=ISTART,NPR
1859                      IN=IPR+J-1
1860                      IR=IRD(I)+ID(IRD(IN))
1861                      IR2=ID(IRD(I))+IRD(IN)
1862                      IF(IR2 .LE. IR ) IR=IR2
1863                      GO TO (160,170,180),IAM(I,2)
1864   160                NAI(I,J)=(PEWCX(IR)-CEN(I,1))*NAI(I,J)-(PEWCX(IR)-P
1865      1OTP1)*      NAI2(I,J)
1866                      IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS(
1867      1IR)*      0.5D0*(PTD(IR)-PF1(IR))
1868                      GO TO 190
1869   170                NAI(I,J)=(PEWCY(IR)-CEN(I,2))*NAI(I,J)-(PEWCY(IR)-P
1870      1OTP2)*      NAI2(I,J)
1871                      IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS(
1872      1IR)*      0.5D0*(PTD(IR)-PF1(IR))
1873                      GO TO 190
1874   180                NAI(I,J)=(PEWCZ(IR)-CEN(I,3))*NAI(I,J)-(PEWCZ(IR)-P
1875      1OTP3)*      NAI2(I,J)
1876                      IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS(
1877      1IR)*      0.5D0*(PTD(IR)-PF1(IR))
1878   190             CONTINUE
1879                ENDIF
1880   200       CONTINUE
1881 C
1882 C     FORM INTEGRALS OVER CONTRACTED FUNCTIONS
1883 C
1884             IPS=IC*ICD-ICD+1
1885             DO 220 I=IC,NC
1886                JPS=I*ICD-ICD+1
1887                ESPI(I,IC)=0.D0
1888                DO 210 J=JPS,JPS+ICD-1
1889                   DO 210 K=IPS,IPS+ICD-1
1890                      ESPI(I,IC)=ESPI(I,IC)+CC(J)*CC(K)*NAI(J,K-IPS+1)
1891   210          CONTINUE
1892                ES(IESP)=ES(IESP)+2.D0*CESPM(INDC(I),INDC(IC))*ESPI(I,IC)
1893   220       CONTINUE
1894             ES(IESP)=ES(IESP)-CESPM(INDC(IC),INDC(IC))*ESPI(IC,IC)
1895   230    CONTINUE
1896 C
1897 C     WRITE OUT RESTART INFORMATION EVERY NESP/10 POINTS
1898 C
1899          IF(MOD(IESP,NESP/IDN) .EQ. 0) THEN
1900             OPEN(UNIT=15,FILE='ESP.DUMP',STATUS='UNKNOWN',FORM='UNFORMAT
1901      1TED')
1902             JSTART=ISC*2
1903             WRITE(15) JSTART,IESP
1904             DO 240 I=1,NESP
1905                WRITE(15) ES(I)
1906   240       CONTINUE
1907             CLOSE(15)
1908             IDC=IDC+1
1909             WRITE(6,'(A,F6.2,A)')
1910      1'NAICAP DUMPED: ',100.D0/IDN*IDC,' PERCENT COMPLETE'
1911          ENDIF
1912   250 CONTINUE
1913       RETURN
1914       END