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!/* Comdeck rev_log */
21!Revision 1.2  2001/01/17 10:41:29  vebjornb
22!Calls to *MPA*B* in arhpack.F have been replaced with DGEMM calls
23!
24!Revision 1.2  2000/05/01 12:13:01  hjj
25!Removed obsolete KAVER test.
26!
27!Written july / august 1992 Hinne Hettema
28!
29!Purpose of routines :
30!=====================
31!1. perform averaging of start vectors and linear transformations in response
32!   calculations (RSPAVE)
33!2. symmetrize the E[2] matrix (RFANTI ... RSPFCA)
34!
35!List of updates:
36!================
37!===========================================================================
38
39C  /* Deck rspave */
40      SUBROUTINE RSPAVE (VECS,NVDIM,NVSIM)
41C
42C 23-Jul-1992 Hinne Hettema
43C
44#include "implicit.h"
45      DIMENSION VECS(NVDIM,NVSIM)
46C
47C Used from common blocks:
48C  INFINP : SUPSYM
49C  INFVAR : NWOPT,JWOPSY
50C  INFRSP : RSPSUP
51C  WRKRSP : KSYMOP, KZWOPT
52C  INFPRI : IPRAVE
53C
54#include "maxorb.h"
55#include "priunit.h"
56#include "infinp.h"
57#include "infvar.h"
58#include "infpri.h"
59C
60#include "infrsp.h"
61#include "wrkrsp.h"
62C
63      CALL QENTER('RSPAVE')
64C
65      IF ( IPRAVE .GT. 20 ) THEN
66         WRITE(LUPRI,'(A)')    ' Test output from RSPAVE : '
67         WRITE(LUPRI,'(A)')    ' ========================= '
68         WRITE(LUPRI,'(A,I8)') ' NVDIM     =  ', NVDIM
69         WRITE(LUPRI,'(A,I8)') ' NVSIM     =  ', NVSIM
70         WRITE(LUPRI,'(A,I8)') ' NWOPT     =  ', NWOPT
71         WRITE(LUPRI,'(A,I8)') ' JWOPSY    =  ', JWOPSY
72         WRITE(LUPRI,'(A,I8)') ' KSYMOP    =  ', KSYMOP
73         WRITE(LUPRI,'(A,I8)') ' KZWOPT    =  ', KZWOPT
74         WRITE(LUPRI,'(A,L8)') ' SUPSYM    =  ', SUPSYM
75         WRITE(LUPRI,'(A,L8)') ' RSPSUP    =  ', RSPSUP
76         WRITE(LUPRI,'(A)')    ' ========================= '
77      END IF
78      IF (.NOT. SUPSYM) GO TO 9999
79      IF (NWOPT  .EQ. 0) GO TO 9999
80      IF (JWOPSY .NE. 1) GO TO 910
81      IF (NVDIM  .LT. NWOPT) GO TO 920
82      IF (JWOPSY .NE. KSYMOP) GO TO 930
83      IF (NWOPT  .NE. KZWOPT) GO TO 940
84      IF (SUPSYM) THEN
85         IF ( IPRAVE .GT. 80 ) THEN
86            WRITE(LUPRI,'(A)')    ' Vector before averaging : '
87            WRITE(LUPRI,'(A)')    ' ========================= '
88            CALL OUTPUT(VECS,1,NVDIM,1,NVSIM,NVDIM,NVSIM,1,LUPRI)
89         END IF
90         CALL AVERSS(VECS,NVDIM,NVSIM)
91         IF ( IPRAVE .GT. 80 ) THEN
92            WRITE(LUPRI,'(A)')    ' Vector after averaging : '
93            WRITE(LUPRI,'(A)')    ' ======================== '
94            CALL OUTPUT(VECS,1,NVDIM,1,NVSIM,NVDIM,NVSIM,1,LUPRI)
95         END IF
96      END IF
97C
98 9999 CONTINUE
99      IF (IPRAVE .GT. 80) THEN
100         WRITE(LUPRI,'(A)')    ' End of RSPAVE'
101      END IF
102      CALL QEXIT('RSPAVE')
103      RETURN
104C
105C *** Error sections
106C
107  910 CONTINUE
108      WRITE (LUERR,9010) JWOPSY
109      CALL QTRACE(LUERR)
110      CALL QUIT('RSPAVE ERROR, operator symmetry JWOPSY .ne. 1')
111 9010 FORMAT(/' ERROR-RSPAVE, operator symmetry .ne. 1; JWOPSY =',I8)
112C
113  920 CONTINUE
114      WRITE (LUERR,9020) NVDIM,NWOPT
115      CALL QTRACE(LUERR)
116      CALL QUIT('RSPAVE ERROR, vector length lt NWOPT')
117 9020 FORMAT(/' ERROR-RSPAVE, vector length',I8,' is less than NWOPT',
118     *   I8)
119C
120  930 CONTINUE
121      WRITE (LUERR,9020) JWOPSY, KSYMOP
122      CALL QTRACE(LUERR)
123      CALL QUIT('RSPAVE ERROR, JWOPSY .NE. KSYMOP')
124 9030 FORMAT(/' ERROR-RSPAVE, JWOPSY = ',I8,' KSYMOP = ', I8)
125C
126  940 CONTINUE
127      WRITE (LUERR,9020) NWOPT, KZWOPT
128      CALL QTRACE(LUERR)
129      CALL QUIT('RSPAVE ERROR, NWOPT .NE. KZWOPT')
130 9040 FORMAT(/' ERROR-RSPAVE, NWOPT = ',I8,' KZWOPT = ', I8)
131C
132      END
133C  /* Deck rfanti */
134      SUBROUTINE RFANTI(NSIM,EVECS,ZYMAT,WRK,LWRK)
135C
136C Subroutine to obtain antisymmetric part of the
137C MCSCF gradient from total MCSCF Fock matrix
138C
139#include "implicit.h"
140C
141      DIMENSION EVECS(*),ZYMAT(NORBT,NORBT,NSIM),WRK(LWRK)
142C
143C Used from common blocks
144C
145#include "maxorb.h"
146#include "priunit.h"
147#include "infinp.h"
148#include "inftap.h"
149#include "infopt.h"
150#include "infpri.h"
151#include "inforb.h"
152#include "infvar.h"
153#include "infrsp.h"
154#include "wrkrsp.h"
155C
156      PARAMETER ( DHALF = 0.50D0 )
157C
158      CALL QENTER('RFANTI')
159C
160C     Layout work space
161C
162      KGRA   = 1
163      KFOCK  = KGRA    + NSIM * N2ORBX
164      KANTI  = KFOCK   + N2ORBT
165      KWRK   = KANTI   + N2ORBX
166      LWRK1  = LWRK    - KWRK
167C
168      IF (LWRK1 .LT. 0 ) THEN
169         CALL QTRACE(LUERR)
170         CALL QUIT(' Not enough work space in RFANTI ')
171      END IF
172C
173C  -- Get the Fock matrix from the interface file:
174C     Comment : we know that RSPCI is False, so no test.
175C
176      REWIND(LUSIFC)
177      CALL MOLLAB(LBSIFC,LUSIFC,LUERR)
178      READ (LUSIFC)
179      READ (LUSIFC)
180      READ (LUSIFC)
181      READ (LUSIFC)
182      READ (LUSIFC)
183      CALL READT(LUSIFC,N2ORBT,WRK(KFOCK))
184      IF (IPRAVE .GT. 80 ) THEN
185         WRITE(LUPRI,'(//A)')  ' Fock matrix for use in RFANTI'
186         WRITE(LUPRI,'(A)')    ' ============================='
187         DO 10 ISYM = 1, NSYM
188            I2ORBI = I2ORB(ISYM)
189            NORBI  = NORB(ISYM)
190            WRITE(LUPRI,'(/A,I4)')  ' Symmetry representation: ', ISYM
191            WRITE(LUPRI,'(A)')      ' =============================='
192            CALL OUTPUT(WRK(KFOCK+I2ORBI),1,NORBI,1,NORBI,NORBI,NORBI,
193     &                  1,LUPRI)
194 10      CONTINUE
195      END IF
196C
197C  -- Get the antisymmetric component of the Fock matrix
198C
199      CALL DZERO(WRK(KANTI), NSIM*N2ORBX)
200      DO 100 ISYM = 1, NSYM
201         NORBI  = NORB(ISYM)
202         IORBI  = IORB(ISYM)
203         I2ORBI = I2ORB(ISYM)
204         CALL RSPANT(WRK(KFOCK+I2ORBI),NORBI,IORBI,NORBT,WRK(KANTI))
205  100 CONTINUE
206C
207      IF (IPRAVE .GT.  80) THEN
208         WRITE(LUPRI,'(//A)') ' Antisymmetry component'
209         WRITE(LUPRI,'(A)')   ' ======================'
210         CALL OUTPUT(WRK(KANTI),1,NORBT,1,NORBT,NORBT,NORBT,
211     &               1,LUPRI)
212         DO 90 ISIM = 1, NSIM
213            WRITE(LUPRI,'(//A)') ' ZYMAT matrix used :'
214            WRITE(LUPRI,'(A)')   ' ==================='
215            CALL OUTPUT(ZYMAT(1,1,ISIM),1,NORBT,1,NORBT,NORBT,NORBT,
216     &                  1,LUPRI)
217  90     CONTINUE
218      END IF
219C
220C  -- Do matrix multiplication with ZYMAT to get correction
221C     0.5 sum(s) ( k(q,s) g(p,s) - g(s,q) k(s,p) )
222C
223      DO 200 ISIM = 1, NSIM
224         CALL DGEMM('N','T',NORBT,NORBT,NORBT,1.D0,
225     &              ZYMAT(1,1,ISIM),NORBT,
226     &              WRK(KANTI),NORBT,0.D0,
227     &              WRK(KGRA+N2ORBX*(ISIM-1)),NORBT)
228         CALL DGEMM('T','N',NORBT,NORBT,NORBT,-1.D0,
229     &              WRK(KANTI),NORBT,
230     &              ZYMAT(1,1,ISIM),NORBT,1.D0,
231     &              WRK(KGRA+N2ORBX*(ISIM-1)),NORBT)
232         IF (IPRAVE .GT.  80) THEN
233            WRITE(LUPRI,'(//A,I4)') ' ISIM = ', ISIM
234            WRITE(LUPRI,'(A)') ' Multiplication on ZYMAT'
235            WRITE(LUPRI,'(A)') ' ======================='
236            CALL OUTPUT(WRK(KGRA+N2ORBX*(ISIM-1)),1,NORBT,1,
237     &                  NORBT,NORBT,NORBT,1,LUPRI)
238         END IF
239  200 CONTINUE
240      CALL DSCAL(NSIM*N2ORBX,DHALF,WRK(KGRA),1)
241C
242C     and correct the gradient
243C
244      CALL RSPCFA(NSIM,WRK(KGRA),EVECS)
245C
246      CALL QEXIT('RFANTI')
247C
248      RETURN
249      END
250C  /* Deck rspant */
251      SUBROUTINE RSPANT(FOCK,NORB,IOFF,NORBT,GMAT)
252C
253C Subroutine to compute antisymmetric part of Fock matrix
254C
255#include "implicit.h"
256      DIMENSION FOCK(NORB,NORB), GMAT(NORBT,NORBT)
257C
258      CALL QENTER('RSPANT')
259C
260      DO 100 I = 1, NORB
261         DO 200 J = 1, NORB
262            GMAT(IOFF+I,IOFF+J) = FOCK(I,J) - FOCK(J,I)
263            GMAT(IOFF+J,IOFF+I) = FOCK(J,I) - FOCK(I,J)
264200      CONTINUE
265100   CONTINUE
266C
267      CALL QEXIT('RSPANT')
268C
269      RETURN
270      END
271C  /* Deck rspcfa */
272      SUBROUTINE RSPCFA(NSIM,GRA,EVECS)
273#include "implicit.h"
274C
275      DIMENSION GRA(NORBT,NORBT,*), EVECS(KZYVAR,*)
276C
277#include "maxorb.h"
278#include "maxash.h"
279#include "infvar.h"
280#include "inforb.h"
281#include "infind.h"
282#include "infdim.h"
283#include "infpri.h"
284#include "infrsp.h"
285#include "wrkrsp.h"
286C
287      CALL QENTER('RSPCFA')
288C
289      KYCONF = KZCONF + KZVAR
290C
291C Distribute gradient matrices into sigma vector in EVECS
292C
293      DO 100 ISIM = 1, NSIM
294         DO 200 IG = 1,KZWOPT
295            K     = JWOP(1,IG)
296            L     = JWOP(2,IG)
297            EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
298     *            + GRA(L,K,ISIM)
299            EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
300     *            + GRA(K,L,ISIM)
301 200    CONTINUE
302 100  CONTINUE
303C
304      CALL QEXIT('RSPCFA')
305C
306      RETURN
307      END
308