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
19c /* deck cc_ijcb */
20*=====================================================================*
21       SUBROUTINE CC_IJCB( XINT1,   ISY1ALBE,
22     &                     XINT2,   ISY2ALBE,
23     &                     IDEL,    IGAM,
24     &                     X1IJCB,  X1CJIB,
25     &                     X2IJCB,  X2CJIB,
26     &                     XLAMDP1, XLAMDH1, ISYML1,
27     &                     XLAMDP2, XLAMDH2, ISYML2,
28     &                     XLAMDP3, XLAMDH3, ISYML3,
29     &                     XLAMDP4, XLAMDH4, ISYML4,
30     &                     WORK,    LWORK,
31     &                     IOPT,    LDERIV,  LRELAX, LZERO,
32     &                     LNEWTA,  LX2ISQ )
33*---------------------------------------------------------------------*
34*
35* Purpose: drive the generalized transformation to
36*          (ij^|cb) + (ij|c^b) integrals and
37*          (cj^|ib) + (c^j|ib) integrals
38*          for the two-index (**|gam del) approach
39*          assumes three-index arrays XIJCB & XCJIB in core
40*
41*          this routine transforms the indices ij and c, the
42*          transformation of the delta index to b has to be done
43*          from the outside.
44*
45*        XINT1 zero order AO integrals
46*        XINT2 derivative AO integrals
47*        XA1IJCB,  XA1CJIB : TA dependent (zero ord) integrals
48*        XA2IJCB,  XA2CJIB : TA dependent deriv & relax integrals
49*
50*        IOPT=0: (ij^|c del) + (ij|c^ del)  only
51*
52*        IOPT=1: (ij^|c del) + (ij|c^ del) and
53*                (cj^|i del) + (c^j|i del) integrals
54*
55*        IF LDERIV=.TRUE. transform also the derivative integrals g[1]:
56*                           with the XLAMD_1 and XLAMD_3
57*           IOPT=0: (ij^|c del) + (ij|c^ del)
58*           IOPT=1: (ij^|c del) + (ij|c^ del) and (cj^|i del) + (c^j|i del)
59*
60*        IF LRELAX=.TRUE. include relaxation contribution to the
61*                         derivative integrals from the transformation
62*                         of g[0] with XLAMD_1 * XLAMD_2 * XLAMD_3 * XLAMD_4
63*                         (or just reorthonormalization if IRELAX = 0)
64*
65*        IF LZERO=.FALSE. skip calculation of zero-order integrals
66*        IF LX2ISQ = .TRUE. the (al bet| part of X2INT is already full matr.
67*
68*    Written by Sonia Coriani, February 1999
69*    based  on Christof's CC_IAJB
70*    Note that no special case is needed for LAO apart LX2ISQ = TRUE
71*
72*=====================================================================*
73#if defined (IMPLICIT_NONE)
74      IMPLICIT NONE
75#else
76#  include "implicit.h"
77#endif
78#include "priunit.h"
79#include "ccorb.h"
80#include "ccsdsym.h"
81#include "maxorb.h"
82#include "ccisao.h"
83
84#if defined (SYS_CRAY)
85      REAL ONE, ZERO
86#else
87      DOUBLE PRECISION ONE, ZERO
88#endif
89      PARAMETER (ONE = 1.0d0, ZERO = 0.0d0)
90
91      LOGICAL LDERIV, LRELAX, LZERO, LNEWTA, LX2ISQ
92      INTEGER IDEL, IGAM, ISY1ALBE, ISY2ALBE
93      INTEGER ISYML1, ISYML2, ISYML3, ISYML4, LWORK, IOPT
94
95#if defined (SYS_CRAY)
96      REAL XLAMDP1(*), XLAMDH1(*)
97      REAL XLAMDP2(*), XLAMDH2(*)
98      REAL XLAMDP3(*), XLAMDH3(*)
99      REAL XLAMDP4(*), XLAMDH4(*)
100      REAL XINT1(*),  XINT2(*)
101      REAL X1IJCB(*), X1CJIB(*)
102      REAL X2IJCB(*), X2CJIB(*)
103      REAL WORK(LWORK)
104#else
105      DOUBLE PRECISION XLAMDP1(*), XLAMDH1(*)
106      DOUBLE PRECISION XLAMDP2(*), XLAMDH2(*)
107      DOUBLE PRECISION XLAMDP3(*), XLAMDH3(*)
108      DOUBLE PRECISION XLAMDP4(*), XLAMDH4(*)
109      DOUBLE PRECISION XINT1(*), XINT2(*)
110      DOUBLE PRECISION X1IJCB(*), X1CJIB(*)
111      DOUBLE PRECISION X2IJCB(*), X2CJIB(*)
112      DOUBLE PRECISION WORK(LWORK)
113#endif
114
115* local
116
117      INTEGER ISYL11, ISYL12, ISYL13, ISYL14, ISYL23
118      INTEGER ISYDEL, ISYGAM
119      INTEGER KSCR1, KSCR3, KEND1, LWRK1, KDUM
120      PARAMETER (KDUM = + 999 999 999)
121
122*---------------------------------------------------------------------*
123*     set some symmetries
124*---------------------------------------------------------------------*
125
126      ISYDEL  = ISAO(IDEL)
127      ISYGAM  = ISAO(IGAM)
128
129*---------------------------------------------------------------------*
130*   work space allocation:
131*
132*   KSCR1  --  I^{del,gam}(alp bet)
133*
134*   KSCR3  --  I^{del,gam}(alp bet)[1]
135*
136*---------------------------------------------------------------------*
137      KSCR1   = 1
138      KEND1   = KSCR1  + N2BST(ISY1ALBE)
139
140c      IF ((LDERIV .OR. LRELAX).AND.(.NOT.LX2ISQ)) THEN
141
142      IF ((LDERIV).AND.(.NOT.LX2ISQ)) THEN
143        KSCR3  = KEND1
144        KEND1  = KSCR3  + N2BST(ISY2ALBE)
145      ELSE
146        KSCR3  = KDUM
147      END IF
148
149      LWRK1   = LWORK - KEND1
150
151      IF ( LWRK1 .LT. 0) THEN
152        CALL QUIT('Insufficient memory in CC_IJCB.')
153      END IF
154
155*---------------------------------------------------------------------*
156*     square zero-order integral matrix up (alp and bet)
157*     and derivative also if (LDERIV).AND.(.NOT.LX2ISQ)
158*---------------------------------------------------------------------*
159
160      CALL CCSD_SYMSQ(XINT1,ISY1ALBE,WORK(KSCR1))
161
162      IF ((LDERIV).AND.(.NOT.LX2ISQ)) THEN
163        CALL CCSD_SYMSQ(XINT2,ISY2ALBE,WORK(KSCR3))
164      END IF
165
166*---------------------------------------------------------------------*
167*     call routine for actual calculation of transformed integrals
168*---------------------------------------------------------------------*
169      IF (.NOT.LX2ISQ) THEN
170         CALL CCIJCB(WORK(KSCR1),ISY1ALBE,WORK(KSCR3),ISY2ALBE,
171     &              IDEL, IGAM, X1IJCB,  X1CJIB, X2IJCB,  X2CJIB,
172     &              XLAMDP1, XLAMDH1, ISYML1, XLAMDP2, XLAMDH2, ISYML2,
173     &              XLAMDP3, XLAMDH3, ISYML3, XLAMDP4, XLAMDH4, ISYML4,
174     &              WORK(KEND1),  LWRK1,
175     &              IOPT,    LDERIV,  LRELAX, LZERO, LNEWTA )
176      ELSE
177         CALL CCIJCB(WORK(KSCR1),ISY1ALBE,XINT2,ISY2ALBE,
178     &              IDEL, IGAM, X1IJCB,  X1CJIB, X2IJCB,  X2CJIB,
179     &              XLAMDP1, XLAMDH1, ISYML1, XLAMDP2, XLAMDH2, ISYML2,
180     &              XLAMDP3, XLAMDH3, ISYML3, XLAMDP4, XLAMDH4, ISYML4,
181     &              WORK(KEND1),  LWRK1,
182     &              IOPT,    LDERIV,  LRELAX, LZERO, LNEWTA )
183      END IF
184
185*---------------------------------------------------------------------*
186*     return
187*---------------------------------------------------------------------*
188
189      RETURN
190      END
191*=====================================================================*
192*                 END OF SUBROUTINE CC_IJCB                           *
193*=====================================================================*
194c /* deck ccijcb */
195*=====================================================================*
196       SUBROUTINE CCIJCB( X1INT,   ISY1ALBE,
197     &                    X2INT,   ISY2ALBE,
198     &                    IDEL,    IGAM,
199     &                    X1IJCB,  X1CJIB,
200     &                    X2IJCB,  X2CJIB,
201     &                    XLAMDP1, XLAMDH1, ISYML1,
202     &                    XLAMDP2, XLAMDH2, ISYML2,
203     &                    XLAMDP3, XLAMDH3, ISYML3,
204     &                    XLAMDP4, XLAMDH4, ISYML4,
205     &                    WORK,    LWORK,
206     &                    IOPT,    LDERIV,  LRELAX, LZERO,
207     &                    LNEWTA )
208*---------------------------------------------------------------------*
209*
210* Purpose: realise the generalized transformation to
211*          (ij^|cb) + (ij|c^b) integrals and
212*          (cj^|ib) + (c^j|ib) integrals
213*          for the two-index (**|gam del) approach
214*          assumes three-index arrays XIJCB & XCJIB in core
215*
216* See CC_IJCB for details
217*
218*    Written by Sonia Coriani, February 1999
219*
220*=====================================================================*
221#if defined (IMPLICIT_NONE)
222      IMPLICIT NONE
223#else
224#  include "implicit.h"
225#endif
226#include "priunit.h"
227#include "ccorb.h"
228#include "ccsdsym.h"
229#include "maxorb.h"
230#include "ccisao.h"
231
232#if defined (SYS_CRAY)
233      REAL ONE, ZERO
234#else
235      DOUBLE PRECISION ONE, ZERO
236#endif
237      PARAMETER (ONE = 1.0d0, ZERO = 0.0d0)
238
239      LOGICAL LDERIV, LRELAX, LZERO, LNEWTA
240      INTEGER IDEL, IGAM, ISY1ALBE, ISY2ALBE, ISYML1, ISYML2
241      INTEGER ISYML3, ISYML4, LWORK, IOPT, IDUMMY
242
243#if defined (SYS_CRAY)
244      REAL XLAMDP1(*), XLAMDH1(*)
245      REAL XLAMDP2(*), XLAMDH2(*)
246      REAL XLAMDP3(*), XLAMDH3(*)
247      REAL XLAMDP4(*), XLAMDH4(*)
248      REAL X1INT(*), X2INT(*)
249      REAL X1IJCB(*), X1CJIB(*)
250      REAL X2IJCB(*), X2CJIB(*)
251      REAL WORK(LWORK)
252#else
253      DOUBLE PRECISION XLAMDP1(*), XLAMDH1(*)
254      DOUBLE PRECISION XLAMDP2(*), XLAMDH2(*)
255      DOUBLE PRECISION XLAMDP3(*), XLAMDH3(*)
256      DOUBLE PRECISION XLAMDP4(*), XLAMDH4(*)
257      DOUBLE PRECISION X1INT(*), X2INT(*)
258      DOUBLE PRECISION X1IJCB(*), X1CJIB(*)
259      DOUBLE PRECISION X2IJCB(*), X2CJIB(*)
260      DOUBLE PRECISION WORK(LWORK)
261#endif
262
263* local
264
265      INTEGER ISYL11, ISYL12, ISYL13, ISYL14, ISYL23
266      INTEGER ISYALP, ISYBET, ISYDEL, ISYGAM
267      INTEGER ISYALP1, ISYBET1
268      INTEGER ISYMI, ISYMJ, ISYMIJ, ISYMC
269      INTEGER KSCR2, KSCR4, KSCR41, KSCR5
270      INTEGER KSCR6, KSCR61, KSCR7, KEND1, LWRK1
271      INTEGER KOFF1, KOFF2, KOFF3, KOFF4, KOFF41, KLAMD
272      INTEGER KOFF5, KOFF6, KOFF61, KOFF7
273      INTEGER NBASA, NBASB, NVIRC, NRHFI
274      INTEGER ISYIJ1,ISYIJ2,ISYIJ3,ISYIJ4
275      INTEGER ISYCJ2,ISYCJ4,ISYCJ3
276
277*---------------------------------------------------------------------*
278*     set some symmetries
279*---------------------------------------------------------------------*
280
281      ISYDEL  = ISAO(IDEL)
282      ISYGAM  = ISAO(IGAM)
283
284*---------------------------------------------------------------------*
285*   work space allocation:
286*
287*   KSCR2  --  I^{del,gam}(alp j); (alp j^); (alp j-bar); (alp j^-bar)
288*   KSCR4  --  I^{del,gam}(i j)
289*   KSCR41 --  I^{del,gam}(i j^)
290*
291*   KSCR6  --  I^{del,gam}([i-bar j + i j-bar]) + I^[1]^{del,gam}(i j)
292*   KSCR61 --  I^{del,gam}([i-bar j^ + i j^-bar]) + I^[1]^{del,gam}(i j^)
293*
294*   KSCR5 --  I^{del,gam}([c^ j  + c j^])
295*   KSCR7 --  I^{del,gam}([c^-bar j + c-bar j^ + c j^-bar + c^ j-bar])
296*             + I^[1]^{del,gam}(c^ j + c j^)
297*
298*---------------------------------------------------------------------*
299      KSCR2   = 1
300      KSCR4   = KSCR2  + NBAST*NRHFT
301      KSCR41  = KSCR4  + NRHFT*NRHFT
302      KEND1   = KSCR41 + NRHFT*NRHFT
303
304      IF (IOPT.EQ.1) THEN
305        KSCR5 = KEND1
306        KEND1 = KSCR5 + NVIRT*NRHFT
307      END IF
308
309      IF (LDERIV .OR. LRELAX) THEN
310        KSCR6  = KEND1
311        KSCR61 = KSCR6  + NRHFT*NRHFT
312        KEND1  = KSCR61 + NRHFT*NRHFT
313
314        IF (IOPT.EQ.1) THEN
315          KSCR7 = KEND1
316          KEND1 = KSCR7 + NVIRT*NRHFT
317        END IF
318      END IF
319
320      LWRK1   = LWORK - KEND1
321
322      IF ( LWRK1 .LT. 0) THEN
323        CALL QUIT('Insufficient memory in CC_IJCB.')
324      END IF
325
326*---------------------------------------------------------------------*
327*     transform beta index to J using XLAMDH1
328*      -- store (alf J|gam del)  in SCR2
329*---------------------------------------------------------------------*
330c**
331c      ISYSCR2 = MULD2H(ISY1ALBE,ISYML1)
332c      CALL DZERO(NT1AO(ISYSCR2),WORK(KSCR2))
333c**
334      KOFF2 = KSCR2
335      DO ISYMJ = 1, NSYM
336
337        ISYBET = MULD2H(ISYML1,ISYMJ)
338        ISYALP = MULD2H(ISYBET,ISY1ALBE)
339
340        KOFF1 = 1 + IAODIS(ISYALP,ISYBET)
341        KLAMD = IGLMRH(ISYBET,ISYMJ) + 1
342
343        NBASA = MAX(NBAS(ISYALP),1)
344        NBASB = MAX(NBAS(ISYBET),1)
345
346        CALL DGEMM('N','N',NBAS(ISYALP),NRHF(ISYMJ),NBAS(ISYBET),
347     *             ONE,X1INT(KOFF1),NBASA,XLAMDH1(KLAMD),
348     *             NBASB,ZERO,WORK(KOFF2),NBASA)
349
350        KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
351
352      END DO
353*---------------------------------------------------------------------*
354*     transform alpha index to I using XLAMDP1
355*      -- store (i j|gam del) in SCR4
356*---------------------------------------------------------------------*
357      KOFF2 = KSCR2
358      DO ISYMJ = 1, NSYM
359
360        ISYBET = MULD2H(ISYML1,ISYMJ)
361        ISYALP = MULD2H(ISYBET,ISY1ALBE)
362        ISYMI  = MULD2H(ISYML1,ISYALP)
363
364        KLAMD = IGLMRH(ISYALP,ISYMI) + 1
365        KOFF4 = KSCR4 + IMATIJ(ISYMI,ISYMJ)
366
367        NBASA = MAX(NBAS(ISYALP),1)
368        NRHFI = MAX(NRHF(ISYMI),1)
369
370        CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYALP),
371     *             ONE,XLAMDP1(KLAMD),NBASA,WORK(KOFF2),
372     *             NBASA,ZERO,WORK(KOFF4),NRHFI)
373
374        KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
375
376      END DO
377
378*---------------------------------------------------------------------*
379* LRELAX   transform alpha index to i-bar using XLAMDP2
380*          -- store (i-bar j|gam del) in SCR6
381*---------------------------------------------------------------------*
382      IF ( LRELAX ) THEN
383
384         KOFF2 = KSCR2
385
386         DO ISYMJ = 1, NSYM
387
388           ISYBET = MULD2H(ISYML1,ISYMJ)
389           ISYALP = MULD2H(ISYBET,ISY1ALBE)
390           ISYMI  = MULD2H(ISYML2,ISYALP)
391
392           KLAMD = IGLMRH(ISYALP,ISYMI) + 1
393           KOFF6 = KSCR6 + IMATIJ(ISYMI,ISYMJ)
394
395           NBASA = MAX(NBAS(ISYALP),1)
396           NRHFI = MAX(NRHF(ISYMI),1)
397
398           CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),
399     *                NBAS(ISYALP),ONE,XLAMDP2(KLAMD),NBASA,
400     *                WORK(KOFF2),NBASA,ZERO,WORK(KOFF6),NRHFI)
401
402           KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
403
404         END DO
405
406      END IF
407
408*---------------------------------------------------------------------*
409*     for IOPT=1 transform alpha index to C^ using XLAMDP3
410*        -- store (c^ j|gam del) in SCR5
411*---------------------------------------------------------------------*
412      IF ( IOPT.EQ.1 ) THEN
413
414        KOFF2 = KSCR2
415
416        DO ISYMJ = 1, NSYM
417
418          ISYBET = MULD2H(ISYML1,ISYMJ)
419          ISYALP = MULD2H(ISYBET,ISY1ALBE)
420          ISYMC  = MULD2H(ISYML3,ISYALP)
421
422          KLAMD = IGLMVI(ISYALP,ISYMC) + 1
423          KOFF5 = KSCR5 + IT1AM(ISYMC,ISYMJ)
424
425          NBASA = MAX(NBAS(ISYALP),1)
426          NVIRC = MAX(NVIR(ISYMC),1)
427
428          CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),NBAS(ISYALP),
429     *               ONE,XLAMDP3(KLAMD),NBASA,WORK(KOFF2),NBASA,
430     *               ZERO,WORK(KOFF5),NVIRC)
431
432          KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
433c
434        END DO
435
436      END IF
437*---------------------------------------------------------------------*
438*     for IOPT=1 and LRELAX transform alpha index to C^(bar)
439*     using XLAMDP4 -- store (c^-bar j|gam del) in SCR7
440*---------------------------------------------------------------------*
441      IF (( IOPT.EQ.1 ).AND.LRELAX) THEN
442
443        KOFF2 = KSCR2
444
445        DO ISYMJ = 1, NSYM
446
447          ISYBET = MULD2H(ISYML1,ISYMJ)
448          ISYALP = MULD2H(ISYBET,ISY1ALBE)
449          ISYMC = MULD2H(ISYML4,ISYALP)
450
451          KLAMD = IGLMVI(ISYALP,ISYMC) + 1
452          KOFF7 = KSCR7 + IT1AM(ISYMC,ISYMJ)
453
454          NBASA = MAX(NBAS(ISYALP),1)
455          NVIRC = MAX(NVIR(ISYMC),1)
456
457          CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),NBAS(ISYALP),
458     *               ONE,XLAMDP4(KLAMD),NBASA,WORK(KOFF2),NBASA,
459     *               ZERO,WORK(KOFF7),NVIRC)
460
461          KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
462
463        END DO
464
465      END IF
466
467*---------------------------------------------------------------------*
468*     transform beta index to J^ using XLAMDH3
469*      -- store (alf J^|gam del)  in SCR2
470*---------------------------------------------------------------------*
471      KOFF2 = KSCR2
472
473      DO ISYMJ = 1, NSYM
474
475        ISYBET = MULD2H(ISYML3,ISYMJ)
476        ISYALP = MULD2H(ISYBET,ISY1ALBE)
477
478        KOFF1 = 1 + IAODIS(ISYALP,ISYBET)
479        KLAMD = IGLMRH(ISYBET,ISYMJ) + 1
480
481        NBASA = MAX(NBAS(ISYALP),1)
482        NBASB = MAX(NBAS(ISYBET),1)
483
484        CALL DGEMM('N','N',NBAS(ISYALP),NRHF(ISYMJ),NBAS(ISYBET),
485     *             ONE,X1INT(KOFF1),NBASA,XLAMDH3(KLAMD),
486     *             NBASB,ZERO,WORK(KOFF2),NBASA)
487
488        KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
489
490      END DO
491
492*---------------------------------------------------------------------*
493*     transform alpha index to I using XLAMDP1
494*      -- store (i j^|gam del) in SCR41
495*---------------------------------------------------------------------*
496      KOFF2 = KSCR2
497
498      DO ISYMJ = 1, NSYM
499
500        ISYBET = MULD2H(ISYML3,ISYMJ)
501        ISYALP = MULD2H(ISYBET,ISY1ALBE)
502        ISYMI  = MULD2H(ISYML1,ISYALP)
503
504        KLAMD  = IGLMRH(ISYALP,ISYMI) + 1
505        KOFF41 = KSCR41 + IMATIJ(ISYMI,ISYMJ)
506
507        NBASA = MAX(NBAS(ISYALP),1)
508        NRHFI = MAX(NRHF(ISYMI),1)
509
510        CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYALP),
511     *             ONE,XLAMDP1(KLAMD),NBASA,WORK(KOFF2),
512     *             NBASA,ZERO,WORK(KOFF41),NRHFI)
513
514        KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
515
516      END DO
517
518*---------------------------------------------------------------------*
519*     if LRELAX transform alpha index to i-bar using XLAMDP2
520*        -- store (i-bar j^| gam del) in SCR61
521*---------------------------------------------------------------------*
522      IF ( LRELAX ) THEN
523
524        KOFF2 = KSCR2
525
526        DO ISYMJ = 1, NSYM
527
528          ISYBET = MULD2H(ISYML3,ISYMJ)
529          ISYALP = MULD2H(ISYBET,ISY1ALBE)
530          ISYMI  = MULD2H(ISYML2,ISYALP)
531
532          KLAMD  = IGLMRH(ISYALP,ISYMI) + 1
533          KOFF61 = KSCR61 + IMATIJ(ISYMI,ISYMJ)
534
535          NBASA = MAX(NBAS(ISYALP),1)
536          NRHFI = MAX(NRHF(ISYMI),1)
537
538          CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYALP),
539     *               ONE,XLAMDP2(KLAMD),NBASA,WORK(KOFF2),NBASA,
540     *               ZERO,WORK(KOFF61),NRHFI)
541
542          KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
543
544        END DO
545
546      END IF
547
548*---------------------------------------------------------------------*
549*     for IOPT=1 transform alpha index to C using XLAMDP1
550*     -- add (c j^|gam del) to (c^ j|gam del) in SCR5
551*---------------------------------------------------------------------*
552      IF ( IOPT.EQ.1 ) THEN
553
554        KOFF2 = KSCR2
555
556        DO ISYMJ = 1, NSYM
557
558          ISYBET = MULD2H(ISYML3,ISYMJ)
559          ISYALP = MULD2H(ISYBET,ISY1ALBE)
560          ISYMC  = MULD2H(ISYML1,ISYALP)
561
562          KLAMD = IGLMVI(ISYALP,ISYMC) + 1
563          KOFF5 = KSCR5 + IT1AM(ISYMC,ISYMJ)
564
565          NBASA = MAX(NBAS(ISYALP),1)
566          NVIRC = MAX(NVIR(ISYMC),1)
567
568          CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),NBAS(ISYALP),
569     *               ONE,XLAMDP1(KLAMD),NBASA,WORK(KOFF2),NBASA,
570     *               ONE,WORK(KOFF5),NVIRC)
571
572          KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
573
574        END DO
575      END IF
576*---------------------------------------------------------------------*
577*     for IOPT=1.and.LRELAX  transform alpha to C-bar using XLAMDP2
578*     -- add (c-bar j^|gam del) to (c^-bar j|gam del) in SCR7
579*---------------------------------------------------------------------*
580      IF (( IOPT.EQ.1 ).AND.(LRELAX)) THEN
581
582        KOFF2 = KSCR2
583
584        DO ISYMJ = 1, NSYM
585
586          ISYBET = MULD2H(ISYML3,ISYMJ)
587          ISYALP = MULD2H(ISYBET,ISY1ALBE)
588          ISYMC  = MULD2H(ISYML2,ISYALP)
589
590          KLAMD = IGLMVI(ISYALP,ISYMC) + 1
591          KOFF7 = KSCR7 + IT1AM(ISYMC,ISYMJ)
592
593          NBASA = MAX(NBAS(ISYALP),1)
594          NVIRC = MAX(NVIR(ISYMC),1)
595
596          CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),NBAS(ISYALP),
597     *               ONE,XLAMDP2(KLAMD),NBASA,WORK(KOFF2),NBASA,
598     *               ONE,WORK(KOFF7),NVIRC)
599
600          KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
601
602        END DO
603
604      END IF
605C
606C Finished with SCR2 again
607C
608*---------------------------------------------------------------------*
609*     for LRELAX add extra contributions from (alp j-bar|gam del),
610*     (alp ^j-bar|gam del) etc
611*---------------------------------------------------------------------*
612      IF ( LRELAX ) THEN
613*---------------------------------------------------------------------*
614*     transform beta to J-bar using XLAMDH2
615*     -- store (alp j-bar|gam del) in SCR2
616*     If (LDERIV) add also (alp j|gam del)[1]
617*---------------------------------------------------------------------*
618
619         KOFF2 = KSCR2
620
621         DO ISYMJ = 1, NSYM
622
623            ISYBET = MULD2H(ISYML2,ISYMJ)
624            ISYALP = MULD2H(ISYBET,ISY1ALBE)
625
626            KOFF1 = 1 + IAODIS(ISYALP,ISYBET)
627            KLAMD = IGLMRH(ISYBET,ISYMJ) + 1
628
629            NBASA = MAX(NBAS(ISYALP),1)
630            NBASB = MAX(NBAS(ISYBET),1)
631
632            CALL DGEMM('N','N',NBAS(ISYALP),NRHF(ISYMJ),NBAS(ISYBET),
633     *                 ONE,X1INT(KOFF1),NBASA,XLAMDH2(KLAMD),
634     *                 NBASB,ZERO,WORK(KOFF2),NBASA)
635
636            IF (LDERIV) THEN
637
638              ISYBET1 = MULD2H(ISYML1,ISYMJ)
639              ISYALP1 = MULD2H(ISYBET1,ISY2ALBE)
640              IF (ISYALP1.NE.ISYALP)
641     *          CALL QUIT('Symmetry mismatch in CC_IJCB')
642
643              KOFF3  = 1 + IAODIS(ISYALP1,ISYBET1)
644              KLAMD  = IGLMRH(ISYBET1,ISYMJ) + 1
645              NBASA = MAX(NBAS(ISYALP1),1)
646              NBASB = MAX(NBAS(ISYBET1),1)
647              CALL DGEMM('N','N',NBAS(ISYALP1),NRHF(ISYMJ),
648     *                  NBAS(ISYBET1),ONE,X2INT(KOFF3),NBASA,
649     *                  XLAMDH1(KLAMD),NBASB,ONE,WORK(KOFF2),NBASA)
650            END IF
651
652
653            KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
654
655          END DO
656
657*---------------------------------------------------------------------*
658*     transform alpha to I using XLAMDP1
659*     -- add (i j-bar|gam del) to (i-bar j|gam del) in SCR6
660*---------------------------------------------------------------------*
661
662          KOFF2 = KSCR2
663
664          DO ISYMJ = 1, NSYM
665
666             ISYBET = MULD2H(ISYML2,ISYMJ)
667             ISYALP = MULD2H(ISYBET,ISY1ALBE)
668             ISYMI  = MULD2H(ISYML1,ISYALP)
669
670             KLAMD = IGLMRH(ISYALP,ISYMI) + 1
671             KOFF6 = KSCR6 + IMATIJ(ISYMI,ISYMJ)
672
673             NBASA = MAX(NBAS(ISYALP),1)
674             NRHFI = MAX(NRHF(ISYMI),1)
675
676             CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),
677     *                  NBAS(ISYALP),ONE,XLAMDP1(KLAMD),NBASA,
678     *                  WORK(KOFF2), NBASA,ONE,WORK(KOFF6),NRHFI)
679
680             KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
681
682           END DO
683*---------------------------------------------------------------------*
684*          for IOPT=1  transform alpha to C^ using XLAMDP3
685*          -- add (c^ j-bar | gam del) in SCR7
686*---------------------------------------------------------------------*
687           IF ( IOPT.EQ.1 ) THEN
688
689           KOFF2 = KSCR2
690
691           DO ISYMJ = 1, NSYM
692
693             ISYBET = MULD2H(ISYML2,ISYMJ)
694             ISYALP = MULD2H(ISYBET,ISY1ALBE)
695             ISYMC  = MULD2H(ISYML3,ISYALP)
696
697             KLAMD = IGLMVI(ISYALP,ISYMC) + 1
698             KOFF7 = KSCR7 + IT1AM(ISYMC,ISYMJ)
699
700             NBASA = MAX(NBAS(ISYALP),1)
701             NVIRC = MAX(NVIR(ISYMC),1)
702
703             CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),
704     *               NBAS(ISYALP),ONE,XLAMDP3(KLAMD),NBASA,
705     *               WORK(KOFF2),NBASA,ONE,WORK(KOFF7),NVIRC)
706
707             KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
708
709           END DO
710           END IF
711
712*---------------------------------------------------------------------*
713*     transform beta to J^-bar using XLAMDH4
714*     -- store (alp j^-bar|gam del) in SCR2
715* if (LDERIV) add derivative contribution (al j^|gam del)(1)
716*---------------------------------------------------------------------*
717
718         KOFF2 = KSCR2
719
720         DO ISYMJ = 1, NSYM
721
722            ISYBET = MULD2H(ISYML4,ISYMJ)
723            ISYALP = MULD2H(ISYBET,ISY1ALBE)
724
725            KOFF1 = 1 + IAODIS(ISYALP,ISYBET)
726            KLAMD = IGLMRH(ISYBET,ISYMJ) + 1
727
728            NBASA = MAX(NBAS(ISYALP),1)
729            NBASB = MAX(NBAS(ISYBET),1)
730
731            CALL DGEMM('N','N',NBAS(ISYALP),NRHF(ISYMJ),NBAS(ISYBET),
732     *                 ONE,X1INT(KOFF1),NBASA,XLAMDH4(KLAMD),
733     *                 NBASB,ZERO,WORK(KOFF2),NBASA)
734
735            IF (LDERIV) THEN
736
737              ISYBET1 = MULD2H(ISYML3,ISYMJ)
738              ISYALP1 = MULD2H(ISYBET1,ISY2ALBE)
739              IF (ISYALP1.NE.ISYALP)
740     *           CALL QUIT('Symmetry mismatch in CC_IJCB')
741
742              KOFF3  = 1 + IAODIS(ISYALP1,ISYBET1)
743              KLAMD  = IGLMRH(ISYBET1,ISYMJ) + 1
744              NBASA = MAX(NBAS(ISYALP1),1)
745              NBASB = MAX(NBAS(ISYBET1),1)
746              CALL DGEMM('N','N',NBAS(ISYALP1),NRHF(ISYMJ),
747     *                  NBAS(ISYBET1),ONE,X2INT(KOFF3),NBASA,
748     *                  XLAMDH3(KLAMD),NBASB,ONE,WORK(KOFF2),NBASA)
749            END IF
750
751            KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
752
753          END DO
754*---------------------------------------------------------------------*
755*     transform alpha to I using XLAMDP1
756*     -- add (i j^-bar|gam del) to (i-bar j^|gam del) in SCR61
757*---------------------------------------------------------------------*
758
759          KOFF2 = KSCR2
760
761          DO ISYMJ = 1, NSYM
762
763             ISYBET = MULD2H(ISYML4,ISYMJ)
764             ISYALP = MULD2H(ISYBET,ISY1ALBE)
765             ISYMI  = MULD2H(ISYML1,ISYALP)
766
767             KLAMD  = IGLMRH(ISYALP,ISYMI) + 1
768             KOFF61 = KSCR61 + IMATIJ(ISYMI,ISYMJ)
769
770             NBASA = MAX(NBAS(ISYALP),1)
771             NRHFI = MAX(NRHF(ISYMI),1)
772
773             CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),
774     *                 NBAS(ISYALP),ONE,XLAMDP1(KLAMD),NBASA,
775     *                 WORK(KOFF2),NBASA,ONE,WORK(KOFF61),NRHFI)
776
777             KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
778
779           END DO
780*---------------------------------------------------------------------*
781*          for IOPT=1  transform alpha to C using XLAMDP1
782*          -- add (c j^-bar | gam del) in SCR7
783*---------------------------------------------------------------------*
784           IF ( IOPT.EQ.1 ) THEN
785
786           KOFF2 = KSCR2
787
788           DO ISYMJ = 1, NSYM
789
790             ISYBET = MULD2H(ISYML4,ISYMJ)
791             ISYALP = MULD2H(ISYBET,ISY1ALBE)
792             ISYMC  = MULD2H(ISYML1,ISYALP)
793
794             KLAMD = IGLMVI(ISYALP,ISYMC) + 1
795             KOFF7 = KSCR7 + IT1AM(ISYMC,ISYMJ)
796
797             NBASA = MAX(NBAS(ISYALP),1)
798             NVIRC = MAX(NVIR(ISYMC),1)
799
800             CALL DGEMM('T','N',NVIR(ISYMC),NRHF(ISYMJ),
801     *               NBAS(ISYALP),ONE,XLAMDP1(KLAMD),NBASA,
802     *               WORK(KOFF2),NBASA,ONE,WORK(KOFF7),NVIRC)
803
804             KOFF2 = KOFF2 + NBAS(ISYALP)*NRHF(ISYMJ)
805
806           END DO
807           END IF
808      END IF
809*---------------------------------------------------------------------*
810*     Add the contribution to the result X1IJCB and X2IJCB vector:
811*     (transform the gamma index to VIRTUAL)
812*     Note that IJCB integrals are sorted as I_cj,i,del
813*---------------------------------------------------------------------*
814      ISYL11 = MULD2H(ISYML1,ISYML1)
815      ISYL12 = MULD2H(ISYML1,ISYML2)
816      ISYL13 = MULD2H(ISYML1,ISYML3)
817      ISYL14 = MULD2H(ISYML1,ISYML4)
818      ISYL23 = MULD2H(ISYML2,ISYML3)
819
820      IF ( LNEWTA ) THEN
821C        -------------------------
822C        add (ij|c^ del) to X1IJCB:
823C        -------------------------
824         ISYIJ1 = MULD2H(ISY1ALBE,ISYL11)
825
826         CALL CC_IJCB2(IGAM, WORK(KSCR4), ISYIJ1, ISYGAM,
827     &                 XLAMDP3, ISYML3, X1IJCB)
828C        -------------------------
829C        add (ij^|c del) to X1IJCB:
830C        -------------------------
831         ISYIJ3 = MULD2H(ISY1ALBE,ISYL13)
832         CALL CC_IJCB2(IGAM, WORK(KSCR41), ISYIJ3, ISYGAM,
833     &                 XLAMDP1, ISYML1, X1IJCB)
834      END IF
835
836
837      IF ( LRELAX ) THEN
838C        ------------------------------
839C        add (ij|c^-bar del) to X2IJCB:
840C        ------------------------------
841         ISYIJ1 = MULD2H(ISY1ALBE,ISYL11)
842
843         CALL CC_IJCB2(IGAM, WORK(KSCR4), ISYIJ1, ISYGAM,
844     &                 XLAMDP4, ISYML4, X2IJCB)
845C        ------------------------------
846C        add (ij^|c-bar del) to X2IJCB:
847C        ------------------------------
848         ISYIJ3 = MULD2H(ISY1ALBE,ISYL13)
849
850         CALL CC_IJCB2(IGAM, WORK(KSCR41), ISYIJ3, ISYGAM,
851     &                 XLAMDP2, ISYML2, X2IJCB)
852      END IF
853      IF ( LDERIV .OR. LRELAX ) THEN
854C        ------------------------------------------------
855C        add ([i-bar j^ + i j^-bar]|c del) to X2IJCB:
856C        ------------------------------------------------
857         ISYIJ4 = MULD2H(ISY1ALBE,ISYL14)
858
859         CALL CC_IJCB2(IGAM, WORK(KSCR61), ISYIJ4, ISYGAM,
860     &                 XLAMDP1, ISYML1, X2IJCB)
861C        ------------------------------------------------
862C        add ([i-bar j + i j-bar]|c^ del) to X2IJCB:
863C        ------------------------------------------------
864         ISYIJ2 = MULD2H(ISY1ALBE,ISYL12)
865
866         CALL CC_IJCB2(IGAM, WORK(KSCR6), ISYIJ2, ISYGAM,
867     &                 XLAMDP3, ISYML3, X2IJCB)
868      END IF
869
870*---------------------------------------------------------------------*
871*     Add the contribution to the result X1CJIB and X2CJIB vector:
872*     transform gamma to OCCUPIED
873*     Integrals sorted as I_cj,i,del
874*---------------------------------------------------------------------*
875      IF ( IOPT.EQ.1 ) THEN
876
877         IF ( LNEWTA ) THEN
878C           -------------------------
879C           add ([c^ j + c j^]|i del) to X1CJIB:
880C           -------------------------
881c***
882            ISYCJ3 = MULD2H(ISY1ALBE,ISYL13)
883
884            CALL CC_IAJB1(IGAM, WORK(KSCR5), ISYCJ3, ISYGAM,
885     &                    XLAMDP1, ISYML1, X1CJIB, .FALSE., IDUMMY)
886         END IF
887
888
889         IF ( LRELAX ) THEN
890C           ------------------------------
891C           add  ([c^ j + c j^]|i-bar del) to X2CJIB:
892C           ------------------------------
893            ISYCJ3 = MULD2H(ISY1ALBE,ISYL13)
894
895            CALL CC_IAJB1(IGAM, WORK(KSCR5), ISYCJ3, ISYGAM,
896     &                    XLAMDP2, ISYML2, X2CJIB, .FALSE., IDUMMY)
897         END IF
898
899         IF ( LDERIV .OR. LRELAX ) THEN
900C           ------------------------------------------------
901C           add ([c-bar j^ + c j^-bar + c^-bar j + c^ j-bar] | i del) to X2IABJ:
902C           ------------------------------------------------
903            ISYCJ4 = MULD2H(ISY1ALBE,ISYL14)
904
905            CALL CC_IAJB1(IGAM, WORK(KSCR7), ISYCJ4, ISYGAM,
906     &                    XLAMDP1, ISYML1, X2CJIB, .FALSE., IDUMMY)
907         END IF
908
909      END IF
910
911      RETURN
912      END
913*=====================================================================*
914*                 END OF SUBROUTINE CCIJCB                            *
915*=====================================================================*
916c /* deck cc_ijcb2 */
917*=====================================================================*
918      SUBROUTINE CC_IJCB2(IGAM, XIJG, ISYMIJ, ISYGAM,
919     &                    XLAMDA, ISYLAM, XIJCB  )
920*---------------------------------------------------------------------*
921*
922*   Purpose: transform (ij|gam del) to (ij|c del), with sorting
923*            I_cj,i,del
924*   Sonia, March 1999
925*---------------------------------------------------------------------*
926#if defined (IMPLICIT_NONE)
927      IMPLICIT NONE
928#else
929#  include "implicit.h"
930#endif
931#include "ccsdsym.h"
932#include "ccorb.h"
933#include "maxorb.h"
934#include "ccisao.h"
935
936#if defined (SYS_CRAY)
937      REAL XIJG(*), XIJCB(*), XLAMDA(*)
938#else
939      DOUBLE PRECISION XIJG(*), XIJCB(*), XLAMDA(*)
940#endif
941
942      INTEGER IGAM, ISYMIJ, ISYGAM, ISYLAM, ISYMC, KLAMD, KRES
943      INTEGER ISYMCJ, ISYMI, ISYMJ, KOFF5, NBASG
944      INTEGER KXIJ
945
946* transform integral batch:
947      ISYMC  = MULD2H(ISYGAM,ISYLAM)
948      G      = IGAM - IBAS(ISYGAM)
949      NBASG  = MAX(NBAS(ISYGAM),1)
950
951      DO ISYMJ = 1, NSYM
952         ISYMI  = MULD2H(ISYMIJ,ISYMJ)
953         ISYMCJ = MULD2H(ISYMC,ISYMJ)
954
955         DO J = 1, NRHF(ISYMJ)
956           DO I = 1, NRHF(ISYMI)
957
958              KLAMD = IGLMVI(ISYGAM,ISYMC) + G               !offs Lambda_gam,c
959              !offset I_ij,gam
960              KXIJ  = IMATIJ(ISYMI,ISYMJ)  + NRHF(ISYMI)*(J-1) + I
961              !offset I_cj,i
962              KRES  = IT2BCD(ISYMCJ,ISYMI) + NT1AM(ISYMCJ)*(I-1)
963     &                + IT1AM(ISYMC,ISYMJ) + NVIR(ISYMC)*(J-1) + 1
964
965              CALL DAXPY(NVIR(ISYMC),XIJG(KXIJ),XLAMDA(KLAMD),
966     &                                        NBASG,XIJCB(KRES),1)
967           END DO
968         END DO
969
970      END DO
971
972      RETURN
973      END
974*=====================================================================*
975*                 END OF SUBROUTINE CC_IJCB2                          *
976*=====================================================================*
977