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      SUBROUTINE CC_R12NO(T1AM,T2AM,WORK,LWORK)
21*-----------------------------------------------------------------------
22* Purpose: Calculate semi-natural virtual orbitals for R12
23*
24* Note: T1AM, T2AM will be overwritten! T2AM is expected to contain
25*       the integrals (ia|jb) on entry!
26*
27* Christian Neiss, autumn 2005
28*-----------------------------------------------------------------------
29C
30      implicit none
31#include "priunit.h"
32#include "ccorb.h"
33#include "ccsdsym.h"
34#include "ccsdinp.h"
35#include "inftap.h"
36#include "r12int.h"
37C
38      INTEGER KFOCKD,KYMAT,KEIVAL,KEIVEC,KT2SQ,KSCR,KFCKR12
39      INTEGER KCMO,KFCKHF,LUFCK,KOFF1,KOFF2,KOFF3,KCMOTR
40      INTEGER LWORK,LWRK1,LWRK2,KEND1,KEND2
41      INTEGER MATZ,KFV1,KFV2,IERR,ISYM,IDUMMY
42      INTEGER NFCKR12,IFCKR12(8),ICOUNT,IVIR1(8),KFOCKX,KCMOX
43      INTEGER NSYMX,NORBTSX,NBASTX,NLAMDSX,NRHFSX(8),NORBSX(8)
44      INTEGER NR12xVIR,IR12xVIR(8),KFCKMIX,LUNIT
45      LOGICAL LOCDBG
46#if defined (SYS_CRAY)
47      REAL WORK(LWORK),T1AM(NT1AMX),T2AM(NT2AMX)
48#else
49      DOUBLE PRECISION WORK(LWORK),T1AM(NT1AMX),T2AM(NT2AMX)
50#endif
51C
52      PARAMETER (LOCDBG = .FALSE.)
53C
54      CALL QENTER('CC_R12NO')
55      IF (LOCDBG) WRITE(LUPRI,*)'Entered CC_R12NO'
56C
57      NFCKR12 = 0
58      ICOUNT  = 0
59      DO ISYM = 1, NSYM
60        IFCKR12(ISYM) = NFCKR12
61        IVIR1(ISYM)   = ICOUNT
62        NFCKR12 = NFCKR12 + NRXR12(ISYM)*NRXR12(ISYM)
63        ICOUNT  = ICOUNT + NVIR(ISYM)
64      END DO
65C
66      KFOCKD = 1
67      KCMO   = KFOCKD + NORBTS
68      KFCKHF = KCMO  + NLAMDS
69      KYMAT  = KFCKHF + N2BAST
70      KEIVAL = KYMAT + NMATAB(1)
71      KEIVEC = KEIVAL + NVIRT
72      KFCKR12= KEIVEC + NMATAB(1)
73      KEND1  = KFCKR12 + NFCKR12
74      LWRK1  = LWORK - KEND1
75      IF (LWRK1.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
76C
77      !read canonical orbital energies and CMO matrix:
78      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
79     &            .FALSE.)
80      REWIND LUSIFC
81      CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
82      READ (LUSIFC)
83      READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS)
84      READ (LUSIFC) (WORK(KCMO+I-1),I=1,NLAMDS)
85      CALL GPCLOSE(LUSIFC,'KEEP')
86      IF (FROIMP .OR. FROEXP)
87     *  CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND1),LWRK1)
88      CALL FOCK_REORDER(WORK(KFOCKD),WORK(KEND1),LWRK1)
89      CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1)
90      IF (LOCDBG) THEN
91        WRITE(LUPRI,*)'vir-vir block of CMO matrix:'
92        DO ISYM = 1, NSYM
93          KOFF1 = KCMO + IGLMVI(ISYM,ISYM)
94          CALL OUTPUT(WORK(KOFF1),1,NBAS(ISYM),1,NVIR(ISYM),
95     &                NBAS(ISYM),NVIR(ISYM),1,LUPRI)
96        END DO
97      END IF
98C
99      !read SCF AO-Fock matrix:
100      LUFCK = -1
101      CALL GPOPEN(LUFCK,'CC_FCKREF','UNKNOWN',' ','UNFORMATTED',
102     *           IDUMMY,.FALSE.)
103      REWIND(LUFCK)
104      READ(LUFCK)(WORK(KFCKHF + I-1),I = 1,N2BST(1))
105      CALL GPCLOSE(LUFCK,'KEEP' )
106      IF (LOCDBG) THEN
107        WRITE(LUPRI,*)'AO-Fock matrix read in:'
108        DO ISYM = 1, NSYM
109          WRITE(LUPRI,*)' Symmetry block: ',ISYM
110          CALL OUTPUT(WORK(KFCKHF+IAODIS(ISYM,ISYM)),1,NBAS(ISYM),
111     &                1,NBAS(ISYM),NBAS(ISYM),NBAS(ISYM),1,LUPRI)
112        END DO
113      END IF
114C
115      !make MP2-Guess:
116      CALL CCSD_GUESS(T1AM,T2AM,WORK(KFOCKD),IPRINT)
117C
118      !make 2C-E combination of T2AM:
119      KT2SQ = KEND1
120      KEND2 = KT2SQ + NT2SQ(1)
121      LWRK2 = LWORK - KEND2
122      IF (LWRK2.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
123      CALL CC_T2SQ(T2AM,WORK(KT2SQ),1)
124      CALL CCRHS_T2TR(WORK(KT2SQ),WORK(KEND2),LWRK2,1)
125C     CALL DSCAL(NT2SQ(1),2.0D0,WORK(KT2SQ),1)
126C
127      !calculate virtual-virtual block of density matrix:
128      CALL CC_YI(WORK(KYMAT),WORK(KT2SQ),1,T2AM,1,
129     &           WORK(KEND2),LWRK2)
130C     CALL CC_PRSQ(IDUMMY,WORK(KT2SQ),1,0,1)
131C     CALL CC_PRP(IDUMMY,T2AM,1,0,1)
132C     WRITE(LUPRI,*)'YMAT in CC_R12NO:'
133C     DO ISYM = 1, NSYM
134C       WRITE(LUPRI,*)'symmetry block ',ISYM
135C       CALL OUTPUT(WORK(KYMAT+IMATAB(ISYM,ISYM)),1,NVIR(ISYM),
136C    &              1,NVIR(ISYM),NVIR(ISYM),NVIR(ISYM),
137C    &              1,LUPRI)
138C     END DO
139C
140      !allocate work space for transformed CMO matrix:
141      KCMOTR = KEND1
142      KEND1  = KCMOTR + NLAMDS
143      LWRK1  = LWORK - KEND1
144      IF (LWRK1.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
145      CALL DCOPY(NLAMDS,WORK(KCMO),1,WORK(KCMOTR),1)
146C
147      !symmetrize and diagonalize:
148C     CALL DSCAL(NMATAB(1),0.5D0,WORK(KYMAT),1)
149      DO ISYM = 1, NSYM
150        KSCR  = KEND1
151        KEND2 = KSCR + NVIR(ISYM)*NVIR(ISYM)
152        LWRK2 = LWORK - KEND2
153        IF (LWRK2.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
154C
155        CALL TRSREC(NVIR(ISYM),NVIR(ISYM),
156     &              WORK(KYMAT+IMATAB(ISYM,ISYM)),WORK(KSCR))
157        CALL DAXPY(NVIR(ISYM)*NVIR(ISYM),1.0D0,WORK(KSCR),1,
158     &             WORK(KYMAT+IMATAB(ISYM,ISYM)),1)
159        IF (LOCDBG) THEN
160          WRITE(LUPRI,*)'vir-vir block of density matrix before diag.:'
161          CALL OUTPUT(WORK(KYMAT+IMATAB(ISYM,ISYM)),1,NVIR(ISYM),
162     &                1,NVIR(ISYM),NVIR(ISYM),NVIR(ISYM),1,LUPRI)
163        END IF
164C
165        MATZ = 1
166        KFV1  = KEND1
167        KFV2  = KFV1 + NVIR(ISYM)
168        KEND2 = KFV2 + NVIR(ISYM)
169        LWRK2 = LWORK - KEND2
170        IF (LWRK2.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
171        CALL RS(NVIR(ISYM),NVIR(ISYM),WORK(KYMAT+IMATAB(ISYM,ISYM)),
172     &          WORK(KEIVAL+IVIR1(ISYM)),
173     &          MATZ,WORK(KEIVEC+IMATAB(ISYM,ISYM)),
174     &          WORK(KFV1),WORK(KFV2),IERR)
175        IF ( IERR.NE.0 ) THEN
176          WRITE(LUPRI,'(/A,I5)')
177     *    ' EIGENVALUE PROBLEM NOT CONVERGED IN RS, IERR =',IERR
178          CALL QUIT(' CCRED: EIGENVALUE EQUATION NOT CONVERGED ')
179        END IF
180        IF (LOCDBG) THEN
181          WRITE(LUPRI,*)'Eigenvalues before reorder: symmetry=',ISYM
182          CALL OUTPUT(WORK(KEIVAL+IVIR1(ISYM)),1,NVIR(ISYM),
183     &                1,1,NVIR(ISYM),1,1,LUPRI)
184          WRITE(LUPRI,*)'Eigenvectors before reorder:'
185          CALL OUTPUT(WORK(KEIVEC+IMATAB(ISYM,ISYM)),1,NVIR(ISYM),
186     &                1,NVIR(ISYM),NVIR(ISYM),NVIR(ISYM),1,LUPRI)
187        END IF
188        CALL RGORD(NVIR(ISYM),NVIR(ISYM),WORK(KEIVAL+IVIR1(ISYM)),
189     &             WORK(KFV1),WORK(KEIVEC+IMATAB(ISYM,ISYM)),.TRUE.)
190        WRITE(LUPRI,*)'R12-NATVIR: Nat. occ. numbers after ordering: '//
191     &                'symmetry=',ISYM
192        IF (NVIR(ISYM).EQ.0) THEN
193          WRITE(LUPRI,*) 'This block is empty'
194        ELSE
195          CALL OUTPUT(WORK(KEIVAL+IVIR1(ISYM)),1,MIN(10,NVIR(ISYM)),
196     &                1,1,NVIR(ISYM),1,1,LUPRI)
197        END IF
198        IF (LOCDBG) THEN
199          WRITE(LUPRI,*)'Eigenvectors after reorder:'
200          CALL OUTPUT(WORK(KEIVEC+IMATAB(ISYM,ISYM)),1,NVIR(ISYM),
201     &                1,NVIR(ISYM),NVIR(ISYM),NVIR(ISYM),1,LUPRI)
202          CALL FLSHFO(LUPRI)
203        END IF
204C
205        !calculate CMO-matrix for natural virtual orbitals:
206        KOFF1 = KCMO + IGLMVI(ISYM,ISYM)
207        KOFF2 = KEIVEC+IMATAB(ISYM,ISYM)
208        KOFF3 = KCMOTR + IGLMVI(ISYM,ISYM)
209        CALL DGEMM('N','N',NBAS(ISYM),NRXR12(ISYM),NVIR(ISYM),1.0D0,
210     &             WORK(KOFF1),MAX(1,NBAS(ISYM)),
211     &             WORK(KOFF2),MAX(1,NVIR(ISYM)),
212     &             0.0D0,WORK(KOFF3),MAX(1,NBAS(ISYM)))
213C
214      END DO
215C
216      IF (LOCDBG) THEN
217        WRITE(LUPRI,*)'vir-vir block of CMO matrix with natural '//
218     &                'virtual R12-orbitals:'
219        DO ISYM = 1, NSYM
220          KOFF1 = KCMOTR + IGLMVI(ISYM,ISYM)
221          CALL OUTPUT(WORK(KOFF1),1,NBAS(ISYM),1,NVIR(ISYM),
222     &                NBAS(ISYM),NVIR(ISYM),1,LUPRI)
223        END DO
224      END IF
225C
226      !transform Fock matrix into subspace of additional r12-orbitals
227      !(which are chosen from the natural vir. orbials just generated):
228      DO ISYM = 1, NSYM
229        KOFF1 = KCMOTR + IGLMVI(ISYM,ISYM)
230        KOFF2 = KFCKHF + IAODIS(ISYM,ISYM)
231        KOFF3 = KEND1
232        CALL DGEMM('N','N',NBAS(ISYM),NRXR12(ISYM),NBAS(ISYM),1.0D0,
233     &             WORK(KOFF2),MAX(1,NBAS(ISYM)),
234     &             WORK(KOFF1),MAX(1,NBAS(ISYM)),0.0D0,
235     &             WORK(KOFF3),MAX(1,NBAS(ISYM)))
236        IF (LOCDBG) THEN
237          WRITE(LUPRI,*)'CMO part used:'
238          WRITE(LUPRI,*)'symmetry block: ',ISYM
239          CALL OUTPUT(WORK(KOFF1),1,NBAS(ISYM),
240     &                1,NRXR12(ISYM),NBAS(ISYM),NVIR(ISYM),
241     &                1,LUPRI)
242          WRITE(LUPRI,*)'AO-Fock part used:'
243          WRITE(LUPRI,*)'symmetry block: ',ISYM
244          CALL OUTPUT(WORK(KOFF2),1,NBAS(ISYM),
245     &                1,NBAS(ISYM),NBAS(ISYM),NBAS(ISYM),
246     &                1,LUPRI)
247          WRITE(LUPRI,*)'half-transformed R12-Fock-part:'
248          WRITE(LUPRI,*)'symmetry block: ',ISYM
249          CALL OUTPUT(WORK(KOFF3),1,NBAS(ISYM),
250     &                1,NRXR12(ISYM),NBAS(ISYM),NRXR12(ISYM),
251     &                1,LUPRI)
252        END IF
253        CALL DGEMM('T','N',NRXR12(ISYM),NRXR12(ISYM),NBAS(ISYM),1.0D0,
254     &             WORK(KOFF3),MAX(1,NBAS(ISYM)),
255     &             WORK(KOFF1),MAX(1,NBAS(ISYM)),0.0D0,
256     &             WORK(KFCKR12+IFCKR12(ISYM)),MAX(1,NRXR12(ISYM)))
257        IF (LOCDBG) THEN
258          WRITE(LUPRI,*)'R12-Fock-part before diagonalization:'
259          WRITE(LUPRI,*)'symmetry block: ',ISYM
260          CALL OUTPUT(WORK(KFCKR12+IFCKR12(ISYM)),1,NRXR12(ISYM),
261     &                1,NRXR12(ISYM),NRXR12(ISYM),NRXR12(ISYM),
262     &                1,LUPRI)
263        END IF
264        !diagonalize:
265        MATZ = 1
266        KFV1  = KEND1
267        KFV2  = KFV1 + NRXR12(ISYM)
268        KEND2 = KFV2 + NRXR12(ISYM)
269        LWRK2 = LWORK - KEND2
270        IF (LWRK2.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
271        CALL RS(NRXR12(ISYM),NRXR12(ISYM),WORK(KFCKR12+IFCKR12(ISYM)),
272     &          WORK(KEIVAL+IVIR1(ISYM)),
273     &          MATZ,WORK(KEIVEC+IMATAB(ISYM,ISYM)),
274     &          WORK(KFV1),WORK(KFV2),IERR)
275        IF ( IERR.NE.0 ) THEN
276          WRITE(LUPRI,'(/A,I5)')
277     *    ' EIGENVALUE PROBLEM NOT CONVERGED IN RS, IERR =',IERR
278          CALL QUIT(' CCRED: EIGENVALUE EQUATION NOT CONVERGED ')
279        END IF
280        IF (LOCDBG) THEN
281          WRITE(LUPRI,*)'Fock-Eigenvalues after RS: symmetry=',ISYM
282          CALL OUTPUT(WORK(KEIVAL+IVIR1(ISYM)),1,NRXR12(ISYM),
283     &                1,1,NVIR(ISYM),1,1,LUPRI)
284          WRITE(LUPRI,*)'Eigenvectors after RS:'
285          CALL OUTPUT(WORK(KEIVEC+IMATAB(ISYM,ISYM)),1,NRXR12(ISYM),
286     &                1,NRXR12(ISYM),NRXR12(ISYM),NRXR12(ISYM),1,LUPRI)
287        END IF
288        CALL RGORD(NRXR12(ISYM),NRXR12(ISYM),WORK(KEIVAL+IVIR1(ISYM)),
289     &             WORK(KFV1),WORK(KEIVEC+IMATAB(ISYM,ISYM)),.FALSE.)
290        WRITE(LUPRI,*)
291        WRITE(LUPRI,*)'R12-NATVIR: Fock-Eigenvalues after ordering: '//
292     &                'symmetry=',ISYM
293        IF (NRXR12(ISYM).EQ.0) THEN
294          WRITE(LUPRI,*) 'This block is empty'
295        ELSE IF (NVIR(ISYM).LT.NRXR12(ISYM)) THEN
296          WRITE(LUPRI,*) 'You specified more virtual R12 orbitals '//
297     &                   'than there are virtual orbitals in this '//
298     &                   'symmery:'
299          WRITE(LUPRI,*) 'NVIR(',ISYM,') = ',NVIR(ISYM)
300          WRITE(LUPRI,*) 'NRXR12(',ISYM,') = ',NRXR12(ISYM)
301          CALL QUIT('Too many virtual R12 orbitals')
302        ELSE
303          CALL OUTPUT(WORK(KEIVAL+IVIR1(ISYM)),1,NRXR12(ISYM),
304     &                1,1,NVIR(ISYM),1,1,LUPRI)
305        END IF
306        IF (LOCDBG) THEN
307          WRITE(LUPRI,*)'Eigenvectors after reorder:'
308          CALL OUTPUT(WORK(KEIVEC+IMATAB(ISYM,ISYM)),1,NRXR12(ISYM),
309     &                1,NRXR12(ISYM),NRXR12(ISYM),NRXR12(ISYM),1,LUPRI)
310          CALL FLSHFO(LUPRI)
311        END IF
312
313        !calculate semi-natural virtual R12-part of CMO-matrix :
314        KOFF1 = KCMOTR + IGLMVI(ISYM,ISYM)
315        KOFF2 = KEIVEC+IMATAB(ISYM,ISYM)
316        KOFF3 = KEND1
317        CALL DGEMM('N','N',NBAS(ISYM),NRXR12(ISYM),NRXR12(ISYM),1.0D0,
318     &             WORK(KOFF1),MAX(1,NBAS(ISYM)),
319     &             WORK(KOFF2),MAX(1,NRXR12(ISYM)),
320     &             0.0D0,WORK(KOFF3),MAX(1,NBAS(ISYM)))
321        CALL DCOPY(NBAS(ISYM)*NRXR12(ISYM),WORK(KOFF3),1,WORK(KOFF1),1)
322      END DO
323C
324      WRITE(LUPRI,*)
325C
326      IF (IPRINT.GE.5) THEN
327        WRITE(LUPRI,*)'R12-NATVIR: MO coefficient matrix of '//
328     &                'semi-natural virtual R12-orbitals:'
329        DO ISYM = 1, NSYM
330          IF (NRXR12(ISYM).GT.0) THEN
331            WRITE(LUPRI,*) 'symmetry=',ISYM
332            KOFF1 = KCMOTR + IGLMVI(ISYM,ISYM)
333            CALL OUTPUT(WORK(KOFF1),1,NBAS(ISYM),1,NRXR12(ISYM),
334     &                  NBAS(ISYM),NRXR12(ISYM),1,LUPRI)
335          END IF
336        END DO
337      END IF
338C
339      !generate mixed r12-vir Fock-matrix and write it to file:
340      NR12xVIR = 0
341      DO ISYM = 1, NSYM
342        IR12xVIR(ISYM) = NR12xVIR
343        NR12xVIR = NR12xVIR + NVIR(ISYM)*NRXR12(ISYM)
344      END DO
345C
346      KFCKMIX = KEND1
347      KEND2   = KFCKMIX + NR12xVIR
348      LWRK2   = LWORK - KEND2
349      IF (LWRK2.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
350      DO ISYM = 1, NSYM
351        KSCR  = KEND2
352        KEND2 = KSCR + NBAS(ISYM)*NRXR12(ISYM)
353        LWRK2   = LWORK - KEND2
354        IF (LWRK2.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
355        KOFF1 = KCMOTR + IGLMVI(ISYM,ISYM)
356        KOFF2 = KFCKHF + IAODIS(ISYM,ISYM)
357        KOFF3 = KSCR
358        CALL DGEMM('N','N',NBAS(ISYM),NRXR12(ISYM),NBAS(ISYM),1.0D0,
359     &             WORK(KOFF2),MAX(1,NBAS(ISYM)),
360     &             WORK(KOFF1),MAX(1,NBAS(ISYM)),0.0D0,
361     &             WORK(KOFF3),MAX(1,NBAS(ISYM)))
362        KOFF1 = KCMO + IGLMVI(ISYM,ISYM)
363        CALL DGEMM('T','N',NVIR(ISYM),NRXR12(ISYM),NBAS(ISYM),1.0D0,
364     &             WORK(KOFF1),MAX(1,NBAS(ISYM)),
365     &             WORK(KOFF3),MAX(1,NBAS(ISYM)),0.0D0,
366     &             WORK(KFCKMIX+IR12xVIR(ISYM)),MAX(1,NVIR(ISYM)))
367      END DO
368      IF (LOCDBG) THEN
369        WRITE(LUPRI,*) 'R12 x VIR FOCK MATRIX:'
370        DO ISYM = 1, NSYM
371          WRITE(LUPRI,*) 'Symmetry block: ',ISYM
372          CALL OUTPUT(WORK(KFCKMIX+IR12xVIR(ISYM)),1,NVIR(ISYM),
373     &                1,NRXR12(ISYM),NVIR(ISYM),NRXR12(ISYM),
374     &                1,LUPRI)
375        END DO
376      END IF
377      LUNIT = -1
378      CALL GPOPEN(LUNIT,'R12FVIR','UNKNOWN',' ','UNFORMATTED',IDUMMY,
379     &            .FALSE.)
380      WRITE (LUNIT) (WORK(KFCKMIX+I-1), I=1, NR12xVIR)
381      CALL GPCLOSE(LUNIT,'KEEP')
382C
383      !merge results to SIRIUS interface file:
384      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
385     &            .FALSE.)
386      REWIND LUSIFC
387      CALL MOLLAB('FULLBAS ',LUSIFC,LUPRI)
388      READ (LUSIFC) NSYMX,NORBTSX,NBASTX,NLAMDSX,(NRHFSX(I),I=1,NSYM),
389     &              (NORBSX(I),I=1,NSYM)
390      KFOCKX = KEND1
391      KCMOX  = KFOCKX + NORBTSX
392      KEND2  = KCMOX + NLAMDSX
393      LWRK2  = LWORK - KEND2
394      IF (LWRK2.LT.0) CALL QUIT ('Insufficient memory in CC_R12NO')
395      READ (LUSIFC) (WORK(KFOCKX+I-1), I=1,NORBTSX)
396      READ (LUSIFC) (WORK(KCMOX+I-1),I=1,NLAMDSX)
397C     CALL GPCLOSE(LUSIFC,'KEEP')
398      IF (LOCDBG) THEN
399        WRITE(LUPRI,*)'Orbital energies read in:'
400        WRITE(LUPRI,*) (WORK(KFOCKX+I-1),I=1,NORBTSX)
401      END IF
402C
403      KOFF2 = KFOCKX
404      DO ISYM = 1, NSYM
405        KOFF1 = KEIVAL + IVIR1(ISYM)
406        KOFF2 = KOFF2 + NRHFSA(ISYM)
407        CALL DCOPY(NRXR12(ISYM),WORK(KOFF1),1,WORK(KOFF2),1)
408        KOFF2 = KOFF2 + NRXR12(ISYM) + NRHFSA(ISYM)
409C       CALL DCOPY(NRXR12(ISYM),WORK(KOFF1),1,WORK(KOFF2),1)
410        KOFF2 = KOFF2 + NVIR(ISYM) + NORB2(ISYM)
411      END DO
412      IF (LOCDBG) THEN
413        WRITE(LUPRI,*)'Orbital energies write out:'
414        WRITE(LUPRI,*) (WORK(KFOCKX+I-1),I=1,NORBTSX)
415      END IF
416      IF (LOCDBG) THEN
417        WRITE(LUPRI,*)'MO-coefficient matrix read in:'
418        KOFF1 = KCMOX
419        DO ISYM = 1, NSYM
420          WRITE(LUPRI,*)' Symmetry block: ',ISYM
421          CALL OUTPUT(WORK(KOFF1),1,MBAS1(ISYM)+MBAS2(ISYM),
422     &                1,NORBSX(ISYM),MBAS1(ISYM)+MBAS2(ISYM),
423     &                NORBSX(ISYM),1,LUPRI)
424          KOFF1 = KOFF1 + NORBSX(ISYM)*(MBAS1(ISYM)+MBAS2(ISYM))
425        END DO
426      END IF
427      KOFF2 = KCMOX
428C     KOFF3 = KCMOX
429      DO ISYM = 1, NSYM
430        KOFF1 = KCMOTR + IGLMVI(ISYM,ISYM)
431        KOFF2 = KOFF2 + NRHFSA(ISYM)*(MBAS1(ISYM)+MBAS2(ISYM))
432C       KOFF3 = KOFF3 + NRHFSB(ISYM)*(MBAS1(ISYM)+MBAS2(ISYM)) +
433C    &          NRHFSA(ISYM)*(MBAS1(ISYM)+MBAS2(ISYM))
434        DO I = 1, NRXR12(ISYM)
435          CALL DCOPY(MBAS1(ISYM),WORK(KOFF1),1,WORK(KOFF2),1)
436          KOFF2 = KOFF2 + MBAS1(ISYM)+MBAS2(ISYM)
437C         CALL DCOPY(MBAS1(ISYM),WORK(KOFF1),1,WORK(KOFF3),1)
438          KOFF1 = KOFF1 + MBAS1(ISYM)
439C         KOFF3 = KOFF3 + MBAS1(ISYM)+MBAS2(ISYM)
440        END DO
441        KOFF2 = KOFF2 + (NORB1(ISYM)+NORB2(ISYM))*
442     &                  (MBAS1(ISYM)+MBAS2(ISYM))
443C       KOFF3 = KOFF3 + NORB2(ISYM)*(MBAS1(ISYM)+MBAS2(ISYM))
444      END DO
445      IF (LOCDBG) THEN
446        WRITE(LUPRI,*)'MO-coefficient matrix write out:'
447        KOFF1 = KCMOX
448        DO ISYM = 1, NSYM
449          WRITE(LUPRI,*)' Symmetry block: ',ISYM
450          CALL OUTPUT(WORK(KOFF1),1,MBAS1(ISYM)+MBAS2(ISYM),
451     &                1,NORBSX(ISYM),MBAS1(ISYM)+MBAS2(ISYM),
452     &                NORBSX(ISYM),1,LUPRI)
453          KOFF1 = KOFF1 + NORBSX(ISYM)*(MBAS1(ISYM)+MBAS2(ISYM))
454        END DO
455      END IF
456C     CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
457C    &            .FALSE.)
458      REWIND LUSIFC
459      CALL MOLLAB('FULLBAS ',LUSIFC,LUPRI)
460      READ (LUSIFC)
461      WRITE (LUSIFC) (WORK(KFOCKX+I-1), I=1,NORBTSX)
462      WRITE (LUSIFC) (WORK(KCMOX+I-1),I=1,NLAMDSX)
463      CALL GPCLOSE(LUSIFC,'KEEP')
464C
465C     Compute new matrix <ij|r12**2|kl>
466      CALL R12AUX1(WORK(KCMOX),WORK(KEND2),LWRK2)
467C
468      IF (LOCDBG) WRITE(LUPRI,*)'Leaving CC_R12NO'
469      CALL QEXIT('CC_R12NO')
470      RETURN
471      END
472*=======================================================================
473
474