1!
2!...   Copyright (c) 2014 by the authors of Dalton (see below).
3!...   All Rights Reserved.
4!...
5!...   The source code in this file is part of
6!...   "Dalton, a molecular electronic structure program,
7!...    Release DALTON2014 (2014), see http://daltonprogram.org"
8!...
9!...   This source code is provided under a written licence and may be
10!...   used, copied, transmitted, or stored only in accord with that
11!...   written licence.
12!...
13!...   In particular, no part of the source code or compiled modules may
14!...   be distributed outside the research group of the licence holder.
15!...   This means also that persons (e.g. post-docs) leaving the research
16!...   group of the licence holder may not take any part of Dalton,
17!...   including modified files, with him/her, unless that person has
18!...   obtained his/her own licence.
19!...
20!...   For further information, including how to get a licence, see:
21!...      http://daltonprogram.org
22!
23!
24
25
26C*****************************************************************************
27      SUBROUTINE GETRHO_OLD(DMAT,GSO,RHO,DMATGAO,DFTHRI,IPRINT)
28C*****************************************************************************
29C
30C     T. Helgaker feb 01
31C     (RHO13 removed Aug 18)
32C
33C Output:
34C    RHO
35C    DMATGAO(i) = sum(j) dmat(i,j) * gso(j)
36C
37C*****************************************************************************
38      implicit none
39#include "priunit.h"
40#include "maxaqn.h"
41#include "maxorb.h"
42#include "mxcent.h"
43#include "symmet.h"
44#include "inforb.h"
45      real*8, parameter :: D0 = 0.0d0
46      real*8, parameter :: D1 = 1.0d0
47      real*8, parameter :: D2 = 2.0d0
48      real*8, parameter :: DP5 = 0.5d0
49      real*8, parameter :: DP3 = D1/3.0d0
50C
51      integer, intent(in) :: IPRINT
52      real*8, intent(in) :: DMAT(NBAST,NBAST), GSO(NBAST), DFTHRI
53      real*8, intent(out) :: RHO, DMATGAO(NBAST)
54C
55      integer :: ISTR, ISYM, NBASI
56      real*8 :: DDOT
57C
58      IF (IPRINT .GT. 100) THEN
59         write(lupri,*) 'GETRHO_OLD NBAST, NSYM',NBAST,NSYM
60         write(lupri,*) 'DFTHRI, DMAT(1,1)',DFTHRI,DMAT(1,1),GSO(1)
61         write(lupri,*) 'IBAS',ibas(1:nsym)
62         write(lupri,*) 'NBAS',nbas(1:nsym)
63      END IF
64      IF (NSYM.EQ.1) THEN
65         CALL DSYMV('U',NBAST,D1,DMAT,NBAST,GSO,1,D0,DMATGAO,1)
66      ELSE
67         CALL DZERO(DMATGAO,NBAST)
68         DO ISYM = 1, NSYM
69            ISTR = IBAS(ISYM) + 1
70            NBASI= NBAS(ISYM)
71            CALL DSYMV('U',NBASI,D1,DMAT(ISTR,ISTR),NBAST,GSO(ISTR),1,
72     &                 D0,DMATGAO(ISTR),1)
73         END DO
74      END IF
75      RHO = DDOT(NBAST,GSO,1,DMATGAO,1)
76      RETURN
77      END
78
79
80C*****************************************************************************
81      SUBROUTINE DFTKSMGGASPIN(EXCMAT, GAO, GAO1, RHG, d1E, DFTHRL)
82C*****************************************************************************
83C
84C     Erik Hedegaard (feb. 2016) based on DFTKSM by T. Helgaker
85C     Revised Apr. 2018 by Hans Joergen Aa. Jensen
86C
87C     RHG(1:3,1:4): grad(rhoc), grad(rhos), grad(rhoa), grad(rhob)
88C     VXC  = d1E(1) = d(e_xc) / d(rhoc)
89C     VSC  = d1E(2) = d(e_xc) / d(rhos)
90C     VXB  = d1E(3) = d(e_xc) / d(grdcc)
91C     VSB  = d1E(4) = d(e_xc) / d(grdss)
92C     VXSB = d1E(5) = d(e_xc) / d(grdcs)
93C
94C*****************************************************************************
95      implicit none
96! inforb.h : NBAST, NSYM
97#include "inforb.h"
98      real*8, intent(in) :: GAO(NBAST), GAO1(NBAST,3), RHG(3,2)
99      real*8, intent(in) :: d1E(5), DFTHRL
100      real*8, intent(inout) :: EXCMAT(NBAST,NBAST,2)
101
102      ! local variables:
103      real*8 :: VXB, VXC, VSB, VSC, VXSB
104      real*8 :: FX, FY, FZ, SX, SY, SZ
105      real*8 :: FC, FCS
106      integer :: I, J, ISYM, ISTR, IEND
107C
108C     Exchange-correlation contribution to Kohn-Sham matrix
109C
110      VXC  = d1E(1)
111      VSC  = d1E(2)
112      VXB  = d1E(3) * 4.0d0
113      VSB  = d1E(4) * 4.0d0
114      VXSB = d1E(5) * 2.0d0
115      FX = VXB*RHG(1,1) + VXSB*RHG(1,2)
116      FY = VXB*RHG(2,1) + VXSB*RHG(2,2)
117      FZ = VXB*RHG(3,1) + VXSB*RHG(3,2)
118      SX = VSB*RHG(1,2) + VXSB*RHG(1,1)
119      SY = VSB*RHG(2,2) + VXSB*RHG(2,1)
120      SZ = VSB*RHG(3,2) + VXSB*RHG(3,1)
121      IF (NSYM.EQ.1) THEN
122         DO J = 1, NBAST
123            FC = VXC*GAO(J)+FX*GAO1(J,1)+FY*GAO1(J,2)+FZ*GAO1(J,3)
124            IF (abs(FC).GT.DFTHRL) THEN
125               CALL DAXPY(NBAST,FC,GAO,1,EXCMAT(1,J,1),1)
126            END IF
127            FCS = VSC*GAO(J)+SX*GAO1(J,1)+SY*GAO1(J,2)+SZ*GAO1(J,3)
128            IF (abs(FCS).GT.DFTHRL) THEN
129               CALL DAXPY(NBAST,FCS,GAO,1,EXCMAT(1,J,2),1)
130            END IF
131         END DO
132      ELSE
133         DO ISYM = 1, NSYM
134            ISTR = IBAS(ISYM) + 1
135            IEND = IBAS(ISYM) + NBAS(ISYM)
136            DO J = ISTR, IEND
137               FC = VXC*GAO(J)+FX*GAO1(J,1)+FY*GAO1(J,2)+FZ*GAO1(J,3)
138               IF (abs(FC).GT.DFTHRL) THEN
139                  DO I = ISTR, IEND
140                     EXCMAT(I,J,1) = EXCMAT(I,J,1) + FC*GAO(I)
141                  END DO
142               END IF
143               FCS = VSC*GAO(J)+SX*GAO1(J,1)+SY*GAO1(J,2)+SZ*GAO1(J,3)
144               IF (abs(FCS).GT.DFTHRL) THEN
145                  DO I = ISTR, IEND
146                     EXCMAT(I,J,2) = EXCMAT(I,J,2) + FCS*GAO(I)
147                  END DO
148               END IF
149            END DO
150         END DO
151      END IF
152      RETURN
153      END
154
155
156C*****************************************************************************
157      SUBROUTINE DFTKSM_MGGA(EXCMAT, GAO, GAO1, GAO2, RHG, d1E,
158     &                           DFTHRL, DO_SPIN)
159C*****************************************************************************
160C
161C     Erik Kjellgren (oct. 2018) based on DFTKSMGGASPIN by E. Hedegaard
162C
163C     RHG(1:3,1:4): grad(rhoc), grad(rhos), grad(rhoa), grad(rhob)
164C     VXC  = d1E(1) = d(e_xc) / d(rhoc)
165C     VSC  = d1E(2) = d(e_xc) / d(rhos)
166C     VXB  = d1E(3) = d(e_xc) / d(grdcc)
167C     VSB  = d1E(4) = d(e_xc) / d(grdss)
168C     VXSB = d1E(5) = d(e_xc) / d(grdcs)
169C     VXT  = d1E(6) = d(e_xc) / d(tauc)
170C     VST  = d1E(7) = d(e_xc) / d(taus)
171C     VXL  = d1E(8) = d(e_xc) / d(laplace(rhoc))
172C     VSL  = d1E(9) = d(e_xc) / d(laplace(rhos))
173C
174C     GAO2 is just a dummy variable for now, since it is not calculated
175C     in the code. It is put in now to make it easier to implement MGGA
176C     functional in the future that depends on laplace(rho). GAO2 is
177C     supposed to be defined as:
178C     GAO2 = laplace(Omega) = lapalce(GAOp)*GAOq + GAOp*laplace(GAOq)
179C                              + 2*GRAD(GAOp)*GRAD(GAOq)
180C     When calling this subroutine just set GAO2 = 0.0d0
181C     Later remember to change it to a vector when actually need in the
182C     declaration.
183C
184C*****************************************************************************
185      implicit none
186! inforb.h : NBAST, NSYM
187#include "inforb.h"
188      logical, intent(in) :: DO_SPIN
189      real*8, intent(in) :: GAO(NBAST), GAO1(NBAST,3), RHG(3,2)
190      real*8, intent(in) :: GAO2, d1E(9), DFTHRL
191      real*8, intent(inout) :: EXCMAT(NBAST,NBAST,2)
192
193      ! local variables:
194      real*8 :: VXB, VXC, VSB, VSC, VXSB, VXT, VST
195      !real*8 :: VXL, VSL
196      real*8 :: FX, FY, FZ, SX, SY, SZ
197      real*8 :: FC, FCS
198      integer :: I, J, ISYM, ISTR, IEND
199C
200C     Exchange-correlation contribution to Kohn-Sham matrix
201C
202      VXC  = d1E(1)
203      VSC  = d1E(2)
204      VXB  = d1E(3) * 4.0d0
205      VSB  = d1E(4) * 4.0d0
206      VXSB = d1E(5) * 2.0d0
207      VXT  = d1E(6) * 0.5d0
208      VST  = d1E(7) * 0.5d0
209      !VXL  = d1E(8), no laplace(rho) implementation yet
210      !VSL  = d1E(9), no laplace(rho) implementation yet
211      FX = VXB*RHG(1,1) + VXSB*RHG(1,2)
212      FY = VXB*RHG(2,1) + VXSB*RHG(2,2)
213      FZ = VXB*RHG(3,1) + VXSB*RHG(3,2)
214      SX = VSB*RHG(1,2) + VXSB*RHG(1,1)
215      SY = VSB*RHG(2,2) + VXSB*RHG(2,1)
216      SZ = VSB*RHG(3,2) + VXSB*RHG(3,1)
217      IF (NSYM.EQ.1) THEN
218         DO J = 1, NBAST
219            FC = VXC*GAO(J)+FX*GAO1(J,1)+FY*GAO1(J,2)+FZ*GAO1(J,3)
220            IF (abs(FC).GT.DFTHRL) THEN
221               CALL DAXPY(NBAST,FC,GAO,1,EXCMAT(1,J,1),1)
222            END IF
223            ! Meta-GGA tau part of FC
224            DO I = 1, NBAST
225               EXCMAT(I,J,1) = EXCMAT(I,J,1)
226     &                         + VXT * (GAO1(J,1)*GAO1(I,1)
227     &                                 + GAO1(J,2)*GAO1(I,2)
228     &                                 + GAO1(J,3)*GAO1(I,3))
229C     &                        + VXL * (GAO2(J)*GAO(I)
230C     &                               + GAO(J)*GAO2(I)
231C     &                               + 2.0d0*GAO1(J,1)*GAO1(I,1)
232C     &                               + 2.0d0*GAO1(J,2)*GAO1(I,2)
233C     &                               + 2.0d0*GAO1(J,3)*GAO1(I,3)
234            END DO
235            IF (DO_SPIN) THEN
236               FCS = VSC*GAO(J)+SX*GAO1(J,1)+SY*GAO1(J,2)+SZ*GAO1(J,3)
237               IF (abs(FCS).GT.DFTHRL) THEN
238                  CALL DAXPY(NBAST,FCS,GAO,1,EXCMAT(1,J,2),1)
239               END IF
240               ! Meta-GGA tau part of FCS
241               DO I = 1, NBAST
242                  EXCMAT(I,J,2) = EXCMAT(I,J,2)
243     &                           + VST * (GAO1(J,1)*GAO1(I,1)
244     &                                   + GAO1(J,2)*GAO1(I,2)
245     &                                   + GAO1(J,3)*GAO1(I,3))
246C     &                           + VSL * (GAO2(J)*GAO(I)
247C     &                                  + GAO(J)*GAO2(I)
248C     &                                  + 2.0d0*GAO1(J,1)*GAO1(I,1)
249C     &                                  + 2.0d0*GAO1(J,2)*GAO1(I,2)
250C     &                                  + 2.0d0*GAO1(J,3)*GAO1(I,3)
251              END DO
252            END IF
253         END DO
254      ELSE
255         DO ISYM = 1, NSYM
256            ISTR = IBAS(ISYM) + 1
257            IEND = IBAS(ISYM) + NBAS(ISYM)
258            DO J = ISTR, IEND
259               FC = VXC*GAO(J)+FX*GAO1(J,1)+FY*GAO1(J,2)+FZ*GAO1(J,3)
260               IF (abs(FC).GT.DFTHRL) THEN
261                  DO I = ISTR, IEND
262                     EXCMAT(I,J,1) = EXCMAT(I,J,1) + FC*GAO(I)
263                  END DO
264               END IF
265               DO I = ISTR, IEND
266                  EXCMAT(I,J,1) = EXCMAT(I,J,1)
267     &                            + VXT * (GAO1(J,1)*GAO1(I,1)
268     &                                    + GAO1(J,2)*GAO1(I,2)
269     &                                    + GAO1(J,3)*GAO1(I,3))
270C     &                           + VXL * (GAO2(J)*GAO(I)
271C     &                                  + GAO(J)*GAO2(I)
272C     &                                  + 2.0d0*GAO1(J,1)*GAO1(I,1)
273C     &                                  + 2.0d0*GAO1(J,2)*GAO1(I,2)
274C     &                                  + 2.0d0*GAO1(J,3)*GAO1(I,3)
275               END DO
276               IF (DO_SPIN) THEN
277                 FCS = VSC*GAO(J)+SX*GAO1(J,1)+SY*GAO1(J,2)+SZ*GAO1(J,3)
278                 IF (abs(FCS).GT.DFTHRL) THEN
279                    DO I = ISTR, IEND
280                       EXCMAT(I,J,2) = EXCMAT(I,J,2) + FCS*GAO(I)
281                    END DO
282                 END IF
283                 DO I = ISTR, IEND
284                   EXCMAT(I,J,2) = EXCMAT(I,J,2)
285     &                             + VST * (GAO1(J,1)*GAO1(I,1)
286     &                                     + GAO1(J,2)*GAO1(I,2)
287     &                                     + GAO1(J,3)*GAO1(I,3))
288C     &                            + VSL * (GAO2(J)*GAO(I)
289C     &                                   + GAO(J)*GAO2(I)
290C     &                                   + 2.0d0*GAO1(J,1)*GAO1(I,1)
291C     &                                   + 2.0d0*GAO1(J,2)*GAO1(I,2)
292C     &                                   + 2.0d0*GAO1(J,3)*GAO1(I,3)
293                 END DO
294               END IF
295            END DO
296         END DO
297      END IF
298      RETURN
299      END
300
301
302C*****************************************************************************
303      SUBROUTINE DFTKSM(EXCMAT,GAO,GAO1,RHG,VXC,VXB,DOGGA,FROMVX,DFTHRL)
304C*****************************************************************************
305C
306C     T. Helgaker oct 2000
307C
308C     VXC = d(e_xc) / d(rhoc)
309C     VXB = d(e_xc) / d(grdcc)
310C
311C*****************************************************************************
312      implicit none
313#include "inforb.h"
314      real*8, parameter :: D2 = 2.0D0
315C
316      logical, intent(in) :: DOGGA, FROMVX
317      real*8, intent(in) :: GAO(NBAST), GAO1(NBAST,3), RHG(3),
318     &                      VXC, VXB, DFTHRL
319      real*8, intent(inout) :: EXCMAT(NBAST,NBAST)
320C
321      integer :: I, IEND, ISTR, ISYM, J
322      real*8 :: FC, FX, FY, FZ, GJ, GVXC
323C
324C     Exchange-correlation contribution to Kohn-Sham matrix
325C
326      IF (DOGGA .AND. .NOT.FROMVX) THEN
327         FX = 4.0d0*VXB*RHG(1)
328         FY = 4.0d0*VXB*RHG(2)
329         FZ = 4.0d0*VXB*RHG(3)
330         IF (NSYM.EQ.1) THEN
331            DO J = 1, NBAST
332               FC = VXC*GAO(J)+FX*GAO1(J,1)+FY*GAO1(J,2)+FZ*GAO1(J,3)
333               IF (abs(FC).GT.DFTHRL) THEN
334                  CALL DAXPY(NBAST,FC,GAO,1,EXCMAT(1,J),1)
335               END IF
336            END DO
337         ELSE
338            DO ISYM = 1, NSYM
339               ISTR = IBAS(ISYM) + 1
340               IEND = IBAS(ISYM) + NBAS(ISYM)
341               DO J = ISTR, IEND
342                  FC = VXC*GAO(J)+FX*GAO1(J,1)+FY*GAO1(J,2)+FZ*GAO1(J,3)
343                  IF (abs(FC).GT.DFTHRL) THEN
344                     DO I = ISTR, IEND
345                        EXCMAT(I,J) = EXCMAT(I,J) + FC*GAO(I)
346                     END DO
347                  END IF
348               END DO
349            END DO
350         END IF
351      ELSE
352         DO ISYM = 1, NSYM
353            ISTR = IBAS(ISYM) + 1
354            IEND = IBAS(ISYM) + NBAS(ISYM)
355            DO J = ISTR, IEND
356               GJ   = GAO(J)
357               GVXC = D2*VXC*GJ
358               IF (abs(GVXC).GT.DFTHRL) THEN
359                  DO I = ISTR, J - 1
360                     EXCMAT(I,J) = EXCMAT(I,J) + GVXC*GAO(I)
361                  END DO
362                  EXCMAT(J,J) = EXCMAT(J,J) + VXC*GJ*GJ
363               END IF
364            END DO
365         END DO
366      END IF
367      RETURN
368      END
369
370
371C*****************************************************************************
372      SUBROUTINE DFTFRC(DMAT,GAO,GAO1,GAO2,VXC,VXB,RHX,RHY,RHZ,DOGGA)
373C*****************************************************************************
374C
375C     Exchange-correlation contribution to molecular gradient
376C
377C     T. Helgaker sep 99/oct 00/feb 01
378C
379C*****************************************************************************
380#include "implicit.h"
381#include "priunit.h"
382#include "maxaqn.h"
383#include "maxorb.h"
384#include "mxcent.h"
385C
386      PARAMETER (D0 = 0.0D0, D2 = 2.0D0)
387C
388      LOGICAL DOGGA
389      REAL*8  DMAT(NBAST,NBAST),
390     &        GAO(NBAST), GAO1(NBAST,3), GAO2(NBAST,6)
391C
392#include "inforb.h"
393#include "energy.h"
394#include "symmet.h"
395#include "shells.h"
396C
397      V2 = D2*VXC
398      DO IX = 1, 3
399         IF (IX.EQ.1) THEN
400            K1 = 1
401            K2 = 2
402            K3 = 3
403         ELSE IF (IX.EQ.2) THEN
404            K1 = 2
405            K2 = 4
406            K3 = 5
407         ELSE
408            K1 = 3
409            K2 = 5
410            K3 = 6
411         END IF
412         DO IREPA = 0, MAXREP
413            ISTR = IBAS(IREPA+1) + 1
414            IEND = IBAS(IREPA+1) + NBAS(IREPA+1)
415            IRPAX = IEOR(IREPA,ISYMAX(IX,1))
416            IORBA = 0
417            DO ISHELA = 1, KMAX
418               ISCOOR = IPTCNT(3*(NCENT(ISHELA) - 1) + IX,0,1)
419               DO ICOMPA = 1, KHKT(ISHELA)
420                  IORBA = IORBA + 1
421                  IA = IPTSYM(IORBA,IREPA)
422                  KA = IPTSYM(IORBA,IRPAX)
423                  IF (KA.GT.0) THEN
424                     IF (DOGGA) THEN
425                        GA  = GAO1(KA,IX)
426                        GAX = RHX*GA
427                        GAY = RHY*GA
428                        GAZ = RHZ*GA
429                        GA2 = RHX*GAO2(KA,K1) + RHY*GAO2(KA,K2)
430     &                                        + RHZ*GAO2(KA,K3)
431                        GD = D0
432                        GF = D0
433                        DO IB = ISTR, IEND
434                           GD = GD + DMAT(IB,IA)*GAO(IB)
435                           GF = GF + DMAT(IB,IA)*(GAO(IB)*GA2
436     &                                          + GAO1(IB,1)*GAX
437     &                                          + GAO1(IB,2)*GAY
438     &                                          + GAO1(IB,3)*GAZ)
439                        END DO
440                        FRC = V2*GD*GA + VXB*GF
441                     ELSE
442                        GD = D0
443                        DO IB = ISTR, IEND
444                           GD = GD + GAO(IB)*DMAT(IB,IA)
445                        END DO
446                        FRC = V2*GD*GAO1(KA,IX)
447                     END IF
448                     GRADFT(ISCOOR) = GRADFT(ISCOOR) - FRC
449                  END IF
450               END DO
451            END DO
452         END DO
453      END DO
454      RETURN
455      END
456
457
458C*****************************************************************************
459      SUBROUTINE DFTLTR(KSYMOP,DTRMAT,EXCMAT,GAO,GAO1,C0,C1,C2,HES,
460     &                  RHO,RHO13,RHOGRD,RHG,WDRC,WVWN,WBCK,WLYP,DODRC,
461     &                  DOVWN,DOBCK,DOLYP,DOGGA,TRPLET,DOHES,DTGAO)
462C*****************************************************************************
463C
464C     T. Helgaker sep 99/oct 00
465C
466C*****************************************************************************
467#include "implicit.h"
468#include "priunit.h"
469#include "mxcent.h"
470C
471      PARAMETER (D0 = 0.0D0, D1 = 1.0D0, DP5 = 0.5D0)
472C
473#include "inforb.h"
474#include "nuclei.h"
475#include "dftcom.h"
476C
477      INTEGER A, B
478      LOGICAL DODRC, DOVWN, DOBCK, DOLYP, DOGGA, TRPLET, DOHES
479      DIMENSION DTRMAT(NBAST,NBAST),
480     &          GAO(NBAST), GAO1(NBAST,3),
481     &          EXCMAT(NBAST,NBAST),
482     &          C0(NORBT), C1(NORBT,3), C2(NORBT),
483     &          HES(NOCCT,NVIRT,NOCCT,NVIRT),RHG(3),DTGAO(NBAST)
484      DIMENSION B3(3)
485C
486      IF (DOHES .AND. NSYM.NE.1) THEN
487         WRITE (LUPRI,'(2X,A,/A)')
488     &      ' Symmetry not implemented for explicit Hessian in DFTLTR.',
489     &      ' Program aborted.'
490         CALL QUIT('Symmetry not implementd for DOHES in DFTLTR')
491      END IF
492      IF (DOHES) THEN
493         B0 = D1
494      ELSE
495         CALL DGEMV('N',NBAST,NBAST,D1,DTRMAT,NBAST,GAO,1,D0,DTGAO,1)
496         B0 = DDOT(NBAST,DTGAO,1,GAO,1)
497      END IF
498C
499C     ***************
500C     ***** LDA *****
501C     ***************
502C
503      IF (.NOT.DOGGA) THEN
504         IF (abs(B0).GT.DFTHRL) THEN
505C
506C           Calculated VT
507C
508            IF (DODRC) THEN
509               CALL V1DRC(VDRC0,VDRC1,RHO,RHO13)
510            ELSE
511               VDRC0 = 0D0
512               VDRC1 = 0D0
513            ENDIF
514            IF (DOVWN) THEN
515               IF (.NOT.TRPLET) THEN
516                  CALL V1VWN(VVWN0,FRRVWN,RHO,RHO13)
517               ELSE
518                  CALL VTVWN(FRRVWN,RHO,RHO13)
519                  FRRVWN = DP5*FRRVWN
520               END IF
521            ELSE
522               VVWN0 = 0D0
523               FRRVWN = 0D0
524            END IF
525            VT = (WDRC*VDRC1 + WVWN*FRRVWN)*B0
526C
527C           Linear transformation
528C
529            IF (.NOT.DOHES) THEN
530               IF (NSYM.EQ.1) THEN
531                  DO I = 1, NBAST
532                     GVI = VT*GAO(I)
533                     DO J = 1, I
534                        EXCMAT(J,I) = EXCMAT(J,I) + GVI*GAO(J)
535                     END DO
536                  END DO
537               ELSE
538                  DO ISYM = 1, NSYM
539                     ISTR = IBAS(ISYM) + 1
540                     IEND = IBAS(ISYM) + NBAS(ISYM)
541                     JSYM = MULD2H(ISYM,KSYMOP)
542                     IF (ISYM.GE.JSYM) THEN
543                        JSTR = IBAS(JSYM) + 1
544                        JEND = IBAS(JSYM) + NBAS(JSYM)
545                        DO I = ISTR, IEND
546                           GVI = VT*GAO(I)
547                           DO J = JSTR, MIN(I,JEND)
548                              EXCMAT(J,I) = EXCMAT(J,I) + GVI*GAO(J)
549                           END DO
550                        END DO
551                     END IF
552                  END DO
553               END IF
554C
555C           Explicit Hessian
556C
557            ELSE
558               CALL DGEMM('T','N',NORBT,1,NBAST,1.D0,
559     &                    DTRMAT,NBAST,
560     &                    GAO,NBAST,0.D0,
561     &                    C0,NORBT)
562               DO 300 B = NOCCT + 1, NORBT
563               DO 300 J = 1, NOCCT
564               DO 300 A = NOCCT + 1, NORBT
565               DO 300 I = 1, NOCCT
566                  IA = A - NOCCT
567                  IB = B - NOCCT
568                  HES(I,IA,J,IB) = HES(I,IA,J,IB)
569     &                           + VT*C0(A)*C0(I)*C0(B)*C0(J)
570  300          CONTINUE
571            END IF
572         END IF
573      ELSE
574C
575C        ***************
576C        ***** GGA *****
577C        ***************
578C
579C        B0, BX=B3(1), BY=B3(2), BZ=B3(3), BMAX
580C
581         IF (DOHES) THEN
582            BMAX = D1
583         ELSE
584C           B3 = GAO1'*DTGAO
585            CALL DGEMV('T',NBAST,3,D1,GAO1,NBAST,DTGAO,1,D0,B3,1)
586C           DTGAO= DTRMAT'*GAO
587            CALL DGEMV('T',NBAST,NBAST,D1,DTRMAT,NBAST,GAO,1,D0,DTGAO,1)
588C           B3 = B3 + GAO1'*DTGAO
589            CALL DGEMV('T',NBAST,3,D1,GAO1,NBAST,DTGAO,1,D1,B3,1)
590            BMAX = MAX(abs(B0),abs(B3(1)),abs(B3(2)),abs(B3(3)))
591         END IF
592C
593         IF (BMAX.GT.DFTHRL) THEN
594C
595C           ZNV, FZ0, FRR, FRZ, FZZ
596C
597            IF (DODRC) CALL V1DRC(VDRC0,FRRDRC,RHO,RHO13)
598            IF (DOBCK) CALL V1BCK(FR0BCK,FZ0BCK,FRRBCK,FRZBCK,
599     &                            FZZBCK,RHO,RHOGRD)
600            IF (DOVWN) THEN
601               IF (.NOT.TRPLET) THEN
602                  CALL V1VWN(VVWN0,FRRVWN,RHO,RHO13)
603               ELSE
604                  CALL VTVWN(FRRVWN,RHO,RHO13)
605                  FRRVWN = DP5*FRRVWN
606               END IF
607            END IF
608            IF (DOLYP) THEN
609               RHOA = DP5*RHO
610               RHGA = (DP5*RHOGRD)**2
611C              CALL V1LYP(FR0LYP,FZ0LYP,FRRLYP,FRZLYP,
612C    &                    FZZLYP,RHO,RHOGRD)
613               CALL GLYPCO(DF1000,DF0100,DF0010,DF0001,
614     &                     DF00001,RHO,RHO13,RHOGRD,.TRUE.)
615               CALL VTLYP (DF2000,DF0200,DF1100,DF1010,
616     &                     DF0101,DF1001,DF0110,DF10001,
617     &                     DF01001,RHOA,RHOA,RHGA,RHGA,RHGA)
618               IF (.NOT.TRPLET) THEN
619                  FZ0LYP = DP5*(DF0010 + DF00001)*RHOGRD
620                  FRRLYP = DP5*(DF2000 + DF1100)
621                  FRZLYP = DP5*(DF1010 + DF1001+DF10001)*RHOGRD
622                  FZZLYP = FZ0LYP/RHOGRD
623               ELSE
624                  FZ0LYP = DP5*(DF0010 - DF00001)*RHOGRD
625                  FRRLYP = DP5*(DF2000 - DF1100)
626                  FRZLYP = DP5*(DF1010 - DF1001)*RHOGRD
627                  FZZLYP = FZ0LYP/RHOGRD
628               END IF
629            ELSE
630               FZ0LYP = 0D0
631               FRRLYP = 0D0
632               FRZLYP = 0D0
633               FZZLYP = 0D0
634            END IF
635C
636            ZNV = D1/RHOGRD
637            FZ0 = ZNV*(WBCK*FZ0BCK + WLYP*FZ0LYP)
638            FRR = WDRC*FRRDRC + WVWN*FRRVWN
639     &          + WBCK*FRRBCK + WLYP*FRRLYP
640            FRZ = WBCK*FRZBCK + WLYP*FRZLYP
641            FZZ = WBCK*FZZBCK + WLYP*FZZLYP
642C
643            RX = ZNV*RHG(1)
644            RY = ZNV*RHG(2)
645            RZ = ZNV*RHG(3)
646C
647C           Linear transformation
648C
649            IF (.NOT.DOHES) THEN
650               BR = B3(1)*RX + B3(2)*RY + B3(3)*RZ
651               FAC0 = FRR*B0 + FRZ*BR
652               FACR = FRZ*B0 + FZZ*BR
653               IF (NSYM.EQ.1) THEN
654                  DO I = 1, NBAST
655                     G0 = GAO(I)
656                     GX = GAO1(I,1)
657                     GY = GAO1(I,2)
658                     GZ = GAO1(I,3)
659                     DO J = 1, I
660                        A0 = G0*GAO(J)
661                        AX = GX*GAO(J) + G0*GAO1(J,1)
662                        AY = GY*GAO(J) + G0*GAO1(J,2)
663                        AZ = GZ*GAO(J) + G0*GAO1(J,3)
664                        AR = AX*RX + AY*RY + AZ*RZ
665                        AB = AX*B3(1) + AY*B3(2) + AZ*B3(3) - AR*BR
666                        EXCMAT(J,I) = EXCMAT(J,I)+FAC0*A0+FACR*AR+FZ0*AB
667                     END DO
668                  END DO
669               ELSE
670                  DO ISYM = 1, NSYM
671                     ISTR = IBAS(ISYM) + 1
672                     IEND = IBAS(ISYM) + NBAS(ISYM)
673                     JSYM = MULD2H(ISYM,KSYMOP)
674                     IF (ISYM.GE.JSYM) THEN
675                        JSTR = IBAS(JSYM) + 1
676                        JEND = IBAS(JSYM) + NBAS(JSYM)
677                        DO I = ISTR, IEND
678                           G0 = GAO(I)
679                           GX = GAO1(I,1)
680                           GY = GAO1(I,2)
681                           GZ = GAO1(I,3)
682                           DO J = JSTR, MIN(I,JEND)
683                              A0 = G0*GAO(J)
684                              AX = GX*GAO(J) + G0*GAO1(J,1)
685                              AY = GY*GAO(J) + G0*GAO1(J,2)
686                              AZ = GZ*GAO(J) + G0*GAO1(J,3)
687                              AR = AX*RX + AY*RY + AZ*RZ
688                              AB = AX*B3(1) + AY*B3(2) + AZ*B3(3) -AR*BR
689                              EXCMAT(J,I) = EXCMAT(J,I) + FAC0*A0
690     &                                    + FACR*AR + FZ0*AB
691                           END DO
692                        END DO
693                     END IF
694                  END DO
695               END IF
696C
697C           Explicit Hessian
698C
699            ELSE
700               CALL DGEMM('T','N',NORBT,1,NBAST,1.D0,
701     &                    DTRMAT,NBAST,
702     &                    GAO,NBAST,0.D0,
703     &                    C0,NORBT)
704               CALL DGEMM('T','N',NORBT,3,NBAST,1.D0,
705     &                    DTRMAT,NBAST,
706     &                    GAO1,NBAST,0.D0,
707     &                    C1,NORBT)
708C
709               DO I = 1, NORBT
710                 C2(I) = RX*C1(I,1)+RY*C1(I,2)+RZ*C1(I,3)
711               END DO
712C
713               DO B = NOCCT + 1, NORBT
714               DO J = 1, NOCCT
715                 IB  = B - NOCCT
716                 GBJ = C0(B)*C0(J)
717                 PBJ = C2(B)*C0(J)   + C2(J)*C0(B)
718                 CBX = C0(B)*C1(J,1) + C1(B,1)*C0(J)
719                 CBY = C0(B)*C1(J,2) + C1(B,2)*C0(J)
720                 CBZ = C0(B)*C1(J,3) + C1(B,3)*C0(J)
721                 DO A = NOCCT + 1, B
722                 DO I = 1, NOCCT
723                   IA  = A - NOCCT
724                   GAI = C0(A)*C0(I)
725                   PAI = C2(A)*C0(I) + C2(I)*C0(A)
726                   CAB = CBX*(C0(A)*C1(I,1) + C1(A,1)*C0(I))
727     &                 + CBY*(C0(A)*C1(I,2) + C1(A,2)*C0(I))
728     &                 + CBZ*(C0(A)*C1(I,3) + C1(A,3)*C0(I))
729                   HES(I,IA,J,IB) = HES(I,IA,J,IB)
730     &                            + FRR*GAI*GBJ
731     &                            + FRZ*(PAI*GBJ + GAI*PBJ)
732     &                            + FZZ*PAI*PBJ
733     &                            + FZ0*(CAB - PAI*PBJ)
734                  END DO
735                  END DO
736               END DO
737               END DO
738            END IF
739         END IF
740      END IF
741      RETURN
742      END
743