1!
2!  Dalton, a molecular electronic structure program
3!  Copyright (C) by the authors of Dalton.
4!
5!  This program is free software; you can redistribute it and/or
6!  modify it under the terms of the GNU Lesser General Public
7!  License version 2.1 as published by the Free Software Foundation.
8!
9!  This program is distributed in the hope that it will be useful,
10!  but WITHOUT ANY WARRANTY; without even the implied warranty of
11!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12!  Lesser General Public License for more details.
13!
14!  If a copy of the GNU LGPL v2.1 was not distributed with this
15!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
16!
17!
18C
19*=====================================================================*
20      SUBROUTINE CC_QMAT(QMATP, QMATH, RMAT, XKAPPA,
21     &                   IREAL, ISYMQ, NOKAPPA, CMO, WORK, LWORK )
22*---------------------------------------------------------------------*
23*
24*     Purpose: build the Q matrices, which are basically the connection
25*              matrix transformed to the MO representation minus
26*              the orbital relaxation matrix:
27*
28*              Q^h :  R   - kappa
29*              Q^p :  R^* + kappa^T
30*
31*              QMATP  : Q^p matrices
32*              QMATH  : Q^h matrices
33*              RMAT   : orbital connection matrix in AO basis
34*              XKAPPA : orbital relaxation vector (packed)
35*              CMO    : orbital coefficient matrix, Sirius ordering
36*
37*              NOKAPPA = .TRUE.  --> skip contribution of kappa
38*                                  (used f.x. for SCF kappa rhs vectors)
39*
40*     Christof Haettig, March 1999
41*     generalized for non-symmetric R and Kappa in November 1999
42*
43*=====================================================================*
44#if defined (IMPLICIT_NONE)
45      IMPLICIT NONE
46#else
47#  include "implicit.h"
48#endif
49#include "priunit.h"
50#include "ccorb.h"
51#include "ccfro.h"
52#include "ccsdsym.h"
53
54      LOGICAL LOCDBG
55      PARAMETER (LOCDBG = .FALSE.)
56
57      INTEGER ISYM0
58      PARAMETER (ISYM0 = 1)
59
60      LOGICAL NOKAPPA
61      INTEGER IREAL, ISYMQ, LWORK
62
63#if defined (SYS_CRAY)
64      REAL QMATP(*), QMATH(*), RMAT(*), XKAPPA(*), CMO(*)
65      REAL WORK(LWORK), ONE, ZERO, SIGN
66#else
67      DOUBLE PRECISION QMATP(*), QMATH(*), RMAT(*), XKAPPA(*), CMO(*)
68      DOUBLE PRECISION WORK(LWORK), ONE, ZERO, SIGN
69#endif
70      PARAMETER(ONE=1.0D0, ZERO=0.0D0)
71
72      INTEGER ISYALP, ISYBET, NBASA, NBASB, NORBSA, IORBI, IORBA
73      INTEGER KOFF1,KOFF2,KOFF3,KOFF4,KKAI,KRAI,KRIA,KSCR1,KSCR2,KEND
74      INTEGER NCMO(8), ICMO(8,8), ISYM, ICOUNT, ISYM2, ISYM1, LEN
75
76*---------------------------------------------------------------------*
77*     set ICMO & NCMO arrays:
78*---------------------------------------------------------------------*
79      DO ISYM = 1, NSYM
80         ICOUNT = 0
81         DO ISYM2 = 1, NSYM
82            ISYM1 = MULD2H(ISYM,ISYM2)
83            ICMO(ISYM1,ISYM2) = ICOUNT
84            ICOUNT = ICOUNT + NBAS(ISYM1)*NORBS(ISYM2)
85         END DO
86         NCMO(ISYM) = ICOUNT
87      END DO
88
89*---------------------------------------------------------------------*
90*     put the orbital relaxation (kappa) vector into the Q matrices:
91*---------------------------------------------------------------------*
92      IF (.NOT. NOKAPPA) THEN
93        CALL CCKAPPASQ(QMATH,XKAPPA,ISYMQ,'N')
94        CALL CCKAPPASQ(QMATP,XKAPPA,ISYMQ,'T')
95      ELSE
96        CALL DZERO(QMATH,N2BST(ISYMQ))
97        CALL DZERO(QMATP,N2BST(ISYMQ))
98      END IF
99
100      IF (LOCDBG) THEN
101        WRITE (LUPRI,*) 'debug CC_QMAT> NOKAPPA :',NOKAPPA
102        WRITE (LUPRI,*) 'debug CC_QMAT> IREAL   :',IREAL
103        WRITE (LUPRI,*) 'debug CC_QMAT> kappa matrix in MO:'
104        CALL CC_PRONELAO(QMATH,ISYMQ)
105        WRITE (LUPRI,*) 'debug CC_QMAT> kappa^T matrix in MO:'
106        CALL CC_PRONELAO(QMATP,ISYMQ)
107        WRITE (LUPRI,*) 'debug CC_QMAT> R matrix in AO:'
108        CALL CC_PRONELAO(RMAT,ISYMQ)
109      END IF
110
111*---------------------------------------------------------------------*
112*     transform the connection matrix R to MO using the CMO vector,
113*     which is not resorted to the CC standard ordering.
114*     put the result into the Q matrices.
115*
116*     THE FOLLOWING HAS TO BE COUNTERCHECKED!
117*     WHAT HAPPENS FOR FROZEN AND DELETED(!) ORBITALS!
118*---------------------------------------------------------------------*
119      DO ISYALP = 1, NSYM
120         ISYBET = MULD2H(ISYALP,ISYMQ)
121
122         KSCR1  = 1
123         KSCR2  = KSCR1 +  NBAS(ISYALP)*NORBS(ISYBET)
124         KEND   = KSCR2 + NORBS(ISYALP)*NORBS(ISYBET)
125
126         IF ( KEND .GT. LWORK ) THEN
127            CALL QUIT('Insufficient memory in CC_QMAT.')
128         END IF
129
130         NBASA  = MAX(NBAS(ISYALP),1)
131         NBASB  = MAX(NBAS(ISYBET),1)
132         NORBSA = MAX(NORBS(ISYALP),1)
133
134         ! transform second index of R(alpha,beta) to MO
135         KOFF1 = IAODIS(ISYALP,ISYBET) + 1
136         KOFF2 = ICMO(ISYBET,ISYBET)   + 1
137         CALL DGEMM('N','N',NBAS(ISYALP),NORBS(ISYBET),NBAS(ISYBET),
138     &              ONE,RMAT(KOFF1),NBASA,CMO(KOFF2),NBASB,
139     &              ZERO,WORK(KSCR1),NBASA)
140
141
142         ! transform leading index of R(alpha,beta) to MO
143         KOFF3 = ICMO(ISYALP,ISYALP)   + 1
144         CALL DGEMM('T','N',NORBS(ISYALP),NORBS(ISYBET),NBAS(ISYALP),
145     &              ONE,CMO(KOFF3),NBASA,WORK(KSCR1),NBASA,
146     &              ZERO,WORK(KSCR2),NORBSA)
147
148         LEN   = NORBS(ISYALP) * NORBS(ISYBET)
149         KOFF4 = IAODIS(ISYALP,ISYBET) + 1
150
151         ! Q^p(a,b) = R(a,b) - kappa(a,b)
152         CALL DSCAL(LEN,-ONE,QMATH(KOFF4),1)
153         CALL DAXPY(LEN,+ONE,WORK(KSCR2),1,QMATH(KOFF4),1)
154
155
156         ! Q^p(a,b) = R^*(a,b) + kappa^T(a,b)
157         SIGN = DBLE(IREAL)
158         CALL DAXPY(LEN,SIGN,WORK(KSCR2),1,QMATP(KOFF4),1)
159
160      END DO
161
162
163      IF (LOCDBG) THEN
164        WRITE (LUPRI,*) 'debug CC_QMAT> Q^h matrix in MO:'
165        CALL CC_PRONELAO(QMATH,ISYMQ)
166        WRITE (LUPRI,*) 'debug CC_QMAT> Q^p matrix in MO:'
167        CALL CC_PRONELAO(QMATP,ISYMQ)
168      END IF
169
170      RETURN
171      END
172*=====================================================================*
173