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 CCX_GRAD(WORK,LWORK)
21*---------------------------------------------------------------------*
22*
23*    Purpose: calculate relaxed excited state expectation values
24*             of one- and two- electron operators the dummy way,
25*             i.e., setting up the O1 vectors and some other vectors
26*             in the MO basis and calculating LE A^(1) RE.
27*
28*
29*     Written by Christof Haettig, July 2002
30*
31*=====================================================================*
32#if defined (IMPLICIT_NONE)
33      IMPLICIT NONE
34#else
35#  include "implicit.h"
36#endif
37#include "priunit.h"
38#include "ccexgr.h"
39#include "ccgr.h"
40#include "maxorb.h"
41#include "cclists.h"
42#include "mxcent.h"
43#include "nuclei.h"
44#include "energy.h"
45#include "ccroper.h"
46#include "ccr1rsp.h"
47#include "cco1rsp.h"
48#include "ccsdinp.h"
49
50      LOGICAL LOCDBG
51      PARAMETER (LOCDBG = .FALSE.)
52
53      INTEGER LWORK
54      INTEGER IXETRAN(MXDIM_XEVEC,NAXGRO*NXGRST)
55      INTEGER IOTRAN(NXGRST)
56      INTEGER IEDOTS(NAXGRO*NXGRST), IODOTS(NAXGRO,NXGRST)
57
58      REAL*8 WORK(LWORK)
59      REAL*8 ECONS(NAXGRO*NXGRST), OCONS(NAXGRO,NXGRST)
60      REAL*8 RDUM, PROPT, PROPE, PROPN, GRDNRM
61
62      CHARACTER CDUM*1
63      INTEGER NXETRAN, IOP, IOPER, IRELAX, IOPTRES, IORDER, IDUM, IX,
64     &        IDXO1, ISYST, IEX, ISCOOR
65
66      INTEGER IR1TAMP, IRHSR1, ILSTSYM, IDXSYM
67      REAL*8 CC_NUCCON, DDOT
68
69! trkoor.h : NCOOR
70#include "trkoor.h"
71      REAL*8 ERGMOL
72      REAL*8, allocatable :: GRDMOL(:), HESMOL(:,:)
73
74      CALL QENTER('CCX_GRAD')
75
76      IF (LOCDBG) WRITE(LUPRI,*) 'entered CCX_GRAD...'
77
78
79      DO IX = 1, NXGRST
80        IOTRAN(IX) = IX
81        DO IOP = 1, NAXGRO
82          IODOTS(IOP,IX) = 0
83          OCONS(IOP,IX)  = 0.0D0
84        END DO
85      END DO
86
87      NXETRAN = 0
88      DO IOP = 1, NAXGRO
89        IOPER = IAXGRO(IOP)
90        IF (LPDBSOP(IOPER) .AND. ISYOPR(IOPER).EQ.1) THEN
91         IRELAX = IR1TAMP(LBLOPR(IOPER),.TRUE.,0.0D0,ISYOPR(IOPER))
92         IDXO1  = IRHSR1(LBLOPR(IOPER),.TRUE.,0.0D0,ISYOPR(IOPER))
93cch
94          write(lupri,*) 'ccx_grad> IOP,IOPER,LBLOPR(IOPER):',
95     &                              IOP,IOPER,LBLOPR(IOPER)
96          write(lupri,*) 'ccx_grad> IRELAX:',IRELAX
97          write(lupri,*) 'ccx_grad> IDXO1 :',IDXO1
98cch
99         DO IX = 1, NXGRST
100cch
101          write(lupri,*) 'ccx_grad> IX,IXGRST(IX):',IX,IXGRST(IX)
102cch
103          NXETRAN = NXETRAN + 1
104          IXETRAN(1,NXETRAN) = IOPER
105          IXETRAN(2,NXETRAN) = IXGRST(IX)
106          IXETRAN(3,NXETRAN) = -1
107          IXETRAN(4,NXETRAN) =  0
108          IXETRAN(5,NXETRAN) = IRELAX
109          IXETRAN(6,NXETRAN) = -1
110          IXETRAN(7,NXETRAN) = -1
111          IXETRAN(8,NXETRAN) = -1
112          IEDOTS(NXETRAN)    = IXGRST(IX)
113          IOTRAN(IX)         = IX
114          IODOTS(IOP,IX)     = IDXO1
115         END DO
116        END IF
117      END DO
118
119      ! check if anything at all to do...
120      IF (NXETRAN.EQ.0) THEN
121        CALL QEXIT('CCX_GRAD')
122        RETURN
123      END IF
124
125      CALL AROUND(' Relaxed and/or two-electron perturbations ')
126      WRITE(LUPRI,'(/10x,A)')
127     & 'Operator   Total       Electronic  Nuclear'
128
129      IF (.NOT.(CCS.OR.CC2.OR.CCSD)) THEN
130        WRITE(LUPRI,*) 'Requested excited states properties not     '
131        WRITE(LUPRI,*) 'available for the current wavefunction model'
132      END IF
133
134      IF (LOCDBG) WRITE(LUPRI,*) 'call cc_xieta...'
135
136      IF (XGROPT) THEN
137        ! set all gradient components to zero
138        CALL DZERO(GRADKE,MXCOOR)
139        CALL DZERO(GRADNA,MXCOOR)
140        CALL DZERO(GRADEE,MXCOOR)
141        CALL DZERO(GRADNN,MXCOOR)
142        CALL DZERO(GRADFS,MXCOOR)
143      END IF
144
145      CALL DZERO(ECONS,NXETRAN)
146      CALL DZERO(OCONS,NAXGRO*NXGRST)
147
148      IOPTRES = 5
149      IORDER  = 1
150      CALL CC_XIETA(IXETRAN,NXETRAN,IOPTRES, IORDER, 'LE ',
151     &              CDUM, IDUM, RDUM, 'RE ', IEDOTS, ECONS,
152     &              1, WORK,LWORK)
153
154
155      IF (LOCDBG) WRITE(LUPRI,*) 'returned from cc_xieta...'
156      IF (LOCDBG) WRITE(LUPRI,*) 'call cc_dotdrv...'
157
158      CALL CC_DOTDRV('E0 ','O1 ',NXGRST,NAXGRO,IOTRAN,IODOTS,OCONS,
159     &               WORK,LWORK)
160
161      IF (LOCDBG) WRITE(LUPRI,*) 'returned from cc_dotdrv...'
162
163
164      ! calculate nuclear repulsion contribution to gradient
165      CALL NUCREP(WORK,WORK(MXCOOR*MXCOOR+1),WORK(2*MXCOOR*MXCOOR+1))
166
167      NXETRAN = 0
168      DO IOP = 1, NAXGRO
169        IOPER = IAXGRO(IOP)
170        IF (LPDBSOP(IOPER)) THEN
171         IF (ISYOPR(IOPER).EQ.1) THEN
172           IRELAX = IR1TAMP(LBLOPR(IOPER),.TRUE.,0.0D0,ISYOPR(IOPER))
173           IDXO1  = IRHSR1(LBLOPR(IOPER),.TRUE.,0.0D0,ISYOPR(IOPER))
174           PROPN  = CC_NUCCON(LBLOPR(IOPER),ISYOPR(IOPER))
175         ELSE
176           PROPN  = 0.0D0
177         END IF
178         DO IX = 1, NXGRST
179           NXETRAN = NXETRAN + 1
180
181           IF (ISYOPR(IOPER).EQ.1) THEN
182             PROPE = ECONS(NXETRAN)+OCONS(IOP,IX)+AVEO1(IDXO1)
183             PROPT = PROPN + PROPE
184           ELSE
185             PROPE = 0.0D0
186             PROPT = PROPN + PROPE
187           END IF
188
189           ISYST = ILSTSYM('LE ',IXGRST(IX))
190           IEX   = IDXSYM('LE ',ISYST,IXGRST(IX))
191
192           IF (XGROPT .AND. ISYST.EQ.IXSTSY .AND. IEX.EQ.IXSTAT .AND.
193     &         LBLOPR(IOPER)(1:5).EQ.'1DHAM') THEN
194             READ(LBLOPR(IOPER)(6:8),'(I3)') ISCOOR
195             GRADEE(ISCOOR) = PROPE
196           END IF
197
198           IF (ISYOPR(IOPER).EQ.1) THEN
199            WRITE(LUPRI,'(10x,A9,2x,F10.6,2x,F10.6,2x,F10.6)')
200     &        LBLOPR(IOPER)//':',PROPT,PROPE,PROPN
201            IF (LOCDBG) THEN
202             WRITE(LUPRI,'(20X,A9,1X,F10.6)') LBLOPR(IOPER)//':', PROPT
203             WRITE(LUPRI,'(20X,A9,1X,F10.6)')'E A^(1) E:',ECONS(NXETRAN)
204             WRITE(LUPRI,'(20X,A9,1X,F10.6)')'E0 O1    :',OCONS(IOP,IX)
205             WRITE(LUPRI,'(20X,A9,1X,F10.6)')'AVERAGE  :',AVEO1(IDXO1)
206             WRITE(LUPRI,'(20X,A9,1X,F10.6)')'NUCCON   :',PROPN
207            END IF
208           ELSE
209             WRITE(LUPRI,'(10X,A9,1X,A)') LBLOPR(IOPER)//':',
210     &         ' zero by symmetry '
211           END IF
212
213         END DO
214        END IF
215      END DO
216
217      WRITE(LUPRI,*) ' '
218
219      IF (XGROPT) THEN
220        NCOOR = 3*NUCDEP ! define NCOOR in trkoor.h, used in ABAREAD_TAYMOL
221        IF ( IPRINT.GT.1) THEN
222          CALL CC_SETDORPS('1DHAM   ',.TRUE.,IPRINT)
223          ! calculate nuclear repulsion contribution to gradient
224          CALL NUCREP(WORK,WORK(MXCOOR*MXCOOR+1),
225     &                     WORK(2*MXCOOR*MXCOOR+1))
226          CALL HEADER('Electronic gradient',-1)
227          CALL PRIGRD(GRADEE,WORK,WORK(MXCOOR*MXCOOR+1))
228        END IF
229
230        CALL ZERGRD
231        CALL ADDGRD(GRADNN)
232        CALL ADDGRD(GRADEE)
233        allocate ( GRDMOL(NCOOR), HESMOL(NCOOR, NCOOR) )
234        CALL ABAREAD_TAYMOL(ERGMOL,GRDMOL,HESMOL,NCOOR)
235
236        CALL HEADER('Molecular gradient',-1)
237        CALL PRIGRD(GRDMOL,WORK,WORK(MXCOOR*MXCOOR+1))
238
239        GRDNRM = DDOT(NCOOR,GRDMOL,1,GRDMOL,1)
240        WRITE (LUPRI,'(/19X,A,1P,E10.2)')
241     *     'Molecular gradient norm:', GRDNRM
242        deallocate ( GRDMOL, HESMOL )
243
244      END IF
245
246
247      CALL QEXIT('CCX_GRAD')
248
249      RETURN
250      END
251
252*---------------------------------------------------------------------*
253*              END OF SUBROUTINE CCX_GRAD                             *
254*---------------------------------------------------------------------*
255