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! 18#if defined (VAR_DEBUG) 19#define HSODEBUG .TRUE. 20#else 21#define HSODEBUG .FALSE. 22#endif 23 24C 25C======================================================================= 26C May 2000 hjj: 27C Use MZCONF(1) instead NCREF for CREF (for triplet response w. CSF) 28C Removed delete of LUMHSO which caused trouble because LUMHSO was neeeded later! 29C======================================================================= 30 31C /* Deck x2hso */ 32 SUBROUTINE X2HSO(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 33 * OPLBL,IOPSYM,IOPSPI, 34 * GRVEC,X2TRS,VEC,UDV,PV,OPMAT,DEN1,DEN2,XINDX,CMO, 35 * MJWOP,WRK,LWRK) 36C 37C HSO[2]*N X2 style, i.e. 38C 39C ( <0|[q,HSO(N)]|0> ) 40C HSO[2]*N = - ( <j|HSO(N)|0> ) (Case 1) 41C ( <0|[q+,HSO(N)]|0> ) 42C ( -<0|HSO(N)|j> ) 43C 44C 45C ( <L|[q,HSO]|0> + <0|[q,HSO]|R> ) 46C - ( <j|HSO|R> ) (Case 2) 47C ( <L|[q+,HSO]|0> + <0|[q+,HSO]|R> ) 48C ( -<L|HSO|j> ) 49C 50C 51C ( 0 ) 52C - <0|HSO|0> ( Sj' ) (Case 3) 53C ( 0 ) 54C ( Sj ) 55C 56C Case 3 does not contribute for orbitally non-degenerate states (assumed) 57C 58#include "implicit.h" 59#include "dummy.h" 60C 61 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, DM1 = -1.0D0 ) 62C 63#include "maxorb.h" 64#include "maxash.h" 65#include "priunit.h" 66#include "inforb.h" 67#include "infind.h" 68#include "infvar.h" 69#include "infdim.h" 70#include "infrsp.h" 71#include "wrkrsp.h" 72#include "rspprp.h" 73#include "infhyp.h" 74#include "qrinf.h" 75#include "infpri.h" 76#include "inftap.h" 77#include "infspi.h" 78#include "infhso.h" 79#include "trhso.h" 80#include "codata.h" 81C 82C Used from common: 83C 84 CHARACTER*8 OPLBL 85 DIMENSION WRK(*) 86 DIMENSION X2TRS(KZYVR), MJWOP(2,MAXWOP,8) 87 DIMENSION GRVEC(KZYVR) 88 DIMENSION VEC(KZYV) 89 DIMENSION DEN1(NASHDI,NASHDI), DEN2(*) 90 DIMENSION OPMAT(NORBT,NORBT) 91 DIMENSION XINDX(*), CMO(*) 92 DIMENSION UDV(NASHDI,NASHDI) 93C 94 LOGICAL NORHO2, TDM, LORB, LCON, LREFST, LREF 95 LOGICAL DO1, DO2 96 CHARACTER*8 LABEL,HSOLBL(3) 97 LOGICAL ACALL 98 COMMON /CB_HSO_ACALL/ACALL 99 ACALL = .FALSE. 100 DATA HSOLBL/'X SPNORB','Y SPNORB','Z SPNORB'/ 101 102 IPRHSO = MAX(IPRRSP, IPRHSO) 103C 104C If expicit one-electron part go normal way 105C 106 HSOFAC=ALPHAC**2/4 107 IF (OPLBL(2:2).EQ.'1') THEN 108 CALL X2DRV(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 109 * OPLBL,IOPSYM,IOPSPI,GRVEC,X2TRS, 110 * VEC,UDV,OPMAT,DEN1, 111 * XINDX,CMO,MJWOP,WRK,LWRK) 112 CALL DSCAL(KZYVR,HSOFAC,X2TRS,1) 113 RETURN 114 END IF 115C 116C Everything from AO integrals 117C 118 IF (TDHF) THEN 119 CALL X2HSOAO( 120 & OPLBL,CMO,MJWOP, 121 & KZYV,ISYMV,ISPINV,VEC, 122 & KZYVR,IGRSYM,IGRSPI,X2TRS, 123 & WRK,LWRK) 124 CALL DSCAL(KZYVR,-HSOFAC,X2TRS,1) 125 RETURN 126 END IF 127 IF (OPLBL(2:2).EQ.'2') THEN 128 DO1 = .FALSE. 129 DO2 = .TRUE. 130 END IF 131 IF (OPLBL(2:2).EQ.' ') THEN 132 DO1 = DOSO1 133 DO2 = DOSO2 134 END IF 135 LABEL = OPLBL 136 LABEL(2:2) = '1' 137 CALL QENTER('X2HSO') 138C 139C Layout some workspace 140C 141 KCREF = 1 142 KZYM = KCREF + MZCONF(1) 143 KX2X = KZYM + N2ORBX 144 KFREE = KX2X + N2ORBX 145 LFREE = LWRK - KFREE + 1 146 IF (LFREE.LT.0) CALL ERRWRK('X2HSO',KFREE-1,LWRK) 147C 148C Transform if necessary: last transformed component: ILXYZ 149C 150 IF (OPLBL(1:1).NE.HSOLBL(ILXYZ)(1:1)) THEN 151 IF (OPLBL(1:1) .EQ. 'X') ILXYZ=1 152 IF (OPLBL(1:1) .EQ. 'Y') ILXYZ=2 153 IF (OPLBL(1:1) .EQ. 'Z') ILXYZ=3 154 KSYMSO = IOPSYM 155 ITRLSO = 2 156 IF (IPRHSO.GT.0) WRITE(LUPRI,'(/A,I3)') 157 & ' X2HSO: Transforming 2-el. spin-orbit component', ILXYZ 158 CALL GPOPEN(LUAHSO,'AO2SOINT','OLD',' ','UNFORMATTED',IDUMMY, 159 & .FALSE.) 160 CALL TRAHSO(ITRLSO,CMO,WRK(KFREE),LFREE) 161 CALL GPCLOSE(LUAHSO,'KEEP') 162 END IF 163 IF (X2TEST) THEN 164 X2TMZC = D0 165 X2TMZO = D0 166 X2TMYC = D0 167 X2TMYO = D0 168 X2TM = D0 169 KZX2O = MZWOPT(IGRSYM) 170 KZX2C = MZCONF(IGRSYM) 171 KZC = 1 172 KZO = KZC + KZX2C 173 KYC = KZO + KZX2O 174 KYO = KYC + KZX2C 175 END IF 176C 177 KSYMOP = IOPSYM 178 TDM = .TRUE. 179 NORHO2 = .NOT.DOSO2 180 IDAX = LUMHSO 181C 182C Get the operator matrix 183C Put the minus sign in the whole term at the end 184C 185 IF (DO1) THEN 186 KSYMP = -1 187 CALL PRPGET(LABEL,CMO,OPMAT,KSYMP,ANTSYM,WRK,LWRK,IPRHSO) 188 IF (KSYMP .NE. KSYMOP) CALL QUIT( 189 & 'ERROR: Unexpected symmetry of property matrix') 190 ELSE 191 CALL DZERO(OPMAT,N2ORBX) 192 END IF 193C 194C Case 1 195C ====== 196 IF (MZWOPT(ISYMV) .EQ. 0 ) GOTO 2000 197C 198 ISPIN = MULSP(ISPINV,IOPSPI) 199 IF (ISPIN.NE.IGRSPI) CALL QUIT('X2HSO: SPIN ERROR') 200C 201C Transform the operator 202C 203 CALL GTZYMT(NSIM,VEC,KZYV,ISYMV,WRK(KZYM),MJWOP) 204 CALL DZERO(WRK(KX2X),N2ORBX) 205 CALL OITH1(ISYMV,WRK(KZYM),OPMAT,WRK(KX2X),IOPSYM) 206C 207C Make copies of the MC density matrix in DEN1 208C 209 CALL DCOPY(NASHT*NASHT,UDV,1,DEN1,1) 210 OVLAP = D1 211 ISYMDN = 1 212 IKLVL = 1 213C 214C Make the gradient 215C 216 ISYMR = IREFSY 217 INTSYM = KSYMOP 218 ISYMST = MULD2H(IGRSYM,IREFSY) 219 IF ( ISYMST .EQ. IREFSY ) THEN 220 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 221 ELSE 222 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 223 END IF 224 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 225 NZYVEC = MZCONF(1) 226 NZCVEC = MZCONF(1) 227 LREFST = .TRUE. 228 CALL IPSET(0,0,1,1) 229 CALL HSO2CR(IGRSYM,IGRSPI,X2TRS,VDUMMY,NZYVEC,NZCVEC,ISYMV,ISYMDN, 230 * XINDX,OVLAP,UDV,PV,WRK(KX2X),WRK(KFREE),LFREE,KZYVR, 231 * LCON,LORB,LREFST,IDAX,INTSYM,IOPSPI,ISPINV, 232 * IKLVL,WRK(KZYM),ISYMV,ZYMC,ISYMC,MJWOP) 233C CALL HSO2CR(IGRSYM,IGRSPI,ETRS,VEC,NZYVEC,NZCVEC,ISYMV,ISYMDN, 234C * XINDX,OVLAP,DEN1,DEN2,FI,WRK,LWRK,KZYVR, 235C * LCON,LORB,LREFST,IDAX,INTSYM,ISPIN1,ISPIN2, 236C * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2,MJWOP) 237C CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,X2TRS, 238C * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,WRK(KX2X), 239C * XINDX,WRK(KFREE),LFREE,LORB,LCON,LREFST) 240C 241 IF (X2TEST) THEN 242 X2ZO = -DDOT(KZX2O,GRVEC(KZO),1,X2TRS(KZO),1) 243 X2ZC = -DDOT(KZX2C,GRVEC(KZC),1,X2TRS(KZC),1) 244 X2YO = -DDOT(KZX2O,GRVEC(KYO),1,X2TRS(KYO),1) 245 X2YC = -DDOT(KZX2C,GRVEC(KYC),1,X2TRS(KYC),1) 246 X2TOT = X2ZO+X2ZC+X2YO+X2YC 247 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 ZO :',X2ZO - X2TMZO 248 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 ZC :',X2ZC - X2TMZC 249 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 YO :',X2YO - X2TMYO 250 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 YC :',X2YC - X2TMYC 251 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 :',X2TOT - X2TM 252 X2TMZO = X2ZO 253 X2TMZC = X2ZC 254 X2TMYO = X2YO 255 X2TMYC = X2YC 256 X2TM = X2TOT 257 END IF 258 IF ( IPRRSP .GT. 100 ) THEN 259 WRITE(LUPRI,'(A)') ' Case 1 for X2 term' 260 CALL OUTPUT(X2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 261 END IF 262C 263C Case 2 264C ====== 265 2000 IF (MZCONF(ISYMV) .EQ. 0 ) GOTO 3000 266C Construct the density matrix <02L|..|0> + <0|..|02R> 267C 268C 269 CALL GETREF(WRK(KCREF),MZCONF(1)) 270C 271 CALL DCOPY(N2ORBX,OPMAT,1,WRK(KX2X),1) 272 CALL DZERO(DEN1,NASHT*NASHT) 273 OVLAP = D0 274 IKLVL = 0 275C 276 ILSYM = IREFSY 277 IRSYM = MULD2H(IREFSY,ISYMV) 278 NCL = MZCONF(1) 279 NCR = MZCONF(ISYMV) 280 KZVARL = MZCONF(1) 281 KZVARR = MZYVAR(ISYMV) 282 LREF = .TRUE. 283 ISYMDN = MULD2H(ILSYM,IRSYM) 284C 285C Densities for orbital gradient: <e(+,-)>,<e(-,+)> for IGRSPI = 0 286C Densities for orbital gradient: <e(-,-)>,<e(+,+)> for IGRSPI = 1 287C 288 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 289 * WRK(KCREF),VEC,OVLAP,DEN1,DEN2,IGRSPI,1, 290 * TDM,NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 291 IF (IGRSPI.EQ.1) THEN 292 IS1 = 0 293 IS2 = 0 294 ELSE 295 IS1 = 1 296 IS2 = IGRSPI 297 END IF 298 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 299 * WRK(KCREF),VEC,OVLAP,DEN1,DEN2(1+N2ASHX*N2ASHX), 300 * IS1,IS2, 301 * TDM,NORHO2,XINDX, 302 * WRK,KFREE,LFREE,LREF) 303 CALL IPSET(IGRSPI,1,IS1,IS2) 304C 305C Make the gradient 306C 307 ISYMR = ISYMV 308 INTSYM = KSYMOP 309 ISYMST = MULD2H(IGRSYM,IREFSY) 310 IF ( ISYMST .EQ. IREFSY ) THEN 311 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 312 ELSE 313 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 314 END IF 315 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 316 NZYVEC = MZYVAR(ISYMV) 317 NZCVEC = MZCONF(ISYMV) 318 LREFST = .FALSE. 319 CALL HSO2CR(IGRSYM,IGRSPI,X2TRS,VEC,NZYVEC,NZCVEC,ISYMV,ISYMDN, 320 * XINDX,OVLAP,DEN1,DEN2,WRK(KX2X), 321 * WRK(KFREE),LFREE,KZYVR, 322 * LCON,LORB,LREFST,IDAX,INTSYM,IOPSPI,0, 323 * IKLVL,ZYMB,ISYMB,ZYMC,ISYMC,MJWOP) 324C CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,X2TRS, 325C * VEC,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,OPMAT, 326C * XINDX,WRK,LWRK,LORB,LCON,LREFST) 327C 328 IF (X2TEST) THEN 329 X2ZO = -DDOT(KZX2O,GRVEC(KZO),1,X2TRS(KZO),1) 330 X2ZC = -DDOT(KZX2C,GRVEC(KZC),1,X2TRS(KZC),1) 331 X2YO = -DDOT(KZX2O,GRVEC(KYO),1,X2TRS(KYO),1) 332 X2YC = -DDOT(KZX2C,GRVEC(KYC),1,X2TRS(KYC),1) 333 X2TOT = X2ZO+X2ZC+X2YO+X2YC 334 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 ZO :',X2ZO - X2TMZO 335 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 ZC :',X2ZC - X2TMZC 336 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 YO :',X2YO - X2TMYO 337 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 YC :',X2YC - X2TMYC 338 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 :',X2TOT - X2TM 339 X2TMZO = X2ZO 340 X2TMZC = X2ZC 341 X2TMYO = X2YO 342 X2TMYC = X2YC 343 X2TM = X2TOT 344 END IF 345 IF ( IPRRSP .GT. 100 ) THEN 346 WRITE(LUPRI,'(A)') ' Case 2 for X2 term' 347 CALL OUTPUT(X2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 348 END IF 349C 350C Skip expectation value of imaginary operator 351C 352C Case 3 353C ====== 354C Do Sj(2) <0 |op| 0> 355C 356 3000 CONTINUE 357 IF (X2TEST) THEN 358 WRITE(LUPRI,'(/A,F24.8)')' X2TEST Final result: ZO', X2ZO 359 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Final result: ZC', X2ZC 360 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Final result: YO', X2YO 361 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Final result: YC', X2YC 362 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Final result: ', X2TOT 363 END IF 364C 365C Take care of the minus sign of the whole term 366C and transform to atomic units 367C 368 CALL DSCAL(KZYVR,-HSOFAC,X2TRS,1) 369 CALL QEXIT('X2HSO') 370 RETURN 371 END 372C /* Deck a2hso */ 373 SUBROUTINE A2HSO(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 374 * OPLBL,IOPSYM,IOPSPI, 375 * GRVEC,A2TRS,VEC,UDV,PV,OPMAT,XINDX,CMO, 376 * MJWOP,WRK,LWRK) 377C 378C 379C HSO[2]*N A2 style, i.e. 380C 381C ( 1/2 <0|[q+,HSO(N)]|0> ) 382C HSO[2]*N = - ( <0|HSO(N)|j> ) (Case 1) 383C ( 1/2 <0|[q,HSO(N)]|0> ) 384C ( <j|HSO(N)|0> ) 385C 386C 387C ( 0 ) 388C - ( - <L|HSO|j> ) (Case 2) 389C ( 0 ) 390C ( <j|HSO|R> ) 391C 392C 393C ( 0 ) 394C - 1/2 <0|HSO|0> ( Sj' ) (Case 3) 395C ( 0 ) 396C ( Sj ) 397C 398C Case 3 does not contribute for orbitally non-degenerate states (assumed) 399C 400C Drive the computation of A[2] times one vector 401C 402#include "implicit.h" 403#include "dummy.h" 404C 405 PARAMETER( D0 = 0.0D0, D1 = 1.0D0, DH = 0.5D0, DM1 = -1.0D0) 406 PARAMETER (D2 = 2.0D0) 407C 408#include "maxorb.h" 409#include "maxash.h" 410#include "priunit.h" 411#include "infdim.h" 412#include "inftap.h" 413#include "inforb.h" 414#include "infrsp.h" 415#include "wrkrsp.h" 416#include "infvar.h" 417#include "infind.h" 418#include "qrinf.h" 419#include "infpri.h" 420#include "infspi.h" 421#include "trhso.h" 422#include "infhso.h" 423#include "codata.h" 424C 425 CHARACTER*8 OPLBL,HSOLBL(3) 426 DATA HSOLBL/'X SPNORB','Y SPNORB','Z SPNORB'/ 427C 428 DIMENSION WRK(*) 429 DIMENSION A2TRS(KZYVR), MJWOP(2,MAXWOP,8) 430 DIMENSION GRVEC(KZYVR) 431 DIMENSION VEC(KZYV) 432 DIMENSION OPMAT(NORBT,NORBT) 433 DIMENSION XINDX(*) 434 DIMENSION UDV(NASHDI,NASHDI), PV(*) 435 DIMENSION CMO(*) 436C 437 CHARACTER*8 LABEL 438 LOGICAL NORHO2, TDM, LORB, LCON, LREFST, LREF 439 LOGICAL DO1, DO2 440 LOGICAL ACALL 441 COMMON /CB_HSO_ACALL/ACALL 442 ACALL = .TRUE. 443C 444C If expicit one-electron part go normal way 445C 446 HSOFAC=ALPHAC**2/4 447 IF (OPLBL(2:2).EQ.'1') THEN 448 CALL A2DRV(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 449 * OPLBL,IOPSYM,IOPSPI, 450 * GRVEC,A2TRS,VEC,UDV,OPMAT,XINDX,CMO,MJWOP, 451 * WRK,LWRK) 452 CALL DSCAL(KZYVR,HSOFAC,A2TRS,1) 453 RETURN 454 END IF 455C 456C Everything from AO integrals 457C 458 IF (TDHF) THEN 459 CALL X2HSOAO( 460 & OPLBL,CMO,MJWOP, 461 & KZYV,ISYMV,ISPINV,VEC, 462 & KZYVR,IGRSYM,IGRSPI,A2TRS, 463 & WRK,LWRK) 464 CALL DSCAL(KZYVR,-HSOFAC/2,A2TRS,1) 465 CALL DSWAP(MZVAR(IGRSYM),A2TRS,1, 466 * A2TRS(MZVAR(IGRSYM)+1),1) 467 RETURN 468 END IF 469 IF (OPLBL(2:2).EQ.'2') THEN 470 DO1 = .FALSE. 471 DO2 = .TRUE. 472 END IF 473 IF (OPLBL(2:2).EQ.' ') THEN 474 DO1 = DOSO1 475 DO2 = DOSO2 476 END IF 477 LABEL = OPLBL 478 LABEL(2:2) = '1' 479 CALL QENTER('A2HSO') 480C 481C Layout some workspace 482C 483 KZYM = 1 484 KA2X = KZYM + N2ORBX 485 KFREE = KA2X + N2ORBX 486 LFREE = LWRK - KFREE 487 IF (LFREE.LT.0) CALL ERRWRK('A2HSO',KFREE-1,LWRK) 488C 489C Transform if necessary: last transformed component: ILXYZ 490C 491 IF (OPLBL(1:1).NE.HSOLBL(ILXYZ)(1:1)) THEN 492 IF (OPLBL(1:1) .EQ. 'X') ILXYZ=1 493 IF (OPLBL(1:1) .EQ. 'Y') ILXYZ=2 494 IF (OPLBL(1:1) .EQ. 'Z') ILXYZ=3 495 KSYMSO = IOPSYM 496 ITRLSO = 2 497 IF (IPRHSO.GT.0) WRITE(LUPRI,'(/A,I3)') 498 & ' A2HSO: Transforming 2-el. spin-orbit component', ILXYZ 499 CALL GPOPEN(LUAHSO,'AO2SOINT','OLD',' ','UNFORMATTED', 500 & IDUMMY,.FALSE.) 501 CALL TRAHSO(ITRLSO,CMO,WRK(KFREE),LFREE) 502 CALL GPCLOSE(LUAHSO,'KEEP') 503 END IF 504C 505 TDM = .TRUE. 506 NORHO2 = .NOT.DO2 507C 508 IF (A2TEST) THEN 509 A2TMZC = D0 510 A2TMZO = D0 511 A2TMYC = D0 512 A2TMYO = D0 513 A2TM = D0 514 KZA2O = MZWOPT(IGRSYM) 515 KZA2C = MZCONF(IGRSYM) 516 KZC = 1 517 KZO = KZC + KZA2C 518 KYC = KZO + KZA2O 519 KYO = KYC + KZA2C 520 END IF 521C 522C Get the operator matrix 523C 910408-hjaaj: transpose the operator matrix to 524C get right sign for imaginary operators 525C 526 KSYMOP = IOPSYM 527 IDAX = LUMHSO 528 IF (DO1) THEN 529 KSYMP = -1 530 CALL PRPGET(LABEL,CMO,OPMAT,KSYMP,ANTSYM,WRK,LWRK,IPRHSO) 531 IF (KSYMP .NE. KSYMOP) CALL QUIT( 532 & 'ERROR: Unexpected symmetry of property matrix') 533 CALL DGETRN(OPMAT,NORBT,NORBT) 534 ELSE 535 CALL DZERO(OPMAT,N2ORBX) 536 END IF 537C 538C Case 1: 539C ======= 540C / 1/2 <0| [qj+,A(K)] |0> \ 541C | - <0| A(K) |j> | 542C | 1/2 <0| [qj ,A(K)] |0> | 543C \ <j| A(K) |0> / 544C 545 IF ( MZWOPT(ISYMV) .EQ. 0 ) GOTO 2000 546C 547 ISPIN = MULSP(ISPINV,IOPSPI) 548 IF (ISPIN.NE.IGRSPI) CALL QUIT('A2HSO: SPIN ERROR') 549C 550C Transform the operator 551C 552 INTSYM = IOPSYM 553 CALL GTZYMT(NSIM,VEC,KZYV,ISYMV,WRK(KZYM),MJWOP) 554 CALL DZERO(WRK(KA2X),N2ORBX) 555 CALL OITH1(ISYMV,WRK(KZYM),OPMAT,WRK(KA2X),IOPSYM) 556C 557C Make the gradient 558C 559 OVLAP = D1 560 IKLVL = 1 561 ISYMDN = 1 562 ISYMR = IREFSY 563 ISYMST = MULD2H(IGRSYM,IREFSY) 564 IF(ISYMST .EQ. IREFSY ) THEN 565 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 566 ELSE 567 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 568 END IF 569 LORB = (MZWOPT(IGRSYM) .GT. 0 ) 570 NZYVEC = MZCONF(1) 571 NZCVEC = MZCONF(1) 572 LREFST = .TRUE. 573 CALL IPSET(0,0,1,1) 574 IF (LCON) THEN 575 CALL HSO2CR(IGRSYM,IGRSPI,A2TRS,VDUMMY,NZYVEC,NZCVEC, 576 * ISYMV,ISYMDN, 577 * XINDX,OVLAP,UDV,PV,WRK(KA2X),WRK(KFREE),LFREE,KZYVR, 578 * LCON,.FALSE.,LREFST,IDAX,INTSYM,IOPSPI,ISPINV, 579 * IKLVL,WRK(KZYM),ISYMV,ZYMC,ISYMC,MJWOP) 580C CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,A2TRS, 581C * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV, 582C * WRK(KA2X),XINDX,WRK(KFREE),LFREE, 583C * .FALSE.,LCON,LREFST) 584 END IF 585 IF (LORB) THEN 586C multiply A(K) with 1/2 for orbital part 587C multiply ZYMAT with 1/2 for orbital part 588 CALL DSCAL(N2ORBX,DH,WRK(KZYM),1) 589 CALL DZERO(WRK(KA2X),N2ORBX) 590 CALL OITH1(ISYMV,WRK(KZYM),OPMAT,WRK(KA2X),IOPSYM) 591 CALL HSO2CR(IGRSYM,IGRSPI,A2TRS,VDUMMY,NZYVEC,NZCVEC, 592 * ISYMV,ISYMDN, 593 * XINDX,OVLAP,UDV,PV,WRK(KA2X),WRK(KFREE),LFREE,KZYVR, 594 * .FALSE.,LORB,LREFST,IDAX,INTSYM,IOPSPI,ISPINV, 595 * IKLVL,WRK(KZYM),ISYMV,ZYMC,ISYMC,MJWOP) 596C CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,A2TRS, 597C * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV, 598C * WRK(KA2X),XINDX,WRK(KFREE),LFREE, 599C * LORB,.FALSE.,LREFST) 600 END IF 601C 602 IF (A2TEST) THEN 603 A2ZO = DDOT(KZA2O,GRVEC(KZO),1,A2TRS(KYO),1) 604 A2ZC = DDOT(KZA2C,GRVEC(KZC),1,A2TRS(KYC),1) 605 A2YO = DDOT(KZA2O,GRVEC(KYO),1,A2TRS(KZO),1) 606 A2YC = DDOT(KZA2C,GRVEC(KYC),1,A2TRS(KZC),1) 607 A2TOT = A2ZO+A2ZC+A2YO+A2YC 608 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 ZO :',A2ZO - A2TMZO 609 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 ZC :',A2ZC - A2TMZC 610 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 YO :',A2YO - A2TMYO 611 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 YC :',A2YC - A2TMYC 612 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 :',A2TOT - A2TM 613 A2TMZO = A2ZO 614 A2TMZC = A2ZC 615 A2TMYO = A2YO 616 A2TMYC = A2YC 617 A2TM = A2TOT 618 END IF 619 IF ( IPRRSP .GT. 30 ) THEN 620 WRITE(LUPRI,'(A)') ' Case 1 for A2 term' 621 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 622 END IF 623C 624C Case 2: 625C ======= 6262000 IF (MZCONF(ISYMV) .EQ. 0 ) GO TO 4000 627C 628C Multiply the factor of one half into the operator matrix 629C for evaluation of Case 2 and 3. 630C 631C I don't like it, but we scale the vector for cases 2 and 3 632C 633C CALL DSCAL(NORBT*NORBT,DH,OPMAT,1) 634 CALL DSCAL(MZCONF(ISYMV),DH,VEC,1) 635 CALL DSCAL(MZCONF(ISYMV),DH,VEC(1+MZVAR(ISYMV)),1) 636C 637 ISPIN = IOPSPI 638C 639C Make the gradient 640C 641 OVLAP = D0 642 ISYMDN = 1 643 INTSYM = IOPSYM 644 IKLVL = 0 645 ISYMR = ISYMV 646 NZYVEC = MZYVAR(ISYMV) 647 NZCVEC = MZCONF(ISYMV) 648 LORB = .FALSE. 649 ISYMST = MULD2H(IGRSYM,IREFSY) 650 IF(ISYMST .EQ. IREFSY ) THEN 651 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 652 ELSE 653 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 654 END IF 655 LREFST = .FALSE. 656 CALL HSO2CR(IGRSYM,IGRSPI,A2TRS,VEC,NZYVEC,NZCVEC,ISYMV,ISYMDN, 657 * XINDX,OVLAP,UDV,PV,OPMAT,WRK(KFREE),LFREE,KZYVR, 658 * LCON,LORB,LREFST,IDAX,INTSYM,IOPSPI,0, 659 * IKLVL,ZYMB,ISYMB,ZYMC,ISYMC,MJWOP) 660C CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,A2TRS, 661C * VEC,NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,OPMAT, 662C * XINDX,WRK,LWRK,LORB,LCON,LREFST) 663C 664 CALL DSCAL(MZCONF(ISYMV),D2,VEC,1) 665 CALL DSCAL(MZCONF(ISYMV),D2,VEC(1+MZVAR(ISYMV)),1) 666 IF (A2TEST) THEN 667 A2ZO = DDOT(KZA2O,GRVEC(KZO),1,A2TRS(KYO),1) 668 A2ZC = DDOT(KZA2C,GRVEC(KZC),1,A2TRS(KYC),1) 669 A2YO = DDOT(KZA2O,GRVEC(KYO),1,A2TRS(KZO),1) 670 A2YC = DDOT(KZA2C,GRVEC(KYC),1,A2TRS(KZC),1) 671 A2TOT = A2ZO+A2ZC+A2YO+A2YC 672 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 ZO :',A2ZO - A2TMZO 673 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 ZC :',A2ZC - A2TMZC 674 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 YO :',A2YO - A2TMYO 675 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 YC :',A2YC - A2TMYC 676 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 :',A2TOT - A2TM 677 A2TMZO = A2ZO 678 A2TMZC = A2ZC 679 A2TMYO = A2YO 680 A2TMYC = A2YC 681 A2TM = A2TOT 682 END IF 683 IF ( IPRRSP .GT. 100 ) THEN 684 WRITE(LUPRI,'(A)') ' Case 2 for A2 term' 685 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 686 END IF 687C 688C Skip expectation value of imaginary operator 689C 690C Case 3: 691C ======= 692C Do Sj(2) <0 |op| 0> 693C 6944000 CONTINUE 695 IF (A2TEST) THEN 696 WRITE(LUPRI,'(/A,F24.8)')' A2TEST Final result: ZO', A2ZO 697 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Final result: ZC', A2ZC 698 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Final result: YO', A2YO 699 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Final result: YC', A2YC 700 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Final result: ', A2TOT 701 END IF 702 IF ( IPRRSP .GT. 150 ) THEN 703 WRITE(LUPRI,'(//A)') ' Gradient before swapping in A2DRV' 704 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 705 END IF 706C 707C Swap the final result to conform with the notation for A[2] 708C 709 CALL DSWAP(MZVAR(IGRSYM),A2TRS,1, 710 * A2TRS(MZVAR(IGRSYM)+1),1) 711C 712 IF ( IPRRSP .GT. 100 ) THEN 713 WRITE(LUPRI,'(//A)') ' Gradient after swapping in A2DRV' 714 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 715 END IF 716 CALL DSCAL(KZYVR,HSOFAC,A2TRS,1) 717 CALL QEXIT('A2HSO') 718C 719 RETURN 720 END 721C /* Deck hso2cr */ 722 SUBROUTINE HSO2CR(IGRSYM,IGRSPI,ETRS,VEC,NZYVEC,NZCVEC,ISYMV, 723 * ISYMDN, 724 * XINDX,OVLAP,DEN1,DEN2,FI,WRK,LWRK,KZYVA, 725 * LCON,LORB,LREFST,IDAX,INTSYM,IOPSPI,ISPINV, 726 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2,MJWOP) 727C 728C Layout the core for RSP2GR 729C 730#include "implicit.h" 731C 732#include "maxorb.h" 733#include "infdim.h" 734#include "inforb.h" 735#include "infvar.h" 736#include "infrsp.h" 737#include "wrkrsp.h" 738#include "qrinf.h" 739#include "infpri.h" 740#include "infspi.h" 741C 742 DIMENSION WRK(*) 743 DIMENSION ETRS(KZYVA), MJWOP(2,MAXWOP,8) 744 DIMENSION XINDX(*) 745 DIMENSION DEN1(*), DEN2(N2ASHX*N2ASHX,*) 746 DIMENSION FI(*) 747 DIMENSION VEC(NZYVEC) 748 DIMENSION ZYMAT1(*), ZYMAT2(*) 749C 750 LOGICAL LCON, LORB, LREFST, HSORUN 751C 752C Layout of WRK: 753C We keep the H2AX array at the beginning, in order to 754C release extra workspace in the gradient routine after processing 755C the orbital part of the linear transformation. 756C 757 KH2AX = 1 758 IF (LCON) THEN 759 IF (DIROIT) THEN 760 KFA = KH2AX + N2ASHX * N2ASHX * 2 761 ELSE 762 KFA = KH2AX + N2ASHX * N2ASHX 763 END IF 764 ELSE 765 KFA = KH2AX 766 END IF 767 KQA = KFA + NORBT * NORBT 768 KQB = KQA + NORBT * NASHDI 769 KPVD = KQB + NORBT * NASHDI 770 KH2 = KPVD + N2ASHX * N2ASHX 771 KH2X = KH2 + N2ORBX 772 KWRKO = KH2X + N2ORBX 773C 774C Get free workspace for orbital part of linear transformation 775C 776 LWRKO = LWRK - KWRKO 777 IF (LWRKO.LT.0) CALL ERRWRK('HSO2CR 1',KWRKO-1,LWRK) 778 CALL DZERO(WRK,KWRKO) 779C 780C Get space that can be used in configuration part 781C We need the arrays H2XAC and FI, so keep them intact 782C 783 KWRKC = KFA 784 LWRKC = LWRK - KWRKC 785 IF (LWRKC.LT.0) CALL ERRWRK('HSO2CR 2',KWRKC-1,LWRK) 786C 787 HSORUN = .TRUE. 788 CALL RSP2GR (IGRSYM,IGRSPI,ETRS,VEC,NZYVEC,NZCVEC,ISYMV,ISYMDN, 789 * FI,WRK(KFA),WRK(KQA),WRK(KQB),WRK(KH2),WRK(KH2X), 790 * OVLAP,DEN1,DEN2,WRK(KPVD),WRK(KH2AX),XINDX, 791 * WRK(KWRKO),LWRKO,WRK(KWRKC),LWRKC,KZYVA, 792 * LCON,LORB,LREFST,IDAX,INTSYM,IOPSPI,ISPINV,IDUM, 793 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2,DUM,IDUM,HSORUN, 794 * DUM,DUM,MJWOP) 795C 796 RETURN 797 END 798C /* Deck hsofxd */ 799 SUBROUTINE HSOFXD (FI,FA,QA,QB,H2AX, 800 * IDAX,INTSYM, ISYMDN,DEN1,DEN2, 801 * PVDEN,H2,H2X,WRK,KFREE,LFREE, 802 * LCON,LORB,IPRFX,IGRSPI,IOPSPI,ISPINV, 803 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2) 804C 805C 806C 920206 Olav Vahtras 807C 808C This routine computes the FX matrices, that is, the Fock 809C type matrices FI, FA, QA, and QB with transformed ("X") 810C integrals. In addition, If LCON true then the H2AX matrix 811C with transformed integrals is extracted for the CI routines. 812C 813C The operator is of the form 814C HSO = h(r,s)T(r,s) + (rs|t^u)(e(+,-) + 2e(-,+))(rstu) 815C 816C 817#include "implicit.h" 818C 819 PARAMETER ( D0 = 0.0D0, DH = 0.50D0, D2 = 2.0D0 ) 820 PARAMETER ( IPLUS = 0, IMINUS = 1 ) 821C 822C Used from common blocks: 823C ? 824C 825#include "maxorb.h" 826#include "maxash.h" 827#include "priunit.h" 828#include "inforb.h" 829#include "infind.h" 830#include "infdim.h" 831#include "infrsp.h" 832#include "infpri.h" 833#include "orbtypdef.h" 834#include "inftap.h" 835#include "rspprp.h" 836#include "infhyp.h" 837#include "infspi.h" 838#include "trhso.h" 839#include "cbgetdis.h" 840#include "infhso.h" 841#include "infdis.h" 842C 843C 844 DIMENSION FI(NORBT,NORBT), FA(NORBT,NORBT) 845 DIMENSION QA(NORBT,NASHDI),QB(NORBT,NASHDI) 846 DIMENSION H2(NORBT,NORBT), H2X(NORBT,NORBT) 847 DIMENSION H2AX(N2ASHX*N2ASHX,*) 848 DIMENSION DEN1(NASHDI,NASHDI), DEN2(*), PVDEN(NASHDI,NASHDI) 849 DIMENSION ZYMAT1(NORBT,NORBT),ZYMAT2(NORBT,NORBT) 850 DIMENSION WRK(*) 851C 852 DIMENSION NEEDED(-4:6) 853C 854 LOGICAL LCON, LORB 855 LOGICAL ACALL 856 COMMON /CB_HSO_ACALL/ACALL 857C 858C Put up the structure for reading the (transformed) integrals: 859C IEND < 0 flags the completeness of the distributions 860C 861 CALL QENTER('HSOFXD') 862 CALL DZERO(H2,N2ORBX) 863 864 IF (IPRFX.GT.200) THEN 865 WRITE(LUPRI,'(/A)') 'ENTER HSOFXD' 866 WRITE(LUPRI,'(A,I5)') 'IDAX =',IDAX 867 WRITE(LUPRI,'(A,I5)') 'INTSYM =',INTSYM 868 WRITE(LUPRI,'(A,I5)') 'ISYMDN =',ISYMDN 869 WRITE(LUPRI,'(A,L1)') 'LCON =',LCON 870 WRITE(LUPRI,'(A,L2)') 'LORB =',LORB 871 WRITE(LUPRI,'(A,I5)') 'IPRFX =',IPRFX 872 WRITE(LUPRI,'(A,I5)') 'IGRSPI =',IGRSPI 873 WRITE(LUPRI,'(A,I5)') 'IOPSPI =',IOPSPI 874 WRITE(LUPRI,'(A,I5)') 'ISPINV =',ISPINV 875 WRITE(LUPRI,'(A,I5)') 'IKLVL =',IKLVL 876 WRITE(LUPRI,'(A,I5)') 'ISYM1 =',ISYM1 877 WRITE(LUPRI,'(A,I5)') 'ISYM2 =',ISYM2 878 WRITE(LUPRI,'(A)') 'One-electron FI' 879 CALL OUTPUT(FI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 880 IF (LORB) THEN 881 WRITE(LUPRI,'(A)') 'One-electron density' 882 CALL OUTPUT(DEN1,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 883 WRITE(LUPRI,'(A)') 'Two-electron density 1' 884 CALL PRIAC2(DEN2,NASHT,LUPRI) 885 WRITE(LUPRI,'(A)') 'Two-electron density 2' 886 CALL PRIAC2(DEN2(1+N2ASHX*N2ASHX),NASHT,LUPRI) 887 END IF 888 END IF 889 IF (IOPSPI.NE.1) THEN 890 WRITE(LUPRI,'(/A)') 'HSOFXD: WRONG VALUE OF IOPSPI' 891 WRITE(LUPRI,'(A,I5)') 'IOPSPI =', IOPSPI 892 CALL QTRACE(0) 893 CALL QUIT('HSOFXD: WRONG VALUE OF IOPSPI') 894 END IF 895 IPERCD = -1 896 IF (NASHT.GT.0) THEN 897 CALL MEMGET('REAL',KH2XAC,N2ASHX,WRK,KFREE,LFREE) 898 CALL DZERO(WRK(KH2XAC),N2ASHX) 899 IF (LORB) THEN 900 LSCR = MAX(NORBT,N2ASHX) 901 ELSE 902 LSCR = 0 903 END IF 904 CALL MEMGET('REAL',KSCR,LSCR,WRK,KFREE,LFREE) 905 ELSE 906 CALL MEMGET('REAL',KH2XAC,0,WRK,KFREE,LFREE) 907 CALL MEMGET('REAL',KSCR ,0,WRK,KFREE,LFREE) 908 END IF 909 910 IF (.NOT.DOSO2) GOTO 2000 911 912 NEEDED(-4:6) = 0 913 IF (IKLVL.EQ.1) THEN 914 NEEDED(1:6) = 1 915 ELSE 916 NEEDED(1:5) = 1 917 END IF 918 IDIST = 0 9191000 CALL NXTHSO(IC,ID,H2,NEEDED,WRK,KFREE,LFREE,IDIST) 920 IF ( IDIST .LT. 0) GO TO 2000 921 IF (IC.EQ.ID) GO TO 1000 922C 923C Transpose the operator for A[2] calcualtion 924C 925 IF (ACALL) THEN 926 ITMP = IC 927 IC = ID 928 ID = ITMP 929 END IF 930C 931C Determine type of distribution 932C 933 ICTYP = IOBTYP(IC) 934 IDTYP = IOBTYP(ID) 935 ICDTYP = IDBTYP(ICTYP,IDTYP) 936C 937C Determine symmetry of distribution 938C 939 ICSYM = ISMO( IC ) 940 IDSYM = ISMO( ID ) 941 ICDSYM = MULD2H(ICSYM,IDSYM) 942C 943C 944C 945 IF ( IPRFX .GT. 300) THEN 946 WRITE(LUPRI,'(//A)') 'Characteristics of distribution:' 947 WRITE(LUPRI,'(A)') '================================' 948 WRITE(LUPRI,'(A,2I5)')'IC ID = ', IC,ID 949 WRITE(LUPRI,'(3A)') 'TYP = ', COBTYP(ICTYP),COBTYP(IDTYP) 950 WRITE(LUPRI,'(A,2I5)')'Symmetry= ', ICSYM, IDSYM 951 WRITE(LUPRI,'(A)') '================================' 952C 953 CALL HEADER('Two-electron integrals from this distribution',3) 954 CALL OUTPUT(H2,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 955 END IF 956C 957 IF (IDAX .EQ. LUMHSO) THEN 958C 959C Regular MO-integrals 960C 961 IF (IKLVL.EQ.0) THEN 962C 963C ( |^)e(+,-) 964C 965 IF (LORB) THEN 966 IPDA = IPDENS(IGRSPI,1) 967 IPDB = IPDENS(0,MULSP(IGRSPI,1)) 968 END IF 969 CALL GETAC1(H2,WRK(KH2XAC)) 970 CALL FQDIS1(INTSYM,IC,ID,ICDTYP,ICDSYM,ICSYM,IDSYM, 971 * FI,FA,QA,QB,H2,WRK(KH2XAC),H2AX, 972 * DEN1,DEN2(IPDA),DEN2(IPDB),ISYMDN,PVDEN, 973 * LORB,LCON, 974 * IGRSPI,IPLUS,IMINUS,IPERCD) 975C 976C 2( |^)e(-,+) 977C 978 IF (LORB) THEN 979 IPDA = IPDENS(MULSP(IGRSPI,1),0) 980 IPDB = IPDENS(1,IGRSPI) 981 END IF 982 CALL DSCAL(N2ORBX,D2,H2,1) 983 CALL DSCAL(N2ASHX,D2,WRK(KH2XAC),1) 984 CALL FQDIS1(INTSYM,IC,ID,ICDTYP,ICDSYM,ICSYM,IDSYM, 985 * FI,FA,QA,QB,H2,WRK(KH2XAC),H2AX, 986 * DEN1,DEN2(IPDA),DEN2(IPDB),ISYMDN,PVDEN, 987 * LORB,.FALSE., 988 * IGRSPI,IMINUS,IPLUS,IPERCD) 989 ELSE IF (IKLVL.EQ.1) THEN 990 IH2SYM = MULD2H(ICDSYM,INTSYM) 991C 992C Transform left (~| ) 993C 994 CALL DZERO(H2X,N2ORBX) 995 CALL OITH1(ISYM1,ZYMAT1,H2,H2X,IH2SYM) 996 INTSY1 = MULD2H(INTSYM,ISYM1) 997 IF (IPRFX.GT.300) THEN 998 WRITE(LUPRI,'(/A,I5)')'Left transform H2 of symmetry', 999 * IH2SYM 1000 CALL OUTPUT(H2,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1001 WRITE(LUPRI,'(/A,I5)') 'with ZYMAT1 of symmetry',ISYM1 1002 CALL OUTPUT(ZYMAT1,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1003 WRITE(LUPRI,'(/A)') 'to H2X ' 1004 CALL OUTPUT(H2X,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1005 END IF 1006 IF (LORB) THEN 1007 IPDA = IPDENS(MULSP(IGRSPI,ISPINV),1) 1008 IPDB = IPDENS(ISPINV,MULSP(IGRSPI,1)) 1009 END IF 1010 CALL GETAC1(H2X,WRK(KH2XAC)) 1011C 1012C (~|^)e(~+,-) 1013C 1014 CALL FQDIS1(INTSY1,IC,ID,ICDTYP,ICDSYM,ICSYM,IDSYM, 1015 * FI,FA,QA,QB,H2X,WRK(KH2XAC),H2AX, 1016 * DEN1,DEN2(IPDA),DEN2(IPDB),ISYMDN,PVDEN, 1017 * LORB,LCON, 1018 * IGRSPI,ISPINV,IMINUS,IPERCD) 1019 IF (LORB) THEN 1020 IPDA = IPDENS(0,0) 1021 IPDB = IPDENS(MULSP(ISPINV,1),IGRSPI) 1022 END IF 1023 CALL DSCAL(N2ORBX,D2,H2X,1) 1024 CALL DSCAL(N2ASHX,D2,WRK(KH2XAC),1) 1025C 1026C 2(~|^)e(~-,+) 1027C 1028 CALL FQDIS1(INTSY1,IC,ID,ICDTYP,ICDSYM,ICSYM,IDSYM, 1029 * FI,FA,QA,QB,H2X,WRK(KH2XAC),H2AX, 1030 * DEN1,DEN2(IPDA),DEN2(IPDB),ISYMDN,PVDEN, 1031 * LORB,.FALSE., 1032 * IGRSPI,MULSP(ISPINV,IMINUS),IPLUS,IPERCD) 1033C 1034C Transform right ( |~^) 1035C 1036 IF (LORB) THEN 1037 IPDA = IPDENS(IGRSPI,MULSP(ISPINV,1)) 1038 IPDB = IPDENS(0,0) 1039 END IF 1040 CALL GETAC1(H2,WRK(KH2XAC)) 1041C 1042C ( |~^)e(+,~-) 1043C 1044 CALL FQDIS2(INTSYM,IC,ID,ICDTYP,ICDSYM,ICSYM,IDSYM, 1045 * FI,FA,QA,QB,H2,WRK(KH2XAC),H2AX(1,2), 1046 * ZYMAT1,ISYM1, 1047 * DEN1,DEN2(IPDA),DEN2(IPDB),ISYMDN,PVDEN, 1048 * LORB,LCON, 1049 * IGRSPI,IPLUS,MULSP(ISPINV,IMINUS), 1050 * IPERCD,WRK(KSCR)) 1051 IF (LORB) THEN 1052 IPDA = IPDENS(MULSP(IGRSPI,1),ISPINV) 1053 IPDB = IPDENS(1,MULSP(IGRSPI,ISPINV)) 1054 END IF 1055 CALL DSCAL(N2ORBX,D2,H2,1) 1056 CALL DSCAL(N2ASHX,D2,WRK(KH2XAC),1) 1057C 1058C 2 ( |~^)e(-,~+) 1059C 1060 CALL FQDIS2(INTSYM,IC,ID,ICDTYP,ICDSYM,ICSYM,IDSYM, 1061 * FI,FA,QA,QB,H2,WRK(KH2XAC),H2AX(1,2), 1062 * ZYMAT1,ISYM1, 1063 * DEN1,DEN2(IPDA),DEN2(IPDB),ISYMDN,PVDEN, 1064 * LORB,.FALSE., 1065 * IGRSPI,IMINUS,ISPINV,IPERCD,WRK(KSCR)) 1066 ELSE 1067 WRITE(LUPRI,'(/A)') ' WRONG VALUE OF IKLVL IN HSOFXD' 1068 WRITE(LUPRI,'(/A,I5)') ' IKLVL =',IKLVL 1069 CALL QUIT(' WRONG VALUE OF IKLVL IN HSOFXD') 1070 END IF 1071 ELSE 1072 WRITE(LUPRI,'(/A,I5)') 'HSOFXD: WRONG VALUE OF UNIT IDAX',IDAX 1073 WRITE(LUPRI,'(/A,2I5)') ' IDAX .NE. LUMHSO',IDAX,LUMHSO 1074 CALL QUIT('HSOFXD: WRONG VALUE OF UNIT IDAX') 1075 END IF 1076C 1077C 1078C 1079 IF( IPRFX .GT. 300 ) THEN 1080 WRITE(LUPRI,'(/A)') ' PARTIAL MATRICES IN HSOFXD' 1081 WRITE(LUPRI,'(A)') ' ===========================' 1082 WRITE(LUPRI,'(//A)') ' Inactive fock matrix' 1083 CALL OUTPUT(FI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1084 WRITE(LUPRI,'(//A)') ' Active fock matrix' 1085 CALL OUTPUT(FA,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1086 WRITE(LUPRI,'(//A)') ' Qa matrix' 1087 CALL OUTPUT(QA,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1088 WRITE(LUPRI,'(//A)') ' Qb matrix' 1089 CALL OUTPUT(QB,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1090 WRITE(LUPRI,'(/A)') ' Active two-electron integrals' 1091 CALL OUTPUT(H2AX,1,N2ASHX,1,N2ASHX,N2ASHX,N2ASHX,1,LUPRI) 1092 IF (IKLVL.EQ.1) THEN 1093 WRITE(LUPRI,'(/A)') ' Active two-electron integrals' 1094 CALL OUTPUT(H2AX(1,2),1,N2ASHX,1,N2ASHX,N2ASHX,N2ASHX,1, 1095 & LUPRI) 1096 END IF 1097 END IF 1098C 1099C 1100C We have now completed processing this load, so get the next 1101C 1102 GO TO 1000 11032000 CONTINUE 1104C 1105C Account for 1/2 in definition of two-electron integrals 1106C 1107C and for the fact that we calculated 2 times FA 1108C 1109 CALL DSCAL(NORBT*NORBT,DH,FA,1) 1110C 1111C Set DISTYP for the cases we may have ci-gradients 1112C 1113 IF (IKLVL.EQ.0) THEN 1114 INFDIS(1) = 23 1115 INFDIS(2) = 24 1116 END IF 1117 IF (IKLVL.EQ.1) THEN 1118 INFDIS(1) = 19 1119 INFDIS(2) = 20 1120 END IF 1121 IF (IPRFX.GT.200) THEN 1122 WRITE(LUPRI,'(/A)') ' Completed matrices in HSOFXD' 1123 WRITE(LUPRI,'(/A)') ' Inactive Fock matrix' 1124 CALL OUTPUT(FI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1125 IF (LORB) THEN 1126 WRITE(LUPRI,'(/A)') ' Active Fock matrix' 1127 CALL OUTPUT(FA,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1128 WRITE(LUPRI,'(/A)') ' QA matrix' 1129 CALL OUTPUT(QA,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1130 WRITE(LUPRI,'(/A)') ' QB matrix' 1131 CALL OUTPUT(QB,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1132 END IF 1133 IF (LCON)THEN 1134 WRITE(LUPRI,'(/A)') ' Active two-electron integrals' 1135 CALL PRIAC2(H2AX,NASHT,LUPRI) 1136 IF (IKLVL.EQ.1) THEN 1137 WRITE(LUPRI,'(/A)') ' Active two-electron integrals' 1138 CALL PRIAC2(H2AX(1,2),NASHT,LUPRI) 1139 END IF 1140 END IF 1141 END IF 1142C 1143 CALL QEXIT('HSOFXD') 1144 RETURN 1145 END 1146C /* Deck ipdens */ 1147 FUNCTION IPDENS(ISPIN1,ISPIN2) 1148C 1149C Get a pointer to the density <e(ISPIN1,ISPIN2)> 1150C The vector INFDEN(4) should be reset when a two-electron density 1151C is constructed. In a triplet calculation we normally have two densities 1152C stored after one another, e.g. ( <e(+,+)> <e(-,-)> ) 1153C Referring to these densities INFDEN should be set to (0,0,1,1) 1154C 1155#include "infden.h" 1156#include "inforb.h" 1157 CALL QENTER('IPDENS') 1158 IF (ISPIN1.EQ.INFDEN(1) .AND. ISPIN2.EQ.INFDEN(2)) THEN 1159 IPDENS = 1 1160 ELSE IF (ISPIN1.EQ.INFDEN(3) .AND. ISPIN2.EQ.INFDEN(4)) THEN 1161 IPDENS = N2ASHX*N2ASHX + 1 1162 ELSE 1163 CALL QTRACE(0) 1164 WRITE(0,'(/A,4I3)') 'INFDEN =',(INFDEN(I), I=1,4) 1165 CALL QUIT('IPDENS: REQUIRED DENSITIES ARE NOT STORED') 1166 END IF 1167 CALL QEXIT('IPDENS') 1168 RETURN 1169 END 1170C /* Deck ipset */ 1171 SUBROUTINE IPSET(I,J,K,L) 1172#include "infden.h" 1173 INFDEN(1) = I 1174 INFDEN(2) = J 1175 INFDEN(3) = K 1176 INFDEN(4) = L 1177 RETURN 1178 END 1179C 1180C /* Deck x2hsoao */ 1181C 1182 SUBROUTINE X2HSOAO( 1183 & OPLBL,CMO,MJWOP, 1184 & LVEC,VECSYM,VECSPIN, 1185 & VEC, 1186 & LGR,GRSYM,GRSPIN, 1187 & GR, 1188 & WORK,LWORK 1189 & ) 1190C 1191C Built gradient of one-index transformed spin-orbit operator 1192C from AO integrals (single determinant) 1193C 1194 IMPLICIT NONE 1195C 1196C Input 1197C 1198 CHARACTER*8 OPLBL 1199 REAL*8 CMO(*) 1200 INTEGER MJWOP(*) 1201 INTEGER LVEC, VECSYM, VECSPIN 1202 REAL*8 VEC(*) 1203 INTEGER LGR, GRSYM, GRSPIN 1204C 1205C Output 1206C 1207 REAL*8 GR 1208C 1209C Work 1210C 1211 INTEGER LWORK 1212 REAL*8 WORK(LWORK) 1213#include "inforb.h" 1214#include "priunit.h" 1215#include "infrsp.h" 1216#include "infhso.h" 1217#include "maxorb.h" 1218#include "infinp.h" 1219#include "aovec.h" 1220#include "dftcom.h" 1221#include "dummy.h" 1222C 1223C Local 1224C 1225 INTEGER KFREE,LFREE,KAO,KMO,KS,KD,KD1,KD2,KF,KF1,KF2 1226 LOGICAL TRIPLET(2) 1227 REAL*8 D1 1228 PARAMETER (D1 = 1.0D0) 1229C 1230 INTEGER IPRINT, NDMAT 1231 INTEGER KJSTRS , KNPRIM , KNCONT , KIORBS , KJORBS , NPAO , KKORBS 1232 REAL*8 TIMEND,TIMSTR,SO1WAL,SO1CPU 1233 REAL*8 SO2WAL,SO2CPU,SOCPU,SOWAL, HFXSAV 1234 INTEGER I2TYP 1235 INTEGER IPRTWO,IPRNTA,IPRNTB,IPRNTC,IPRNTD 1236 INTEGER IRNTYP , NUMDIS, IFCTYP(2), ISYMDM(2) 1237 LOGICAL TKTIME,RTNTWO 1238C 1239 IPRINT=MAX(IPRRSP,IPRHSO) 1240 NDMAT=2 1241 IF (NASHT.GT.0) 1242 & CALL QUIT('X2HSOAO not implemented for open shells') 1243 1244 CALL QENTER('X2HSOAO') 1245 KFREE=1 1246 LFREE=LWORK 1247 CALL MEMGET('REAL',KAO,N2BASX,WORK,KFREE,LFREE) 1248 CALL MEMGET('REAL',KMO,N2ORBX,WORK,KFREE,LFREE) 1249C 1250C Unpack vectors 1251C 1252 CALL DZERO(WORK(KMO),N2ORBX) 1253 CALL GTZYMT(1,VEC,LVEC,VECSYM,WORK(KMO),MJWOP) 1254C 1255C Ao-tranformed kappa 1256C 1257 CALL DZERO(WORK(KAO),N2BASX) 1258 CALL MOTOAO(WORK(KMO),WORK(KAO),CMO,VECSYM,WORK(KFREE),LFREE) 1259 CALL MEMREL('MOTOAO',WORK,KMO,KMO,KFREE,LFREE) 1260 IF (IPRHSO .gt. 100 .or. HSODEBUG) THEN 1261 CALL MAOPRI(WORK(KAO),'X2HSOAO:kappa(ao)') 1262 END IF 1263C 1264C Get overlap 1265C 1266 CALL MEMGET('REAL',KS,N2BASX,WORK,KFREE,LFREE) 1267 CALL GET_H1(WORK(KS),'OVERLAP',WORK(KFREE),LFREE) 1268C 1269C Allocate densities and Fock 1270C 1271 CALL MEMGET('REAL',KF,2*N2BASX,WORK,KFREE,LFREE) 1272 CALL MEMGET('REAL',KD,2*N2BASX,WORK,KFREE,LFREE) 1273 CALL DZERO(WORK(KF),2*N2BASX) 1274 CALL DZERO(WORK(KD),2*N2BASX) 1275 KD1=KD 1276 KD2=KD+N2BASX 1277 KF1=KF 1278 KF2=KF+N2BASX 1279C 1280C Transform density 1281C 1282 CALL GTDMSO(WORK(KFREE),CMO,WORK(KD1),WORK(KFREE),WORK(KFREE)) 1283 CALL D_K(NBAST,WORK(KAO),WORK(KD1),WORK(KD2),WORK(KS), 1284 & WORK(KFREE),LFREE) 1285 IF (IPRHSO .gt. 100 .or. HSODEBUG) THEN 1286 CALL MAOPRI(WORK(KD1),'X2HSOAO:Density 1') 1287 CALL MAOPRI(WORK(KD2),'X2HSOAO:Density 2') 1288 END IF 1289C 1290C Build fock 1291C 1292 IF (DIRFCK) THEN 1293 IF (OPLBL(1:1).EQ.'X') IFCTYP(1)=1 1294 IF (OPLBL(1:1).EQ.'Y') IFCTYP(1)=2 1295 IF (OPLBL(1:1).EQ.'Z') IFCTYP(1)=3 1296 IFCTYP(2) = IFCTYP(1) 1297 IF (GRSPIN.NE.VECSPIN) IFCTYP(1) = -IFCTYP(1) 1298 IF (GRSPIN.EQ.1) IFCTYP(2) = -IFCTYP(2) 1299C 1300C Transform density to AO basis 1301C 1302 CALL DSOTAO(WORK(KD1),WORK(KF1),NBAST,0,IPRINT) 1303 CALL DSOTAO(WORK(KD2),WORK(KF2),NBAST,VECSYM-1,IPRINT) 1304 CALL DCOPY(NDMAT*N2BASX,WORK(KF),1,WORK(KD),1) 1305C 1306C Setup for TWOINT 1307C 1308 NPAO = MXSHEL*MXAOVC 1309 CALL MEMGET('INTE',KJSTRS,NPAO*2,WORK,KFREE,LFREE) 1310 CALL MEMGET('INTE',KNPRIM,NPAO*2,WORK,KFREE,LFREE) 1311 CALL MEMGET('INTE',KNCONT,NPAO*2,WORK,KFREE,LFREE) 1312 CALL MEMGET('INTE',KIORBS,NPAO ,WORK,KFREE,LFREE) 1313 CALL MEMGET('INTE',KJORBS,NPAO ,WORK,KFREE,LFREE) 1314 CALL MEMGET('INTE',KKORBS,NPAO ,WORK,KFREE,LFREE) 1315 CALL PAOVEC(WORK(KJSTRS),WORK(KNPRIM),WORK(KNCONT), 1316 & WORK(KIORBS),WORK(KJORBS),WORK(KKORBS),0, 1317 & .FALSE.,IPRRSP) 1318 CALL MEMREL('HERINT.PAOVEC',WORK,1,KJORBS,KFREE,LFREE) 1319 CALL TIMER('START ',TIMSTR,TIMEND) 1320 CALL GETTIM(SO1CPU,SO1WAL) 1321 I2TYP = 0 1322 IRNTYP = -20 1323 IPRTWO = 0 1324 TKTIME = .FALSE. 1325 CALL DZERO(WORK(KF),NDMAT*N2BASX) 1326C Always full exchange in spin-orbit integrals: 1327 HFXSAV=HFXFAC 1328 HFXFAC=D1 1329 CALL TWOINT(WORK(KFREE),LFREE,WORK(KFREE), 1330 & WORK(KF),WORK(KD),NDMAT, 1331 & IDUMMY, IFCTYP, 1332 & DUMMY,IDUMMY,NUMDIS,1,IRNTYP,0,0,.TRUE.,.TRUE., 1333 & .FALSE.,TKTIME,IPRTWO,IPRNTA,IPRNTB,IPRNTC, 1334 & IPRNTD,RTNTWO,IDUMMY,I2TYP,WORK(KJSTRS), 1335 & WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 1336 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 1337 & .FALSE.,.false.) 1338 CALL MEMREL('HSOFCKAO.TWOINT',WORK,KF,KJSTRS,KFREE,LFREE) 1339 HFXFAC=HFXSAV 1340 IF (IPRHSO .gt. 100 .or. HSODEBUG) THEN 1341 CALL MAOPRI(WORK(KF1),'X2HSOAO(TWOINT):Fock 1') 1342 CALL MAOPRI(WORK(KF2),'X2HSOAO(TWOINT):Fock 2') 1343 END IF 1344C 1345C Transform to SO 1346C 1347 ISYMDM(1)=0 1348 ISYMDM(2)=VECSYM-1 1349 CALL SKLFCK(WORK(KF),VDUMMY,WORK(KFREE),LFREE, 1350 & IPRTWO,.FALSE., 1351 & .FALSE.,.FALSE.,.FALSE.,.TRUE.,IDUMMY,.TRUE.,NDMAT, 1352 & ISYMDM,IFCTYP,IDUMMY,.TRUE.) 1353C 1354 CALL MEMCHK('HSOFCKAO.SKLFCK',WORK,KF) 1355 CALL GETTIM(SO2CPU,SO2WAL) 1356 SOCPU=SO2CPU-SO1CPU 1357 SOWAL=SO2WAL-SO1WAL 1358 CALL TIMER('TWOINT',TIMSTR,TIMEND) 1359 CALL FLSHFO(LUPRI) 1360 WRITE(LUPRI,'(2(/A),2(/A,F5.0,A))') 1361 & ' Two-electron spin-orbit integrals', 1362 & ' =================================', 1363 & ' Spin-orbit 2-electron CPU time ',SOCPU,' seconds', 1364 & ' Spin-orbit 2-electron wall time ',SOWAL,' seconds' 1365 ELSE 1366 TRIPLET(1) = GRSPIN.NE.VECSPIN 1367 TRIPLET(2) = GRSPIN.EQ.1 1368 CALL GET_FSO_AO(OPLBL,TRIPLET,WORK(KF),WORK(KD),2) 1369 END IF 1370C 1371C One-electron 1372C 1373 IF (OPLBL(2:2).EQ.' ') THEN 1374 CALL GET_P(OPLBL(1:1)//'1'//OPLBL(3:8),WORK(KD1)) 1375 CALL DAXPY(N2BASX,D1,WORK(KD1),1,WORK(KF1),1) 1376 END IF 1377 IF (IPRHSO .gt. 100 .or. HSODEBUG) THEN 1378 CALL MAOPRI(WORK(KF1),'X2HSOAO:Fock 1') 1379 CALL MAOPRI(WORK(KF2),'X2HSOAO:Fock 2') 1380 END IF 1381C 1382C Final transform (adds [k,f1] to f2 1383C 1384 CALL F_K(NBAST,WORK(KAO),WORK(KF1),WORK(KF2),WORK(KS), 1385 & WORK(KFREE),LFREE) 1386 IF (IPRHSO .gt. 100 .or. HSODEBUG) THEN 1387 CALL MAOPRI(WORK(KF2),'X2HSOAO:Final fock') 1388 END IF 1389C 1390C Mo basis (add to one-electron in input) 1391C 1392 CALL DZERO(WORK(KF1),N2ORBX) 1393 CALL AOTOMO(WORK(KF2),WORK(KF1),CMO,GRSYM,WORK(KFREE),LFREE) 1394 IF (IPRHSO .gt. 100 .or. HSODEBUG) THEN 1395 CALL MAOPRI(WORK(KF1),'X2HSOAO:Final fock (mo)') 1396 END IF 1397C 1398C Final gradient 1399C 1400 TRPLET=GRSPIN.EQ.VECSPIN 1401 CALL ORBEX(GRSYM,1,LGR, 1402 & WORK(KF1),WORK(KFREE),WORK(KFREE),WORK(KFREE), 1403 & GR,D1,WORK(KFREE),MJWOP) 1404 IF (IPRHSO .gt. 100 .or. HSODEBUG) THEN 1405 CALL OUTPUT(GR,1,LGR/2,1,2,LGR/2,2,1,LUPRI) 1406 END IF 1407 CALL MEMREL('X2HSOAO',WORK,1,1,KFREE,LFREE) 1408 CALL QEXIT('X2HSOAO') 1409 END 1410C 1411C /* Deck get_p */ 1412C 1413 SUBROUTINE GET_P(LABEL,P) 1414C 1415C Simplified PRPGET 1416C 1417#include "implicit.h" 1418 CHARACTER*8 LABEL 1419 DIMENSION P(*) 1420C 1421C External 1422C 1423#include "inforb.h" 1424#include "inftap.h" 1425#include "dummy.h" 1426 LOGICAL FNDLB2, CLOSED 1427 EXTERNAL FNDLB2 1428C 1429C Local 1430C 1431 DIMENSION TMP(N2BASX) 1432 CHARACTER*8 RTNLBL(2) 1433C 1434 CALL QENTER('GET_P') 1435C 1436C Leave AOPROPER in the same state (open/closed) 1437C 1438 CLOSED=LUPROP.LT.0 1439 IF (CLOSED) THEN 1440 CALL GPOPEN(LUPROP,'AOPROPER','OLD',' ','UNFORMATTED',IDUMMY, 1441 & .FALSE.) 1442 END IF 1443 REWIND LUPROP 1444 IF ( FNDLB2(LABEL,RTNLBL,LUPROP)) THEN 1445 IF (RTNLBL(2).EQ.'SQUARE ') THEN 1446 CALL READT(LUPROP,N2BASX,P) 1447 ELSE IF (RTNLBL(2).EQ.'SYMMETRI') THEN 1448 CALL READT(LUPROP,NNBASX,TMP) 1449 CALL DSPTSI(NBAST,TMP,P) 1450 ELSE IF (RTNLBL(2).EQ.'ANTISYMM') THEN 1451 CALL READT(LUPROP,NNBASX,TMP) 1452 CALL DAPTGE(NBAST,TMP,P) 1453 ELSE 1454 CALL QUIT('Error: No antisymmetry label on LUPROP') 1455 END IF 1456 ELSE 1457 CALL QUIT('GET_P: '//LABEL//' not found on LUPROP') 1458 END IF 1459 IF (CLOSED) CALL GPCLOSE(LUPROP,'KEEP') 1460 CALL QEXIT('GET_P') 1461 END 1462