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_PQI(CTR2,ISYCTR,T2AM,ISYTAMP,
21     &                  FILNAM,LUPQIM,IADRPQ,IADR,IV,WORK,LWORK)
22*----------------------------------------------------------------------
23*
24*     Purpose: Calculate the P and Q intermediates from the
25*              Lagrangian multipiers CTR2
26*              and the amplitude vector T2AM
27*              and write them to a file direct access file FILNAM
28*
29*   P{aik,del} = sum_{dl} (2 t(dl,k;del) - t(dk,l;del)) Zeta(dl;ai)
30*   P{aik,del} = sum_{dl} t(dl,k;del) (2Zeta(di;al) + Zeta(dl;ai))
31*
32*   Christof Haettig & Asger Halkier, August 1998
33*
34*======================================================================
35#if defined (IMPLICIT_NONE)
36      IMPLICIT NONE
37#else
38#  include "implicit.h"
39#endif
40#include "priunit.h"
41#include "ccorb.h"
42#include "maxorb.h"
43#include "ccsdsym.h"
44#include "ccsdinp.h"
45#include "iratdef.h"
46#include "second.h"
47
48      LOGICAL LOCDBG
49      PARAMETER (LOCDBG = .FALSE.)
50      INTEGER LWORK, LUPQIM
51      CHARACTER*(*) FILNAM
52
53#if defined (SYS_CRAY)
54      REAL ONE, TWO, THREE, HALF, XNORM, DDOT, ZERO
55      REAL WORK(*), CTR2(*), T2AM(*)
56      REAL TIMET, TIMIO, TIMZWV, TIMTAM, DTIME, TIMTI
57      REAL TIMZET,TIMSCL
58#else
59      DOUBLE PRECISION ONE, TWO, THREE, HALF, DDOT, ZERO
60      DOUBLE PRECISION WORK(*), CTR2(*), T2AM(*)
61      DOUBLE PRECISION TIMET,TIMIO,TIMTI,TIMZWV,TIMTAM,DTIME
62      DOUBLE PRECISION TIMZET,TIMSCL
63#endif
64      PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,HALF=0.5D0,ZERO=0.0D0)
65
66      INTEGER IV, IADR, IADRQ
67      INTEGER IADRPQ(MXCORB_CC,IV)
68      INTEGER ISYTAMP, ISYCTR, ISYMD, ISYTIN, ISYAIK
69      INTEGER KEND1, LWRK1, KEND2, LWRK2
70      INTEGER IOPT, IDEL, ILLL
71      INTEGER KTINT, KTJNT, KPINT, KQINT, KLAMDH, KLAMDP, KT1AM
72
73
74      CALL QENTER('CC_PQI')
75*----------------------------------------------------------------------
76* set symmetries and allocate work space:
77*----------------------------------------------------------------------
78      KLAMDH = 1
79      KLAMDP = KLAMDH + NLAMDT
80      KT1AM  = KLAMDP + NLAMDT
81      KEND1  = KT1AM  + NT1AMX
82      LWRK1  = LWORK - KEND1
83
84      IF ( LWRK1 .LT. 0 ) THEN
85         WRITE (LUPRI,*) 'Insufficient memory in CC_PQI.'
86         CALL QUIT('Insufficient memory in CC_PQI.')
87      END IF
88
89      TIMET  = SECOND()
90      TIMIO  = ZERO
91      TIMTI  = ZERO
92      TIMZWV = ZERO
93      TIMTAM = ZERO
94      TIMZET = ZERO
95      TIMSCL = ZERO
96
97      IF (LOCDBG) THEN
98         WRITE (LUPRI,*) 'Norm of T2AMP:',
99     &        DDOT(NT2AM(ISYTAMP),T2AM,1,T2AM,1)
100         WRITE (LUPRI,*) 'T2AMP:'
101         CALL CC_PRP(WORK,T2AM,ISYTAMP,0,1)
102         WRITE (LUPRI,*) 'Norm of CTR2 :',
103     &        DDOT(NT2SQ(ISYCTR),CTR2,1,CTR2,1)
104         WRITE (LUPRI,*) 'CTR2:'
105         CALL CC_PRSQ(WORK,CTR2,ISYCTR,0,1)
106      END IF
107*----------------------------------------------------------------------
108* get XLAMD matrices from zero T1 amplitudes:
109*----------------------------------------------------------------------
110      CALL DZERO(WORK(KT1AM),NT1AMX)
111
112      CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),
113     &            WORK(KEND1),LWRK1)
114
115*----------------------------------------------------------------------
116* calculate the P intermediate in a loop over AO index delta:
117*----------------------------------------------------------------------
118      DO ISYMD = 1, NSYM
119      DO ILLL = 1, NBAS(ISYMD)
120        IDEL = IBAS(ISYMD) + ILLL
121
122        ISYTIN = MULD2H(ISYMD,ISYTAMP)
123        ISYAIK = MULD2H(ISYTIN,ISYCTR)
124
125        KTINT = KEND1
126        KTJNT = KTINT + NT2BCD(ISYTIN)
127        KPINT = KTJNT + NT2BCD(ISYTIN)
128        KQINT = KPINT + NT2BCD(ISYAIK)
129        KEND2 = KQINT + NT2BCD(ISYAIK)
130        LWRK2 = LWORK  - KEND2
131
132        IF  (LWRK2 .LT. 0) THEN
133           WRITE (LUPRI,*) 'Insufficient memory in CC_PQIM.'
134           CALL QUIT('Insufficient memory in CC_PQIM.')
135        END IF
136C
137C       calculate delta batch of backtransformed amplitudes:
138C
139        DTIME = SECOND()
140        CALL CC_TI(WORK(KTINT),ISYTIN,T2AM,ISYTAMP,WORK(KLAMDH),1,
141     &             WORK(KEND2),LWRK2,IDEL,ISYMD)
142        TIMTI = TIMTI + SECOND() - DTIME
143C
144C       calculate 2 x Coul - Exc combination of backtransformed T:
145C
146        DTIME = SECOND()
147        CALL DCOPY(NT2BCD(ISYTIN),WORK(KTINT),1,WORK(KTJNT),1)
148
149        CALL DSCAL(NT2BCD(ISYTIN),-ONE,WORK(KTJNT),1)
150
151        CALL CCLT_P21I(WORK(KTJNT),ISYTIN,WORK(KEND2),LWRK2,
152     &                 IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR)
153
154        CALL DAXPY(NT2BCD(ISYTIN),TWO,WORK(KTINT),1,WORK(KTJNT),1)
155        TIMTAM = TIMTAM + SECOND() - DTIME
156C
157C       calculate the P intermediate:
158C
159        IOPT = 1
160        DTIME = SECOND()
161        CALL CC_ZWVI(WORK(KPINT),CTR2,ISYCTR,WORK(KTJNT),
162     &               ISYTIN,WORK(KEND2),LWRK2,IOPT)
163        TIMZWV = TIMZWV + SECOND() - DTIME
164
165        DTIME = SECOND()
166        CALL DSCAL(NT2BCD(ISYAIK),HALF,WORK(KPINT),1)
167        TIMSCL = TIMSCL + SECOND() - DTIME
168C
169C       write the intermediate to file:
170C
171        IADRPQ(IDEL,IV) = IADR
172
173        DTIME = SECOND()
174        CALL  PUTWA2(LUPQIM,FILNAM,WORK(KPINT),IADR,NT2BCD(ISYAIK))
175        TIMIO = TIMIO + SECOND() - DTIME
176
177        IADR = IADR + 2*NT2BCD(ISYAIK)
178
179        IF (LOCDBG) THEN
180           WRITE (LUPRI,*) 'CC_PQI> P interm. in AO for IDEL = ',IDEL
181           WRITE (LUPRI,'(5(F12.8))')
182     &          (WORK(KPINT+I),I=0,NT2BCD(ISYAIK)-1)
183        END IF
184
185      END DO
186      END DO
187
188*----------------------------------------------------------------------
189*     calculate (2 x Exc + Coul)/3 combination of ZETA:
190*     (we interrupt here the loop over AO to calculate the modified
191*      ZETA and accept that we recalculate the backtransformed
192*      amplitude since it the transformation of the amplitudes is
193*      less expansive than the transformation and restruction of ZETA
194*      for each delta. both the transformation of TAMP and the
195*      transformation/restruction of ZETA scale with N V^2 O^2. )
196*----------------------------------------------------------------------
197      DTIME = SECOND()
198      CALL CCRHS_T2BT(CTR2,WORK(KEND1),LWRK1,ISYCTR)
199      CALL CCSD_T2TP(CTR2,WORK(KEND1),LWRK1,ISYCTR)
200      TIMZET = TIMZET + SECOND() - DTIME
201
202*----------------------------------------------------------------------
203* calculate the Q intermediate in a loop over AO index delta:
204*----------------------------------------------------------------------
205      DO ISYMD = 1, NSYM
206      DO ILLL = 1, NBAS(ISYMD)
207        IDEL = IBAS(ISYMD) + ILLL
208
209        ISYTIN = MULD2H(ISYMD,ISYTAMP)
210        ISYAIK = MULD2H(ISYTIN,ISYCTR)
211
212        KTINT = KEND1
213        KTJNT = KTINT + NT2BCD(ISYTIN)
214        KPINT = KTJNT + NT2BCD(ISYTIN)
215        KQINT = KPINT + NT2BCD(ISYAIK)
216        KEND2 = KQINT + NT2BCD(ISYAIK)
217        LWRK2 = LWORK  - KEND2
218
219        IF (LWRK2 .LT. 0) THEN
220           WRITE (LUPRI,*) 'Insufficient memory in CC_PQIM.'
221           CALL QUIT('Insufficient memory in CC_PQIM.')
222        END IF
223C
224C       calculate delta batch of backtransformed amplitudes:
225C
226        DTIME = SECOND()
227        CALL CC_TI(WORK(KTINT),ISYTIN,T2AM,ISYTAMP,WORK(KLAMDH),1,
228     &             WORK(KEND2),LWRK2,IDEL,ISYMD)
229        TIMTI = TIMTI + SECOND() - DTIME
230C
231C       calculate Q intermediate:
232C
233        IOPT = 2
234        DTIME = SECOND()
235        CALL CC_ZWVI(WORK(KQINT),CTR2,ISYCTR,WORK(KTINT),
236     &               ISYTIN,WORK(KEND2),LWRK2,IOPT)
237        TIMZWV = TIMZWV + SECOND() - DTIME
238
239        DTIME = SECOND()
240        CALL DSCAL(NT2BCD(ISYAIK),THREE*HALF,WORK(KQINT),1)
241        TIMSCL = TIMSCL + SECOND() - DTIME
242C
243C       write the intermediate to file:
244C
245        IADRQ = IADRPQ(IDEL,IV) + NT2BCD(ISYAIK)
246
247        DTIME = SECOND()
248        CALL  PUTWA2(LUPQIM,FILNAM,WORK(KQINT),IADRQ,NT2BCD(ISYAIK))
249        TIMIO = TIMIO + SECOND() - DTIME
250
251        IF (LOCDBG) THEN
252           WRITE (LUPRI,*) 'CC_PQI> Q interm. in AO for IDEL = ',IDEL
253           WRITE (LUPRI,'(5(F12.8))')
254     &          (WORK(KQINT+I),I=0,NT2BCD(ISYAIK)-1)
255        END IF
256
257      END DO
258      END DO
259
260*----------------------------------------------------------------------
261* recover Zeta vector:
262*----------------------------------------------------------------------
263      DTIME = SECOND()
264      CALL CCSD_T2TP(CTR2,WORK(KEND1),LWRK1,ISYCTR)
265      CALL CCRHS_T2TR(CTR2,WORK(KEND1),LWRK1,ISYCTR)
266      TIMZET = TIMZET + SECOND() - DTIME
267
268*----------------------------------------------------------------------
269* return:
270*----------------------------------------------------------------------
271      IF (IPRINT.GE.10) THEN
272         TIMET = SECOND() - TIMET
273         WRITE(LUPRI,*) '  Timings of CC_PQIM:'
274         WRITE(LUPRI,1) 'I/O             ', TIMIO
275         WRITE(LUPRI,1) 'CC_TI           ', TIMTI
276         WRITE(LUPRI,1) '2C-E of T2      ', TIMTAM
277         WRITE(LUPRI,1) '2E+C of L2      ', TIMZET
278         WRITE(LUPRI,1) 'CC_ZWV          ', TIMZWV
279         WRITE(LUPRI,1) 'scaling etc.    ', TIMSCL
280         WRITE(LUPRI,1) 'total CC_PQIM:  ', TIMET
281      END IF
282
283   1  FORMAT(1x,'Time used for',2x,A18,2x,': ',f10.2,' seconds')
284
285      CALL QEXIT('CC_PQI')
286      RETURN
287      END
288C
289*======================================================================
290      SUBROUTINE CC_PQI3(CTR2,ISYCTR,T2AM,ISYTAMP,
291     &                    FILNAM,LUPQIM,IADRPQ,IADR,IV,
292     &                    XLAMDH,ICAL,WORK,LWORK)
293*----------------------------------------------------------------------
294*
295*     Purpose: Calculate the P intermediates from the
296*              Lagrangian multipiers CTR2
297*              and the amplitude vector T2AM
298*              and write them to a file direct access file FILNAM
299*
300*   P{aik,del} = sum_{dl} (2 t(dl,k;del) - t(dk,l;del)) Ctr2(dl;ai)
301C   where in the triplet case, ctr2 should contain
302C   (+)L(dl;ai) - (-)L(dl;ai)
303*
304C   Based on CCPQI31 by
305*   Christof Haettig & Asger Halkier, August 1998
306*
307*======================================================================
308      IMPLICIT NONE
309#include "priunit.h"
310#include "ccorb.h"
311#include "maxorb.h"
312#include "ccsdsym.h"
313#include "ccsdinp.h"
314#include "iratdef.h"
315#include "second.h"
316
317      CHARACTER, PARAMETER :: MYNAME = 'CC_PQI3'
318      INTEGER, PARAMETER :: NCAL = 3 ! Number of times this is called
319      LOGICAL LOCDBG
320      PARAMETER (LOCDBG = .FALSE.)
321
322
323      INTEGER, INTENT(IN) :: LWORK, LUPQIM
324      INTEGER, INTENT(INOUT) :: IADRPQ(MXCORB_CC,IV), IADR
325      INTEGER, INTENT(IN) :: ISYCTR, ISYTAMP, IV, ICAL
326      CHARACTER*(*),INTENT(IN) :: FILNAM
327
328
329      DOUBLE PRECISION, INTENT(IN) :: WORK(*), CTR2(*), T2AM(*),
330     &                                XLAMDH(*)
331      DOUBLE PRECISION TIMET,TIMIO,TIMTI,TIMZWV,TIMTAM,DTIME
332      DOUBLE PRECISION TIMZET,TIMSCL
333
334      DOUBLE PRECISION ONE, TWO, THREE, HALF, DDOT, ZERO
335      PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,HALF=0.5D0,ZERO=0.0D0)
336
337      INTEGER :: IADRQ
338      INTEGER ISYMD, ISYTIN, ISYAIK
339      INTEGER KEND1, LWRK1, KEND2, LWRK2
340      INTEGER IOPT, IDEL, ILLL
341      INTEGER KTINT, KTJNT, KPINT, KQINT, KLAMDH, KLAMDP, KT1AM
342
343
344      CALL QENTER(myname)
345*----------------------------------------------------------------------
346* set symmetries and allocate work space:
347*----------------------------------------------------------------------
348      KEND1  = 1
349      LWRK1  = LWORK
350
351      IF ( LWRK1 .LT. 0 ) THEN
352         WRITE (LUPRI,*) 'Insufficient memory in '//myname//'.'
353         CALL QUIT('Insufficient memory in '//myname//'.')
354      END IF
355
356      TIMET  = SECOND()
357      TIMIO  = ZERO
358      TIMTI  = ZERO
359      TIMZWV = ZERO
360      TIMTAM = ZERO
361      TIMZET = ZERO
362      TIMSCL = ZERO
363
364      IF (LOCDBG) THEN
365         WRITE (LUPRI,*) 'Norm of T2AMP:',
366     &        DDOT(NT2AM(ISYTAMP),T2AM,1,T2AM,1)
367         WRITE (LUPRI,*) 'T2AMP:'
368         CALL CC_PRP(WORK,T2AM,ISYTAMP,0,1)
369         WRITE (LUPRI,*) 'Norm of CTR2 :',
370     &        DDOT(NT2SQ(ISYCTR),CTR2,1,CTR2,1)
371         WRITE (LUPRI,*) 'CTR2:'
372         CALL CC_PRSQ(WORK,CTR2,ISYCTR,0,1)
373      END IF
374
375*----------------------------------------------------------------------
376* calculate the P intermediate in a loop over AO index delta:
377*----------------------------------------------------------------------
378      DO ISYMD = 1, NSYM
379      DO ILLL = 1, NBAS(ISYMD)
380        IDEL = IBAS(ISYMD) + ILLL
381
382        ISYTIN = MULD2H(ISYMD,ISYTAMP)
383        ISYAIK = MULD2H(ISYTIN,ISYCTR)
384
385        KTINT = KEND1
386        KTJNT = KTINT + NT2BCD(ISYTIN)
387        KPINT = KTJNT + NT2BCD(ISYTIN)
388        KEND2 = KPINT + NT2BCD(ISYAIK)
389        LWRK2 = LWORK  - KEND2
390
391        IF  (LWRK2 .LT. 0) THEN
392           WRITE (LUPRI,*) 'Insufficient memory in '//myname//'.'
393           CALL QUIT('Insufficient memory in '//myname//'.')
394        END IF
395C
396C       calculate delta batch of backtransformed amplitudes:
397C
398        DTIME = SECOND()
399        CALL CC_TI(WORK(KTINT),ISYTIN,T2AM,ISYTAMP,XLAMDH,1,
400     &             WORK(KEND2),LWRK2,IDEL,ISYMD)
401        TIMTI = TIMTI + SECOND() - DTIME
402C
403C       calculate 2 x Coul - Exc combination of backtransformed T:
404C
405        IF (ICAL .EQ. 0) THEN
406            DTIME = SECOND()
407            CALL DCOPY(NT2BCD(ISYTIN),WORK(KTINT),1,WORK(KTJNT),1)
408
409            CALL DSCAL(NT2BCD(ISYTIN),TWO,WORK(KTINT),1)
410
411            CALL CCLT_P21I(WORK(KTJNT),ISYTIN,WORK(KEND2),LWRK2,
412     &                 IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR)
413
414            CALL DAXPY(NT2BCD(ISYTIN),-ONE,WORK(KTJNT),1,WORK(KTINT),1)
415
416            IADRPQ(IDEL,IV) = IADR
417            IADR = IADR + NCAL*NT2BCD(ISYAIK)
418        ELSE
419            CALL CCLT_P21I(WORK(KTINT),ISYTIN,WORK(KEND2),LWRK2,
420     &                 IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR)
421        END IF
422        IADRQ = IADRPQ(IDEL,IV) + NT2BCD(ISYAIK)*ICAL
423
424        TIMTAM = TIMTAM + SECOND() - DTIME
425C
426C       calculate the P intermediate:
427C
428        IOPT = 1
429        DTIME = SECOND()
430        CALL CC_ZWVI3(WORK(KPINT),CTR2,ISYCTR,WORK(KTINT),
431     &               ISYTIN,WORK(KEND2),LWRK2)
432        TIMZWV = TIMZWV + SECOND() - DTIME
433
434        DTIME = SECOND()
435        CALL DSCAL(NT2BCD(ISYAIK),HALF,WORK(KPINT),1)
436        TIMSCL = TIMSCL + SECOND() - DTIME
437C
438C       write the intermediate to file:
439C
440
441        DTIME = SECOND()
442        CALL  PUTWA2(LUPQIM,FILNAM,WORK(KPINT),IADRQ,NT2BCD(ISYAIK))
443        TIMIO = TIMIO + SECOND() - DTIME
444
445
446        IF (LOCDBG) THEN
447           WRITE (LUPRI,*) 'CC_PQI> P interm. in AO for IDEL = ',IDEL
448           WRITE (LUPRI,'(5(F12.8))')
449     &          (WORK(KPINT+I),I=0,NT2BCD(ISYAIK)-1)
450        END IF
451
452      END DO
453      END DO
454
455      CALL QEXIT(myname)
456      RETURN
457      END
458
459