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      subroutine cc_r12mkxajkl(kvajkl,cmo,work,lwork,fel)
20c---------------------------------------------------------------------
21c     purpose: calculate X_{aj}^{kl}, X = V,T necessary for the
22c     calculation of MP2-R12 contributions to the right hand side
23c     of Z-Vector Equations (cc_r12grad)
24c---------------------------------------------------------------------
25
26      implicit none
27#include "priunit.h"
28#include "dummy.h"
29#include "ccorb.h"
30#include "ccsdsym.h"
31#include "r12int.h"
32#if defined (SYS_CRAY)
33      real zero,one
34#else
35      double precision zero,one
36#endif
37      logical ldum,test,fel
38      parameter (test = .false.,zero = 0.0d0, one = 1.0d0)
39
40      integer  lwork,luvajkl,lwrk1,lwrk2,isymai
41      integer isymkl,isymij,isymmu,isyml,isymk,isymak,isymaj,icount7,
42     &     ntota,isymv,isymmuj,koffl,koffv,kvmujkl,
43     &     isymj,isyma,idxkl,xajkl,iglmrv(8,8)
44      integer kend1,kend2,koffc,isymi,ntklaj,yajkl
45      integer icou2,isymmua,ntotmu,nr1bas(8),ir1bas(8,8),icount2
46      integer noffro,ntota1,noffset,isym,nvir0(8)
47      integer idxij,idxji,kcmoa,idum,ntotgu,ntotj,lu43
48      integer nklaj(8),nr1vir(8),iglmro(8,8),icou3
49      character*10 fnelena
50      data fnelena/'DDR12xxxxx'/
51      real*8 cmo(*),work(*),kvajkl(*)
52c
53      call qenter('cc_r12mkxajkl')
54      fnelena = fnelena(1:5)//fnvajkl(6:10)
55
56      if (fnvajkl .eq. 'CCR12QAJKL' .and. fel .eqv. .true.) then
57         lu43 = -43
58         call dzero(cmo,nlamds)
59         if (R12HYB) then
60             call gpopen(lu43,'GUMAT.1','UNKNOWN',' ','UNFORMATTED',
61     &                      IDUM,LDUM)
62         else
63             call quit('Sorry, calculation of R12 corrections to'//
64     &                  ' first order properties not implemented'//
65     &                  ' with .NO HYB')
66         endif
67         REWIND(LU43)
68         READ(LU43) NTOTGU
69         READ(LU43) (cmo(I), I = 1, NTOTGU)
70         CALL GPCLOSE(LU43,'KEEP')
71      endif
72c
73      if (fnvajkl .eq. 'CCR12QAJKL' .and. fel .eqv. .false.) then
74         kend1 = 1
75         lwrk1 = lwork - kend1
76         lu43 = -43
77         call dzero(cmo,nlamds)
78         call gpopen(lu43,'SIRIFC','OLD',' ','UNFORMATTED',
79     &                      IDUMMY,.FALSE.)
80         REWIND(LU43)
81         CALL MOLLAB(LABEL,LU43,LUPRI)
82         READ(LU43)
83         READ(LU43)
84         READ(LU43) (cmo(I), I = 1, NLAMDS)
85         CALL GPCLOSE(LU43,'KEEP')
86         CALL CMO_REORDER(CMO,WORK(KEND1),LWRK1)
87
88      endif
89c
90
91      do isymak = 1, nsym
92        nr1bas(isymak) = 0
93        nr1vir(isymak) = 0
94        do isymk = 1, nsym
95           isyma = muld2h(isymak,isymk)
96           nvir0(isyma)    = norb1(isyma) - nrhfs(isyma)
97           nr1bas(isymak)  = nr1bas(isymak) +mbas1(isyma)*nrhf(isymk)
98           nr1vir(isymak)  = nr1vir(isymak) +nvir0(isyma)*nrhf(isymk)
99        end do
100      end do
101
102      do isymak = 1, nsym
103         nvajkl(isymak) = 0
104         icount7 = 0
105         icount2 = 0
106         do isymk = 1, nsym
107            isyma = muld2h(isymak,isymk)
108            ivajkl(isyma,isymk) = icount7
109            nvajkl(isymak) = nvajkl(isymak)+ nr1bas(isyma)*nmatij(isymk)
110            icount7 = icount7 + nr1bas(isyma)*nmatij(isymk)
111            ir1bas(isyma,isymk)  = icount2
112            icount2 = icount2 + nrhf(isymk)*mbas1(isyma)
113         end do
114      end do
115
116      ntklaj = 0
117      do isymkl = 1,nsym
118         isymaj = isymkl
119         nklaj(isymaj) = 0
120         do isymj = 1,nsym
121            isyma = muld2h(isymaj,isymj)
122            nklaj(isymaj) = nklaj(isymaj) + nmatij(isymj)*nr1vir(isyma)
123            ntklaj = ntklaj + nmatij(isymkl)*nr1vir(isymj)
124         enddo
125      enddo
126
127      do isymmua = 1,nsym
128         icou2 = 0
129         icou3 = 0
130         do isyma = 1,nsym
131            isymmu = muld2h(isymmua,isyma)
132            iglmrv(isymmu,isyma) = icou2
133            iglmro(isymmu,isyma) = icou3
134            icou2 = icou2 + nbas(isymmu)*(nvir(isyma)-nrhffr(isyma))
135            icou3 = icou3 + nbas(isymmu)*(norb1(isyma)-nrhffr(isyma))
136         enddo
137      enddo
138
139
140      noffro = 0
141      do isym = 1,nsym
142         noffro = noffro + nbas(isym)*nrhf(isym)
143      enddo
144
145
146      kend1 = 1
147      yajkl = kend1
148      kcmoa =  yajkl + ntklaj
149      kend2  = kcmoa + nlamds
150      lwrk2  = lwork - kend2
151
152      if (lwrk2 .lt.0) then
153         call quit('Insufficient work space in cc_r12mkxajkl')
154      end if
155c
156c
157      call dzero(work(yajkl),ntklaj)
158      xajkl = yajkl
159      do isymkl = 1, nsym
160         isymaj  = isymkl
161         isymmuj = isymkl
162         do isyml =1, nsym
163            isymk = muld2h(isymkl,isyml)
164            do k = 1, nrhf(isymk)
165               do l = 1, nrhf(isyml)
166                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
167                  noffset = 0
168                  do isymmu =1, nsym
169                     isyma = muld2h(1,isymmu)
170                     isymj  = muld2h(isymmuj,isymmu)
171                     nvir0(isyma) = norb1(isyma) - nrhfs(isyma)
172                     if (fel) then
173                        koffc = 1 + iglmro(isymmu,isyma)
174     &                          +nrhf(isymmu)*nbas(isymmu)
175                     else
176                        noffset = noffset+ nbas(isymmu)*nrhffr(isymmu)
177                        koffc = 1 + nrhf(isymmu)*nbas(isymmu)
178     &                        +iglmrv(isymmu,isyma)+noffset
179     &                        +noffro
180                     endif
181                     kvmujkl = 1+ivajkl(isymmuj,isymkl)+
182     &                         nr1bas(isymmuj)*(idxkl-1)
183     &                         +ir1bas(isymmu,isymj)
184                     ntotmu = max(1,mbas1(isymmu))
185                     ntota1 = max(1,nbas(isymmu))
186                     ntota = max(1,nvir0(isyma))
187                     ntotj = max(1,nrhf(isymj))
188c
189                     call dgemm('T','N',nvir0(isyma),
190     &                 nrhf(isymj),mbas1(isymmu),one,cmo(koffc),
191     &                 ntota1,kvajkl(kvmujkl),ntotmu,zero,
192     &                 work(xajkl),ntota)
193
194                     xajkl = xajkl + nvir0(isyma)*nrhf(isymj)
195                  end do
196               end do
197            end do
198         end do
199      end do
200
201      if (r12cbs) call dscal(nvajkl(1),1.00d0,work(yajkl),1)
202
203
204      luvajkl = -1
205      call gpopen(luvajkl,fnelena,'unknown',' ','unformatted',
206     &               idummy,.false.)
207      rewind(luvajkl)
208      write(luvajkl) (work(yajkl+i-1), i = 1,ntklaj)
209      call gpclose(luvajkl,'KEEP')
210
211
212      call qexit('cc_r12mkxajkl')
213      end
214C=====================================================================*
215C  /* Deck cc_r12grad */
216      subroutine cc_r12grad(ipdd,orb,v2am,work,lwork)
217C     Calculate MP2-R12 contributions to the right hand side of Z-Vector
218C     Equations
219      implicit none
220#include "ccsdsym.h"
221#include "ccorb.h"
222#include "dummy.h"
223#include "cbirea.h"
224#include "r12int.h"
225#include "priunit.h"
226#include "ccfro.h"
227      logical locdbg
228      logical lauxd,ldum
229      parameter (locdbg = .false.)
230      character*10 CCR12YAJIJ
231      integer kend1,kccr12,kend2,lwrk1,lwork,nrhftria,isymij
232      integer lu43,idum,n2,lwrk2,ipdd
233      integer krsing,krtrip,kend3,lwrk3,nr1bas(8)
234      integer isymkl,isyml,isymk,idxij,idxji,idxklij,idxklji
235      integer isymai,isymv,isymi,isymj,idxkl,kctil,ntklaj
236      integer kbajkl,kbijal,kvajkl,kvijal,isymaj,isyma,luvajkl
237      integer lcvai,koffc,nvir0(8),isymmu,ntota,ntotmu,kcvai
238      integer ir1bas(8,8),isymmuj
239      integer dklij,ntotk,lcvia,kvmujkl,ntoti
240      integer koffd,bdajij,lccr12
241      integer dklijz,kajij,kcvia,ntotj,imataj(8,8),icou3
242      integer nr1orb(8),nr1xbas(8),nrgkl,n2bst1(8),ir1orb(8,8)
243      integer ir2bas(8,8),ir2xbas(8,8),irgkl,iaodis1(8,8),nr2bas(8)
244      integer ir1xbas(8,8),nijkl(8)
245      integer koffe,kctil2,idxijkl,lctil,ncvai
246      integer res,kscr1
247
248#if defined (SYS_CRAY)
249      real work(lwork),one,two,zero
250      real orb(*),v2am(*)
251#else
252      double precision work(lwork),one,zero,two,orb(*),v2am(*)
253#endif
254      parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0)
255
256
257      call qenter('cc_r12grad')
258
259
260      nrhftria = nrhft*(nrhft+1)/2
261      n2 = nrhftria * nrhftria
262
263      kccr12  = 1
264      krsing  = kccr12  + ngamsq(1)
265      krtrip  = krsing  + n2
266      kctil   = krtrip  + n2
267      kctil2  = kctil   + ngamsq(1)
268      kend1   = kctil2  + ngamsq(1)
269      lwrk1   = lwork   - kend1
270
271      call dzero(work(kccr12),ngamsq(1))
272      call dzero(work(krsing),n2)
273      call dzero(work(krtrip),n2)
274      call dzero(work(kctil),ngamsq(1))
275      call dzero(work(kctil2),ngamsq(1))
276
277      isymv = 1
278
279      if (lwrk1 .lt. 0) then
280         call quit('Insufficient work space in cc_r12grad')
281      END IF
282
283c     ----------------------------------------
284c     read R12 amplitudes
285c     ----------------------------------------
286
287
288      lu43 = -43
289      call gpopen(lu43,'CCR12_C','unknown',' ','formatted',
290     &                   idummy,.false.)
291      read(lu43,'(4e30.20)') (work(krsing+i), i = 0, n2-1 )
292      read(lu43,'(4e30.20)') (work(krtrip+i), i = 0, n2-1 )
293      call gpclose(lu43,'KEEP')
294
295      call dscal(nrhftria*nrhftria,0.5d0,work(krsing),1)
296      call dscal(nrhftria*nrhftria,0.5d0,work(krtrip),1)
297
298      call cc_r12vpack(work(kccr12),isymv,work(krsing),work(krtrip),
299     &     nrhf,nrhfa,nijkl)
300
301c     ---------------------------------
302C     get c_tilde = 2*c_ij^kl - c_ji^kl
303c     ---------------------------------
304
305      do isymkl = 1, nsym
306         isymij = muld2h(isymkl,isymv)
307         do isyml =1, nsym
308            isymk = muld2h(isymkl,isyml)
309            do k = 1, nrhf(isymk)
310               do l = 1, nrhf(isyml)
311                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
312                  do isymi = 1, nsym
313                     isymj = muld2h(isymij,isymi)
314                     do j = 1, nrhf(isymj)
315                        do i = 1, nrhf(isymi)
316                           idxij = imatij(isymi,isymj)+
317     &                          nrhf(isymi)*(j-1) + i
318                           idxji = imatij(isymj,isymi)+
319     &                          nrhf(isymj)*(i-1) + j
320                           idxklij =  igamsq(isymkl,isymij)+
321     &                          nmatij(isymkl)*(idxij-1)+idxkl
322                           idxklji =igamsq(isymkl,isymij)+
323     &                          nmatij(isymkl)*(idxji-1)+idxkl
324                           work(kctil-1+idxklij) =
325     &                          two*work(kccr12-1+idxklij) -
326     &                          work(kccr12-1+idxklji)
327                        end do
328                     end do
329
330                  end do
331               end do
332            end do
333         end do
334      end do
335
336      if (locdbg) then
337         write(lupri,*) 'resorted R12 amplitudes (ctilde):'
338         do isymkl = 1, nsym
339            isymij = muld2h(isymkl,1)
340            write(lupri,*) 'symmetry block ',isymkl,isymij
341            call output(work(kctil+igamsq(isymkl,isymij)),
342     &           1,nmatij(isymkl),1,nmatij(isymij),
343     &            nmatij(isymkl),nmatij(isymij),1,lupri)
344         end do
345      endif
346
347c     ---------------------------------
348C     Transpone c_tilde
349c     ---------------------------------
350
351      do isymkl = 1, nsym
352         isymij = muld2h(isymkl,isymv)
353         do isyml =1, nsym
354            isymk = muld2h(isymkl,isyml)
355            do k = 1, nrhf(isymk)
356               do l = 1, nrhf(isyml)
357                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
358                  do isymi = 1, nsym
359                     isymj = muld2h(isymij,isymi)
360                     do j = 1, nrhf(isymj)
361                        do i = 1, nrhf(isymi)
362                           idxij = imatij(isymi,isymj)+
363     &                          nrhf(isymi)*(j-1) + i
364                           idxji = imatij(isymj,isymi)+
365     &                          nrhf(isymj)*(i-1) + j
366                           idxklij =  igamsq(isymij,isymkl)+
367     &                          nmatij(isymkl)*(idxij-1)+idxkl
368                           idxklji =igamsq(isymij,isymkl)+
369     &                          nmatij(isymkl)*(idxji-1)+idxkl
370                           idxijkl = igamsq(isymkl,isymij)+
371     &                          nmatij(isymij)*(idxkl-1) +idxij
372                           work(kctil2-1+idxklij) =
373     &                          work(kctil-1+idxijkl)
374                        end do
375                     end do
376
377                  end do
378               end do
379            end do
380         end do
381      end do
382
383
384c
385C---------------------------------------------
386C     Calculate some offsets and dimensions
387C---------------------------------------------
388      ntklaj = 0
389      do isymkl = 1, nsym
390         isymaj  = isymkl
391         do isyml =1, nsym
392            isymk = muld2h(isymkl,isyml)
393            do k = 1, nrhf(isymk)
394               do l = 1, nrhf(isyml)
395                  do isyma =1, nsym
396                     isymj  = muld2h(isymaj,isyma)
397                     isymi  = isyma
398                     ntklaj = ntklaj + nrhf(isymj) *
399     &                        (norb1(isyma) - nrhfs(isyma))
400                  enddo
401               enddo
402            enddo
403         enddo
404      enddo
405
406
407      icou3 = 0
408      do isymai = 1,nsym
409         icou3 = 0
410         do isyma = 1,nsym
411            isymi = muld2h(isymai,isyma)
412            imataj(isyma,isymi)= icou3
413            icou3 = icou3+nrhf(isymi)*(norb1(isyma) - nrhfs(isyma))
414         enddo
415      enddo
416
417      do isymaj = 1,nsym
418         isymij = isymaj
419        ncvai = 0
420         do isyma = 1,nsym
421            isymj = muld2h(isymaj,isyma)
422            isymi = muld2h(isymij,isymj)
423            ncvai = ncvai + (norb1(isyma) - nrhfs(isyma))*nrhf(isymi)
424         enddo
425      enddo
426
427
428      kbajkl  = kend1
429      kbijal  = kbajkl  + ntklaj
430      kvajkl  = kbijal  + ntklaj
431      kvijal  = kvajkl  + ntklaj
432      kcvai   = kvijal  + ntklaj
433      kcvia   = kcvai   + ncvai
434      kajij   = kcvia   + ncvai
435      dklij   = kajij   + ncvai
436      res     = dklij   + ngamsq(1)
437      kend2   = res     + ncvai
438      lwrk2   = lwork   - kend2
439
440c     ----------------------------------------------------
441c     Read V_{aj}^{kl},V_{kl}^{aj},B_{aj}^{kl},B_{kl}^{aj}
442c               and transform it with Ctilde
443c     -----------------------------------------------------
444
445      if (lwrk2 .lt.0) then
446         call quit('Insufficient work space in cc_r12grad')
447      end if
448
449      call dzero(work(kbajkl),ntklaj)
450      call dzero(work(kbijal),ntklaj)
451      luvajkl = -1
452      call gpopen(luvajkl,'DDR12BAJKL','unknown',' ','unformatted',
453     &                   idummy,.false.)
454      rewind(luvajkl)
455      read(luvajkl) (work(kbajkl+i-1), i = 1,ntklaj)
456      if (ipdd .eq. 5) then
457          call gpclose(luvajkl,'DELETE')
458      else
459          call gpclose(luvajkl,'KEEP')
460      endif
461
462
463      call gpopen(luvajkl,'DDR12BIJAL','unknown',' ','unformatted',
464     &                   idummy,.false.)
465      rewind(luvajkl)
466      read(luvajkl) (work(kbijal+i-1), i = 1,ntklaj)
467      if (ipdd .eq. 5) then
468          call gpclose(luvajkl,'DELETE')
469      else
470          call gpclose(luvajkl,'KEEP')
471      endif
472c
473      call daxpy(ntklaj,one,work(kbajkl),1,work(kbijal),1)
474
475      call dzero(work(kvajkl),ntklaj)
476      call dzero(work(kvijal),ntklaj)
477      call gpopen(luvajkl,'DDR12VAJKL','unknown',' ','unformatted',
478     &                   idummy,.false.)
479      rewind(luvajkl)
480      read(luvajkl) (work(kvajkl+i-1), i = 1,ntklaj)
481      if (ipdd .eq. 5) then
482          call gpclose(luvajkl,'DELETE')
483      else
484          call gpclose(luvajkl,'KEEP')
485      endif
486
487      call gpopen(luvajkl,'DDR12VIJAL','unknown',' ','unformatted',
488     &                   idummy,.false.)
489      rewind(luvajkl)
490      read(luvajkl) (work(kvijal+i-1), i = 1,ntklaj)
491      if (ipdd .eq. 5) then
492          call gpclose(luvajkl,'DELETE')
493      else
494          call gpclose(luvajkl,'KEEP')
495      endif
496
497
498
499c     Calculate d_{kl}^{ij} = c_{kl}^{mn}*ctilde_{ij}^{mn}
500
501      call dzero(work(dklij),ngamsq(1))
502      do isymkl = 1,nsym
503         isymij  = isymkl
504         lccr12 = kccr12 + igamsq(isymij,isymkl)
505         lctil  = kctil + igamsq(isymij,isymkl)
506         dklijz = dklij + igamsq(isymij,isymkl)
507         ntoti = max(1,nmatij(isymij))
508         ntotk = max(1,nmatij(isymkl))
509         call dgemm('T','N',nmatij(isymij),nmatij(isymkl),
510     &              nmatij(isymkl),one,work(lccr12),
511     &              ntotk,work(lctil),ntoti,zero,
512     &              work(dklijz),ntotk)
513      end do
514
515
516      call dzero(work(kcvai),ncvai)
517      lcvai = kcvai
518      call dzero(work(kcvia),ncvai)
519      lcvia = kcvia
520      call dzero(work(kajij),ncvai)
521      bdajij = kajij
522
523c     Calculate V_{aj}^{kl}*ctilde_{kl}^{ij}
524
525
526      do isymkl = 1, nsym
527         isymaj  = isymkl
528         isymij  = isymkl
529         do isyml =1, nsym
530            isymk = muld2h(isymkl,isyml)
531           do k = 1, nrhf(isymk)
532             do l = 1, nrhf(isyml)
533                idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
534                 do isyma =1, nsym
535                    isymj  = muld2h(isymaj,isyma)
536                    isymi  = muld2h(isymij,isymj)
537                    nvir0(isyma)  = norb1(isyma) - nrhfs(isyma)
538                    koffc  = kctil + igamsq(isymij,isymkl)
539     &                 + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
540                    ntotmu = max(1,nvir0(isyma))
541                    ntotj = max(1,nrhf(isymj))
542                    ntoti = max(1,nrhf(isymi))
543                    lcvai  = kcvai + imataj(isyma,isymi)
544                    call dgemm('N','T',nvir0(isyma),nrhf(isymi),
545     &                    nrhf(isymj),one,work(kvajkl),
546     &                    ntotmu,work(koffc),ntoti,one,
547     &                    work(lcvai),ntotmu)
548
549                    kvajkl = kvajkl + nvir0(isyma)*nrhf(isymj)
550
551c     Calculate V_{kl}^{aj}*ctilde_{ij}^{kl}
552
553                    lcvia  = kcvia + imataj(isyma,isymi)
554                    koffe = kctil2+ igamsq(isymkl,isymij)
555     &              +  nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
556                    call dgemm('N','T',nvir0(isyma),nrhf(isymi),
557     &                         nrhf(isymj),one,work(kvijal),
558     &                         ntotmu,work(koffe),ntoti,one,
559     &                         work(lcvia),ntotmu)
560                    kvijal = kvijal + nvir0(isyma)*nrhf(isymj)
561c
562c     Calculate B_{aj}^{kl}*d_{kl}^{ij}
563c
564                    koffd = dklij + igamsq(isymij,isymkl)
565     &              +  nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
566                    bdajij  = kajij + imataj(isyma,isymi)
567                    call dgemm('N','T',nvir0(isyma),nrhf(isymi),
568     &                         nrhf(isymj),one,work(kbijal),
569     &                         ntotmu,work(koffd),ntoti,one,
570     &                         work(bdajij),ntotmu)
571                    kbijal  = kbijal + nvir0(isyma)*nrhf(isymj)
572                end do
573              end do
574            end do
575         end do
576      end do
577      call dscal(ncvai,two,work(kcvai),1)
578      call dscal(ncvai,two,work(kcvia),1)
579      call daxpy(ncvai,one,work(kcvai),1,work(kcvia),1)
580      call daxpy(ncvai,one,work(kcvia),1,work(kajij),1)
581      call dscal(ncvai,two,work(kajij),1)
582      if (.not. lmulbs) call dscal(ncvai,-one,work(kajij),1)
583
584      if (locdbg) then
585         write(lupri,*) 'Matrixelemente MP2-R12 Gradient',ipdd
586         do isymaj = 1,nsym
587            do isyma = 1,nsym
588               isymi =isyma
589               write(lupri,*) 'symmetry block',imataj(isyma,isymi)
590               call output(work(kajij+imataj(isyma,isymi)),
591     &              1,nvir0(isyma),1,nrhf(isymi),
592     &              nvir0(isyma),nrhf(isymi),1,lupri)
593
594            enddo
595         enddo
596      end if
597
598      if (ipdd .eq. 2) then
599        luvajkl = -1
600        call gpopen(luvajkl,'CCR12YAJIJ','unknown',' ','unformatted',
601     &               idummy,.false.)
602        write(luvajkl) (work(kajij+i-1), i = 1,ncvai)
603        call gpclose(luvajkl,'KEEP')
604      else if (ipdd .eq. 3) then
605        call dzero(work(res),ncvai)
606        call r12gradap(ipdd,work(res),orb,work(kend2),v2am,lwrk2)
607        call dscal(ncvai,2.0D0,work(res),1)
608        call daxpy(ncvai,one,work(kajij),1,work(res),1)
609        luvajkl = -1
610        call gpopen(luvajkl,'CCR12ZAJIJ','unknown',' ','unformatted',
611     &               idummy,.false.)
612        write(luvajkl) (work(res+i-1), i = 1,ncvai)
613        call gpclose(luvajkl,'KEEP')
614      else if (ipdd .eq. 5) then
615        call dzero(work(res),ncvai)
616        call r12gradap(ipdd,work(res),orb,work(kend2),v2am,lwrk2)
617        call dscal(ncvai,2.00D0,work(res),1)
618        call daxpy(ncvai,one,work(kajij),1,work(res),1)
619        luvajkl = -1
620        call gpopen(luvajkl,'CCR12XAJIJ','unknown',' ','unformatted',
621     &               idummy,.false.)
622        write(luvajkl) (work(res+i-1), i = 1,ncvai)
623        call gpclose(luvajkl,'KEEP')
624
625      endif
626c
627c
628      call qexit('cc_r12grad')
629      end
630c------------------------------------------------------------------------
631C=====================================================================*
632C  /* Deck r12gradap */
633      subroutine r12gradap(ipdd,res,orb,work,v2am,lwork)
634      implicit none
635#include "ccsdsym.h"
636#include "ccorb.h"
637#include "dummy.h"
638#include "r12int.h"
639#include "priunit.h"
640#include "ccfro.h"
641#include "ccsdinp.h"
642      integer kccr12,krsing,krtrip,kctil,kend1,lwrk1,nrhftria,n2
643      integer lwork,lu43,idxkl,isymv,isymi,isymj,isymk,idxij
644      integer isyml,isymkl,isymij,idxji,idxklij,idxijkl,idxklji
645      integer kxajkl,lunit,luvajkl,ntklaj,isyma,isymaj,kend2
646      integer nmatan(8),isyman,isymn,ntotij,ntotan,nvir0(8)
647      integer lctil,klijan,ntota,ntoti,ntotn,lcvai,isymai
648      integer kijan,dklij,dklijz,ntotk,koffc,lccr12,koffd
649      integer idxijlk,idxlk,isymin,kbajkl,ncvai,imatan(8,8),icou3
650      integer lwrk2,lctil1,klij,kcotil,kdotil,koccr,koctil
651      integer idxlkji,idxlkij,idxjilk,krklaj,krajkl,ntaj
652      integer kxsing,kxtrip,kxxr12,klijmn,lxxr12,ij,al,idx
653      integer isym,nvir0t,blajkl,klijz,krajklz,klijanz,kxajklz
654      integer ipdd,kxtil,koffk,koffl,koffj,koffi,isymmu,isymmuj
655      integer icou4,imataj(8,8),blajklf,kl,aj,ajkl,nrhfst
656      integer imatmn(8,8),icou6,nmatkl1(8),imatkl1(8,8),icou5
657      integer imatka(8,8),icou7,nmatak(8),nijkl(8)
658#if defined (SYS_CRAY)
659      real res(*),work(*),orb(*),v2am(*),one,two,zero
660#else
661      double precision work(*),orb(*),v2am(*),one,zero,two
662      double precision res(*)
663#endif
664      parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0)
665
666      call qenter('r12gradap')
667
668      isymv = 1
669
670      nrhftria = nrhft*(nrhft+1)/2
671      n2 = nrhftria * nrhftria
672
673      nrhfst = 0
674      do isymi = 1,nsym
675         nrhfst = nrhfst + nrhfs(isymi)
676      enddo
677
678      kccr12  = 1
679      krsing  = kccr12  + ngamsq(1)
680      krtrip  = krsing  + n2
681      kctil   = krtrip  + n2
682      kend1   = kctil  + ngamsq(1)
683      lwrk1   = lwork   - kend1
684
685      call dzero(work(krsing),2*n2)
686      call dzero(work(kccr12),ngamsq(1))
687      call dzero(work(kctil),ngamsq(1))
688
689      if (lwrk1 .lt. 0) then
690         call quit('Insufficient work space in r12gradap')
691      end if
692
693c     ----------------------------------------
694c     read R12 amplitudes
695c     ----------------------------------------
696      lu43 = -43
697      call gpopen(lu43,'CCR12_C','unknown',' ','formatted',
698     &                   idummy,.false.)
699      read(lu43,'(4e30.20)') (work(krsing+i), i = 0, n2-1 )
700      read(lu43,'(4e30.20)') (work(krtrip+i), i = 0, n2-1 )
701      call gpclose(lu43,'KEEP')
702
703
704      call dscal(nrhftria*nrhftria,0.5d0,work(krsing),1)
705      call dscal(nrhftria*nrhftria,0.5d0,work(krtrip),1)
706
707      call cc_r12vpack(work(kccr12),isymv,work(krsing),work(krtrip),
708     &     nrhf,nrhfa,nijkl)
709
710c     ---------------------------------
711C     get c_tilde = 2*c_ij^kl - c_ji^kl
712c     ---------------------------------
713
714      do isymkl = 1, nsym
715         isymij = muld2h(isymkl,isymv)
716         do isyml =1, nsym
717            isymk = muld2h(isymkl,isyml)
718            do k = 1, nrhf(isymk)
719               do l = 1, nrhf(isyml)
720                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
721                  do isymi = 1, nsym
722                     isymj = muld2h(isymij,isymi)
723                     do j = 1, nrhf(isymj)
724                        do i = 1, nrhf(isymi)
725                           idxij = imatij(isymi,isymj)+
726     &                          nrhf(isymi)*(j-1) + i
727                           idxji = imatij(isymj,isymi)+
728     &                          nrhf(isymj)*(i-1) + j
729                           idxklij =  igamsq(isymij,isymkl)+
730     &                          nmatij(isymkl)*(idxij-1)+idxkl
731                           idxklji =igamsq(isymkl,isymij)+
732     &                          nmatij(isymkl)*(idxji-1)+idxkl
733                           idxijkl = igamsq(isymkl,isymij)+
734     &                          nmatij(isymij)*(idxkl-1) +idxij
735                           work(kctil-1+idxklij) =
736     &                          two*work(kccr12-1+idxklij) -
737     &                          work(kccr12-1+idxklji)
738                        end do
739                     end do
740
741                  end do
742               end do
743            end do
744         end do
745      end do
746
747
748
749      ntklaj = 0
750      do isymkl = 1, nsym
751         isymaj  = isymkl
752         do isyml =1, nsym
753            isymk = muld2h(isymkl,isyml)
754            do k = 1, nrhf(isymk)
755               do l = 1, nrhf(isyml)
756                  do isyma =1, nsym
757                     isymj  = muld2h(isymaj,isyma)
758                     isymi  = isyma
759                     ntklaj = ntklaj + nrhf(isymj) *
760     &                        (norb1(isyma) - nrhfs(isyma))
761                  enddo
762               enddo
763            enddo
764         enddo
765      enddo
766
767      do isyman = 1,nsym
768         icou3 = 0
769         icou4 = 0
770         icou5 = 0
771         do isyma = 1,nsym
772            isymn = muld2h(isyman,isyma)
773            nvir0(isyma)  = norb1(isyma) - nrhfs(isyma)
774            icou3 = icou3 + nvir0(isyma)*nrhfs(isymn)
775            icou4 = icou4 + nrhfs(isyma)*nrhfs(isymn)
776            icou5 = icou5 + nvir0(isyma)*nrhf(isymn)
777          enddo
778          nmatan(isyman) = icou3
779          nmatkl1(isyman) = icou4
780          nmatak(isyman) = icou5
781      enddo
782
783      icou3 = 0
784      do isyman = 1,nsym
785         icou3 = 0
786         icou4 = 0
787         icou5 = 0
788         icou6 = 0
789         icou7 = 0
790         do isyma = 1,nsym
791            isymn = muld2h(isyman,isyma)
792            imatan(isyma,isymn)= icou3
793            imataj(isymn,isyma)= icou4
794            imatkl1(isyma,isymn)= icou5
795            imatmn(isymn,isyma)= icou6
796            imatka(isymn,isyma)= icou7
797            icou3 = icou3+nrhf(isymn)*(norb1(isyma) - nrhfs(isyma))
798            icou4 = icou4+nrhfs(isymn)*(norb1(isyma) - nrhfs(isyma))
799            icou5 = icou5+nmatan(isymn)* nmatkl1(isyma)
800            icou6 = icou6+nrhfs(isyma)*nrhfs(isymn)
801            icou7 = icou7+nmatak(isyma)*nmatij(isymn)
802         enddo
803      enddo
804      do isyman = 1,nsym
805         do isyma = 1,nsym
806            isymn = muld2h(isyma,isyman)
807         enddo
808      enddo
809
810      do isymaj = 1,nsym
811         isymij = isymaj
812         ncvai = 0
813         do isyma = 1,nsym
814            isymj = muld2h(isymaj,isyma)
815            isymi = muld2h(isymij,isymj)
816            ncvai = ncvai + (norb1(isyma) - nrhfs(isyma))*nrhf(isymi)
817         enddo
818      enddo
819
820      ntaj = 0
821      do isymkl = 1,nsym
822         isymij = isymkl
823         do isymj = 1,nsym
824            isyma = muld2h(isymij,isymj)
825             ntaj = ntaj + nrhf(isymj)
826     &            *(norb1(isyma) - nrhfs(isyma))
827         enddo
828      enddo
829
830      nvir0t = 0
831      do isym = 1, nsym
832         nvir0(isym) = norb1(isym) - nrhfs(isym)
833         nvir0t      = nvir0t + nvir0(isym)
834      end do
835
836
837      kxajkl  = kend1
838      klijan  = kxajkl + ntklaj
839      dklij   = klijan + ncvai
840      koctil  = dklij  + ngamsq(1)
841      kdotil  = koctil + ngamsq(1)
842      klij    = kdotil + ngamsq(1)
843      krajkl  = klij   + ngamsq(1)
844      kxxr12  = krajkl + nvir0t*nrhfst**3
845      kxsing  = kxxr12 + ngamsq(1)
846      kxtrip  = kxsing + n2
847      klijmn  = kxtrip + n2
848      blajkl  = klijmn + ngamsq(1)
849      kxtil   = blajkl + nvir0t*nrhfst**3
850      blajklf = kxtil  + ngamsq(1)
851      kend2   = blajklf + ntklaj
852      lwrk2   = lwork  - kend2
853
854
855
856
857      call dzero(work(koctil),ngamsq(1))
858       do isymkl = 1, nsym
859         isymij = muld2h(isymkl,isymv)
860         do isyml =1, nsym
861            isymk = muld2h(isymkl,isyml)
862            do k = 1, nrhf(isymk)
863               do l = 1, nrhf(isyml)
864                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
865                  idxlk=imatij(isyml,isymk)+nrhf(isyml)*(k-1)+l
866                  do isymi = 1, nsym
867                     isymj = muld2h(isymij,isymi)
868                     do j = 1, nrhf(isymj)
869                        koffj = irhf(isymj)+j
870                        do i = 1, nrhf(isymi)
871                           ij = (i-1)*nrhf(isymj)+j
872                           koffi = irhf(isymi)+i
873                           idxij = imatij(isymi,isymj)+
874     &                          nrhf(isymi)*(j-1) + i
875                           idxji = imatij(isymj,isymi)+
876     &                          nrhf(isymj)*(i-1) + j
877                           idxklij =  igamsq(isymij,isymkl)+
878     &                          nmatij(isymkl)*(idxij-1)+idxkl
879                           idxklji =igamsq(isymkl,isymij)+
880     &                          nmatij(isymkl)*(idxji-1)+idxkl
881                           idxijkl = igamsq(isymkl,isymij)+
882     &                          nmatij(isymij)*(idxkl-1) +idxij
883                           idxijlk = igamsq(isymij,isymkl)+
884     &                          nmatij(isymij)*(idxlk-1) +idxij
885                           work(koctil-1+idxijlk) =
886     &                           work(koctil-1+idxijlk) -two*
887     &                          work(kctil-1+idxijlk) * orb(koffi)
888                           work(koctil-1+idxijlk) =
889     &                          work(koctil-1+idxijlk) -two*
890     &                          work(kctil-1+idxijlk) * orb(koffj)
891                        end do
892                     end do
893
894                  end do
895               end do
896            end do
897         end do
898      end do
899
900c     Calculate d_{kl}^{ij} = c_{kl}^{mn}*ctilde_{ij}^{mn}
901
902      call dzero(work(dklij),ngamsq(1))
903      call dzero(work(klij),ngamsq(1))
904      do isymkl = 1,nsym
905         isymij  = isymkl
906         lccr12 = kccr12 + igamsq(isymij,isymkl)
907         lctil  = koctil + igamsq(isymij,isymkl)
908         lctil1  = kctil + igamsq(isymij,isymkl)
909         dklijz = dklij + igamsq(isymij,isymkl)
910         klijz = klij + igamsq(isymij,isymkl)
911         ntoti = max(1,nmatij(isymij))
912         ntotk = max(1,nmatij(isymkl))
913         call dgemm('T','N',nmatij(isymij),nmatij(isymkl),
914     &              nmatij(isymkl),one,work(lccr12),
915     &              ntotk,work(lctil),ntoti,zero,
916     &              work(dklijz),ntotk)
917         call dgemm('T','N',nmatij(isymij),nmatij(isymkl),
918     &              nmatij(isymkl),one,work(lccr12),
919     &              ntotk,work(lctil1),ntoti,zero,
920     &              work(klijz),ntotk)
921      end do
922
923
924
925      call dzero(work(kdotil),ngamsq(1))
926      do isymkl = 1, nsym
927        isymij = muld2h(isymkl,isymv)
928         do isyml =1, nsym
929            isymk = muld2h(isymkl,isyml)
930            do k = 1, nrhf(isymk)
931               koffk = irhf(isymk)+k
932               do l = 1, nrhf(isyml)
933                  koffl = irhf(isyml)+l
934                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
935                  idxlk=imatij(isyml,isymk)+nrhf(isyml)*(k-1)+l
936                  do isymi = 1, nsym
937                     isymj = muld2h(isymij,isymi)
938                     do j = 1, nrhf(isymj)
939                        koffj = irhf(isymj)+j
940                        do i = 1, nrhf(isymi)
941                           koffi = irhf(isymi)+i
942                           idxij = imatij(isymi,isymj)+
943     &                          nrhf(isymi)*(j-1) + i
944                           idxji = imatij(isymj,isymi)+
945     &                          nrhf(isymj)*(i-1) + j
946                           idxklij =  igamsq(isymij,isymkl)+
947     &                          nmatij(isymkl)*(idxij-1)+idxkl
948                           idxklji =igamsq(isymkl,isymij)+
949     &                          nmatij(isymkl)*(idxji-1)+idxkl
950                           idxijkl = igamsq(isymkl,isymij)+
951     &                          nmatij(isymij)*(idxkl-1) +idxij
952                           idxijlk = igamsq(isymkl,isymij)+
953     &                          nmatij(isymij)*(idxlk-1) +idxij
954                           work(kdotil-1+idxijlk) =
955     &                          work(kdotil-1+idxijlk) +
956     &                          work(klij-1+idxijlk) * orb(koffi)
957                           work(kdotil-1+idxijlk) =
958     &                          work(kdotil-1+idxijlk) +
959     &                          work(klij-1+idxijlk) * orb(koffj)
960                           work(kdotil-1+idxijlk) =
961     &                          work(kdotil-1+idxijlk) +
962     &                          work(klij-1+idxijlk) * orb(koffk)
963                           work(kdotil-1+idxijlk) =
964     &                          work(kdotil-1+idxijlk) +
965     &                          work(klij-1+idxijlk) * orb(koffl)
966                        end do
967                     end do
968
969                  end do
970               end do
971            end do
972         end do
973      end do
974
975      call daxpy(ngamsq(1),one,work(dklij),1,work(kdotil),1)
976
977
978      luvajkl = -1
979      call dzero(work(kxajkl),ntklaj)
980      call gpopen(luvajkl,'DDR12XAJKL','unknown',' ','unformatted',
981     &                   idummy,.false.)
982      rewind(luvajkl)
983      read(luvajkl) (work(kxajkl+i-1), i = 1,ntklaj)
984      call gpclose(luvajkl,'KEEP')
985
986
987
988      call dzero(work(krajkl),nvir0t*nrhfst**3)
989      luvajkl = -1
990      call gpopen(luvajkl,'AUXQA12','unknown',' ','formatted',
991     &                   idummy,.false.)
992      rewind(luvajkl)
993      read(luvajkl,'(4E30.20)') (work(krajkl+i-1), i=1,nvir0t*nrhfst**3)
994      call gpclose(luvajkl,'KEEP')
995
996      call dzero(work(blajkl),ntklaj)
997      call cc_r12xsort(work(krajkl),work(blajkl),nvir0)
998
999      if (froimp) then
1000         call dzero(work(blajklf),ntklaj)
1001         idx = 0
1002         do isymkl = 1, nsym
1003            isymaj  = isymkl
1004            isymij  = isymkl
1005            do isyml =1, nsym
1006               isymk = muld2h(isymkl,isyml)
1007               do k = nrhffr(isymk)+1, nrhfs(isymk)
1008                  do l = nrhffr(isyml)+1, nrhfs(isyml)
1009                     do isyma =1, nsym
1010                        isymj  = muld2h(isymaj,isyma)
1011                        nvir0(isyma)  = norb1(isyma) - nrhfs(isyma)
1012                        do j = nrhffr(isymj)+1, nrhfs(isymj)
1013                           do a = 1, nvir0(isyma)
1014                              idx = idx + 1
1015                              kl =imatmn(isymk,isyml)+
1016     &                            (k-1)*nrhfs(isyml) + l
1017                              aj = imataj(isymj,isyma)+
1018     &                            (j-1)*nvir0(isyma) + a
1019                              ajkl = imatkl1(isymaj,isymkl)+
1020     &                          aj +(kl-1)*nmatan(isymaj)
1021                              work(blajklf+idx-1) = work(blajkl+ajkl-1)
1022                           enddo
1023                        enddo
1024                     enddo
1025                  enddo
1026               enddo
1027            enddo
1028         enddo
1029      endif
1030
1031
1032      if (froimp) then
1033         call daxpy(ntklaj,one,work(blajklf),1,work(kxajkl),1)
1034      else
1035         call daxpy(ntklaj,one,work(blajkl),1,work(kxajkl),1)
1036      endif
1037
1038      call dzero(work(klijan),ncvai)
1039      do isymkl = 1, nsym
1040         isyman  = isymkl
1041         isymij  = isymkl
1042         do isyml =1, nsym
1043            isymk = muld2h(isymkl,isyml)
1044           do k = 1, nrhf(isymk)
1045             do l = 1, nrhf(isyml)
1046                idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
1047                 do isyma =1, nsym
1048                    isymn  = muld2h(isyman,isyma)
1049                    isymi  = muld2h(isymij,isymn)
1050                    nvir0(isyma)  = norb1(isyma) - nrhfs(isyma)
1051                    ntota = max(1,nvir0(isyma))
1052                    ntotn = max(1,nrhf(isymn))
1053                    ntoti = max(1,nrhf(isymi))
1054                    koffd = kdotil+ igamsq(isymij,isymkl)
1055     &              +  nmatij(isymij)*(idxkl-1)+imatij(isymi,isymn)
1056                    klijanz = klijan + imatan(isyma,isymi)
1057                    call dgemm('N','T',nvir0(isyma),nrhf(isymi),
1058     &                         nrhf(isymn),one,work(kxajkl),
1059     &                         ntota,work(koffd),ntoti,one,
1060     &                         work(klijanz),ntota)
1061                    kxajkl = kxajkl + nvir0(isyma)*nrhf(isymn)
1062                    blajklf = blajklf + nvir0(isyma)*nrhf(isymn)
1063
1064                 enddo
1065              enddo
1066           enddo
1067         enddo
1068      enddo
1069
1070
1071
1072      call r12gradaxd(ipdd,work(kccr12),work(kctil),work(klij),
1073     &                work(kend2),
1074     &                v2am,res,lwrk2)
1075      call daxpy(ncvai,one,work(klijan),1,res(1),1)
1076
1077
1078      call qexit('r12gradap')
1079      end
1080c--------------------------------------------------------------------------
1081C=====================================================================*
1082C  /* Deck r12gradaxd */
1083      subroutine r12gradaxd(ipdd,kccr12,kctil,dklij,work,v2am,res,
1084     &                      lwork)
1085      implicit none
1086#include "ccsdsym.h"
1087#include "ccorb.h"
1088#include "dummy.h"
1089#include "r12int.h"
1090#include "priunit.h"
1091#include "ccfro.h"
1092      character*8 CCR12YIJ
1093      integer krsing,krtrip,kend1,lwrk1,nrhftria,n2
1094      integer lwork,lu43,idxkl,isymv,isymi,isymj,isymk,idxij
1095      integer isyml,isymkl,isymij,idxji,idxklij,idxijkl,idxklji
1096      integer isymia,luvajkl,ntklaj,isymaj,kend2
1097      integer nmatan(8),isyman,isyma,isymn,ntotij,ntotan,nvir0(8)
1098      integer lctil,klijan,ntota,ntoti,ntotn,lcvai,isymai
1099      integer kijan,ntotk,koffc,lccr12,koffd
1100      integer idxijlk,idxlk,isymin,kbajkl,ncvai,imatan(8,8),icou3
1101      integer lwrk2,lctil1,kcotil,kdotil,koccr,koctil
1102
1103      integer kxijkl,kaklij,nak,isymik,nai,naiki,nkj,naikj,isymjk
1104      integer ntotmu,dklijz,eklij,eklijz,lxxr12,kxxr12,kxsing,kxtrip
1105      integer kiakj,nia,njk,niajk,idx,ntotj,kgiakj,kgaijk,kaxdia
1106      integer fklij,fklijz,kxcijkl,kxcijklz,ntklij,ntaj
1107      integer naj,nki,najki,nji,naijk,nakji,kxcijklt,nka
1108      integer rklij,koffh,dijkl,hklij,kaxcia,kaxdai,kiajk,kgaikj
1109      integer dlkmnz,dlkmn,kctil1,kcxijklt,kcxijkl,krcc12,lunit
1110      integer kcxklmnz,kcxmnklz,kcxmnkl,kcxklmn,kxxr12n,idxjikl,idxjilk
1111      integer aijkl,hijkl,hklijz,lxxr12n,dijklz,eijkl,eijklz
1112      integer njaik,nik,nja,najik,nikja,aklij,kctil2,aijklz,lrcc12,bijkl
1113      integer isymja,isymka,ntoto,isymoj,isymo,idxlkij,kres,nikaj
1114      integer kiakjz,aj,ajkl,kl,isymkj,kiajkz,rklijz,index
1115      integer ipdd,resb,kxtil,kresz,icoun4,nmataj(8),ncvij,nijkl(8)
1116#if defined (SYS_CRAY)
1117      real work(*),dklij(*),v2am(*),kctil(*),kccr12(*),one,two,zero
1118      real res(*)
1119#else
1120      double precision work(*),dklij(*),v2am(*),kctil(*),one,zero,two
1121      double precision kccr12(*),res(*)
1122#endif
1123      parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0)
1124
1125      index(i,j) = max(i,j)*(max(i,j)-3)/2 +i+j
1126
1127      call qenter('r12gradaxd')
1128
1129C     Calculate some offsets
1130
1131      isymv = 1
1132
1133      nrhftria = nrhft*(nrhft+1)/2
1134      n2 = nrhftria * nrhftria
1135
1136
1137
1138      do isymaj = 1,nsym
1139         isymij = isymaj
1140         ncvij = 0
1141         do isyma = 1,nsym
1142            isymj = muld2h(isymaj,isyma)
1143            isymi = muld2h(isymij,isymj)
1144            ncvij = ncvij +  nrhf(isyma)*nrhf(isymi)
1145         enddo
1146      enddo
1147
1148
1149
1150
1151      do isymaj = 1,nsym
1152         isymij = isymaj
1153        ntaj = 0
1154         do isyma = 1,nsym
1155            isymj = muld2h(isymaj,isyma)
1156            isymi = muld2h(isymij,isymj)
1157            ntaj = ntaj + (norb1(isyma) - nrhfs(isyma))*nrhf(isymi)
1158         enddo
1159      enddo
1160
1161
1162      ntklaj = 0
1163      do isymkl = 1,nsym
1164         isymaj = isymkl
1165         do isymj = 1,nsym
1166            isyma = muld2h(isymaj,isymj)
1167            ntklaj = ntklaj + nmatij(isymkl)*nrhf(isymj)
1168     &               *(norb1(isyma) - nrhfs(isyma))
1169         enddo
1170      enddo
1171
1172
1173
1174      kxxr12  = 1
1175      kxsing  = kxxr12  + ngamsq(1)
1176      kxtrip  = kxsing  + n2
1177      eklij   = kxtrip  + n2
1178      kgiakj  = eklij   + ncvij
1179      kgaijk  = kgiakj  + ntklaj
1180      kxcijkl = kgaijk  + ntklaj
1181      fklij   = kxcijkl + ngamsq(1)
1182      kiajk   = fklij   + ncvij
1183      kgaikj  = kiajk   + ntklaj
1184      rklij   = kgaikj  + ntklaj
1185      kcxijkl = rklij   + ncvij
1186      eijkl   = kcxijkl + ngamsq(1)
1187      kxtil   = eijkl   + ncvij
1188      resb    = kxtil   + ngamsq(1)
1189      kend2   = resb    + ntaj
1190      lwrk2   = lwork   - kend2
1191
1192C     Read matrix X
1193
1194      lu43 = -43
1195      call gpopen(lu43,'CCR12_XP','unknown',' ','formatted',
1196     &                   idummy,.false.)
1197      read(lu43,*) i
1198      read(lu43,'(4e30.20)') (work(kxsing+i), i = 0, n2-1 )
1199      read(lu43,'(4e30.20)') (work(kxtrip+i), i = 0, n2-1 )
1200      if (ipdd.eq.3) then
1201          call gpclose(lu43,'KEEP')
1202      else
1203          call gpclose(lu43,'DELETE')
1204      endif
1205      call dscal(n2,2d0,work(kxsing),1)
1206      call dscal(n2,2d0/3d0,work(kxtrip),1)
1207
1208      call cc_r12vpack(work(kxxr12),isymv,work(kxsing),work(kxtrip),
1209     &     nrhf,nrhfa,nijkl)
1210
1211
1212C     Calculate XC_{ij}^{kl} = X_{ij}^{mn} * Ctilde_{mn}^{kl}
1213
1214      call dzero(work(kxcijkl),ngamsq(1))
1215      do isymkl = 1,nsym
1216         isymij  = isymkl
1217         lxxr12   = kxxr12 + igamsq(isymij,isymkl)
1218         lctil    = 1 + igamsq(isymij,isymkl)
1219         kxcijklz = kxcijkl + igamsq(isymij,isymkl)
1220         ntoti = max(1,nmatij(isymij))
1221         ntotk = max(1,nmatij(isymkl))
1222         call dgemm('N','T',nmatij(isymij),nmatij(isymkl),
1223     &              nmatij(isymkl),one,work(lxxr12),
1224     &              ntotk,kctil(lctil),ntoti,zero,
1225     &              work(kxcijklz),ntotk)
1226      end do
1227c     ---------------------------------
1228C     Transpone CX_{ij}^{kl}
1229c     ---------------------------------
1230
1231      call dzero(work(kcxijkl),ngamsq(1))
1232      do isymkl = 1, nsym
1233         isymij = muld2h(isymkl,isymv)
1234         do isyml =1, nsym
1235            isymk = muld2h(isymkl,isyml)
1236            do k = 1, nrhf(isymk)
1237               do l = 1, nrhf(isyml)
1238                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
1239                  do isymi = 1, nsym
1240                     isymj = muld2h(isymij,isymi)
1241                     do j = 1, nrhf(isymj)
1242                        do i = 1, nrhf(isymi)
1243                           idxij = imatij(isymi,isymj)+
1244     &                          nrhf(isymi)*(j-1) + i
1245                           idxklij =  igamsq(isymkl,isymij)+
1246     &                          nmatij(isymkl)*(idxij-1)+idxkl
1247                           idxijkl = igamsq(isymij,isymkl)+
1248     &                          nmatij(isymij)*(idxkl-1) +idxij
1249                           work(kcxijkl-1+idxijkl) =
1250     &                         work(kxcijkl-1+idxklij)
1251                        end do
1252                     end do
1253
1254                  end do
1255               end do
1256            end do
1257         end do
1258      end do
1259
1260
1261
1262      call dzero(work(eklij),ncvij)
1263      call dzero(work(eijkl),ncvij)
1264      call dzero(work(fklij),ncvij)
1265      do isymkl = 1, nsym
1266         isymoj  = isymkl
1267         isymij = isymkl
1268         do isyml =1, nsym
1269            isymk = muld2h(isymkl,isyml)
1270           do k = 1, nrhf(isymk)
1271             do l = 1, nrhf(isyml)
1272                idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
1273                 do isymo =1, nsym
1274                    isymj  = muld2h(isymoj,isymo)
1275                    isymi  = muld2h(1,isymo)
1276
1277C     Calculate X_{ol}^{ij} * D_{ij}^{kl}
1278
1279                    koffd  = 1 + igamsq(isymij,isymkl)
1280     &                 + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
1281                    lxxr12 = kxxr12 + igamsq(isymij,isymkl)
1282     &                 + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
1283                    ntotj = max(1,nrhf(isymj))
1284                    ntoti = max(1,nrhf(isymi))
1285                    ntoto = max(1,nrhf(isymo))
1286                    eklijz = eklij + imatij(isymo,isymi)
1287                    call dgemm('N','T',nrhf(isymo),nrhf(isymi),
1288     &                    nrhf(isymj),one,work(lxxr12),
1289     &                    ntoto,dklij(koffd),ntoti,one,
1290     &                    work(eklijz),ntoto)
1291                    eijklz = eijkl + imatij(isymi,isymo)
1292                    call dgemm('N','T',nrhf(isymi),nrhf(isymo),
1293     &                    nrhf(isymj),one,dklij(koffd),
1294     &                    ntoti,work(lxxr12),ntoto,one,
1295     &                    work(eijklz),ntoti)
1296
1297C     Calculate   F_{}^{} = C_{ij}^{kl} * XC_{ij}^{ol}
1298
1299                    lccr12 = 1+ igamsq(isymij,isymkl)
1300     &                 + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
1301                    kcxijklt = kcxijkl + igamsq(isymij,isymkl)
1302     &                 + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
1303                    fklijz = fklij + imatij(isymo,isymi)
1304                    call dgemm('N','T',nrhf(isymi),nrhf(isymo),
1305     &                    nrhf(isymj),one,kccr12(lccr12),
1306     &                    ntoti,work(kcxijklt),ntoti,one,
1307     &                    work(fklijz),ntoti)
1308
1309                 enddo
1310              enddo
1311            enddo
1312         enddo
1313      end do
1314
1315      call dzero(work(rklij),ncvij)
1316      call daxpy(ncvij,one,work(eklij),1,work(rklij),1)
1317      call daxpy(ncvij,one,work(eijkl),1,work(rklij),1)
1318      call daxpy(ncvij,-two,work(fklij),1,work(rklij),1)
1319      call dscal(ncvij,two,work(rklij),1)
1320      lunit = -1
1321      if (ipdd .eq. 3) then
1322         call gpopen(lunit,'CCR12YIJ','unknown',' ','unformatted',
1323     &               idummy,.false.)
1324      else if (ipdd .eq. 5) then
1325         call gpopen(lunit,'CCR12YIJB','unknown',' ','unformatted',
1326     &               idummy,.false.)
1327      endif
1328      write(lunit)(work(rklij+i-1),i=1,ncvij)
1329      call gpclose(lunit,'KEEP')
1330
1331
1332      call dzero(work(kgiakj),ntklaj)
1333      call dzero(work(kgaijk),ntklaj)
1334      call dzero(work(kgaikj),ntklaj)
1335      call dzero(work(kiajk),ntklaj)
1336
1337      idx = 0
1338      do isymka = 1,nsym
1339         isymij = isymka
1340         do isyma = 1, nsym
1341            nvir0(isyma) = norb1(isyma) - nrhfs(isyma)
1342            isymk = muld2h(isymka,isyma)
1343            do isymj = 1,nsym
1344               isymi = muld2h(isymij,isymj)
1345               do i = 1,nrhf(isymi)
1346                  do j = 1,nrhf(isymj)
1347                     isymia = muld2h(isymi,isyma)
1348                     isymkj = muld2h(isymk,isymj)
1349                     isymja = muld2h(isymj,isyma)
1350                     isymik = muld2h(isymi,isymk)
1351                     do k = 1, nrhf(isymk)
1352                        do a = 1,nvir0(isyma)
1353                          idx = idx + 1
1354                          nai = it1am(isyma,isymi)
1355     *                        + nvir(isyma)*(i-1)
1356     &                         + a + nrhfs(isyma)
1357                          njk = it1am(isymj,isymk)
1358     *                         + nvir(isymj)*(k-1) + j+nrhffr(isymj)
1359                          naijk = it2am(isymia,isymkj)
1360     *                           + index(nai,njk)
1361                          nak = it1am(isyma,isymk)
1362     *                         + nvir(isyma)*(k-1)
1363     &                         + a + nrhfs(isyma)
1364                          nji = it1am(isymj,isymi)
1365     *                          + nvir(isymj)*(i-1) + j+nrhffr(isymj)
1366                          nakji = it2am(isymka,isymij)
1367     *                            + index(nak,nji)
1368                          naj = it1am(isyma,isymj)
1369     *                          + nvir(isyma)*(j-1)
1370     &                          + a + nrhfs(isyma)
1371                          nki = it1am(isymk,isymi)
1372     *                          + nvir(isymk)*(i-1) + k +nrhffr(isymk)
1373                          nik = it1am(isymi,isymk)
1374     *                          + nvir(isymi)*(k-1) + i
1375                          najki = it2am(isymja,isymik)
1376     *                            + index(naj,nki)
1377                          work(kgiakj+idx-1) = v2am(naijk)
1378                          work(kgaijk+idx-1) = v2am(nakji)
1379                          work(kgaikj+idx-1) = v2am(najki)
1380                          work(kiajk+idx-1)  = 4.0d0*work(kgaijk+idx-1)
1381     &                        - work(kgaikj+idx-1) - work(kgiakj+idx-1)
1382                      enddo
1383                    enddo
1384                  enddo
1385               enddo
1386            enddo
1387         enddo
1388      enddo
1389
1390
1391
1392      do isymi = 1,nsym
1393          nvir0(isymi) = norb1(isymi)-nrhfs(isymi)
1394          nmataj(isymi) = nvir0(isymi)*nrhf(isymi)
1395      enddo
1396
1397
1398
1399      call dzero(res(1),ntaj)
1400      kres = 1
1401      do isymij = 1,nsym
1402         isymkl = isymij
1403         isymaj = isymij
1404         do isyma =1, nsym
1405            isymj  = muld2h(isymaj,isyma)
1406            isymi  = muld2h(1,isyma)
1407            nvir0(isyma) = norb1(isyma)-nrhfs(isyma)
1408            ntotk = max(1,nvir0(isyma)*nrhf(isymj))
1409            ntoti = max(1,nmatij(isymkl))
1410            call dgemm('N','N',nrhf(isymj)*nvir0(isyma),1,
1411     &                 nmatij(isymkl),one,work(kiajk),ntotk,
1412     &                 work(rklij),ntoti,zero,res(kres),ntotk)
1413            kres = kres + nvir0(isyma)*nrhf(isymj)
1414            kiajk = kiajk + nrhf(isymj)*nvir0(isyma)
1415     &                      * nmatij(isymkl)
1416         enddo
1417      enddo
1418
1419      call dscal(ntaj,0.5d0,res(1),1)
1420
1421
1422      if (ipdd .eq. 5) then
1423         call dscal(ntaj,1.0d0,res(1),1)
1424         call dzero(work(resb),ntaj)
1425         call r12gradb(work(resb),dklij,work(kend2),ntklaj,ntaj,lwrk2)
1426
1427         call daxpy(ntaj,one,work(resb),1,res(1),1)
1428      endif
1429      call qexit('r12gradaxd')
1430      end
1431c--------------------------------------------------------------------------
1432C=====================================================================*
1433C  /* Deck r12gradb */
1434      subroutine r12gradb(res,dklij,work,ntklaj,ncvai,lwork)
1435      implicit none
1436#include "ccsdsym.h"
1437#include "ccorb.h"
1438#include "dummy.h"
1439#include "r12int.h"
1440#include "priunit.h"
1441#include "ccfro.h"
1442      logical locdbg
1443      parameter (locdbg = .false.)
1444      integer kqajkl,kqijal,kuijal,krajkl,kend1
1445      integer isym,luvajkl,nvir0t,lwork,lwrk1,ntklaj
1446      integer icou3,imataj(8,8),koffd,ksajkl,nvir0(8),ncvai
1447      integer kcqia,kcqai,lcqia,lcqai,isymai,isyma,isymi
1448      integer isymaj,isymkl,isymij,isyml,isymk,isymj,idxkl,ntotj
1449      integer isymv,idxlk,idxklij,idxlkij,idxij,idxji
1450      integer idxjikl,idxijkl,idxijlk,koff,ntota,ntoti
1451      integer kuajkl,ktajkl,ktijal,kcuia,kcuai,lcuai,lcuia
1452      integer kssijal,ksrajkl,kstijal,idxklji,koff1,kdjikl
1453      integer kuijals,kqijals,kssajkl
1454
1455      integer kccr12,krsing,krtrip,kctil,lccr12,lctil
1456      integer nrhftria,n2,lu43,ntotk
1457#if defined (SYS_CRAY)
1458      real dklij(*),res(*),work(*),one,two,zero
1459#else
1460      double precision dklij(*),res(*),work(*),one,zero,two
1461#endif
1462      parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0)
1463
1464      call qenter('r12gradb')
1465
1466      isymv = 1
1467
1468      nrhftria = nrhft*(nrhft+1)/2
1469      n2 = nrhftria * nrhftria
1470
1471
1472      nvir0t = 0
1473      do isym = 1, nsym
1474         nvir0(isym) = norb1(isym) - nrhfs(isym)
1475         nvir0t      = nvir0t + nvir0(isym)
1476      end do
1477
1478      icou3 = 0
1479      do isymai = 1,nsym
1480         icou3 = 0
1481         do isyma = 1,nsym
1482            isymi = muld2h(isymai,isyma)
1483            imataj(isyma,isymi)= icou3
1484            icou3 = icou3+nrhf(isymi)*(norb1(isyma) - nrhfs(isyma))
1485         enddo
1486      enddo
1487
1488
1489      kqajkl = 1
1490      kqijal = kqajkl + ntklaj
1491      kuijal = kqijal + ntklaj
1492      krajkl = kuijal + ntklaj
1493      ksajkl = krajkl + nvir0t*nrhft**3
1494      kcqai  = ksajkl + nvir0t*nrhft**3
1495      kcqia  = kcqai  + ncvai
1496      ktijal = kcqia  + ncvai
1497      kcuai  = ktijal + nvir0t*nrhft**3
1498      kuajkl = kcuai  + ncvai
1499      kcuia  = kuajkl + ntklaj
1500      kstijal= kcuia  + ncvai
1501      ksrajkl= kstijal+ ntklaj
1502      kssijal= ksrajkl+ ntklaj
1503      kssajkl= kssijal+ ntklaj
1504      kuijals= kssajkl+ ntklaj
1505      kqijals= kuijals+ ntklaj
1506      kend1  = kqijals+ ntklaj
1507      lwrk1  = lwork  - kend1
1508
1509
1510      luvajkl = -1
1511      call dzero(work(kqajkl),ntklaj)
1512      call gpopen(luvajkl,'DDR12QAJKL','unknown',' ','unformatted',
1513     &                   idummy,.false.)
1514      rewind(luvajkl)
1515      read(luvajkl) (work(kqajkl+i-1), i = 1,ntklaj)
1516      call gpclose(luvajkl,'DELETE')
1517
1518
1519      luvajkl = -1
1520      call dzero(work(kqijal),ntklaj)
1521      call gpopen(luvajkl,'DDR12QIJAL','unknown',' ','unformatted',
1522     &                   idummy,.false.)
1523      rewind(luvajkl)
1524      read(luvajkl) (work(kqijal+i-1), i = 1,ntklaj)
1525      call gpclose(luvajkl,'DELETE')
1526
1527      call dzero(work(kqijals),ntklaj)
1528
1529      luvajkl = -1
1530      call dzero(work(kuijal),ntklaj)
1531      call gpopen(luvajkl,'DDR12UIJAL','unknown',' ','unformatted',
1532     &                   idummy,.false.)
1533      rewind(luvajkl)
1534      read(luvajkl) (work(kuijal+i-1), i = 1,ntklaj)
1535          call gpclose(luvajkl,'DELETE')
1536
1537      call dzero(work(kuijals),ntklaj)
1538
1539
1540      luvajkl = -1
1541      call dzero(work(kuajkl),ntklaj)
1542      call gpopen(luvajkl,'DDR12UAJKL','unknown',' ','unformatted',
1543     &                   idummy,.false.)
1544      rewind(luvajkl)
1545      read(luvajkl) (work(kuajkl+i-1), i = 1,ntklaj)
1546      call gpclose(luvajkl,'DELETE')
1547
1548      call dzero(work(krajkl),nvir0t*nrhft**3)
1549      luvajkl = -1
1550      call gpopen(luvajkl,'AUXQSA12','unknown',' ','formatted',
1551     &                   idummy,.false.)
1552      rewind(luvajkl)
1553      read(luvajkl,'(4E30.20)') (work(krajkl+i-1), i=1,nvir0t*nrhft**3)
1554      call gpclose(luvajkl,'DELETE')
1555
1556      call dzero(work(ksrajkl),ntklaj)
1557      call cc_r12xsort2(work(krajkl),work(ksrajkl),nvir0)
1558
1559
1560      call dzero(work(ksajkl),nvir0t*nrhft**3)
1561      luvajkl = -1
1562      call gpopen(luvajkl,'AUXQSAJ12','unknown',' ','formatted',
1563     &                   idummy,.false.)
1564      rewind(luvajkl)
1565      read(luvajkl,'(4E30.20)') (work(ksajkl+i-1), i=1,nvir0t*nrhft**3)
1566      call gpclose(luvajkl,'DELETE')
1567
1568      call dzero(work(kssijal),ntklaj)
1569      call cc_r12xsort2(work(ksajkl),work(kssijal),nvir0)
1570
1571      call dzero(work(kssajkl),ntklaj)
1572      call cc_r12xsort3(work(ksajkl),work(kssajkl),nvir0)
1573
1574
1575
1576      call dzero(work(ktijal),nvir0t*nrhft**3)
1577      luvajkl = -1
1578      call gpopen(luvajkl,'AUXQSAI12','unknown',' ','formatted',
1579     &                   idummy,.false.)
1580      rewind(luvajkl)
1581      read(luvajkl,'(4E30.20)') (work(ktijal+i-1), i=1,nvir0t*nrhft**3)
1582      call gpclose(luvajkl,'DELETE')
1583
1584      call dzero(work(kstijal),ntklaj)
1585      call cc_r12xsort3(work(ktijal),work(kstijal),nvir0)
1586
1587
1588      call daxpy(ntklaj,one,work(kssijal),1,work(kqijal),1)
1589      call daxpy(ntklaj,one,work(kssajkl),1,work(kuijal),1)
1590      call daxpy(ntklaj,one,work(kstijal),1,work(kuajkl),1)
1591      call daxpy(ntklaj,one,work(ksrajkl),1,work(kqajkl),1)
1592
1593
1594
1595      call dzero(work(kcqai),ncvai)
1596      lcqai = kcqai
1597      call dzero(work(kcqia),ncvai)
1598      lcqia = kcqia
1599      call dzero(work(kcuai),ncvai)
1600      lcuai = kcuai
1601      call dzero(work(kcuia),ncvai)
1602      lcuia = kcuia
1603
1604c     Calculate Q_{aj}^{kl}*D_{kl}^{ij}
1605
1606
1607      do isymkl = 1, nsym
1608         isymaj  = isymkl
1609         isymij  = isymkl
1610         do isyml =1, nsym
1611            isymk = muld2h(isymkl,isyml)
1612           do k = 1, nrhf(isymk)
1613             do l = 1, nrhf(isyml)
1614                idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
1615                 do isyma =1, nsym
1616                    isymj  = muld2h(isymaj,isyma)
1617                    isymi  = muld2h(isymij,isymj)
1618                    nvir0(isyma)  = norb1(isyma) - nrhfs(isyma)
1619                    koffd  = 1 + igamsq(isymij,isymkl)
1620     &                 + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
1621                    ntota = max(1,nvir0(isyma))
1622                    ntotj = max(1,nrhf(isymj))
1623                    ntoti = max(1,nrhf(isymi))
1624                    lcqai  = kcqai + imataj(isyma,isymi)
1625                    call dgemm('N','T',nvir0(isyma),nrhf(isymi),
1626     &                    nrhf(isymj),one,work(kqajkl),
1627     &                    ntota,dklij(koffd),ntoti,one,
1628     &                    work(lcqai),ntota)
1629
1630                    kqajkl = kqajkl + nvir0(isyma)*nrhf(isymj)
1631                    ksrajkl = ksrajkl + nvir0(isyma)*nrhf(isymj)
1632
1633                    ntota = max(1,nvir0(isyma))
1634                    ntotj = max(1,nrhf(isymj))
1635                    ntoti = max(1,nrhf(isymi))
1636                    lcqia  = kcqia + imataj(isyma,isymi)
1637                    call dgemm('N','T',nvir0(isyma),nrhf(isymi),
1638     &                    nrhf(isymj),one,work(kqijal),
1639     &                    ntota,dklij(koffd),ntoti,one,
1640     &                    work(lcqia),ntota)
1641
1642                    kqijal = kqijal + nvir0(isyma)*nrhf(isymj)
1643
1644                    lcuai  = kcuai + imataj(isyma,isymi)
1645                    call dgemm('N','T',nvir0(isyma),nrhf(isymi),
1646     &                    nrhf(isymj),one,work(kuijal),
1647     &                    ntota,dklij(koffd),ntoti,one,
1648     &                    work(lcuai),ntota)
1649
1650                    kuijal = kuijal + nvir0(isyma)*nrhf(isymj)
1651                    kssijal = kssijal + nvir0(isyma)*nrhf(isymj)
1652                    kssajkl = kssajkl + nvir0(isyma)*nrhf(isymj)
1653
1654                    lcuia  = kcuia + imataj(isyma,isymi)
1655                    call dgemm('N','T',nvir0(isyma),nrhf(isymi),
1656     &                    nrhf(isymj),one,work(kuajkl),
1657     &                    ntota,dklij(koffd),ntoti,one,
1658     &                    work(lcuia),ntota)
1659
1660                    kuajkl  = kuajkl + nvir0(isyma)*nrhf(isymj)
1661                    kstijal = kstijal + nvir0(isyma)*nrhf(isymj)
1662                    ktijal  = ktijal + nvir0(isyma)*nrhf(isymj)
1663
1664                 enddo
1665              enddo
1666           enddo
1667         enddo
1668      enddo
1669
1670
1671      call r12gradbd(res,dklij,work(kend1),lwrk1)
1672
1673      call daxpy(ncvai,one,work(kcqai),1,res,1)
1674      call daxpy(ncvai,one,work(kcuai),1,res,1)
1675      call daxpy(ncvai,one,work(kcqia),1,res,1)
1676      call daxpy(ncvai,one,work(kcuia),1,res,1)
1677
1678      if (locdbg) then
1679         write(lupri,*) 'Ergebnisse r12gradb'
1680         do isymaj = 1,nsym
1681            do isyma = 1,nsym
1682               isymi =isyma
1683               write(lupri,*) 'symmetry block',isyma
1684               call output(res(1+imataj(isyma,isymi)),
1685     &                 1,nvir0(isyma),1,nrhf(isymi),
1686     &                 nvir0(isyma),nrhf(isymi),1,lupri)
1687            enddo
1688         enddo
1689       endif
1690
1691
1692      call qexit('r12gradb')
1693      end
1694c------------------------------------------------------------------------
1695C=====================================================================*
1696C  /* Deck r12gradbd */
1697      subroutine r12gradbd(res,dklij,work,lwork)
1698      implicit none
1699#include "ccsdsym.h"
1700#include "ccorb.h"
1701#include "dummy.h"
1702#include "r12int.h"
1703#include "priunit.h"
1704#include "ccfro.h"
1705#include "ccsdinp.h"
1706      integer lxdmmu,nr1bast,isymak,isymk,isyma,nr1bas(8)
1707      integer lwork,kxmujkl,kxdmmu,kend1,nmum,luvajkl
1708      integer isymkl,isymij,isymi,isymj,idxkl,koffd,isymaj
1709      integer isymai,icou3,imatmuj(8,8),ntota,ntoti,ntotj,isyml
1710      integer smujkl,koffc,lusifc,lwrk1,kcmo,ntota1,kxdnumu
1711      integer kxpjkl,ntklpj,isymmua,iglmrv(8,8),isymmu
1712      integer kxajkl,ntotmu,isymmuj,kxpjklz,icount2,icou2,icount7
1713      integer nr1xbas(8),koff,idxji,idxij,kzmunu,kxdpmu,ir1bas(8,8)
1714      integer nrhftria,kxsing,kxtrip,n2,isymv,lu43,nxajkl(8),kcmoa
1715      integer isym,noffset,kspjkl,nvir0t,nvir0(8),ntklaj,krajkl
1716      integer imbas1(8),kdens,lxdnumu,icou4,imatmunu(8,8),lxdpmu
1717      integer nrhfst,ksqijklz,ksqjkl,ksajkl,sajkl,icoun1,norb1t
1718      integer imataj(8,8),idx,kl,aj,ajkl,kspjklf,ksajklf
1719      integer imatmn(8,8),icou6,nmatkl1(8),imatkl1(8,8),icou5
1720      integer nmatan(8),icou7
1721      integer imatoi(8,8),icou8,isyman,isymn
1722      integer iglmrv1(8,8),noffset1,icou9,nrhffrt
1723      integer ksijklf,ij,ijkl,ksijkl,sijklf
1724#if defined (SYS_CRAY)
1725      real res(*),work(*),dklij(*),zero,one,two
1726#else
1727      double precision res(*),work(*),dklij(*),zero,one,two
1728#endif
1729      parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0)
1730
1731      call qenter('r12gradbd')
1732
1733      isymv = 1
1734
1735      nrhftria = nrhft*(nrhft+1)/2
1736      n2 = nrhftria * nrhftria
1737
1738
1739      nvir0t = 0
1740      norb1t = 0
1741      nrhfst = 0
1742      nrhffrt = 0
1743      do isymi = 1, nsym
1744         nvir0(isymi) = norb1(isymi) - nrhfs(isymi)
1745         nvir0t      = nvir0t + nvir0(isymi)
1746         norb1t      = norb1t + norb1(isymi)
1747         nrhfst      = nrhfst + nrhfs(isymi)
1748         nrhffrt     = nrhffrt+ nrhffr(isymi)
1749      end do
1750
1751
1752      do isymmua = 1,nsym
1753         icou2 = 0
1754         icou3 = 0
1755         do isyma = 1,nsym
1756            isymmu = muld2h(isymmua,isyma)
1757            iglmrv(isymmu,isyma) = icou2
1758            iglmrv1(isymmu,isyma) = icou3
1759            icou2 = icou2 + nbas(isymmu)*(nvir(isyma))
1760            icou3 = icou3 + nbas(isymmu)*(nvir(isyma)-nrhffr(isyma))
1761         enddo
1762      enddo
1763
1764      do isyman = 1,nsym
1765         icou3 = 0
1766         icou4 = 0
1767         icou7 = 0
1768         nmatkl1(isyman) = 0
1769         do isyma = 1,nsym
1770            isymn = muld2h(isyman,isyma)
1771            nvir0(isyma)  = norb1(isyma) - nrhfs(isyma)
1772            icou3 = icou3 + nvir0(isyma)*nrhfs(isymn)
1773            icou4 = icou4 + nrhfs(isyma)*nrhfs(isymn)
1774          enddo
1775          nmatan(isyman) = icou3
1776          nmatkl1(isyman) = icou4
1777      enddo
1778
1779      do isyman = 1,nsym
1780         icou4 = 0
1781         icou5 = 0
1782         icou6 = 0
1783         icou8 = 0
1784         do isyma = 1,nsym
1785            isymn = muld2h(isyman,isyma)
1786            imataj(isyma,isymn)= icou4
1787            imatkl1(isyma,isymn)= icou5
1788            imatmn(isyma,isymn)= icou6
1789            imatoi(isyma,isymn)= icou8
1790            icou4 = icou4+nrhfs(isymn)*(norb1(isyma) - nrhfs(isyma))
1791            icou5 = icou5+nmatan(isymn)* nmatkl1(isyma)
1792            icou6 = icou6+nrhfs(isyma)*nrhfs(isymn)
1793            icou8 = icou8+nmatkl1(isyma)*nmatkl1(isymn)
1794         enddo
1795      enddo
1796
1797
1798
1799      do isymak = 1, nsym
1800        nr1bas(isymak) = 0
1801        nr1xbas(isymak) = 0
1802        do isymk = 1, nsym
1803           isyma = muld2h(isymak,isymk)
1804           nr1bas(isymak)  = nr1bas(isymak) +mbas1(isyma)*nrhf(isymk)
1805           nr1xbas(isymak) = nr1xbas(isymak)+mbas2(isyma)*nrhf(isymk)
1806        end do
1807      end do
1808
1809      nr1bast = 0
1810      do isyma = 1,nsym
1811         nr1bast = nr1bast + nr1bas(isyma)
1812      enddo
1813
1814      ntklpj = 0
1815      do isymkl = 1,nsym
1816         isymaj = isymkl
1817         do isymj = 1,nsym
1818            isyma = muld2h(isymaj,isymj)
1819            ntklpj = ntklpj + nmatij(isymkl)*nrhf(isymj)*norb1(isyma)
1820         enddo
1821      enddo
1822
1823
1824
1825      do isymak = 1, nsym
1826         nvajkl(isymak) = 0
1827         icount7        = 0
1828         icount2 = 0
1829         do isymk = 1, nsym
1830            isyma = muld2h(isymak,isymk)
1831            ivajkl(isyma,isymk) = icount7
1832            nvajkl(isymak) = nvajkl(isymak)+ nr1bas(isyma)*nmatij(isymk)
1833            icount7 = icount7 + nr1bas(isyma)*nmatij(isymk)
1834            ir1bas(isyma,isymk)  = icount2
1835            icount2 = icount2 + nrhf(isymk)*mbas1(isyma)
1836         end do
1837      end do
1838
1839      icou3 = 0
1840      do isymai = 1,nsym
1841         icou3 = 0
1842         icou4 = 0
1843         do isyma = 1,nsym
1844            isymi = muld2h(isyma,isymai)
1845            imatmuj(isyma,isymi)= icou3
1846            icou3 = icou3+nrhf(isymi)*norb1(isyma)
1847            imatmunu(isymi,isyma) = icou4
1848            icou4 = icou4+mbas1(isymi)*mbas1(isyma)
1849         enddo
1850      enddo
1851
1852      ntklaj = 0
1853      do isymkl = 1, nsym
1854         isymaj  = isymkl
1855         do isyml =1, nsym
1856            isymk = muld2h(isymkl,isyml)
1857            do k = 1, nrhf(isymk)
1858               do l = 1, nrhf(isyml)
1859                  do isyma =1, nsym
1860                     isymj  = muld2h(isymaj,isyma)
1861                     isymi  = isyma
1862                     ntklaj = ntklaj + nrhf(isymj) *
1863     &                        (norb1(isyma) - nrhfs(isyma))
1864                  enddo
1865               enddo
1866            enddo
1867         enddo
1868      enddo
1869
1870
1871
1872      kxdmmu  = 1
1873      kxdnumu = kxdmmu  + nr1bast
1874      kcmo    = kxdnumu + 10*mbas1t*mbas1t
1875      smujkl  = kcmo    + nlamds
1876      kxpjkl  = smujkl  + nrhfst**4
1877      kxdpmu  = kxpjkl  + ntklpj
1878      kzmunu  = kxdpmu  + mbas1t*nrhft
1879      kcmoa   = kzmunu  + n2basx
1880      kxmujkl = kcmoa   + nlamds
1881      kspjkl  = kxmujkl + nvajkl(1)
1882      kdens   = kspjkl  + ntklpj
1883      sajkl   = kdens   + n2basx
1884      ksajkl  = sajkl   + nvir0t*nrhfst**3
1885      ksqjkl  = ksajkl  + nvir0t*nrhfst**3
1886      ksajklf = ksqjkl  + ntklpj
1887      kspjklf = ksajklf + nvir0t*nrhfst**3
1888      ksijklf = kspjklf + nrhfst**4
1889      ksijkl  = ksijklf + nrhfst*nrhfst**3
1890      kend1   = ksijkl  + nrhfst*nrhfst**3
1891      lwrk1   = lwork   - kend1
1892
1893
1894      call dzero(work(kxmujkl),nvajkl(1))
1895      luvajkl = -1
1896      call gpopen(luvajkl,'CCR12QAJKL','unknown',' ','unformatted',
1897     &               idummy,.false.)
1898      rewind(luvajkl)
1899      read(luvajkl) (work(kxmujkl+i-1), i = 1,nvajkl(1))
1900      call gpclose(luvajkl,'KEEP')
1901c
1902      call dzero(work(smujkl),nrhfst**4)
1903      luvajkl = -1
1904
1905      call gpopen(luvajkl,'AUXQ12','unknown',' ','formatted',
1906     &               idummy,.false.)
1907      rewind(luvajkl)
1908      read(luvajkl,'(4E30.20)') (work(smujkl+i-1), i = 1,
1909     &                        nrhfst**4)
1910      call gpclose(luvajkl,'KEEP')
1911c
1912      luvajkl = -1
1913      call gpopen(luvajkl,'AUXQA12','unknown',' ','formatted',
1914     &               idummy,.false.)
1915      rewind(luvajkl)
1916      read(luvajkl,'(4E30.20)') (work(sajkl+i-1), i = 1,
1917     &                        nvir0t*nrhfst**3)
1918      call gpclose(luvajkl,'DELETE')
1919
1920
1921
1922      call dzero(work(kspjkl),ntklpj)
1923      call dzero(work(ksijkl),nrhfst**4)
1924      call cc_r12xsort1(work(smujkl),work(kspjkl))
1925      call cc_r12xsort(work(sajkl),work(ksajkl),nvir0)
1926      call cc_r12xsort4(work(smujkl),work(ksijkl))
1927
1928      if (froimp) then
1929         call dzero(work(ksijklf),nrhfst**4)
1930         idx = 0
1931         do isymkl = 1, nsym
1932            isymij  = isymkl
1933            do isyml =1, nsym
1934               isymk = muld2h(isymkl,isyml)
1935              do k = nrhffr(isymk)+1, nrhfs(isymk)
1936                 do l = nrhffr(isyml)+1, nrhfs(isyml)
1937                    idxkl=imatij(isymk,isyml)+nrhfs(isymk)*(l-1)+k
1938                    do isymi =1, nsym
1939                       isymj  = muld2h(isymij,isymi)
1940                       do j = nrhffr(isymj)+1, nrhfs(isymj)
1941                          do i = 1, nrhffr(isymi)
1942                            idx = idx + 1
1943                             kl = imatmn(isymk,isyml)+
1944     &                             (l-1)*nrhfs(isymk) + k
1945                             ij = imatmn(isymi,isymj)+
1946     &                             (j-1)*nrhfs(isymi) + i
1947                             ijkl = imatoi(isymij,isymkl)+
1948     &                             (kl-1)*nmatkl1(isymij) + ij
1949                             work(ksijklf+idx-1) = work(ksijkl+ijkl-1)
1950                         enddo
1951                       enddo
1952                    enddo
1953                 enddo
1954              enddo
1955            enddo
1956         enddo
1957
1958
1959
1960         call dzero(work(ksajklf),ntklaj)
1961         idx = 0
1962         do isymkl = 1, nsym
1963            isymaj  = isymkl
1964            isymij  = isymkl
1965            do isyml =1, nsym
1966               isymk = muld2h(isymkl,isyml)
1967              do k = nrhffr(isymk)+1, nrhfs(isymk)
1968                do l = nrhffr(isyml)+1, nrhfs(isyml)
1969                   idxkl=imatij(isymk,isyml)+nrhfs(isymk)*(l-1)+k
1970                    do isyma =1, nsym
1971                       isymj  = muld2h(isymaj,isyma)
1972                       nvir0(isyma)  = norb1(isyma) - nrhfs(isyma)
1973                       do j = nrhffr(isymj)+1, nrhfs(isymj)
1974                          do a = 1, nvir0(isyma)
1975                             idx = idx + 1
1976                             kl =imatmn(isyml,isymk)+
1977     &                            (k-1)*nrhfs(isyml) + l
1978                              aj = imataj(isyma,isymj)+
1979     &                            (j-1)*nvir0(isyma) + a
1980                              ajkl = imatkl1(isymaj,isymkl)+
1981     &                          aj +(kl-1)*nmatan(isymaj)
1982                             work(ksajklf+idx-1) = work(ksajkl+ajkl-1)
1983                          enddo
1984                       enddo
1985                    enddo
1986                 enddo
1987              enddo
1988            enddo
1989         enddo
1990
1991        call dzero(work(kspjklf),nrhfst**4)
1992         idx = 0
1993         do isymkl = 1, nsym
1994            isymaj  = isymkl
1995            isymij  = isymkl
1996            do isyml =1, nsym
1997               isymk = muld2h(isymkl,isyml)
1998              do k = nrhffr(isymk)+1, nrhfs(isymk)
1999                do l = nrhffr(isyml)+1, nrhfs(isyml)
2000                    do isyma =1, nsym
2001                       isymj  = muld2h(isymaj,isyma)
2002                       nvir0(isyma)  = norb1(isyma) - nrhfs(isyma)
2003                       do j = nrhffr(isymj)+1, nrhfs(isymj)
2004                          do a = nrhffr(isyma)+1, nrhfs(isyma)
2005                             idx = idx + 1
2006                             kl =imatmn(isyml,isymk)+
2007     &                            (k-1)*nrhfs(isyml) + l
2008                             aj = imatmn(isyma,isymj)+
2009     &                            (j-1)*nrhfs(isyma) + a
2010                             ajkl = imatoi(isymaj,isymkl)+
2011     &                          aj +(kl-1)*nmatkl1(isymaj)
2012                             work(kspjklf+idx-1) = work(kspjkl+ajkl-1)
2013                          enddo
2014                       enddo
2015                    enddo
2016                 enddo
2017              enddo
2018            enddo
2019         enddo
2020
2021      endif
2022
2023c     Read MO-coefficients
2024
2025      call dzero(work(kcmo),nlamds)
2026      lusifc = -1
2027      call gpopen(lusifc,'SIRIFC','OLD',' ','UNFORMATTED',
2028     &            idummy,.false.)
2029      rewind(lusifc)
2030      call mollab(label,lusifc,lupri)
2031      read(lusifc)
2032      read(lusifc)
2033      read(lusifc) (work(kcmo+i-1),i=1,nlamds)
2034      call gpclose(lusifc,'KEEP')
2035
2036!     reorder to occupied, symmetry; virtual, symmetry
2037!             from symmetry, (occupied,virtual)
2038
2039      call cmo_reorder(work(kcmo),work(kend1),lwrk1)
2040
2041      noffset  = 0
2042      noffset1 = 0
2043      do isym = 1,nsym
2044         noffset  = noffset  + nbas(isym)*nrhf(isym)
2045         noffset1 = noffset1 + nbas(isym)*nrhfs(isym)
2046      enddo
2047
2048c     Transform into orbital basis
2049
2050      call dzero(work(kxpjkl),ntklpj)
2051      kxpjklz = kxpjkl
2052      do isymkl = 1, nsym
2053         isymaj  = isymkl
2054         isymmuj = isymkl
2055         do isyml =1, nsym
2056            isymk = muld2h(isymkl,isyml)
2057            do k = 1, nrhf(isymk)
2058               do l = 1, nrhf(isyml)
2059                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
2060                  do isymmu =1, nsym
2061                     isyma = muld2h(1,isymmu)
2062                     isymj  = muld2h(isymmuj,isymmu)
2063                     koffc = kcmo + noffset+
2064     &                      iglmrv(isymmu,isyma)
2065                     kxajkl = kxmujkl+ivajkl(isymmuj,isymkl)+
2066     &                         nr1bas(isymmuj)*(idxkl-1)
2067     &                         +ir1bas(isymmu,isymj)
2068                     ntotmu = max(1,mbas1(isymmu))
2069                     ntota1 = max(1,nbas(isymmu))
2070                     ntota = max(1,norb1(isyma))
2071                     ntotj = max(1,nrhf(isymj))
2072c
2073                     call dgemm('T','N',norb1(isyma),
2074     &                 nrhf(isymj),mbas1(isymmu),one,work(koffc),
2075     &                 ntota1,work(kxajkl),ntotmu,zero,
2076     &                 work(kxpjkl),ntota)
2077                     kxpjkl = kxpjkl + norb1(isyma)*nrhf(isymj)
2078                  end do
2079               end do
2080            end do
2081         end do
2082      end do
2083
2084
2085      if (froimp) then
2086         ksajkl = ksajklf
2087         kspjkl = kspjklf
2088         ksijkl = ksijklf
2089      endif
2090      ksqijklz = ksqjkl
2091      do isymkl = 1, nsym
2092         isymaj  = isymkl
2093         do isyml =1, nsym
2094            isymk = muld2h(isymkl,isyml)
2095            do k = 1, nrhf(isymk)
2096               do l = 1, nrhf(isyml)
2097                  idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
2098                  do isym =1, nsym
2099                     isyma = muld2h(1,isym)
2100                     isymj  = muld2h(isymkl,isym)
2101                     do i = 1,nrhf(isymj)
2102                        call dcopy(nrhffr(isym),work(ksijklf),1,
2103     &                             work(ksqjkl),1)
2104                        ksqjkl = ksqjkl + nrhffr(isym)
2105                        ksijklf = ksijklf+ nrhffr(isym)
2106                        call dcopy(nrhf(isym),work(kspjkl),1,
2107     &                             work(ksqjkl),1)
2108                        ksqjkl = ksqjkl + nrhf(isym)
2109                        kspjkl = kspjkl + nrhf(isym)
2110                        call dcopy(nvir0(isym),work(ksajkl),1,
2111     &                             work(ksqjkl),1)
2112                        ksqjkl = ksqjkl + nvir0(isym)
2113                        ksajkl = ksajkl + nvir0(isym)
2114                      enddo
2115                  enddo
2116               enddo
2117             enddo
2118          enddo
2119      enddo
2120
2121
2122      call daxpy(ntklpj,one,work(ksqijklz),1,work(kxpjklz),1)
2123c
2124
2125c--------------------------------------------------------------------
2126c     Calculate X_{pj}^{kl}*D_{kl}^{ij}
2127c---------------------------------------------------------------------
2128
2129      call dzero(work(kxdmmu),nr1bast)
2130      do isymkl = 1, nsym
2131         isymaj  = isymkl
2132         isymij  = isymkl
2133         do isyml =1, nsym
2134            isymk = muld2h(isymkl,isyml)
2135            do k = 1, nrhf(isymk)
2136              do l = 1, nrhf(isyml)
2137                 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k
2138                 do isyma =1, nsym
2139                    isymj  = muld2h(isymaj,isyma)
2140                    isymi  = muld2h(isymij,isymj)
2141                    koffd  = 1 + igamsq(isymij,isymkl)
2142     &                 + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj)
2143                    ntota = max(1,norb1(isyma))
2144                    ntotj = max(1,nrhf(isymj))
2145                    ntoti = max(1,nrhf(isymi))
2146                    ntota1= max(1,nbas(isyma))
2147                    lxdmmu= kxdmmu + imatmuj(isyma,isymi)
2148                    call dgemm('N','T',norb1(isyma),nrhf(isymi),
2149     &                    nrhf(isymj),one,work(kxpjklz),
2150     &                    ntota,dklij(koffd),ntoti,one,
2151     &                    work(lxdmmu),ntota)
2152
2153                    kxpjklz = kxpjklz + norb1(isyma)*nrhf(isymj)
2154                    ksqijklz= ksqijklz+ norb1(isyma)*nrhf(isymj)
2155                    ksijkl  = ksijkl+ nrhffr(isyma)*nrhf(isymj)
2156                 enddo
2157              enddo
2158           enddo
2159         enddo
2160      enddo
2161
2162
2163
2164
2165      call dzero(work(kxdpmu),mbas1t*nrhft)
2166      call dzero(work(kxdnumu),10*mbas1t**2)
2167      do isymaj = 1, nsym
2168         isymkl  = isymaj
2169         isymij  = isymaj
2170         noffset1 = 0
2171         do isyma =1, nsym
2172            isymmu = muld2h(1,isyma)
2173            isymj  = muld2h(isymaj,isyma)
2174            isymi  = muld2h(isymij,isymj)
2175            ntota = max(1,mbas1(isymi))
2176            ntotj = max(1,norb1(isyma))
2177            ntoti = max(1,mbas1(isymi))
2178            ntota1= max(1,nbas(isymi))
2179            koffc = kcmo + noffset +
2180     &              iglmrv(isyma,isymi)
2181            lxdmmu = kxdmmu + imatmuj(isyma,isymi)
2182            lxdpmu = kxdpmu + imatmuj(isyma,isymi)
2183            lxdnumu = kxdnumu + imatmunu(isymi,isymi)
2184
2185            call dgemm('N','N',mbas1(isymi),
2186     &           nrhf(isyma),norb1(isyma),one,
2187     &           work(koffc),
2188     &           ntota1,work(lxdmmu),ntotj,zero,
2189     &           work(lxdpmu),ntota)
2190
2191            noffset1 = noffset1 + nbas(isymi)*nrhffr(isymi)
2192            koffd = kcmo + noffset + noffset1 +
2193     &              iglmrv1(isyma,isymi)
2194
2195            call dgemm('N','T',mbas1(isymi),
2196     &           mbas1(isymi),nrhf(isyma),one,work(koffd),
2197     &           ntota1,work(lxdpmu),ntoti,zero,
2198     &           work(lxdnumu),ntota)
2199
2200         enddo
2201      enddo
2202
2203
2204
2205      call dzero(work(kzmunu),n2basx)
2206      call dzero(work(kdens),n2basx)
2207      koff = 0
2208      do isymij = 1,nsym
2209         isymj = isymij
2210         isymi = isymij
2211         koff = imatmunu(isymi,isymj)
2212         do j = 1,mbas1(isymj)
2213            do i = 1,mbas1(isymi)
2214               koff = koff+1
2215               idxij = iaodis(isymi,isymj)+nbas(isymi)*
2216     &                 (j-1)+i
2217               work(kzmunu+idxij-1) = work(kxdnumu+koff-1)
2218             enddo
2219          enddo
2220      enddo
2221
2222
2223      call cc_r12fsym(imatmunu,work(kzmunu),work(kdens))
2224
2225
2226      call r12exdens(res,work(kcmo),work(kdens),work(kend1),lwrk1)
2227
2228
2229      call qexit('r12gradbd')
2230      end
2231c------------------------------------------------------------------------
2232C=====================================================================*
2233C  /* Deck r12exdens*/
2234      subroutine r12exdens(res,cmo,dens,work,lwork)
2235      implicit none
2236#include "ccsdsym.h"
2237#include "ccorb.h"
2238#include "dummy.h"
2239#include "r12int.h"
2240#include "priunit.h"
2241#include "ccfro.h"
2242#include "inftap.h"
2243#include "ccsdinp.h"
2244ccc#include "inforb.h"
2245      integer kcmo,lwork,kfrsav,kfree,lfree,kfvao,kfcao,kdcao,kdvao,kdv
2246      integer pq,kcao,isym,ioff,ioff1,kexia,isymi,NNORB1
2247      integer kctao,koff,index,ISYMAI,isyma,nnbast
2248      integer idxkl,idxlk,idxij,idxji,idxklij,idxijkl,idxijlk
2249      integer idxjikl,isymj,isymij,isymkl,isymk,isyml,idxklji
2250      integer koffc,kexiaz,noffset,nvir0(8),kextia,kfctao
2251      integer icoun1,imatpq(8,8),icoun2,imatsq(8,8),ij,ji
2252      integer kajfr,isymaj,luvajkl,ncfai,kexfa,kexfaz
2253      logical onlyfc
2254c
2255      integer icoun3,imatfq(8,8),kxij,kxijij,ncvij,nvirf(8)
2256#include "maxorb.h"
2257#include "infinp.h"
2258      real*8 res(*),cmo(*),emcmy,zero,one,two,work(*),dens(*),densx
2259      parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0)
2260
2261      index(i,j) = max(i,j)*(max(i,j)-3)/2 +i+j
2262
2263      call qenter('r12exdens')
2264
2265
2266      nrhft = 0
2267      do isymi = 1, nsym
2268         nvir0(isymi) = norb1(isymi) - nrhfs(isymi)
2269         nvirf(isymi) = nvir(isymi) - nrhffr(isymi)
2270         nrhft        = nrhft + nrhf(isymi)
2271      end do
2272
2273      do isymaj = 1,nsym
2274         isymij = isymaj
2275         ncfai = 0
2276         do isyma = 1,nsym
2277            isymj = muld2h(isymaj,isyma)
2278            isymi = muld2h(isymij,isymj)
2279            ncfai = ncfai + (norb1(isyma) - nrhfs(isyma))*nrhffr(isymi)
2280         enddo
2281      enddo
2282
2283
2284      do isymi = 1,nsym
2285         icoun1 = 0
2286         icoun2 = 0
2287         icoun3 = 0
2288         do isymj = 1,nsym
2289            isymk = muld2h(isymj,isymi)
2290            imatpq(isymk,isymj) = icoun1
2291            imatsq(isymk,isymj) = icoun2
2292            imatfq(isymk,isymj) = icoun3
2293            icoun1 = icoun1 + nvir(isymk)*(nvir(isymj)+1)/2
2294            icoun2 = icoun2 + nrhf(isymk)*nvir0(isymj)
2295            icoun3 = icoun3 + nrhffr(isymk)*nvir0(isymj)
2296         enddo
2297      enddo
2298
2299      do isymaj = 1,nsym
2300         isymij = isymaj
2301         ncvij = 0
2302         do isyma = 1,nsym
2303            isymj = muld2h(isymaj,isyma)
2304            isymi = muld2h(isymij,isymj)
2305            ncvij = ncvij +  nrhffr(isyma)*nrhf(isymi)
2306         enddo
2307      enddo
2308
2309
2310
2311
2312      NBAST  = 0
2313      NNBAST = 0
2314      DO I = 1,NSYM
2315         NBAST  = NBAST  + NBAS(I)
2316         NNBAST = NNBAST + NBAS(I)*(NBAS(I)+1)/2
2317      ENDDO
2318
2319      noffset = 0
2320      do isym = 1,nsym
2321         noffset = noffset+nbas(isym)*nrhf(isym)
2322      enddo
2323
2324      kfrsav = 1
2325      KFREE  = 1
2326      LFREE  = lwork
2327      CALL MEMGET('REAL',kfcao,N2BASX,WORK,KFREE,LFREE)
2328      CALL MEMGET('REAL',kfvao,     0,WORK,KFREE,LFREE)
2329      CALL MEMGET('REAL',kdcao,N2BASX,WORK,KFREE,LFREE)
2330      CALL MEMGET('REAL',kdvao,     0,WORK,KFREE,LFREE)
2331      CALL MEMGET('REAL',kdv  ,     0,WORK,KFREE,LFREE)
2332      CALL MEMGET('REAL',kcmo, nlamds,WORK,KFREE,LFREE)
2333      CALL MEMGET('REAL',kexia,nnbast,WORK,KFREE,LFREE)
2334      CALL MEMGET('REAL',kexfa,nnbast,WORK,KFREE,LFREE)
2335      CALL MEMGET('REAL',kajfr,ncfai ,WORK,KFREE,LFREE)
2336      CALL MEMGET('REAL',kxijij,ncvij,WORK,KFREE,LFREE)
2337      CALL MEMGET('REAL',kxij  ,ncvij,WORK,KFREE,LFREE)
2338
2339      call erisfk(dens,nbast,1)
2340      call dscal(n2basx,0.50d0,dens,1)
2341
2342
2343      DIRFCK = .TRUE.
2344      R12TRA = .TRUE.
2345      R12INT = .FALSE.
2346      V12INT = .TRUE.
2347      NORXR = NORXR .OR. R12HYB
2348      mbsmax = 4
2349
2350      LUSUPM = -1
2351      call dzero(work(kfcao),n2basx)
2352      onlyfc = .TRUE.
2353      call fckmao(onlyfc,emcmy,work(kfcao),work(kfvao),
2354     &     dens,work(kdvao),work(kdv),cmo(1+noffset),work(kfree),lfree)
2355
2356!     Note: dens is destroyed in fckmao because of dirfck, but that doesn't
2357!     matter as it is not used any more /Jan 2011
2358
2359c
2360      WRITE(LUPRI,'(/A,I1,A)')
2361     &   ' Computation of exchange matrix done [',MBSMAX,']'
2362
2363
2364      R12INT = .TRUE.
2365      V12INT = .FALSE.
2366      call dzero(WORK(KDCAO),N2BASX)
2367      call dzero(WORK(kexia),NNBAST)
2368      CALL DCOPY(N2BASX,WORK(KFCAO),1,WORK(KDCAO),1)
2369c N2BASX = NBAST*NBAST!
2370      CALL DGETSP(NBAST,WORK(KDCAO),WORK(KFCAO))
2371      IF (NSYM .GT. 1) CALL PKSYM1(WORK(KFCAO),WORK(KFCAO),NBAS,NSYM,2)
2372      CALL DZERO(WORK(kexia),NNBAST)
2373      kexiaz = kexia
2374      koffc= 1
2375      CALL UTHUB(WORK(KFCAO),WORK(kexia),
2376     *           CMO(1+noffset),WORK(kfree),NSYM,nbas,nvir)
2377
2378
2379      koff = 0
2380      do isym = 1,nsym
2381         do p = 1+nrhfs(isym),norb1(isym)
2382            do q = 1+nrhffr(isym),nrhfs(isym)
2383               koff = imatsq(isym,isym)+nvir0(isym)*(q-nrhffr(isym)-1)
2384     &                        +p-nrhfs(isym)
2385               pq = imatpq(isym,isym)+index(p,q)
2386               res(koff) = 4.00d0*work(kexia+pq-1)
2387            enddo
2388         enddo
2389      enddo
2390
2391      if (froimp) then
2392
2393          koff = 0
2394          do isym = 1,nsym
2395             do p = 1+nrhfs(isym),norb1(isym)
2396                do q = 1,nrhffr(isym)
2397                   koff = imatfq(isym,isym)+nvir0(isym)
2398     &                       *(q-1)+p - nrhfs(isym)
2399                   pq = imatpq(isym,isym)+index(p,q)
2400                   work(kajfr+koff-1) = 8.00d0*work(kexia+pq-1)
2401                enddo
2402            enddo
2403          enddo
2404
2405          luvajkl = -1
2406          call gpopen(luvajkl,'CCR12YAIFR','unknown',' ','unformatted',
2407     &                idummy,.false.)
2408          write(luvajkl) (work(kajfr+i-1), i = 1,ncfai)
2409          call gpclose(luvajkl,'KEEP')
2410
2411
2412      endif
2413
2414
2415      CALL MEMREL('r12exdens',WORK,KFRSAV,KFRSAV,KFREE,LFREE)
2416      call qexit('r12exdens')
2417      end
2418c--------------------------------------------------------------------------
2419