From breneman@vax1.chem.rpi.edu Mon Feb 11 09:52:13 1991 Received: from quant.chem.rpi.edu by rpi.edu (4.1/RPI-ITS-SM55); id AA04379; Mon, 11 Feb 91 09:47:46 EST for bbesler@vela.acs.oakland.edu Received: from VAX1.CHEM.RPI.EDU by quant.chem.rpi.edu (4.0/ST22); id AA04057; Mon, 11 Feb 91 09:41:34 EST for taylor@vax1.chem.rpi.edu Date: Mon, 11 Feb 1991 09:45:20 EST From: breneman@vax1.chem.rpi.edu To: sybyl@quant.chem.rpi.edu Message-Id: <0094410d.9918d9a0.2596@VAX1.CHEM.RPI.EDU> Subject: SYBYL Users Community Bulletin Board System. Status: OR Dear Charter Members of the SYBYL BBS, This is the inaugural message from the recently-installed SYBYL Users Community BBS. You may distribute messages to the registered recipients of the SYBYL BBS by sending them to sybyl@quant.chem.rpi.edu. Please encourage others with an interest in SYBYL to join our group. Messages concerning the operation of the BBS and new member registrations should be sent to breneman@quant.chem.rpi.edu. Happy Modeling, Curt Breneman Asst. Professor of Chemistry Rensselaer Polytechnic Institute Troy, NY 12180 (518) 276-2678 From chemistry-request@ccl.net Fri May 3 19:21:55 1991 Return-Path: Received: by oscsunb.ccl.net (5.64+IDA+kva.1/901126.13) id AA28929; Fri, 3 May 91 17:40:08 -0400 Received: from rpi.edu by oscsunb.ccl.net (5.64+IDA+kva.1/901126.13) id AA28921; Fri, 3 May 91 17:40:02 -0400 From: breneman@quant.chem.rpi.edu Received: from quant.chem.rpi.edu by rpi.edu (4.1/RPI-ITS-SM55); id AA27594; Fri, 3 May 91 17:41:26 EDT for chemistry@ccl.net Received: by quant.chem.rpi.edu (4.0/ST22); id AA08744; Fri, 3 May 91 17:39:54 EDT for chemistry@ccl.net Date: Fri, 3 May 91 17:39:54 EDT Message-Id: <9105032139.AA08744@quant.chem.rpi.edu> Apparently-To: chemistry@ccl.net Sender: chemistry-request@ccl.net Errors-To: owner-chemistry@ccl.net Precedence: bulk Status: OR Folks, For computing electrostatic potential-derived charges from ab-initio wavefunctions generated by one of the Gaussian 86/88/90 packages, the CHELPG program is helpful. This program is a rotationally-invariant modification of the QCPE program CHELP by Chirlian and Francl. For a description of the CHELPG program and some results, see: Breneman, C. and Wiberg, K, J. Comput. Chem. 1990, 11, 361. The program source code is available for Unix and VMS via internet mail. Curt Breneman breneman@quant.chem.rpi.edu From chemistry-request@ccl.net Sun May 12 17:17:20 1991 Return-Path: Received: by oscsunb.ccl.net (5.64+IDA+kva.1/901126.13) id AA24121; Sun, 12 May 91 15:55:32 -0400 Received: from rpi.edu by oscsunb.ccl.net (5.64+IDA+kva.1/901126.13) id AA24113; Sun, 12 May 91 15:55:20 -0400 From: breneman@quant.chem.rpi.edu Received: from quant.chem.rpi.edu by rpi.edu (4.1/RPI-ITS-SM55); id AA00810; Sun, 12 May 91 15:56:54 EDT for chemistry@ccl.net Received: by quant.chem.rpi.edu (4.0/ST22); id AA13534; Sun, 12 May 91 15:54:18 EDT for chemistry@ccl.net Date: Sun, 12 May 91 15:54:18 EDT Message-Id: <9105121954.AA13534@quant.chem.rpi.edu> Apparently-To: chemistry@ccl.net Sender: chemistry-request@ccl.net Errors-To: owner-chemistry@ccl.net Precedence: bulk Status: OR To all those who requested the CHELPG source code: It is now on its way to you via Internet mail. You will need to compile and link it with the Gaussian utility library, and with the l401 library. Please contact me if there are difficulties. Curt Breneman breneman@quant.chem.rpi.edu (518) 276-2678 (Thanks for the great response.) From breneman@quant.chem.rpi.edu Sun May 12 21:45:58 1991 Received: from quant.chem.rpi.edu by rpi.edu (4.1/RPI-ITS-SM55); id AA02467; Sun, 12 May 91 21:45:53 EDT for bbesler@vela.acs.oakland.edu Received: by quant.chem.rpi.edu (4.0/ST22); id AA13580; Sun, 12 May 91 21:40:34 EDT for bbesler@vela.acs.oakland.edu Date: Sun, 12 May 91 21:40:34 EDT Message-Id: <9105130140.AA13580@quant.chem.rpi.edu> Apparently-To: bbesler@vela.acs.oakland.edu Status: OR Brent, Here comes the CHELPG source code. The following three messages are: 1) the CHELPG source code in generic Unix fortran77; 2) A sample h2o input file; 3) A sample h2o output file. The program should work with G88 or G90 in its present form. Merely compile it, and then link with l401.a and util.a. You may need to run c8690 on your checkpoint files if they were written with an early version of G88. Let me know if I can be of any assitance. Curt Breneman RPI Chemistry breneman@quant.chem.rpi.edu From breneman@quant.chem.rpi.edu Sun May 12 21:46:21 1991 Received: from quant.chem.rpi.edu by rpi.edu (4.1/RPI-ITS-SM55); id AA02474; Sun, 12 May 91 21:46:11 EDT for bbesler@vela.acs.oakland.edu Received: by quant.chem.rpi.edu (4.0/ST22); id AA13586; Sun, 12 May 91 21:44:16 EDT for bbesler@vela.acs.oakland.edu Date: Sun, 12 May 91 21:44:16 EDT Message-Id: <9105130144.AA13586@quant.chem.rpi.edu> Apparently-To: bbesler@vela.acs.oakland.edu Status: OR h2o.chk Water, STO-3G STO-3G G88/90 ChelpG TEST 1 0 0 0 0 6 1.7,1.45,1.45 0 1. 2.8,1.2 From breneman@quant.chem.rpi.edu Sun May 12 21:46:30 1991 Received: from quant.chem.rpi.edu by rpi.edu (4.1/RPI-ITS-SM55); id AA02479; Sun, 12 May 91 21:46:20 EDT for bbesler@vela.acs.oakland.edu Received: by quant.chem.rpi.edu (4.0/ST22); id AA13589; Sun, 12 May 91 21:44:24 EDT for bbesler@vela.acs.oakland.edu Date: Sun, 12 May 91 21:44:24 EDT Message-Id: <9105130144.AA13589@quant.chem.rpi.edu> Apparently-To: bbesler@vela.acs.oakland.edu Status: O Warning: nonstandard floating-point arithmetic mode was enabled in this program and was never cleared. COORDINATES 0.000000 0.000000 0.240356 0.000000 1.432694 -0.961424 0.000000 -1.432694 -0.961424 NATOMS = 3 ICHARG = 0 MULTIP = 1 NAE = 5 NBE = 5 NE = 10 NBASIS = 7 NSHELL = 4 MAXTYP = 1 IAN 8 1 1 CENTER TYPE SHELLA 1 0 1 1 1 4 2 0 7 3 0 10 0 EXPON EXPCOF(S) EXPCOF(P) 0.130709321E+03 0.425194328E+01 0.000000000E+00 0.238088661E+02 0.411229442E+01 0.000000000E+00 0.644360831E+01 0.128162255E+01 0.000000000E+00 0.503315132E+01 -0.239413005E+00 0.167545020E+01 EXPCOF(D) EXPCOF(F) 0.000000000E+00 0.000000000E+00 0.000000000E+00 0.000000000E+00 0.000000000E+00 0.000000000E+00 0.000000000E+00 0.000000000E+00 ATOM # V.D.W. RADII (MULTIPLIED BY FACTOR) 1 1.70 2 1.45 3 1.45 Wavefunction is RHF. RMAX = 2.8000000000000 (ANGS), DELR = 0.30000000000000 (ANGS). THERE ARE 3 ATOMS TO CONSIDER. XMAX = 0. (AU), XMIN = 0. (AU). YMAX = 1.4326936135682 (AU), YMIN = -1.4326936135682 (AU). ZMAX = 0.24035589486613 (AU), ZMIN = -0.96142357946451 (AU). NUMBER OF X POINTS REQUIRED = 18 NUMBER OF Y POINTS REQUIRED = 23 NUMBER OF Z POINTS REQUIRED = 20 TOTAL NUMBER OF POINTS CONSIDERED = 8280 NUMBER OF POINTS SELECTED FOR FITTING : 3975 CHARGES FROM ELECTROSTATIC POTENTIAL GRID CHELPGrid Grid Modification - C. Breneman, Yale Chem, 1989 Water, STO-3G STO-3G G88/90 ChelpG TEST CHECKPOINT FILE: h2o.chk 0- 0- 0 MOLECULAR GEOMETRY ATOMIC NUMBER X Y Z 8 0.0000000 0.0000000 0.2403559 1 0.0000000 1.4326936 -0.9614236 1 0.0000000 -1.4326936 -0.9614236 THE TOTAL CHARGE IS CONSTRAINED TO: 0 NET CHARGES ATOMIC NUMBER CHARGE 8 -0.5637 1 0.2818 1 0.2818 THE DIPOLE MOMENT OF THESE CHARGES IS: 1.72184 FIT TO ELECTROSTATIC POTENTIAL AT 3975 POINTS ROOT MEAN SQUARE DEVIATION IS 0.0087 KCAL/MOLE 172.6u 0.9s 6:54 41% 0+1240k 8+1io 52pf+0w From breneman@quant.chem.rpi.edu Sun May 12 21:46:57 1991 Received: from quant.chem.rpi.edu by rpi.edu (4.1/RPI-ITS-SM55); id AA02473; Sun, 12 May 91 21:46:03 EDT for bbesler@vela.acs.oakland.edu Received: by quant.chem.rpi.edu (4.0/ST22); id AA13583; Sun, 12 May 91 21:44:06 EDT for bbesler@vela.acs.oakland.edu Date: Sun, 12 May 91 21:44:06 EDT Message-Id: <9105130144.AA13583@quant.chem.rpi.edu> Apparently-To: bbesler@vela.acs.oakland.edu Status: O C CHELP-NET ATOMIC CHARGES FROM AB INITIO WAVE FUNCTIONS C Modified for Grid Operations by Curt Breneman, Yale University C Department of Chemistry, 3/88 (Currently of Rensselaer C Polytechnic Institute, Troy, NY 12180.) C C CHELPG C C (NET ATOMIC) CHARGES FIT TO ELECTROSTATIC POTENTIALS C C Original CHELP code by: C M.M. FRANCL C L.E. CHIRLIAN C C OCTOBER 1985 C PRINCETON CHEMISTRY DEPARTMENT VAX 11/780 C VMS 3.7 C C FEBRUARY 1988 C MODIFIED TO USE GAUSSIAN86 CHECKPOINT FILES C Modified to use G88/90 checkpoint files 1/89 C YALE UNIVERSITY DEPARTMENT OF CHEMISTRY C WIBERG GROUP VMS 4.5 C CURT BRENEMAN C IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4 HANDLE1 HANDLE1=0 C*** TRACE-7 C ISTAT1=LIB$INIT_TIMER(HANDLE1) C*** C C READ IN DATA FROM CHECKPOINT FILE C CALL READIN C C C SELECT POINTS FOR FITTING, BEGIN WITH SHELL OF RADIUS 2A AND C INCREASING BY .5A SELECT POINTS FROM THE ROUGHLY RADIAL DISTRIBUTION C WHICH ARE NOT ENCLOSED BY THE VAN DER WAALS ENVELOPE OF THE MOLECULE C UNTIL A PREDETERMINED MAXIMUM NUMBER OF POINTS HAS BEEN REACHED C CALL BALL C C CALCULATE THE ELECTROSTATIC POTENTIAL USING FIRST ORDER HARTREE-FOCK C PERTURBATION THEORY C CALL EP C C USING METHOD OF LAGRANGE MULTIPLIERS, FIT BY LEAST SQUARES THE CHARGES C TO THE ELECTROSTATIC POTENTIAL, CONSTRAINING THE FIT TO REPRODUCE THE C TOTAL MOLECULAR CHARGE C CALL FIT C C PRINT OUT TABLE OF RESULTS C CALL OUTPUT C C*** TRACE-7 C ISTAT1=LIB$SHOW_TIMER(HANDLE1) C*** END C C SUBROUTINE BALL C C ROUTINE TO SELECT POINTS FOR FITTING TO THE ELECTROSTATIC POTENTIAL. C C POINTS WHICH LIE WITHIN THE VAN DER WAALS ENVELOPE OF THE MOLECULE C ARE EXCLUDED. C C POINTS ARE INITIALLY SELECTED IN A CUBE AROUND THE MOLECULE WHICH C IS SCALED TO THE SIZE OF THE MOLECULE+RMAX. THIS IS PRESENTLY AN INPUT C PARAMETER. POINTS ARE THEN EXCLUDED IF THEY FALL WITHIN THE INPUT C VDW RADIUS OF ANY OF THE ATOMS, OR, IF THEY FALL OUTSIDE C A DESIGNATED DISTANCE (RMAX) FROM ALL OF THE ATOMS. THE REMAINING C POINTS ARE PACKED IN A SET OF THREE (X,Y,Z) VECTORS, AND SENT TO THE C LAGRANGE LEAST-SQUARES FITTING ROUTINE. THE ORIGINAL CHELP INPUT C DECK IS AUGMENTED BY ADDING TWO FREE-FORMAT VARIABLES AT THE END. C THE TWO NEW INPUT VARIABLES ARE 'RMAX' AND 'DELR', WHERE RMAX C IS THE MAXIMUM DISTANCE A POINT CAN BE FROM ANY ATOM AND STILL C BE CONSIDERED IN THE FIT, AND DELR IS THE DISTANCE BETWEEN POINTS C IN THE GRID. BOTH RMAX AND DELR ARE IN ANGSTROMS. C C CURT BRENEMAN AND TERESA LEPAGE C YALE UNIVERSITY DEPARTMENT OF CHEMISTRY 3/88 C C ORIGINAL CODE BY: C C L.E. CHIRLIAN C M.M. FRANCL C APRIL 1985 C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (NPOINTS = 50000) COMMON /IO/ IN,IOUT C+++ COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS, $ IAN(401),ATMCHG(400),C(3,400) C+++ COMMON /IPO/ IPO(5) COMMON /SPHERE/ RADII(400),NTOTP COMMON /POINTS/ P(3,NPOINTS), MAXPNTS C DATA ANG2AU /1.889726878D0/ C C*** READ IN THE THE RMAX AND DELR VALUES IN ANGSTROMS. C read(IN,*) RMAX, DELR write(IOUT,*) ' RMAX = ',RMAX,' (ANGS), DELR = ',DELR,' (ANGS).' C*** C C CONVERT RADII TO AU C DELR = DELR * ANG2AU RMAX = RMAX * ANG2AU C C WHILE CONVERTING THE VDW RADII TO AU, FIND THE EXTREMA OF THE C MOLECULAR GEOMETRY. C XMAX=-50.0D0 XMIN=50.0D0 YMAX=-50.0D0 YMIN=50.0D0 ZMAX=-50.0D0 ZMIN=50.0D0 C WRITE(IOUT,*) ' THERE ARE ',NATOMS,' ATOMS TO CONSIDER.' DO 10 I=1,NATOMS RADII(I) = RADII(I) * ANG2AU C IF (C(1,I) .GT. XMAX) XMAX = C(1,I) IF (C(1,I) .LT. XMIN) XMIN = C(1,I) IF (C(2,I) .GT. YMAX) YMAX = C(2,I) IF (C(2,I) .LT. YMIN) YMIN = C(2,I) IF (C(3,I) .GT. ZMAX) ZMAX = C(3,I) IF (C(3,I) .LT. ZMIN) ZMIN = C(3,I) 10 CONTINUE C WRITE(IOUT,*) ' XMAX = ',XMAX,' (AU), XMIN = ',XMIN,' (AU).' WRITE(IOUT,*) ' YMAX = ',YMAX,' (AU), YMIN = ',YMIN,' (AU).' WRITE(IOUT,*) ' ZMAX = ',ZMAX,' (AU), ZMIN = ',ZMIN,' (AU).' C C DETERMINE THE MINIMUM CUBE DIMENSIONS REQUIRED TO CONTAIN C THE MOLECULE, INCLUDING THE MAXIMUM SELECTION RADIUS (RMAX) C ON BOTH SIDES. C XRANGE = XMAX - XMIN + 2.0D0 * RMAX YRANGE = YMAX - YMIN + 2.0D0 * RMAX ZRANGE = ZMAX - ZMIN + 2.0D0 * RMAX C NXPTS = INT(XRANGE/DELR) NYPTS = INT(YRANGE/DELR) NZPTS = INT(ZRANGE/DELR) C WRITE(IOUT,*) ' NUMBER OF X POINTS REQUIRED = ',NXPTS WRITE(IOUT,*) ' NUMBER OF Y POINTS REQUIRED = ',NYPTS WRITE(IOUT,*) ' NUMBER OF Z POINTS REQUIRED = ',NZPTS MAXPOSS = NXPTS * NYPTS * NZPTS WRITE(IOUT,*) ' TOTAL NUMBER OF POINTS CONSIDERED = ',MAXPOSS C C C RESET POINT COUNTER FOR NUMBER OF SELECTED POINTS C IPOINT = 0 C C LOOP OVER POSSIBLE POINTS C DO 200 II = 1,NXPTS + 1 C P1 = XMIN - RMAX + DBLE(II-1) * DELR C DO 200 JJ = 1,NYPTS + 1 C P2 = YMIN - RMAX + DBLE(JJ-1) * DELR C DO 200 KK = 1,NZPTS + 1 C P3 = ZMIN - RMAX + DBLE(KK-1) * DELR C C C IS THIS POINT WITHIN A VAN DER WAALS SPHERE OR OUTSIDE THE C RMAX DISTANCE FROM ALL ATOMS? C RADMIN=50.0D0 DO 100 I=1,NATOMS VRAD = RADII(I) DIST = (P1 - C(1,I))**2 + (P2 - C(2,I))**2 + (P3 - C(3,I))**2 DIST = DSQRT(DIST) IF (DIST .LT. VRAD) GOTO 210 IF (DIST .LT. RADMIN) RADMIN = DIST 100 CONTINUE IF (RADMIN .GT. RMAX) GOTO 210 C C STORE POINTS (IN ATOMIC UNITS) C IPOINT = IPOINT + 1 P(1,IPOINT) = P1 P(2,IPOINT) = P2 P(3,IPOINT) = P3 IF (IPO(2) .EQ. 1) $ WRITE(IOUT,*) 'POINT ',IPOINT,' X,Y,Z ',P1,P2,P3 210 CONTINUE 200 CONTINUE C MAXPNTS = IPOINT WRITE(IOUT,*) ' NUMBER OF POINTS SELECTED FOR FITTING : ',MAXPNTS RETURN END C SUBROUTINE EP C C ROUTINE TO CALCULATE THE ELECTROSTATIC POTENTIAL FROM FIRST ORDER C PERTURBATION THEORY C C M.M. FRANCL APRIL 1985 C MODIFIED VERSION OF A MEPHISTO ROUTINE C RESTRICTED TO CLOSED SHELL MOLECULES C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NPOINTS = 50000) INTEGER*4 SHELLA,SHELLN,SHELLT,AOS,SHELLC,AON,HANDLE C COMMON /IO/ IN,IOUT COMMON /IPO/ IPO(5) C+++ COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS, $ IAN(401),ATMCHG(400),C(3,400) C C=== Gaussian88 Modification for enlarged common /b/. Common/B/EXX(6000),C1(6000),C2(6000),C3(6000),X(2000),Y(2000), $Z(2000),JAN(2000),ShellA(2000),ShellN(2000),ShellT(2000), $ShellC(2000),AOS(2000),AON(2000),NShell,MaxTyp C==== Old G86 Version of common /b/ c COMMON/B/EXX(1200),C1(1200),C2(1200),C3(1200), c $ X(400),Y(400),Z(400),JAN(400),SHELLA(400),SHELLN(400), c $ SHELLT(400),SHELLC(400),AOS(400),AON(400),NSHELL,MAXTYP C+++ C COMMON /B/ EXX(240),C1(240),C2(240),C3(240),X(80),Y(80),Z(80), C $ JAN(80),SHELLA(80),SHELLN(80),SHELLT(80),SHELLC(80) C $ ,AOS(80),AON(80),NSHELL,MAXTYP C COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101), C $ ATMCHG(100),C(3,100) COMMON /POINTS/ P(3,NPOINTS),MAXPNTS COMMON /ELP/ ELECP(NPOINTS) COMMON /CHARGE/ COEF_ALPHA(100000),COEF_BETA(100000),IUHF COMMON /OUT/ NTITLE(20,3),CHKFIL,Q(400),I6TO5,RMS,PERCENT, 1 NLIN,NEND(3) C DIMENSION HPERT(100000),INDEX(1280) C DATA IPTCHG/1.0/ DATA ZERO/0.0/, TWO/2.0/, VNUCMAX/30.0/ C DIVERT TO ROUTINE UEP IF WAVEFUNCTION IS UNRESTRICTED C HARTREE-FOCK WAVEFUNCTION C IF (IUHF .EQ. 1) THEN CALL UEP RETURN END IF C HANDLE = 0 C C SET UP THE INDEXING TABLE FOR HPERT C DO 100 I=1,NBASIS INDEX(I) = (I-1)*I/2 100 CONTINUE C C BEGIN LOOP TO CALCULATE ELECTROSTATIC POTENTIAL C NOCC = NEL / 2 MVIR = NOCC + 1 C C START OF LOOP C DO 200 NPNT=1,MAXPNTS X1 = P(1,NPNT) X2 = P(2,NPNT) X3 = P(3,NPNT) C C CALCULATE THE ONE-ELECTRON INTEGRALS C IF (IPO(5).EQ.1) THEN WRITE(IOUT,3010) 3010 FORMAT(1X,'TIME FOR INTEGRALS') C*** C ISTAT = LIB$INIT_TIMER(HANDLE) C*** END IF C CALL INTGRL (HPERT,X1,X2,X3,IPTCHG,I6TO5) C C*** C IF (IPO(5).EQ.1) ISTAT = LIB$SHOW_TIMER(HANDLE) C*** C IF (IPO(4).EQ.1) CALL LINOUT (HPERT,NBASIS,0,0) C IF (IPO(5).EQ.1) THEN WRITE(IOUT,3000) 3000 FORMAT(1X,'TIME FOR TRANSFORM') C*** C ISTAT = LIB$INIT_TIMER(HANDLE) C*** END IF C C FORM THE HPERT MATRIX ELEMENTS C E = ZERO ICOEFI = -NBASIS C C SUM OVER OCCUPIED MOS C DO 220 II=1,NOCC ICOEFI = ICOEFI + NBASIS C C CALCULATE ELECTROSTATIC POTENTIAL C DO 221 IP=1,NBASIS CPI = COEF_ALPHA(ICOEFI+IP) IPDEX = INDEX(IP) C DO 222 IQ=1,IP E = E + CPI * COEF_ALPHA(ICOEFI+IQ) * HPERT(IPDEX+IQ) 222 CONTINUE DO 223 IQ=IP+1,NBASIS E = E + CPI * COEF_ALPHA(ICOEFI+IQ) * HPERT(IP+INDEX(IQ)) 223 CONTINUE C 221 CONTINUE 220 CONTINUE C C*** C IF (IPO(5).EQ.1) ISTAT = LIB$SHOW_TIMER(HANDLE) C*** C C CALCULATE NUCLEAR PART OF ELECTROSTATIC POTENTIAL C VNUC = ZERO DO 300 IATOM=1,NATOMS DEL1 = C(1,IATOM) - X1 DEL2 = C(2,IATOM) - X2 DEL3 = C(3,IATOM) - X3 RA = DSQRT(DEL1*DEL1 + DEL2*DEL2 + DEL3*DEL3) IF (RA.EQ.ZERO) THEN VNUC=VNUCMAX GOTO 310 END IF VNUC = VNUC + IAN(IATOM) / RA 300 CONTINUE 310 CONTINUE C ELECP(NPNT) = (E * TWO + VNUC * IPTCHG) IF (IPO(5) .EQ. 1) WRITE(IOUT,*) 'E(',NPNT,') = ',E 200 CONTINUE RETURN END SUBROUTINE FIT C C ROUTINE TO USE METHOD OF LAGRANGE MULTIPLIERS TO OBTAIN BEST C LEAST SQUARE FIT WITH CONSTRAINTS C C M.M. FRANCL C APRIL 1985 C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NPOINTS = 50000) INTEGER*4 WHICH1 CHARACTER*40 CHKFIL C COMMON /IO/ IN,IOUT COMMON /IPO/ IPO(5) COMMON /ELP/ E(NPOINTS) COMMON /POINTS/ P(3,NPOINTS),MAXPNTS COMMON /OUT/ NTITLE(20,3),CHKFIL,X(400),I6TO5,RMS,PERCENT, 1 NLIN,NEND(3) C+++ COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS, $ IAN(401),ATMCHG(400),C(3,400) C+++ C COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101), C $ ATMCHG(100),C(3,100) C DIMENSION A(400,400),Y(400),IS(2,400),IAD1(400),IAD2(400) DIMENSION D(400),WHICH1(3) C C DEBYE = CONVERSION FROM DEBYES TO AU C DATA ONE/1.0/, ZERO/0.0/, DEBYE/0.393427328/, MAXDIM/400/ DATA AU2CAL/627.51/, HALF/0.5/, HUNDRED/100.0/,NCONSTR/1/ C C SET UP MATRIX OF LINEAR COEFFICIENTS, A C C BEGIN LOOP OVER ROWS C DO 100 K=1,NATOMS C C BEGIN LOOP OVER COLUMNS C DO 200 MU=1,NATOMS C SUM = ZERO DO 400 I=1,MAXPNTS RIK = (P(1,I)-C(1,K))**2 + (P(2,I)-C(2,K))**2 + (P(3,I)-C(3,K))**2 RIK = DSQRT(RIK) RIMU = (P(1,I)-C(1,MU))**2 + (P(2,I)-C(2,MU))**2 + $ (P(3,I)-C(3,MU))**2 RIMU = DSQRT(RIMU) SUM = SUM + ONE / (RIK * RIMU) 400 CONTINUE C A(K,MU) = SUM 200 CONTINUE C C FILL OUT COLUMNS CORRESPONDING TO LAGRANGE MULTIPLIERS C A(K,NATOMS+1) = HALF C C 100 CONTINUE C C FILL OUT THE ROWS CORRESPONDING TO CONSTRAINTS C DO 500 MU=1,NATOMS A(NATOMS+1,MU) = ONE C 500 CONTINUE C C FILL OUT THE BLOCK WHICH CONNECTS LAGRANGE MULTIPLIERS TO C CONSTRAINTS C DO 600 K=NATOMS+1,NATOMS+NCONSTR DO 600 MU=NATOMS+1,NATOMS+NCONSTR A(K,MU) = ZERO 600 CONTINUE C C****DEBUG***** C IF (IPO(3) .EQ. 1) THEN WRITE(IOUT,*) 'A MATRIX' DO 699 K=1,NATOMS+NCONSTR WRITE(IOUT,1699) (A(K,MU),MU=1,NATOMS+NCONSTR) 1699 FORMAT(1X,10F10.4) 699 CONTINUE END IF C*************** C C CONSTRUCT COLUMN VECTOR, Y C DO 700 K=1,NATOMS SUM = ZERO DO 710 I=1,MAXPNTS RIK = (P(1,I)-C(1,K))**2 + (P(2,I)-C(2,K))**2 + $ (P(3,I)-C(3,K))**2 RIK = DSQRT(RIK) SUM = SUM + E(I) / RIK 710 CONTINUE Y(K) = SUM IF (IPO(3) .EQ. 1) WRITE(IOUT,*) K,Y(K) 700 CONTINUE C C CONSTRUCT THE PORTION OF Y CORRESPONDING TO LAGRANGE MULTIPLIERS C C Y(NATOMS+1) = DFLOAT(ICHARG) C C IF (IPO(3) .EQ. 1) $ WRITE(IOUT,*) 'COL VECTR Y', (Y(KK),KK=1,NATOMS+NCONSTR) C C SOLVE MATRIX EQUATION AX = Y; C WHERE X = (Q1,Q2, ... QN,L1,L2, ... ,LN) C C X = A(INV)Y C C INVERT A C CALL INV(A,NATOMS+NCONSTR,IS,IAD1,IAD2,D,MAXDIM) C C****DEBUG***** C IF (IPO(3) .EQ. 1) THEN WRITE(IOUT,*) 'A INVERSE' DO 799 K=1,NATOMS+NCONSTR WRITE(IOUT,1699) (A(K,MU),MU=1,NATOMS+NCONSTR) 799 CONTINUE END IF C************** C C PERFORM MATRIX MULTIPLICATION A(INV)Y C CALL MULTAY(A,Y,X,NATOMS+NCONSTR,MAXDIM) C IF (IPO(3) .EQ. 1) THEN WRITE(IOUT,*) 'CHARGES: ' DO 899 I=1,NATOMS WRITE(IOUT,*) IAN(I),X(I) 899 CONTINUE END IF C C COMPUTE RMS DEVIATION AND MEAN ABSOLUTE % DEVIATION C RMS = ZERO PERCENT = ZERO DO 800 I=1,MAXPNTS EQ = ZERO DO 810 J=1,NATOMS DIST = (P(1,I)-C(1,J))**2 + (P(2,I)-C(2,J))**2 + $ (P(3,I)-C(3,J))**2 DIST = DSQRT(DIST) EQ = EQ + X(J) / DIST 810 CONTINUE RMS = RMS + (E(I) - EQ)**2 PERCENT = PERCENT + DABS((E(I) - EQ) / E(I) * HUNDRED) IF (IPO(3) .EQ. 1) WRITE(IOUT,*) 'ACTUAL,CALC ',E(I),EQ 800 CONTINUE IF (IPO(3) .EQ. 1) WRITE(IOUT,*) 'SUM OF SQUARES ',RMS RMS = DSQRT(RMS) * AU2CAL / MAXPNTS PERCENT = PERCENT / MAXPNTS IF (IPO(3) .EQ. 1) WRITE(IOUT,*) 'RMS, %',RMS,PERCENT RETURN END SUBROUTINE FMGEN(F,T,M) C IMPLICIT REAL*8 (A-H,O-Z) COMMON/IO/IN,IOUT C DIMENSION F(M) DIMENSION GA(35) C EQUIVALENCE (APPROX,OLDSUM) C DATA ZERO/0.0E0/, HALF/0.5E0/, ONE/1.0E0/, TWO/2.0E0/, TEN/10.0E0/ $ ,PI/3.14159265358979E0/, F42/42.0E0/, F80/80.0E0/ C 2001 FORMAT(42H1FAILURE IN FMGEN FOR SMALL T: IX.GT.50, / $ 6H IX = ,I3,7H, T = ,E20.14) 2002 FORMAT(37H1FAILURE IN FMGEN FOR INTERMEDIATE T,/ $ 6H T = ,E20.14) C TEXP=ZERO IF(T-F80)2,3,3 2 TEXP=EXP(-T) 3 CONTINUE IF(T-TEN)10,70,70 C*********************************************************************** C 0 .LT. T .LT. 10 C*********************************************************************** 10 TERM=HALF*GA(M)*TEXP TX=ONE IX=M+1 SUM=TX/GA(IX) OLDSUM=SUM 20 IX=IX+1 TX=TX*T IF(IX - 35) 40,40,30 30 WRITE(IOUT,2001)IX,T STOP 'FMGEN' 40 SUM=SUM+TX/GA(IX) IF(TOL-ABS(OLDSUM/SUM-ONE))50,60,60 50 OLDSUM=SUM GO TO 20 60 F(M)=SUM*TERM GO TO 160 C 70 IF(T-F42)80,150,150 C*********************************************************************** C 10 .LE. T .LT. 42 C*********************************************************************** 80 A=FLOAT(M-1) B=A+HALF A=A-HALF TX=ONE/T MM1=M-1 APPROX=RPITWO*SQRT(TX)*(TX**MM1) IF(MM1)90,110,90 90 DO 100 IX=1,MM1 B=B-ONE 100 APPROX=APPROX*B 110 FIMULT=HALF*TEXP*TX SUM=ZERO IF(FIMULT)120,140,120 120 FIPROP=FIMULT/APPROX TERM=ONE SUM =ONE NOTRMS=INT(T)+MM1 DO 130 IX=2,NOTRMS TERM=TERM*A*TX SUM=SUM+TERM IF(ABS(TERM*FIPROP/SUM)-TOL)140,140,130 130 A=A-ONE WRITE(IOUT,2002)T STOP 'FMGEN' 140 F(M)=APPROX-FIMULT*SUM GO TO 160 C*********************************************************************** C T .GE. 42 C*********************************************************************** 150 TX=FLOAT(M)-HALF F(M)=HALF*GA(M)/(T**TX) C*********************************************************************** C RECUR DOWNWARDS TO F(1) C*********************************************************************** 160 TX=T+T SUM=FLOAT(M+M-3) MM1=M-1 IF(MM1)170,190,170 170 DO 180 IX=1,MM1 F(M-IX)=(TX*F(M-IX+1)+TEXP)/SUM 180 SUM=SUM-TWO 190 RETURN C ENTRY FMSET C GA(1)=SQRT(PI) TOL=HALF DO 200 I=2,35 GA(I)=GA(I-1)*TOL 200 TOL=TOL+ONE TOL = 5.0E-09 RPITWO=HALF*GA(1) RETURN END SUBROUTINE INTGRL (H,X1,X2,X3,ICHARG,I6TO5) C C ROUTINE TO CALCULATE THE ELECTRON-CHARGE MATRIX ELEMENTS FOR THE C POLARIZATION POTENTIAL. CODE REVISED FROM THE ONE ELECTRON PACKAGE C AS IT EXISTED AUGUST, 1983. C C C REVISED BY M.M. FRANCL JANUARY 1984 FOR PRINCETON CHEMISTRY C DEPARTMENT VAX 11/780 C C REVISED TO BE COMPATIBLE WITH COMMON /B/ FROM GAUSSIAN 82 C MAY 1984 M.M. FRANCL C C REVISED TO USE ** BASIS SETS AND THOSE HAVING P ONLY SHELLS C JANUARY 1986 M.M. FRANCL C C REVISED FOR GAUSSIAN 86 CHECKPOINT FILES FOR YALE UNIVERSITY C FEBRUARY 1988 CURT BRENEMAN C C IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4 SHELLA,SHELLN,SHELLT,SHELLC,AOS,AON,SHLADF C C+++ COMMON /MOL/ NATOMS,JCHARG,MULTIP,NAE,NBE,NEL,NBASIS, $ IAN(401),ATMCHG(400),C(3,400) C C=== Gaussian88 Modification. New Common /b/ size. Common/B/EXX(6000),C1(6000),C2(6000),C3(2000),CF(2000), $SHLADF(4000),X(2000),Y(2000), $Z(2000),JAN(2000),ShellA(2000),ShellN(2000),ShellT(2000), $ShellC(2000),AOS(2000),AON(2000),NShell,MaxTyp C C=== Old G86 common /b/ c COMMON/B/EXX(1200),C1(1200),C2(1200),C3(400),CF(400),SHLADF(800), c $ X(400),Y(400),Z(400),JAN(400),SHELLA(400),SHELLN(400), c $ SHELLT(400),SHELLC(400),AOS(400),AON(400),NSHELL,MAXTYP c C+++ C COMMON /B/ EXX(240),C1(240),C2(240),C3(80),CF(80),SHLADF(160), C $ X(80),Y(80),Z(80), C $ JAN(80),SHELLA(80),SHELLN(80),SHELLT(80),SHELLC(80) C $ ,AOS(80),AON(80),NSHELL,MAXTYP C COMMON /MOL/ NATOMS,JCHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101), C $ ATMCHG(100),C(3,100) COMMON /IPO/ IPO(10) COMMON/IO/ IN,IOUT C DIMENSION H(1) DIMENSION RENORM(10) DIMENSION OF(9),OX(9),TX(13),ABX(5),ABY(5),ABZ(5),ABSQ(5), *A(5),B(5),F(5),APB(5),CPX(5),CPY(5),CPZ(5),FM(5) DIMENSION EPN(100) C COMMON/H100/ $EP00,EP10,EP20,EP30,EP40,EP50,EP60,EP70,EP80,EP90, SPDSTV $EP01,EP11,EP21,EP31,EP41,EP51,EP61,EP71,EP81,EP91, SPDSTV $EP02,EP12,EP22,EP32,EP42,EP52,EP62,EP72,EP82,EP92, SPDSTV $EP03,EP13,EP23,EP33,EP43,EP53,EP63,EP73,EP83,EP93, SPDSTV $EP04,EP14,EP24,EP34,EP44,EP54,EP64,EP74,EP84,EP94, SPDSTV $EP05,EP15,EP25,EP35,EP45,EP55,EP65,EP75,EP85,EP95, SPDSTV $EP06,EP16,EP26,EP36,EP46,EP56,EP66,EP76,EP86,EP96, SPDSTV $EP07,EP17,EP27,EP37,EP47,EP57,EP67,EP77,EP87,EP97, SPDSTV $EP08,EP18,EP28,EP38,EP48,EP58,EP68,EP78,EP88,EP98, SPDSTV $EP09,EP19,EP29,EP39,EP49,EP59,EP69,EP79,EP89,EP99 SPDSTV C DIMENSION EEP(100) SPDSTV DIMENSION MAX(6) C SPDSTV C LOCAL VARIABLES. SPDSTV C DIMENSION AG(6),CSA(6),CPA(6),CDA(6), SPDSTV $ BG(6),CSB(6),CPB(6),CDB(6), SPDSTV $ DPP(9) SPDSTV EQUIVALENCE(OF0,OF(1)),(OF1,OF(2)),(OF2,OF(3)), SPDSTV $ (OF3,OF(4)),(OF4,OF(5)),(OF5,OF(6)), SPDSTV $ (OF6,OF(7)),(OF7,OF(8)),(OF8,OF(9)) SPDSTV EQUIVALENCE(OX0,OX(1)),(OX1,OX(2)),(OX2,OX(3)), SPDSTV $ (OX3,OX(4)),(OX4,OX(5)),(OX5,OX(6)), SPDSTV $ (OX6,OX(7)),(OX7,OX(8)),(OX8,OX(9)) SPDSTV EQUIVALENCE(A1,A(2)),(A2,A(3)),(A3,A(4)),(A4,A(5)) SPDSTV EQUIVALENCE(B1,B(2)),(B2,B(3)),(B3,B(4)),(B4,B(5)) SPDSTV EQUIVALENCE(T01,T0),(T02,T1),(T03,T2), SPDSTV $ (T04,T3),(T05,T4),(T06,T5), SPDSTV $ (T07,T6),(T08,T7),(T09,T8) SPDSTV EQUIVALENCE(T10,TX(10)),(T11,TX(11)),(T12,TX(12)),(T13,TX(13)) SPDSTV EQUIVALENCE(T0,TX(1)),(T1,TX(2)),(T2,TX(3)), SPDSTV $ (T3,TX(4)),(T4,TX(5)),(T5,TX(6)), SPDSTV $ (T6,TX(7)),(T7,TX(8)),(T8,TX(9)) SPDSTV EQUIVALENCE(C001,T01),(C050,T02),(C054,T09), SPDSTV $ (C067,T13),(C068,T08),(C074,T03) SPDSTV EQUIVALENCE(ABX1,ABX(2)),(ABX2,ABX(3)), SPDSTV $ (ABX3,ABX(4)),(ABX4,ABX(5)) SPDSTV EQUIVALENCE(AB004,ABX1),(AB006,ABX2),(AB023,ABX3),(AB029,ABX4) SPDSTV EQUIVALENCE(ABY1,ABY(2)),(ABY2,ABY(3)), SPDSTV $ (ABY3,ABY(4)),(ABY4,ABY(5)) SPDSTV EQUIVALENCE(AB007,ABY1),(AB010,ABY2),(AB032,ABY3),(AB035,ABY4) SPDSTV EQUIVALENCE(ABZ1,ABZ(2)),(ABZ2,ABZ(3)), SPDSTV $ (ABZ3,ABZ(4)),(ABZ4,ABZ(5)) SPDSTV EQUIVALENCE(AB002,ABZ1),(AB003,ABZ2),(AB011,ABZ3),(AB017,ABZ4) SPDSTV EQUIVALENCE(ABSQ1,ABSQ(2)),(ABSQ2,ABSQ(3)), SPDSTV $ (ABSQ3,ABSQ(4)),(ABSQ4,ABSQ(5)) SPDSTV EQUIVALENCE(APB1,APB(2)),(APB2,APB(3)), SPDSTV $ (APB3,APB(4)),(APB4,APB(5)) SPDSTV EQUIVALENCE(CPX1,CPX(2)),(CPX2,CPX(3)), SPDSTV $ (CPX3,CPX(4)),(CPX4,CPX(5)) SPDSTV EQUIVALENCE(CPY1,CPY(2)),(CPY2,CPY(3)), SPDSTV $ (CPY3,CPY(4)),(CPY4,CPY(5)) SPDSTV EQUIVALENCE(CPZ1,CPZ(2)),(CPZ2,CPZ(3)), SPDSTV $ (CPZ3,CPZ(4)),(CPZ4,CPZ(5)) SPDSTV EQUIVALENCE(F1,F(2)),(F2,F(3)),(F3,F(4)),(F4,F(5)) SPDSTV EQUIVALENCE(FM0,FM(1)),(FM1,FM(2)),(FM2,FM(3)),(FM3,FM(4)), SPDSTV $ (FM4,FM(5)) SPDSTV EQUIVALENCE (D001,FM0) SPDSTV EQUIVALENCE(EP00,EEP(1)) SPDSTV C DATA MAX/1,4,9,1,4,10/ SPDSTV DATA TX/1.0E0,0.5E0,0.25E0,0.125E0,0.375E0,0.625E-01,0.1875E0, $ 0.75E0,1.5E0,2.25E0,1.125E0,0.0E0,3.0E0/ DATA ZERO/0.0/,HALF/0.5/,ONE/1.0/,ONEPT5/1.5/,TWO/2.0/,THREE/3.0/, *ROOT3/1.732050808/,PI/3.14159265358979/ DATA ANTOAU /1.889726878D0/ C 2010 FORMAT(/1X,'ELECTRON-CHARGE MATRIX ELEMENTS'/) C C CALL ROUTINE TO MODIFY COMMON /B/ IF P ONLY SHELLS ARE PRESENT C CALL STAR (NBASIS,SHELLT,SHELLC,AOS,NSHELL,NOSTAR) C C*********************************************************************** SPDSTV C INITIALIZE THIS SEGMENT. SPDSTV C*********************************************************************** SPDSTV C SPDSTV C ****************************************************************** SPDSTV C COMPUTE SIZE OF S T AND V ARRAYS C ****************************************************************** SPDSTV NTT=(NBASIS*(NBASIS+1))/2 I5OR6=3 CC IF(IGO(4) .NE. 0) I5OR6 = 0 C ****************************************************************** SPDSTV C INITIALIZE RENORM USED TO NORMALIZE D FUNCTIONS C ****************************************************************** SPDSTV DO 100 I=1,10 SPDSTV 100 RENORM(I)=ONE SPDSTV RENORM(5)=ROOT3 SPDSTV RENORM(8)=ROOT3 SPDSTV RENORM(9)=ROOT3 SPDSTV C ****************************************************************** SPDSTV C CLEAR H ARRAY C ****************************************************************** SPDSTV DO 50 I=1,NTT SPDSTV 50 H(I)=ZERO C ****************************************************************** SPDSTV C * INITIALIZE THE VARIABLES USED BY ROUTINE FMGEN. * SPDSTV C ****************************************************************** SPDSTV CALL FMSET DO 95 I=1,5 95 FM(I)=ZERO ABX(1)=ONE SPDSTV ABY(1)=ONE SPDSTV ABZ(1)=ONE SPDSTV A(1)=ONE SPDSTV B(1)=ONE SPDSTV F(1)=ONE SPDSTV CPX(1)=ONE SPDSTV CPY(1)=ONE SPDSTV CPZ(1)=ONE SPDSTV APB(1)=ONE SPDSTV ABSQ(1)=ONE SPDSTV C*********************************************************************** SPDSTV C LOOP OVER SHELLS ISHELL AND JSHELL. SPDSTV C*********************************************************************** SPDSTV DO 1000 ISHELL=1,NSHELL SPDSTV DO 1000 JSHELL=1,ISHELL SPDSTV SYMFAC = ONE C ****************************************************************** SPDSTV C ZERO LOCATIONS SPDSTV C ****************************************************************** SPDSTV 80 CONTINUE DO 9447 JI=1,100 SPDSTV EPN(JI)=ZERO SPDSTV 9447 CONTINUE SPDSTV IF(SHELLT(ISHELL)-SHELLT(JSHELL))120,120,110 SPDSTV 110 INEW=JSHELL SPDSTV JNEW=ISHELL SPDSTV LA=SHELLT(JSHELL) SPDSTV LB=SHELLT(ISHELL) SPDSTV GO TO 200 SPDSTV 120 INEW=ISHELL SPDSTV JNEW=JSHELL SPDSTV LA=SHELLT(ISHELL) SPDSTV LB=SHELLT(JSHELL) SPDSTV 200 CONTINUE SPDSTV LAP1=LA+1 SPDSTV LBP1=LB+1 SPDSTV LAMAX=MAX(LAP1+I5OR6) SPDSTV LBMAX=MAX(LBP1+I5OR6) SPDSTV ITYPE=3*LB+LA SPDSTV M=LA+LB+1 SPDSTV NGA=SHELLN(INEW) SPDSTV NGB=SHELLN(JNEW) SPDSTV AX=X(INEW) SPDSTV BX=X(JNEW) SPDSTV AY=Y(INEW) SPDSTV BY=Y(JNEW) SPDSTV AZ=Z(INEW) SPDSTV BZ=Z(JNEW) SPDSTV ISHA=SHELLA(INEW) SPDSTV ISHB=SHELLA(JNEW) SPDSTV ISHAD = SHLADF(INEW) ISHBD = SHLADF(JNEW) IAOS=AOS(INEW) SPDSTV JAOS=AOS(JNEW) SPDSTV C ****************************************************************** SPDSTV C OBTAIN INFORMATION ABOUT SHELLS INEW AND JNEW SPDSTV C ****************************************************************** SPDSTV DO 101 I=1,NGA SPDSTV N=ISHA+I-1 SPDSTV ND = ISHAD + I -1 IF (MAXTYP .LE. 1) ND=1 AG(I)=EXX(N) SPDSTV CSA(I)=C1(N) SPDSTV CPA(I)=C2(N) SPDSTV 101 CDA(I)=C3(ND) SPDSTV DO 102 I=1,NGB SPDSTV N=ISHB+I-1 SPDSTV ND = ISHBD + I -1 BG(I)=EXX(N) SPDSTV CSB(I)=C1(N) SPDSTV CPB(I)=C2(N) SPDSTV 102 CDB(I)=C3(ND) SPDSTV ABX(2)=BX-AX SPDSTV ABY(2)=BY-AY SPDSTV ABZ(2)=BZ-AZ SPDSTV RABSQ=ABX(2)*ABX(2)+ABY(2)*ABY(2)+ABZ(2)*ABZ(2) SPDSTV ABSQ(2)=RABSQ SPDSTV DO 103 I=3,5 SPDSTV ABX(I)=ABX(I-1)*ABX(2) SPDSTV ABY(I)=ABY(I-1)*ABY(2) SPDSTV ABZ(I)=ABZ(I-1)*ABZ(2) SPDSTV 103 ABSQ(I)=ABSQ(I-1)*ABSQ(2) SPDSTV AB001=ONE SPDSTV AB005=ABX1*ABZ1 SPDSTV AB008=ABY1*ABZ1 SPDSTV AB009=ABX1*ABY1 SPDSTV AB012=ABX1*ABZ2 SPDSTV AB013=ABX2*ABZ1 SPDSTV AB014=ABY1*ABZ2 SPDSTV AB015=ABX1*ABY1*ABZ1 SPDSTV AB016=ABY2*ABZ1 SPDSTV AB018=ABX1*ABZ3 SPDSTV AB019=ABX2*ABZ2 SPDSTV AB020=ABY1*ABZ3 SPDSTV AB021=ABX1*ABY1*ABZ2 SPDSTV AB022=ABY2*ABZ2 SPDSTV AB024=ABX2*ABY1 SPDSTV AB025=ABX1*ABY2 SPDSTV AB026=ABX3*ABZ1 SPDSTV AB027=ABX2*ABY1*ABZ1 SPDSTV AB028=ABX1*ABY2*ABZ1 SPDSTV AB030=ABX3*ABY1 SPDSTV AB031=ABX2*ABY2 SPDSTV AB033=ABY3*ABZ1 SPDSTV AB034=ABX1*ABY3 SPDSTV C*********************************************************************** SPDSTV C LOOP OVER GAUSSIANS (CONTRACTION LOOP). SPDSTV C*********************************************************************** SPDSTV DO 105 IGAUSS=1,NGA SPDSTV AA=AG(IGAUSS) SPDSTV DO 105 JGAUSS=1,NGB SPDSTV BB=BG(JGAUSS) SPDSTV AAPBB=AA+BB SPDSTV APBB=ONE/AAPBB SPDSTV F2=TWO*AA*BB*APBB SPDSTV PX=(AA*AX+BB*BX)*APBB SPDSTV PY=(AA*AY+BB*BY)*APBB SPDSTV PZ=(AA*AZ+BB*BZ)*APBB SPDSTV A(2)=ONE/AA SPDSTV B(2)=ONE/BB SPDSTV F(2)=F2 SPDSTV APB(2)=APBB SPDSTV YX=PI*APBB SPDSTV EXX1=HALF*F2*RABSQ SPDSTV IF(EXX1-80.0E0)4172,4173,4173 4173 EXX1=ZERO GO TO 4714 SPDSTV 4172 EXX1=EXP(-EXX1) 4714 CONTINUE SPDSTV OV=(YX**ONEPT5)*EXX1 SPDSTV OVEK=THREE*AA*BB*APBB SPDSTV EK=F2*AA*BB*APBB*OV SPDSTV EP=TWO*YX*EXX1 SPDSTV DO 119 I=3,5 SPDSTV A(I)=A(I-1)*A(2) SPDSTV B(I)=B(I-1)*B(2) SPDSTV APB(I)=APB(I-1)*APB(2) SPDSTV 119 F(I)=F(I-1)*F(2) SPDSTV DPP(1)=CSA(IGAUSS)*CSB(JGAUSS) SPDSTV DPP(2)=CPA(IGAUSS)*CSB(JGAUSS) SPDSTV DPP(3)=CDA(IGAUSS)*CSB(JGAUSS) SPDSTV DPP(4)=CSA(IGAUSS)*CPB(JGAUSS) SPDSTV DPP(5)=CPA(IGAUSS)*CPB(JGAUSS) SPDSTV DPP(6)=CDA(IGAUSS)*CPB(JGAUSS) SPDSTV DPP(7)=CSA(IGAUSS)*CDB(JGAUSS) SPDSTV DPP(8)=CPA(IGAUSS)*CDB(JGAUSS) SPDSTV DPP(9)=CDA(IGAUSS)*CDB(JGAUSS) SPDSTV DO 2132 I=1,9 SPDSTV OF(I)=DPP(I)*OV SPDSTV 2132 OX(I)=DPP(I)*EK SPDSTV DO 2139 I=1,100 2139 EEP(I)=ZERO C002=T02*A1*F1 SPDSTV C006=T02*B1*F1 SPDSTV C007=T03*A1*B1*F2 SPDSTV C008=T03*A1*B1*F1 SPDSTV C027=T01*A1 SPDSTV C031=T01*A1*B1*F1 SPDSTV C032=T02*A1*B1 SPDSTV C051=T02*A1*B1*F2 SPDSTV C012=T02*B1 SPDSTV C013=T03*B2*F2 SPDSTV C014=T03*B2*F1 SPDSTV C036=T01*B2*F1 SPDSTV C037=T02*B2 SPDSTV C056=T01*B1*F1 SPDSTV C030=T01*B1 SPDSTV C018=T04*A1*B2*F2 SPDSTV IF(ITYPE-7)3060,3040,3041 SPDSTV 3041 CONTINUE SPDSTV C003=T02*A1 SPDSTV C004=T03*A2*F2 SPDSTV C005=T03*A2*F1 SPDSTV C009=T04*A2*B1*F3 SPDSTV C010=T05*A2*B1*F2 SPDSTV C011=T04*A2*B1*F2 SPDSTV C017=T03*A1*B1 SPDSTV C019=T04*A1*B2*F1 SPDSTV C020=T04*A2*B1*F1 SPDSTV C021=T06*A2*B2*F4 SPDSTV C022=T05*A2*B2*F3 SPDSTV C023=T07*A2*B2*F2 SPDSTV C024=T07*A2*B2*F3 SPDSTV C025=T06*A2*B2*F3 SPDSTV C026=T06*A2*B2*F2 SPDSTV C028=T01*A2*F1 SPDSTV C029=T02*A2 SPDSTV C033=T08*A2*B1*F2 SPDSTV C034=T09*A2*B1*F1 SPDSTV C035=T02*A2*B1*F1 SPDSTV C040=T02*A1*B2*F1 SPDSTV C041=T03*A1*B2 SPDSTV C042=T03*A2*B1 SPDSTV C043=T02*A2*B2*F3 SPDSTV C044=T10*A2*B2*F2 SPDSTV C045=T08*A2*B2*F1 SPDSTV C046=T11*A2*B2*F2 SPDSTV C047=T05*A2*B2*F2 SPDSTV C048=T03*A2*B2*F1 SPDSTV C049=T01*A1*F1 SPDSTV C057=T12*A1*B1*F1 SPDSTV C058=T03*A1 SPDSTV C059=T03*B1 SPDSTV C060=T03*A1*B2*F3 SPDSTV C061=T04*B2*F2 SPDSTV C062=T04*B2*F1 SPDSTV C063=T03*A2*B1*F3 SPDSTV C064=T01*A1*B1*F2 SPDSTV C065=T09*B1*F1 SPDSTV C066=T09*A1*F1 SPDSTV C069=T04*A2*F2 SPDSTV C070=T04*A2*F1 SPDSTV C071=T03*A2*B1*F2 SPDSTV C072=T08*A1*F1 SPDSTV C073=T03*A1*B2*F2 SPDSTV C075=T08*B1*F1 SPDSTV C076=T04*A1*B1*F2 SPDSTV 3040 CONTINUE SPDSTV C015=T04*A1*B2*F3 SPDSTV C016=T05*A1*B2*F2 SPDSTV C038=T08*A1*B2*F2 SPDSTV C039=T09*A1*B2*F1 SPDSTV C040=T02*A1*B2*F1 SPDSTV C052=T02*A1*B1*F1 SPDSTV C053=T03*B1*F1 SPDSTV C055=T03*A1*F1 SPDSTV 3060 CONTINUE CX=X1 CY=X2 CZ=X3 CPX(2)=PX-CX SPDSTV CPY(2)=PY-CY SPDSTV CPZ(2)=PZ-CZ SPDSTV CP2=CPX(2)*CPX(2)+CPY(2)*CPY(2)+CPZ(2)*CPZ(2) SPDSTV CALL FMGEN(FM,AAPBB*CP2,M) SPDSTV DO 108 I=3,5 SPDSTV CPX(I)=CPX(I-1)*CPX(2) SPDSTV CPY(I)=CPY(I-1)*CPY(2) SPDSTV 108 CPZ(I)=CPZ(I-1)*CPZ(2) SPDSTV EPAN=EP*FLOAT(-ICHARG) DO 2136 I=1,9 SPDSTV 2136 OF(I)=DPP(I)*EPAN SPDSTV D002=CPZ1*FM1 SPDSTV D003=CPZ2*FM2 SPDSTV D004=APB1*FM1 SPDSTV D005=CPX1*FM1 SPDSTV D006=CPX1*CPZ1*FM2 SPDSTV D007=CPX2*FM2 SPDSTV D008=CPY1*FM1 SPDSTV D009=CPY1*CPZ1*FM2 SPDSTV D010=CPX1*CPY1*FM2 SPDSTV D011=CPY2*FM2 SPDSTV D012=CPZ3*FM3 SPDSTV D013=APB1*CPZ1*FM2 SPDSTV D014=CPX1*CPZ2*FM3 SPDSTV D015=APB1*CPX1*FM2 SPDSTV D016=CPX2*CPZ1*FM3 SPDSTV D017=CPY1*CPZ2*FM3 SPDSTV D018=APB1*CPY1*FM2 SPDSTV D019=CPX1*CPY1*CPZ1*FM3 SPDSTV D020=CPY2*CPZ1*FM3 SPDSTV D034=CPX3*FM3 SPDSTV D035=CPX2*CPY1*FM3 SPDSTV D036=CPX1*CPY2*FM3 SPDSTV D043=CPY3*FM3 SPDSTV C ****************************************************************** SPDSTV C * SS * SPDSTV C ****************************************************************** SPDSTV EP00=OF0*(+C001*AB001*D001) SPDSTV IF(ITYPE)3230,3262,3230 SPDSTV C ****************************************************************** SPDSTV C * SP * SPDSTV C ****************************************************************** SPDSTV 3230 CONTINUE SPDSTV EP01=OF3*(-C006*AB002*D001-C001*AB001*D002) SPDSTV EP03=OF3*(-C006*AB004*D001-C001*AB001*D005) SPDSTV EP06=OF3*(-C006*AB007*D001-C001*AB001*D008) SPDSTV IF(ITYPE-7)3240,3242,3241 SPDSTV 3240 IF(ITYPE-4)3262,3261,3260 SPDSTV C ****************************************************************** SPDSTV C * DD * SPDSTV C ****************************************************************** SPDSTV 3241 CONTINUE SPDSTV D021=CPZ4*FM4 SPDSTV D022=APB1*CPZ2*FM3 SPDSTV D023=APB2*FM2 SPDSTV D024=CPX1*CPZ3*FM4 SPDSTV D025=APB1*CPX1*CPZ1*FM3 SPDSTV D026=CPX2*CPZ2*FM4 SPDSTV D027=APB1*CPX2*FM3 SPDSTV D028=CPY1*CPZ3*FM4 SPDSTV D029=APB1*CPY1*CPZ1*FM3 SPDSTV D030=CPX1*CPY1*CPZ2*FM4 SPDSTV D031=APB1*CPX1*CPY1*FM3 SPDSTV D032=CPY2*CPZ2*FM4 SPDSTV D033=APB1*CPY2*FM3 SPDSTV D037=CPX3*CPZ1*FM4 SPDSTV D038=CPX2*CPY1*CPZ1*FM4 SPDSTV D039=CPX1*CPY2*CPZ1*FM4 SPDSTV D040=CPX4*FM4 SPDSTV D041=CPX3*CPY1*FM4 SPDSTV D042=CPX2*CPY2*FM4 SPDSTV D044=CPY3*CPZ1*FM4 SPDSTV D045=CPX1*CPY3*FM4 SPDSTV D046=CPY4*FM4 SPDSTV EP20=OF2*(+C003*AB001*D001+C004*AB003*D001-C005*AB001*D001-C049*AB SPDSTV $002*D002+C001*AB001*D003-C050*AB001*D004) SPDSTV EP40=OF2*(+C004*AB005*D001-C002*AB004*D002-C002*AB002*D005+C001*AB SPDSTV $001*D006) SPDSTV EP50=OF2*(+C003*AB001*D001+C004*AB006*D001-C005*AB001*D001-C049*AB SPDSTV $004*D005+C001*AB001*D007-C050*AB001*D004) SPDSTV EP70=OF2*(+C004*AB008*D001-C002*AB007*D002-C002*AB002*D008+C001*AB SPDSTV $001*D009) SPDSTV EP80=OF2*(+C004*AB009*D001-C002*AB004*D008-C002*AB007*D005+C001*AB SPDSTV $001*D010) SPDSTV EP90=OF2*(+C003*AB001*D001+C004*AB010*D001-C005*AB001*D001-C049*AB SPDSTV $007*D008+C001*AB001*D011-C050*AB001*D004) SPDSTV EP21=OF5*(-C008*AB002*D001-C003*AB001*D002-C009*AB011*D001+C010*AB SPDSTV $002*D001+C051*AB003*D002-C052*AB001*D002-C006*AB002*D003+C053*AB00 SPDSTV $2*D004-C004*AB003*D002+C005*AB001*D002+C049*AB002*D003-C002*AB002* SPDSTV $D004-C001*AB001*D012+C054*AB001*D013) SPDSTV EP41=OF5*(-C009*AB012*D001+C011*AB004*D001+C007*AB005*D002+C007*AB SPDSTV $003*D005-C008*AB001*D005-C006*AB002*D006-C004*AB005*D002+C002*AB00 SPDSTV $4*D003-C055*AB004*D004+C002*AB002*D006-C001*AB001*D014+C050*AB001* SPDSTV $D015) SPDSTV EP51=OF5*(-C008*AB002*D001-C003*AB001*D002-C009*AB013*D001+C011*AB SPDSTV $002*D001+C051*AB005*D005-C006*AB002*D007+C053*AB002*D004-C004*AB00 SPDSTV $6*D002+C005*AB001*D002+C049*AB004*D006-C001*AB001*D016+C050*AB001* SPDSTV $D013) SPDSTV EP71=OF5*(-C009*AB014*D001+C011*AB007*D001+C007*AB008*D002+C007*AB SPDSTV $003*D008-C008*AB001*D008-C006*AB002*D009-C004*AB008*D002+C002*AB00 SPDSTV $7*D003-C055*AB007*D004+C002*AB002*D009-C001*AB001*D017+C050*AB001* SPDSTV $D018) SPDSTV EP81=OF5*(-C009*AB015*D001+C007*AB005*D008+C007*AB008*D005-C006*AB SPDSTV $002*D010-C004*AB009*D002+C002*AB004*D009+C002*AB007*D006-C001*AB00 SPDSTV $1*D019) SPDSTV EP91=OF5*(-C008*AB002*D001-C003*AB001*D002-C009*AB016*D001+C011*AB SPDSTV $002*D001+C051*AB008*D008-C006*AB002*D011+C053*AB002*D004-C004*AB01 SPDSTV $0*D002+C005*AB001*D002+C049*AB007*D009-C001*AB001*D020+C050*AB001* SPDSTV $D013) SPDSTV EP22=OF8*(+C017*AB001*D001+C018*AB003*D001-C019*AB001*D001+C057*AB SPDSTV $002*D002+C003*AB001*D003-C058*AB001*D004+C011*AB003*D001-C020*AB00 SPDSTV $1*D001+C012*AB001*D003-C059*AB001*D004+C021*AB017*D001-C022*AB003* SPDSTV $D001-C060*AB011*D002+C023*AB001*D001+C038*AB002*D002+C013*AB003*D0 SPDSTV $03-C061*AB003*D004-C014*AB001*D003+C062*AB001*D004+C063*AB011*D002 SPDSTV $-C033*AB002*D002-C064*AB003*D003+C051*AB003*D004+C031*AB001*D003-C SPDSTV $052*AB001*D004+C056*AB002*D012-C065*AB002*D013+C004*AB003*D003-C00 SPDSTV $5*AB001*D003-C049*AB002*D012+C066*AB002*D013+C001*AB001*D021-C067* SPDSTV $AB001*D022+C068*AB001*D023-C069*AB003*D004+C070*AB001*D004) SPDSTV EP42=OF8*(+C011*AB005*D001-C008*AB004*D002-C008*AB002*D005+C012*AB SPDSTV $001*D006+C021*AB018*D001-C024*AB005*D001-C015*AB012*D002-C015*AB01 SPDSTV $1*D005+C016*AB002*D005+C013*AB003*D006+C018*AB004*D002-C014*AB001* SPDSTV $D006+C063*AB012*D002-C071*AB004*D002-C051*AB005*D003+C007*AB005*D0 SPDSTV $04-C051*AB003*D006+C052*AB001*D006+C056*AB002*D014-C006*AB002*D015 SPDSTV $+C004*AB005*D003-C002*AB004*D012+C072*AB004*D013-C002*AB002*D014+C SPDSTV $001*AB001*D024-C054*AB001*D025-C069*AB005*D004+C055*AB002*D015) SPDSTV EP52=OF8*(+C017*AB001*D001+C018*AB003*D001-C019*AB001*D001+C052*AB SPDSTV $002*D002+C003*AB001*D003-C058*AB001*D004+C011*AB006*D001-C020*AB00 SPDSTV $1*D001-C052*AB004*D005+C012*AB001*D007-C059*AB001*D004+C021*AB019* SPDSTV $D001-C025*AB003*D001-C060*AB012*D005+C013*AB003*D007-C061*AB003*D0 SPDSTV $04-C025*AB006*D001+C026*AB001*D001+C073*AB004*D005-C014*AB001*D007 SPDSTV $+C062*AB001*D004+C063*AB013*D002-C071*AB002*D002-C064*AB005*D006+C SPDSTV $056*AB002*D016-C006*AB002*D013+C004*AB006*D003-C005*AB001*D003-C04 SPDSTV $9*AB004*D014+C001*AB001*D026-C050*AB001*D022-C069*AB006*D004+C070* SPDSTV $AB001*D004+C002*AB004*D015-C050*AB001*D027+C074*AB001*D023) SPDSTV EP72=OF8*(+C011*AB008*D001-C008*AB007*D002-C008*AB002*D008+C012*AB SPDSTV $001*D009+C021*AB020*D001-C024*AB008*D001-C015*AB014*D002-C015*AB01 SPDSTV $1*D008+C016*AB002*D008+C013*AB003*D009+C018*AB007*D002-C014*AB001* SPDSTV $D009+C063*AB014*D002-C071*AB007*D002-C051*AB008*D003+C007*AB008*D0 SPDSTV $04-C051*AB003*D009+C052*AB001*D009+C056*AB002*D017-C006*AB002*D018 SPDSTV $+C004*AB008*D003-C002*AB007*D012+C072*AB007*D013-C002*AB002*D017+C SPDSTV $001*AB001*D028-C054*AB001*D029-C069*AB008*D004+C055*AB002*D018) SPDSTV EP82=OF8*(+C011*AB009*D001-C008*AB004*D008-C008*AB007*D005+C012*AB SPDSTV $001*D010+C021*AB021*D001-C015*AB012*D008-C015*AB014*D005+C013*AB00 SPDSTV $3*D010-C025*AB009*D001+C018*AB004*D008+C018*AB007*D005-C014*AB001* SPDSTV $D010+C063*AB015*D002-C051*AB005*D009-C051*AB008*D006+C056*AB002*D0 SPDSTV $19+C004*AB009*D003-C002*AB004*D017-C002*AB007*D014+C001*AB001*D030 SPDSTV $-C069*AB009*D004+C055*AB004*D018+C055*AB007*D015-C050*AB001*D031) SPDSTV EP92=OF8*(+C017*AB001*D001+C018*AB003*D001-C019*AB001*D001+C052*AB SPDSTV $002*D002+C003*AB001*D003-C058*AB001*D004+C011*AB010*D001-C020*AB00 SPDSTV $1*D001-C052*AB007*D008+C012*AB001*D011-C059*AB001*D004+C021*AB022* SPDSTV $D001-C025*AB003*D001-C060*AB014*D008+C013*AB003*D011-C061*AB003*D0 SPDSTV $04-C025*AB010*D001+C026*AB001*D001+C073*AB007*D008-C014*AB001*D011 SPDSTV $+C062*AB001*D004+C063*AB016*D002-C071*AB002*D002-C064*AB008*D009+C SPDSTV $056*AB002*D020-C006*AB002*D013+C004*AB010*D003-C005*AB001*D003-C04 SPDSTV $9*AB007*D017+C001*AB001*D032-C050*AB001*D022-C069*AB010*D004+C070* SPDSTV $AB001*D004+C002*AB007*D018-C050*AB001*D033+C074*AB001*D023) SPDSTV EP23=OF5*(-C008*AB004*D001-C003*AB001*D005-C009*AB012*D001+C011*AB SPDSTV $004*D001+C051*AB005*D002-C006*AB004*D003+C053*AB004*D004-C004*AB00 SPDSTV $3*D005+C005*AB001*D005+C049*AB002*D006-C001*AB001*D014+C050*AB001* SPDSTV $D015) SPDSTV EP43=OF5*(-C009*AB013*D001+C007*AB006*D002+C011*AB002*D001-C008*AB SPDSTV $001*D002+C007*AB005*D005-C006*AB004*D006-C004*AB005*D005+C002*AB00 SPDSTV $4*D006+C002*AB002*D007-C001*AB001*D016-C055*AB002*D004+C050*AB001* SPDSTV $D013) SPDSTV EP53=OF5*(-C008*AB004*D001-C003*AB001*D005-C009*AB023*D001+C010*AB SPDSTV $004*D001+C051*AB006*D005-C052*AB001*D005-C006*AB004*D007+C053*AB00 SPDSTV $4*D004-C004*AB006*D005+C005*AB001*D005+C049*AB004*D007-C002*AB004* SPDSTV $D004-C001*AB001*D034+C054*AB001*D015) SPDSTV EP73=OF5*(-C009*AB015*D001+C007*AB009*D002+C007*AB005*D008-C006*AB SPDSTV $004*D009-C004*AB008*D005+C002*AB007*D006+C002*AB002*D010-C001*AB00 SPDSTV $1*D019) SPDSTV EP83=OF5*(-C009*AB024*D001+C007*AB006*D008+C011*AB007*D001-C008*AB SPDSTV $001*D008+C007*AB009*D005-C006*AB004*D010-C004*AB009*D005+C002*AB00 SPDSTV $4*D010+C002*AB007*D007-C001*AB001*D035-C055*AB007*D004+C050*AB001* SPDSTV $D018) SPDSTV EP93=OF5*(-C008*AB004*D001-C003*AB001*D005-C009*AB025*D001+C011*AB SPDSTV $004*D001+C051*AB009*D008-C006*AB004*D011+C053*AB004*D004-C004*AB01 SPDSTV $0*D005+C005*AB001*D005+C049*AB007*D010-C001*AB001*D036+C050*AB001* SPDSTV $D015) SPDSTV EP24=OF8*(+C018*AB005*D001+C008*AB004*D002+C008*AB002*D005+C003*AB SPDSTV $001*D006+C021*AB018*D001-C024*AB005*D001-C060*AB012*D002+C073*AB00 SPDSTV $4*D002+C013*AB005*D003-C061*AB005*D004+C009*AB012*D002-C011*AB004* SPDSTV $D002-C051*AB005*D003+C007*AB005*D004+C006*AB004*D012-C075*AB004*D0 SPDSTV $13+C009*AB011*D005-C010*AB002*D005-C051*AB003*D006+C052*AB001*D006 SPDSTV $+C006*AB002*D014-C053*AB002*D015+C004*AB003*D006-C005*AB001*D006-C SPDSTV $049*AB002*D014+C002*AB002*D015+C001*AB001*D024-C054*AB001*D025) SPDSTV EP44=OF8*(+C021*AB019*D001-C025*AB006*D001-C015*AB013*D002-C025*AB SPDSTV $003*D001+C026*AB001*D001+C018*AB002*D002-C015*AB012*D005+C018*AB00 SPDSTV $4*D005+C013*AB005*D006+C009*AB013*D002-C007*AB006*D003+C076*AB006* SPDSTV $D004-C011*AB002*D002+C008*AB001*D003-C008*AB001*D004-C051*AB005*D0 SPDSTV $06+C006*AB004*D014-C053*AB004*D015+C009*AB012*D005-C011*AB004*D005 SPDSTV $-C007*AB003*D007+C008*AB001*D007+C006*AB002*D016+C076*AB003*D004-C SPDSTV $053*AB002*D013+C004*AB005*D006-C002*AB004*D014+C055*AB004*D015-C00 SPDSTV $2*AB002*D016+C001*AB001*D026-C050*AB001*D027+C055*AB002*D013-C050* SPDSTV $AB001*D022+C074*AB001*D023) SPDSTV EP54=OF8*(+C018*AB005*D001+C008*AB004*D002+C008*AB002*D005+C003*AB SPDSTV $001*D006+C021*AB026*D001-C024*AB005*D001-C060*AB013*D005+C073*AB00 SPDSTV $2*D005+C013*AB005*D007-C061*AB005*D004+C009*AB023*D002-C010*AB004* SPDSTV $D002-C051*AB006*D006+C052*AB001*D006+C006*AB004*D016-C053*AB004*D0 SPDSTV $13+C009*AB013*D005-C011*AB002*D005-C051*AB005*D007+C007*AB005*D004 SPDSTV $+C006*AB002*D034-C075*AB002*D015+C004*AB006*D006-C005*AB001*D006-C SPDSTV $049*AB004*D016+C002*AB004*D013+C001*AB001*D037-C054*AB001*D025) SPDSTV EP74=OF8*(+C021*AB021*D001-C025*AB009*D001-C015*AB015*D002-C015*AB SPDSTV $012*D008+C018*AB004*D008+C013*AB005*D009+C009*AB015*D002-C007*AB00 SPDSTV $9*D003+C076*AB009*D004-C007*AB005*D009+C006*AB004*D017-C053*AB004* SPDSTV $D018+C009*AB014*D005-C011*AB007*D005-C007*AB008*D006-C007*AB003*D0 SPDSTV $10+C008*AB001*D010+C006*AB002*D019+C004*AB008*D006-C002*AB007*D014 SPDSTV $+C055*AB007*D015-C002*AB002*D019+C001*AB001*D030-C050*AB001*D031) SPDSTV EP84=OF8*(+C021*AB027*D001-C015*AB013*D008-C025*AB008*D001+C018*AB SPDSTV $002*D008-C015*AB015*D005+C013*AB005*D010+C009*AB024*D002-C007*AB00 SPDSTV $6*D009-C011*AB007*D002+C008*AB001*D009-C007*AB009*D006+C006*AB004* SPDSTV $D019+C009*AB015*D005-C007*AB005*D010-C007*AB008*D007+C006*AB002*D0 SPDSTV $35+C076*AB008*D004-C053*AB002*D018+C004*AB009*D006-C002*AB004*D019 SPDSTV $-C002*AB007*D016+C001*AB001*D038+C055*AB007*D013-C050*AB001*D029) SPDSTV EP94=OF8*(+C018*AB005*D001+C008*AB004*D002+C008*AB002*D005+C003*AB SPDSTV $001*D006+C021*AB028*D001-C025*AB005*D001-C060*AB015*D008+C013*AB00 SPDSTV $5*D011-C061*AB005*D004+C009*AB025*D002-C011*AB004*D002-C051*AB009* SPDSTV $D009+C006*AB004*D020-C053*AB004*D013+C009*AB016*D005-C011*AB002*D0 SPDSTV $05-C051*AB008*D010+C006*AB002*D036-C053*AB002*D015+C004*AB010*D006 SPDSTV $-C005*AB001*D006-C049*AB007*D019+C001*AB001*D039-C050*AB001*D025) SPDSTV EP25=OF8*(+C017*AB001*D001+C018*AB006*D001-C019*AB001*D001+C052*AB SPDSTV $004*D005+C003*AB001*D007-C058*AB001*D004+C011*AB003*D001-C020*AB00 SPDSTV $1*D001-C052*AB002*D002+C012*AB001*D003-C059*AB001*D004+C021*AB019* SPDSTV $D001-C025*AB006*D001-C060*AB013*D002+C013*AB006*D003-C061*AB006*D0 SPDSTV $04-C025*AB003*D001+C026*AB001*D001+C073*AB002*D002-C014*AB001*D003 SPDSTV $+C062*AB001*D004+C063*AB012*D005-C071*AB004*D005-C064*AB005*D006+C SPDSTV $056*AB004*D014-C006*AB004*D015+C004*AB003*D007-C005*AB001*D007-C04 SPDSTV $9*AB002*D016+C001*AB001*D026-C050*AB001*D027-C069*AB003*D004+C070* SPDSTV $AB001*D004+C002*AB002*D013-C050*AB001*D022+C074*AB001*D023) SPDSTV EP45=OF8*(+C011*AB005*D001-C008*AB004*D002-C008*AB002*D005+C012*AB SPDSTV $001*D006+C021*AB026*D001-C015*AB023*D002-C024*AB005*D001+C016*AB00 SPDSTV $4*D002-C015*AB013*D005+C013*AB006*D006+C018*AB002*D005-C014*AB001* SPDSTV $D006+C063*AB013*D005-C051*AB006*D006-C071*AB002*D005+C052*AB001*D0 SPDSTV $06-C051*AB005*D007+C056*AB004*D016+C007*AB005*D004-C006*AB004*D013 SPDSTV $+C004*AB005*D007-C002*AB004*D016-C002*AB002*D034+C001*AB001*D037+C SPDSTV $072*AB002*D015-C054*AB001*D025-C069*AB005*D004+C055*AB004*D013) SPDSTV EP55=OF8*(+C017*AB001*D001+C018*AB006*D001-C019*AB001*D001+C057*AB SPDSTV $004*D005+C003*AB001*D007-C058*AB001*D004+C011*AB006*D001-C020*AB00 SPDSTV $1*D001+C012*AB001*D007-C059*AB001*D004+C021*AB029*D001-C022*AB006* SPDSTV $D001-C060*AB023*D005+C023*AB001*D001+C038*AB004*D005+C013*AB006*D0 SPDSTV $07-C061*AB006*D004-C014*AB001*D007+C062*AB001*D004+C063*AB023*D005 SPDSTV $-C033*AB004*D005-C064*AB006*D007+C051*AB006*D004+C031*AB001*D007-C SPDSTV $052*AB001*D004+C056*AB004*D034-C065*AB004*D015+C004*AB006*D007-C00 SPDSTV $5*AB001*D007-C049*AB004*D034+C066*AB004*D015+C001*AB001*D040-C067* SPDSTV $AB001*D027+C068*AB001*D023-C069*AB006*D004+C070*AB001*D004) SPDSTV EP75=OF8*(+C011*AB008*D001-C008*AB007*D002-C008*AB002*D008+C012*AB SPDSTV $001*D009+C021*AB027*D001-C015*AB024*D002-C015*AB013*D008+C013*AB00 SPDSTV $6*D009-C025*AB008*D001+C018*AB007*D002+C018*AB002*D008-C014*AB001* SPDSTV $D009+C063*AB015*D005-C051*AB009*D006-C051*AB005*D010+C056*AB004*D0 SPDSTV $19+C004*AB008*D007-C002*AB007*D016-C002*AB002*D035+C001*AB001*D038 SPDSTV $-C069*AB008*D004+C055*AB007*D013+C055*AB002*D018-C050*AB001*D029) SPDSTV EP85=OF8*(+C011*AB009*D001-C008*AB004*D008-C008*AB007*D005+C012*AB SPDSTV $001*D010+C021*AB030*D001-C015*AB023*D008-C024*AB009*D001+C016*AB00 SPDSTV $4*D008-C015*AB024*D005+C013*AB006*D010+C018*AB007*D005-C014*AB001* SPDSTV $D010+C063*AB024*D005-C051*AB006*D010-C071*AB007*D005+C052*AB001*D0 SPDSTV $10-C051*AB009*D007+C056*AB004*D035+C007*AB009*D004-C006*AB004*D018 SPDSTV $+C004*AB009*D007-C002*AB004*D035-C002*AB007*D034+C001*AB001*D041+C SPDSTV $072*AB007*D015-C054*AB001*D031-C069*AB009*D004+C055*AB004*D018) SPDSTV EP95=OF8*(+C017*AB001*D001+C018*AB006*D001-C019*AB001*D001+C052*AB SPDSTV $004*D005+C003*AB001*D007-C058*AB001*D004+C011*AB010*D001-C020*AB00 SPDSTV $1*D001-C052*AB007*D008+C012*AB001*D011-C059*AB001*D004+C021*AB031* SPDSTV $D001-C025*AB006*D001-C060*AB024*D008+C013*AB006*D011-C061*AB006*D0 SPDSTV $04-C025*AB010*D001+C026*AB001*D001+C073*AB007*D008-C014*AB001*D011 SPDSTV $+C062*AB001*D004+C063*AB025*D005-C071*AB004*D005-C064*AB009*D010+C SPDSTV $056*AB004*D036-C006*AB004*D015+C004*AB010*D007-C005*AB001*D007-C04 SPDSTV $9*AB007*D035+C001*AB001*D042-C050*AB001*D027-C069*AB010*D004+C070* SPDSTV $AB001*D004+C002*AB007*D018-C050*AB001*D033+C074*AB001*D023) SPDSTV EP26=OF5*(-C008*AB007*D001-C003*AB001*D008-C009*AB014*D001+C011*AB SPDSTV $007*D001+C051*AB008*D002-C006*AB007*D003+C053*AB007*D004-C004*AB00 SPDSTV $3*D008+C005*AB001*D008+C049*AB002*D009-C001*AB001*D017+C050*AB001* SPDSTV $D018) SPDSTV EP46=OF5*(-C009*AB015*D001+C007*AB009*D002+C007*AB008*D005-C006*AB SPDSTV $007*D006-C004*AB005*D008+C002*AB004*D009+C002*AB002*D010-C001*AB00 SPDSTV $1*D019) SPDSTV EP56=OF5*(-C008*AB007*D001-C003*AB001*D008-C009*AB024*D001+C011*AB SPDSTV $007*D001+C051*AB009*D005-C006*AB007*D007+C053*AB007*D004-C004*AB00 SPDSTV $6*D008+C005*AB001*D008+C049*AB004*D010-C001*AB001*D035+C050*AB001* SPDSTV $D018) SPDSTV EP76=OF5*(-C009*AB016*D001+C007*AB010*D002+C011*AB002*D001-C008*AB SPDSTV $001*D002+C007*AB008*D008-C006*AB007*D009-C004*AB008*D008+C002*AB00 SPDSTV $7*D009+C002*AB002*D011-C001*AB001*D020-C055*AB002*D004+C050*AB001* SPDSTV $D013) SPDSTV EP86=OF5*(-C009*AB025*D001+C011*AB004*D001+C007*AB009*D008+C007*AB SPDSTV $010*D005-C008*AB001*D005-C006*AB007*D010-C004*AB009*D008+C002*AB00 SPDSTV $4*D011-C055*AB004*D004+C002*AB007*D010-C001*AB001*D036+C050*AB001* SPDSTV $D015) SPDSTV EP96=OF5*(-C008*AB007*D001-C003*AB001*D008-C009*AB032*D001+C010*AB SPDSTV $007*D001+C051*AB010*D008-C052*AB001*D008-C006*AB007*D011+C053*AB00 SPDSTV $7*D004-C004*AB010*D008+C005*AB001*D008+C049*AB007*D011-C002*AB007* SPDSTV $D004-C001*AB001*D043+C054*AB001*D018) SPDSTV EP27=OF8*(+C018*AB008*D001+C008*AB007*D002+C008*AB002*D008+C003*AB SPDSTV $001*D009+C021*AB020*D001-C024*AB008*D001-C060*AB014*D002+C073*AB00 SPDSTV $7*D002+C013*AB008*D003-C061*AB008*D004+C009*AB014*D002-C011*AB007* SPDSTV $D002-C051*AB008*D003+C007*AB008*D004+C006*AB007*D012-C075*AB007*D0 SPDSTV $13+C009*AB011*D008-C010*AB002*D008-C051*AB003*D009+C052*AB001*D009 SPDSTV $+C006*AB002*D017-C053*AB002*D018+C004*AB003*D009-C005*AB001*D009-C SPDSTV $049*AB002*D017+C002*AB002*D018+C001*AB001*D028-C054*AB001*D029) SPDSTV EP47=OF8*(+C021*AB021*D001-C025*AB009*D001-C015*AB015*D002-C015*AB SPDSTV $014*D005+C018*AB007*D005+C013*AB008*D006+C009*AB015*D002-C007*AB00 SPDSTV $9*D003+C076*AB009*D004-C007*AB008*D006+C006*AB007*D014-C053*AB007* SPDSTV $D015+C009*AB012*D008-C011*AB004*D008-C007*AB005*D009-C007*AB003*D0 SPDSTV $10+C008*AB001*D010+C006*AB002*D019+C004*AB005*D009-C002*AB004*D017 SPDSTV $+C055*AB004*D018-C002*AB002*D019+C001*AB001*D030-C050*AB001*D031) SPDSTV EP57=OF8*(+C018*AB008*D001+C008*AB007*D002+C008*AB002*D008+C003*AB SPDSTV $001*D009+C021*AB027*D001-C025*AB008*D001-C060*AB015*D005+C013*AB00 SPDSTV $8*D007-C061*AB008*D004+C009*AB024*D002-C011*AB007*D002-C051*AB009* SPDSTV $D006+C006*AB007*D016-C053*AB007*D013+C009*AB013*D008-C011*AB002*D0 SPDSTV $08-C051*AB005*D010+C006*AB002*D035-C053*AB002*D018+C004*AB006*D009 SPDSTV $-C005*AB001*D009-C049*AB004*D019+C001*AB001*D038-C050*AB001*D029) SPDSTV EP77=OF8*(+C021*AB022*D001-C025*AB010*D001-C015*AB016*D002-C025*AB SPDSTV $003*D001+C026*AB001*D001+C018*AB002*D002-C015*AB014*D008+C018*AB00 SPDSTV $7*D008+C013*AB008*D009+C009*AB016*D002-C007*AB010*D003+C076*AB010* SPDSTV $D004-C011*AB002*D002+C008*AB001*D003-C008*AB001*D004-C051*AB008*D0 SPDSTV $09+C006*AB007*D017-C053*AB007*D018+C009*AB014*D008-C011*AB007*D008 SPDSTV $-C007*AB003*D011+C008*AB001*D011+C006*AB002*D020+C076*AB003*D004-C SPDSTV $053*AB002*D013+C004*AB008*D009-C002*AB007*D017+C055*AB007*D018-C00 SPDSTV $2*AB002*D020+C001*AB001*D032-C050*AB001*D033+C055*AB002*D013-C050* SPDSTV $AB001*D022+C074*AB001*D023) SPDSTV EP87=OF8*(+C021*AB028*D001-C025*AB005*D001-C015*AB015*D008-C015*AB SPDSTV $016*D005+C018*AB002*D005+C013*AB008*D010+C009*AB025*D002-C011*AB00 SPDSTV $4*D002-C007*AB009*D009-C007*AB010*D006+C008*AB001*D006+C006*AB007* SPDSTV $D019+C009*AB015*D008-C007*AB005*D011+C076*AB005*D004-C007*AB008*D0 SPDSTV $10+C006*AB002*D036-C053*AB002*D015+C004*AB009*D009-C002*AB004*D020 SPDSTV $+C055*AB004*D013-C002*AB007*D019+C001*AB001*D039-C050*AB001*D025) SPDSTV EP97=OF8*(+C018*AB008*D001+C008*AB007*D002+C008*AB002*D008+C003*AB SPDSTV $001*D009+C021*AB033*D001-C024*AB008*D001-C060*AB016*D008+C073*AB00 SPDSTV $2*D008+C013*AB008*D011-C061*AB008*D004+C009*AB032*D002-C010*AB007* SPDSTV $D002-C051*AB010*D009+C052*AB001*D009+C006*AB007*D020-C053*AB007*D0 SPDSTV $13+C009*AB016*D008-C011*AB002*D008-C051*AB008*D011+C007*AB008*D004 SPDSTV $+C006*AB002*D043-C075*AB002*D018+C004*AB010*D009-C005*AB001*D009-C SPDSTV $049*AB007*D020+C002*AB007*D013+C001*AB001*D044-C054*AB001*D029) SPDSTV EP28=OF8*(+C018*AB009*D001+C008*AB004*D008+C008*AB007*D005+C003*AB SPDSTV $001*D010+C021*AB021*D001-C025*AB009*D001-C060*AB015*D002+C013*AB00 SPDSTV $9*D003-C061*AB009*D004+C009*AB012*D008-C011*AB004*D008-C051*AB005* SPDSTV $D009+C006*AB004*D017-C053*AB004*D018+C009*AB014*D005-C011*AB007*D0 SPDSTV $05-C051*AB008*D006+C006*AB007*D014-C053*AB007*D015+C004*AB003*D010 SPDSTV $-C005*AB001*D010-C049*AB002*D019+C001*AB001*D030-C050*AB001*D031) SPDSTV EP48=OF8*(+C021*AB027*D001-C015*AB024*D002-C025*AB008*D001+C018*AB SPDSTV $007*D002-C015*AB015*D005+C013*AB009*D006+C009*AB013*D008-C007*AB00 SPDSTV $6*D009-C011*AB002*D008+C008*AB001*D009-C007*AB005*D010+C006*AB004* SPDSTV $D019+C009*AB015*D005-C007*AB009*D006-C007*AB008*D007+C006*AB007*D0 SPDSTV $16+C076*AB008*D004-C053*AB007*D013+C004*AB005*D010-C002*AB004*D019 SPDSTV $-C002*AB002*D035+C001*AB001*D038+C055*AB002*D018-C050*AB001*D029) SPDSTV EP58=OF8*(+C018*AB009*D001+C008*AB004*D008+C008*AB007*D005+C003*AB SPDSTV $001*D010+C021*AB030*D001-C024*AB009*D001-C060*AB024*D005+C073*AB00 SPDSTV $7*D005+C013*AB009*D007-C061*AB009*D004+C009*AB023*D008-C010*AB004* SPDSTV $D008-C051*AB006*D010+C052*AB001*D010+C006*AB004*D035-C053*AB004*D0 SPDSTV $18+C009*AB024*D005-C011*AB007*D005-C051*AB009*D007+C007*AB009*D004 SPDSTV $+C006*AB007*D034-C075*AB007*D015+C004*AB006*D010-C005*AB001*D010-C SPDSTV $049*AB004*D035+C002*AB004*D018+C001*AB001*D041-C054*AB001*D031) SPDSTV EP78=OF8*(+C021*AB028*D001-C015*AB025*D002-C025*AB005*D001+C018*AB SPDSTV $004*D002-C015*AB015*D008+C013*AB009*D009+C009*AB015*D008-C007*AB00 SPDSTV $9*D009-C007*AB005*D011+C006*AB004*D020+C076*AB005*D004-C053*AB004* SPDSTV $D013+C009*AB016*D005-C007*AB010*D006-C011*AB002*D005+C008*AB001*D0 SPDSTV $06-C007*AB008*D010+C006*AB007*D019+C004*AB008*D010-C002*AB007*D019 SPDSTV $-C002*AB002*D036+C001*AB001*D039+C055*AB002*D015-C050*AB001*D025) SPDSTV EP88=OF8*(+C021*AB031*D001-C025*AB006*D001-C015*AB024*D008-C025*AB SPDSTV $010*D001+C026*AB001*D001+C018*AB007*D008-C015*AB025*D005+C018*AB00 SPDSTV $4*D005+C013*AB009*D010+C009*AB024*D008-C007*AB006*D011+C076*AB006* SPDSTV $D004-C011*AB007*D008+C008*AB001*D011-C008*AB001*D004-C051*AB009*D0 SPDSTV $10+C006*AB004*D036-C053*AB004*D015+C009*AB025*D005-C011*AB004*D005 SPDSTV $-C007*AB010*D007+C008*AB001*D007+C006*AB007*D035+C076*AB010*D004-C SPDSTV $053*AB007*D018+C004*AB009*D010-C002*AB004*D036+C055*AB004*D015-C00 SPDSTV $2*AB007*D035+C001*AB001*D042-C050*AB001*D027+C055*AB007*D018-C050* SPDSTV $AB001*D033+C074*AB001*D023) SPDSTV EP98=OF8*(+C018*AB009*D001+C008*AB004*D008+C008*AB007*D005+C003*AB SPDSTV $001*D010+C021*AB034*D001-C024*AB009*D001-C060*AB025*D008+C073*AB00 SPDSTV $4*D008+C013*AB009*D011-C061*AB009*D004+C009*AB025*D008-C011*AB004* SPDSTV $D008-C051*AB009*D011+C007*AB009*D004+C006*AB004*D043-C075*AB004*D0 SPDSTV $18+C009*AB032*D005-C010*AB007*D005-C051*AB010*D010+C052*AB001*D010 SPDSTV $+C006*AB007*D036-C053*AB007*D015+C004*AB010*D010-C005*AB001*D010-C SPDSTV $049*AB007*D036+C002*AB007*D015+C001*AB001*D045-C054*AB001*D031) SPDSTV EP29=OF8*(+C017*AB001*D001+C018*AB010*D001-C019*AB001*D001+C052*AB SPDSTV $007*D008+C003*AB001*D011-C058*AB001*D004+C011*AB003*D001-C020*AB00 SPDSTV $1*D001-C052*AB002*D002+C012*AB001*D003-C059*AB001*D004+C021*AB022* SPDSTV $D001-C025*AB010*D001-C060*AB016*D002+C013*AB010*D003-C061*AB010*D0 SPDSTV $04-C025*AB003*D001+C026*AB001*D001+C073*AB002*D002-C014*AB001*D003 SPDSTV $+C062*AB001*D004+C063*AB014*D008-C071*AB007*D008-C064*AB008*D009+C SPDSTV $056*AB007*D017-C006*AB007*D018+C004*AB003*D011-C005*AB001*D011-C04 SPDSTV $9*AB002*D020+C001*AB001*D032-C050*AB001*D033-C069*AB003*D004+C070* SPDSTV $AB001*D004+C002*AB002*D013-C050*AB001*D022+C074*AB001*D023) SPDSTV EP49=OF8*(+C011*AB005*D001-C008*AB004*D002-C008*AB002*D005+C012*AB SPDSTV $001*D006+C021*AB028*D001-C015*AB025*D002-C015*AB016*D005+C013*AB01 SPDSTV $0*D006-C025*AB005*D001+C018*AB004*D002+C018*AB002*D005-C014*AB001* SPDSTV $D006+C063*AB015*D008-C051*AB009*D009-C051*AB008*D010+C056*AB007*D0 SPDSTV $19+C004*AB005*D011-C002*AB004*D020-C002*AB002*D036+C001*AB001*D039 SPDSTV $-C069*AB005*D004+C055*AB004*D013+C055*AB002*D015-C050*AB001*D025) SPDSTV EP59=OF8*(+C017*AB001*D001+C018*AB010*D001-C019*AB001*D001+C052*AB SPDSTV $007*D008+C003*AB001*D011-C058*AB001*D004+C011*AB006*D001-C020*AB00 SPDSTV $1*D001-C052*AB004*D005+C012*AB001*D007-C059*AB001*D004+C021*AB031* SPDSTV $D001-C025*AB010*D001-C060*AB025*D005+C013*AB010*D007-C061*AB010*D0 SPDSTV $04-C025*AB006*D001+C026*AB001*D001+C073*AB004*D005-C014*AB001*D007 SPDSTV $+C062*AB001*D004+C063*AB024*D008-C071*AB007*D008-C064*AB009*D010+C SPDSTV $056*AB007*D035-C006*AB007*D018+C004*AB006*D011-C005*AB001*D011-C04 SPDSTV $9*AB004*D036+C001*AB001*D042-C050*AB001*D033-C069*AB006*D004+C070* SPDSTV $AB001*D004+C002*AB004*D015-C050*AB001*D027+C074*AB001*D023) SPDSTV EP79=OF8*(+C011*AB008*D001-C008*AB007*D002-C008*AB002*D008+C012*AB SPDSTV $001*D009+C021*AB033*D001-C015*AB032*D002-C024*AB008*D001+C016*AB00 SPDSTV $7*D002-C015*AB016*D008+C013*AB010*D009+C018*AB002*D008-C014*AB001* SPDSTV $D009+C063*AB016*D008-C051*AB010*D009-C071*AB002*D008+C052*AB001*D0 SPDSTV $09-C051*AB008*D011+C056*AB007*D020+C007*AB008*D004-C006*AB007*D013 SPDSTV $+C004*AB008*D011-C002*AB007*D020-C002*AB002*D043+C001*AB001*D044+C SPDSTV $072*AB002*D018-C054*AB001*D029-C069*AB008*D004+C055*AB007*D013) SPDSTV EP89=OF8*(+C011*AB009*D001-C008*AB004*D008-C008*AB007*D005+C012*AB SPDSTV $001*D010+C021*AB034*D001-C024*AB009*D001-C015*AB025*D008-C015*AB03 SPDSTV $2*D005+C016*AB007*D005+C013*AB010*D010+C018*AB004*D008-C014*AB001* SPDSTV $D010+C063*AB025*D008-C071*AB004*D008-C051*AB009*D011+C007*AB009*D0 SPDSTV $04-C051*AB010*D010+C052*AB001*D010+C056*AB007*D036-C006*AB007*D015 SPDSTV $+C004*AB009*D011-C002*AB004*D043+C072*AB004*D018-C002*AB007*D036+C SPDSTV $001*AB001*D045-C054*AB001*D031-C069*AB009*D004+C055*AB007*D015) SPDSTV EP99=OF8*(+C017*AB001*D001+C018*AB010*D001-C019*AB001*D001+C057*AB SPDSTV $007*D008+C003*AB001*D011-C058*AB001*D004+C011*AB010*D001-C020*AB00 SPDSTV $1*D001+C012*AB001*D011-C059*AB001*D004+C021*AB035*D001-C022*AB010* SPDSTV $D001-C060*AB032*D008+C023*AB001*D001+C038*AB007*D008+C013*AB010*D0 SPDSTV $11-C061*AB010*D004-C014*AB001*D011+C062*AB001*D004+C063*AB032*D008 SPDSTV $-C033*AB007*D008-C064*AB010*D011+C051*AB010*D004+C031*AB001*D011-C SPDSTV $052*AB001*D004+C056*AB007*D043-C065*AB007*D018+C004*AB010*D011-C00 SPDSTV $5*AB001*D011-C049*AB007*D043+C066*AB007*D018+C001*AB001*D046-C067* SPDSTV $AB001*D033+C068*AB001*D023-C069*AB010*D004+C070*AB001*D004) SPDSTV C ****************************************************************** SPDSTV C * PD * SPDSTV C ****************************************************************** SPDSTV 3242 CONTINUE SPDSTV EP12=OF7*(+C008*AB002*D001-C012*AB001*D002+C015*AB011*D001-C016*AB SPDSTV $002*D001-C013*AB003*D002+C014*AB001*D002+C051*AB003*D002-C052*AB00 SPDSTV $1*D002-C056*AB002*D003+C006*AB002*D004+C002*AB002*D003-C001*AB001* SPDSTV $D012+C054*AB001*D013-C055*AB002*D004) SPDSTV EP32=OF7*(+C008*AB004*D001-C012*AB001*D005+C015*AB012*D001-C013*AB SPDSTV $003*D005-C018*AB004*D001+C014*AB001*D005+C051*AB005*D002-C056*AB00 SPDSTV $2*D006+C002*AB004*D003-C001*AB001*D014-C055*AB004*D004+C050*AB001* SPDSTV $D015) SPDSTV EP62=OF7*(+C008*AB007*D001-C012*AB001*D008+C015*AB014*D001-C013*AB SPDSTV $003*D008-C018*AB007*D001+C014*AB001*D008+C051*AB008*D002-C056*AB00 SPDSTV $2*D009+C002*AB007*D003-C001*AB001*D017-C055*AB007*D004+C050*AB001* SPDSTV $D018) SPDSTV EP14=OF7*(+C015*AB012*D001-C018*AB004*D001-C013*AB005*D002+C007*AB SPDSTV $005*D002-C006*AB004*D003+C053*AB004*D004+C007*AB003*D005-C008*AB00 SPDSTV $1*D005-C006*AB002*D006+C002*AB002*D006-C001*AB001*D014+C050*AB001* SPDSTV $D015) SPDSTV EP34=OF7*(+C015*AB013*D001-C018*AB002*D001-C013*AB005*D005+C007*AB SPDSTV $006*D002-C008*AB001*D002-C006*AB004*D006+C007*AB005*D005-C006*AB00 SPDSTV $2*D007+C053*AB002*D004+C002*AB004*D006-C001*AB001*D016+C050*AB001* SPDSTV $D013) SPDSTV EP64=OF7*(+C015*AB015*D001-C013*AB005*D008+C007*AB009*D002-C006*AB SPDSTV $004*D009+C007*AB008*D005-C006*AB002*D010+C002*AB007*D006-C001*AB00 SPDSTV $1*D019) SPDSTV EP15=OF7*(+C008*AB002*D001-C012*AB001*D002+C015*AB013*D001-C013*AB SPDSTV $006*D002-C018*AB002*D001+C014*AB001*D002+C051*AB005*D005-C056*AB00 SPDSTV $4*D006+C002*AB002*D007-C001*AB001*D016-C055*AB002*D004+C050*AB001* SPDSTV $D013) SPDSTV EP35=OF7*(+C008*AB004*D001-C012*AB001*D005+C015*AB023*D001-C016*AB SPDSTV $004*D001-C013*AB006*D005+C014*AB001*D005+C051*AB006*D005-C052*AB00 SPDSTV $1*D005-C056*AB004*D007+C006*AB004*D004+C002*AB004*D007-C001*AB001* SPDSTV $D034+C054*AB001*D015-C055*AB004*D004) SPDSTV EP65=OF7*(+C008*AB007*D001-C012*AB001*D008+C015*AB024*D001-C013*AB SPDSTV $006*D008-C018*AB007*D001+C014*AB001*D008+C051*AB009*D005-C056*AB00 SPDSTV $4*D010+C002*AB007*D007-C001*AB001*D035-C055*AB007*D004+C050*AB001* SPDSTV $D018) SPDSTV EP17=OF7*(+C015*AB014*D001-C018*AB007*D001-C013*AB008*D002+C007*AB SPDSTV $008*D002-C006*AB007*D003+C053*AB007*D004+C007*AB003*D008-C008*AB00 SPDSTV $1*D008-C006*AB002*D009+C002*AB002*D009-C001*AB001*D017+C050*AB001* SPDSTV $D018) SPDSTV EP37=OF7*(+C015*AB015*D001-C013*AB008*D005+C007*AB009*D002-C006*AB SPDSTV $007*D006+C007*AB005*D008-C006*AB002*D010+C002*AB004*D009-C001*AB00 SPDSTV $1*D019) SPDSTV EP67=OF7*(+C015*AB016*D001-C018*AB002*D001-C013*AB008*D008+C007*AB SPDSTV $010*D002-C008*AB001*D002-C006*AB007*D009+C007*AB008*D008-C006*AB00 SPDSTV $2*D011+C053*AB002*D004+C002*AB007*D009-C001*AB001*D020+C050*AB001* SPDSTV $D013) SPDSTV EP18=OF7*(+C015*AB015*D001-C013*AB009*D002+C007*AB005*D008-C006*AB SPDSTV $004*D009+C007*AB008*D005-C006*AB007*D006+C002*AB002*D010-C001*AB00 SPDSTV $1*D019) SPDSTV EP38=OF7*(+C015*AB024*D001-C018*AB007*D001-C013*AB009*D005+C007*AB SPDSTV $006*D008-C008*AB001*D008-C006*AB004*D010+C007*AB009*D005-C006*AB00 SPDSTV $7*D007+C053*AB007*D004+C002*AB004*D010-C001*AB001*D035+C050*AB001* SPDSTV $D018) SPDSTV EP68=OF7*(+C015*AB025*D001-C018*AB004*D001-C013*AB009*D008+C007*AB SPDSTV $009*D008-C006*AB004*D011+C053*AB004*D004+C007*AB010*D005-C008*AB00 SPDSTV $1*D005-C006*AB007*D010+C002*AB007*D010-C001*AB001*D036+C050*AB001* SPDSTV $D015) SPDSTV EP19=OF7*(+C008*AB002*D001-C012*AB001*D002+C015*AB016*D001-C013*AB SPDSTV $010*D002-C018*AB002*D001+C014*AB001*D002+C051*AB008*D008-C056*AB00 SPDSTV $7*D009+C002*AB002*D011-C001*AB001*D020-C055*AB002*D004+C050*AB001* SPDSTV $D013) SPDSTV EP39=OF7*(+C008*AB004*D001-C012*AB001*D005+C015*AB025*D001-C013*AB SPDSTV $010*D005-C018*AB004*D001+C014*AB001*D005+C051*AB009*D008-C056*AB00 SPDSTV $7*D010+C002*AB004*D011-C001*AB001*D036-C055*AB004*D004+C050*AB001* SPDSTV $D015) SPDSTV EP69=OF7*(+C008*AB007*D001-C012*AB001*D008+C015*AB032*D001-C016*AB SPDSTV $007*D001-C013*AB010*D008+C014*AB001*D008+C051*AB010*D008-C052*AB00 SPDSTV $1*D008-C056*AB007*D011+C006*AB007*D004+C002*AB007*D011-C001*AB001* SPDSTV $D043+C054*AB001*D018-C055*AB007*D004) SPDSTV C ****************************************************************** SPDSTV C * SD * SPDSTV C ****************************************************************** SPDSTV 3260 CONTINUE SPDSTV EP02=OF6*(+C012*AB001*D001+C013*AB003*D001-C014*AB001*D001+C056*AB SPDSTV $002*D002+C001*AB001*D003-C050*AB001*D004) SPDSTV EP04=OF6*(+C013*AB005*D001+C006*AB004*D002+C006*AB002*D005+C001*AB SPDSTV $001*D006) SPDSTV EP05=OF6*(+C012*AB001*D001+C013*AB006*D001-C014*AB001*D001+C056*AB SPDSTV $004*D005+C001*AB001*D007-C050*AB001*D004) SPDSTV EP07=OF6*(+C013*AB008*D001+C006*AB007*D002+C006*AB002*D008+C001*AB SPDSTV $001*D009) SPDSTV EP08=OF6*(+C013*AB009*D001+C006*AB004*D008+C006*AB007*D005+C001*AB SPDSTV $001*D010) SPDSTV EP09=OF6*(+C012*AB001*D001+C013*AB010*D001-C014*AB001*D001+C056*AB SPDSTV $007*D008+C001*AB001*D011-C050*AB001*D004) SPDSTV IF(ITYPE-6)3261,3262,3261 SPDSTV C ****************************************************************** SPDSTV C * PP * SPDSTV C ****************************************************************** SPDSTV 3261 CONTINUE SPDSTV EP10=OF1*(+C002*AB002*D001-C001*AB001*D002) SPDSTV EP30=OF1*(+C002*AB004*D001-C001*AB001*D005) SPDSTV EP60=OF1*(+C002*AB007*D001-C001*AB001*D008) SPDSTV EP11=OF4*(-C007*AB003*D001+C008*AB001*D001+C006*AB002*D002-C002*AB SPDSTV $002*D002+C001*AB001*D003-C050*AB001*D004) SPDSTV EP31=OF4*(-C007*AB005*D001+C006*AB002*D005-C002*AB004*D002+C001*AB SPDSTV $001*D006) SPDSTV EP61=OF4*(-C007*AB008*D001+C006*AB002*D008-C002*AB007*D002+C001*AB SPDSTV $001*D009) SPDSTV EP13=OF4*(-C007*AB005*D001+C006*AB004*D002-C002*AB002*D005+C001*AB SPDSTV $001*D006) SPDSTV EP33=OF4*(-C007*AB006*D001+C008*AB001*D001+C006*AB004*D005-C002*AB SPDSTV $004*D005+C001*AB001*D007-C050*AB001*D004) SPDSTV EP63=OF4*(-C007*AB009*D001+C006*AB004*D008-C002*AB007*D005+C001*AB SPDSTV $001*D010) SPDSTV EP16=OF4*(-C007*AB008*D001+C006*AB007*D002-C002*AB002*D008+C001*AB SPDSTV $001*D009) SPDSTV EP36=OF4*(-C007*AB009*D001+C006*AB007*D005-C002*AB004*D008+C001*AB SPDSTV $001*D010) SPDSTV EP66=OF4*(-C007*AB010*D001+C008*AB001*D001+C006*AB007*D008-C002*AB SPDSTV $007*D008+C001*AB001*D011-C050*AB001*D004) SPDSTV 3262 CONTINUE DO 2137 I=1,100 SPDSTV 2137 EPN(I)=EPN(I)+EEP(I) SPDSTV 105 CONTINUE SPDSTV C ****************************************************************** SPDSTV C END OF LOOP OVER GAUSSIANS SPDSTV C STORE IN ARRAYS SPDSTV C ****************************************************************** SPDSTV INTC=0 SPDSTV DO 500 J=1,10 SPDSTV R3B=RENORM(J) SPDSTV DO 500 I=1,10 SPDSTV R3A=R3B*RENORM(I) SPDSTV INTC=INTC+1 SPDSTV 500 EPN(INTC) = ( EPN(INTC) )*R3A*SYMFAC CALL REDUC1(EPN,LAMAX,LBMAX,I6TO5) CALL MATFIL(H,EPN,AOS,SHELLT,INEW,JNEW,LAMAX,LBMAX,LA,LB) 1000 CONTINUE SPDSTV C C REFORMAT COMMON /B/ AND THE H ARRAY IF THIS BASIS CONTAINS C P ONLY SHELLS C IF (IPO(4) .EQ. 0) GOTO 1285 WRITE(IOUT,*) 'DEBUG OF UNSTAR' CALL LINOUT (H,NBASIS,0,0) 1285 CONTINUE C CALL UNSTAR (NBASIS,SHELLT,SHELLC,AOS,NSHELL,H,NOSTAR) C IF(IPO(4).EQ.0) GOTO 1500 WRITE(IOUT,2010) CALL LINOUT(H,NBASIS,0,0) 1500 CONTINUE RETURN END SPDSTV SUBROUTINE INV(A,N,IS,IAD1,IAD2,D,MDM) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C ****************************************************************** C INVERSION OF SQUARE MATRIX A BY MEANS OF THE GAUSS-JORDAN C ALGORITHM C C APRIL 72/RS9B C ****************************************************************** DIMENSION A(MDM,MDM),IS(2,MDM),IAD1(MDM),IAD2(MDM),D(MDM) C COMMON/IO/IN,IOUT,IPUNCH C DATA ZERO/0.0D0/, ONE/1.0D0/, SMALL/1.0D-20/ C 2000 FORMAT(' WARNING FROM INV: MATRIX IS SINGULAR') C ****************************************************************** DO 1 L=1,N IS(1,L)=0 1 IS(2,L)=0 DO 9 IMA=1,N B= ZERO DO 2 L=1,N DO 2 M=1,N IF(IS(1,L).EQ.1.OR.IS(2,M).EQ.1) GOTO 2 E=DABS(A(L,M)) IF(E.LT.B) GOTO 8 I=L K=M 8 B=DMAX1(B,E) 2 CONTINUE IS(1,I)=1 IS(2,K)=1 IAD1(K)=I IAD2(I)=K B=A(I,K) C.....PIVOT IF(DABS(B).LT. SMALL) GOTO 20 A(I,K)=ONE/B DO 6 L=1,N IF(L.EQ.K) GOTO 6 C.....KELLERZEILE A(I,L)=-A(I,L)/B 6 CONTINUE DO 5 L=1,N DO 5 M=1,N IF(L.EQ.I.OR.M.EQ.K) GOTO 5 C.....RECHTECK-REGEL A(L,M)=A(L,M)+A(L,K)*A(I,M) 5 CONTINUE DO 11 L=1,N IF(L.EQ.I) GOTO 11 C.....PIVOT-SPALTE A(L,K)=A(L,K)/B 11 CONTINUE 9 CONTINUE C.....PERMUTATION DER ZEILEN, UM DIE NATUERLICHE ORDNUNG WIEDER HERZUSTE DO 15 L=1,N DO 13 J=1,N K=IAD1(J) 13 D(J)=A(K,L) DO 14 J=1,N 14 A(J,L)=D(J) 15 CONTINUE C.....PERMUTATION DER SPALTEN DO 16 L=1,N DO 17 J=1,N K=IAD2(J) 17 D(J)=A(L,K) DO 18 J=1,N 18 A(L,J)=D(J) 16 CONTINUE RETURN C C ERROR EXIT: MATRIX IS SINGULAR 20 WRITE(IOUT,2000) STOP 'INV IN POLAR' END SUBROUTINE LINOUT(X,N,KEY,IZERO) C C GENERAL LINEAR MATRIX OUTPUT ROUTINE C C KEY=0 MATRIX SYMMETRIC C KEY=1 MATRIX SQUARE ASYMMETRIC C C IZERO=0 ZERO MATRIX ELEMENTS LESS THAN CUTOFF C IZERO=1 DO NOT ZERO MATRIX ELEMENTS C C CUTOFF=1.0E-06 C IMPLICIT REAL*8 (A-H,O-Z) COMMON/IO/IN,IOUT C DIMENSION S(9),X(1) C DATA CUTOFF/1.0E-06/ DATA ZERO/0.0E0/ C IA(I)=(I*(I-1))/2 C C ILOWER=1 100 IUPPER=MIN0(ILOWER+8,N) IRANGE=MIN0(IUPPER-ILOWER+1,9) WRITE (IOUT,9000) (J,J=ILOWER,IUPPER) WRITE (IOUT,9010) DO 160 I=1,N K=1 DO 150 J=ILOWER,IUPPER IF(KEY)110,120,110 110 IJ=N*(J-1)+I GO TO 140 120 IJ=IA(I)+J IF(I-J)130,140,140 130 IJ=IA(J)+I 140 S(K)=X(IJ) IF(IZERO.EQ.0.AND.ABS(S(K)).LE.CUTOFF) S(K)=ZERO 150 K=K+1 160 WRITE (IOUT,9020) I,(S(J),J=1,IRANGE) WRITE (IOUT,9010) ILOWER=ILOWER+9 IF(N-IUPPER)170,170,100 170 RETURN 9000 FORMAT(12X,8(I3,11X),I3) 9010 FORMAT(/) 9020 FORMAT(1X,I3,2X,9E14.6) END SUBROUTINE MATFIL(A,AA,AOS,SHELLT,INEW,JNEW,LAMAX,LBMAX,LA,LB) C C GAUSSIAN 77/UCI C IMPLICIT REAL*8 (A-H,O-Z) INTEGER AOS(1), SHELLT(1) C DIMENSION A(1),AA(1) C LIND(I)=(I*(I-1))/2 C ISTART=AOS(INEW) MATFIL JSTART=AOS(JNEW) MATFIL IAL = 0 IAU = 5 IBL = 0 IBU = 5 IMA = 0 IMB = 0 IF(SHELLT(INEW) .EQ. 2) IMA = 1 IF(SHELLT(JNEW) .EQ. 2) IMB = 1 C MATFIL 120 INTC=0 MATFIL DO 170 J=1,LBMAX MATFIL DO 170 I=1,LAMAX MATFIL INTC=INTC+1 MATFIL IF( LA.GT.1 .AND. I.GT.IAL .AND. I.LT.IAU )GO TO 170 MATFIL IF( LA.EQ.1 .AND. I.EQ.IMA )GO TO 170 MATFIL IF( LB.GT.1 .AND. J.GT.IBL .AND. J.LT.IBU )GO TO 170 MATFIL IF( LB.EQ.1 .AND. J.EQ.IMB )GO TO 170 MATFIL IND=ISTART+I-1 MATFIL JND=JSTART+J-1 MATFIL IF(IND-JND)130,140,150 MATFIL 130 IJ=LIND(JND)+IND MATFIL GO TO 160 MATFIL 140 IJ=LIND(IND+1) MATFIL GO TO 160 MATFIL 150 IJ=LIND(IND)+JND MATFIL 160 A(IJ)=AA(INTC) MATFIL 170 CONTINUE MATFIL RETURN MATFIL END MATFIL SUBROUTINE MULTAY(A,Y,X,N,MAXDIM) C C MATRIX MULTIPLICATION ROUTINE C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION A(MAXDIM,MAXDIM),Y(MAXDIM),X(MAXDIM) C DATA ZERO/0.0/ C DO 200 IROW=1,N SUM = ZERO DO 100 JCOL=1,N SUM = SUM + A(IROW,JCOL) * Y(JCOL) 100 CONTINUE X(IROW) = SUM 200 CONTINUE RETURN END C SUBROUTINE OUTPUT C C C L.E. CHIRLIAN C APRIL 1985 C C A SUBROUTINE TO OUTPUT THE CHARGES AND OTHER PERTINANT C INFORMATION FROM THE CHELP PROGRAM C C Slightly Modified for CHELPG operations by Curt Breneman C Yale University Department of Chemistry, 2/88 C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NPOINTS = 50000) INTEGER*4 SHELLA,SHELLN,SHELLT,SHELLC,AOS,AON,SHLADF CHARACTER*40 CHKFIL C COMMON /IO/ IN,IOUT COMMON /IPO/ IPO(5) C+++ COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NE,NBASIS, $ IAN(401),ATMCHG(400),C(3,400) C+++ C COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NE,NBASIS,IAN(101), C 1 ATMCHG(100),C(3,100) COMMON /OUT/ NTITLE(20,3),CHKFIL,Q(400),ND,RMS,PD, 1 NLIN,NEND(3) COMMON /POINTS/ P(3,NPOINTS), NP DATA DEB/0.393427328/ c c c write(6,*) 'Debug --', nlin,nend(1),nend(2),nend(3),nwords C C CALCULATE THE DIPOLE MOMENT FROM THE FITTED CHARGES C C DIPX=0. DIPY=0. DIPZ=0. DO 99 I=1,NATOMS DIPX=DIPX+(Q(I)*C(1,I)) DIPY=DIPY+(Q(I)*C(2,I)) DIPZ=DIPZ+(Q(I)*C(3,I)) 99 CONTINUE DIPX=DIPX/DEB DIPY=DIPY/DEB DIPZ=DIPZ/DEB C C CALCULATE TOTAL DIPOLE MOMENT C DIPTOT = DSQRT(DIPX**2+DIPY**2+DIPZ**2) C C CREATE OUTPUT C WRITE (IOUT,100) 100 FORMAT (/,/,17X,'CHARGES FROM ELECTROSTATIC POTENTIAL GRID') WRITE (IOUT,110) WRITE (IOUT,111) 110 FORMAT (/,36X,'CHELPGrid',/) 111 FORMAT (/,15X,'Grid Modification.') DO 24 I=1,NLIN WRITE (6,1200)(NTITLE(J,I),J=1,NEND(I)) 1200 FORMAT(2X,19A4) 24 CONTINUE C C WRITE CHECKPOINT FILE NAME C WRITE (IOUT,150)CHKFIL 150 FORMAT(/2X,'CHECKPOINT FILE: ',A40) C C PRINT DATE C C*** Take out for Trace-7 C CALL FOR$JDATE(IMONTH,IDATE,IYEAR) C*** WRITE (IOUT,170)IMONTH,IDATE,IYEAR 170 FORMAT (/2X,I2,'-'I2,'-'I2) C C************************************************************************** C WRITE GEOMETRY C************************************************************************** C WRITE (IOUT,180) 180 FORMAT (/2X,36X,'MOLECULAR GEOMETRY') WRITE (IOUT,190) 190 FORMAT (/,/,17X,'ATOMIC NUMBER',8X,'X',12X,'Y',12X,'Z') DO 30 I=1,NATOMS WRITE (IOUT,200)IAN(I),C(1,I),C(2,I),C(3,I) 200 FORMAT (/,23X,I2,8X,F10.7,3X,F10.7,3X,F10.7) 30 CONTINUE WRITE (IOUT,210)ICHARG 210 FORMAT (/2X,'THE TOTAL CHARGE IS CONSTRAINED TO: ',I3) WRITE (IOUT,240) 240 FORMAT (/,36X,'NET CHARGES') WRITE (IOUT,250) 250 FORMAT (/,28X,'ATOMIC NUMBER',5X,'CHARGE') WRITE (IOUT,260)(IAN(I),Q(I),I=1,NATOMS) WRITE (6,101) DIPTOT 101 FORMAT (/,2X,'THE DIPOLE MOMENT OF THESE CHARGES IS: ',F8.5) 260 FORMAT (/,34X,I2,10X,F8.4) WRITE (IOUT,270)NP 270 FORMAT(/,2X,'FIT TO ELECTROSTATIC POTENTIAL AT ',I6,' POINTS') WRITE (IOUT,280)RMS 280 FORMAT (/,2X,'ROOT MEAN SQUARE DEVIATION IS ',F6.4,' KCAL/MOLE') RETURN END SUBROUTINE READIN C C WRITTEN BY M.M. FRANCL FOR THE C PRINCETON CHEMISTRY DEPARTMENT VAX 11/780 UNDER VMS 3.4. C MODIFIED BY L.E. CHIRLIAN UNDER VMS 3.7. C C MODIFIED FOR GAUSSIAN 86 BY CURT BRENEMAN (VMS 4.5) C Modified for G88/90 by Curt Breneman, 2/89 C YALE UNIVERSITY DEPARTMENT OF CHEMISTRY. C C THIS VERSION IS COMPATIBLE WITH GAUSSIAN 82 FROM CARNEGIE- C MELLON UNIVERSITY AND IS DESIGNED FOR THE INPUT OF MO AND C BASIS INFORMATION FROM CHECKPOINT FILES. THIS VERSION IS C TO BE USED FOR THE DETERMINATION OF ATOMIC CHARGES FROM C ELECTROSTATIC POTENTIALS DETERMINED BY FIRST ORDER HARTREE C FOCK PERURBATION THEORY. C C OLD LIMITATIONS: NO MORE THAN 256 BASIS FUNCTIONS C NO MORE THAN 80 SHELLS C C NEW LIMITATIONS: NO MORE THAN 1280 BASIS FUNCTIONS C NO MORE THAN 400 SHELLS C IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4 SHELLA,SHELLN,SHELLT,SHELLC,AOS,AON,SHLADF,FILNUM CHARACTER*40 CHKFIL PARAMETER (NUMPTS = 20) DIMENSION LINE(20) C COMMON /IO/ IN,IOUT c+++ c Change for G86 : New Commons /MOL/ and /B/ c COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NE,NBASIS, $ IAN(401),ATMCHG(400),C(3,400) C C=== Gaussian88 Modification. New Common /b/ size. Common/B/EXX(6000),C1(6000),C2(6000),C3(2000),CF(2000), $SHLADF(4000),X(2000),Y(2000), $Z(2000),JAN(2000),ShellA(2000),ShellN(2000),ShellT(2000), $ShellC(2000),AOS(2000),AON(2000),NShell,MaxTyp C=== Old G86 Common /b/ c COMMON/B/EXX(1200),C1(1200),C2(1200),C3(400),CF(400),SHLADF(800), c $ X(400),Y(400),Z(400),JAN(400),SHELLA(400),SHELLN(400), c $ SHELLT(400),SHELLC(400),AOS(400),AON(400),NSHELL,MAXTYP c c+++ C%%% c Original CHELP common /B/ c c COMMON /B/ EXX(240),C1(240),C2(240),C3(80),CF(80),SHLADF(160), c $ X(80),Y(80),Z(80), c $ JAN(80),SHELLA(80),SHELLN(80),SHELLT(80),SHELLC(80), c $ AOS(80),AON(80),NSHELL,MAXTYP C COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NE,NBASIS,IAN(101), C $ ATMCHG(100),C(3,100) C%%% COMMON /IPO/ IPO(5) COMMON /OUT/ NTITLE(20,3),CHKFIL,Q(400),ND,RMS,PD, 1 NLIN,NEND(3) COMMON /CHARGE/ COEF_ALPHA(100000),COEF_BETA(100000),IUHF COMMON /SPHERE/ VDWR(400), NPTS C VECT(3,NUMPTS) C DATA MAXBAS /1280/ DATA MAXPTS/ 50000/ DATA NPO/ 5/ DATA IUNIT/ 3/, IREAD/ 2/, IFIND/11/, IBLNK/4H / C DATA VECT/0.00000000, 0.00000000, 1.00000000, C $ 0.23807485, -0.08801514, -0.96725059, C $ -0.46055608, 0.17026540, 0.87114740, C $ 0.65287142, -0.24136347, -0.71798508, C $ -0.80242446, 0.29665251, 0.51779559, C $ 0.00000000, 0.00000000, 1.00000000, C $ -0.20401674, -0.15100817, -0.96725059, C $ 0.39467063, 0.29212549, 0.87114740, C $ -0.55947405, -0.41410893, -0.71798508, C $ 0.68763258, 0.50896872, 0.51779559, C $ -0.94978689, -0.28553486, -0.12796369, C $ 0.88757697, 0.26683266, 0.37550960, C $ -0.76723180, -0.23065324, -0.59846007, C $ 0.59663385, 0.17936630, 0.78221211, C $ -0.38695708, -0.11633108, -0.91473018, C $ 0.28119199, 0.95108168, -0.12796369, C $ -0.26277424, -0.88878695, 0.37550960, C $ 0.22714509, 0.76827772, -0.59846007, C $ -0.17663821, -0.59744720, 0.78221211, C $ 0.11456173, 0.38748460, -0.91473018/ C C Old "spherical" unit vectors (14 of them) C c 0.5773502691896258,0.5773502691896258, c $ 0.5773502691896258, c 1 -0.5773502691896258,-0.5773502691896258,0.5773502691896258, c 2 0.5773502691896258,-0.5773502691896258,-0.5773502691896258, c 3 -0.5773502691896258,0.5773502691896258,-0.5773502691896258, c 4 0.0000000000000000E+00,-1.000000000000000, c $ 0.0000000000000000E+00, c 5 0.0000000000000000E+00,0.0000000000000000E+00, c $ -1.000000000000000, c 6 -1.000000000000000,0.0000000000000000E+00, c $ 0.0000000000000000E+00, c 7 -0.5773502691896258,-0.5773502691896258,-0.5773502691896258, c 8 1.000000000000000,0.0000000000000000E+00, c $ 0.0000000000000000E+00, c 9 0.0000000000000000E+00,1.000000000000000, c $ 0.0000000000000000E+00, c $ 0.5773502691896258,0.5773502691896258,-0.5773502691896258, c $ 0.0000000000000000E+00,0.0000000000000000E+00, c $ 1.000000000000000, c $ -0.5773502691896258,0.5773502691896258,0.5773502691896258, c $ 0.5773502691896258,-0.5773502691896258,0.5773502691896258, IN = 5 IOUT = 6 C C CHECKPOINT FILE NAME C READ(IN,1000) CHKFIL 1000 FORMAT(A40) C C INPUT INFORMATION C 1010 FORMAT(20A4) DO 1921 I=1,3 NLIN= I READ (5,1010) LINE IF (LINE(1) .EQ. IBLNK) THEN NLIN=NLIN-1 GOTO 192 END IF DO 1922 J=1,20 L=20-J IF (LINE(L) .NE. IBLNK) THEN NEND(I) = L DO 1923 K=1,NEND(I) NTITLE(K,I) = LINE(K) 1923 CONTINUE GOTO 1921 END IF 1922 CONTINUE 1921 CONTINUE NLIN=NLIN+1 192 CONTINUE C C C READ IN PRINT OPTIONS C 3000 READ(IN,*) (IPO(I),I=1,NPO) C C READ IN # OF D FUNCTIONS C NOTE: IF THE BASIS SET USES 5 D FUNCTION, ND MUST BE C SET EQUAL TO 1 TO ACCOMODATE THE INTEGRAL PACKAGE. C IF THE BASIS SET USES 6 D FUNCTIONS, ND IS SET EQUAL TO C 0. C READ(IN,*) ND IF (ND .NE. 5 .AND. ND .NE. 6) THEN STOP '# OF D FUNCTIONS MUST BE 5 OR 6' END IF IF (ND .EQ. 5) THEN ND = 1 GOTO 15 END IF ND = 0 15 CONTINUE C C C INITIATE FILEIO C C*** Different Fopen statement for Trace-7 CALL FOPEN (IUNIT,5,CHKFIL,IALLOC,junk) c CALL FOPEN (IUNIT,'old',CHKFIL(1:linend(chkfil))//char(0)) C IWWRIT = IPO(1) C*********************************************************************** C C READ IN COMMON /MOL/ C c NWORDS = 1804 c IFILENO = 30997 c CALL FILEIO (IREAD,IFILENO,NWORDS,NATOMS,IALLOC) IRwMol=997 MaxAtm=400 LenMol = 4*MaxAtm + InToWP(8+MaxAtm) Call FileIO(2,-FilNum(IRwMol,IUnit),LenMol,NAtoms,0) C C*********************************************************************** C C READ IN SPHERE DATA (VAN DER WAALS RADII, # OF POINTS TO C FIT C C READ (IN,*)(VDWR(I),I=1,NATOMS) READ (IN,*)NPTS IF (NPTS .GT. MAXPTS) THEN STOP 'MAXIMUM NUMBER OF POINTS MUST BE LESS THAN 50000' END IF C C SET MAXIMUM VAN DER WAALS RADII TO 4 C VMAX=4. DO 20 I=1,NATOMS IF (VDWR(I) .GT. VMAX) THEN WRITE (IOUT, 2500) I 2500 FORMAT (3X, 'THE VAN DER WAALS RADII OF ATOM', I3, 1 'IS OUT OF RANGE') STOP END IF 20 CONTINUE C READ (IN,*) VFACT DO 21 I=1,NATOMS VDWR(I)=VDWR(I)*VFACT 21 CONTINUE C*********************************************************************** C C READ IN BASIS SET INFORMATION (COMMON /B/) C IFILENO = 30506 CALL FILEIO (IFIND,IFILENO,NWORDS,EXX,0) CALL FILEIO (IREAD,-IFILENO,NWORDS,EXX,0) C*********************************************************************** IF(IWWRIT .NE. 1) GOTO 170 WRITE(IOUT,8000)(C(1,I),C(2,I),C(3,I),I=1,NATOMS) 8000 FORMAT(/1X,'COORDINATES'/(1X,3F12.6)) WRITE(IOUT,8020) NATOMS,ICHARG,MULTIP,NAE,NBE,NE,NBASIS $ ,NSHELL,MAXTYP 8020 FORMAT(/1X,'NATOMS = ',I3 $/1X,'ICHARG = ',I3 $/1X,'MULTIP = ',I3 $/1X,'NAE = ',I3 $/1X,'NBE = ',I3 $/1X,'NE = ',I3 $/1X,'NBASIS = ',I3 $/1X,'NSHELL = ',I3 $/1X,'MAXTYP = ',I3) WRITE(IOUT,8030) (IAN(I),I=1,NATOMS) 8030 FORMAT(/1X,'IAN'/(1X,20I3)) WRITE(IOUT,8050) (JAN(I),SHELLT(I),SHELLA(I),I=1,NSHELL) 8050 FORMAT(/1X,'CENTER TYPE SHELLA'/(1X,3I7)) WRITE(IOUT,8055) SHELLA(NSHELL+1) 8055 FORMAT(1X,14X,I7) WRITE(IOUT,8060) (EXX(I),C1(I),C2(I),I=1,NSHELL) 8060 FORMAT(/1X,12X,'EXPON',8X,'EXPCOF(S)',8X,'EXPCOF(P)', $/(1X,3E17.9)) WRITE(IOUT,8070) (C3(I),CF(I),I=1,NSHELL) 8070 FORMAT(/1X,12X,'EXPCOF(D)',8X,'EXPCOF(F)', $/(1X,2E17.9)) WRITE (IOUT,3575) 3575 FORMAT(/,15X,'ATOM #',5X,'V.D.W. RADII (MULTIPLIED BY FACTOR)') DO 3500 I=1,NATOMS WRITE (IOUT,4000)I,VDWR(I) 4000 FORMAT (15X,I5,20X,F5.2) 3500 CONTINUE c WRITE (IOUT,4500)NPTS 4500 FORMAT(/2X,'NUMBER OF POINTS TO FIT',I6) 170 CONTINUE C*********************************************************************** C C READ IN ALPHA MO COEFFICIENTS C IFILENO = 30524 CALL FILEIO (IFIND,-IFILENO,NWORDS,COEF_ALPHA,0) CALL FILEIO (IREAD,-IFILENO,NWORDS,COEF_ALPHA,0) C*********************************************************************** C C READ IN THE BETA MO COEFFICIENTS C IFILENO = 30526 CALL FILEIO (IFIND,IFILENO,NWORDS,COEF_BETA,0) IF (NWORDS.EQ.0) THEN IUHF = 0 GOTO 300 END IF IUHF = 1 CALL FILEIO (IREAD,-IFILENO,NWORDS,COEF_BETA,0) C*********************************************************************** 300 CONTINUE C*********************************************************************** RETURN END SUBROUTINE REDUC1(X,LAMAX,LBMAX,I6TO5) C C MODIFIED FOR POLARIZATION POTENTIAL CALCULATIONS C M.M. FRANCL FEBRUARY 1984 C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION X(100),S(100),IND5(9),IND6(10) C DATA PT5/0.5/ DATA R3OV2/0.8660254040/ DATA IND5/1, $ 4,7,2, $ 3,6,9,5,8/ DATA IND6/1, $ 4,7,2, $ 6,10,3,9,5,8/ C C ****************************************************************** C ROUTINE REORDERS FROM ARRANGEMENT: S,Z,ZZ,X,XZ,XX,Y,YZ,XY,YY C TO S,X,Y,Z,XX,YY,ZZ,XY,XZ,YZ C OR FROM S,Z,X,Y TO S,X,Y,Z C THIS ENSURES LABELING COMPATIBILITY BETWEEN THE SP AND D PACKAGES C AT THE SAME TIME THE INTEGRALS ARE MOVED TO THE FIRST 1,4,10,16, C 40 OR 100 LOCATIONS, DEPENDING ON THE SHELL QUANTUM NUMBERS C ****************************************************************** NWORD=LAMAX*LBMAX IF(NWORD-1)5,180,5 5 INTC=0 IF(I6TO5 .EQ. 1) GOTO 40 10 DO 20 I=1,LBMAX ISB=10*(IND6(I)-1) DO 20 J=1,LAMAX ISA=ISB+IND6(J) INTC=INTC+1 20 S(INTC)=X(ISA) GO TO 160 C ****************************************************************** C ROUTINE TO REDUCE SIX D FUNCTIONS TO FIVE C ALSO REORDERS FROM : S,Z,ZZ,X,XZ,XX,Y,YZ,YX,YY C TO S,X,Y,Z,ZZ,XX-YY,XY,XZ,YZ C OR FROM S,Z,X,Y TO S,X,Y,Z C FOR COMPATIBILITY WITH SP PACKAGE C ****************************************************************** 40 DO 150 I=1,LBMAX ISB=10*(IND5(I)-1) C IFB=0 FOR S,X,Y,Z,XY,XZ,YZ, IFB=1 FOR ZZ-RR, IFB=2 FOR XX-YY IFB = 0 IF(I .EQ. 5) IFB = 1 IF(I .EQ. 6) IFB = 2 80 DO 150 J=1,LAMAX ISA=ISB+IND5(J) IFA = 0 IF(J .EQ. 5) IFA = 1 IF(J .EQ. 6) IFA = 2 120 IHOP = 3*IFB + IFA + 1 GOTO(130,122,123,124,125,126,127,128,129),IHOP C C ****************************************************************** C * (F,O,ZA2) * C ****************************************************************** 122 XX=ZZ1(X,ISA,3,7) GO TO 140 C C ****************************************************************** C * (F,O,XA2-YA2) * C ****************************************************************** 123 XX=XY1(X,ISA,4) GO TO 140 C C ****************************************************************** C * (ZB2,O,F) * C ****************************************************************** 124 XX=ZZ1(X,ISA,30,70) GO TO 140 C C ****************************************************************** C * (ZB2,O,ZA2) * C ****************************************************************** 125 XX=ZZ1(X,ISA,30,70)-PT5*(ZZ1(X,ISA+3,30,70)+ZZ1(X,ISA+7,30,70)) GO TO 140 C C ****************************************************************** C * (ZB2,O,XA2-YA2) * C ****************************************************************** 126 XX=R3OV2*(ZZ1(X,ISA,30,70)-ZZ1(X,ISA+4,30,70)) GO TO 140 C C ****************************************************************** C * (XB2-YB2,O,F) * C ****************************************************************** 127 XX=XY1(X,ISA,40) GO TO 140 C C ****************************************************************** C * (XB2-YB2,O,ZA2) * C ****************************************************************** 128 XX=R3OV2*(ZZ1(X,ISA,3,7)-ZZ1(X,ISA+40,3,7)) GO TO 140 C C ****************************************************************** C * (XB2-YB2,O,XA2-YA2) * C ****************************************************************** 129 XX=R3OV2*(XY1(X,ISA,4)-XY1(X,ISA+40,4)) GO TO 140 C C ****************************************************************** C * (F,O,F) * C ****************************************************************** 130 XX=X(ISA) 140 INTC=INTC+1 150 S(INTC)=XX 160 DO 170 I=1,NWORD 170 X(I)=S(I) 180 RETURN END SUBROUTINE STAR(NBASIS,SHELLT,SHELLC,AOS,NSHELL,NOSTAR) C C ROUTINE TO MODIFY COMMON /B/ TO THE EXPECTED FORMAT FOR INTGRL C FOR BASIS SETS HAVING P ONLY SHELLS, SUCH AS THE 6-31G** BASIS C IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4 SHELLC,SHELLT,AOS C DIMENSION SHELLC(2000),SHELLT(2000),AOS(2000) C C LOOP OVER SHELLS C DO 100 I=1,NSHELL IF (SHELLT(I).EQ.1 .AND. SHELLC(I).EQ.1) THEN NBASIS = NBASIS + 1 NOSTAR = 1 DO 200 J=I,NSHELL AOS(J) = AOS(J) + 1 200 CONTINUE END IF 100 CONTINUE RETURN END SUBROUTINE UEP C C ROUTINE TO CALCULATE THE ELECTROSTATIC POTENTIAL FROM FIRST ORDER C PERTURBATION THEORY C C M.M. FRANCL JULY 1985 C MODIFIED VERSION OF A MEPHISTO ROUTINE C RESTRICTED TO UNRESTRICTED HARTREE-FOCK WAVEFUNCTIONS C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NPOINTS = 50000) INTEGER*4 SHELLA,SHELLN,SHELLT,AOS,SHELLC,AON,HANDLE C COMMON /IO/ IN,IOUT COMMON /IPO/ IPO(5) C+++ COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS, $ IAN(401),ATMCHG(400),C(3,400) C C=== Gaussian88 Modification for enlarged common /b/. Common/BB/EXX(6000),C1(6000),C2(6000),C3(6000),X(2000),Y(2000), $Z(2000),JAN(2000),ShellA(2000),ShellN(2000),ShellT(2000), $ShellC(2000),AOS(2000),AON(2000),NShell,MaxTyp C C=== Old G86 Common /b/ c COMMON/B/EXX(1200),C1(1200),C2(1200),C3(1200), c $ X(400),Y(400),Z(400),JAN(400),SHELLA(400),SHELLN(400), c $ SHELLT(400),SHELLC(400),AOS(400),AON(400),NSHELL,MAXTYP C+++ C COMMON /B/ EXX(240),C1(240),C2(240),C3(240),X(80),Y(80),Z(80), C $ JAN(80),SHELLA(80),SHELLN(80),SHELLT(80),SHELLC(80) C $ ,AOS(80),AON(80),NSHELL,MAXTYP C COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101), C $ ATMCHG(100),C(3,100) COMMON /POINTS/ P(3,NPOINTS),MAXPNTS COMMON /ELP/ ELECP(NPOINTS) COMMON /CHARGE/ COEF_ALPHA(100000),COEF_BETA(100000),IUHF COMMON /OUT/ NTITLE(20,3),CHKFIL,Q(400),I6TO5,RMS,PERCENT, 1 NLIN,NEND(3) C DIMENSION HPERT(100000),INDEX(1280) C DATA IPTCHG/1.0/ DATA ZERO/0.0/, TWO/2.0/, VNUCMAX/30.0/ C C INTIALIZE TIMING C HANDLE = 0 C C SET UP THE INDEXING TABLE FOR HPERT C DO 100 I=1,NBASIS INDEX(I) = (I-1)*I/2 100 CONTINUE C C BEGIN LOOP TO CALCULATE ELECTROSTATIC POTENTIAL C DO 200 NPNT=1,MAXPNTS X1 = P(1,NPNT) X2 = P(2,NPNT) X3 = P(3,NPNT) C C CALCULATE THE ONE-ELECTRON INTEGRALS C IF (IPO(5).EQ.1) THEN WRITE(IOUT,3010) 3010 FORMAT(1X,'TIME FOR INTEGRALS') C*** C ISTAT = LIB$INIT_TIMER(HANDLE) C*** END IF C CALL INTGRL (HPERT,X1,X2,X3,IPTCHG,I6TO5) C C*** C IF (IPO(5).EQ.1) ISTAT = LIB$SHOW_TIMER(HANDLE) C*** C IF (IPO(4).EQ.1) CALL LINOUT (HPERT,NBASIS,0,0) C IF (IPO(5).EQ.1) THEN WRITE(IOUT,3000) 3000 FORMAT(1X,'TIME FOR TRANSFORM') C*** C ISTAT = LIB$INIT_TIMER(HANDLE) C*** END IF C C FORM THE HPERT MATRIX ELEMENTS C C ALPHA CODE C E = ZERO ICOEFI = -NBASIS C C SUM OVER OCCUPIED ALPHA MOS C DO 220 II=1,NAE ICOEFI = ICOEFI + NBASIS C C CALCULATE ELECTROSTATIC POTENTIAL C DO 221 IP=1,NBASIS CPI = COEF_ALPHA(ICOEFI+IP) IPDEX = INDEX(IP) C DO 222 IQ=1,IP E = E + CPI * COEF_ALPHA(ICOEFI+IQ) * HPERT(IPDEX+IQ) 222 CONTINUE DO 223 IQ=IP+1,NBASIS E = E + CPI * COEF_ALPHA(ICOEFI+IQ) * HPERT(IP+INDEX(IQ)) 223 CONTINUE C 221 CONTINUE 220 CONTINUE C C BETA CODE C ICOEFI = -NBASIS C C SUM OVER OCCUPIED BETA MOS C DO 420 II=1,NBE ICOEFI = ICOEFI + NBASIS C C CALCULATE ELECTROSTATIC POTENTIAL C DO 421 IP=1,NBASIS CPI = COEF_BETA(ICOEFI+IP) IPDEX = INDEX(IP) C DO 422 IQ=1,IP E = E + CPI * COEF_BETA(ICOEFI+IQ) * HPERT(IPDEX+IQ) 422 CONTINUE DO 423 IQ=IP+1,NBASIS E = E + CPI * COEF_BETA(ICOEFI+IQ) * HPERT(IP+INDEX(IQ)) E = E + CPI * COEF_* HPERT(IP+INDEX(IQ)) 423 CONTINUE C 421 CONTINUE 420 CONTINUE C C*** C IF (IPO(5) .EQ. 1) ISTAT = LIB$SHOW_TIMER(HANDLE) C*** C C CALCULATE NUCLEAR PART OF ELECTROSTATIC POTENTIAL C VNUC = ZERO DO 300 IATOM=1,NATOMS DEL1 = C(1,IATOM) - X1 DEL2 = C(2,IATOM) - X2 DEL3 = C(3,IATOM) - X3 RA = DSQRT(DEL1*DEL1 + DEL2*DEL2 + DEL3*DEL3) IF (RA.EQ.ZERO) THEN VNUC=VNUCMAX GOTO 310 END IF VNUC = VNUC + IAN(IATOM) / RA 300 CONTINUE 310 CONTINUE C ELECP(NPNT) = (E + VNUC * IPTCHG) IF (IPO(5) .EQ. 1) WRITE(IOUT,*) 'E(',NPNT,') = ',E 200 CONTINUE RETURN END SUBROUTINE UNSTAR(NBASIS,SHELLT,SHELLC,AOS,NSHELL,H,NOSTAR) C C ROUTINE TO REFORMAT COMMON/B/ AND THE H ARRAY FOR BASIS C SETS HAVING P ONLY SHELLS C IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4 SHELLC,SHELLT,AOS C DIMENSION SHELLC(2000),SHELLT(2000),AOS(2000),H(1) C IF (NOSTAR.EQ.0) RETURN C C LOOP OVER SHELLS C DO 100 I=1,NSHELL IF (SHELLT(I).EQ.1 .AND. SHELLC(I).EQ.1) THEN C C REMOVE EXTRA ROWS AND COLUMNS C C LOOP OVER ROWS C IBASIS = AOS(I) ITEM = (IBASIS-1) * IBASIS / 2 DO 200 J=IBASIS+1,NBASIS C C LOOP OVER COLUMNS C NEWITEM = (J-1) * J /2 DO 250 K=1,IBASIS-1 ITEM = ITEM + 1 NEWITEM = NEWITEM + 1 H(ITEM) = H(NEWITEM) 250 CONTINUE C C SKIP THE VALUE IN THE IBASIS TH COLUMN C NEWITEM = NEWITEM + 1 C DO 260 K=IBASIS+1,J ITEM = ITEM + 1 NEWITEM = NEWITEM + 1 H(ITEM) = H(NEWITEM) 260 CONTINUE 200 CONTINUE NBASIS = NBASIS - 1 C C RESTRUCTURE AOS TO ACCOUNT FOR THE S SHELL REMOVED C DO 300 IAOS = I,NSHELL AOS(IAOS) = AOS(IAOS) - 1 300 CONTINUE C END IF 100 CONTINUE RETURN END FUNCTION XY1(X,I,IY) C DIMENSION X(100) C DATA HALFR3/0.8660254040/ C XY1=HALFR3*(X(I)-X(I+IY)) RETURN END FUNCTION ZZ1(X,I,IX,IY) C DIMENSION X(100) C DATA HALF/0.5/ C ZZ1=X(I)-HALF*(X(I+IX)+X(I+IY)) RETURN END