(****************************************************************************
* PPP-mo calculation program. The original program, written in Algol, came  *
* from the Prof. John Griffiths group at the University of Leeds. This is   *
* the pascal translation of that program. This is the Apple II UCSD Pascal  *
* version.                                                                  *
*                                                                           *
* Copyright (C) <1985> <Giorgio Castellani>                                 *
* This program is free software: you can redistribute it and/or             *
* modify it under the terms of the GNU General Public License               *
* as published by the Free Software Foundation,                             *
* either version 3 of the License, or (at your option) any later version.   *
* This program is distributed in the hope that it will be useful,           *
* but WITHOUT ANY WARRANTY; without even the implied                        *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.          *
* See the GNU General Public License for more details.                      *
* You should have received a copy of the GNU General Public License         *
* along with this program. If not, see <http://www.gnu.org/licenses/>.      *
*****************************************************************************)

(*L PRINTER:*)
PROGRAM PPP;
  FUNCTION SIN(X:REAL):REAL;
  FUNCTION COS(X:REAL):REAL;
  FUNCTION EXP(X:REAL):REAL;
  FUNCTION ATAN(X:REAL):REAL;
  FUNCTION LN(X:REAL):REAL;
  FUNCTION LOG(X:REAL):REAL;
  FUNCTION SQRT(X:REAL):REAL;

USES TRASCEND;
CONST M=30;
TYPE VECTN= ARRAY  [1..M] OF REAL;
     ARRN= ARRAY [1..M] OF ARRAY [1..M] OF REAL;
VAR I,J,N,NF,PM,NIT,CI,NU,NL,NG,STI,STJ:INTEGER;
    C,PRC,P,Q,BETA,GAMMA:ARRN;
    EIG,CX,CY,CZ,AI,E,Z:VECTN;
    CONV:REAL;
    OUT:CHAR;
    FO:FILE OF CHART;

PROCEDURE PLOTT;
VAR I,J,K,L:INTEGER;
    YMIN,ZMIN,YMAX,ZMAX:REAL;
    CYP,CZP:VECTN;

(*$i Pictur.pas*)

BEGIN
  YMAX:=0;
  ZMAX:=0;
  YMIN:=999;
  ZMIN:=999;
  FOR I:=1 TO N DO BEGIN
   IF CY[I] < YMIN THEN
      YMIN:=CY[I];
   IF CZ[I] < ZMIN THEN
      ZMIN:=CZ[I]
  END;
  FOR I:=1 TO N DO BEGIN
   CYP[I]:=CY[I]-YMIN;
   CZP[I]:=CZ[I]-ZMIN
  END;
  FOR I:=1 TO N DO BEGIN
   IF CYP[I] > YMAX THEN
      YMAX:=CYP[I];
   IF CZP[I] > ZMAX THEN
      ZMAX:=CZP[I]
  END;
  IF YMAX > ZMAX THEN
     PICTURE(CYP,CZP,YMAX,ZMAX)
  ELSE
     PICTURE(CZP,CYP,ZMAX,YMAX);
END;

PROCEDURE NEWLINE(NLIN:INTEGER);FORWARD;

SEGMENT PROCEDURE HEADER;
BEGIN
REPEAT
 WRITE(CHR(12));
 WRITE('OUTPUT:C(ONSOLE,P(RINTER,R(EMOUT');
 READLN(OUT);
 CASE OUT OF
   'C':REWRITE(FO,'CONSOLE:');
   'P':REWRITE(FO,'PRINTER:');
   'R':REWRITE(FO,'REMOUT:');
 END;
 UNTIL OUT IN ['C','P','R'];
 WRITELN(FO,'PPP PROGRAM');
 WRITELN(FO,'PASCAL VERSION BY DR.G.CASTELLANI');
 NEWLINE(2);
END;

SEGMENT PROCEDURE REEDDATA;
VAR I,J,K,CONA,CONB: INTEGER;
    A,X,Y,Z1,RD:REAL;

PROCEDURE REEDSM(VAR A:VECTN);
VAR I:INTEGER;
BEGIN
 FOR I:=1 TO N DO
   READ(A[I])
 END;

PROCEDURE SHAPE(VAR R:ARRN);
VAR I,J:INTEGER;
    ALPHA,R1:REAL;
    ANG:VECTN;
BEGIN
 FOR I:=1 TO N DO BEGIN
  CX[I]:=0;
   FOR J:=I TO N DO BEGIN
    R[I,J]:=0;
    R[J,I]:=0
   END
  END;
  CZ[1]:=0;
  CY[1]:=0;
  ANG[1]:=0;
  REPEAT
   READ(I);
   IF I <> 0 THEN BEGIN
      READ(J);
      READ(R1);
      R[I,J]:=R1;
      R[J,I]:=R1;
      READ(ALPHA);
      ALPHA:=ALPHA+ANG[I];
      IF ALPHA >= 360.0 THEN
         ALPHA:=ALPHA-360.0;
      ANG[J]:=ALPHA;
      ALPHA:=ALPHA*0.01745327;
      CZ[J]:=CZ[I]+R1*COS(ALPHA);
      CY[J]:=CY[I]+R1*SIN(ALPHA);
   END;
  UNTIL I = 0;
  FOR I:=1 TO N-1 DO
   FOR J:=I+1 TO N DO
    IF R[I,J] <= 0.001 THEN BEGIN
       R1:=CZ[J]-CZ[I];
       ALPHA:=CY[J]-CY[I];
       R[I,J]:=SQRT(SQR(ALPHA)+SQR(R1));
       R[J,I]:=R[I,J]
    END
END;

PROCEDURE COORDS;
VAR I,J:INTEGER;
    A,B,C,D,MAX : REAL;
BEGIN
  CX[3]:=0;
  A:= GAMMA[1,3];
  WRITELN(A);
  A:= SQR(A);
  B:=GAMMA[2,3];
  C:=CZ[2];
  C:=SQR(C);
  B:=(A-SQR(B)+C)/CZ[2]*0.5;
  CZ[3]:=B;
  B:=A-SQR(B);
  CY[3]:=SQRT(B);
  IF N = 3 THEN
     EXIT(COORDS);
  FOR I:=4 TO N DO BEGIN
   A:=GAMMA[1,I];
   D:=GAMMA[2,I];
   CZ[I]:=(SQR(A)-SQR(D)+C)/CZ[2]*0.5;
   A:=CZ[2]-CZ[I];
   D:=CZ[3]-CZ[I];
   A:=SQR(D)-SQR(A);
   D:=GAMMA[3,I];
   A:=A-SQR(D);
   D:=GAMMA[2,I];
   CY[I]:=(A+SQR(D)+B)/CY[3]*0.5
  END;
  IF PM = 0 THEN BEGIN
     FOR I:=4 TO N DO
      CX[I]:=0;
     EXIT(COORDS);
  END;
  MAX:=0;
  FOR I:=4 TO N DO BEGIN
   A:=GAMMA[1,I];
   B:=CZ[I];
   C:=CY[I];
   D:=SQR(A)-SQR(B)-SQR(C);
   IF D < 0 THEN BEGIN
      D:=-D;
      WRITELN(FO,'INACCURATE GEOMETRY ')
   END;
   CX[I]:=SQRT(D);
   IF CX[I] > MAX THEN BEGIN
      MAX:=CX[I];
      J:=I
   END
  END;
  IF N = 4 THEN
     EXIT(COORDS);
  FOR I:=4 TO N DO
   IF I <> J THEN BEGIN
      B:=CY[J]-CY[I];
      C:=CZ[J]-CZ[I];
      D:=GAMMA[I,J];
      D:=SQR(D)-SQR(B)-SQR(C);
      A:=CX[J]-CX[I];
      A:=D-SQR(A);
      B:=CX[J]+CX[I];
      B:=D-SQR(B);
      IF SQR(A) > SQR(B) THEN
         CX[I]:=-CX[I]
   END;
END;

BEGIN   {REEDDATA}
 READ(N);
 READ (CONA);
 READ (CONB);
 READ (PM);
 READ (CI);
 IF (CI = 1) OR (CI = 0) THEN BEGIN
    NU:=N;
    NL:=1
 END ELSE BEGIN
    READ(NU);
    READ(NL)
 END;
 IF CONA = 2 THEN
    NF:= N DIV 2
 ELSE
    READ (NF);
 IF CONA <> 3 THEN BEGIN
    CONV:=0.0001;
    NIT:=15
 END ELSE BEGIN
    READ(NIT);
    READ(CONV)
 END;
 IF CONA = 2 THEN BEGIN
    READ(X);
    READ(Y);
    FOR I:=1 TO N DO BEGIN
     AI[I]:=X;
     E[I]:=Y;
     Z[I]:=1
    END
 END ELSE BEGIN
    REEDSM(AI);
    REEDSM(E);
    REEDSM(Z)
 END;
 FOR I:=1 TO N DO begin
  CX[i]:=0;
  CY[i]:=0;
  CZ[i]:=0;
  FOR J:=I TO N DO BEGIN
   BETA[I,J]:=0;
   BETA[J,I]:=0;
  END;
 end;
 READ(K);
 WHILE K<>0 DO BEGIN
  IF K = -1 THEN BEGIN
     READ(A);
     REPEAT
      READ(K);
      IF K <> -1 THEN BEGIN
        READ(FI,J);
        BETA[K,J]:=A;
        BETA[J,K]:=A
      END;
     UNTIL K = -1
  END ELSE BEGIN
     READ(J);
     READ(BETA[K,J]);
     BETA[J,K]:=BETA[K,J]
  END;
  READ(K)
 END;
 IF CONB = 1 THEN
  SHAPE(GAMMA);
 IF CONB = 2 THEN BEGIN
    FOR I:=1 TO N DO BEGIN
     GAMMA[I,I]:=0;
     J:=I+1;
     WHILE J <= N DO BEGIN
      READ(RD);
      GAMMA[I,J]:=RD;
      GAMMA[I,J]:=RD;
      J:=J+1
     END
 END;
 CX[1]:=0;
 CY[1]:=0;
 CZ[1]:=0;
 CX[2]:=0;
 CY[2]:=0;
 CZ[2]:=GAMMA[1,2];
 IF N <> 2 THEN
    COORDS
 END;
 IF CONB > 2 THEN BEGIN
    IF PM =1 THEN
       REEDSM(CX)
    ELSE
     FOR I:=1 TO N DO
      CX[I]:=0;
    REEDSM (CY);
    REEDSM(CZ);
    FOR I:= 1 TO N DO BEGIN
     GAMMA[I,I]:=0;
     J:=I+1;
     WHILE J <= N DO BEGIN
      IF PM = 0 THEN
         X:=0
      ELSE
         X:=CX[I]-CX[J];
      Y:=CY[I]-CY[J];
      Z1:=CZ[I]-CZ[J];
      GAMMA[I,J]:=SQRT(SQR(X)+SQR(Y)+SQR(Z1));
      GAMMA[J,I]:=GAMMA[I,J];
      J:=J+1
     END
    END
  END;
  NG:=(NF+1-NL)*(NU-NF)
END;

PROCEDURE RITEA(VAR X:ARRN;S:ST80;N:INTEGER);FORWARD;

PROCEDURE RITEV(VAR X:VECTN;S:ST80);FORWARD;

SEGMENT PROCEDURE RITEDATA;
BEGIN
 IF NG > M THEN BEGIN
   WRITELN('NUMBER OF TRANSITIONS EXCEED MAX. ALLOWED (',M:2,')');
   EXIT(PROGRAM);
 END ELSE
 WRITELN(FO,'NUMBER OF CENTRES  =  ',N);
 WRITELN(FO,'NUMBER OF DOUBLY OCCUPIED SHELLS  =  ',NF);
 IF CI <> 0 THEN BEGIN
    WRITELN(FO,'UPPER LIMIT OF CI  ',NU);
    WRITELN(FO,'LOWER LIMIT OF CI  ',NL);
 END;
 WRITELN(FO,'MAX NUMBER OF SCF ITERATIONS  ',NIT);
 NEWLINE(2);
 WRITE(FO,'CONVERGENCE CRITERION =  ',CONV:10:7);
 RITEA(BETA,'RESONANCE INTEGRAL MATRIX',N);
 RITEA(GAMMA,'INTER-ATOMIC DISTANCE MATRIX',N);
 WRITELN(FO,'ATOMIC COORDINATES');
 WRITELN(FO);
 RITEV(CX,'X-COORDS');
 RITEV(CY,'Y-COORDS');
 RITEV(CZ,'Z-COORDS');
 RITEV(AI,'ATOMIC IONISATION POTENTIALS');
 RITEV(E,'ATOMIC ELECTRON AFFINITIES');
 RITEV(Z,'ATOMIC VIRTUAL CHARGES')
END;

SEGMENT PROCEDURE DIAG(VAR A,S:ARRN;VAR EIG:VECTN;N:INTEGER);
VAR I,P,Q,J:INTEGER;
    APP,AQQ,APQ,API,AQI,SPI,SQI,CO2,SI2,SICO,SICO2,
    SINE,COSINE,OMEGA,NU,NUFINAL,MU:REAL;
    NOTOFFDIAG:BOOLEAN;

 PROCEDURE REORD(VAR A,S:ARRN);
 VAR I,J,P:INTEGER;
     APP:REAL;
 BEGIN
  FOR I:=1 TO N DO WRITE(CHR(7));
  FOR I:=1 TO N DO
   EIG[I]:=A[I,I];
  FOR I:=1 TO N-1 DO BEGIN
   APP:=EIG[I];
   P:=I;
    FOR J:=I+1 TO N DO
     IF APP >= EIG[J] THEN BEGIN
        APP:=EIG[J];
        P:=J;
     END;
   EIG[P]:=EIG[I];
   EIG[I]:=APP;
   FOR J:=1 TO N DO BEGIN
    APP:=S[I,J];
    S[I,J]:=S[P,J];
    S[P,J]:=APP
   END
  END;
  FOR I:=1 TO N DO
   FOR J:=I TO N DO BEGIN
    A[I,J]:=S[I,J];
    A[J,I]:=S[J,I]
   END
 END;

BEGIN
 NU:=0;
 FOR I:=1 TO N DO BEGIN
  S[I,I]:=1;
  J:=I+1;
  WHILE J <= N DO BEGIN
   S[I,J]:=0;
   S[J,I]:=0;
   NU:=NU+SQR(A[I,J])+SQR(A[J,I]);
   J:=J+1
  END
 END;
 NU:=SQRT(NU+NU);
 NUFINAL:=NU*5.0E-9/N;
 REPEAT
  NU:=NU/N;
  REPEAT
   NOTOFFDIAG:=TRUE;
   FOR Q:=1 TO N-1 DO
    FOR P:=Q+1 TO N DO BEGIN
     APQ:=A[P,Q];
     IF ABS(APQ) > NU THEN BEGIN
        NOTOFFDIAG:=FALSE;
        APP:=A[P,P];
        AQQ:=A[Q,Q];
        MU:=(APP-AQQ)*0.5;
        IF MU = 0 THEN BEGIN
           SINE:=-0.70710678;
           COSINE:=0.70710678;
           SI2:=0.50
        END ELSE BEGIN
           OMEGA:=-APQ/SQRT(SQR(APQ)+SQR(MU));
           IF MU < 0 THEN OMEGA := -OMEGA;
           SINE:=OMEGA/SQRT(2.0*(1.0+SQRT(1.0-SQR(OMEGA))));
           SI2:=SQR(SINE);
           COSINE:=SQRT(1.0-SI2)
        END;
        CO2:=SQR(COSINE);
        FOR I := 1 TO N DO BEGIN
         IF (I <> P) AND (I <> Q) THEN BEGIN
            API:=A[P,I];
            AQI:=A[Q,I];
            A[Q,I]:=API*SINE+AQI*COSINE;
            A[I,Q]:=A[Q,I];
            A[P,I]:=API*COSINE-AQI*SINE;
            A[I,P]:=A[P,I]
         END;
         SPI:=S[P,I];
         SQI:=S[Q,I];
         S[Q,I]:=SPI*SINE+SQI*COSINE;
         S[P,I]:=SPI*COSINE-SQI*SINE
        END;
        SICO:=SINE*COSINE;
        SICO2:=2.0*SICO*APQ;
        A[P,P]:=APP*CO2+AQQ*SI2-SICO2;
        A[Q,Q]:=APP*SI2+AQQ*CO2+SICO2;
        A[P,Q]:=(APP-AQQ)*SICO+APQ*(CO2-SI2);
        A[Q,P]:=A[P,Q]
     END
    END;
  UNTIL NOTOFFDIAG;
 UNTIL NU <= NUFINAL;
 REORD(A,S)
END;
(*$R+*)
SEGMENT PROCEDURE BONDER(VAR P:ARRN);
VAR EN:REAL;
  I,J,K:INTEGER;
 BEGIN
  FOR I:=1 TO N DO
   FOR J:=I TO N DO BEGIN
    EN:=0;
    FOR K:=1 TO NF DO
     EN:=EN+C[K,I]*C[K,J];
    P[I,J]:=EN+EN;
    P[J,I]:=P[I,J]
   END
  END;

PROCEDURE RITEOUT(VAR P:ARRN);FORWARD;

SEGMENT PROCEDURE HUCKEL(VAR PRC:ARRN);
VAR I,J:INTEGER;
BEGIN
 FOR I:=1 TO N DO
  FOR J:=I TO N DO BEGIN
   C[I,J]:=BETA[I,J];
   C[J,I]:=C[I,J]
  END;
  WRITELN('HUCKEL APPROXIMATION');
  DIAG(C,PRC,EIG,N);
  BONDER(P);
  RITEOUT(P);
  FOR I:=1 TO N DO BEGIN
   GAMMA[I,I]:=AI[I]-E[I];
   J:=I+1;
   WHILE J <= N DO BEGIN
    GAMMA[I,J]:=1.0/((GAMMA[I,J]/14.41)+(2/(AI[I]-E[I]+AI[J]-E[J])));
    GAMMA[J,I]:=GAMMA[I,J];
    J:=J+1
   END
  END
END;

SEGMENT PROCEDURE SCF;
LABEL 100;
VAR I,J,L,M,K:INTEGER;
    SUM:REAL;
BEGIN
 FOR K:=1 TO NIT DO BEGIN
  FOR I:=1 TO N DO BEGIN
   SUM:=0;
   FOR L:=1 TO N DO
    SUM:=SUM+(P[L,L]-Z[L])*GAMMA[I,L];
   C[I,I]:=-AI[I]-GAMMA[I,I]*(0.5*P[I,I]-Z[I])+SUM;
   J:=I+1;
   WHILE J <= N DO BEGIN
    C[I,J]:=BETA[I,J]-0.5*P[I,J]*GAMMA[I,J];
    C[J,I]:=C[I,J];
    J:=J+1
   END
  END;
  DIAG(C,PRC,EIG,N);
  BONDER(Q);
(* CONVERGENCE TEST *);
  FOR I:=1 TO N DO
   FOR J:=I TO N DO
    IF ABS(P[I,J]-Q[I,J]) > CONV THEN BEGIN
       FOR M:=1 TO N DO
        FOR L:=M TO N DO BEGIN
         P[M,L]:=Q[M,L];
         P[L,M]:=Q[M,L]
        END;
       GOTO 100
    END;
  WRITE(FO,'BOND MATRIX CONVERGED TO REQUIRED');
  WRITELN(FO,' EXTENT,NUMBER OF ITERATIONS  ',K);
  EXIT(SCF);
 100:
 END;
 WRITE(FO,'BOND MATRIX NOT CONVERGED TO REQUIRED');
 WRITELN(FO,' EXTENT,NUMBER OF ITERATIONS  ',NIT);
END;

SEGMENT FUNCTION PATAN(A,B:REAL):REAL;
VAR C:REAL;
BEGIN
 IF ABS(A) < 1E-6 THEN BEGIN
    PATAN:=0;
    EXIT(PATAN);
 END;
 IF ABS(B) < 1E-6 THEN
    C:=1.570796
 ELSE
    C:= ATAN(ABS(A / B));
 IF B < 0 THEN BEGIN
    IF A > 0 THEN
       C:=3.141593-C
    ELSE
       C:=3.141593+C
 END ELSE
    IF A < 0 THEN
       C:=6.283186-C;
 PATAN:=C*57.29578;
END;

SEGMENT PROCEDURE SINGTRIP;
(*$R PATAN*)
VAR M,P,I,MA,J,K,L:INTEGER;
    A,B,Q,X,X1,F:REAL;
    SING,TRIP:ARRN;
    QX,QY,QZ:VECTN;

 PROCEDURE HEAD;
  BEGIN
   NEWLINE(2);
   WRITE(FO,'TRANS.  SINGLET  TRIPLET  TRANS.');
   WRITE(FO,' DIPOLE  ALPHA  BETA  OSC. STRENGTH');
   WRITELN(FO,'  ABSORPTION');
   NEWLINE(2)
  END;

 PROCEDURE SPEC(S,T,Q,A,B,F:REAL);
 BEGIN
  WRITE(FO,S:10:4,T:9:4,Q:12:4);
  WRITE(FO,A:10:2,B:6:2,F:12:4);
  WRITELN(FO,1240 / S:14:4)
 END;

 PROCEDURE COMPSNTP(K,L:INTEGER;VAR C,GAMMA:ARRN);
 VAR F,C1,X:REAL;
     P,R:INTEGER;
 BEGIN
  A:=0;
  B:=0;
  FOR P:=1 TO N DO BEGIN
   Q:=C[I,P];
   C1:=0;
   F:=0;
   FOR R:=1 TO N DO BEGIN
    X:=C[L,R]*GAMMA[P,R];
    C1:=C1+X*C[J,R];
    F:=F+X*C[K,R];
   END;
   A:=A+C1*Q*C[K,P];
   B:=B+F*Q*C[J,P];
  END;
 END;

PROCEDURE MAKECI(VAR SING,TRIP:ARRN;VAR QX,QY,QZ:VECTN);
VAR J,P:INTEGER;
    A,B,Q,C1,X,F:REAL;
    ES,ET:VECTN;
BEGIN
 WRITELN(FO);
 WRITELN(FO);
 WRITELN(FO,'NUMBER OF STATES USED IN CI', NG);
 DIAG(SING,PRC,ES,NG);
 DIAG(TRIP,PRC,ET,NG);
 NEWLINE(2);
 WRITELN(FO,'TRANSITIONS AFTER CI');
 HEAD;
 FOR J:=1 TO NG DO BEGIN
  A:=0;
  IF PM <> 0 THEN
     FOR P:=1 TO NG DO
      A:=A+SING[J,P]*QX[P];
  B:=0;
  C1:=0;
  FOR P:=1 TO NG DO BEGIN
   X:=SING[J,P];
   B:=B+X*QY[P];
   C1:=C1+X*QZ[P];
  END;
  F:=SQR(B)+SQR(C1);
  Q:=2*(SQR(A)+F);
  B:=PATAN(B,C1);
  IF PM = 0 THEN
     A:=0
  ELSE
     A:=PATAN(A,SQRT(F));
  F:=0.0875161*Q*ES[J];
  WRITE(FO,J:3,'  ');
  SPEC(ES[J],ET[J],Q,B,A,F);
 END;
 NEWLINE(2);
 RITEA(SING,'SINGLET CONFIGURATION VECTORS  ',NG);
 RITEA(TRIP,'TRIPLET CONFIGURATION VECTORS  ',NG);
END;

BEGIN
 NEWLINE(2);
 WRITELN(FO,'SINGLE ELECTRON EXCITATIONS');
 NEWLINE(2);
 WRITELN(FO,'TRANSITIONS BEFORE CI');
 HEAD;
 M:=0;
 X1:=1E10;
 FOR I:=NL TO NF DO
  FOR J:=NF+1 TO NU DO BEGIN
   M:=M+1;
   MA:=0;
   FOR K:=NL TO NF DO
    FOR L:=NF+1 TO NU DO BEGIN
     MA:=MA+1;
     IF MA >= M THEN
        IF (CI <> 0) OR (MA = M) THEN BEGIN
           COMPSNTP(K,L,C,GAMMA);
           SING[M,MA]:=-A+B+B;
           SING[MA,M]:=SING[M,MA];
           TRIP[M,MA]:=-A;
           TRIP[MA,M]:=-A;
           IF MA = M THEN BEGIN
              A:=EIG[J]-EIG[I];
              SING[M,M]:=SING[M,M]+A;
              TRIP[M,M]:=TRIP[M,M]+A;
              Q:=0;
              IF PM <> 0 THEN BEGIN
                 FOR P:=1 TO N DO
                  Q:=Q+C[I,P]*C[J,P]*CX[P];
                 QX[M]:=Q
              END;
              A:=0;
              B:=0;
              FOR P:=1 TO N DO BEGIN
               X:=C[I,P]*C[J,P];
               A:=A+X*CY[P];
               B:=B+X*CZ[P]
              END;
              QY[M]:=A;
              QZ[M]:=B;
              B:=SQR(A)+SQR(B);
              Q:=2*(SQR(Q)+B);
              A:=PATAN(QY[M],QZ[M]);
              IF PM = 0 THEN
                 B:=0
              ELSE
                 B:=PATAN(QX[M],SQRT(B));
              F:=0.0875161*Q*SING[M,M];
              WRITE(FO,I:2,',',J:2);
              SPEC(SING[M,M],TRIP[M,M],Q,A,B,F);
              IF X1 > SING[M,M] THEN BEGIN
                 X1:=SING[M,M];
                 STI:=I;
                 STJ:=J
              END;
           END
        END
    END
  END;
  FOR I:=1 TO N DO
   EIG[I]:=SQR(PRC[STJ,I])-SQR(PRC[STI,I]);
  IF CI <> 0 THEN
     MAKECI(SING,TRIP,QX,QY,QZ);
END;

SEGMENT PROCEDURE DENSITY;
 VAR I:INTEGER;
  BEGIN
   WRITELN(FO,'ELECTRON DENSITY CHANGES');
   NEWLINE(2);
   WRITELN(FO,'      ROWS USED = ',STI:5,',',STJ:5);
   NEWLINE(2);
   FOR I:=1 TO N DO
    BEGIN
     WRITE(FO,EIG[I]:10:6);
     IF I MOD 8 = 0 THEN WRITELN(FO)
    END;
   WRITELN(FO)
  END;

SEGMENT PROCEDURE TOTEN;
VAR I,J,L:INTEGER;
    EN,SUM,SUM1:REAL;

 PROCEDURE DIPOLE;
(*$R PATAN*)
 VAR I:INTEGER;
     A,B,C,D,E:REAL;
 BEGIN
  A:=0;
  B:=0;
  C:=0;
  FOR I:=1 TO N DO BEGIN
   D:=Z[I]-P[I,I];
   A:=A+CX[I]*D;
   B:=B+CY[I]*D;
   C:=C+CZ[I]*D;
  END;
  A:=A*4.8;
  B:=B*4.8;
  C:=C*4.8;
  D:=SQR(B)+SQR(C);
  NEWLINE(2);
  WRITELN(FO,'PI-CONTRIBUTION TO DIPOLE MOMENT');
  NEWLINE(2);
  WRITE(FO,'X-MOMENT= ',A:8:2,'     Y-MOMENT = ',B:8:2);
  WRITELN(FO,'     Z-MOMENT = ',C:8:2);
  NEWLINE(2);
  E:=D+SQR(A);
  C:=PATAN(B,C);
  B:=PATAN(A,SQRT(D));
  WRITE(FO,'TOTAL MOMENT = ',E:8:2);
  WRITELN(FO,'     ALPHA = ',C:8:2,'     BETA = ',B:8:2);
  NEWLINE(3)
 END;

BEGIN
 FOR I:=1 TO N DO BEGIN
  SUM:=0;
  FOR L:=1 TO N DO
   SUM:=SUM+(Q[L,L]-Z[L])*GAMMA[I,L];
  C[I,I]:=-AI[I]-GAMMA[I,I]*(0.5*Q[I,I]-Z[I])+SUM;
  J:=I+1;
  WHILE J <= N DO BEGIN
   C[I,J]:=BETA[I,J]-0.5*Q[I,J]*GAMMA[I,J];
   C[J,I]:=C[I,J];
   J:=J+1
  END
 END;
 EN:=0;
 FOR I:=1 TO N DO BEGIN
  SUM1:=0;
  SUM:=0;
  FOR J:=1 TO N DO BEGIN
   SUM:=SUM+Z[J]*GAMMA[I,J];
   SUM1:=SUM1+Q[I,J]*(BETA[I,J]+C[I,J])
  END;
  EN:=EN+SUM1+Q[I,I]*(-BETA[I,I]-AI[I]+
      Z[I]*GAMMA[I,I]-SUM)
 END;
 EN:=EN*0.5;
 NEWLINE(2);
 WRITELN(FO,'TOTAL ENERGY = ',EN:12:5);
 DIPOLE
END;

PROCEDURE NEWLINE(NLIN:INTEGER);
  VAR I:INTEGER;
  BEGIN
   FOR I:=1 TO NLIN DO
    WRITELN(FO)
  END;

PROCEDURE RITEA;
  VAR I,J:INTEGER;
BEGIN
 NEWLINE (2);
 WRITELN(FO,S);
 FOR I:=1 TO N DO BEGIN
  FOR J:=1 TO N DO BEGIN
   WRITE(FO,X[I,J] :10:6);
   IF J MOD 8 = 0 THEN WRITELN(FO)
  END;
   WRITELN(FO)
 END;
 NEWLINE (2)
END;

PROCEDURE RITEV;
  VAR I:INTEGER;
BEGIN
 NEWLINE(2);
 WRITELN(FO,S);
 FOR I:=1 TO N DO BEGIN
  WRITE(FO,X[I] :10:6);
  IF I MOD 8 = 0 THEN WRITELN(FO);
 END;
 NEWLINE(2)
END;

PROCEDURE RITEOUT(VAR P:ARRN);
BEGIN
 NEWLINE(2);
 RITEV(EIG,'EIGENVALUES');
 RITEA(C,'EIGENVECTOR',N);
 RITEA(P,'BOND-ORDER MATRIX',N);
END;

(* MAIN PROGRAM *)

BEGIN
 HEADER;
 REEDDATA;
 RITEDATA;
 HUCKEL(PRC);
 SCF;
 RITEOUT(Q);
 SINGTRIP;
 DENSITY;
 TOTEN
END.