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_FIND_SO_OP(LABELA,LABELB,LABSOP,ISYSOP,ISIGN,INUM,
21     &                         WORK,LWORK)
22*---------------------------------------------------------------------*
23*
24*     Purpose: find second-order operator from the labels of the
25*              corresponding first-order perturbations and try to
26*              read the operator from the AONEINT file
27*              if the labels is identified and the operator is found
28*              on file, it is added to the IROPER2 list
29*
30*
31*     Input:   LABELA,LABELB  --  labels of first-order perturbations
32*
33*     Output:  LABSOP    --  labels for second-order operator
34*              ISYSOP    --  point group symmetry
35*              INUM      --  index on IROPER2 list, -1 if not found
36*              ISIGN     --  flag for sign, since for some operators
37*                            the sign conventions differs from the
38*                            one for the first-order operators
39*
40*     Christof Haettig, March 1999
41*
42*=====================================================================*
43#if defined (IMPLICIT_NONE)
44      IMPLICIT NONE
45#else
46#  include "implicit.h"
47#endif
48#include "priunit.h"
49#include "ccorb.h"
50#include "ccsdsym.h"
51#include "ccroper.h"
52#include "ccropr2.h"
53
54      LOGICAL LOCDBG
55      PARAMETER (LOCDBG = .FALSE.)
56
57      CHARACTER*8 LABELA, LABELB, LABSOP
58      LOGICAL FOUND, LOPNSAVE
59      INTEGER LWORK, ISYSOP, INUM, IMATRIX, IERR, KPROPAO
60      INTEGER KEND1, LWRK1, ISIGN
61      INTEGER ISYMA, ISYMB, INUMA, INUMB, INUMAB
62
63      REAL*8   WORK(LWORK)
64
65* external functions:
66      INTEGER IROPER2
67      INTEGER IROPER
68
69*---------------------------------------------------------------------*
70* check if already known (if so, return without any further action):
71*---------------------------------------------------------------------*
72      LOPNSAVE = LOPR2OPN
73      LOPR2OPN = .FALSE.
74      LQUIET   = .TRUE.
75
76      INUM = IROPER2(LABELA,LABELB,LABSOP,ISIGN,ISYSOP)
77
78      LQUIET   = .FALSE.
79      LOPR2OPN = LOPNSAVE
80
81      IF (INUM.GT.0) RETURN
82
83*---------------------------------------------------------------------*
84* for some second-order operators the sign convention differs from the
85* one for the first-order operators and we need to turn the sign
86* (see below), but default is not to turn the sign.
87*---------------------------------------------------------------------*
88      ISIGN = +1
89
90*---------------------------------------------------------------------*
91* Hessian (geometric second derivatives): not yet available
92*---------------------------------------------------------------------*
93      IF    (LABELA(1:5).EQ.'1DHAM'.AND.LABELB(1:5).EQ.'1DHAM') THEN
94
95        FOUND = .FALSE.
96
97*---------------------------------------------------------------------*
98* Dipole gradient:
99*---------------------------------------------------------------------*
100      ELSEIF(LABELA(2:8).EQ.'DIPLEN '.AND.LABELB(1:5).EQ.'1DHAM') THEN
101
102        WRITE(LABSOP,'(A3,A4,A1)') LABELB(6:8), 'DPG ', LABELA(1:1)
103        FOUND = .TRUE.
104        ISIGN = -1
105
106      ELSEIF(LABELB(2:8).EQ.'DIPLEN '.AND.LABELA(1:5).EQ.'1DHAM') THEN
107
108        WRITE(LABSOP,'(A3,A4,A1)') LABELA(6:8), 'DPG ', LABELB(1:1)
109        FOUND = .TRUE.
110        ISIGN = -1
111
112*---------------------------------------------------------------------*
113* Second moment of charge:
114*---------------------------------------------------------------------*
115      ELSEIF(LABELA(3:8).EQ.'SECMOM'.AND.LABELB(1:5).EQ.'1DHAM') THEN
116
117        WRITE(LABSOP,'(A3,A3,A2)') LABELB(6:8), 'QDG', LABELA(1:2)
118        FOUND = .TRUE.
119        ISIGN = -1
120
121      ELSEIF(LABELB(3:8).EQ.'SECMOM'.AND.LABELA(1:5).EQ.'1DHAM') THEN
122
123        WRITE(LABSOP,'(A3,A3,A2)') LABELA(6:8), 'QDG', LABELB(1:2)
124        FOUND = .TRUE.
125        ISIGN = -1
126
127*---------------------------------------------------------------------*
128* Third moment of charge:
129*---------------------------------------------------------------------*
130      ELSEIF(LABELA(5:8).EQ.'3MOM'.AND.LABELB(1:5).EQ.'1DHAM') THEN
131
132        WRITE(LABSOP,'(A2,A3,A3)') LABELB(7:8), 'OCG', LABELA(1:3)
133        FOUND = .TRUE.
134        ISIGN = -1
135
136      ELSEIF(LABELB(5:8).EQ.'3MOM'.AND.LABELA(1:5).EQ.'1DHAM') THEN
137
138        WRITE(LABSOP,'(A2,A3,A3)') LABELA(7:8), 'OCG', LABELB(1:3)
139        FOUND = .TRUE.
140        ISIGN = -1
141
142*---------------------------------------------------------------------*
143* Quadrupole gradient: not yet available
144*---------------------------------------------------------------------*
145
146*---------------------------------------------------------------------*
147* Octupole gradient: not yet available
148*---------------------------------------------------------------------*
149
150*---------------------------------------------------------------------*
151* Magnetizabilities (magnetic second derivatives): not yet available
152*---------------------------------------------------------------------*
153      ELSEIF(LABELA(1:5).EQ.'dh/dB'.AND.LABELB(1:5).EQ.'dh/dB') THEN
154
155        FOUND = .FALSE.
156
157*---------------------------------------------------------------------*
158* geometric derivatives of magnetic properties:
159*---------------------------------------------------------------------*
160      ELSEIF(LABELA(1:5).EQ.'1DHAM'.AND.LABELB(1:5).EQ.'dh/dB') THEN
161
162        FOUND = .FALSE.
163
164      ELSEIF(LABELA(1:5).EQ.'dh/dB'.AND.LABELB(1:5).EQ.'1DHAM') THEN
165
166        FOUND = .FALSE.
167
168*---------------------------------------------------------------------*
169* mixed electric/magnetic properties with London orbitals:
170*---------------------------------------------------------------------*
171      ELSEIF(LABELA(2:7).EQ.'DIPLEN'.AND.LABELB(1:5).EQ.'dh/dB') THEN
172
173        WRITE(LABSOP,'(A1,A5,A2)') LABELA(1:1), '-CM1 ', LABELB(6:7)
174        FOUND = .TRUE.
175        ISIGN = +1
176
177      ELSEIF(LABELA(1:5).EQ.'dh/dB'.AND.LABELB(2:7).EQ.'DIPLEN') THEN
178
179        WRITE(LABSOP,'(A1,A5,A2)') LABELB(1:1), '-CM1 ', LABELA(6:7)
180        FOUND = .TRUE.
181        ISIGN = +1
182
183*---------------------------------------------------------------------*
184* mixed electric/magnetic properties with usual (not London) orbitals:
185*---------------------------------------------------------------------*
186      ELSEIF(LABELA(2:7).EQ.'DIPLEN'.AND.LABELB(2:7).EQ.'ANGMOM') THEN
187
188        WRITE(LABSOP,'(A6,A1,A1)') '-> DxL',LABELA(1:1),LABELB(1:1)
189        FOUND = .TRUE.
190        ISIGN = 0
191
192      ELSEIF(LABELA(2:7).EQ.'ANGMOM'.AND.LABELB(2:7).EQ.'DIPLEN') THEN
193
194        WRITE(LABSOP,'(A6,A1,A1)') '-> DxL',LABELB(1:1),LABELA(1:1)
195        FOUND = .TRUE.
196        ISIGN = 0
197
198*---------------------------------------------------------------------*
199* two dipole operators --> second-order operator is zero
200*---------------------------------------------------------------------*
201      ELSEIF(LABELA(2:7).EQ.'DIPLEN'.AND.LABELB(2:7).EQ.'DIPLEN') THEN
202
203        WRITE(LABSOP,'(A6,A1,A1)') '-> DxD',LABELA(1:1),LABELB(1:1)
204        FOUND = .TRUE.
205        ISIGN = 0
206
207*---------------------------------------------------------------------*
208* nuclear shieldings:
209*---------------------------------------------------------------------*
210      ELSEIF(LABELA(1:5).EQ.'dh/dB'.AND.LABELB(1:4).EQ.'PSO ') THEN
211
212        WRITE(LABSOP,'(A3,A4,A1)') LABELB(5:7),' NST',LABELA(6:6)
213        FOUND = .TRUE.
214        ISIGN = +1
215
216      ELSEIF(LABELB(1:5).EQ.'dh/dB'.AND.LABELA(1:4).EQ.'PSO ') THEN
217
218        WRITE(LABSOP,'(A3,A4,A1)') LABELA(5:7),' NST',LABELB(6:6)
219        FOUND = .TRUE.
220        ISIGN = +1
221
222*---------------------------------------------------------------------*
223* default: no second-order operator in the Hamiltonian
224*---------------------------------------------------------------------*
225      ELSE
226        FOUND = .FALSE.
227      END IF
228
229*---------------------------------------------------------------------*
230* check if the LABSOP is available on the AONEINT file:
231*---------------------------------------------------------------------*
232      IF (FOUND .AND. LABSOP(1:2).NE.'->') THEN
233
234         KPROPAO = 1
235         KEND1   = KPROPAO + N2BASX
236         LWRK1   = LWORK   - KEND1
237
238         CALL CCPRPAO(LABSOP,.TRUE.,WORK(KPROPAO),ISYSOP,IMATRIX,IERR,
239     &                WORK(KEND1),LWRK1)
240
241         IF (IERR.GT.0) FOUND = .FALSE.
242
243         IF (IERR.LT.0) THEN
244            INUMA   = IROPER(LABELA,ISYMA)
245            INUMB   = IROPER(LABELB,ISYMB)
246            ISYSOP  = MULD2H(ISYMA,ISYMB)
247            IMATRIX = ISYMAT(INUMA) * ISYMAT(INUMB)
248            IF (LOCDBG) THEN
249              WRITE(LUPRI,*) 'LABELA,ISYMA,INUMA:',LABELA,ISYMA,INUMA
250              WRITE(LUPRI,*) 'LABELB,ISYMB,INUMB:',LABELB,ISYMB,INUMB
251              WRITE(LUPRI,*) 'IMAT:',ISYMAT(INUMA),ISYMAT(INUMB)
252            END IF
253         END IF
254
255         IF (LOCDBG) THEN
256            WRITE (LUPRI,*) '"',LABSOP,'" integrals found on file :',
257     &           FOUND
258         END IF
259
260      END IF
261
262      IF (FOUND .AND. LABSOP(1:2).EQ.'->') THEN
263         INUMA   = IROPER(LABELA,ISYMA)
264         INUMB   = IROPER(LABELB,ISYMB)
265         ISYSOP  = MULD2H(ISYMA,ISYMB)
266         IMATRIX = 0
267      END IF
268
269*---------------------------------------------------------------------*
270* add operator to the IROPER2 list:
271*---------------------------------------------------------------------*
272      IF (FOUND) THEN
273         LOPNSAVE = LOPR2OPN
274         LOPR2OPN = .TRUE.
275
276         INUM = IROPER2(LABELA,LABELB,LABSOP,ISIGN,ISYSOP)
277
278         ! save symmetry of integral matrix on common blocks
279         ISYMAT2(INUM) = IMATRIX
280
281         ! get index of operator on IROPER list and
282         ! set symmetry of integral matrix correct for this list:
283         INUMAB = IROPER(LABSOP,ISYSOP)
284         ISYMAT(INUMAB) = IMATRIX
285
286         LOPR2OPN = LOPNSAVE
287      ELSE
288         INUM = -1
289      END IF
290
291      IF (LOCDBG) THEN
292        WRITE (LUPRI,*) 'CC_FIND_SO_OP> LABELA :',LABELA
293        WRITE (LUPRI,*) 'CC_FIND_SO_OP> LABELB :',LABELB
294        WRITE (LUPRI,*) 'CC_FIND_SO_OP> LABSOP :',LABSOP
295        WRITE (LUPRI,*) 'CC_FIND_SO_OP> ISYSOP :',ISYSOP
296        WRITE (LUPRI,*) 'CC_FIND_SO_OP> IMATRIX:',IMATRIX
297      END IF
298
299      RETURN
300      END
301*=====================================================================*
302