1!
2!  Dalton, a molecular electronic structure program
3!  Copyright (C) by the authors of Dalton.
4!
5!  This program is free software; you can redistribute it and/or
6!  modify it under the terms of the GNU Lesser General Public
7!  License version 2.1 as published by the Free Software Foundation.
8!
9!  This program is distributed in the hope that it will be useful,
10!  but WITHOUT ANY WARRANTY; without even the implied warranty of
11!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12!  Lesser General Public License for more details.
13!
14!  If a copy of the GNU LGPL v2.1 was not distributed with this
15!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
16!
17!
18C
19!===========================================================================
20! look for HJMAERKE; (1) symort of UMAT should be done
21!
22!961213-hjaaj:
23!AVESET: disable SUPSYM if no supersymmetry found
24!        HJMAERKE: this means that debug SUPSYM for abelian group
25!        cannot be performed without commenting new code in AVESET!
26!940815-hjaaj:
27!AVESE1,AVEDEG,AVEGTU,AVEPHA: use THRSSY from INFINP (new in INFINP)
28!   instead of THREQL parameter; include INFINP
29!AVESE1: Check and reset THRSSY if needed; more print if disentangling
30!940427-hjaaj: SCALAR directives for Cray (as SCALSUB does not work
31!  for Cray because directive has to be after SUBROUTINE statement)
32!930715-hjaaj
33!AVESE0: define MXDGSS and NINFSS(*,3)
34!AVEDEG: define MXDGSS
35!===========================================================================
36C  /* Deck averag */
37      SUBROUTINE AVERAG (VECS,NVDIM,NVSIM)
38C
39C 7-Jul-1992 hjaaj+hh
40C
41#include "implicit.h"
42#include "priunit.h"
43      DIMENSION VECS(NVDIM,NVSIM)
44C
45C Used from common blocks:
46C  INFINP : SUPSYM
47C  INFVAR : NWOPT,JWOPSY
48C
49#include "maxorb.h"
50#include "infinp.h"
51#include "infvar.h"
52#include "infpri.h"
53C
54      CALL QENTER('AVERAG')
55      IF (.NOT. SUPSYM) GO TO 9999
56      IF (NWOPT  .EQ. 0) GO TO 9999
57      IF (JWOPSY .NE. 1) GO TO 910
58      IF (NVDIM  .LT. NWOPT) GO TO 920
59      IF (SUPSYM) THEN
60         CALL AVERSS(VECS,NVDIM,NVSIM)
61      END IF
62 9999 CALL QEXIT('AVERAG')
63      RETURN
64C
65C *** Error sections
66C
67  910 CONTINUE
68      WRITE (LUERR,9010) JWOPSY
69      CALL QTRACE(LUERR)
70      CALL QUIT('AVERAG ERROR, operator symmetry JWOPSY .ne. 1')
71 9010 FORMAT(/' ERROR-AVERAG, operator symmetry .ne. 1; JWOPSY =',I8)
72C
73  920 CONTINUE
74      WRITE (LUERR,9020) NVDIM,NWOPT
75      CALL QTRACE(LUERR)
76      CALL QUIT('AVERAG ERROR, vector length lt NWOPT')
77 9020 FORMAT(/' ERROR-AVERAG, vector length',I8,' is less than NWOPT',
78     *   I8)
79C
80      END
81C  /* Deck averss */
82      SUBROUTINE AVERSS (VECS,NVDIM,NVSIM)
83C
84C 7-Jul-1992 Hans Joergen Aa. Jensen + Hinne Hettema
85C
86C Purpose: Average degenerate super symmetries
87C          (which are different components of a degenerate irrep).
88C          This corresponds to doing a microcanonical average.
89C
90#include "implicit.h"
91      DIMENSION VECS(NVDIM,NVSIM)
92C
93C Used from common blocks:
94C   INFORB : NSSYM,NORBSS(),IORBSS(),NINFSS(*,3), NOCCT,NORBT
95C   INFIND : ISW(*),ISSORD(*)
96C   INFVAR : NWOPT
97C
98#include "maxash.h"
99#include "maxorb.h"
100#include "priunit.h"
101#include "infind.h"
102#include "inforb.h"
103#include "infvar.h"
104#include "infpri.h"
105C
106      PARAMETER (MAXDEG = 20, D1 = 1.0D0)
107      DIMENSION ISSDEG(MAXDEG),ISSOFF(MAXDEG),IDEGV(MAXDEG)
108C
109Chjaaj-may2000:
110C     Dynamic allocation of KLWOP with adjustable run-time dimensions
111      DIMENSION KLWOP(NOCCT,NORBT)
112C
113      CALL QENTER('AVERSS')
114      CALL MAKE_KLWOP(KLWOP)
115      IF (IPRAVE .GT. 100) THEN
116         WRITE(LUPRI,'(/A)') ' AVERSS, Incoming orbital vector(s):'
117         CALL OUTPUT(VECS,1,NWOPT,1,NVSIM,NVDIM,NVSIM,-1,LUPRI)
118      END IF
119C
120      DO 800 ISSYMR = 1,NSSYM
121C        Skip this supsym if not degenerate or not "root supsym"
122         IF (NINFSS(ISSYMR,1) .EQ. 1) GO TO 800
123         IF (NINFSS(ISSYMR,2) .NE. ISSYMR) GO TO 800
124         NDEG = 0
125         DO 100 ISSYM = ISSYMR,NSSYM
126            IF (NINFSS(ISSYM,2) .EQ. ISSYMR) THEN
127               NDEG = NDEG + 1
128               ISSDEG(NDEG) = ISSYM
129               ISSOFF(NDEG) = IORBSS(ISSYM) - IORBSS(ISSYMR)
130            END IF
131  100    CONTINUE
132         IF (IPRAVE .GT. 5 .OR. NDEG .NE. NINFSS(ISSYMR,1)) THEN
133            WRITE(LUPRI,'(/A,I3,A/A,(T35,8I5))')
134     &         ' AVERSS averaging ',NDEG,' degenerate super symmetries',
135     &         ' namely sup.sym.s no.',(ISSDEG(IDEG),IDEG=1,NDEG)
136         END IF
137         IF (NDEG .NE. NINFSS(ISSYMR,1)) THEN
138            WRITE(LUPRI,'(/A)')
139     &      ' AVERSS-software error: inconsistent degeneracy in NINFSS'
140            CALL AVEDMP(LUPRI)
141            CALL QUIT('AVERSS: inconsistent degeneracy in NINFSS')
142         END IF
143C
144         AVFAC = NDEG
145         AVFAC = D1 / AVFAC
146C
147         IOFFSR = IORBSS(ISSYMR)
148         DO 480 IMOSS = IOFFSR+1,IOFFSR+NORBSS(ISSYMR)
149            IMO = ISSORD(IMOSS)
150            IMOW = ISW(IMO)
151            IF (IMOW .GT. NOCCT) GO TO 480
152            DO 470 JMOSS = IMOSS+1,IOFFSR+NORBSS(ISSYMR)
153               JMO = ISSORD(JMOSS)
154               JMOW = ISW(JMO) - NISHT
155               IF (JMOW .LE. 0) GO TO 470
156Chjaaj         ... this is an inactive orbital
157               IF (KLWOP(IMOW,JMOW) .EQ. 0) GO TO 470
158               IDEGV(1) = KLWOP(IMOW,JMOW)
159               NERR = 0
160               DO 410 IDEG = 2,NDEG
161                  KMO = ISSORD(IMOSS + ISSOFF(IDEG))
162                  LMO = ISSORD(JMOSS + ISSOFF(IDEG))
163                  LMOW = ISW(LMO) - NISHT
164                  IF (LMOW .LE. 0) THEN
165                     IDEGV(IDEG) = 0
166                     NERR = NERR + 1
167                     WRITE(LUPRI,'(/A)')
168     &       'AVERSS error: degenerate orbitals not same orbital class.'
169                  ELSE
170                     IDEGV(IDEG) = KLWOP( ISW(KMO) , LMOW )
171                     IF (IDEGV(IDEG) .EQ. 0) NERR = NERR + 1
172                  END IF
173  410          CONTINUE
174               IF (IPRAVE .GT. 10 .OR. NERR .GT. 0) THEN
175                  IF (NERR .GT. 0) THEN
176                     WRITE(LUPRI,'(/A/A/)')
177     &         ' AVERSS error: not all degenerate rotations included',
178     &         ' Probable error: some but not all'//
179     &            ' of the degenerate orbitals have been frozen.'
180                  END IF
181                  WRITE(LUPRI,'(A,(T35,4I10))')
182     &               ' AVERSS averaging orb.rot. no.s',
183     &               (IDEGV(IDEG),IDEG=1,NDEG)
184                  IF (NERR .GT. 0) CALL AVEDMP(LUPRI)
185               END IF
186               IF (NERR .GT. 0) THEN
187                  CALL QTRACE(LUPRI)
188                  CALL QUIT(
189     &            'AVERSS error: not all degenerate rotations included')
190               END IF
191               DO 460 IVSIM = 1,NVSIM
192                  XKAPIJ = VECS(IDEGV(1),IVSIM)
193                  DO 420 IDEG = 2,NDEG
194                     XKAPIJ = XKAPIJ + VECS(IDEGV(IDEG),IVSIM)
195  420             CONTINUE
196                  XKAPIJ = AVFAC*XKAPIJ
197                  DO 430 IDEG = 1,NDEG
198                     VECS(IDEGV(IDEG),IVSIM) = XKAPIJ
199  430             CONTINUE
200  460          CONTINUE
201  470       CONTINUE
202  480    CONTINUE
203  800 CONTINUE
204C
205      IF (IPRAVE .GT. 100) THEN
206         WRITE(LUPRI,'(/A)') ' AVERSS, Averaged orbital vector(s):'
207         CALL OUTPUT(VECS,1,NWOPT,1,NVSIM,NVDIM,NVSIM,-1,LUPRI)
208      END IF
209C
210      CALL QEXIT('AVERSS')
211      RETURN
212      END
213C  /* Deck aveset */
214      SUBROUTINE AVESET(CMO,WRK,KFREE,LFREE)
215C
216C 2-Jul-1992 hjaaj : core allocation for AVESE1
217C
218#include "implicit.h"
219      REAL*8  CMO(*), WRK(*)
220C
221C Used from common blocks:
222C  INFINP: SUPSYM
223C  INFORB: N2ORBX, MXSSYM
224C
225#include "maxorb.h"
226#include "priunit.h"
227#include "infinp.h"
228#include "inforb.h"
229#include "infpri.h"
230C
231      CALL QENTER('AVESET')
232      IF (SUPSYM) THEN
233         IF (IPRAVE .GE. 3) THEN
234            WRITE(LUPRI,'(//A)') ' TEST OUTPUT FROM AVESET'
235         END IF
236         KFRSAV = KFREE
237         LFRSAV = LFREE
238         CALL MEMGET('REAL',KTKMO,N2ORBX,WRK,KFREE,LFREE)
239         CALL MEMGET('INTE',KNDEGS,MXSSYM,WRK,KFREE,LFREE)
240         CALL AVESE1(CMO,WRK(KTKMO),WRK(KNDEGS),
241     *               WRK,KFREE,LFREE)
242         CALL MEMREL('AVESET',WRK,KFRSAV,KFRSAV,KFREE,LFREE)
243C
244C        Check if supersymmetry is identical to "D2h" symmetry
245C
246         IF (NSSYM .LE. NSYM .AND. MXDGSS .EQ. 1) THEN
247            DO 200 ISYM = 1,NSYM
248               JSYM = 0
249               DO 100 ISSYM = 1,NSSYM
250                  IF (NINFSS(ISSYM,3) .EQ. ISYM) JSYM = JSYM + 1
251  100          CONTINUE
252               IF (JSYM .GT. 1) GO TO 300
253  200       CONTINUE
254            IF (IPRAVE .GE. 1) THEN
255               WRITE(LUPRI,'(/A)')' AVESET: No "super symmetry" found.'
256            END IF
257            SUPSYM = .FALSE.
258            GO TO 300
259         END IF
260  300    CONTINUE
261      END IF
262      IF (.NOT.SUPSYM) THEN
263         CALL AVESE0
264      END IF
265      CALL QEXIT('AVESET')
266      RETURN
267      END
268C  /* Deck avese0 */
269      SUBROUTINE AVESE0
270C
271C 7-Jul-1992 hjaaj
272C Set super symmetry information as symmetry information
273C (when super symmetry is not to be used()
274C
275#include "implicit.h"
276C
277C Used from common blocks:
278C  INFORB : NSSYM, NORBSS(), IORBSS(); NSYM, NORB(), IORB(), NORBT
279C  INFIND : ISSMO(), ISSORD(); ISMO()
280C
281#include "maxash.h"
282#include "maxorb.h"
283#include "inforb.h"
284#include "infind.h"
285C
286      CALL QENTER('AVESE0')
287      NSSYM = NSYM
288      DO 100 ISYM = 1,NSYM
289         NORBSS(ISYM) = NORB(ISYM)
290         IORBSS(ISYM) = IORB(ISYM)
291         NINFSS(ISYM,1) = 1
292         NINFSS(ISYM,2) = ISYM
293         NINFSS(ISYM,3) = ISYM
294  100 CONTINUE
295      DO 200 IMO = 1,NORBT
296         ISSMO(IMO)  = ISMO(IMO)
297         ISSORD(IMO) = IMO
298  200 CONTINUE
299      MXDGSS = 1
300      CALL QEXIT('AVESE0')
301      RETURN
302      END
303C  /* Deck avese1 */
304      SUBROUTINE AVESE1(CMO,TKMO,NDEGSS,WRK,KFREE,LFREE)
305C
306C Written 29-Jun-1992 Hans Joergen Aa. Jensen and Hinne Hettema
307C
308C Purpose:
309C  Find super symmetry and degeneracy of molecular orbitals in CMO
310C  Disentangle mixed degenerate orbitals
311C  Coordinate relative phases of degenerate orbitals
312C
313C Input:
314C  CMO; molecular orbitals
315C
316C Output:
317C
318#include "implicit.h"
319      REAL*8    CMO(NCMOT), TKMO(NORBT,NORBT), WRK(*)
320      INTEGER   NDEGSS(*) ! DIMENSION NDEGSS(MXSSYM)
321C
322#include "threql.h"
323      PARAMETER (TMXSSY = 1.0D-4)
324C
325C Used from common blocks:
326C   INFINP : THRSSY
327C   INFORB : NSYM,NSSYM,NORBT,NORB(),?
328C   INFIND : ISSMO(),ISSORD(),?
329C   INFPRI : IPRAVE,NINFO,?
330C
331#include "maxash.h"
332#include "maxorb.h"
333#include "priunit.h"
334#include "infinp.h"
335#include "inforb.h"
336#include "infind.h"
337#include "infpri.h"
338C
339C
340      CALL QENTER('AVESE1')
341      IF (THRSSY .GT. TMXSSY) THEN
342         NWARN = NWARN + 1
343         WRITE (LUPRI,'(/A/2(A,1P,D10.2))') ' AVESE1-WARNING: SUPSYM'//
344     &      ' threshold for degeneracy and "zero" elements too big;',
345     &      '                 value reset from',THRSSY,' to',TMXSSY
346         THRSSY = TMXSSY
347      END IF
348      IF (THRSSY .LT. THREQL) THEN
349C        ... THREQL must be .gt. accuracy of diagonalization routine
350         NINFO = NINFO + 1
351         WRITE (LUPRI,'(/A/2(A,1P,D10.2))') ' AVESE1-INFO:'//
352     &      ' SUPSYM threshold for degeneracy and "zero" elements',
353     &      '              is reset from',THRSSY,' to',THREQL
354         THRSSY = THREQL
355      END IF
356      KFRSAV = KFREE
357      ITURN = 0
358  100 ITURN = ITURN + 1
359      KFREE1 = KFREE
360C
361C Step 0: Obtain TKMO matrix for finding super symmetries
362C         TKMO must be the matrix of a totally symmetric
363C         operator which is very unlikely to have
364C         accidental zeroes (we use the kinetic energy).
365C
366      CALL AVETK(TKMO,CMO,WRK,KFREE,LFREE)
367C     Note: this can give problems for FINITE FIELD (NFIELD.gt.0)
368C     because the TKMO will reflect the symmetry without the finite field(s)!
369C
370C Step 1: Find super symmetry of each MO
371C         Define ISSMO(i)  = super sym. of m.o. no. i
372C                ISSORD(i) = sup. sym. reordering array
373C                NSSYM     = no. of super symmetries
374C                NORBSS(*) = no. of mo's in each sup.sym.
375C                IORBSS(*) = pointer array to ISSORD
376C
377      CALL IZERO(ISSMO,NORBT)
378      CALL IZERO(ISSORD,NORBT)
379      NSSYMX = 0
380      DO 180 ISYM = 1,NSYM
381         IF (NORB(ISYM) .EQ. 0) GO TO 180
382         IOFFMO = IORB(ISYM)
383         DO 170 IMO = IOFFMO+1,IOFFMO+NORB(ISYM)
384            IF (ISSMO(IMO) .EQ. 0) THEN
385               NSSYMX = NSSYMX + 1
386               IF (NSSYMX .GT. MXSSYM) GO TO 901
387C     v-------- error exit
388               ISSMO(IMO) = NSSYMX
389               NORBSS(NSSYMX) = 1
390            END IF
391            ISSYM = ISSMO(IMO)
392            DO 160 JMO = IMO+1,IOFFMO+NORB(ISYM)
393            IF (ABS(TKMO(JMO,IMO)) .GT. THRSSY) THEN
394            IF (ISSMO(JMO) .EQ. 0) THEN
395C -- We have a new orbital in the supersymmetry rep.
396               ISSMO(JMO) = ISSYM
397               NORBSS(ISSYM) = NORBSS(ISSYM) + 1
398            ELSE IF (ISSMO(JMO) .NE. ISSYM) THEN
399C -- We have a conflict and must collapse the two supersym.s
400               JSSYMH = MAX(ISSMO(JMO),ISSYM)
401               JSSYML = MIN(ISSMO(JMO),ISSYM)
402               IF (IPRAVE .GE. 4) THEN
403                  WRITE(LUPRI,'(/A,2I4)')
404     *             ' Collapsing super symmetries :',JSSYML,JSSYMH
405               END IF
406               DO 140 KMO = IOFFMO+1,IOFFMO+NORB(ISYM)
407                  IF (ISSMO(KMO) .EQ. JSSYMH) ISSMO(KMO) = JSSYML
408  140          CONTINUE
409               ISSYM  = JSSYML
410               NORBSS(JSSYML) = NORBSS(JSSYML) + NORBSS(JSSYMH)
411               NORBSS(JSSYMH) = 0
412            END IF
413            END IF
414  160       CONTINUE
415  170    CONTINUE
416  180 CONTINUE
417C
418C  in case of debugging, get out some print
419C
420      IF ( IPRAVE .GE. 10) THEN
421         NSSYM = NSSYMX
422         CALL IZERO(IORBSS,NSSYM)
423         WRITE(LUPRI,'(/A/A/)')
424     *      ' After initializing in AVESE1 (step 1)',
425     *      ' ====================================='
426         CALL AVEDMP(LUPRI)
427      END IF
428C
429C Step 2: Find NSSYM, IORBSS(*) and ISSORD(*)
430C
431      NSSYM = 0
432      IMOSS = 0
433      IORBSS(1) = 0
434      DO 260 ISSYMX = 1,NSSYMX
435      IF (NORBSS(ISSYMX) .GT. 0) THEN
436         NSSYM = NSSYM + 1
437         IF (NSSYM .NE. ISSYMX) THEN
438            NORBSS(NSSYM) = NORBSS(ISSYMX)
439            DO 220 IMO = 1,NORBT
440               IF (ISSMO(IMO) .EQ. ISSYMX) ISSMO(IMO) = NSSYM
441  220       CONTINUE
442         END IF
443         IORBSS(NSSYM) = IMOSS
444         DO 240 IMO = 1,NORBT
445            IF (ISSMO(IMO) .EQ. NSSYM) THEN
446               IMOSS = IMOSS + 1
447               ISSORD(IMOSS) = IMO
448            END IF
449  240    CONTINUE
450         IF ( (IMOSS - IORBSS(NSSYM)) .NE. NORBSS(NSSYM) ) THEN
451            GO TO 902
452         END IF
453      END IF
454  260 CONTINUE
455C
456C  in case of debugging, get out some print
457C
458      IF ( IPRAVE .GE. 3) THEN
459         WRITE(LUPRI,'(/A/A/)')
460     *      ' After clearing up in AVESE1 (step 2)',
461     *      ' ===================================='
462         CALL AVEDMP(LUPRI)
463      END IF
464C
465C Step 3: Find degeneracies
466C
467      CALL MEMGET('INTE',KINDX1,NORBT,WRK,KFREE,LFREE)
468      CALL MEMGET('INTE',KINDX2,NORBT,WRK,KFREE,LFREE)
469      CALL AVEDEG(TKMO,WRK(KINDX1),WRK(KINDX2),NDEGSS)
470C
471C Step 4: Disentangle degenerate super symmetries
472C
473      NDIS = 0
474      DO 480 ISSYM = 1,NSSYM
475         NDEGI = NDEGSS(ISSYM)
476         IF (NDEGI .GT. 1) THEN
477            NDIS = NDIS + 1
478         IF (ITURN .EQ. 1) THEN
479            ISYM  = NINFSS(ISSYM,3)
480            NORBI = NORB(ISYM)
481            NBASI = NBAS(ISYM)
482            JCMOI = ICMO(ISYM) + 1
483            NSSMOI = NORBSS(ISSYM)
484            NSSMON = NSSMOI / NDEGI
485            IF (NDEGI*NSSMON .NE. NSSMOI) THEN
486               CALL QUIT('AVESE1 error: NDEGI*NSSMON .ne. NSSMOI')
487            END IF
488            KFREE4 = KFREE
489            CALL MEMGET('REAL',KUMAT,NDEGI*NDEGI,WRK,KFREE,LFREE)
490            CALL MEMGET('INTE',KIMODE,NSSMOI,WRK,KFREE,LFREE)
491            CALL MEMGET('REAL',KWRK,NBASI*NDEGI,WRK,KFREE,LFREE)
492            CALL AVEDIS(ISSYM,ISYM,NDEGI,NORBI,NBASI,NSSMOI,NSSMON,
493     *                  CMO(JCMOI),TKMO,WRK(KINDX1),
494     *                  WRK(KUMAT),WRK(KIMODE),WRK(KWRK))
495            CALL MEMREL('AVESE1.AVEDIS',WRK,1,KFREE4,KFREE,LFREE)
496         END IF
497         END IF
498  480 CONTINUE
499C
500      CALL MEMREL('AVESE1.step 4',WRK,1,KFREE1,KFREE,LFREE)
501C
502      IF (NDIS .GT. 0) THEN
503         IF (ITURN .EQ. 1) THEN
504            IF (IPRAVE .GE. 3) WRITE(LUPRI,'(/A)')
505     &         ' AVESE1: new super symmetry analysis after'//
506     &         ' disentangling degenerate super symmetries.'
507            GO TO 100
508         END IF
509         WRITE(LUPRI,'(/A)')
510     &      ' AVESE1 error: disentangling failed, sorry!'
511         CALL QUIT('AVESE1 error: disentangling failed, sorry!')
512      END IF
513C
514C Step 5: Check (and change) relative phases of MO's
515C
516      IF (IPRAVE .GE. 4) THEN
517         WRITE(LUPRI,'(/A)')
518     &      ' AVESE1 : checking relative phases of degenerate MOs'
519      END IF
520      IPHCHA = 0
521      DO 580 ISSYM = 1, NSSYM
522         IDEG   = NINFSS(ISSYM,1)
523         ISSYMR = NINFSS(ISSYM,2)
524         IF (IDEG .GT. 1 .AND. ISSYM .EQ. ISSYMR) THEN
525C        --- we have a degeneracy to other supersym.(s)
526            CALL AVEPHA(ISSYM,CMO,TKMO,IPHCHA)
527         END IF
528  580 CONTINUE
529      IF (IPRAVE .GE. 4 .AND. IPHCHA .GT. 0) THEN
530         WRITE(LUPRI,'(A,I5)')
531     &      ' AVESE1 : # of phase shifts of degenerate MOs :',IPHCHA
532      END IF
533C
534C *** end of subroutine AVESE1
535C
536      CALL QEXIT('AVESE1')
537      RETURN
538C
539  901 CONTINUE
540      CALL QUIT('AVESE1: NSSYM .gt. MXSSYM in AVESET code 901')
541  902 CONTINUE
542      CALL QUIT('AVESE1: Inconsistent NORBSS(NSSYM) code 902')
543      END
544C  /* Deck aveord */
545      SUBROUTINE AVEORD()
546C
547C Aug. 2004 hjaaj - create ISSORD() from ISSMO()
548C
549C The order of the super symmetry orbitals may be changed
550C in ORDRSS and ORD2SS calls, thus we need to remake ISSORD()
551C to be sure it is correct.
552C
553#include "implicit.h"
554#include "priunit.h"
555#include "maxash.h"
556#include "maxorb.h"
557#include "infind.h"
558#include "inforb.h"
559#include "infpri.h"
560C
561      CALL QENTER('AVEORD')
562      IMOSS = 0
563      IORBSS(1) = 0
564      DO ISSYM = 1,NSSYM
565         IF (IORBSS(ISSYM) .NE. IMOSS) THEN
566            CALL AVEDMP(LUPRI)
567            CALL QUIT('AVEORD: Inconsistent IORBSS(ISSYM)')
568         END IF
569         DO IMO = 1,NORBT
570            IF (ISSMO(IMO) .EQ. ISSYM) THEN
571               IMOSS = IMOSS + 1
572               ISSORD(IMOSS) = IMO
573            END IF
574         END DO
575         IF ( (IMOSS - IORBSS(ISSYM)) .NE. NORBSS(ISSYM) ) THEN
576            CALL AVEDMP(LUPRI)
577            CALL QUIT('AVEORD: Inconsistent NORBSS(ISSYM)')
578         END IF
579      END DO
580C
581C  In case of debugging, get out some print
582C
583      IF ( IPRAVE .GE. 3) THEN
584         WRITE(LUPRI,'(/A/A/)')
585     *      ' After remake of ISSORD() in AVEORD  ',
586     *      ' ===================================='
587         CALL AVEDMP(LUPRI)
588      END IF
589      CALL QEXIT('AVEORD')
590      RETURN
591      END
592C  /* Deck avetk */
593      SUBROUTINE AVETK(TKMO,CMO,WRK,KFREE,LFREE)
594C
595C 920629-hjaaj:
596C Purpose: obtain TKMO(NORBT,NORBT)
597C
598#include "implicit.h"
599      DIMENSION CMO(*), TKMO(NORBT,NORBT), WRK(*)
600C
601      CHARACTER*8 TKLABEL
602      PARAMETER (TKLABEL = 'KINETINT')
603C
604C Used from common blocks:
605C   INFORB : NORBT,NNBAST,...
606C   INFPRI : IPRAVE
607C
608#include "priunit.h"
609#include "inforb.h"
610#include "infpri.h"
611C
612      LOGICAL FOUND
613C
614      CALL QENTER('AVETK ')
615      KFRSAV = KFREE
616      CALL MEMGET('REAL',KTKAO,NNBAST,WRK,KFREE,LFREE)
617      FOUND = .FALSE.
618      CALL RDONEL(TKLABEL,FOUND,WRK(KTKAO),NNBAST)
619      IF (.NOT. FOUND) THEN
620         CALL QUIT('AVETK error: TK integrals not found')
621      END IF
622C
623C
624      CALL MEMGET('REAL',KTKPK,NNORBT,WRK,KFREE,LFREE)
625      CALL MEMGET('REAL',KSCR1,NNORBX,WRK,KFREE,LFREE)
626      DO 200 ISYM = 1,NSYM
627         NORBI = NORB(ISYM)
628      IF (NORBI.EQ.0) GO TO 200
629         NBASI = NBAS(ISYM)
630         JTKAO = KTKAO + IIBAS(ISYM)
631         JTKPK = KTKPK + IIORB(ISYM)
632         JCMO  = 1     + ICMO(ISYM)
633C
634         CALL UTHU(WRK(JTKAO),WRK(JTKPK),CMO(JCMO),WRK(KSCR1),
635     *             NBASI,NORBI)
636C        CALL UTHU(H,HT,U,WRK,NAO,NMO)
637  200 CONTINUE
638C
639      CALL PKSYM1(WRK(KSCR1),WRK(KTKPK),NORB,NSYM,-1)
640C     CALL PKSYM1(ASP,APK,NDIM,NBLK,IWAY)
641C
642      CALL DSPTSI(NORBT,WRK(KSCR1),TKMO)
643C     CALL DSPTSI(N,ASP,ASI)
644C
645      CALL MEMREL('AVETK',WRK,1,KFRSAV,KFREE,LFREE)
646C
647      IF ( IPRAVE .GE. 15 ) THEN
648         WRITE(LUPRI,'(3(/A))')
649     &      ' Matrix used for determining supersymmetries',
650     &      ' (the kinetic energy matrix)',
651     &      ' ==========================================='
652         CALL OUTPUT(TKMO,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
653      END IF
654C
655      CALL QEXIT('AVETK ')
656      RETURN
657      END
658C  /* Deck avedeg */
659      SUBROUTINE AVEDEG(TKMO,INDX1,INDX2,NDEGSS)
660C
661C   920630 hjaaj
662C   Purpose : find degeneracies, and construct
663C             INDX1(*) = labeling degenerate orbitals
664C             INDX2(*) = degeneracies of labels in INDX1
665C
666#include "implicit.h"
667      DIMENSION TKMO(NORBT,NORBT),INDX1(NORBT),INDX2(NORBT)
668      DIMENSION NDEGSS(*)
669C     DIMENSION NDEGSS(MXSSYM)
670C
671C Used from common blocks:
672C  INFINP : THRSSY
673C  INFORB : NORBT, NSSYM, NINFSS(MXSSYM,3),MXDGSS
674C  INFIND : ISSMO()
675C
676#include "maxash.h"
677#include "maxorb.h"
678#include "priunit.h"
679#include "infinp.h"
680#include "inforb.h"
681#include "infind.h"
682#include "infpri.h"
683C
684      CALL QENTER('AVEDEG')
685      CALL IZERO(INDX1,NORBT)
686      CALL IZERO(INDX2,NORBT)
687      IMOX = 0
688      NDEGMX = 0
689      DO 300 IMO = 1,NORBT
690      IF (INDX1(IMO) .EQ. 0) THEN
691         IMOX = IMOX + 1
692         INDX1(IMO) = IMOX
693         INDX2(IMOX) = 1
694         TKMOI = TKMO(IMO,IMO)
695         DO 200 JMO = IMO+1,NORBT
696            IF (ABS(TKMO(JMO,JMO)-TKMOI) .LT. THRSSY) THEN
697               INDX1(JMO)  = IMOX
698               INDX2(IMOX) = INDX2(IMOX) + 1
699            END IF
700  200    CONTINUE
701         NDEGMX = MAX(NDEGMX,INDX2(IMOX))
702      END IF
703  300 CONTINUE
704      MXDGSS = NDEGMX
705      NORXT = IMOX
706C
707C     Define NINFSS(i,1) = degeneracy of MO's in supsym_i
708C            NINFSS(i,2) = "root" supsym of this degeneracy
709C            NINFSS(i,3) = Abelian symmetry of supsym_i
710C            NDEGSS(i  ) = no. of deg. MO's in this supsym
711C
712      CALL IZERO(NINFSS,MXSSYM*3)
713      CALL IZERO(NDEGSS,MXSSYM)
714      IDGERR = 0
715      DO 500 IMO = 1,NORBT
716         ISSYM = ISSMO(IMO)
717         IMOX  = INDX1(IMO)
718         NDEGI = INDX2(IMOX)
719         IF (NINFSS(ISSYM,1) .EQ. 0) THEN
720            NINFSS(ISSYM,1) = NDEGI
721            NINFSS(ISSYM,2) = ISSYM
722            NINFSS(ISSYM,3) = ISMO(IMO)
723            NDEGSS(ISSYM)   = 1
724            DO 400 JMO = IMO+1,NORBT
725            IF (INDX1(JMO) .EQ. IMOX) THEN
726               JSSYM = ISSMO(JMO)
727               NINFSS(JSSYM,1) = NDEGI
728               NINFSS(JSSYM,2) = ISSYM
729               NINFSS(JSSYM,3) = ISMO(JMO)
730               NDEGSS(JSSYM)   = NDEGSS(JSSYM) + 1
731            END IF
732  400       CONTINUE
733         ELSE IF (NINFSS(ISSYM,1) .NE. NDEGI) THEN
734            WRITE(LUPRI,'(A,4I4)')
735     &           ' Error for IMO,ISSYM, NINFSS(ISSYM,1),NDEGI  = ',
736     &           IMO,ISSYM, NINFSS(ISSYM,1),NDEGI
737            IDGERR = IDGERR + 1
738         END IF
739  500 CONTINUE
740C
741C Print NINFSS for debugging
742C
743      IF ( IPRAVE .GE. 3 .OR. IDGERR .GT. 0) THEN
744         IF (IDGERR .GT. 0) THEN
745            WRITE(LUPRI,'(7(/A))')
746     &      ' AVEDEG: Degeneracy error, dump of NINFSS follows below.',
747     &      ' Probable cause:'//
748     &         ' nuclear geometry has only approximate super symmetry.',
749     &      ' Some possible actions:',
750     &' (1) include more digits in nuclear coordinates',
751     &' (2) attempt to force super symmetry by making .THRSSY larger',
752     &' (3) attempt to avoid degeneracy error by making .THRSSY smaller'
753     &,' (4) disable supersymmetry by removing .SUPSYM'
754         END IF
755         WRITE(LUPRI,'(4(/A))')
756     &   ' NINFSS(i,1) = degeneracy of MOs in supsym_i',
757     &   ' NINFSS(i,2) = root supsym of this degeneracy',
758     &   ' NINFSS(i,3) = Abelian symmetry of supsym_i',
759     &   ' NDEGSS(i)   = no. of deg. supsyms currently in this supsym'
760         WRITE(LUPRI,'(/A/A)')
761     &      '  ISS         NINFSS 1 to 3            NDEGSS',
762     &      ' ============================================'
763         DO 710 ISSYM = 1, NSSYM
764            WRITE(LUPRI,'(1X,I4,4I10)')
765     &         ISSYM,(NINFSS(ISSYM,J),J=1,3),NDEGSS(ISSYM)
766  710    CONTINUE
767         WRITE(LUPRI,'(A/)')
768     &      ' ============================================'
769      END IF
770C
771C     print for debuggers
772C
773      IF ( IPRAVE .GE. 10 .OR. IDGERR .GT. 0 ) THEN
774         WRITE(LUPRI,'(A,F15.12/)')
775     &      ' Degeneracy pattern found with THRSSY =',THRSSY
776         WRITE(LUPRI,'("    I  ISSMO INDX1 INDX2       TKMO(I,I)")')
777       WRITE(LUPRI,'(" ============================================")')
778         DO 700 IMO = 1, NORBT
779            WRITE(LUPRI,'(I5,3I6,F25.15)')
780     &         IMO, ISSMO(IMO), INDX1(IMO), INDX2(IMO), TKMO(IMO,IMO)
781  700    CONTINUE
782       WRITE(LUPRI,'(" ============================================")')
783      END IF
784      IF (IDGERR .GT. 0) THEN
785         CALL QUIT('AVEDEG: Degeneracy error ...')
786      END IF
787C
788      CALL QEXIT('AVEDEG')
789      RETURN
790      END
791C  /* Deck avedis */
792      SUBROUTINE AVEDIS(ISSYM,ISYM,NDEGI,NORBI,NBASI,NSSMOI,NSSMON,
793     *                  CMOI,TKMO,INDX1,UMAT,IMODEG,WRK)
794C
795#include "implicit.h"
796C
797      DIMENSION CMOI(NBASI,NORBI), TKMO(NORBT,NORBT)
798      DIMENSION INDX1(NORBT),UMAT(NDEGI,NDEGI)
799      DIMENSION IMODEG(NDEGI,NSSMON), WRK(NBASI,NDEGI)
800C
801C Used from common blocks
802C  INFORB: NSSYM, NORBT, NORBSS(), IORBSS()
803C  INFIND: ISMO(),ISSORD()
804C  INFPRI: IPRAVE
805C
806#include "maxash.h"
807#include "maxorb.h"
808#include "priunit.h"
809#include "inforb.h"
810#include "infind.h"
811#include "infpri.h"
812C
813      CALL QENTER('AVEDIS')
814      IF (IPRAVE .GE. 3) THEN
815         WRITE(LUPRI,'(/A,I5/A,I5)')
816     *      ' AVEDIS disentangling mixed degenerate components in'//
817     *      ' supersymmetry:',ISSYM,' Abelian symmetry:',ISYM
818      END IF
819      IMOXP = 0
820      KMOSS = 0
821      IMOSSE = IORBSS(ISSYM)+NORBSS(ISSYM)
822      DO 180 IMOSS = IORBSS(ISSYM)+1,IMOSSE
823         IMO = ISSORD(IMOSS)
824         IMOX = INDX1(IMO)
825         IF (IMOX .GT. IMOXP) THEN
826            IMOXP = IMOX
827            KMOSS = KMOSS + 1
828            IDEG = 1
829            IMODEG(1,KMOSS) = IMO
830            DO 120 JMOSS = IMOSS+1,IMOSSE
831               JMO = ISSORD(JMOSS)
832               IF (INDX1(JMO) .EQ. IMOX) THEN
833                  IDEG = IDEG + 1
834                  IMODEG(IDEG,KMOSS) = JMO
835               END IF
836  120       CONTINUE
837         END IF
838  180 CONTINUE
839      IF (KMOSS .NE. NSSMON) THEN
840         CALL QUIT('AVEDIS error: KMOSS .ne. NSSMON')
841      END IF
842C
843      IOFFMO = IORB(ISYM)
844C     check that first block is diagonal
845      CALL AVEGTU(NDEGI,NORBT,1,1,
846     *            IMODEG(1,1),IMODEG(1,1),UMAT,TKMO,WRK)
847      DO 280 KMOSS = 2,NSSMON
848C        a) check that KMOSS block is diagonal
849         CALL AVEGTU(NDEGI,NORBT,KMOSS,KMOSS,
850     *               IMODEG(1,KMOSS),IMODEG(1,KMOSS),UMAT,TKMO,WRK)
851C        b) find UMAT to transform the KMOSS degenerate orbital set to
852C           the same irrep rows as the first deg. orbital set
853         CALL AVEGTU(NDEGI,NORBT,1,KMOSS,
854     *               IMODEG(1,1),IMODEG(1,KMOSS),UMAT,TKMO,WRK)
855         CALL DZERO(WRK,NDEGI*NBASI)
856C        c) do the transformation
857         DO 240 IDEG = 1,NDEGI
858            IMO = IMODEG(IDEG,KMOSS) - IOFFMO
859            DO 220 JDEG = 1,NDEGI
860               CALL DAXPY(NBASI,UMAT(JDEG,IDEG),
861     *                    CMOI(1,IMO),1,WRK(1,JDEG),1)
862  220       CONTINUE
863  240    CONTINUE
864         DO 260 IDEG = 1,NDEGI
865            IMO = IMODEG(IDEG,KMOSS) - IOFFMO
866            CALL DCOPY(NBASI,WRK(1,IDEG),1,CMOI(1,IMO),1)
867  260    CONTINUE
868  280 CONTINUE
869C
870      CALL QEXIT('AVEDIS')
871      RETURN
872      END
873C  /* Deck avegtu */
874      SUBROUTINE AVEGTU(NDEGI,NORBT,IMOSS,KMOSS,IMOI,IMOK,
875     *                  UMAT,TKMO,UTUMAT)
876C
877#include "implicit.h"
878      DIMENSION IMOI(NDEGI), IMOK(NDEGI)
879      DIMENSION UMAT(NDEGI,NDEGI), TKMO(NORBT,NORBT)
880      DIMENSION UTUMAT(NDEGI,NDEGI)
881      PARAMETER (D1 = 1.0D0)
882C
883C Used from common blocks:
884C   INFINP : THRSSY
885C   INFPRI : IPRAVE
886C
887#include "maxorb.h"
888#include "priunit.h"
889#include "infinp.h"
890#include "infpri.h"
891C
892      CALL QENTER('AVEGTU')
893      DO 120 K = 1,NDEGI
894         DO 110 I = 1,NDEGI
895            UMAT(I,K) = TKMO( IMOI(I) , IMOK(K) )
896  110    CONTINUE
897  120 CONTINUE
898C
899C  -- Do some printout for debuggers
900C
901      IF (IPRAVE .GE. 10) THEN
902         WRITE(LUPRI,'(/A/A/A,3I5)')
903     *      ' Test output from AVEGTU',
904     *      ' =======================',
905     *      ' IMOSS, KMOSS, NDEGI =',IMOSS,KMOSS,NDEGI
906         WRITE(LUPRI,'(A,20I4)') ' IMOI :',(IMOI(I),I=1,NDEGI)
907         WRITE(LUPRI,'(A,20I4)') ' IMOK :',(IMOK(I),I=1,NDEGI)
908         WRITE(LUPRI,'(/A/A/)') ' UMAT found in AVEGTU',
909     *                         ' ===================='
910         CALL OUTPUT(UMAT,1,NDEGI,1,NDEGI,NDEGI,NDEGI,1,LUPRI)
911      END IF
912C
913      IF (IMOSS .EQ. KMOSS) THEN
914C     ... check that UMAT is diagonal
915         NERR = NOFDIA(NDEGI,NDEGI,UMAT,THRSSY)
916         IF (NERR .GT. 0) THEN
917            WRITE(LUPRI,'(2A,I5)') ' AVEGTU error, ',
918     *         'no. of non-zero off-diagonal elements in U =',NERR
919            WRITE(LUPRI,'(A,I4)') ' -- IMOSS = ', IMOSS
920            WRITE(LUPRI,'(A,I4)') ' -- KMOSS = ', KMOSS
921            WRITE(LUPRI,'(/A)') ' The U matrix :'
922            CALL OUTPUT(UMAT,1,NDEGI,1,NDEGI,NDEGI,NDEGI,1,LUPRI)
923            CALL QUIT('AVEGTU error for IMOSS .eq. KMOSS')
924         END IF
925      ELSE
926C     ... IMOSS .ne. KMOSS
927C         a) check that matrix is x*orthogonal matrix
928         ITURN = 0
929  300    ITURN = ITURN + 1
930         CALL DGEMM('T','N',NDEGI,NDEGI,NDEGI,1.D0,
931     &              UMAT,NDEGI,
932     &              UMAT,NDEGI,0.D0,
933     &              UTUMAT,NDEGI)
934         NERR = NOFDIA(NDEGI,NDEGI,UTUMAT,THRSSY)
935         IF (NERR .GT. 0 .OR. IPRAVE .GE. 10) THEN
936            IF (NERR .GT. 0) WRITE(LUPRI,'(2A,I5)') ' AVEGTU error, ',
937     *         'no. of non-zero off-diagonal elements in U^(T)U =',NERR
938            WRITE(LUPRI,'(A,I4)') ' -- IMOSS = ', IMOSS
939            WRITE(LUPRI,'(A,I4)') ' -- KMOSS = ', KMOSS
940            WRITE(LUPRI,'(/A)') ' The U matrix :'
941            CALL OUTPUT(UMAT,1,NDEGI,1,NDEGI,NDEGI,NDEGI,1,LUPRI)
942            WRITE(LUPRI,'(/A)') ' The UTU matrix :'
943            CALL OUTPUT(UTUMAT,1,NDEGI,1,NDEGI,NDEGI,NDEGI,1,LUPRI)
944            IF (NERR .GT. 0)
945     &           CALL QUIT('AVEGTU error for IMOSS .ne. KMOSS')
946         END IF
947C
948C HJMAERKE 940818-hjaaj: we ought to force UMAT orthogonal with
949C   symmetric orthonormalization; with the current code we may
950C   get deviations from orthonormality up to THRSSY
951         X2 = UTUMAT(1,1)
952         DO 310 I = 2,NDEGI
953            X2 = X2 * UTUMAT(I,I)
954  310    CONTINUE
955C        det(U) = sqrt(X2) = X2 ** ( .5 )
956         XEXP = 2*NDEGI
957         XEXP = - D1 / XEXP
958         XINV = X2 ** ( XEXP )
959         CALL DSCAL(NDEGI*NDEGI,XINV,UMAT,1)
960         IF (ABS(D1-XINV) .GT. THRSSY) THEN
961            IF (IPRAVE .GE. 10 .OR. ITURN .GE. 3) THEN
962               WRITE(LUPRI,'(A,I4)') ' -- ITURN = ', ITURN
963               WRITE(LUPRI,'(A,1D20.10)') ' -- XINV  = ', XINV
964               WRITE(LUPRI,'(A,I4)') ' -- IMOSS = ', IMOSS
965               WRITE(LUPRI,'(A,I4)') ' -- KMOSS = ', KMOSS
966               WRITE(LUPRI,'(/A)') ' The U matrix :'
967               CALL OUTPUT(UMAT,1,NDEGI,1,NDEGI,NDEGI,NDEGI,1,LUPRI)
968            END IF
969            IF (ITURN .LT. 3) GO TO 300
970            WRITE(LUPRI,'(A)')
971     *      ' AVEGTU error : failed to normalize UMAT'
972            CALL QUIT('AVEGTU error: could not normalize UMAT')
973         END IF
974      END IF
975      CALL QEXIT('AVEGTU')
976      RETURN
977      END
978C  /* Deck avepha */
979      SUBROUTINE AVEPHA(ISSYMR,CMO,TKMO,IPHCHA)
980C
981C  2-Jul-1992 hjaaj+hh
982C  Purpose: make relative phases the same in each block
983C           of a set of degenerate mo.s.
984C           A block corresponds to orbitals following
985C           a row of an irrep.
986C
987#include "implicit.h"
988      DIMENSION CMO(NCMOT), TKMO(NORBT,NORBT)
989C
990C Used from common blocks:
991C   INFINP : THRSSY
992C   INFORB : NCMOT,NORBT,NORBSS(), IORBSS(),...,NINFSS(*,3),...
993C   INFIND : ISSORD()
994C
995#include "maxash.h"
996#include "maxorb.h"
997#include "priunit.h"
998#include "infinp.h"
999#include "inforb.h"
1000#include "infind.h"
1001#include "infpri.h"
1002C
1003      CALL QENTER('AVEPHA')
1004      IF (NINFSS(ISSYMR,2) .NE. ISSYMR) THEN
1005         WRITE(LUPRI,'(/A)')
1006     *      ' AVEPHA error: ISSYMR is not a "root" symmetry'
1007         WRITE(LUPRI,'(A,I4)') ' ISSYMR           = ', ISSYMR
1008         WRITE(LUPRI,'(A,I4)') ' NINFSS(ISSYMR,2) = ', NINFSS(ISSYMR,2)
1009         CALL QUIT('AVEPHA error: ISSYMR is not a "root" symmetry')
1010      END IF
1011      IOSSI = IORBSS(ISSYMR)
1012      NOSSI = NORBSS(ISSYMR)
1013      DO 800 JSSYM = 1,NSSYM
1014         JSSYMR = NINFSS(JSSYM,2)
1015         IF (JSSYMR .EQ. ISSYMR) THEN
1016            IF (JSSYM .LT. ISSYMR) THEN
1017               CALL QUIT('AVEPHA error: JSSYM .lt ISSYMR')
1018            END IF
1019            IF (NORBSS(JSSYM) .NE. NOSSI) THEN
1020               WRITE(LUPRI,'(/A/A,2I5/A)')
1021     *            ' AVEHPA error: Diff. no. of orb.s in deg. supsym.s',
1022     *            ' JSSYMR, JSSYM =',JSSYMR,JSSYM,
1023     *            ' Dump of super symmetry information :'
1024               CALL AVEDMP(LUPRI)
1025               CALL QUIT(
1026     &         'AVEHPA error: Diff. no. of orb.s in deg. supsym.s')
1027            END IF
1028            ISYMJ = NINFSS(JSSYM,3)
1029            ICMOJ = ICMO(ISYMJ)
1030            NBASJ = NBAS(ISYMJ)
1031            IORBJ = IORB(ISYMJ)
1032            IOSSJ = IORBSS(JSSYM)
1033            IMO1  = ISSORD(IOSSI + 1)
1034            JMO1  = ISSORD(IOSSJ + 1)
1035            DO 700 K = 1,NOSSI
1036               IMO = ISSORD(IOSSI + K)
1037               JMO = ISSORD(IOSSJ + K)
1038               IF (ABS(TKMO(IMO,IMO) - TKMO(JMO,JMO)) .GT. THRSSY) THEN
1039                  WRITE (LUPRI,'(/A,2I5)')
1040     &               ' AVEPHA error: different ordering'//
1041     &               ' of degenerate orbitals in supsyms:',ISSYMR,JSSYM
1042                  WRITE(LUPRI,'(/A,2I5/A,D20.10/A,D20.10)')
1043     &            ' IMO, JMO       =',IMO,JMO,
1044     &            ' TKMO(IMO,IMO)  =',TKMO(IMO,IMO),
1045     &            ' TKMO(JMO,JMO)  =',TKMO(JMO,JMO)
1046                  CALL QUIT(
1047     &            'AVEPHA error: different ordering of deg. orb.s')
1048               END IF
1049               IF      (ABS(TKMO(IMO,IMO1)) .LT. THRSSY) THEN
1050                  WRITE(LUPRI,'(/A/A,2I5/A,D10.2)')
1051     &            ' AVEPHA error: TKMO(IMO,IMO1) is zero',
1052     &            ' IMO, IMO1      =',IMO,IMO1,
1053     &            ' TKMO(IMO,IMO1) =',TKMO(IMO,IMO1)
1054                  CALL QTRACE(LUERR)
1055                  CALL QUIT('AVEPHA error: TKMO(IMO,IMO1) is zero')
1056               ELSE IF (ABS(TKMO(IMO,IMO1)+TKMO(JMO,JMO1)) .LT. THRSSY)
1057     *         THEN
1058                  IPHCHA = IPHCHA + 1
1059                  IF (IPRAVE .GE. 4) THEN
1060                     WRITE(LUPRI,'(A,I4)')
1061     *                  ' AVEPHA changing phase of mo no.',JMO
1062                  END IF
1063                  JOFF = ICMOJ + NBASJ*(JMO-IORBJ-1)
1064                  DO 600 J = 1,NBASJ
1065                     CMO(JOFF+J) = -CMO(JOFF+J)
1066  600             CONTINUE
1067               ELSE IF (ABS(TKMO(IMO,IMO1)-TKMO(JMO,JMO1)) .GT. THRSSY)
1068     *         THEN
1069                  WRITE(LUPRI,'(/A,2(/A,2I5/A,D20.10))')
1070     &            ' AVEPHA error: TKMO(IMO,IMO1) not TKMO(JMO,JMO1)',
1071     &            ' IMO, IMO1      =',IMO,IMO1,
1072     &            ' TKMO(IMO,IMO1) =',TKMO(IMO,IMO1),
1073     &            ' JMO, JMO1      =',JMO,JMO1,
1074     &            ' TKMO(JMO,JMO1) =',TKMO(JMO,JMO1)
1075                  CALL QTRACE(LUERR)
1076                  CALL QUIT('AVEPHA error: error in deg. orb.s code 2')
1077               END IF
1078  700       CONTINUE
1079         END IF
1080  800 CONTINUE
1081      CALL QEXIT('AVEPHA')
1082      RETURN
1083      END
1084#if defined (VAR_NOTUSED  )
1085 AVEUDV is not used
1086C  /* Deck aveudv */
1087      SUBROUTINE AVEUDV(UDV)
1088C
1089C     09-Jul-1992 Hinne Hettema
1090C
1091C     Purpose : perform an averaging of the one-electron density
1092C     matrix. First we identify an IMO which is degenerate with
1093C     other MO's, then we average the matrix rows and colums.
1094C
1095#include "implicit.h"
1096      DIMENSION UDV(NASHT,*)
1097C
1098C     Used from common blocks
1099C     INFORB    : NSYM ,NSSYM, NORBSS(*), IORBSS(*), NINFSS(*,3)
1100C     INFIND    : ICH(*), ISSMO(*), ISSORD(*)
1101C
1102#include "maxorb.h"
1103#include "maxash.h"
1104#include "inforb.h"
1105#include "infind.h"
1106#include "infpri.h"
1107C
1108      PARAMETER ( MAXDEG = 20, D0 = 0.0D0, D1 = 1.0D0 )
1109      DIMENSION ISSOFF(MAXDEG),IDEGV(MAXDEG,2)
1110C
1111      CALL QENTER('AVEUDV')
1112C
1113      DO 800 ISSYMR = 1, NSSYM
1114C     -- skip if supersym not degenerate or not 'root' supersym
1115         IF ( NINFSS(ISSYMR,1) .EQ. 1) GO TO 800
1116         IF ( NINFSS(ISSYMR,2) .NE. ISSYMR) GO TO 800
1117C
1118         NDEG  = 0
1119         DO 100 JSSYM = ISSYMR, NSSYM
1120C        -- skip if not same 'root' supersymmetry
1121C           (because then not degenerate at all)
1122            IF ( NINFSS(JSSYM,2) .EQ. ISSYMR) THEN
1123               NDEG = NDEG + 1
1124               ISSOFF(NDEG) = IORBSS(JSSYM) - IORBSS(ISSYMR)
1125            END IF
1126  100    CONTINUE
1127         IF ( NDEG .NE. NINFSS(ISSYMR,1) ) THEN
1128            WRITE(LUERR,'(A)') ' AVEUDV error : inconsistent'//
1129     &            ' degeneracies'
1130            WRITE(LUERR,'(A,I4)') ' - ISSYMR           = ', ISSYMR
1131            WRITE(LUERR,'(A,I4)') ' - NDEG             = ', NDEG
1132            WRITE(LUERR,'(A,I4)') ' - NINFSS(ISSYMR,1) = ',
1133     &                 NINFSS(ISSYMR,1)
1134            CALL QUIT(' AVEUDV : inconsistent degeneracy in NINFSS')
1135         END IF
1136C
1137C     -- Compute averaging factor
1138         AVFAC = NDEG
1139         AVFAC = D1 / AVFAC
1140C
1141C     -- in the following, use supersymmetry info to get IMO
1142         IRSYM  = NINFSS(ISSYMR,3)
1143         NASHI  = NASH(IRSYM)
1144         IOFFSR = IORBSS(ISSYMR)
1145         IOFFDV = IASH(IRSYM) + 1
1146         DO 700 IMOSS = IOFFSR+1, IOFFSR+NORBSS(ISSYMR)
1147            IMO  = ISSORD(IMOSS)
1148            IMOW = ICH(IMO)
1149            IF (IMOW .LE. 0) GO TO 700
1150            DO 600 JMOSS = IMOSS, IOFFSR+NORBSS(ISSYMR)
1151               JMO = ISSORD(JMOSS)
1152               JMOW = ICH(JMO)
1153               IF ( JMOW .LE. 0 ) GO TO 600
1154               DO 300 IDEG = 2, NDEG
1155                  KMO = ISSORD(IMOSS + ISSOFF(IDEG))
1156                  KMOW = ICH(KMO)
1157                  IF (KMOW .LE. 0 ) THEN
1158                     CALL QUIT('AVEUDV code 300.1: ICH(KMO).le.0')
1159                  END IF
1160                  LMO = ISSORD(JMOSS + ISSOFF(IDEG))
1161                  LMOW = ICH(LMO)
1162                  IF (LMOW .LE. 0 ) THEN
1163                     CALL QUIT('AVEUDV code 300.2: ICH(LMO).le.0')
1164                  END IF
1165                  IDEGV(IDEG,1) = KMOW
1166                  IDEGV(IDEG,2) = LMOW
1167  300          CONTINUE
1168C           -- now collect all elements which should be made equal
1169               XUDVIJ = UDV(IMOW,JMOW)
1170               DO 400 IDEG = 2, NDEG
1171                  KMOW = IDEGV(IDEG,1)
1172                  LMOW = IDEGV(IDEG,2)
1173                  XUDVIJ = XUDVIJ + UDV(KMOW,LMOW)
1174  400          CONTINUE
1175               XUDVIJ = XUDVIJ * AVFAC
1176C           -- now distribute all elements which should be made equal
1177               UDV(IMOW,JMOW) = XUDVIJ
1178               IF (IMOW .NE. JMOW) UDV(JMOW,IMOW) = XUDVIJ
1179               DO 500 IDEG = 2, NDEG
1180                  KMOW = IDEGV(IDEG,1)
1181                  LMOW = IDEGV(IDEG,2)
1182                  UDV(KMOW,LMOW) = XUDVIJ
1183                  IF(KMOW .NE. LMOW) UDV(LMOW,KMOW) = XUDVIJ
1184  500          CONTINUE
1185  600       CONTINUE
1186  700    CONTINUE
1187C
1188  800 CONTINUE
1189C
1190C     Zero elements which are not totally symmetric
1191C     with respect to the supersymmetry.
1192C
1193      DO 980 ISYM = 1,NSYM
1194         IASHI = IASH(ISYM)
1195         NASHI = NASH(ISYM)
1196         DO 960 IMOW = IASHI+1,IASHI+NASHI
1197            IMO = ISX(NISHT+IMOW)
1198            DO 950 JMOW = IMOW+1,IASHI+NASHI
1199               JMO = ISX(NISHT+JMOW)
1200               IF (ISSMO(JMO) .NE. ISSMO(IMO)) THEN
1201                  UDV(IMOW,JMOW) = D0
1202                  UDV(JMOW,IMOW) = D0
1203               END IF
1204  950       CONTINUE
1205  960    CONTINUE
1206  980 CONTINUE
1207C
1208      CALL QEXIT('AVEUDV')
1209C
1210      RETURN
1211      END
1212C  /* Deck aveh1 */
1213 AVEH1 is not used
1214      SUBROUTINE AVEH1(H1)
1215C
1216C     12-Jul-1992 Hinne Hettema
1217C
1218C     Purpose : perform an averaging of a one-electron
1219C     matrix. First we identify an IMO which is degenerate with
1220C     other MO's, then we average the matrix rows and colums.
1221C
1222#include "implicit.h"
1223      DIMENSION H1(NORBT,*)
1224C
1225C     Used from common blocks
1226C     INFORB    : NSYM ,NSSYM, NORBSS(*), IORBSS(*), NINFSS(*,3)
1227C     INFIND    : ICH(*), ISSMO(*), ISSORD(*)
1228C
1229#include "maxorb.h"
1230#include "maxash.h"
1231#include "priunit.h"
1232#include "inforb.h"
1233#include "infind.h"
1234#include "infpri.h"
1235C
1236      PARAMETER ( MAXDEG = 20, D0 = 0.0D0, D1 = 1.0D0 )
1237      DIMENSION ISSOFF(MAXDEG),IDEGV(MAXDEG,2)
1238C
1239      CALL QENTER('AVEH1 ')
1240C
1241      IF ( IPRAVE .GE. 10 ) THEN
1242         WRITE(LUPRI,'(//A)') ' H1 matrix before averaging'
1243         WRITE(LUPRI,'(A)')   ' =========================='
1244         CALL OUTPUT(H1,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
1245      END IF
1246C
1247      DO 800 ISSYMR = 1, NSSYM
1248C     -- skip if supersym not degenerate or not 'root' supersym
1249         IF ( NINFSS(ISSYMR,1) .EQ. 1) GO TO 800
1250         IF ( NINFSS(ISSYMR,2) .NE. ISSYMR) GO TO 800
1251C
1252         NDEG  = 0
1253         DO 100 JSSYM = ISSYMR, NSSYM
1254C        -- skip if not same 'root' supersymmetry
1255C           (because then not degenerate at all)
1256            IF ( NINFSS(JSSYM,2) .EQ. ISSYMR) THEN
1257               NDEG = NDEG + 1
1258               ISSOFF(NDEG) = IORBSS(JSSYM) - IORBSS(ISSYMR)
1259            END IF
1260  100    CONTINUE
1261         IF ( NDEG .NE. NINFSS(ISSYMR,1) ) THEN
1262            WRITE(LUERR,'(A)') ' AVEH1 error : inconsistent'//
1263     &            ' degeneracies'
1264            WRITE(LUERR,'(A,I4)') ' - ISSYMR           = ', ISSYMR
1265            WRITE(LUERR,'(A,I4)') ' - NDEG             = ', NDEG
1266            WRITE(LUERR,'(A,I4)') ' - NINFSS(ISSYMR,1) = ',
1267     &                 NINFSS(ISSYMR,1)
1268            CALL QUIT(' AVEH1 : inconsistent degeneracy in NINFSS')
1269         END IF
1270C
1271C     -- Compute averaging factor
1272         AVFAC = NDEG
1273         AVFAC = D1 / AVFAC
1274C
1275C     -- in the following, use supersymmetry info to get IMO
1276         IRSYM  = NINFSS(ISSYMR,3)
1277         NORBI  = NORB(IRSYM)
1278         IOFFSR = IORBSS(ISSYMR)
1279         DO 700 IMOSS = IOFFSR+1, IOFFSR+NORBSS(ISSYMR)
1280            IMO  = ISSORD(IMOSS)
1281            DO 600 JMOSS = IMOSS, IOFFSR+NORBSS(ISSYMR)
1282               JMO = ISSORD(JMOSS)
1283               IDEGV(1,1) = IMO
1284               IDEGV(1,2) = JMO
1285               DO 300 IDEG = 2, NDEG
1286                  IDEGV(IDEG,1) = ISSORD(IMOSS + ISSOFF(IDEG))
1287                  IDEGV(IDEG,2) = ISSORD(JMOSS + ISSOFF(IDEG))
1288  300          CONTINUE
1289C           -- now collect all elements which should be made equal
1290               XH1IJ = H1(IMO,JMO)
1291               DO 400 IDEG = 2, NDEG
1292                  KMO = IDEGV(IDEG,1)
1293                  LMO = IDEGV(IDEG,2)
1294                  XH1IJ = XH1IJ + H1(KMO,LMO)
1295  400          CONTINUE
1296               XH1IJ = XH1IJ * AVFAC
1297C           -- now distribute all elements which should be made equal
1298               H1(IMO,JMO) = XH1IJ
1299               IF (IMO .NE. JMO) H1(JMO,IMO) = XH1IJ
1300               DO 500 IDEG = 2, NDEG
1301                  KMO = IDEGV(IDEG,1)
1302                  LMO = IDEGV(IDEG,2)
1303                  H1(KMO,LMO) = XH1IJ
1304                  IF(KMO .NE. LMO) H1(LMO,KMO) = XH1IJ
1305  500          CONTINUE
1306  600       CONTINUE
1307  700    CONTINUE
1308C
1309  800 CONTINUE
1310C
1311C     Zero elements which are not totally symmetric
1312C     with respect to the supersymmetry.
1313C
1314      DO 980 ISYM = 1,NSYM
1315         IORBI = IORB(ISYM)
1316         NORBI = NORB(ISYM)
1317         DO 960 IMO = IORBI+1,IORBI+NORBI
1318            DO 950 JMO = IMO+1,IORBI+NORBI
1319               IF (ISSMO(JMO) .NE. ISSMO(IMO)) THEN
1320                  H1(IMO,JMO) = D0
1321                  H1(JMO,IMO) = D0
1322               END IF
1323  950       CONTINUE
1324  960    CONTINUE
1325  980 CONTINUE
1326C
1327      IF ( IPRAVE .GE. 10 ) THEN
1328         WRITE(LUPRI,'(//A)') ' H1 matrix after averaging'
1329         WRITE(LUPRI,'(A)')   ' ========================='
1330         CALL OUTPUT(H1,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
1331      END IF
1332C
1333      CALL QEXIT('AVEH1 ')
1334C
1335      RETURN
1336      END
1337#endif
1338C  /* Deck avedv */
1339      SUBROUTINE AVEDV(DV)
1340C
1341C     09-Jul-1992 Hinne Hettema
1342C
1343C     Purpose : perform an averaging of the one-electron density
1344C     matrix. First we identify an IMO which is degenerate with
1345C     other MO's, then we average the matrix rows and colums.
1346C
1347#include "implicit.h"
1348      DIMENSION DV(NNASHX)
1349C
1350C     Used from common blocks
1351C     INFORB    : NSYM ,NSSYM, NORBSS(*), IORBSS(*), NINFSS(*,3)
1352C     INFIND    : IROW(*), ICH(*), ISSMO(*), ISSORD(*)
1353C
1354#include "maxorb.h"
1355#include "maxash.h"
1356#include "priunit.h"
1357#include "inforb.h"
1358#include "infind.h"
1359#include "infpri.h"
1360C
1361      PARAMETER ( MAXDEG = 20, D0 = 0.0D0, D1 = 1.0D0 )
1362      DIMENSION ISSOFF(MAXDEG),IDEGV(MAXDEG)
1363C
1364      CALL QENTER('AVEDV ')
1365C
1366      IF (IPRAVE .GE. 5) THEN
1367         WRITE(LUPRI,'(/A)')
1368     *      ' AVEDV test output: DV matrix before averaging'
1369         CALL OUTPAK(DV,NASHT,1,LUPRI)
1370      END IF
1371C
1372      DO 800 ISSYMR = 1, NSSYM
1373C     -- skip if supersym not degenerate or not 'root' supersym
1374         IF ( NINFSS(ISSYMR,1) .EQ. 1) GO TO 800
1375         IF ( NINFSS(ISSYMR,2) .NE. ISSYMR) GO TO 800
1376C
1377         NDEG  = 0
1378         DO 100 JSSYM = ISSYMR, NSSYM
1379C        -- skip if not same 'root' supersymmetry
1380C           (because then not degenerate at all)
1381            IF ( NINFSS(JSSYM,2) .EQ. ISSYMR) THEN
1382               NDEG = NDEG + 1
1383               ISSOFF(NDEG) = IORBSS(JSSYM) - IORBSS(ISSYMR)
1384            END IF
1385  100    CONTINUE
1386         IF ( NDEG .NE. NINFSS(ISSYMR,1) ) THEN
1387            WRITE(LUERR,'(A)') ' AVEDV error : inconsistent'//
1388     &            ' degeneracies'
1389            WRITE(LUERR,'(A,I4)') ' - ISSYMR           = ', ISSYMR
1390            WRITE(LUERR,'(A,I4)') ' - NDEG             = ', NDEG
1391            WRITE(LUERR,'(A,I4)') ' - NINFSS(ISSYMR,1) = ',
1392     &                 NINFSS(ISSYMR,1)
1393            CALL QUIT(' AVEDV : inconsistent degeneracy in NINFSS')
1394         END IF
1395C
1396C     -- Compute averaging factor
1397         AVFAC = NDEG
1398         AVFAC = D1 / AVFAC
1399C
1400C     -- in the following, use supersymmetry info to get IMO
1401         IRSYM  = NINFSS(ISSYMR,3)
1402         NASHI  = NASH(IRSYM)
1403         IOFFSR = IORBSS(ISSYMR)
1404         IOFFDV = IASH(IRSYM) + 1
1405         DO 700 IMOSS = IOFFSR+1, IOFFSR+NORBSS(ISSYMR)
1406            IMO  = ISSORD(IMOSS)
1407            IMOW = ICH(IMO)
1408            IF (IMOW .LE. 0) GO TO 700
1409            DO 600 JMOSS = IOFFSR+1, IMOSS
1410               JMO = ISSORD(JMOSS)
1411               JMOW = ICH(JMO)
1412               IF ( JMOW .LE. 0 ) GO TO 600
1413               IDEGV(1) = IROW(IMOW) + JMOW
1414               DO 300 IDEG = 2, NDEG
1415                  KMO = ISSORD(IMOSS + ISSOFF(IDEG))
1416                  KMOW = ICH(KMO)
1417                  IF (KMOW .LE. 0 ) THEN
1418                     CALL QUIT('AVEDV code 300.1')
1419                  END IF
1420                  LMO = ISSORD(JMOSS + ISSOFF(IDEG))
1421                  LMOW = ICH(LMO)
1422                  IF (LMOW .LE. 0 ) THEN
1423                     CALL QUIT('AVEDV code 300.2')
1424                  END IF
1425                  IDEGV(IDEG) = IROW(KMOW) + LMOW
1426  300          CONTINUE
1427C           -- now collect all elements which should be made equal
1428               XDVIJ = DV(IDEGV(1))
1429               DO 400 IDEG = 2, NDEG
1430                  XDVIJ = XDVIJ + DV(IDEGV(IDEG))
1431  400          CONTINUE
1432               XDVIJ = XDVIJ * AVFAC
1433C           -- now distribute all elements which should be made equal
1434               DO 500 IDEG = 1, NDEG
1435                  DV(IDEGV(IDEG)) = XDVIJ
1436  500          CONTINUE
1437  600       CONTINUE
1438  700    CONTINUE
1439C
1440  800 CONTINUE
1441C
1442C     Zero elements which are not totally symmetric
1443C     with respect to the supersymmetry.
1444C
1445      DO 980 ISYM = 1,NSYM
1446         IASHI = IASH(ISYM)
1447         NASHI = NASH(ISYM)
1448         DO 960 IMOW = IASHI+2,IASHI+NASHI
1449            IMO = ISX(NISHT+IMOW)
1450            DO 950 JMOW = IASHI+1,IMOW-1
1451               JMO = ISX(NISHT+JMOW)
1452               IF (ISSMO(JMO) .NE. ISSMO(IMO)) THEN
1453                  IDV = IROW(IMOW) + JMOW
1454                  DV(IDV) = D0
1455               END IF
1456  950       CONTINUE
1457  960    CONTINUE
1458  980 CONTINUE
1459C
1460      IF (IPRAVE .GE. 5) THEN
1461         WRITE(LUPRI,'(/A)')
1462     *      ' AVEDV test output: DV matrix after averaging'
1463         CALL OUTPAK(DV,NASHT,1,LUPRI)
1464      END IF
1465      CALL QEXIT('AVEDV ')
1466C
1467      RETURN
1468      END
1469#if defined (VAR_NOTUSED  )
1470C  /* Deck aveh1p */
1471      SUBROUTINE AVEH1P(H1P)
1472C
1473C     12-Jul-1992 Hinne Hettema
1474C
1475C     Purpose : perform an averaging of a one-electron
1476C     matrix. First we identify an IMO which is degenerate with
1477C     other MO's, then we average the matrix rows and colums.
1478C
1479#include "implicit.h"
1480      DIMENSION H1P(NNORBX)
1481C
1482C     Used from common blocks
1483C     INFORB    : NSYM ,NSSYM, NORBSS(*), IORBSS(*), NINFSS(*,3)
1484C     INFIND    : ISSMO(*), ISSORD(*)
1485C
1486#include "maxorb.h"
1487#include "maxash.h"
1488#include "priunit.h"
1489#include "inforb.h"
1490#include "infind.h"
1491#include "infpri.h"
1492C
1493      PARAMETER ( MAXDEG = 20, D0 = 0.0D0, D1 = 1.0D0 )
1494      DIMENSION ISSOFF(MAXDEG),IDEGV(MAXDEG)
1495C
1496      CALL QENTER('AVEH1P')
1497C
1498      IF (IPRAVE .GE. 5) THEN
1499         WRITE(LUPRI,'(/A)')
1500     *      ' AVEH1P test output: H1P matrix before averaging'
1501         CALL OUTPAK(H1P,NORBT,LUPRI)
1502      END IF
1503C
1504      DO 800 ISSYMR = 1, NSSYM
1505C     -- skip if supersym not degenerate or not 'root' supersym
1506         IF ( NINFSS(ISSYMR,1) .EQ. 1) GO TO 800
1507         IF ( NINFSS(ISSYMR,2) .NE. ISSYMR) GO TO 800
1508C
1509         NDEG  = 0
1510         DO 100 JSSYM = ISSYMR, NSSYM
1511C        -- skip if not same 'root' supersymmetry
1512C           (because then not degenerate at all)
1513            IF ( NINFSS(JSSYM,2) .EQ. ISSYMR) THEN
1514               NDEG = NDEG + 1
1515               ISSOFF(NDEG) = IORBSS(JSSYM) - IORBSS(ISSYMR)
1516            END IF
1517  100    CONTINUE
1518         IF ( NDEG .NE. NINFSS(ISSYMR,1) ) THEN
1519            WRITE(LUERR,'(A)') ' AVEH1P error : inconsistent'//
1520     &            ' degeneracies'
1521            WRITE(LUERR,'(A,I4)') ' - ISSYMR           = ', ISSYMR
1522            WRITE(LUERR,'(A,I4)') ' - NDEG             = ', NDEG
1523            WRITE(LUERR,'(A,I4)') ' - NINFSS(ISSYMR,1) = ',
1524     &                 NINFSS(ISSYMR,1)
1525            CALL QUIT(' AVEH1P : inconsistent degeneracy in NINFSS')
1526         END IF
1527C
1528C     -- Compute averaging factor
1529         AVFAC = NDEG
1530         AVFAC = D1 / AVFAC
1531C
1532C     -- in the following, use supersymmetry info to get IMO
1533         IRSYM  = NINFSS(ISSYMR,3)
1534         NORBI  = NORB(IRSYM)
1535         IOFFSR = IORBSS(ISSYMR)
1536         DO 700 IMOSS = IOFFSR+1, IOFFSR+NORBSS(ISSYMR)
1537            IMO  = ISSORD(IMOSS)
1538            DO 600 JMOSS = IOFFSR+1, IMOSS
1539               JMO = ISSORD(JMOSS)
1540               IDEGV(1) = IROW(IMO) + JMO
1541               DO 300 IDEG = 2, NDEG
1542C                 920714-hjaaj: KMO var. is used to be safe
1543C                 xlf version 2.2 has made errors in a similar statem.
1544C                 to the version with IROW(IMO+ISSOFF(IDEG))
1545                  KMO = ISSORD(IMOSS + ISSOFF(IDEG))
1546                  IDEGV(IDEG) = IROW(KMO) + ISSORD(JMOSS + ISSOFF(IDEG))
1547  300          CONTINUE
1548C           -- now collect all elements which should be made equal
1549               XH1PIJ = H1P(IDEGV(1))
1550               DO 400 IDEG = 2, NDEG
1551                  XH1PIJ = XH1PIJ + H1P(IDEGV(IDEG))
1552  400          CONTINUE
1553               XH1PIJ = XH1PIJ * AVFAC
1554C           -- now distribute all elements which should be made equal
1555               DO 500 IDEG = 1, NDEG
1556                  H1P(IDEGV(IDEG)) = XH1PIJ
1557  500          CONTINUE
1558  600       CONTINUE
1559  700    CONTINUE
1560C
1561  800 CONTINUE
1562C
1563C     Zero elements which are not totally symmetric
1564C     with respect to the supersymmetry.
1565C
1566      DO 980 ISYM = 1,NSYM
1567         IORBI = IORB(ISYM)
1568         NORBI = NORB(ISYM)
1569         DO 960 IMO = IORBI+2,IORBI+NORBI
1570            DO 950 JMO = IORBI+1,IMO-1
1571               IF (ISSMO(JMO) .NE. ISSMO(IMO)) THEN
1572                  IJH1P = IROW(IMO) + JMO
1573                  H1P(IJH1P) = D0
1574               END IF
1575  950       CONTINUE
1576  960    CONTINUE
1577  980 CONTINUE
1578C
1579      IF (IPRAVE .GE. 5) THEN
1580         WRITE(LUPRI,'(/A)')
1581     *      ' AVEH1P test output: H1P matrix after averaging'
1582         CALL OUTPAK(H1P,NORBT,1,LUPRI)
1583      END IF
1584C
1585      CALL QEXIT('AVEH1P')
1586C
1587      RETURN
1588      END
1589#endif
1590C  /* Deck avefck */
1591      SUBROUTINE AVEFCK(FCK)
1592C
1593C     12-Jul-1992 Hinne Hettema
1594C
1595C     Purpose : perform an averaging of a one-electron
1596C     matrix. First we identify an IMO which is degenerate with
1597C     other MO's, then we average the matrix rows and colums.
1598C
1599#include "implicit.h"
1600      DIMENSION FCK(NNORBT)
1601C
1602C     Used from common blocks
1603C     INFORB    : NSYM ,NSSYM, NORBSS(*), IORBSS(*), NINFSS(*,3)
1604C     INFIND    : ISSMO(*), ISSORD(*)
1605C
1606#include "maxorb.h"
1607#include "maxash.h"
1608#include "priunit.h"
1609#include "inforb.h"
1610#include "infind.h"
1611#include "infpri.h"
1612C
1613      PARAMETER ( MAXDEG = 20, D0 = 0.0D0, D1 = 1.0D0 )
1614      DIMENSION ISMDEG(MAXDEG),ISSOFF(MAXDEG),IDEGV(MAXDEG)
1615C
1616      CALL QENTER('AVEFCK')
1617C
1618      IF (IPRAVE .GE. 7) THEN
1619         WRITE(LUPRI,'(/A)')
1620     *      ' AVEFCK test output: FCK matrix before averaging'
1621         CALL OUTPKB(FCK,NORB,NSYM,1,LUPRI)
1622      END IF
1623C
1624      DO 800 ISSYMR = 1, NSSYM
1625C     -- skip if supersym not degenerate or not 'root' supersym
1626         IF ( NINFSS(ISSYMR,1) .EQ. 1) GO TO 800
1627         IF ( NINFSS(ISSYMR,2) .NE. ISSYMR) GO TO 800
1628C
1629         NDEG  = 0
1630         DO 100 JSSYM = ISSYMR, NSSYM
1631C        -- skip if not same 'root' supersymmetry
1632C           (because then not degenerate at all)
1633            IF ( NINFSS(JSSYM,2) .EQ. ISSYMR) THEN
1634               NDEG = NDEG + 1
1635               ISMDEG(NDEG) = NINFSS(JSSYM,3)
1636               ISSOFF(NDEG) = IORBSS(JSSYM) - IORBSS(ISSYMR)
1637            END IF
1638  100    CONTINUE
1639         IF ( NDEG .NE. NINFSS(ISSYMR,1) ) THEN
1640            WRITE(LUERR,'(A)') ' AVEFCK error : inconsistent'//
1641     &            ' degeneracies'
1642            WRITE(LUERR,'(A,I4)') ' - ISSYMR           = ', ISSYMR
1643            WRITE(LUERR,'(A,I4)') ' - NDEG             = ', NDEG
1644            WRITE(LUERR,'(A,I4)') ' - NINFSS(ISSYMR,1) = ',
1645     &                 NINFSS(ISSYMR,1)
1646            CALL QUIT(' AVEFCK : inconsistent degeneracy in NINFSS')
1647         END IF
1648C
1649C     -- Compute averaging factor
1650         AVFAC = NDEG
1651         AVFAC = D1 / AVFAC
1652C
1653C     -- in the following, use supersymmetry info to get IMO
1654         IRSYM  = ISMDEG(1)
1655         IORBI  = IORB(IRSYM)
1656         NORBI  = NORB(IRSYM)
1657         IIORBI = IIORB(IRSYM)
1658         IOFFSR = IORBSS(ISSYMR)
1659         DO 700 IMOSS = IOFFSR+1, IOFFSR+NORBSS(ISSYMR)
1660            IMO  = ISSORD(IMOSS) - IORBI
1661            DO 600 JMOSS = IOFFSR+1, IMOSS
1662               JMO = ISSORD(JMOSS) - IORBI
1663               IDEGV(1) = IIORBI + IROW(IMO) + JMO
1664               DO 300 IDEG = 2, NDEG
1665                  JRSYM  = ISMDEG(IDEG)
1666                  IORBJ  = IORB(JRSYM)
1667                  IIORBJ = IIORB(JRSYM)
1668C                 920714-hjaaj: KMO var. is used to be safe
1669C                 xlf version 2.2 has made errors in a similar statem.
1670C                 to the version with IROW(IMO+ISSOFF(IDEG))
1671                  KMO = ISSORD(IMOSS + ISSOFF(IDEG)) - IORBJ
1672                  LMO = ISSORD(JMOSS + ISSOFF(IDEG)) - IORBJ
1673                  IDEGV(IDEG) = IIORBJ + IROW(KMO) + LMO
1674  300          CONTINUE
1675C           -- now collect all elements which should be made equal
1676               XFCKIJ = FCK(IDEGV(1))
1677               DO 400 IDEG = 2, NDEG
1678                  XFCKIJ = XFCKIJ + FCK(IDEGV(IDEG))
1679  400          CONTINUE
1680               XFCKIJ = XFCKIJ * AVFAC
1681C           -- now distribute all elements which should be made equal
1682               DO 500 IDEG = 1, NDEG
1683                  FCK(IDEGV(IDEG)) = XFCKIJ
1684  500          CONTINUE
1685  600       CONTINUE
1686  700    CONTINUE
1687C
1688  800 CONTINUE
1689C
1690C     Zero elements which are not totally symmetric
1691C     with respect to the supersymmetry.
1692C
1693      DO 980 ISYM = 1,NSYM
1694         IORBI  = IORB(ISYM)
1695         NORBI  = NORB(ISYM)
1696         IIORBI = IIORB(ISYM)
1697         DO 960 IMO = 2,NORBI
1698            DO 950 JMO = 1,IMO-1
1699               IF (ISSMO(IORBI+JMO) .NE. ISSMO(IORBI+IMO)) THEN
1700                  IJFCK = IIORBI + IROW(IMO) + JMO
1701                  FCK(IJFCK) = D0
1702               END IF
1703  950       CONTINUE
1704  960    CONTINUE
1705  980 CONTINUE
1706C
1707      IF (IPRAVE .GE. 7) THEN
1708         WRITE(LUPRI,'(/A)')
1709     *      ' AVEFCK test output: FCK matrix after averaging'
1710         CALL OUTPKB(FCK,NORB,NSYM,1,LUPRI)
1711      END IF
1712C
1713      CALL QEXIT('AVEFCK')
1714C
1715      RETURN
1716      END
1717C  /* Deck aveprt */
1718      SUBROUTINE AVEPRT(LUPRI)
1719C
1720C     7-Jul-92 Hinne Hettema and Hans Joergen Aa. Jensen.
1721C     Revised 15-Jul-1993 hjaaj
1722C     Purpose : write supersymmetry information
1723C
1724#include "implicit.h"
1725C
1726#include "maxash.h"
1727#include "maxorb.h"
1728#include "inforb.h"
1729#include "infind.h"
1730C
1731      PARAMETER ( ISTEP = 10 )
1732C
1733      CALL QENTER('AVEPRT')
1734C
1735C     Check if supersymmetry is identical to "D2h" symmetry
1736C
1737      IF (NSSYM .LE. NSYM .AND. MXDGSS .EQ. 1) THEN
1738         DO 200 ISYM = 1,NSYM
1739            JSYM = 0
1740            DO 100 ISSYM = 1,NSSYM
1741               IF (NINFSS(ISSYM,3) .EQ. ISYM) JSYM = JSYM + 1
1742  100       CONTINUE
1743            IF (JSYM .GT. 1) GO TO 300
1744  200    CONTINUE
1745         WRITE(LUPRI,'(/5X,A)') 'No supersymmetry found.'
1746         GO TO 9999
1747      END IF
1748  300 CONTINUE
1749C
1750      WRITE(LUPRI,'(/5X,A)') 'Supersymmetry specification'
1751      WRITE(LUPRI,'( 5X,A)') '==========================='
1752      IF (MXDGSS .EQ. 1)
1753     &WRITE(LUPRI,'(  5X,A)') '(only one dimensional irreps)'
1754      ICOUNT = 0
1755      DO 500 ISSYM = 1, NSSYM, ISTEP
1756         ICOUNT = ICOUNT + 1
1757         ISTRT = ISSYM
1758         IEND  = MIN(NSSYM, ICOUNT * ISTEP)
1759         ILEN = IEND - ISTRT + 1
1760C
1761         WRITE(LUPRI,1100)
1762     &            'Supersymmetry number',
1763     &            (I, I=ISTRT, IEND)
1764         IF (MXDGSS .GT. 1) WRITE(LUPRI,1100)
1765     &            'degenerate with sup.sym.',
1766     &            (NINFSS(I,2), I = ISTRT, IEND)
1767         WRITE(LUPRI,1100)
1768     &            'and located in abelian symmetry ',
1769     &            (NINFSS(I,3), I=ISTRT, IEND)
1770         WRITE(LUPRI,1110)
1771     &            ('  --', I = 1, ILEN)
1772         WRITE(LUPRI,1100)
1773     &            'Total number of orbitals',
1774     &            (NORBSS(I), I = ISTRT, IEND)
1775         WRITE(LUPRI,'()')
1776  500 CONTINUE
1777 9999 CALL QEXIT('AVEPRT')
1778      RETURN
1779 1100 FORMAT(5X,A,T36,10I4)
1780 1110 FORMAT(5X,  T36,10A4)
1781      END
1782C  /* Deck avecph */
1783      SUBROUTINE AVECPH(IPHCHA,CMO,WRK,KFRSAV,LFRSAV)
1784C
1785C 920715-hjaaj+hh
1786C
1787C Purpose: change phases such that degenerate MO's have
1788C          the same relative phase.
1789C
1790#include "implicit.h"
1791      DIMENSION CMO(NCMOT), WRK(*)
1792C
1793C Used from common blocks:
1794C   INFORB : NCMOT,N2ORBX,NINFSS(*,4), NSSYM
1795C
1796#include "inforb.h"
1797C
1798      CALL QENTER('AVECPH')
1799      KFREE = KFRSAV
1800      LFREE = LFRSAV
1801      CALL MEMGET('REAL',KTKMO,N2ORBX,WRK,KFREE,LFREE)
1802      CALL AVETK(WRK(KTKMO),CMO,WRK,KFREE,LFREE)
1803      IPHCHA = 0
1804      DO 100 ISSYM = 1, NSSYM
1805         IDEG   = NINFSS(ISSYM,1)
1806         ISSYMR = NINFSS(ISSYM,2)
1807         IF (IDEG .GT. 1 .AND. ISSYM .EQ. ISSYMR) THEN
1808C        --- we have a degeneracy to other supersym.(s)
1809            CALL AVEPHA(ISSYMR,CMO,WRK(KTKMO),IPHCHA)
1810         END IF
1811  100 CONTINUE
1812      CALL MEMREL('AVECPH',WRK,1,KFRSAV,KFREE,LFREE)
1813      CALL QEXIT('AVECPH')
1814      RETURN
1815      END
1816C  /* Deck avechk */
1817      SUBROUTINE AVECHK(INPERR)
1818C
1819C 14-Jul-1992 Hinne Hettema and Hans Joergen Aa. Jensen.
1820C
1821C Purpose: Check supersymmetry information which has been built up
1822C          up till now.
1823C
1824#include "implicit.h"
1825C
1826C Used from common blocks:
1827C   INFINP : SUPSYM
1828C   INFORB : NSSYM,NORBSS(),IORBSS(),NINFSS(*,3)
1829C   INFIND : ISW(*),ISSORD(*)
1830C
1831#include "maxash.h"
1832#include "maxorb.h"
1833#include "priunit.h"
1834#include "infinp.h"
1835#include "infind.h"
1836#include "inforb.h"
1837#include "infpri.h"
1838#include "orbtypdef.h"
1839C
1840      PARAMETER (MAXDEG = 20, D1 = 1.0D0)
1841      DIMENSION ISSDEG(MAXDEG),ISSOFF(MAXDEG)
1842C
1843      CALL QENTER('AVECHK')
1844      IF (.NOT.SUPSYM) GO TO 9999
1845C
1846      DO 800 ISSYMR = 1,NSSYM
1847C        Skip this supsym if not degenerate or not "root supsym"
1848         IF (NINFSS(ISSYMR,1) .EQ. 1) GO TO 800
1849         IF (NINFSS(ISSYMR,2) .NE. ISSYMR) GO TO 800
1850         NDEG = 0
1851         DO 100 ISSYM = ISSYMR,NSSYM
1852            IF (NINFSS(ISSYM,2) .EQ. ISSYMR) THEN
1853               NDEG = NDEG + 1
1854               ISSDEG(NDEG) = ISSYM
1855               ISSOFF(NDEG) = IORBSS(ISSYM) - IORBSS(ISSYMR)
1856            END IF
1857  100    CONTINUE
1858         IF (NDEG .NE. NINFSS(ISSYMR,1)) THEN
1859            INPERR = INPERR + 1
1860            WRITE(LUPRI,'(/A/A,I4/A,I4,A,I4)')
1861     &      ' Super symmetry error: inconsistent degeneracy in NINFSS',
1862     &      '    for "root" super symmetry',ISSYMR,
1863     &      '    Found:',NDEG,';   Expected:',NINFSS(ISSYMR,1)
1864         END IF
1865C
1866         IOFFSR = IORBSS(ISSYMR)
1867         ITYPA = 0
1868         DO 480 IMOSS = IOFFSR+1,IOFFSR+NORBSS(ISSYMR)
1869            IMO = ISSORD(IMOSS)
1870            ITYP = IOBTYP(IMO)
1871            IF (ITYP .EQ. JTACT) THEN
1872               ITYPA = IACTYP( ICH(IMO) )
1873            END IF
1874            DO 410 IDEG = 2,NDEG
1875               KMO = ISSORD(IMOSS + ISSOFF(IDEG))
1876               KTYP = IOBTYP(KMO)
1877               IF ( ITYP .NE. KTYP ) THEN
1878                  INPERR = INPERR + 1
1879                  WRITE (LUPRI,'(/A,2I5/A,2X,A,2X,A/A,I4)')
1880     &         ' Super symmetry error: degenerate orbitals',IMO,KMO,
1881     &         '    are not of same type',COBTYP(ITYP),COBTYP(KTYP),
1882     &         '    for "root" super symmetry',ISSYMR
1883               ELSE IF (ITYP .EQ. JTACT) THEN
1884                  KTYPA = IACTYP( ICH(KMO) )
1885                  IF (ITYPA .NE. KTYPA) THEN
1886                     INPERR = INPERR + 1
1887                     WRITE (LUPRI,'(/A,2I5/A,2I5/A,I4)')
1888     &      ' Super symmetry error: degenerate active orbitals',IMO,KMO,
1889     &      '    are not in same RAS shell type',ITYPA,KTYPA,
1890     &      '    for "root" super symmetry',ISSYMR
1891                  END IF
1892               END IF
1893  410       CONTINUE
1894  480    CONTINUE
1895  800 CONTINUE
1896C
1897 9999 CONTINUE
1898      CALL QEXIT('AVECHK')
1899      RETURN
1900      END
1901C  /* Deck avedmp */
1902      SUBROUTINE AVEDMP(LUWDMP)
1903C
1904C     2-Jul-92 Hinne Hettema
1905C     Purpose : dump supersymmetry information to output for debugging
1906C
1907#include "implicit.h"
1908C
1909#include "maxash.h"
1910#include "maxorb.h"
1911#include "inforb.h"
1912#include "infind.h"
1913C
1914      CALL QENTER('AVEDMP')
1915      WRITE(LUWDMP,'(/A/A)')
1916     *   '    Dump of super symmetry information',
1917     *   ' ========================================'
1918      WRITE(LUWDMP,'(/A,I4)')
1919     *   ' Number of supersymmetries (NSSYM) is ',NSSYM
1920      WRITE(LUWDMP,'(/A/A)')
1921     *   '  ISS    NORBSS    IORBSS',
1922     *   ' =============================='
1923      DO 100 I = 1, NSSYM
1924         WRITE(LUWDMP,'(1X,I4,2I10)') I, NORBSS(I), IORBSS(I)
1925  100 CONTINUE
1926      WRITE(LUWDMP,'(A)') ' =============================='
1927      WRITE(LUWDMP,'(/A/A)')
1928     *   '    I     ISSMO    ISSORD      ISMO',
1929     *   ' ========================================'
1930      DO 200 I = 1, NORBT
1931         WRITE(LUWDMP,'(1X,I4,3I10)') I, ISSMO(I), ISSORD(I), ISMO(I)
1932  200 CONTINUE
1933      WRITE(LUWDMP,'(A/)') ' ========================================'
1934      CALL QEXIT('AVEDMP')
1935      RETURN
1936      END
1937