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 19C=========================================================================== 20CRevision 1.2 2000/05/24 19:02:25 hjj 21Cnew GETREF calls, fixing CSF problem for triplet 22C=========================================================================== 23 24 SUBROUTINE A3INIT(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 25 * OPLBL,IOPSYM,VEC1,VEC2,A3TRS,XINDX,UDV, 26 * CMO,MJWOP,WRK,LWRK) 27C 28C Layout the core for the calculation of A3 times two vectors 29C 30#include "implicit.h" 31#include "priunit.h" 32#include "infdim.h" 33#include "inforb.h" 34#include "maxorb.h" 35#include "maxash.h" 36#include "infvar.h" 37#include "infrsp.h" 38#include "wrkrsp.h" 39#include "infpri.h" 40C 41 CHARACTER*8 OPLBL 42C 43 DIMENSION WRK(*) 44 DIMENSION A3TRS(KZYVR), MJWOP(2,MAXWOP,8) 45 DIMENSION VEC1(KZYV1) 46 DIMENSION VEC2(KZYV2) 47 DIMENSION XINDX(*) 48 DIMENSION CMO(*) 49 DIMENSION UDV(NASHDI,NASHDI) 50C 51C Initialise the gradient to zero. 52C 53 CALL DZERO(A3TRS,KZYVR) 54C 55 KOPMAT = 1 56 KDEN1 = KOPMAT + NORBT * NORBT 57 KFREE = KDEN1 + NASHT * NASHT 58 LWRKF = LWRK - KFREE + 1 59 IF (LWRKF.LT.0) CALL ERRWRK('A3INIT',KFREE-1,LWRK) 60C 61 IF (IPRRSP .GT. 50 ) THEN 62 WRITE(LUPRI,'(//A)') ' Characteristics of A3 gradient' 63 WRITE(LUPRI,'(A)') ' ==============================' 64 WRITE(LUPRI,'(A,I8)') ' Gradient symmetry :',IGRSYM 65 WRITE(LUPRI,'(A,I8)') ' Gradient length :',KZYVR 66 WRITE(LUPRI,'(A,I8)') ' Vector1 symmetry :',ISYMV1 67 WRITE(LUPRI,'(A,I8)') ' Vector1 length :',KZYV1 68 WRITE(LUPRI,'(A,I8)') ' Vector2 symmetry :',ISYMV2 69 WRITE(LUPRI,'(A,I8)') ' Vector2 length :',KZYV2 70 WRITE(LUPRI,'(A,I8)') ' Operator symmetry :',IOPSYM 71 WRITE(LUPRI,'(A,A8)') ' Operator label :',OPLBL 72 WRITE(LUPRI,'(A//)') ' ==============================' 73 END IF 74C 75 CALL A3DRV(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 76 * OPLBL,IOPSYM,VEC1,VEC2,A3TRS,WRK(KOPMAT), 77 * WRK(KDEN1),UDV,XINDX,CMO,MJWOP,WRK(KFREE),LWRKF) 78C 79 RETURN 80 END 81 SUBROUTINE A3DRV(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 82 * OPLBL,IOPSYM,VEC1,VEC2,A3TRS,OPMAT, 83 * DEN1,UDV,XINDX,CMO,MJWOP,WRK,LWRK) 84C 85C Purpose: 86C Outer driver routine for X[3] times two vectors. This subroutine 87C calls the setup of the operator transformation, constructs a den- 88C sity matrix if necessary and calls RSP1GR to compute the gradient. 89C 90#include "implicit.h" 91C 92#include "maxorb.h" 93#include "infdim.h" 94#include "inforb.h" 95#include "infvar.h" 96#include "qrinf.h" 97#include "infrsp.h" 98C 99 CHARACTER*8 OPLBL 100C 101 DIMENSION A3TRS(KZYVR), MJWOP(2,MAXWOP,8) 102 DIMENSION XINDX(*) 103 DIMENSION OPMAT(NORBT,NORBT) 104 DIMENSION DEN1(NASHDI,NASHDI) 105 DIMENSION UDV(NASHDI,NASHDI) 106 DIMENSION WRK(*) 107 DIMENSION VEC1(KZYV1) 108 DIMENSION VEC2(KZYV2) 109 DIMENSION CMO(*) 110C 111C Get the operator matrix 112C 113 KSYMP = -1 114 CALL PRPGET(OPLBL,CMO,OPMAT,KSYMP,ANTSYM,WRK,LWRK,IPRRSP) 115 IF (KSYMP.NE.IOPSYM) CALL QUIT( 116 & 'A3DRV: unexpected symmetry of operator matrix') 117 CALL DGETRN(OPMAT,NORBT,NORBT) 118C 119 CALL ACASE1(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 120 * IOPSYM,VEC1,VEC2,A3TRS,OPMAT, 121 * DEN1,UDV,XINDX,WRK,LWRK) 122C 123 CALL ACASE2(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 124 * IOPSYM,VEC1,VEC2,A3TRS,OPMAT, 125 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 126C 127 CALL ACASE3(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 128 * IOPSYM,VEC1,VEC2,A3TRS,OPMAT, 129 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 130C 131 CALL ACASE1(KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1, 132 * IOPSYM,VEC2,VEC1,A3TRS,OPMAT, 133 * DEN1,UDV,XINDX,WRK,LWRK) 134C 135 CALL ACASE2(KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1, 136 * IOPSYM,VEC2,VEC1,A3TRS,OPMAT, 137 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 138C 139 CALL ACASE3(KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1, 140 * IOPSYM,VEC2,VEC1,A3TRS,OPMAT, 141 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 142C 143C Swap the final result to conform with the notation for A[3] 144C 145 CALL DSWAP(MZVAR(IGRSYM),A3TRS,1, 146 * A3TRS(MZVAR(IGRSYM)+1),1) 147C 148 RETURN 149 END 150 SUBROUTINE ACASE1(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 151 * IOPSYM,VEC1,VEC2,A3TRS,OPMAT, 152 * DEN1,UDV,XINDX,WRK,LWRK) 153C 154#include "implicit.h" 155C 156 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, D3 = 3D0, D6 = 6.0D0 ) 157C 158#include "maxorb.h" 159#include "maxash.h" 160#include "inforb.h" 161#include "infind.h" 162#include "infvar.h" 163#include "infdim.h" 164#include "infrsp.h" 165#include "wrkrsp.h" 166#include "rspprp.h" 167#include "infhyp.h" 168#include "qrinf.h" 169#include "infpri.h" 170#include "infspi.h" 171C 172 DIMENSION A3TRS(KZYVR) 173 DIMENSION XINDX(*) 174 DIMENSION OPMAT(NORBT,NORBT) 175 DIMENSION DEN1(NASHDI,NASHDI) 176 DIMENSION UDV(NASHDI,NASHDI) 177 DIMENSION WRK(*) 178 DIMENSION VEC1(KZYV1) 179 DIMENSION VEC2(KZYV2) 180C 181 LOGICAL LCON, LORB 182 LOGICAL TDM, NORHO2, LREF 183C 184 IF (MZCONF(ISYMV1) .EQ. 0 .OR. MZCONF(ISYMV2) .EQ. 0) RETURN 185C 186C Initialise variables and layout some workspace 187C 188 TDM = .TRUE. 189 NORHO2 = .TRUE. 190 NSIM = 1 191 INTSYM = 1 192 ISPIN = 0 193 IPRONE = 75 194C 195 KCREF = 1 196C KRES = KCREF + NCREF 197 KRES = KCREF + MZCONF(1) 198 KFREE = KRES + KZYVR 199 LFREE = LWRK - KFREE 200 IF (LFREE.LT.0) CALL ERRWRK('ACASE1',KFREE-1,LWRK) 201C 202 CALL GETREF(WRK(KCREF),MZCONF(1)) 203C 204C Case 1a corresponding to X3 does not exist 205C ======= 206C 207C Case 1b 208C ======= 209C / 0 \ 210C | { -1/3<02L|A|0> - 1/6<0|A|02R> }*Sj(1)' | 211C | 0 | 212C \ { -1/6<02L|A|0> - 1/3<0|A|02R> }*Sj(1) / 213C 214C F2R=<0|A|-02R> 215C 216 ILSYM = IREFSY 217 IRSYM = MULD2H(IREFSY,ISYMV2) 218 NCL = MZCONF(1) 219 NCR = MZCONF(ISYMV2) 220 IOFF = 1 221 CALL DZERO(DEN1,NASHT*NASHT) 222 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,WRK(KCREF),VEC2(IOFF), 223 * DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2,XINDX,WRK, 224 * KFREE,LFREE) 225 OVLAP = D0 226 IF (ILSYM.EQ.IRSYM) 227 * OVLAP = DDOT(NCL,WRK(KCREF),1,VEC2(IOFF),1) 228C 229 CALL MELONE(OPMAT,IOPSYM,DEN1,OVLAP,F2R,IPRONE,'F2R in ACASE1') 230C 231C F2L=<02L|A|0> 232C 233 ILSYM = MULD2H(IREFSY,ISYMV2) 234 IRSYM = IREFSY 235 NCL = MZCONF(ISYMV2) 236 NCR = MZCONF(1) 237 IOFF = MZVAR(ISYMV2) + 1 238 CALL DZERO(DEN1,NASHT*NASHT) 239 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,VEC2(IOFF),WRK(KCREF), 240 * DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2,XINDX,WRK, 241 * KFREE,LFREE) 242 OVLAP = D0 243 IF (ILSYM.EQ.IRSYM) 244 * OVLAP = DDOT(NCL,WRK(KCREF),1,VEC2(IOFF),1) 245C 246 CALL MELONE(OPMAT,IOPSYM,DEN1,OVLAP,F2L,IPRONE,'F2L in ACASE1') 247C 248 NZCONF = MZCONF(IGRSYM) 249 NZVAR = MZVAR(IGRSYM) 250 IF (IGRSYM.EQ.ISYMV1) THEN 251 FACT = - F2L/D3 + F2R/D6 252 CALL DAXPY(NZCONF,FACT,VEC1,1,A3TRS,1) 253 FACT = - F2L/D6 + F2R/D3 254 CALL DAXPY(NZCONF,FACT,VEC1(NZVAR+1),1,A3TRS(NZVAR+1),1) 255 END IF 256C 257C Case 1c 258C ======= 259C / 0 \ 260C | -<0| A |j> | * -1/6 * S(1)S(2)' 261C | 0 | 262C \ <j| A |0> / 263C 264 IF (ISYMV1.NE.ISYMV2) RETURN 265C 266 NZCONF = MZCONF(ISYMV1) 267 NZVAR = MZVAR(ISYMV1) 268 FACT = DDOT(NZCONF,VEC1,1,VEC2(NZVAR+1),1) 269 FACT = -FACT/D6 270C 271 ISYMST = MULD2H(IGRSYM,IREFSY) 272 IF(ISYMST .EQ. IREFSY ) THEN 273 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 274 ELSE 275 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 276 END IF 277 IF (LCON) THEN 278 LREF = .TRUE. 279 CALL DZERO(WRK(KRES),KZYVR) 280 CALL CONSX(NSIM,KZYVR,IGRSYM,OPMAT,WRK(KCREF),MZCONF(1), 281 * MZCONF(1),IREFSY,MZCONF(IGRSYM),ISYMST,LREF, 282 * WRK(KRES),XINDX,WRK(KFREE),LFREE) 283 CALL DAXPY(KZYVR,FACT,WRK(KRES),1,A3TRS,1) 284 END IF 285C 286 RETURN 287 END 288 SUBROUTINE ACASE2(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 289 * IOPSYM,VEC1,VEC2,A3TRS,OPMAT, 290 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 291C 292#include "implicit.h" 293C 294 PARAMETER ( D1 = 1.0D0, DMH = -0.5D0 ) 295C 296#include "maxorb.h" 297#include "maxash.h" 298#include "inforb.h" 299#include "infind.h" 300#include "infvar.h" 301#include "infdim.h" 302#include "infrsp.h" 303#include "wrkrsp.h" 304#include "rspprp.h" 305#include "infhyp.h" 306#include "qrinf.h" 307#include "infpri.h" 308#include "infspi.h" 309C 310 DIMENSION A3TRS(KZYVR), MJWOP(2,MAXWOP,8) 311 DIMENSION XINDX(*) 312 DIMENSION OPMAT(NORBT,NORBT) 313 DIMENSION DEN1(NASHDI,NASHDI) 314 DIMENSION UDV(NASHDI,NASHDI) 315 DIMENSION WRK(*) 316 DIMENSION VEC1(KZYV1), VEC2(KZYV2) 317C 318 LOGICAL LCON, LREF 319C 320C Initialise variables and layout some workspace 321C 322 NSIM = 1 323 IPRONE = 75 324C 325 KZYM = 1 326 KA3X = KZYM + N2ORBX 327 KFREE = KA3X + N2ORBX 328 LFREE = LWRK - KFREE 329 IF (LFREE.LT.0) CALL ERRWRK('ACASE1',KFREE-1,LWRK) 330C 331 IF (MZCONF(ISYMV2) .EQ. 0) RETURN 332 IF ( MZWOPT(ISYMV1) .EQ. 0 ) RETURN 333C 334C Transform the operator with kappa and put the factor -1/2 335C present in these terms into the operator 336C 337 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,WRK(KZYM),MJWOP) 338 CALL DZERO(WRK(KA3X),N2ORBX) 339 CALL OITH1(ISYMV1,WRK(KZYM),OPMAT,WRK(KA3X),IOPSYM) 340 CALL DSCAL(NORBT*NORBT,DMH,WRK(KA3X),1) 341C 342C Case 2a 343C ======= 344C / 0 \ 345C | 1/2<02L| A(k1) |j> | 346C | 0 | 347C \ -1/2<j| A(k1) |02R> / 348C 349C Make the gradient 350C 351 ISYMST = MULD2H(IGRSYM,IREFSY) 352 IF ( ISYMST .EQ. IREFSY ) THEN 353 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 354 ELSE 355 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 356 END IF 357 LREF = .FALSE. 358 NZYVEC = MZYVAR(ISYMV2) 359 NZCVEC = MZCONF(ISYMV2) 360 IF (LCON) THEN 361 CALL CONSX(NSIM,KZYVR,IGRSYM,WRK(KA3X),VEC2,NZYVEC,NZCVEC, 362 * ISYMV2,MZCONF(IGRSYM),ISYMST,LREF,A3TRS, 363 * XINDX,WRK(KFREE),LFREE) 364 END IF 365C 366C Case 2b 367C ======= 368C / 0 \ 369C | Sj(2)' | * -1/2<0| A(k1) |0> 370C | 0 | 371C \ Sj(2) / 372C 373 IF (IGRSYM.EQ.ISYMV2) THEN 374 OVLAP = D1 375 CALL MELONE(WRK(KA3X),1,UDV,OVLAP,FACT,IPRONE,'FACT in ACASE2') 376 NZCONF = MZCONF(IGRSYM) 377 NZVAR = MZVAR(IGRSYM) 378 CALL DAXPY(NZCONF,FACT,VEC2,1,A3TRS,1) 379 CALL DAXPY(NZCONF,FACT,VEC2(NZVAR+1),1,A3TRS(NZVAR+1),1) 380 END IF 381C 382 RETURN 383 END 384 SUBROUTINE ACASE3(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 385 * IOPSYM,VEC1,VEC2,A3TRS,OPMAT, 386 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 387C 388#include "implicit.h" 389C 390 PARAMETER ( D0 = 0.0D0, DMH = -0.5D0 , D1 = 1.0D0, D3 = 3.0D0 ) 391C 392#include "maxorb.h" 393#include "maxash.h" 394#include "inforb.h" 395#include "infind.h" 396#include "infvar.h" 397#include "infdim.h" 398#include "infrsp.h" 399#include "wrkrsp.h" 400#include "rspprp.h" 401#include "infhyp.h" 402#include "qrinf.h" 403#include "infpri.h" 404#include "infspi.h" 405C 406 DIMENSION A3TRS(KZYVR), MJWOP(2,MAXWOP,8) 407 DIMENSION XINDX(*) 408 DIMENSION OPMAT(NORBT,NORBT) 409 DIMENSION DEN1(NASHDI,NASHDI) 410 DIMENSION UDV(NASHDI,NASHDI) 411 DIMENSION WRK(*) 412 DIMENSION VEC1(KZYV1), VEC2(KZYV2) 413C 414 LOGICAL LCON, LORB, LREF 415C 416C Initialise variables and layout some workspace 417C 418 ISPIN = 0 419C 420 KCREF = 1 421 KRES = KCREF + MZCONF(1) 422 KZYM = KRES + KZYVR 423 KA3X = KZYM + N2ORBX 424 KA3XX = KA3X + N2ORBX 425 KFREE = KA3XX + N2ORBX 426 LFREE = LWRK - KFREE 427 IF (LFREE.LT.0) CALL ERRWRK('ACASE3',KFREE-1,LWRK) 428C 429 IF ((MZWOPT(ISYMV2) .EQ. 0).OR.(MZWOPT(ISYMV1) .EQ. 0)) RETURN 430C 431C / 1/3<0| [qj+,A(k1,k2)] |0> \ 432C | -<0| A(k1,k2) |j> | * -1/2 433C | 1/3<0| [qj ,A(k1,k2)] |0> | 434C \ <j| A(k1,k2) |0> / 435C 436C Transform the operator with 2k,1k 437C 438C Put the factor of minus one half present in this term into the 439C ZY matrix used for transforming the integrals 440C 441 NSIM = 1 442 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,WRK(KZYM),MJWOP) 443 CALL DSCAL(NORBT*NORBT,DMH,WRK(KZYM),1) 444 CALL DZERO(WRK(KA3X),N2ORBX) 445 CALL OITH1(ISYMV1,WRK(KZYM),OPMAT,WRK(KA3X),IOPSYM) 446 CALL GTZYMT(NSIM,VEC2,KZYV2,ISYMV2,WRK(KZYM),MJWOP) 447 CALL DZERO(WRK(KA3XX),N2ORBX) 448 ISYMX = MULD2H(IOPSYM,ISYMV1) 449 CALL OITH1(ISYMV2,WRK(KZYM),WRK(KA3X),WRK(KA3XX),ISYMX) 450C 451 CALL GETREF(WRK(KCREF),MZCONF(1)) 452C 453C We have the density matrix over the reference state already 454C 455 ISYMDN = 1 456 OVLAP = D1 457C 458C Make the gradient 459C 460 ISYMST = MULD2H(IGRSYM,IREFSY) 461 IF ( ISYMST .EQ. IREFSY ) THEN 462 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 463 ELSE 464 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 465 END IF 466 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 467 LREF = .TRUE. 468 IF (LCON) THEN 469 CALL DZERO(WRK(KRES),KZYVR) 470 CALL CONSX(NSIM,KZYVR,IGRSYM,WRK(KA3XX),WRK(KCREF),MZCONF(1), 471 * MZCONF(1),IREFSY,MZCONF(IGRSYM),ISYMST,LREF, 472 * WRK(KRES),XINDX,WRK(KFREE),LFREE) 473 CALL DAXPY(KZYVR,D1,WRK(KRES),1,A3TRS,1) 474 END IF 475 IF (LORB) THEN 476 CALL DZERO(WRK(KRES),KZYVR) 477 CALL ORBSX(NSIM,IGRSYM,KZYVR,WRK(KRES),WRK(KA3XX),OVLAP, 478 * ISYMDN,UDV,MJWOP,WRK(KFREE),LFREE) 479 CALL DAXPY(KZYVR,1/D3,WRK(KRES),1,A3TRS,1) 480 END IF 481C 482 RETURN 483 END 484 SUBROUTINE X3INIT(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 485 * OPLBL,IOPSYM,VEC1,VEC2,X3TRS,XINDX,UDV, 486 * CMO,MJWOP,WRK,LWRK) 487C 488C Layout the core for the calculation of X3 times two vectors 489C 490#include "implicit.h" 491#include "priunit.h" 492#include "infdim.h" 493#include "inforb.h" 494#include "maxorb.h" 495#include "maxash.h" 496#include "infvar.h" 497#include "infrsp.h" 498#include "wrkrsp.h" 499#include "rspprp.h" 500#include "infhyp.h" 501#include "qrinf.h" 502#include "infpri.h" 503#include "infspi.h" 504#include "infcr.h" 505C 506 CHARACTER*8 OPLBL 507C 508 DIMENSION WRK(*) 509 DIMENSION X3TRS(KZYVR), MJWOP(2,MAXWOP,8) 510 DIMENSION VEC1(KZYV1), VEC2(KZYV2) 511 DIMENSION XINDX(*) 512 DIMENSION CMO(*) 513 DIMENSION UDV(NASHDI,NASHDI) 514C 515C Initialise the gradient to zero. 516C 517 CALL DZERO(X3TRS,KZYVR) 518C 519 KOPMAT = 1 520 KDEN1 = KOPMAT + NORBT * NORBT 521 KFREE = KDEN1 + NASHT * NASHT 522 LWRKF = LWRK - KFREE + 1 523 IF (LWRKF.LT.0) CALL ERRWRK('X3INIT',KFREE-1,LWRK) 524C 525 IF (IPRRSP .GT. 100 ) THEN 526 WRITE(LUPRI,'(//A)') ' Characteristics of X3 gradient' 527 WRITE(LUPRI,'(A)') ' ==============================' 528 WRITE(LUPRI,'(A,I8)') ' Gradient symmetry :',IGRSYM 529 WRITE(LUPRI,'(A,I8)') ' Gradient length :',KZYVR 530 WRITE(LUPRI,'(A,I8)') ' Vector1 symmetry :',ISYMV1 531 WRITE(LUPRI,'(A,I8)') ' Vector1 length :',KZYV1 532 WRITE(LUPRI,'(A,I8)') ' Vector2 symmetry :',ISYMV2 533 WRITE(LUPRI,'(A,I8)') ' Vector2 length :',KZYV2 534 WRITE(LUPRI,'(A,I8)') ' Operator symmetry :',IOPSYM 535 WRITE(LUPRI,'(A,A8)') ' Operator label :',OPLBL 536 WRITE(LUPRI,'(A//)') ' ==============================' 537 END IF 538C 539 CALL X3DRV(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 540 * OPLBL,IOPSYM,VEC1,VEC2,X3TRS,WRK(KOPMAT), 541 * WRK(KDEN1),UDV,XINDX,CMO,MJWOP,WRK(KFREE),LWRKF) 542C 543 RETURN 544 END 545 SUBROUTINE X3DRV(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 546 * OPLBL,IOPSYM,VEC1,VEC2,X3TRS,OPMAT, 547 * DEN1,UDV,XINDX,CMO,MJWOP,WRK,LWRK) 548C 549C Purpose: 550C Outer driver routine for X[3] times two vectors. This subroutine 551C calls the setup of the operator transformation, constructs a den- 552C sity matrix if necessary and calls RSP1GR to compute the gradient. 553C 554#include "implicit.h" 555C 556#include "maxorb.h" 557#include "infdim.h" 558#include "inforb.h" 559#include "infvar.h" 560#include "infrsp.h" 561C 562 CHARACTER*8 OPLBL 563C 564 DIMENSION X3TRS(KZYVR), MJWOP(2,MAXWOP,8) 565 DIMENSION XINDX(*) 566 DIMENSION OPMAT(NORBT,NORBT) 567 DIMENSION DEN1(NASHDI,NASHDI), UDV(NASHDI,NASHDI) 568 DIMENSION WRK(*) 569 DIMENSION VEC1(KZYV1), VEC2(KZYV2) 570 DIMENSION CMO(*) 571C 572C 573C Get the operator matrix 574C 575 KSYMP = -1 576 CALL PRPGET(OPLBL,CMO,OPMAT,KSYMP,ANTSYM,WRK,LWRK,IPRRSP) 577 IF (KSYMP.NE.IOPSYM) CALL QUIT( 578 & 'X3DRV: unexpected symmetry of operator matrix') 579C 580 CALL XCASE1(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 581 * IOPSYM,VEC1,VEC2,X3TRS,OPMAT, 582 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 583C 584 CALL XCASE2(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 585 * IOPSYM,VEC1,VEC2,X3TRS,OPMAT, 586 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 587C 588 CALL XCASE3(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 589 * IOPSYM,VEC1,VEC2,X3TRS,OPMAT, 590 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 591C 592 CALL XCASE2(KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1, 593 * IOPSYM,VEC2,VEC1,X3TRS,OPMAT, 594 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 595C 596 CALL XCASE3(KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1, 597 * IOPSYM,VEC2,VEC1,X3TRS,OPMAT, 598 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 599C 600 RETURN 601 END 602 SUBROUTINE XCASE1(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 603 * IOPSYM,VEC1,VEC2,X3TRS,OPMAT, 604 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 605C 606#include "implicit.h" 607C 608 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, DH = 0.5D0 ) 609C 610#include "maxorb.h" 611#include "maxash.h" 612#include "inforb.h" 613#include "infind.h" 614#include "infvar.h" 615#include "infdim.h" 616#include "infrsp.h" 617#include "wrkrsp.h" 618#include "rspprp.h" 619#include "infhyp.h" 620#include "qrinf.h" 621#include "infpri.h" 622#include "infspi.h" 623C 624 DIMENSION X3TRS(KZYVR), MJWOP(2,MAXWOP,8) 625 DIMENSION XINDX(*) 626 DIMENSION OPMAT(NORBT,NORBT) 627 DIMENSION DEN1(NASHDI,NASHDI) 628 DIMENSION UDV(NASHDI,NASHDI) 629 DIMENSION WRK(*) 630 DIMENSION VEC1(KZYV1), VEC2(KZYV2) 631C 632 LOGICAL LCON, LORB, TDM, NORHO2, LREF 633C 634C Initialise variables and layout some workspace 635C 636 TDM = .TRUE. 637 NORHO2 = .TRUE. 638 NSIM = 1 639 ISPIN = 0 640 IPRONE = 75 641C 642 KCREF = 1 643 KRES = KCREF + MZCONF(1) 644 KFREE = KRES + KZYVR 645 LFREE = LWRK - KFREE 646 IF (LFREE.LT.0) CALL ERRWRK('XCASE1',KFREE-1,LWRK) 647C 648 IF (MZCONF(ISYMV1) .EQ. 0 .OR. MZCONF(ISYMV2) .EQ. 0) RETURN 649C 650 CALL GETREF(WRK(KCREF),MZCONF(1)) 651C 652C Case 1a 653C ======= 654C / <01L| [qj,X] |02R> + <02L| [qj,X] |01R> \ 655C | 0 | 656C | <01L| [qj+,X] |02R> + <02L| [qj+,X] |01R> | 657C \ 0 / 658C 659C Construct <01L|..|02R> + <02L|..|01R> density 660C 661 ILSYM = MULD2H(IREFSY,ISYMV1) 662 IRSYM = MULD2H(IREFSY,ISYMV2) 663 NCL = MZCONF(ISYMV1) 664 NCR = MZCONF(ISYMV2) 665 KZVARL = MZYVAR(ISYMV1) 666 KZVARR = MZYVAR(ISYMV2) 667 LREF = .FALSE. 668 ISYMDN = MULD2H(ILSYM,IRSYM) 669 CALL DZERO(DEN1,N2ASHX) 670 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 671 * VEC1,VEC2,OVLAP,DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2, 672 * XINDX,WRK,KFREE,LFREE,LREF) 673C 674C Make the gradient 675C 676C 677 IF ( MZWOPT(IGRSYM) .GT. 0 ) THEN 678 CALL ORBSX(NSIM,IGRSYM,KZYVR,X3TRS,OPMAT,OVLAP,ISYMDN, 679 * DEN1,MJWOP,WRK(KFREE),LFREE) 680 END IF 681C 682C Case 1b 683C ======= 684C / 0 \ 685C | { 1/2<02L|X|0> + <0|X|02R> }*Sj(1) | 686C | 0 | + permutation 687C \ { <02L|X|0> + 1/2<0|X|02R> }*Sj(1)' / 688C 689C F1R=<0|X|-01R> 690C 691 ILSYM = IREFSY 692 IRSYM = MULD2H(IREFSY,ISYMV1) 693 NCL = MZCONF(1) 694 NCR = MZCONF(ISYMV1) 695 IOFF = 1 696 CALL DZERO(DEN1,NASHT*NASHT) 697 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,WRK(KCREF),VEC1(IOFF), 698 * DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2,XINDX,WRK, 699 * KFREE,LFREE) 700 OVLAP = D0 701 IF (ILSYM.EQ.IRSYM) 702 * OVLAP = DDOT(NCL,WRK(KCREF),1,VEC1(IOFF),1) 703C 704 CALL MELONE(OPMAT,IOPSYM,DEN1,OVLAP,F1R,IPRONE,'F1R in XCASE1') 705C 706C F2R=<0|X|-02R> 707C 708 IRSYM = MULD2H(IREFSY,ISYMV2) 709 NCR = MZCONF(ISYMV2) 710 CALL DZERO(DEN1,NASHT*NASHT) 711 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,WRK(KCREF),VEC2(IOFF), 712 * DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2,XINDX,WRK, 713 * KFREE,LFREE) 714 OVLAP = D0 715 IF (ILSYM.EQ.IRSYM) 716 * OVLAP = DDOT(NCL,WRK(KCREF),1,VEC2(IOFF),1) 717C 718 CALL MELONE(OPMAT,IOPSYM,DEN1,OVLAP,F2R,IPRONE,'F2R in XCASE1') 719C 720C F1L=<01L|X|0> 721C 722 ILSYM = MULD2H(IREFSY,ISYMV1) 723 IRSYM = IREFSY 724 NCL = MZCONF(ISYMV1) 725 NCR = MZCONF(1) 726 IOFF = MZVAR(ISYMV1) + 1 727 CALL DZERO(DEN1,NASHT*NASHT) 728 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,VEC1(IOFF),WRK(KCREF), 729 * DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2,XINDX,WRK, 730 * KFREE,LFREE) 731 OVLAP = D0 732 IF (ILSYM.EQ.IRSYM) 733 * OVLAP = DDOT(NCL,WRK(KCREF),1,VEC1(IOFF),1) 734C 735 CALL MELONE(OPMAT,IOPSYM,DEN1,OVLAP,F1L,IPRONE,'F1L in XCASE1') 736C 737C F2L=<02L|X|0> 738C 739 ILSYM = MULD2H(IREFSY,ISYMV2) 740 NCL = MZCONF(ISYMV2) 741 IOFF = MZVAR(ISYMV2) + 1 742 CALL DZERO(DEN1,NASHT*NASHT) 743 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,VEC2(IOFF),WRK(KCREF), 744 * DEN1,DUMMY,ISPIN,ISPIN,TDM,NORHO2,XINDX,WRK, 745 * KFREE,LFREE) 746 OVLAP = D0 747 IF (ILSYM.EQ.IRSYM) 748 * OVLAP = DDOT(NCL,WRK(KCREF),1,VEC2(IOFF),1) 749C 750 CALL MELONE(OPMAT,IOPSYM,DEN1,OVLAP,F2L,IPRONE,'F2L in XCASE1') 751C 752 NZCONF = MZCONF(IGRSYM) 753 NZVAR = MZVAR(IGRSYM) 754 IF (IGRSYM.EQ.ISYMV1) THEN 755 FACT = DH*F2L - F2R 756 CALL DAXPY(NZCONF,FACT,VEC1,1,X3TRS,1) 757 FACT = -DH*F2R + F2L 758 CALL DAXPY(NZCONF,FACT,VEC1(NZVAR+1),1,X3TRS(NZVAR+1),1) 759 END IF 760 IF (IGRSYM.EQ.ISYMV2) THEN 761 FACT = DH*F1L - F1R 762 CALL DAXPY(NZCONF,FACT,VEC2,1,X3TRS,1) 763 FACT = -DH*F1R + F1L 764 CALL DAXPY(NZCONF,FACT,VEC2(NZVAR+1),1,X3TRS(NZVAR+1),1) 765 END IF 766C 767C Case 1c 768C ======= 769C / <0| [qj,X] |0> \ 770C | 1/2<j| X |0> | * ( S(1)S(2)' + S(1)'S(2) ) 771C | <0| [qj+ ,X] |0> | 772C \ -1/2<0| X |j> / 773C 774 IF (ISYMV1.NE.ISYMV2) RETURN 775C 776 NZCONF = MZCONF(ISYMV1) 777 NZVAR = MZVAR(ISYMV1) 778 FACT = DDOT(NZCONF,VEC1,1,VEC2(NZVAR+1),1) + 779 * DDOT(NZCONF,VEC2,1,VEC1(NZVAR+1),1) 780C 781 ISYMDN = 1 782 OVLAP = D1 783 ISYMST = MULD2H(IGRSYM,IREFSY) 784 IF(ISYMST .EQ. IREFSY ) THEN 785 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 786 ELSE 787 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 788 END IF 789 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 790 LREF = .TRUE. 791 IF (LCON) THEN 792 CALL DZERO(WRK(KRES),KZYVR) 793 CALL CONSX(NSIM,KZYVR,IGRSYM,OPMAT,WRK(KCREF), 794 * MZCONF(1),MZCONF(1),IREFSY,MZCONF(IGRSYM),ISYMST, 795 * LREF,WRK(KRES),XINDX,WRK(KFREE),LFREE) 796 CALL DSCAL(KZYVR,DH,WRK(KRES),1) 797 CALL DAXPY(KZYVR,FACT,WRK(KRES),1,X3TRS,1) 798 END IF 799 IF (LORB) THEN 800 CALL DZERO(WRK(KRES),KZYVR) 801 CALL ORBSX(NSIM,IGRSYM,KZYVR,WRK(KRES),OPMAT,OVLAP,ISYMDN, 802 * UDV,MJWOP,WRK(KFREE),LFREE) 803 CALL DAXPY(KZYVR,FACT,WRK(KRES),1,X3TRS,1) 804 END IF 805C 806 RETURN 807 END 808 SUBROUTINE XCASE2(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 809 * IOPSYM,VEC1,VEC2,X3TRS,OPMAT, 810 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 811C 812#include "implicit.h" 813C 814 PARAMETER ( D1 = 1.0D0 ) 815C 816#include "maxorb.h" 817#include "maxash.h" 818#include "inforb.h" 819#include "infind.h" 820#include "infvar.h" 821#include "infdim.h" 822#include "infrsp.h" 823#include "wrkrsp.h" 824#include "rspprp.h" 825#include "infhyp.h" 826#include "qrinf.h" 827#include "infpri.h" 828#include "infspi.h" 829C 830 DIMENSION X3TRS(KZYVR), MJWOP(2,MAXWOP,8) 831 DIMENSION XINDX(*) 832 DIMENSION OPMAT(NORBT,NORBT) 833 DIMENSION DEN1(NASHDI,NASHDI) 834 DIMENSION UDV(NASHDI,NASHDI) 835 DIMENSION WRK(*) 836 DIMENSION VEC1(KZYV1) 837 DIMENSION VEC2(KZYV2) 838C 839 LOGICAL LCON, LORB 840 LOGICAL TDM, NORHO2, LREF 841C 842C Initialise variables and layout some workspace 843C 844 ISPIN = 0 845 TDM = .TRUE. 846 NORHO2 = .TRUE. 847 IPRONE = 75 848C 849 KCREF = 1 850 KZYM = KCREF + MZCONF(1) 851 KX3X = KZYM + N2ORBX 852 KFREE = KX3X + N2ORBX 853 LFREE = LWRK - KFREE 854 IF (LFREE.LT.0) CALL ERRWRK('XCASE1',KFREE-1,LWRK) 855C 856 IF ( MZCONF(ISYMV2) .EQ. 0 ) RETURN 857 IF ( MZWOPT(ISYMV1) .EQ. 0 ) RETURN 858C 859C Transform the operator with kappa 860C 861 NSIM = 1 862 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,WRK(KZYM),MJWOP) 863 CALL DZERO(WRK(KX3X),N2ORBX) 864 CALL OITH1(ISYMV1,WRK(KZYM),OPMAT,WRK(KX3X),IOPSYM) 865C 866 CALL GETREF(WRK(KCREF),MZCONF(1)) 867C 868C Case 2a 869C ======= 870C / <0| [qj,X(k1)] |02R> + <02L| [qj,X(k1)] |0> \ 871C | <j| X(k1) |02R> | 872C | <0| [qj+,X(k1)] |02R> + <02L| [qj+,X(k1)] |0> | 873C \ -<02L| X(k1) |j> / 874C 875C Construct the density matrix <0L|..|0> + <0|..|0R> 876C 877 ILSYM = IREFSY 878 IRSYM = MULD2H(IREFSY,ISYMV2) 879 NCL = MZCONF(1) 880 NCR = MZCONF(ISYMV2) 881 KZVARL = MZCONF(1) 882 KZVARR = MZYVAR(ISYMV2) 883 LREF = .TRUE. 884 ISYMDN = MULD2H(ILSYM,IRSYM) 885 CALL DZERO(DEN1,NASHT*NASHT) 886 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 887 * WRK(KCREF),VEC2,OVLAP,DEN1,DUMMY,ISPIN,ISPIN,TDM, 888 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 889C 890C Make the gradient 891C 892 ISYMST = MULD2H(IGRSYM,IREFSY) 893 IF ( ISYMST .EQ. IREFSY ) THEN 894 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 895 ELSE 896 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 897 END IF 898 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 899 LREF = .FALSE. 900 NZYVEC = MZYVAR(ISYMV2) 901 NZCVEC = MZCONF(ISYMV2) 902 CALL RSP1GR(NSIM,KZYVR,IDUM,ISPIN,IGRSYM,ISPIN,ISYMV2,X3TRS, 903 * VEC2,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,WRK(KX3X), 904 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 905C 906C Case 2b 907C ======= 908C / 0 \ 909C | Sj(2) | * <0| X(k1) |0> 910C | 0 | 911C \ Sj(2)' / 912C 913 IF (IGRSYM.EQ.ISYMV2) THEN 914 OVLAP = D1 915 CALL MELONE(WRK(KX3X),1,UDV,OVLAP,FACT,IPRONE,'FACT in XCASE2') 916 NZCONF = MZCONF(IGRSYM) 917 NZVAR = MZVAR(IGRSYM) 918 CALL DAXPY(NZCONF,FACT,VEC2,1,X3TRS,1) 919 CALL DAXPY(NZCONF,FACT,VEC2(NZVAR+1),1,X3TRS(NZVAR+1),1) 920 END IF 921C 922 RETURN 923 END 924 SUBROUTINE XCASE3(KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2, 925 * IOPSYM,VEC1,VEC2,X3TRS,OPMAT, 926 * DEN1,UDV,XINDX,MJWOP,WRK,LWRK) 927C 928#include "implicit.h" 929C 930 PARAMETER ( D0 = 0.0D0, DH = 0.5D0 , D1 = 1.0D0 ) 931C 932#include "maxorb.h" 933#include "maxash.h" 934#include "inforb.h" 935#include "infind.h" 936#include "infvar.h" 937#include "infdim.h" 938#include "infrsp.h" 939#include "wrkrsp.h" 940#include "rspprp.h" 941#include "infhyp.h" 942#include "qrinf.h" 943#include "infpri.h" 944#include "infspi.h" 945C 946 DIMENSION X3TRS(KZYVR), MJWOP(2,MAXWOP,8) 947 DIMENSION XINDX(*) 948 DIMENSION OPMAT(NORBT,NORBT) 949 DIMENSION DEN1(NASHDI,NASHDI) 950 DIMENSION UDV(NASHDI,NASHDI) 951 DIMENSION WRK(*) 952 DIMENSION VEC1(KZYV1), VEC2(KZYV2) 953C 954 LOGICAL LCON, LORB, TDM, LREF 955C 956C Initialise variables and layout some workspace 957C 958 ISPIN = 0 959 TDM = .TRUE. 960C 961 KCREF = 1 962 KZYM = KCREF + MZCONF(1) 963 KX3X = KZYM + N2ORBX 964 KX3XX = KX3X + N2ORBX 965 KFREE = KX3XX + N2ORBX 966 LFREE = LWRK - KFREE 967 IF (LFREE.LT.0) CALL ERRWRK('XCASE1',KFREE-1,LWRK) 968C 969 IF ((MZWOPT(ISYMV2) .EQ. 0).OR.(MZWOPT(ISYMV1) .EQ. 0)) RETURN 970C 971C / <0| [qj ,X(k1,k2)] |0> \ 972C | <j| X(k1,k2) |0> | * 1/2 973C | <0| [qj+,X(k1,k2)] |0> | 974C \ -<0| X(k1,k2) |j> / 975C 976C Transform the operator with 2k,1k 977C 978C Put the factor of one half present in this term into the 979C ZY matrix used for transforming the integrals 980C 981 NSIM = 1 982 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,WRK(KZYM),MJWOP) 983 CALL DSCAL(NORBT*NORBT,DH,WRK(KZYM),1) 984 CALL DZERO(WRK(KX3X),N2ORBX) 985 CALL OITH1(ISYMV1,WRK(KZYM),OPMAT,WRK(KX3X),IOPSYM) 986 CALL GTZYMT(NSIM,VEC2,KZYV2,ISYMV2,WRK(KZYM),MJWOP) 987 CALL DZERO(WRK(KX3XX),N2ORBX) 988 ISYMX = MULD2H(IOPSYM,ISYMV1) 989 CALL OITH1(ISYMV2,WRK(KZYM),WRK(KX3X),WRK(KX3XX),ISYMX) 990C 991 CALL GETREF(WRK(KCREF),MZCONF(1)) 992C 993C We have the density matrices over the reference state already 994C 995 ISYMDN = 1 996 OVLAP = D1 997C 998C Make the gradient 999C 1000 ISYMV = IREFSY 1001 ISYMST = MULD2H(IGRSYM,IREFSY) 1002 IF ( ISYMST .EQ. IREFSY ) THEN 1003 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 1004 ELSE 1005 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 1006 END IF 1007 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 1008 LREF = .TRUE. 1009 NZYVEC = MZCONF(1) 1010 NZCVEC = MZCONF(1) 1011 CALL RSP1GR(NSIM,KZYVR,IDUM,ISPIN,IGRSYM,ISPIN,ISYMV,X3TRS, 1012 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,WRK(KX3XX), 1013 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 1014C 1015 RETURN 1016 END 1017