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!
18
19
20C  /* Deck sirloc */
21      SUBROUTINE SIRLOC(CMO,WORK,LWORK)
22C     Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003).
23
24#include "implicit.h"
25      PARAMETER (THRES = 1.0D-10)
26#include "priunit.h"
27#include "maxorb.h"
28#include "maxash.h"
29#include "infinp.h"
30#include "inforb.h"
31#include "infind.h"
32#include "infvar.h"
33#include "r12int.h"
34      DIMENSION CMO(*), WORK(LWORK)
35      LOGICAL ANTIS
36
37      CALL TITLER('BOYS LOCALIZATION','*',125)
38      CALL NUMLOC(NOC,NBA,NB1,NB2)
39      CALL FROLOC(NCO)
40C
41      WRITE(LUPRI,*) 'Number of frozen orbitals    :',NCO
42      WRITE(LUPRI,*) 'Number of localized orbitals :',NOC-NCO
43      WRITE(LUPRI,*) 'Number of molecular orbitals :',NB1
44      WRITE(LUPRI,*) 'Number of auxilliary orbitals:',NB2
45      WRITE(LUPRI,*) 'Number of basis functions    :',NBA
46C
47C Parameters are set for valence orbitals only
48
49      NOC = NOC - NCO
50      NTR  = ((NOC-1)*NOC)/2
51      NN   = NBA*NBA
52      NOC3 = NOC*3
53      NTR3 = NTR*3
54      NCMO = NCMOT - NCO*NBA
55      NKB   = NBA*NOC
56      NKCL  = NOC*NOC
57      NNBA1 = NBA*(NBA+1)/2
58
59C     *****   memory allocation  *****
60
61      KINI   = 1
62      KCMO   = KINI
63      KDIPX  = KCMO   + NCMOT + NBA*(NOC + NCO)
64      KDIPY  = KDIPX  + NNBASX
65      KDIPZ  = KDIPY  + NNBASX
66      KILIFQ = KDIPZ  + NNBASX
67      KIKY   = KILIFQ + NOC
68      KMAXI0  = KIKY   + NBA
69      KMINI0  = KMAXI0  + NN
70      KRI    = KMINI0  + NN
71      KRIJ   = KRI    + NN
72      KCL    = KRIJ   + NTR3
73      KB     = KCL    + NKCL
74      KQPIX  = KB     + NKB
75      KQPJX  = KQPIX  + NOC
76      KPAO   = KQPJX  + NOC
77      KEVC   = KPAO   + NB1*NBA
78      KMOC   = KEVC   + NB1*NBA
79      KOVER1 = KMOC   + (NOC+NCO)*NBA
80      KOVER2 = KOVER1 + NNBAST
81      KEND   = KOVER2 + NN
82      LFREE  = LWORK  - KEND
83
84      KCMOC  = KCMO  + NCO*NBA
85
86C molecular coefficients and dipole moment integrals
87
88      JRDMO = 9
89      CALL READMO(WORK(KCMO),JRDMO)
90C     CALL AROUND('Canonical molecular orbitals')
91C     CALL LOC_PRINTORB(NB1,NBA,WORK(KCMO))
92
93      CALL RDPROP('XDIPLEN ',WORK(KDIPX),ANTIS)
94      CALL RDPROP('YDIPLEN ',WORK(KDIPY),ANTIS)
95      CALL RDPROP('ZDIPLEN ',WORK(KDIPZ),ANTIS)
96
97      CALL BOYLOC(WORK(KCMOC),WORK(KDIPX),WORK(KILIFQ),WORK(KIKY),
98     &            WORK(KMAXI0),WORK(KMINI0),WORK(KRI),WORK(KRIJ),
99     &            WORK(KCL),WORK(KB),WORK(KQPIX),WORK(KQPJX),NCMO,
100     &            NNBASX,NTR,NOC,NBA,WORK(KEND),LFREE)
101
102Cccms CALL DCOPY(NCMOT,WORK(KCMO),1,CMO,1)
103
104C Projected atomic orbitals
105c
106c     NOC = NOC + NCO
107c     CALL RDONEL('OVERLAP ',.TRUE.,WORK(KOVER1),NNBAST)
108c     CALL PROAO(WORK(KPAO),WORK(KEVC),WORK(KCMO),WORK(KMOC),
109c    *           WORK(KOVER1),WORK(KOVER2),NOC,NBA,NB1,NNBA1)
110c
111c     KSHI = KCMO + NBA*NB1
112c     LSHI = NOC*NBA
113c     CALL DZERO(WORK(KSHI),LSHI)
114c
115Cccms KSHI = KCMO + NBA*NOC
116C     LSHI = NB1*NBA
117C     CALL DZERO(WORK(KSHI),LSHI)
118C     DO I = 1, NB1
119C        KSHIP = KPAO + NB1*(I-1)
120C        KSHIC = KCMO + NBA*(NOC+I-1)
121C        CALL DCOPY(NB1,WORK(KSHIP),1,WORK(KSHIC),1)
122C     END DO
123C
124C     CALL AROUND('Localized orbitals')
125C     CALL LOC_PRINTORB(NB1,NBA,WORK(KCMO))
126C
127C *** IF R12EIN COMBINED WITH LOCAL MOS
128
129CWMK  IF (R12EOR) THEN
130
131C *** Memory allocation
132
133         KFCAO  = KEND
134         KFVAO  = KFCAO  + N2BASX
135         KFCMO  = KFVAO  + N2BASX
136         KDCAO  = KFCMO  + N2BASX
137         KDVAO  = KDCAO  + N2BASX
138         KFPAO  = KDVAO  + N2BASX
139         KVEC   = KFPAO  + NNBAST
140         KWRK   = KVEC   + NN
141         LWRK   = LFREE  - KWRK
142
143C *** Tranform AO fock matrix in local MO fock matrix
144
145         CALL FCKDEN(.TRUE.,.FALSE.,WORK(KDCAO),WORK(KDVAO),
146     &                                 WORK(KCMO),DUMMY,WORK(KWRK),LWRK)
147         CALL FCKMAO(.TRUE.,EMCMY,WORK(KFCAO),WORK(KFVAO),
148     &               WORK(KDCAO),WORK(KDVAO),DV_DUMMY,WORK(KCMO),
149     &               WORK(KWRK),LWRK)
150         CALL DCOPY(N2BASX,WORK(KFCAO),1,WORK(KDCAO),1)
151         CALL DGETSP(NBAST,WORK(KDCAO),WORK(KFCAO))
152         CALL UTHU(WORK(KFCAO),WORK(KFCMO),WORK(KCMO),
153     &             WORK(KDCAO),NBA,NB1)
154C ccms diag pao
155C         NTFMO = NOBL*(NOBL+1)/2
156C         NB1T  = NB1 * (NB1+1)/2
157C      DO I = 1, NB1
158C         KSHIP = KCMO + NBA*(NOC+I-1)
159C         KSHIC = KVEC + NB1*(I-1)
160C         CALL DCOPY(NB1,WORK(KSHIP),1,WORK(KSHIC),1)
161C      END DO
162C         CALL LOC_DIAPAO(WORK(KFCMO),WORK(KFPAO),WORK(KVEC),
163C     *                                           NTFMO,NB1T,NB1,NOC)
164C      DO I = 1, NB1
165C         KSHIP = KCMO + NBA*(NOC+I-1)
166C         KSHIC = KVEC + NB1*(I-1)
167C         CALL DCOPY(NB1,WORK(KSHIC),1,WORK(KSHIP),1)
168C      END DO
169C      CALL DCOPY(NBA*NOC,CMO,1,WORK(KCMO),1)
170
171      CALL WRLOCMO(WORK(KCMO),WORK(KFCMO),NOC,NCO,NBA,NB1)
172
173C        WRITE(LUPRI,'(A)') 'SUBROUTINE FCKMAT'
174C        CALL FCKMAT(.TRUE.,DUMMY,WORK(KCMO),EMY,WORK(KFCAO),DUMMY,
175C    *               WORK(KWRK),LWRK)
176C        CALL OUTPAK(WORK(KFCAO),NBAST,1,LUPRI)
177C        CALL WRLOCMO(WORK(KCMO),WORK(KFCAO),NOC,NCO,NBA,NB1)
178c
179c     CALL AROUND('Molecular orbitals and Canonical PAOs')
180c     CALL LOC_PRINTORB(NB1+NOC,NBA,WORK(KCMO))
181c
182CWMK  END IF
183
184      CALL TITLER('END BOYS LOCALIZATION','*',103)
185      RETURN
186      END
187
188#ifdef CURRENTLY_NOT_USED
189C  /* Deck LOC_diapao  */
190      SUBROUTINE LOC_DIAPAO(FCMO,FPAO,VEC,NTFMO,NB1T,NB1,NOC)
191C     Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003).
192#include "implicit.h"
193#include "priunit.h"
194      DIMENSION FCMO(NTFMO),FPAO(NB1T),VEC(NB1,NB1)
195
196      NOBL = NOC + NB1
197      IJ = 0
198      KL = 0
199      DO I = 1, NOBL
200         DO J = 1, I
201            IJ = IJ + 1
202            IF (J.GT.NOC) THEN
203               KL = KL +1
204               FPAO(KL)=FCMO(IJ)
205            END IF
206         END DO
207      END DO
208      CALL JACOBI(FPAO,VEC,NB1,NB1)
209      IJ = 0
210      KL = 0
211      DO I = 1, NOBL
212         DO J = 1, I
213            IJ = IJ + 1
214            IF (J.GT.NOC) THEN
215               KL = KL +1
216               FCMO(IJ)=FPAO(KL)
217            END IF
218         END DO
219      END DO
220
221      END
222#endif   /* CURRENTLY_NOT_USED */
223
224C  /* Deck numloc  */
225      SUBROUTINE NUMLOC(NOC,NBA,NB1,NB2)
226C     Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003).
227C     Revised by hjaaj Aug 2004: set NOC = NISHT instead of NOC = NOCCT
228
229#include "implicit.h"
230#include "inforb.h"
231#include "r12int.h"
232
233C INFORMATION ON CANONICAL OCCUPIED ORBITALS
234C Symmetry is switched off
235
236      NOC = 0
237      NBA = 0
238      NB1 = 0
239      NB2 = 0
240      DO ISYM = 1,NSYM
241         NOC  = NOC + NISH(ISYM)
242C        ... localize all doubly occupied orbitals
243C            but not the active orbitals. /hjaaj aug 2004
244         NBA  = NBA + NBAS(ISYM)
245         NB1  = NB1 + NORB1(ISYM)
246         NB2  = NB2 + NORB2(ISYM)
247      END DO
248
249      RETURN
250      END
251
252C  /* Deck froloc  */
253      SUBROUTINE FROLOC(NCO)
254C     Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003).
255C
256C     Revised by hjaaj Aug 2004 to use NFRO(ISYM) from inforb.h instead
257C     of NRHFFR(ISYM) from ccorb.h which is not defined yet in Sirius!
258C     NOTE that NRHFFR(isym) frozen orbitals from CC from ccorb.h not
259C     necessarily are the first, thus also for that reason it cannot be used.
260C     This will also work for MCSCF.
261C
262C     Revised by wmk Aug 2005 to use LOCFRO(ISYM).
263C
264#include "implicit.h"
265#include "inforb.h"
266#include "r12int.h"
267
268C     LGLO = .TRUE.
269C     hjaaj aug 2004: is in ccorb.h but is never used in cc/, thus
270C     I just disabled the definition.
271      NCO = 0
272      DO ISYM = 1, NSYM
273cwmk     NCO = NCO + NFRO(ISYM)
274         NCO = NCO + LOCFRO(ISYM)
275      END DO
276
277      RETURN
278      END
279
280C  /* Deck boyloc */
281      SUBROUTINE BOYLOC(CMO,DIPXYZ,ILIFQ,IKY,MAXI0,MINI0,RI,RIJ,CL,B,
282     &                  QPIX,QPJX,NCMO,NDIP,NTR,NOC,NBA,WORK,LWORK)
283C Purpose :
284C     To localize molecular orbitals with Boys procedure
285C     Adapted from the GAMESS-UK code
286C
287C     Written by Claire C.M. Samson (University of Utrecht, 30 November 2002).
288
289#include "implicit.h"
290       PARAMETER (D0 = 0.0D0, D1 = 1.0D0, D2 = 2.0D0, D4 = 4.0D0,
291     &            D1M = -1.0D0, D2M = -2.0D0, DP4 = 0.25D0,
292     &            THR10 = 1.0D-10, THR8 = 1.0D-8, THR6= 1.0D-6,
293     &            THR8M = -1.0D-8)
294       DIMENSION WORK(LWORK),CMO(NCMO),DIPXYZ(NDIP,3),ILIFQ(NOC),
295     &           IKY(NBA),MAXI0(NBA,NBA),MINI0(NBA,NBA),RI(NBA,NBA),
296     &           RIJ(NTR,3),CL(NOC,NOC),QPIX(NOC),QPJX(NOC),B(NBA,NOC)
297#include "priunit.h"
298
299C PRINT THE DIPOLE MOMENT MATRICES
300
301      IPRT = 0
302      IF (IPRT .GE. 1) THEN
303         CALL AROUND('DIPX')
304         CALL OUTPAK(DIPXYZ(1,1),NBA,1,LUPRI)
305         CALL AROUND('DIPY')
306         CALL OUTPAK(DIPXYZ(1,2),NBA,1,LUPRI)
307         CALL AROUND('DIPZ')
308         CALL OUTPAK(DIPXYZ(1,3),NBA,1,LUPRI)
309      END IF
310
311C END OF THE PRINT
312
313      IF (NOC.EQ.0) GO TO 500
314
315C Tools to read the triangular density matrix
316
317      DO I= 1, NOC
318         ILIFQ(I)=(I-1)*NBA
319      END DO
320      DO I= 1, NBA
321         IKY(I)=((I-1)*I)/2
322         DO J= 1, NBA
323            IF (J.le.I) THEN
324               MAXI0(I,J)=I
325               MINI0(I,J)=J
326            ELSE
327               MAXI0(I,J)=J
328               MINI0(I,J)=I
329            END IF
330         END DO
331      END DO
332
333
334C     Calculate dipole moment coordinates for the molecular orbitals I,J
335C     (I/COOR/J) ; COOR=X,Y,Z
336C     RI are diagonal elements
337C     RIJ are the off diagonal
338
339      DO IC = 1, 3
340         MM=0
341         DO I = 1, NOC
342            DO J = 1, I
343               IF (I.NE.J) MM = MM + 1
344               SUM = 0.0D0
345               DO K = 1, NBA
346                  KI = K + ILIFQ(I)
347                  CCMO = CMO(KI)
348                  DO L = 1, NBA
349                     LJ = L + ILIFQ(J)
350                     KL = IKY(MAXI0(K,L)) + MINI0(K,L)
351                     SUM = SUM + CCMO*CMO(LJ)*DIPXYZ(KL,IC)
352                  END DO
353               END DO
354               IF (I.EQ.J) THEN
355                  RI(I,IC) = SUM
356               ELSE
357                  RIJ(MM,IC) = SUM
358               END IF
359            END DO
360         END DO
361      END DO
362
363C Now do the rotations
364C ---------------------
365C Initialize the array CL
366
367      DO I = 1 , NOC
368         DO J = 1 , I
369            CL(I,J) = D0
370            CL(J,I) = D0
371         END DO
372         CL(I,I) = D1
373      END DO
374      ITER = 0
375      SHIFT = DATAN(D1)
376 100  CHANGE = D0
377      ITER = ITER + 1
378
379C     2X2 unitary transformation
380C         psi prime(i)=   cos(t)*psi(i)+sin(t)*psi(j)
381C         psi prime (j)= -sin(t)*psi(i)+cos(t)*psi(j).
382C     BOYS LOCALIZATION : maximize the sum of the squares of molecular
383C     dipole moment integrals.
384
385      DO 150 I = 1 , NOC
386         IM = I - 1
387         JM = 1
388         IJM = IM*(IM-1)/2 + 1
389         RM = D0
390         TM = D0
391         SM = D0
392         CM = D1
393         DO 120 J = 1 , NOC
394            IF (I.LT.J) THEN
395               IJ = (J-1)*(J-2)/2 + I
396            ELSE IF (I.EQ.J) THEN
397               GO TO 120
398            ELSE
399               IJ = IM*(IM-1)/2 + J
400            END IF
401            T = D0
402            TX = D0
403            DO KK = 1 , 3
404               T = T + D4*RIJ(IJ,KK)**2 - RI(I,KK)**2 - RI(J,KK)**2
405     &               + D2*RI(I,KK)*RI(J,KK)
406               TX = TX + RIJ(IJ,KK)*(RI(J,KK)-RI(I,KK))
407            END DO
408            IF ((DABS(T).LE.THR10).AND.(DABS(TX).LE.THR10))
409     &          GO TO 120
410            TX = D4*TX
411            T = DATAN2(TX,T)
412            T = T*DP4
413            SIGN = D1
414            IF (T.GT.D0) SIGN = D1M
415            T = T + SIGN*SHIFT
416            ITIM = 0
417 110        S = DSIN(T)
418            ITIM = ITIM + 1
419            CO = DCOS(T)
420            RIN = D0
421            DO KK = 1 , 3
422               QPI = CO*CO*RI(I,KK) + S*S*RI(J,KK) +
423     &               D2*CO*S*RIJ(IJ,KK)
424               QPJ = CO*CO*RI(J,KK) + S*S*RI(I,KK) +
425     &               D2M*CO*S*RIJ(IJ,KK)
426               RIN = RIN + QPI*QPI + QPJ*QPJ - RI(I,KK)**2 - RI(J,KK)**2
427            END DO
428            TTEST = DABS(T) - SHIFT
429            IF ((DABS(T).GT.THR8).AND.(DABS(TTEST).GT.THR8)) THEN
430               IF (RIN.LT.THR8M) THEN
431                  IF (ITIM.LE.1) THEN
432                     SIGN = D1
433                     IF (T.GT.D0) SIGN = D1M
434                     T = T + SHIFT*SIGN
435                     GO TO 110
436                  ELSE
437                     GO TO 160
438                  END IF
439               END IF
440            END IF
441            IF (RIN.GT.RM) THEN
442               IJM = IJ
443               RM = RIN
444               JM = J
445               SM = S
446               CM = CO
447               TM = T
448            END IF
449 120     CONTINUE
450         T = TM
451         RIN = RM
452         S = SM
453         CO = CM
454         J = JM
455         IJ = IJM
456         CHANGE = CHANGE + T*T
457         DO KK = 1 , 3
458            QPI = CO*CO*RI(I,KK) + S*S*RI(J,KK) + D2*CO*S*RIJ(IJ,KK)
459            QPJ = CO*CO*RI(J,KK) + S*S*RI(I,KK) + D2M*CO*S*RIJ(IJ,KK)
460            QPIJ = (CO*CO-S*S)*RIJ(IJ,KK) + CO*S*(RI(J,KK)-RI(I,KK))
461            DO 130 K = 1 , NOC
462               IF (I.LT.K) THEN
463                  IK = (K-1)*(K-2)/2 + I
464               ELSE IF (I.EQ.K) THEN
465                  GO TO 130
466               ELSE
467                  IK = (I-1)*(I-2)/2 + K
468               END IF
469               IF (J.LT.K) THEN
470                  JK = (K-1)*(K-2)/2 + J
471               ELSE IF (J.EQ.K) THEN
472                  GO TO 130
473               ELSE
474                  JK = (J-1)*(J-2)/2 + K
475               END IF
476               QPIX(K) = CO*RIJ(IK,KK) + S*RIJ(JK,KK)
477               QPJX(K) = CO*RIJ(JK,KK) - S*RIJ(IK,KK)
478 130        CONTINUE
479            DO 140 K = 1 , NOC
480               IF (I.LT.K) THEN
481                  IK = (K-1)*(K-2)/2 + I
482               ELSE IF (I.EQ.K) THEN
483                  GO TO 140
484               ELSE
485                  IK = (I-1)*(I-2)/2 + K
486               END IF
487               IF (J.LT.K) THEN
488                  JK = (K-1)*(K-2)/2 + J
489               ELSE IF (J.EQ.K) THEN
490                  GO TO 140
491               ELSE
492                  JK = (J-1)*(J-2)/2 + K
493               END IF
494               RIJ(IK,KK) = QPIX(K)
495               RIJ(JK,KK) = QPJX(K)
496 140        CONTINUE
497            RIN = RIN + QPI + QPJ - RI(I,KK) - RI(J,KK)
498            RI(I,KK) = QPI
499            RI(J,KK) = QPJ
500            RIJ(IJ,KK) = QPIJ
501         END DO
502         DO K = 1,NOC
503            C1 = CO*CL(K,I)+S*CL(K,J)
504            C2 = -S*CL(K,I)+CO*CL(K,J)
505            CL(K,I) = C1
506            CL(K,J) = C2
507         END DO
508 150  CONTINUE
509
510C      if convergence has not been reached start another series
511C      of two center rotations.
512
513      CHANGE = DSQRT(D2*CHANGE/(NOC*(NOC-1)))
514      IF (ITER.LE.75) THEN
515         IF (CHANGE.GE.THR10) GO TO 100
516      END IF
517 160  WRITE(LUPRI,2200) ITER
518      IF (ITER.GE.75 .OR. CHANGE.GT.THR6) THEN
519      WRITE(LUPRI,2300)
520      END IF
521
522C *** load localized molecular orbitals in CMO array
523
524      DO I = 1 , NOC
525         CALL DCOPY(NOC,CL(1,I),1,B(1,I),1)
526      END DO
527
528      DO I = 1 , NOC
529C
530         CALL DZERO(RI,NBA)
531C
532         DO J = 1 , NOC
533            IIJ = ILIFQ(J)
534            CALL DAXPY(NBA,B(J,I),CMO(IIJ+1),1,RI(1,1),1)
535         END DO
536         CALL DCOPY(NBA,RI(1,1),1,B(1,I),1)
537      END DO
538
539      DO I = 1 , NOC
540         II = ILIFQ(I)
541         CALL DCOPY(NBA,B(1,I),1,CMO(II+1),1)
542      END DO
543
544C End of the Localization
545
546 500  CONTINUE
547      RETURN
548 2200 FORMAT(/9x,'** LOCALIZATION CONVERGED AFTER',i3,' ITERATIONS **'/)
549 2300 FORMAT(/'** LOCALIZATION HAS BEEN UNSUCESSFUL **')
550
551      END
552
553C  /* Deck proao  */
554      SUBROUTINE PROAO(PAO,EVC,CMO,OMOC,OVLAT,OVLA,NOC,NBA,NB1,NNBA)
555C     Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003).
556#include "implicit.h"
557      PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, D1M = -1.0D0 )
558      DIMENSION CMO(NBA,NOC),OMOC(NB1,NOC)
559      DIMENSION OVLAT(NNBA),OVLA(NB1,NB1)
560      DIMENSION EVC(NB1,NB1),PAO(NB1,NB1)
561
562      DO J = 1, NOC
563         DO I = 1, NB1
564            OMOC(I,J)=CMO(I,J)
565         END DO
566      END DO
567
568      CALL DZERO(PAO,NB1*NB1)
569      IJ = 0
570      DO I = 1, NB1
571         PAO(I,I) = D1
572         DO J = 1, I
573            IJ = IJ + 1
574            OVLA(I,J) = OVLAT(IJ)
575            OVLA(J,I) = OVLAT(IJ)
576         END DO
577      END DO
578
579      CALL DGEMM('N','T',NB1,NB1,NOC,D1,OMOC,NB1,OMOC,NB1,D0,EVC,NB1)
580      CALL DGEMM('N','N',NB1,NB1,NB1,D1M,EVC,NB1,OVLA,NB1,D1,PAO,NB1)
581      DO I=1,NB1
582       CALL DGEMV('N',NB1,NB1,D1,OVLA,NB1,PAO(1,I),1,D0,EVC,1)
583       XNO=DDOT(NB1,PAO(1,I),1,EVC,1)
584       XNO=D1/DSQRT(XNO)
585       CALL DSCAL(NB1,XNO,PAO(1,I),1)
586      ENDDO
587C test
588C      CALL DGEMM('T','N',NB1,NB1,NB1,D1,PAO,NB1,OVLA,NB1,D0,EVC,NB1)
589C      CALL DGEMM('N','N',NB1,NOC,NB1,D1,EVC,NB1,OMOC,NB1,D0,OMOC,NB1)
590
591      RETURN
592      END
593
594C  /* Deck wrlocmo  */
595      SUBROUTINE WRLOCMO(CMO,FKMAT,NOC,NCO,NBA,NB1)
596C     Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003).
597#include "implicit.h"
598#include "priunit.h"
599      DIMENSION FKMAT(*), CMO(*)
600
601      IJ = 0
602      OPEN(99,FILE='FLOCA')
603      DO I = 1, NB1+NOC
604         DO J = 1, I
605            IJ = IJ + 1
606            IF (J.GT.NCO) WRITE(99,'(D30.20)') FKMAT(IJ)
607         END DO
608      END DO
609      CLOSE(99)
610C
611      OPEN(99,FILE='LOCMO',FORM='FORMATTED')
612      WRITE(99,'(D30.20)') (CMO(IJ), IJ = 1, NBA * (NB1+NOC))
613      CLOSE(99)
614      RETURN
615      END
616C  /* Deck LOC_printorb  */
617      SUBROUTINE LOC_PRINTORB(NENDI,NBASI,ARRA)
618C     Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003).
619#include "implicit.h"
620#include "priunit.h"
621#include "maxorb.h"
622#include "maxash.h"
623#include "infinp.h"
624#include "inforb.h"
625      DIMENSION ARRA(*)
626C
627      ISTBAS = 0
628      IF (NENDI.EQ.0) GO TO 20
629      ICMOI = ICMO(1)
630      ISTORB = IORB(1)
631      IEND=0
632   10 IST =IEND+1
633      ISTMO=IEND*NBASI+ICMOI
634      IEND=IEND+7
635      IF(IEND.GT.NENDI) IEND=NENDI
636      IEMO=NBASI*(IEND-1)+ICMOI
637      WRITE(LUPRI,3100) (I,I=IST,IEND)
638      DO I=1,NBASI
639         JSMO=ISTMO+I
640         JEMO=IEMO+I
641         WRITE(LUPRI,3200) I,CENT(I+ISTBAS),TYPE(I+ISTBAS),
642     *                  (ARRA(J),J=JSMO,JEMO,NBASI)
643      END DO
644      IF (IEND.NE.NENDI) GO TO 10
645   20 CONTINUE
646      RETURN
647 3100 FORMAT(/' Orbital  ',7I9)
648 3200 FORMAT(1X,I3,2X,2A4,7F9.4)
649      END
650