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
19c /* deck cc_cphf */
20*=====================================================================*
21      SUBROUTINE CC_CPHF(TYPE,LABEL,ISYMS,ISTAT,EIGV,
22     &                   ISYMO,FREQS,ICAU,NVEC,MAXVEC,
23     &                   WORK,LWORK)
24*---------------------------------------------------------------------*
25*
26*    Purpose: solve CPHF equations needed in CC program
27*
28*    implemented types:  R1
29*
30*    Written by Christof Haettig, november 1998
31*
32*=====================================================================*
33#if defined (IMPLICIT_NONE)
34      IMPLICIT NONE
35#else
36#include "implicit.h"
37#endif
38#include "priunit.h"
39#include "dummy.h"
40#include "exeinf.h"
41#include "maxorb.h"
42#include "infvar.h"
43#include "inftap.h"
44#include "iratdef.h"
45#include "ccsdinp.h"
46#include "ccsdsym.h"
47#include "ccorb.h"
48#include "cclr.h"
49#include "ccfro.h"
50#include "inflin.h"
51Cholesky
52#include "ccdeco.h"
53Cholesky
54
55* local parameters:
56      CHARACTER*(18) MSGDBG
57      PARAMETER (MSGDBG = '[debug] CC_CPHF>  ')
58      LOGICAL LOCDBG
59      PARAMETER (LOCDBG = .true. )
60      CHARACTER*8 FNGDVE, FNSOVE, FNREVE, FILCPHF
61      PARAMETER(FNGDVE='CCPHFRHS',FNSOVE='CCPHFSOL')
62      PARAMETER(FNREVE='CCPHFRES')
63      INTEGER MXFRVAL
64      PARAMETER (MXFRVAL = 100)
65C     PARAMETER (MXFRVAL = 1)  ! switch off simultaneous equations
66
67      CHARACTER TYPE*(3)
68
69      INTEGER NVEC, MAXVEC, LWORK
70      INTEGER ISYMS(MAXVEC,*), ISYMO(MAXVEC,*)
71      INTEGER ISTAT(MAXVEC,*), ICAU(MAXVEC,*)
72
73      CHARACTER*8 LABEL(MAXVEC,*)
74
75      REAL*8  FREQS(MAXVEC,*), EIGV(MAXVEC,*)
76      REAL*8  WORK(LWORK)
77      REAL*8  ZERO, RDUM
78      REAL*8  XNORM, DDOT, DSQRT, GDNORM
79
80      PARAMETER (ZERO = 0.0d0)
81
82      CHARACTER MODEL*(10), MODWR*(10)
83
84      LOGICAL MAYBE_MORE, RELAX, CICLC, HFCLC, TRPCLC,OOTV
85      LOGICAL EXCLC, NEWCMO_SAVE, OPTST, EX
86      INTEGER IOPT, ISYM, IVEC, NSTAT, ORDER, IDUM, IPERT,LWRK0,LWRK1
87      INTEGER NCMOT,NASHT,N2ASHX,LCINDX,KCMO,KUDV,KXINDX,KEND0,KEND1
88      INTEGER LUGDVE, LUSOVE, MXRM, MXPHP, NCOSAV, IOPTWR, KSLVEC
89      INTEGER IREAL, NFRVAL, KFRVAL, NABAOP, NABATY, IFRVAL, IDX
90      INTEGER IPRABA, LUREVE, LUCPHF, MALLAI, INUM
91
92* external functions:
93      INTEGER IR1TAMP
94
95*---------------------------------------------------------------------*
96* do some checks:
97*---------------------------------------------------------------------*
98      CALL QENTER('CC_CPHF')
99
100      IF (NVEC.LT.1) THEN
101        CALL QEXIT('CC_CPHF')
102        RETURN
103      END IF
104
105      IF (LOCDBG) THEN
106         WRITE(LUPRI,*) 'Entered CC_CPHF. NVEC =',NVEC
107         WRITE(LUPRI,*) 'TYPE ',TYPE
108         CALL FLSHFO(LUPRI)
109      END IF
110
111* check vector type:
112      IF (TYPE(1:3).EQ.'R1 ') THEN
113        ORDER = 1
114        NSTAT = 0
115        MODWR = 'SCF?      '
116      ELSE
117        WRITE (LUPRI,*) 'CPHF vectors ',TYPE(1:3),' not implemented.'
118        CALL QUIT('required CPHF vectors not implemented.')
119      END IF
120
121* Refuse to do anything for Cholesky.
122      IF (CHOINT) THEN
123         NWARN = NWARN + 1
124         WRITE(LUPRI,'(/A/A/A/A/)') '*** WARNING from CC_CPHF:',
125     &      ' WARNING: Refusing to do CPHF for Cholesky, because',
126     &      ' WARNING: ABACUS has not been modified yet!',
127     &      ' WARNING: Program continues nevertheless...'
128         RETURN
129      ENDIF
130
131* get some variables from SIRIUS common blocks
132      CALL CC_SIRINF(NCMOT,NASHT,N2ASHX,LCINDX)
133
134* check number of active shells from SIRIUS:
135      IF (NASHT.NE.0) THEN
136        WRITE (LUPRI,*) 'non-zero number of active shells:',NASHT
137        CALL QUIT('Non-zero number of active shells in CC_CPHF.')
138      END IF
139
140*---------------------------------------------------------------------*
141* print header for rhs vector section
142*---------------------------------------------------------------------*
143      WRITE (LUPRI,'(7(/1X,2A),/)')
144     & '------------------------------------',
145     &                               '-------------------------------',
146     & '|                   OUTPUT FROM CPHF',
147     &                               ' SECTION                      |',
148     & '------------------------------------',
149     &                               '-------------------------------'
150
151* print some debug/info output
152      IF (IPRINT .GT. 10 .OR. LOCDBG) THEN
153        WRITE(LUPRI,*) 'CC_CPHF Workspace:',LWORK
154        WRITE(LUPRI,*) '1 MODWR ',MODWR
155      END IF
156
157      CALL FLSHFO(LUPRI)
158
159*---------------------------------------------------------------------*
160* some initilizations:
161*---------------------------------------------------------------------*
162
163* maximum of nallai over isym (used as fixed record lengths for LUCPHF:
164      MALLAI = NALLAI(1)
165      DO ISYM = 2, NSYM
166        MALLAI = MAX(MALLAI,NALLAI(ISYM))
167      END DO
168
169* arrays for GETGPV and ABARSP:
170      KCMO   = 1
171      KUDV   = KCMO   + NCMOT
172      KXINDX = KUDV   + N2ASHX
173      KFRVAL = KXINDX + LCINDX
174      KEND0  = KFRVAL + MXFRVAL
175      LWRK0  = LWORK  - KEND0
176
177      IF (LWRK0 .LT. 0) THEN
178         CALL QUIT('Insufficient memory in CC_CPHF.')
179      END IF
180
181* read MO coefficients from file:
182      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
183     &            .FALSE.)
184      REWIND LUSIFC
185      CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI)
186      READ (LUSIFC)
187      READ (LUSIFC)
188      CALL READI(LUSIFC,IRAT*NCMOT,WORK(KCMO))
189      CALL GPCLOSE(LUSIFC,'KEEP')
190
191* flags for ABARSP:
192      CICLC  = .FALSE.
193      HFCLC  = .TRUE.
194      TRPCLC = .FALSE.
195      OOTV   = .FALSE.
196      EXCLC  = .FALSE.
197      MXRM   = 40
198      MXPHP  = 1
199
200      NEWCMO_SAVE = NEWCMO
201      NCOSAV = NCONF
202
203      NEWCMO = .TRUE.
204      NCONF  = 1
205
206      IF (DIRECT) CALL CCDFFOP
207
208* flags for CC_WRRSP routine:
209      IOPTWR = 4
210
211* open file for right hand side and solution vectors:
212
213      LUGDVE = -1
214      LUSOVE = -1
215      LUREVE = -1
216      CALL GPOPEN(LUGDVE,FNGDVE,'UNKNOWN',' ','UNFORMATTED',IDUMMY,
217     &            .FALSE.)
218      CALL GPOPEN(LUSOVE,FNSOVE,'UNKNOWN',' ','UNFORMATTED',IDUMMY,
219     &            .FALSE.)
220      CALL GPOPEN(LUREVE,FNREVE,'UNKNOWN',' ','UNFORMATTED',IDUMMY,
221     &            .FALSE.)
222* open property integral file:
223      IF (LUPROP .LE. 0) CALL GPOPEN(LUPROP,'AOPROPER','UNKNOWN',' ',
224     &                               'UNFORMATTED',IDUMMY,.FALSE.)
225
226* close twoelectron integral file 'AOTWOINT':
227      IF (LUINTA.GT.0) CALL GPCLOSE(LUINTA,'KEEP')
228
229* file for CPHF vectors:
230      WRITE(FILCPHF,'(A5,A3)') 'CPHF_',TYPE(1:3)
231      DO I = 6,8
232        IF (FILCPHF(I:I).EQ.' ') FILCPHF(I:I) = '_'
233      END DO
234
235      LUCPHF = -1
236      CALL WOPEN2(LUCPHF,FILCPHF,64,0)
237
238*---------------------------------------------------------------------*
239* loop over vectors, set up rhs vectors and call ABARSP to solve eqs.:
240*---------------------------------------------------------------------*
241
242      IVEC = 0
243      DO WHILE(IVEC.LT.NVEC)
244        IVEC = IVEC + 1
245
246        ISYM = 1
247        DO IDX = 1, NSTAT
248          ISYM = MULD2H(ISYM,ISYMS(IVEC,IDX))
249        END DO
250        DO IDX = 1, ORDER
251          ISYM = MULD2H(ISYM,ISYMO(IVEC,IDX))
252        END DO
253
254        IF (LWRK0 .LT. NALLAI(ISYM)) THEN
255           CALL QUIT('Insufficient memory in CC_CPHF.')
256        END IF
257
258        IF (LOCDBG) THEN
259          WRITE (LUPRI,*) 'CC_CPHF> GP vector, label, symmetry:',
260     &                     LABEL(IVEC,1), ISYM
261          WRITE(LUPRI,*) '2 MODWR ',MODWR
262          CALL FLSHFO(LUPRI)
263        END IF
264
265C       ---------------------------
266C       get right hand side vector:
267C       ---------------------------
268        CALL CC_GETHFGD(IVEC,TYPE,LABEL,ISYMS,ISTAT,EIGV,ISYMO,
269     &                  FREQS,ICAU,NVEC,MAXVEC,IREAL,
270     &                  WORK(KCMO),WORK(KUDV),WORK(KXINDX),
271     &                  WORK(KFRVAL),WORK(KEND0),LWRK0)
272
273        IF (LOCDBG) THEN
274          WRITE (LUPRI,'(5X,I5,F12.8)') (I,WORK(KEND0-1+I),
275     &                     I=1,NALLAI(ISYM))
276          WRITE(LUPRI,*) '3 MODWR ',MODWR
277        END IF
278
279
280C       ------------------------------
281C       write gradient vector to file:
282C       ------------------------------
283        REWIND LUGDVE
284        CALL WRITT(LUGDVE,NALLAI(ISYM),WORK(KEND0))
285
286
287C       ----------------------------------
288C       check norm of the gradient vector:
289C       ----------------------------------
290        GDNORM=DSQRT(DDOT(NALLAI(ISYM),WORK(KEND0),1,WORK(KEND0),1))
291        IF (LOCDBG) WRITE (LUPRI,*) 'GDNORM:',GDNORM
292        IF (LOCDBG) WRITE (LUPRI,*) '4 MODWR ',MODWR
293
294
295C       --------------------------------------------------------
296C       for 'R1' check if several frequencies for same operator:
297C       --------------------------------------------------------
298        NFRVAL     = 1
299        MAYBE_MORE = .TRUE.
300        DO WHILE (MAYBE_MORE)
301           IF (TYPE(1:3).EQ.'R1 '
302     &         .AND. IVEC.LT.NVEC .AND. NFRVAL.LT.MXFRVAL
303     &         .AND. LABEL(IVEC,1).EQ.LABEL(IVEC+1,1) ) THEN
304             WORK(KFRVAL+NFRVAL) = FREQS(IVEC+1,1)
305             IVEC   = IVEC + 1
306             NFRVAL = NFRVAL + 1
307           ELSE
308             MAYBE_MORE = .FALSE.
309           END IF
310        END DO
311
312
313        IF (GDNORM.GT.THRLDPHF) THEN
314C
315          CALL SETSIR(.FALSE.,WORK(KEND0),LWRK0)
316
317          IF (LOCDBG) WRITE (LUPRI,*) 'going to call ABARSP'
318          IF (LOCDBG) WRITE (LUPRI,*) '5 MODWR ',MODWR
319C         -----------------------
320C         solve CPHF equation(s):
321C         -----------------------
322          IPRABA = IPRINT
323          NABAOP = 1
324          NABATY = IREAL  ! flag for real/imaginary operators
325          CALL ABARSP(CICLC,HFCLC,TRPCLC,OOTV,ISYM,EXCLC,
326     &                WORK(KFRVAL),NFRVAL,NABATY,NABAOP,
327     &                LABEL(IVEC,1),LUGDVE,LUSOVE,LUREVE,THRLEQ,
328     &                MAXITE,IPRABA,MXRM,MXPHP,WORK(KEND0),LWRK0)
329
330          IF (LOCDBG) WRITE (LUPRI,*) 'returned from ABARSP'
331          IF (LOCDBG) WRITE (LUPRI,*) '6 MODWR ',MODWR
332
333C         --------------------------------------------
334C         clean up `left overs' from ABARSP:
335C         --------------------------------------------
336          CALL GPCLOSE(LUSIFC,'KEEP')
337          IF (LUONEL .GT. 0) CALL GPCLOSE(LUONEL,'KEEP')
338          IF (LUPROP .GT. 0) CALL GPCLOSE(LUPROP,'KEEP')
339C
340C         --------------------------------------------
341C         read solution vector(s) and save on CC file:
342C         --------------------------------------------
343          REWIND LUSOVE
344
345        END IF
346
347        DO IFRVAL = 1, NFRVAL
348           KSLVEC = KEND0
349           KEND1  = KSLVEC + 2*MALLAI
350           LWRK1  = LWORK  - KEND1
351           IF (LWRK1.LT.0) THEN
352              CALL QUIT('Insufficient memory in CC_CPHF.')
353           END IF
354
355           CALL DZERO(WORK(KSLVEC),2*MALLAI)
356           IF (GDNORM.GT.THRLDPHF) THEN
357             CALL  READT(LUSOVE,2*NALLAI(ISYM),WORK(KSLVEC))
358           END IF
359
360
361           ! save on CPHF vector file
362           IDX = IVEC - NFRVAL + IFRVAL
363           CALL PUTWA2(LUCPHF,FILCPHF,WORK(KSLVEC),
364     &                 1+2*MALLAI*(IDX-1),2*MALLAI)
365
366
367           ! check if a corresponding CC vector exists
368           INUM = -1
369           IF (TYPE(1:3).EQ.'R1 ') THEN
370            INUM=IR1TAMP(LABEL(IDX,1),.TRUE.,WORK(KFRVAL-1+IFRVAL),ISYM)
371           END IF
372           IF (LOCDBG) WRITE (LUPRI,*) 'IFRVAL, TYPE, MODWR ',
373     &         IFRVAL,TYPE,MODWR
374
375           ! if yes put the CPHF also on the CC vector file
376           IF (INUM.GT.0) THEN
377             CALL CC_WRRSP(TYPE,INUM,ISYM,IOPTWR,MODWR,
378     &                   WORK(KSLVEC),DUMMY,DUMMY,WORK(KEND1),LWRK1)
379           END IF
380
381           IF (LOCDBG) THEN
382              WRITE (LUPRI,*)
383     &           'CC_CPHF> solution vector, label, freq:',
384     &           LABEL(IDX,1),WORK(KFRVAL-1+IFRVAL)
385              WRITE (LUPRI,'(5X,I5,F12.8)')
386     &           (I,WORK(KSLVEC-1+I),I=1,2*NALLAI(ISYM))
387              WRITE (LUPRI,*)
388     &           'CC_CPHF> saved CPHF solution for ',TYPE,
389     &           ' equation nb. ',IDX,INUM
390           END IF
391           IF (LOCDBG) WRITE (LUPRI,*) '7 MODWR ',MODWR
392
393        END DO
394
395      END DO
396*---------------------------------------------------------------------*
397* that's it: close files, restore variables and return
398*---------------------------------------------------------------------*
399      CALL WCLOSE2(LUCPHF,FILCPHF,'KEEP')
400
401      IF (LUINTM .GT. 0) CALL GPCLOSE(LUINTM,'DELETE')
402      CALL GPCLOSE(LUGDVE,'DELETE')
403      CALL GPCLOSE(LUSOVE,'DELETE')
404      CALL GPCLOSE(LUREVE,'DELETE')
405
406      IF (LUINTA .LE. 0) THEN
407        CALL MAKE_AOTWOINT(WORK,LWORK)
408        CALL GPOPEN(LUINTA,'AOTWOINT','UNKNOWN',' ','UNFORMATTED',
409     *            IDUMMY,.FALSE.)
410      END IF
411
412      NEWCMO = NEWCMO_SAVE
413      NCONF  = NCOSAV
414
415      IF (LOCDBG) THEN
416         WRITE(LUPRI,*) 'leave CC_CPHF'
417         CALL FLSHFO(LUPRI)
418      END IF
419           IF (LOCDBG) WRITE (LUPRI,*) '8 MODWR ',MODWR
420
421      CALL QEXIT('CC_CPHF')
422      RETURN
423      END
424*=====================================================================*
425*              END OF SUBROUTINE CC_CPHF                              *
426*=====================================================================*
427c /* deck cc_sirinf */
428*=====================================================================*
429      SUBROUTINE CC_SIRINF(NCMOT1,NASHT1,N2ASHX1,LCINDX1)
430*---------------------------------------------------------------------*
431*   Purpose: read some variables from SIRIUS common blocks
432*---------------------------------------------------------------------*
433#include "implicit.h"
434#include "inforb.h"
435#include "inftap.h"
436#include "infdim.h"
437
438      NCMOT1  = NCMOT
439      NASHT1  = NASHT
440      N2ASHX1 = N2ASHX
441      LCINDX1 = LCINDX
442      RETURN
443      END
444*=====================================================================*
445c /* deck cc_rdhfrsp */
446*=====================================================================*
447      SUBROUTINE CC_RDHFRSP(LIST,IDXLST,ISYM,XKAPPA)
448*---------------------------------------------------------------------*
449C
450C   Purpose:  Read a CPHF response vector from file
451C             for explanation of LIST, IDXLIST & MODFIL see CC_RDRSP
452C
453C  Christof Haettig, summer 2003
454*---------------------------------------------------------------------*
455      IMPLICIT NONE
456#include "priunit.h"
457#include "dummy.h"
458#include "ccorb.h"
459#include "ccsdsym.h"
460#include "ccfro.h"
461
462#if defined (SYS_CRAY)
463      REAL XKAPPA(*)
464#else
465      DOUBLE PRECISION XKAPPA(*)
466#endif
467
468      CHARACTER FILCPHF*8, LIST*3
469      INTEGER MALLAI, ISYM, JSYM, LUCPHF, IDXLST
470
471      CALL QENTER('CCRDHFRSP')
472
473* set file name:
474      WRITE(FILCPHF,'(A5,A3)') 'CPHF_',LIST(1:3)
475      DO I = 6,8
476        IF (FILCPHF(I:I).EQ.' ') FILCPHF(I:I) = '_'
477      END DO
478
479* calculate record lengths:
480      MALLAI = NALLAI(1)
481      DO JSYM = 2, NSYM
482        MALLAI = MAX(MALLAI,NALLAI(JSYM))
483      END DO
484
485* open direct access file:
486      LUCPHF = -1
487      CALL WOPEN2(LUCPHF,FILCPHF,64,0)
488
489* read vector number IDXLST from file:
490      CALL GETWA2(LUCPHF,FILCPHF,XKAPPA,
491     &                 1+2*MALLAI*(IDXLST-1),2*NALLAI(ISYM))
492
493* close file and return:
494      CALL WCLOSE2(LUCPHF,FILCPHF,'KEEP')
495
496      CALL QEXIT('CCRDHFRSP')
497      RETURN
498      END
499*=====================================================================*
500*              END OF SUBROUTINE CC_RDHFRSP                           *
501*=====================================================================*
502