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