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!Revision 1.2 2000/05/24 19:04:06 hjj 21!new getref calls with appropriate NDREF instead of NCREF 22!(fixing error for triplet with CSF) 23!=========================================================================== 24 25 SUBROUTINE MELONE(PRPMO,IOPSYM,DEN1,OVLAP,ONETOT,IPRONE,PRPLBL) 26C 27C CALCULATE AVERAGE VALUE OF ONE ELECTRON OPERATOR 28C 29#include "implicit.h" 30C 31 CHARACTER*(*)PRPLBL 32C 33 DIMENSION PRPMO(NORBT,*),DEN1(NASHDI,*) 34C 35#include "priunit.h" 36#include "wrkrsp.h" 37#include "infrsp.h" 38#include "inforb.h" 39#include "infdim.h" 40#include "infpri.h" 41C 42 PARAMETER ( D2 = 2.0D0 , D0 = 0.0D0 ) 43 PARAMETER ( BIGLIM = 100000.0D0, SMLLIM = 0.01D0 ) 44C 45 ONEACT = D0 46 ONEINA = D0 47 DO 50 ISYM = 1,NSYM 48 JSYM = MULD2H(IOPSYM,ISYM) 49 NASHI = NASH(ISYM) 50 NISHI = NISH(ISYM) 51 IORBI = IORB(ISYM) 52 IASHI = IASH(ISYM) 53 NASHJ = NASH(JSYM) 54 NISHJ = NISH(JSYM) 55 IORBJ = IORB(JSYM) 56 IASHJ = IASH(JSYM) 57 DO 80 IINAC = 1,NISHI 58 ONEINA = ONEINA + PRPMO(IORBI+IINAC,IORBI+IINAC) 59 80 CONTINUE 60 DO 60 IA = 1,NASHI 61 DO 70 JA = 1,NASHJ 62 ONEACT = ONEACT + DEN1(IASHI+IA,IASHJ+JA) * 63 * PRPMO(IORBI+NISHI+IA,IORBJ+NISHJ+JA) 64 70 CONTINUE 65 60 CONTINUE 66 50 CONTINUE 67C 68 ONEINA = ONEINA * D2 * OVLAP 69 ONETOT = ONEINA + ONEACT 70 IF (IPRRSP.GE.IPRONE) THEN 71 IF (ABS(ONETOT) .GT. SMLLIM .AND. ABS(ONETOT) .LT. BIGLIM) 72 * THEN 73 WRITE(LUPRI,'(3(/5X,2A,F15.8))') 74 * PRPLBL,' INACTIVE PART:',ONEINA, 75 * PRPLBL,' ACTIVE PART :',ONEACT, 76 * PRPLBL,' TOTAL :',ONETOT 77 ELSE 78 WRITE(LUPRI,'(3(/5X,2A,1P,D15.8))') 79 * PRPLBL,' INACTIVE PART:',ONEINA, 80 * PRPLBL,' ACTIVE PART :',ONEACT, 81 * PRPLBL,' TOTAL :',ONETOT 82 END IF 83 END IF 84C 85C END OF MELONE 86C 87 RETURN 88 END 89 SUBROUTINE MELTWO(H1,DEN1,DEN2,OVLAP,ISYMDN,IDAX,IOPSYM, 90 * ISPIN1,ISPIN2,ISPIN3,IKLVL, 91 * ZYMAT1,ISYM1,ZYMAT2,ISYM2,ZYMAT3,ISYM3, 92 * RES,CMO,WRK,LWRK) 93C 94C Calculation of the two-electron matrix element <L|K|R>. 95C The generalized hamiltonian K can be one-index transformed 0, 1 or 2 96C times according to IKLVL. (0 gives the ordinary untransformed hamiltonian.) 97C 98#include "implicit.h" 99C 100#include "priunit.h" 101#include "wrkrsp.h" 102#include "infrsp.h" 103#include "inforb.h" 104#include "infdim.h" 105#include "infpri.h" 106#include "maxash.h" 107#include "maxorb.h" 108#include "infind.h" 109C 110 PARAMETER (D1=1.0D0, D2=2.0D0) 111 DIMENSION WRK(*) 112 DIMENSION H1(NORBT,NORBT) 113 DIMENSION DEN1(NASHDI,NASHDI),DEN2(N2ASHX,N2ASHX) 114 DIMENSION ZYMAT1(NORBT,NORBT),ZYMAT2(NORBT,NORBT) 115 DIMENSION ZYMAT3(NORBT,NORBT) 116 DIMENSION CMO(*) 117C 118 LOGICAL LCON,LORB 119C 120 CALL QENTER('MELTWO') 121C 122 IGRSPI = 0 123C 124 IF (IPRRSP.GT.200) THEN 125 WRITE(LUPRI,'(/A)') 'ENTER MELTWO' 126 WRITE(LUPRI,'(A,I5)') 'IDAX =',IDAX 127 WRITE(LUPRI,'(A,I5)') 'IOPSYM =',IOPSYM 128 WRITE(LUPRI,'(A,I5)') 'ISYMDN =',ISYMDN 129 WRITE(LUPRI,*) 'OVLAP =',OVLAP 130 WRITE(LUPRI,'(A,I5)') 'IGRSPI =',IGRSPI 131 WRITE(LUPRI,'(A,I5)') 'ISPIN1 =',ISPIN1 132 WRITE(LUPRI,'(A,I5)') 'ISPIN2 =',ISPIN2 133 WRITE(LUPRI,'(A,I5)') 'ISPIN3 =',ISPIN3 134 WRITE(LUPRI,'(A,I5)') 'IKLVL =',IKLVL 135 WRITE(LUPRI,'(A,I5)') 'ISYM1 =',ISYM1 136 WRITE(LUPRI,'(A,I5)') 'ISYM2 =',ISYM2 137 WRITE(LUPRI,'(A,I5)') 'ISYM3 =',ISYM3 138 WRITE(LUPRI,'(A)') 'One-electron H1' 139 CALL OUTPUT(H1,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 140 WRITE(LUPRI,'(A)') 'One-electron density' 141 CALL OUTPUT(DEN1,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 142 WRITE(LUPRI,'(A)') 'Two-electron density' 143 CALL PRIAC2(DEN2,NASHT,LUPRI) 144 END IF 145C 146 KH2AX = 1 147 LH2AX = N2ASHX * N2ASHX 148 IF (DIROIT) THEN 149 IF (IKLVL.EQ.2) LH2AX = LH2AX * 2 150 IF (IKLVL.EQ.3) LH2AX = LH2AX * 4 151 END IF 152 KFA = KH2AX + LH2AX 153 KFI = KFA + N2ORBX 154 KQA = KFI + N2ORBX 155 KQB = KQA + NORBT * NASHDI 156 KPVD = KQB + NORBT * NASHDI 157 KH2 = KPVD + N2ASHX * N2ASHX 158 KH2X = KH2 + N2ORBX 159 KWRK1 = KH2X + N2ORBX 160 LWRK1 = LWRK - KWRK1 161 IF (LWRK1.LT.0) CALL ERRWRK('MELTWO',KWRK1-1,LWRK) 162C 163 CALL DCOPY(N2ORBX,H1,1,WRK(KFI),1) 164 CALL DZERO(WRK(KFA),N2ORBX) 165 IF (NASHT .GT. 0) THEN 166 CALL DZERO(WRK(KQA),NASHT*NORBT) 167 CALL DZERO(WRK(KQB),NASHT*NORBT) 168 CALL DZERO(WRK(KH2AX),LH2AX) 169 END IF 170C 171 IF (DIROIT) THEN 172 INTSYM = 1 173 ELSE 174 INTSYM = IOPSYM 175 END IF 176C 177C We only need the transformed inactive Fock matrix and the 178C tranformed active two-electron integrals. 179C However, FA, QA, and QB must also be initialized because 180C they are referenced in RSPFXD and this has given 181C floating point exception on DEC alpha. /July 2001 HJAaJ 182C 183 LCON = .TRUE. 184 LORB = .FALSE. 185 KFREE = 1 186 LFREE = LWRK1 187 IF (DIROIT) THEN 188 CALL RSPFXD(WRK(KFI),WRK(KFA),WRK(KQA),WRK(KQB),WRK(KH2AX), 189 * IDAX,INTSYM,ISYMDN,DEN1,DEN2,WRK(KPVD), 190 * WRK(KH2),WRK(KH2X),WRK(KWRK1),KFREE,LFREE, 191 * LCON,LORB,IPRRSP,IGRSPI,ISPIN1,ISPIN2,ISPIN3, 192 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2,ZYMAT3,ISYM3) 193 ELSE 194 CALL RSPFX(WRK(KFI),WRK(KFA),WRK(KQA),WRK(KQB),WRK(KH2AX),IDAX, 195 * INTSYM,ISYMDN,DEN1,DEN2,WRK(KPVD),WRK(KH2X), 196 * WRK(KWRK1),KFREE,LFREE,LCON,LORB,IPRRSP) 197 END IF 198C 199C The transformed two-electron integrals do not have the factor 1/2 200C that MEL2 assumes. If IKLVL = 2 we have two densities. For the 201C spin free case they are equal and we just add them. 202C 203 IF (DIROIT.AND.(IKLVL.GT.0)) THEN 204 IF (IKLVL.EQ.2) THEN 205 CALL DAXPY(N2ASHX*N2ASHX,D1,WRK(KH2AX+N2ASHX*N2ASHX),1, 206 * WRK(KH2AX),1) 207 END IF 208 CALL DSCAL(N2ASHX*N2ASHX,D2,WRK(KH2AX),1) 209 END IF 210C 211 CALL MEL2(H1,WRK(KFI),WRK(KH2AX),IOPSYM,DEN1,DEN2,OVLAP,RES) 212C 213 CALL QEXIT('MELTWO') 214C 215 RETURN 216 END 217 SUBROUTINE MEL2(H1,FI,H2A,IOPSYM,DEN1,DEN2,OVLAP,RES) 218C 219C Calculate the two-electron matrix element <L|H|R>. 220C 221C RES = OVLAP * sum(i) [FI(i,i) + H1(i,i)] + sum(x,y) FI(x,y) DEN1(x,y) 222C + 1/2 * sum(x,y,u,v) H2A(x,y,u,v) DEN2(x,y,u,v) 223C 224#include "implicit.h" 225C 226#include "priunit.h" 227#include "wrkrsp.h" 228#include "infrsp.h" 229#include "inforb.h" 230#include "infdim.h" 231#include "infpri.h" 232#include "maxash.h" 233#include "maxorb.h" 234#include "infind.h" 235C 236 PARAMETER ( D0 = 0.0D0, DH = 0.5D0 ) 237C 238 DIMENSION H1(NORBT,NORBT),FI(NORBT,NORBT) 239 DIMENSION H2A(N2ASHX*N2ASHX) 240 DIMENSION DEN1(NASHDI,NASHDI),DEN2(N2ASHX*N2ASHX) 241C 242 IF (IPRRSP.GT.200) THEN 243 WRITE(LUPRI,'(/A)') 'ENTER MEL2' 244 WRITE(LUPRI,'(A,I5)') 'IOPSYM =',IOPSYM 245 WRITE(LUPRI,*) 'OVLAP =',OVLAP 246 WRITE(LUPRI,'(A)') 'One-electron H1' 247 CALL OUTPUT(H1,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 248 WRITE(LUPRI,'(A)') 'Inactive Fock matrix FI' 249 CALL OUTPUT(FI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 250 WRITE(LUPRI,'(A)') 'One-electron density' 251 CALL OUTPUT(DEN1,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 252 WRITE(LUPRI,'(A)') 'Two-electron density' 253 CALL PRIAC2(DEN2,NASHT,LUPRI) 254 WRITE(LUPRI,'(A)') 'Two-electron H2A' 255 CALL PRIAC2(H2A,NASHT,LUPRI) 256 IF (DIROIT) THEN 257 CALL PRIAC2(H2A(N2ASHX*N2ASHX+1),NASHT,LUPRI) 258 END IF 259 END IF 260C 261 RES = D0 262 DO 10 ISYM = 1,NSYM 263 NISHI = NISH(ISYM) 264 IORBI = IORB(ISYM) 265 DO 20 IINAC = 1,NISHI 266 RES = RES + FI(IORBI+IINAC,IORBI+IINAC) 267 RES = RES + H1(IORBI+IINAC,IORBI+IINAC) 268 20 CONTINUE 269 10 CONTINUE 270 RES = RES * OVLAP 271C 272 DO 30 ISYM = 1,NSYM 273 JSYM = MULD2H(IOPSYM,ISYM) 274 NASHI = NASH(ISYM) 275 NISHI = NISH(ISYM) 276 IORBI = IORB(ISYM) 277 IASHI = IASH(ISYM) 278 NASHJ = NASH(JSYM) 279 NISHJ = NISH(JSYM) 280 IORBJ = IORB(JSYM) 281 IASHJ = IASH(JSYM) 282 DO 40 IAC = 1,NASHI 283 DO 50 JAC = 1,NASHJ 284 RES = RES + FI(IORBI+NISHI+IAC,IORBJ+NISHJ+JAC)* 285 * DEN1(IASHI+IAC,IASHJ+JAC) 286 50 CONTINUE 287 40 CONTINUE 288 30 CONTINUE 289C 290 RES = RES + DH*DDOT(N2ASHX*N2ASHX,H2A,1,DEN2,1) 291C 292 IF (IPRRSP.GT.10) THEN 293 WRITE(LUPRI,'(/A,F15.8)') ' Result in MEL2:',RES 294 END IF 295C 296 RETURN 297 END 298 SUBROUTINE S4DRV(KZYVR,KZYV1,KZYV2,KZYV3, 299 * IGRSYM,ISYMV1,ISYMV2,ISYMV3, 300 * S4TRS,VEC1,VEC2,VEC3, 301 * UDV,ZYMAT,DEN1,XINDX,MJWOP,WRK,LWRK) 302C 303C Drive the computation of S[4] times three vectors 304C 305#include "implicit.h" 306#include "infdim.h" 307#include "inforb.h" 308#include "maxorb.h" 309#include "maxash.h" 310#include "infrsp.h" 311#include "wrkrsp.h" 312#include "rspprp.h" 313#include "infhyp.h" 314#include "infvar.h" 315#include "infind.h" 316#include "infpri.h" 317#include "qrinf.h" 318C 319 PARAMETER( D1= 1.0D0, D6=6.0D0, DH=0.5D0 ) 320C 321 DIMENSION WRK(*) 322 DIMENSION S4TRS(KZYVR), MJWOP(2,MAXWOP,8) 323 DIMENSION VEC1(KZYV1), VEC2(KZYV2), VEC3(KZYV3) 324 DIMENSION ZYMAT(NORBT,NORBT) 325 DIMENSION XINDX(*) 326 DIMENSION UDV(NASHDI,NASHDI) 327 DIMENSION DEN1(NASHDI,NASHDI) 328C 329 LOGICAL LORB, LCON, LREF, TDM, NORHO2 330C 331C Layout some workspace 332C 333 KCREF = 1 334 KRES = KCREF + MZCONF(1) 335 KFREE = KRES + KZYVR 336 LFREE = LWRK - KFREE + 1 337 IF (LFREE.LT.0) CALL ERRWRK('S4DRV 1',KFREE-1,LWRK) 338C 339C 340C Initialize variables 341C 342 TDM = .TRUE. 343 NORHO2 = .TRUE. 344 NSIM = 1 345 ISPIN = 0 346C 347 CALL GETREF(WRK(KCREF),MZCONF(1)) 348C 349 CALL SCASE1(KZYVR,KZYV1,KZYV2,KZYV3, 350 * IGRSYM,ISYMV1,ISYMV2,ISYMV3, 351 * S4TRS,VEC1,VEC2,VEC3,WRK(KCREF), 352 * UDV,DEN1,XINDX,MJWOP,WRK(KFREE),LFREE) 353C 354 CALL SCASE1(KZYVR,KZYV2,KZYV1,KZYV3, 355 * IGRSYM,ISYMV2,ISYMV1,ISYMV3, 356 * S4TRS,VEC2,VEC1,VEC3,WRK(KCREF), 357 * UDV,DEN1,XINDX,MJWOP,WRK(KFREE),LFREE) 358C 359 CALL SCASE1(KZYVR,KZYV3,KZYV1,KZYV2, 360 * IGRSYM,ISYMV3,ISYMV1,ISYMV2, 361 * S4TRS,VEC3,VEC1,VEC2,WRK(KCREF), 362 * UDV,DEN1,XINDX,MJWOP,WRK(KFREE),LFREE) 363C 364 CALL SCASE2(KZYVR,KZYV1,KZYV2,KZYV3, 365 * IGRSYM,ISYMV1,ISYMV2,ISYMV3, 366 * S4TRS,VEC1,VEC2,VEC3,WRK(KCREF), 367 * UDV,ZYMAT,DEN1,XINDX,MJWOP,WRK(KFREE),LFREE) 368C 369 CALL SCASE2(KZYVR,KZYV1,KZYV3,KZYV2, 370 * IGRSYM,ISYMV1,ISYMV3,ISYMV2, 371 * S4TRS,VEC1,VEC3,VEC2,WRK(KCREF), 372 * UDV,ZYMAT,DEN1,XINDX,MJWOP,WRK(KFREE),LFREE) 373C 374 IF (MZWOPT(ISYMV1).EQ.0 .OR. 375 * MZCONF(ISYMV2).EQ.0 .OR. 376 * MZCONF(ISYMV3).EQ.0) GOTO 4000 377C 378C / <02L| [qj,K(k1)] |03R> + <03L| [qj,K(k1)] |02R> \ 379C | 0 | 380C | <02L| [qj+,K(k1)] |03R> + <03L| [qj+,K(k1)] |02R> | 381C \ 0 / 382C 383C Construct <03L|..|02R> + <02L|..|03R> density 384C 385 ILSYM = MULD2H(IREFSY,ISYMV2) 386 IRSYM = MULD2H(IREFSY,ISYMV3) 387 NCL = MZCONF(ISYMV2) 388 NCR = MZCONF(ISYMV3) 389 KZVARL = MZYVAR(ISYMV2) 390 KZVARR = MZYVAR(ISYMV3) 391 LREF = .FALSE. 392 ISYMDN = MULD2H(ILSYM,IRSYM) 393 CALL DZERO(DEN1,NASHT*NASHT) 394 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 395 * VEC2,VEC3,OVLAP,DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2, 396 * XINDX,WRK,KFREE,LFREE,LREF) 397C 398C Make the gradient 399C 400 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,ZYMAT,MJWOP) 401C 402 IF ( MZWOPT(IGRSYM) .GT. 0 ) THEN 403 CALL ORBSX(NSIM,IGRSYM,KZYVR,S4TRS,ZYMAT,OVLAP,ISYMDN, 404 * DEN1,MJWOP,WRK(KFREE),LFREE) 405 END IF 406 IF (ISYMV2.NE.ISYMV3) GOTO 4000 407C 408C / <0| [qj ,K(k1)] |0> \ 409C | 1/2<j| K(k1) |0> | * ( S(2)S(3)' + S(2)'S(3) ) 410C | <0| [qj+,K(k1)] |0> | 411C \ -1/2<0| K(k1) |j> / 412C 413 NZCONF = MZCONF(ISYMV2) 414 NZVAR = MZVAR(ISYMV2) 415 F1 = DDOT(NZCONF,VEC2,1,VEC3(NZVAR+1),1) + 416 * DDOT(NZCONF,VEC3,1,VEC2(NZVAR+1),1) 417C 418 ISYMDN = 1 419 OVLAP = D1 420 ISYMST = MULD2H(IGRSYM,IREFSY) 421 IF(ISYMST .EQ. IREFSY ) THEN 422 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 423 ELSE 424 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 425 END IF 426 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 427 LREF = .TRUE. 428 IF (LCON) THEN 429 CALL DZERO(WRK(KRES),KZYVR) 430 CALL CONSX(NSIM,KZYVR,IGRSYM,ZYMAT,WRK(KCREF), 431 * MZCONF(1),MZCONF(1),IREFSY,MZCONF(IGRSYM),ISYMST, 432 * LREF,WRK(KRES),XINDX,WRK(KFREE),LFREE) 433 CALL DSCAL(KZYVR,DH,WRK(KRES),1) 434 CALL DAXPY(KZYVR,F1,WRK(KRES),1,S4TRS,1) 435 END IF 436 IF (LORB) THEN 437 CALL DZERO(WRK(KRES),KZYVR) 438 CALL ORBSX(NSIM,IGRSYM,KZYVR,WRK(KRES),ZYMAT,OVLAP,ISYMDN, 439 * UDV,MJWOP,WRK(KFREE),LFREE) 440 CALL DAXPY(KZYVR,F1,WRK(KRES),1,S4TRS,1) 441 END IF 442C 443 4000 CONTINUE 444C 445 CALL SCASE3(KZYVR,KZYV1,KZYV2,KZYV3, 446 * IGRSYM,ISYMV1,ISYMV2,ISYMV3, 447 * S4TRS,VEC1,VEC2,VEC3,WRK(KCREF), 448 * UDV,ZYMAT,DEN1,XINDX,MJWOP,WRK(KFREE),LFREE) 449C 450 CALL SCASE3(KZYVR,KZYV1,KZYV3,KZYV2, 451 * IGRSYM,ISYMV1,ISYMV3,ISYMV2, 452 * S4TRS,VEC1,VEC3,VEC2,WRK(KCREF), 453 * UDV,ZYMAT,DEN1,XINDX,MJWOP,WRK(KFREE),LFREE) 454C 455 IF (MZWOPT(ISYMV1).EQ.0 .OR. 456 * MZWOPT(ISYMV2).EQ.0 .OR. 457 * MZWOPT(ISYMV3).EQ.0) RETURN 458C 459C / <0| [qj ,K(k1,k2,k3)+K(k1,k3,k2)] |0> \ 460C | <j| K(k1,k2,k3)+K(k1,k3,k2) |0> | * 1/6 461C | <0| [qj+,K(k1,k2,k3)+K(k1,k3,k2)] |0> | 462C \ -<0| K(k1,k2,k3)+K(k1,k3,k2) |j> / 463C 464 ISYMDN = 1 465 OVLAP = D1 466C 467C Make the gradient 468C 469 ISYMV = IREFSY 470 ISYMST = MULD2H(IGRSYM,IREFSY) 471 IF ( ISYMST .EQ. IREFSY ) THEN 472 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 473 ELSE 474 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 475 END IF 476 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 477 LREF = .TRUE. 478 NZYVEC = MZCONF(1) 479 NZCVEC = MZCONF(1) 480 CALL TRZYM2(VEC1,VEC2,VEC3,KZYV1,KZYV2,KZYV3, 481 * ISYMV1,ISYMV2,ISYMV3,ZYMAT,MJWOP, 482 * WRK(KFREE),LFREE) 483 CALL DSCAL(NORBT*NORBT,1/D6,ZYMAT,1) 484 CALL RSP1GR(NSIM,KZYVR,IDUM,ISPIN,IGRSYM,ISPIN,ISYMV,S4TRS, 485 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,ZYMAT, 486 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 487C 488 CALL TRZYM2(VEC1,VEC3,VEC2,KZYV1,KZYV3,KZYV2, 489 * ISYMV1,ISYMV3,ISYMV2,ZYMAT,MJWOP, 490 * WRK(KFREE),LFREE) 491 CALL DSCAL(NORBT*NORBT,1/D6,ZYMAT,1) 492 CALL RSP1GR(NSIM,KZYVR,IDUM,ISPIN,IGRSYM,ISPIN,ISYMV,S4TRS, 493 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,ZYMAT, 494 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 495C 496 RETURN 497 END 498 SUBROUTINE SCASE1(KZYVR,KZYV1,KZYV2,KZYV3, 499 * IGRSYM,ISYMV1,ISYMV2,ISYMV3, 500 * S4TRS,VEC1,VEC2,VEC3,CREF, 501 * UDV,DEN1,XINDX,MJWOP,WRK,LWRK) 502C 503C 504C / <01L| qj |0> + <0| qj |01R> \ 505C | 0 | * -2/3*(S(2)S(3)' + S(2)'S(3)) 506C | <01L| qj+|0> + <0| qj+|01R> | 507C \ 0 / 508C 509C 510#include "implicit.h" 511#include "infdim.h" 512#include "inforb.h" 513#include "maxorb.h" 514#include "maxash.h" 515#include "infrsp.h" 516#include "wrkrsp.h" 517#include "rspprp.h" 518#include "infhyp.h" 519#include "infvar.h" 520#include "infind.h" 521#include "infpri.h" 522#include "qrinf.h" 523C 524 PARAMETER( D2=2.0D0, D3=3.0D0 ) 525C 526 DIMENSION WRK(*) 527 DIMENSION S4TRS(KZYVR) 528 DIMENSION VEC1(KZYV1), VEC2(KZYV2), VEC3(KZYV3) 529 DIMENSION XINDX(*) 530 DIMENSION UDV(NASHDI,NASHDI) 531 DIMENSION DEN1(NASHDI,NASHDI) 532 DIMENSION CREF(*), MJWOP(2,MAXWOP,8) 533C 534 LOGICAL LORB, LCON, LREF, TDM, NORHO2 535C 536 IF (MZCONF(ISYMV1).EQ.0 .OR. 537 * MZCONF(ISYMV2).EQ.0 .OR. 538 * MZCONF(ISYMV3).EQ.0) RETURN 539C 540 IF (ISYMV2.NE.ISYMV3) RETURN 541C 542C Initialize variables 543C 544 TDM = .TRUE. 545 NORHO2 = .TRUE. 546 NSIM = 1 547 ISPIN = 0 548 KONE = 1 549C 550C Construct the density matrix <01L|..|0> + <0|..|01R> 551C 552 ILSYM = IREFSY 553 IRSYM = MULD2H(IREFSY,ISYMV1) 554 NCL = MZCONF(1) 555 NCR = MZCONF(ISYMV1) 556 KZVARL = MZCONF(1) 557 KZVARR = MZYVAR(ISYMV1) 558 LREF = .TRUE. 559 ISYMDN = MULD2H(ILSYM,IRSYM) 560 CALL DZERO(DEN1,NASHT*NASHT) 561 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 562 * CREF,VEC1,OVLAP,DEN1,DUMMY,ISPIN,ISPIN,TDM, 563 * NORHO2,XINDX,WRK,KONE,LWRK,LREF) 564C 565C Put the factor into the density matrix 566C 567 NZCONF = MZCONF(ISYMV2) 568 NZVAR = MZVAR(ISYMV2) 569 F1 = DDOT(NZCONF,VEC2,1,VEC3(NZVAR+1),1) + 570 * DDOT(NZCONF,VEC3,1,VEC2(NZVAR+1),1) 571 F1 = -F1*D2/D3 572 CALL DSCAL(NASHT*NASHT,F1,DEN1,1) 573C 574C Make the gradient on the fly: 575C 576 NZCONF = MZCONF(IGRSYM) 577 NYCONF = MZCONF(IGRSYM) + MZVAR(IGRSYM) 578C 579 DO 500 IC = 1, MZWOPT(IGRSYM) 580 K = MJWOP(1,IC,IGRSYM) 581 L = MJWOP(2,IC,IGRSYM) 582 ITYPK = IOBTYP(K) 583 ITYPL = IOBTYP(L) 584 IF(ITYPK.EQ.JTACT .AND. ITYPL. EQ. JTACT) THEN 585 NWK = ISW(K) - NISHT 586 NWL = ISW(L) - NISHT 587 S4TRS(NZCONF+IC) = S4TRS(NZCONF+IC) 588 * + DEN1(NWK,NWL) 589 S4TRS(NYCONF+IC) = S4TRS(NYCONF+IC) 590 * + DEN1(NWL,NWK) 591 END IF 592500 CONTINUE 593C 594 RETURN 595 END 596 SUBROUTINE SCASE2(KZYVR,KZYV1,KZYV2,KZYV3, 597 * IGRSYM,ISYMV1,ISYMV2,ISYMV3, 598 * S4TRS,VEC1,VEC2,VEC3,CREF, 599 * UDV,ZYMAT,DEN1,XINDX,MJWOP,WRK,LWRK) 600C 601C 602#include "implicit.h" 603#include "infdim.h" 604#include "inforb.h" 605#include "maxorb.h" 606#include "maxash.h" 607#include "infrsp.h" 608#include "wrkrsp.h" 609#include "rspprp.h" 610#include "infhyp.h" 611#include "infvar.h" 612#include "infind.h" 613#include "infpri.h" 614#include "qrinf.h" 615C 616 PARAMETER( DH=0.5D0, D0=0.0D0, D1=1.0D0, D6=6.0D0 ) 617C 618 DIMENSION WRK(*) 619 DIMENSION S4TRS(KZYVR), MJWOP(2,MAXWOP,8) 620 DIMENSION VEC1(KZYV1), VEC2(KZYV2), VEC3(KZYV3) 621 DIMENSION XINDX(*) 622 DIMENSION ZYMAT(NORBT,NORBT) 623 DIMENSION UDV(NASHDI,NASHDI) 624 DIMENSION DEN1(NASHDI,NASHDI) 625 DIMENSION CREF(*) 626C 627 LOGICAL LORB, LCON, LREF, TDM, NORHO2 628C 629C Initialize variables 630C 631 TDM = .TRUE. 632 NORHO2 = .TRUE. 633 NSIM = 1 634 ISPIN = 0 635 IPRONE = 75 636 F = D0 637 F1 = D0 638 F2 = D0 639 KONE = 1 640C 641 IF (MZCONF(ISYMV1).EQ.0 .OR. 642 * MZCONF(ISYMV2).EQ.0 .OR. 643 * MZCONF(ISYMV3).EQ.0) GOTO 1000 644C 645C / 0 \ 646C | Sj(1) | * 1/6*S(2)S(3)' 647C | 0 | 648C \ -Sj(1)' / 649C 650 IF (ISYMV2.EQ.ISYMV3) THEN 651 NZCONF = MZCONF(ISYMV2) 652 NZVAR = MZVAR(ISYMV2) 653 FACT = 1/D6*DDOT(NZCONF,VEC2,1,VEC3(NZVAR+1),1) 654 NZCONF = MZCONF(ISYMV1) 655 NZVAR = MZVAR(ISYMV1) 656 CALL DAXPY(NZCONF,FACT,VEC1,1,S4TRS,1) 657 CALL DAXPY(NZCONF,-FACT,VEC1(NZVAR+1),1,S4TRS(NZVAR+1),1) 658 END IF 659C 660 IF (ISYMV1.NE.ISYMV3) RETURN 661C 662C / 0 \ 663C | (F+F1)*Sj(2) | 6*F1 = S(3)'S(1) - 2*S(1)'S(3) 664C | 0 | 665C \ (F+F2)*Sj(2)' / 6*F2 = 2*S(3)'S(1) - S(1)'S(3) 666C 667C F = 1/2*<03L|K(k1)|0> + <0|K(k1)|03R> + 1/2*<0|K(k1,k3)|0> 668C 669 NZCONF = MZCONF(ISYMV1) 670 NZVAR = MZVAR(ISYMV1) 671 F1 = DDOT(NZCONF,VEC1,1,VEC3(NZVAR+1),1) - 672 * 2*DDOT(NZCONF,VEC3,1,VEC1(NZVAR+1),1) 673 F2 = 2*DDOT(NZCONF,VEC1,1,VEC3(NZVAR+1),1) - 674 * DDOT(NZCONF,VEC3,1,VEC1(NZVAR+1),1) 675 F1 = F1/D6 676 F2 = F2/D6 677C 678 1000 CONTINUE 679 IF (MZCONF(ISYMV3).EQ.0 .OR. IGRSYM.NE.ISYMV2) RETURN 680C 681 IF (MZWOPT(ISYMV1).GT.0 .AND. MZCONF(ISYMV3).GT.0) THEN 682C 683 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,ZYMAT,MJWOP) 684C 685C F3R = <0|K(k1)|-03R> 686C 687 ILSYM = IREFSY 688 IRSYM = MULD2H(IREFSY,ISYMV3) 689 NCL = MZCONF(1) 690 NCR = MZCONF(ISYMV3) 691 IOFF = 1 692 CALL DZERO(DEN1,NASHT*NASHT) 693 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,CREF,VEC3(IOFF), 694 * DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2,XINDX,WRK, 695 * KONE,LWRK) 696 OVLAP = D0 697 IF (ILSYM.EQ.IRSYM) 698 * OVLAP = DDOT(NCL,CREF,1,VEC3(IOFF),1) 699C 700 CALL MELONE(ZYMAT,ISYMV1,DEN1,OVLAP,F3R,IPRONE,'F3R in SCASE2') 701C 702C F3L = <03L|K(k1)|0> 703C 704 ILSYM = MULD2H(IREFSY,ISYMV3) 705 IRSYM = IREFSY 706 NCL = MZCONF(ISYMV3) 707 NCR = MZCONF(1) 708 IOFF = MZVAR(ISYMV3) + 1 709 CALL DZERO(DEN1,NASHT*NASHT) 710 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,VEC3(IOFF),CREF, 711 * DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2,XINDX,WRK, 712 * KONE,LWRK) 713 OVLAP = D0 714 IF (ILSYM.EQ.IRSYM) 715 * OVLAP = DDOT(NCL,CREF,1,VEC1(IOFF),1) 716C 717 CALL MELONE(ZYMAT,ISYMV1,DEN1,OVLAP,F3L,IPRONE,'F3L in SCASE2') 718 F = DH*F3L - F3R 719 END IF 720C 721 IF (MZWOPT(ISYMV1).GT.0 .AND. MZWOPT(ISYMV3).GT.0) THEN 722C 723C FACT = <0|K(k1,k3)|0> 724C 725 CALL TRZYMT(NSIM,VEC3,VEC1,KZYV3,KZYV1,ISYMV3,ISYMV1,ZYMAT, 726 * MJWOP,WRK,LWRK) 727 OVLAP = D1 728 CALL MELONE(ZYMAT,1,UDV,OVLAP,FACT,IPRONE,'FACT in SCASE2') 729 F = F + DH*FACT 730 END IF 731C 732 NZCONF = MZCONF(IGRSYM) 733 NZVAR = MZVAR(IGRSYM) 734 CALL DAXPY(NZCONF,F+F1,VEC2,1,S4TRS,1) 735 CALL DAXPY(NZCONF,F+F2,VEC2(NZVAR+1),1,S4TRS(NZVAR+1),1) 736C 737 RETURN 738 END 739 SUBROUTINE SCASE3(KZYVR,KZYV1,KZYV2,KZYV3, 740 * IGRSYM,ISYMV1,ISYMV2,ISYMV3, 741 * S4TRS,VEC1,VEC2,VEC3,CREF, 742 * UDV,ZYMAT,DEN1,XINDX,MJWOP,WRK,LWRK) 743C 744C / <0| [qj,K(k1,k3)] |02R> + <02L| [qj,K(k1,k3)] |0> \ 745C | <j| K(k1,k3) |02R> |*1/2 746C | <0| [qj+,K(k1,k3)] |02R> + <02L| [qj+,K(k1,k3)] |0> | 747C \ -<02L| K(k1,k3) |j> / 748C 749#include "implicit.h" 750#include "infdim.h" 751#include "inforb.h" 752#include "maxorb.h" 753#include "maxash.h" 754#include "infrsp.h" 755#include "wrkrsp.h" 756#include "rspprp.h" 757#include "infhyp.h" 758#include "infvar.h" 759#include "infind.h" 760#include "infpri.h" 761#include "qrinf.h" 762C 763 PARAMETER( DH=0.5D0 ) 764C 765 DIMENSION WRK(*) 766 DIMENSION S4TRS(KZYVR), MJWOP(2,MAXWOP,8) 767 DIMENSION VEC1(KZYV1), VEC2(KZYV2), VEC3(KZYV3) 768 DIMENSION XINDX(*) 769 DIMENSION ZYMAT(NORBT,NORBT) 770 DIMENSION UDV(NASHDI,NASHDI) 771 DIMENSION DEN1(NASHDI,NASHDI) 772 DIMENSION CREF(*) 773C 774 LOGICAL LORB, LCON, LREF, TDM, NORHO2 775C 776 IF (MZWOPT(ISYMV1).EQ.0 .OR. 777 * MZCONF(ISYMV2).EQ.0 .OR. 778 * MZWOPT(ISYMV3).EQ.0) RETURN 779C 780C Initialize variables 781C 782 TDM = .TRUE. 783 NORHO2 = .TRUE. 784 NSIM = 1 785 ISPIN = 0 786 KONE = 1 787C 788C Construct the density matrix <02L|..|0> + <0|..|02R> 789C 790 ILSYM = IREFSY 791 IRSYM = MULD2H(IREFSY,ISYMV2) 792 NCL = MZCONF(1) 793 NCR = MZCONF(ISYMV2) 794 KZVARL = MZCONF(1) 795 KZVARR = MZYVAR(ISYMV2) 796 LREF = .TRUE. 797 ISYMDN = MULD2H(ILSYM,IRSYM) 798 CALL DZERO(DEN1,NASHT*NASHT) 799 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 800 * CREF,VEC2,OVLAP,DEN1,DUMMY,ISPIN,ISPIN,TDM, 801 * NORHO2,XINDX,WRK,KONE,LWRK,LREF) 802C 803C Put the factor into the operator 804C 805 CALL TRZYMT(NSIM,VEC3,VEC1,KZYV3,KZYV1,ISYMV3,ISYMV1,ZYMAT,MJWOP, 806 * WRK,LWRK) 807 CALL DSCAL(NORBT*NORBT,DH,ZYMAT,1) 808C 809C Make the gradient 810C 811 ISYMST = MULD2H(IGRSYM,IREFSY) 812 IF ( ISYMST .EQ. IREFSY ) THEN 813 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 814 ELSE 815 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 816 END IF 817 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 818 LREF = .FALSE. 819 NZYVEC = MZYVAR(ISYMV2) 820 NZCVEC = MZCONF(ISYMV2) 821 CALL RSP1GR(NSIM,KZYVR,IDUM,ISPIN,IGRSYM,ISPIN,ISYMV2,S4TRS, 822 * VEC2,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,ZYMAT, 823 * XINDX,MJWOP,WRK,LWRK,LORB,LCON,LREF) 824C 825 RETURN 826 END 827