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!
18      SUBROUTINE D0FCDIAG(D0DIAG,FCDIAG,FC)
19#include "implicit.h"
20#include "priunit.h"
21#include "infrsp.h"
22#include "inforb.h"
23
24      DIMENSION D0DIAG(*),FCDIAG(*),FC(*)
25
26      DO 300 ISYM=1, NSYM
27        IOFF=IIORB(ISYM)+1
28
29        DO 100 I=1, NISH(ISYM)
30           D0DIAG(IORB(ISYM)+I) =2.0D0
31 100    CONTINUE
32
33
34        DO 200 I=1,NORB(ISYM)
35          FCDIAG(IORB(ISYM)+I)=FC(IOFF)
36          IOFF = IOFF + I + 1
37 200   CONTINUE
38 300  CONTINUE
39
40      IF (IPRRSP.GE.10) THEN
41        CALL PRINTVEC('FCDIAG:        ',1,NORBT,FCDIAG)
42        CALL PRINTVEC('D0DIAG:        ',1,NORBT,D0DIAG)
43
44      END IF
45      RETURN
46      END
47C
48C END OF D0FCDIAG
49C
50
51      SUBROUTINE GETH1MO(H1MO,CMO,WRK,LWRK)
52#include "implicit.h"
53#include "priunit.h"
54#include "infrsp.h"
55#include "inforb.h"
56
57      DIMENSION H1MO(NORBT,NORBT),CMO(*),WRK(*)
58
59C
60C     Calls sirh1 to get the one electron integrals
61C     in AO basis (h1AO) and transforms the triangular
62C     output to the full symmetrical matrix NORBT*NORBT
63C
64
65      KH1AO  = 1
66      KH1AOT = KH1AO + N2BASX
67
68      KWRK1 = KH1AOT + NNORBX+1
69      LWRK1 = LWRK - KWRK1
70
71      CALL SIRH1(WRK(KH1AOT),WRK(KWRK1),LWRK1)
72
73      K    = KH1AOT
74      KOFF = KH1AO
75
76      DO 300 ISYM=1,NSYM
77        DO 200 I=1,NBAS(ISYM)
78          DO 100 J=1,I
79             WRK(KOFF+(I-1)*NBAST+J-1)=WRK(K)
80             WRK(KOFF+(J-1)*NBAST+I-1)=WRK(K)
81             K=K+1
82 100      CONTINUE
83 200    CONTINUE
84        KOFF = KOFF + NBAS(ISYM)*NBAST + NBAS(ISYM)
85 300  CONTINUE
86
87      IF (IPRRSP .GE. 10) THEN
88        CALL PRINTMAT('H1AO            ',1,NBAST,WRK(KH1AO))
89      END IF
90
91C     Transform the integrals to MO basis
92
93      CALL AO2MO(WRK(KH1AO),H1MO,CMO,WRK(KWRK1),LWRK1)
94
95      RETURN
96      END
97
98      SUBROUTINE GETESG_DENMAT_FOCMAT(D,F,NBAST,NNBAST)
99#include "implicit.h"
100#include "rspprp.h"
101#include "esg.h"
102#include "dummy.h"
103#include "priunit.h"
104
105      DIMENSION D(NNBAST), F(NNBAST)
106
107      WRITE (LUPRI,'(A,/)')
108     &        'Reading density and fock matrices in 1-el. part'
109
110      CALL FLSHFO(LUPRI)
111
112      READ (LUESG) D
113      READ (LUESG) F
114
115      CALL GPCLOSE(LUESG,'KEEP')
116
117      RETURN
118      END
119
120      SUBROUTINE GETESG_DENMAT(DMAT,NDMAT,WRK,LWRK)
121#include "implicit.h"
122#include "mxcent.h"
123#include "rspprp.h"
124#include "esg.h"
125#include "dummy.h"
126#include "inforb.h"
127#include "abainf.h"
128#include "priunit.h"
129
130      DIMENSION DMAT(NBAST,NBAST,NDMAT),WRK(*)
131      PARAMETER (MXDMAT=50)
132      DIMENSION ISYMDM(MXDMAT)
133
134      WRITE(LUPRI,'(A,/)') 'Reading density matrices in 2-el. part'
135
136      KDSO = 1
137      KWRK = KDSO + N2BASX
138
139      DO 100 I=1,NDMAT
140        ISYMDM(I) = 0
141 100  CONTINUE
142      ISYMDM(3) = ISYME - 1
143      ISYMDM(4) = ISYME - 1
144
145      IF (KWRK.GT.LWRK) CALL STOPIT('GETESG_DENMAT',' ',KWRK,LWRK)
146
147      CALL GPOPEN(LUESG2,'ESG_DMAT',' ',' ',' ',IDUMMY,.FALSE.)
148      REWIND (LUESG2)
149
150      DO IDMAT=2,NDMAT
151         CALL READT(LUESG2,N2BASX,WRK)
152         CALL DSOTAO(WRK,DMAT(1,1,IDMAT),NBAST,ISYMDM(IDMAT),IPRESG)
153      END DO
154
155      CALL GPCLOSE(LUESG2,'KEEP')
156
157      IF ( IPRESG .GE. 10 ) THEN
158        WRITE(LUPRI,*) 'Reading density matrices in 2-el. part'
159        CALL PRINTMAT('D0AO:          ',1,NBAST,DMAT(1,1,1))
160        CALL PRINTMAT('D2AO:          ',1,NBAST,DMAT(1,1,2))
161        CALL PRINTMAT('DXAO:          ',1,NBAST,DMAT(1,1,3))
162        CALL PRINTMAT('DXSAO:         ',1,NBAST,DMAT(1,1,4))
163      END IF
164
165      RETURN
166      END
167
168      SUBROUTINE INVE2VEC(CMO,UDV,PV,FC,FV,FCAC,
169     & H2AC,NSIM,INPVECS,OUTVECS,XINDX,WRK,LWRK)
170C
171C CALCULATES INV( E2 ) * VECTOR
172C
173#include "implicit.h"
174#include "dummy.h"
175#include "codata.h"
176#include "priunit.h"
177#include "infopt.h"
178#include "infrsp.h"
179#include "wrkrsp.h"
180#include "rspprp.h"
181#include "infpp.h"
182#include "inflr.h"
183#include "inforb.h"
184#include "infdim.h"
185#include "infpri.h"
186#include "inftap.h"
187#include "infsop.h"
188
189      DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*)
190      DIMENSION XINDX(*),WRK(*)
191      DOUBLE PRECISION INPVECS(*),OUTVECS(*)
192      CHARACTER*8 BLANK
193C
194      PARAMETER ( MAXSIM = 15, BLANK = '        ' )
195
196      PARAMETER ( D2 = 2.0D0, D100 = 100.0D0, D0 = 0.0D0)
197      PARAMETER ( D2R3 = (2.0D0/3.0D0))
198C
199C     space allocation for reduced E(2) and reduced S(2)
200C
201      KREDE  = 1
202      KREDS  = KREDE  + MAXRM*MAXRM
203      KIBTYP = KREDS  + MAXRM*MAXRM
204      KEIVAL = KIBTYP + MAXRM
205      KRESID = KEIVAL + MAXRM
206      KEIVEC = KRESID + MAXRM
207      KREDGD = KEIVEC + MAXRM*MAXRM
208      KWRK1  = KREDGD + MAXRM*MAXRM
209      LWRK1  = LWRK + 1 - KWRK1
210
211      CALL DZERO(WRK,KWRK1)
212
213      IF (IPRPP .GT. 2 .OR. LWRK1 .LT. 2*KZYVAR) THEN
214         WRITE(LUPRI,'(/A)') ' --- IN INVE2VEC:'
215         WRITE(LUPRI,*)' THCPP,MAXRM ',THCPP,MAXRM
216         WRITE(LUPRI,*)' KSYMOP,NGPPP(KSYMOP) ',KSYMOP,NGPPP(KSYMOP)
217         WRITE(LUPRI,*)' LWRK ,LWRK1 ',LWRK,LWRK1
218      END IF
219C
220C ALLOCATE WORK SPACE FOR EIGENVECTORS AND TRANSITION MOMENTS
221C
222      LNEED = 100 + 2*KZYVAR
223C
224C MAXIMUM NUMBER OF SIMULTANEOUS SOLUTION VECTORS
225C
226C      NSIM = MIN(KEXCNV, MAXSIM, (LWRK1-LNEED)/KZYVAR )
227      IF (IPRPP .GT. 2 .OR. NSIM .LE. 0) THEN
228         LWRK2 = KWRK1 + LNEED + KZYVAR
229C        ... need at least space for one KZYVAR (NSIM = 1)
230         WRITE (LUPRI,*) ' KEXCNV,NSIM,LWRK2 ',KEXCNV,NSIM,LWRK2
231         IF (NSIM.LE.0) CALL ERRWRK('RSPPP work space',-LWRK2,LWRK)
232      END IF
233
234      KWRK2  = KWRK1 + NSIM*KZYVAR
235
236      LWRK2  = LWRK   - KWRK2
237
238      THCRSP = THCPP
239      IPRRSP = IPRPP
240      MAXIT  = MAXITP
241
242      CALL DZERO(WRK(KEIVAL),MAXRM)
243      CALL DZERO(WRK(KEIVEC),MAXRM)
244
245
246      DO 1100 ISIM=1,NSIM
247        CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC,.TRUE.,BLANK,
248     *            BLANK,INPVECS((ISIM-1)*KZYVAR+1),WRK(KREDGD),
249     *            WRK(KREDE),WRK(KREDS),
250     *            WRK(KIBTYP),WRK(KEIVAL),WRK(KRESID),WRK(KEIVEC),
251     *            XINDX,WRK(KWRK1),LWRK1)
252        WRITE(LUPRI,*) 'Finished RSPCTL for ISIM=', ISIM
253        CALL RSPEVE(WRK(KIBTYP),WRK(KEIVAL),WRK(KEIVEC),
254     &            OUTVECS((ISIM-1)*KZYVAR+1),
255     &            WRK(KWRK1),NSIM,ISIM-1)
256 1100  CONTINUE
257
258      RETURN
259      END
260C
261C *** END OF INVE2VEC
262C
263
264
265      SUBROUTINE FOLD(AIN,AOUT,N)
266C
267C  This subroutine folds the square matrix AIN into the triangular
268C   matrix AOUT
269C
270#include "implicit.h"
271#include "priunit.h"
272      DIMENSION AIN(N,N), AOUT(*)
273
274      IJ = 0
275      DO I=1,N
276       DO J=1,I-1
277        IJ = IJ + 1
278        AOUT(IJ) = AIN(I,J) + AIN(J,I)
279       END DO
280       IJ = IJ + 1
281       AOUT(IJ) = AIN(I,I)
282      END DO
283
284      RETURN
285      END
286
287      SUBROUTINE FMAT2VEC(NSIM,ZYVEC,ZYMAT)
288#include "implicit.h"
289#include "priunit.h"
290      DIMENSION ZYVEC(KZYWOP,*), ZYMAT(NORBT,NORBT,*)
291C
292#include "maxorb.h"
293#include "maxash.h"
294#include "infvar.h"
295#include "inforb.h"
296#include "infind.h"
297#include "infrsp.h"
298#include "wrkrsp.h"
299C
300      PARAMETER ( D1 = 1.0D0 )
301      CALL DZERO(ZYVEC,NSIM*KZYVAR)
302      DO 100 ISIM=1,NSIM
303         DO 200 IG=1,KZWOPT
304            I=JWOP(1,IG)
305            J=JWOP(2,IG)
306            ZYVEC(IG,ISIM)=ZYMAT(J,I,ISIM)-ZYMAT(I,J,ISIM)
307            ZYVEC(IG+KZWOPT,ISIM)=-ZYVEC(IG,ISIM)
308  200    CONTINUE
309  100 CONTINUE
310
311      RETURN
312      END
313
314      SUBROUTINE SYMMETRIZE(MINP,MOUT,NSIM,NORBT)
315#include "implicit.h"
316#include "priunit.h"
317
318      DIMENSION MINP(NORBT,NORBT,*),MOUT(NORBT,NORBT,*)
319      DOUBLE PRECISION MINP, MOUT
320
321      DO 300 ISIM=1,NSIM
322        DO 200 I=1,NORBT
323          DO 100 J=1,I
324            MOUT(I,J,ISIM) = 0.5D0*(MINP(I,J,ISIM)
325     &                             +MINP(J,I,ISIM))
326            MOUT(J,I,ISIM) = 0.5D0*(MINP(I,J,ISIM)
327     &                             +MINP(J,I,ISIM))
328 100      CONTINUE
329 200    CONTINUE
330 300  CONTINUE
331
332      RETURN
333      END
334
335C
336C  Old AO2MO for inputs without the symmetry
337C
338C      SUBROUTINE AO2MO(MAO,MMO,CMO,NORBT,WRK,LWRK)
339C#include "implicit.h"
340C#include "priunit.h"
341C
342C
343C     Purpose: transforms a two-index matrix from
344C     AO to MO basis (using WRK as scratch )
345C
346C     MMO = CMO^T * MAO * CMO
347C
348C
349C      DIMENSION MAO(*),MMO(*),CMO(*),WRK(*)
350C
351C      CALL DGEMM('T','N',NORBT,NORBT,NORBT,1.D0,
352C     &                   CMO,NORBT,MAO,NORBT,0.D0,
353C     &                   WRK,NORBT)
354C
355C      CALL DGEMM('N','N',NORBT,NORBT,NORBT,1.D0,
356C     &                   WRK,NORBT,CMO,NORBT,0.D0,
357C     &                   MMO,NORBT)
358C
359C      RETURN
360C      END
361C
362C
363C     End of AO2MO_OLD
364C
365
366      SUBROUTINE AO2MO(AAO,BMO,CMO,WRK,LWRK)
367#include "implicit.h"
368#include "priunit.h"
369#include "inforb.h"
370
371C
372C     Purpose: transforms a two-index matrix from
373C     AO to MO basis (using WRK as scratch )
374C
375C     BMO = CMO^T * AAO * CMO
376C
377C  CMO is in rectangular symmetry blocks,
378C  AAO, BMO are full without any symmetry reduction
379C
380
381      DIMENSION AAO(*),BMO(*),CMO(*),WRK(*)
382
383      CALL DZERO(WRK,N2ORBX)
384      CALL DZERO(BMO,N2ORBX)
385
386C      CALL PRINTMAT('AO2MO: beg.    ',1,NORBT,AAO)
387
388      ICOFF = 1
389      DO 100 ISYM=1,NSYM
390        IF ( NBAS(ISYM) .GT. 0 ) THEN
391C          IOFF =  IORB(ISYM)*NORBT+IORB(ISYM) +1
392          IOFF =  IORB(ISYM) + 1
393          CALL DGEMM('T','N',
394C     &      NORB(ISYM),NORB(ISYM),NBAS(ISYM),1.D0,
395     &      NORB(ISYM),NORBT,NBAS(ISYM),1.D0,
396     &      CMO(ICOFF),NBAS(ISYM),
397     &      AAO(IOFF),NORBT,1.0D0,
398     &      WRK(IOFF),NORBT)
399          ICOFF = ICOFF + N2ORB(ISYM)
400        END IF
401 100  CONTINUE
402
403      ICOFF = 1
404      DO 200 ISYM=1,NSYM
405        IF ( NBAS(ISYM) .GT. 0 ) THEN
406C          IOFF =  IORB(ISYM)*NORBT+IORB(ISYM)+1
407          IOFF =  IORB(ISYM)*NORBT+1
408          CALL DGEMM('N','N',
409C     &           NORB(ISYM),NORB(ISYM),NBAS(ISYM),1.0D0,
410     &           NORBT,NORB(ISYM),NBAS(ISYM),1.0D0,
411     &           WRK(IOFF),NORBT,
412     &           CMO(ICOFF),NBAS(ISYM),1.0D0,
413     &           BMO(IOFF),NORBT)
414          ICOFF = ICOFF + N2ORB(ISYM)
415        END IF
416 200  CONTINUE
417
418C      CALL PRINTMAT('AO2MO: final.  ',1,NORBT,BMO)
419
420      RETURN
421      END
422
423C
424C     End of AO2MO
425C
426      SUBROUTINE MO2AO(AMO,AAO,CMO,WRK,LWRK)
427#include "implicit.h"
428#include "priunit.h"
429#include "inforb.h"
430
431C
432C     Purpose: transforms a two-index matrix from
433C     AO to MO basis (using WRK as scratch )
434C
435C     AAO = CMO * AMO * CMO^T
436C
437C  AMO, AAO are full matrix without a symmetry reduction
438C
439C We make a loop in which we multiply CMO block
440C of certain symmetry with the MO matrix
441C
442      DIMENSION AAO(*),AMO(*),CMO(*),WRK(*)
443      INTEGER COFF
444
445      CALL DZERO(WRK,NORBT**2)
446      CALL DZERO(AAO,NORBT**2)
447
448      DO ISYM=1,NSYM
449        IF ( NBAS(ISYM) .GT. 0 ) THEN
450          CALL DGEMM('N','N',
451     &          NBAS(ISYM),NORBT,NORB(ISYM),1.D0,
452     &          CMO(ICMO(ISYM) + 1),NBAS(ISYM),
453     &          AMO(IORB(ISYM) + 1),NORBT,0.D0,
454     &          WRK(IBAS(ISYM)+1),NBAST)
455        END IF
456      END DO
457
458      DO ISYM=1,NSYM
459        IF ( NBAS(ISYM) .GT. 0 ) THEN
460          CALL DGEMM('N','T',
461     &             NORBT,NBAS(ISYM),NORB(ISYM),1.D0,
462     &             WRK(IORB(ISYM)*NBAST+1),NBAST,
463     &             CMO(ICMO(ISYM) + 1),NBAS(ISYM),0.D0,
464     &             AAO(IBAS(ISYM)*NBAST+1),NBAST)
465        END IF
466      END DO
467
468      RETURN
469      END
470
471C
472C This version was working for non-symmetry calc
473C
474C
475C      SUBROUTINE MO2AO_OLD(MMO,MAO,CMO,NORBT,WRK,LWRK)
476C#include "implicit.h"
477C#include "priunit.h"
478C
479C
480C     Purpose: transforms a two-index matrix from
481C     AO to MO basis (using WRK as scratch )
482C
483C     MAO = CMO * MMO * CMO^T
484C
485C
486C      DIMENSION MAO(*),MMO(*),CMO(*),WRK(*)
487C      DOUBLE PRECISION MAO,MMO,CMO,WRK
488C
489C      CALL DGEMM('N','N',NORBT,NORBT,NORBT,1.D0,
490C     &                   CMO,NORBT,MMO,NORBT,0.D0,
491C     &                   WRK,NORBT)
492C
493C      CALL DGEMM('N','T',NORBT,NORBT,NORBT,1.D0,
494C     &                   WRK,NORBT,CMO,NORBT,0.D0,
495C     &                   MAO,NORBT)
496C
497C      RETURN
498C      END
499
500C
501C     End of MO2AO
502C
503
504      SUBROUTINE ESG_PVAL(DODFT,HFXFAC,PVAL,DMAT,A,B,C,D,F)
505#include "implicit.h"
506#include "inforb.h"
507      INTEGER A, B, C, D
508      PARAMETER (DP25 = 0.25 D00, DP5 = 0.5 D00,
509     &           D1 = 1.0 D00, D2 = 2.0 D00, ZERADD = 1.D-15,
510     &           D0 = 0.0 D00, DP125 = 0.125D00, D3=3.0D0, D4=4.0D0)
511      DIMENSION PVAL(*),DMAT(NBAST,NBAST,*)
512      DOUBLE PRECISION PVAL,DMAT
513      LOGICAL DODFT
514
515      IF (.NOT. DODFT) THEN
516         FACTOR = 1.0D0
517      ELSE
518         FACTOR = HFXFAC
519      END IF
520      PVAL(1) = D2*F*(DMAT(A,B,2)*DMAT(D,C,1)+DMAT(A,B,1)*DMAT(D,C,2)
521     & - DP25*FACTOR*(DMAT(C,A,2)*DMAT(D,B,1)+DMAT(D,A,2)*DMAT(C,B,1)+
522     &         DMAT(C,A,1)*DMAT(D,B,2)+DMAT(D,A,1)*DMAT(C,B,2))
523     & + ( DMAT(A,B,4)*DMAT(D,C,4)
524     & - DP125*FACTOR*(DMAT(A,D,3)*DMAT(B,C,3)+DMAT(D,A,3)*DMAT(C,B,3)+
525     &            DMAT(A,C,3)*DMAT(B,D,3)+DMAT(C,A,3)*DMAT(D,B,3))))
526
527      RETURN
528      END
529
530C----------------------------------------
531C     TEMPORARY SUBROUTINES FOR DEBUGING:
532C----------------------------------------
533      SUBROUTINE PRINTVEC( TEXT, NSIM, KZYVAR, V )
534#include "implicit.h"
535#include "priunit.h"
536
537      DIMENSION V(*)
538      CHARACTER*15 TEXT
539
540      WRITE (LUPRI,'(A)') TEXT
541
542      DO 200 I=1,NSIM
543        WRITE( LUPRI, * ) 'VECTOR NUMBER : ', I
544        DO 100 J=1,KZYVAR
545          WRITE( LUPRI, 300 ) I, J,
546     &           V( (I-1)*KZYVAR + J )
547 100    CONTINUE
548 200  CONTINUE
549
550 300  FORMAT(3X,' V',I1,'(',I2,')=',F12.6)
551
552      RETURN
553      END
554C
555C *** END OF PRINTVEC
556C
557      SUBROUTINE PRINTVEC2( TEXT, NSIM, KZVAR, V )
558#include "implicit.h"
559#include "priunit.h"
560
561      DIMENSION V(*)
562      CHARACTER*15 TEXT
563
564      WRITE (LUPRI,'(A)') TEXT
565
566      DO ISIM=1,NSIM
567         IF (NSIM.GT.1) WRITE( LUPRI, * ) 'MATRIX NUMBER : ', ISIM
568         CALL OUTPUT(V,1,KZVAR,1,2,KZVAR,2,1,LUPRI)
569      END DO
570
571c      DO 200 I=1,NSIM
572c        WRITE( LUPRI, * ) 'VECTOR NUMBER : ', I
573c        DO 100 J=1,KZYVAR
574c          WRITE( LUPRI, 300 ) I, J,
575c     &           V( (I-1)*KZYVAR + J )
576c 100    CONTINUE
577c 200  CONTINUE
578
579c 300  FORMAT(3X,' V',I1,'(',I2,')=',F12.6)
580
581      RETURN
582      END
583C
584C *** END OF PRINTVEC
585C
586
587
588      SUBROUTINE PRINTMAT( TEXT, NSIM, NORBT, A )
589#include "implicit.h"
590#include "priunit.h"
591
592
593      DIMENSION A(NORBT,NORBT,*)
594      CHARACTER*15 TEXT
595
596      WRITE (LUPRI,'(A)') TEXT
597      DO ISIM=1,NSIM
598         IF (NSIM.GT.1) WRITE( LUPRI, * ) 'MATRIX NUMBER : ', ISIM
599         CALL OUTPUT(A,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
600      END DO
601c      DO 200 ISIM=1,NSIM
602c        WRITE( LUPRI, * ) 'MATRIX NUMBER : ', ISIM
603c          DO 100 I=1,NORBT
604c            WRITE( LUPRI, 400 )
605c     &             (A(I,J,ISIM),J=1,NORBT)
606c 100      CONTINUE
607c 200  CONTINUE
608c
609c 400  FORMAT(6(3X,F12.6))
610
611      RETURN
612      END
613C
614C *** END OF PRINTVEC
615C
616