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#ifdef REVLOG 20=========================================================================== 21Revision 1.2 2000/05/24 19:09:06 hjj 22inserted Dalton header 23some changes for triplet response with CSF 24=========================================================================== 25#endif 26C 27 SUBROUTINE E3SOL(VECA, VEC1, VEC2,ETRS,XINDX,ZYM1,ZYM2, 28 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 29 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP,ISYRLM) 30C 31C Purpose: 32C Outer driver routine for solvent contribution 33C to E[3] times two vectors. 34C 35#include "implicit.h" 36C 37#include "maxorb.h" 38#include "priunit.h" 39#include "inforb.h" 40#include "infdim.h" 41#include "infinp.h" 42#include "infvar.h" 43#include "infrsp.h" 44#include "infpri.h" 45#include "rspprp.h" 46#include "infcr.h" 47C 48 DIMENSION ETRS(KZYVR),XINDX(*) 49 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 50 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 51 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 52 DIMENSION MJWOP(2,MAXWOP,8) 53 INTEGER DBGFLG(10) 54 55C 56 NSIM = 1 57 KFREE = 1 58 CALL MEMGET('REAL',KTA ,N2ORBX,WRK,KFREE,LFREE) 59 CALL MEMGET('REAL',KTB ,N2ORBX,WRK,KFREE,LFREE) 60 CALL MEMGET('REAL',KTB1,N2ORBX,WRK,KFREE,LFREE) 61 CALL MEMGET('REAL',KTB2,N2ORBX,WRK,KFREE,LFREE) 62 CALL MEMGET('REAL',KTC1,N2ORBX,WRK,KFREE,LFREE) 63 CALL MEMGET('REAL',KTC2,N2ORBX,WRK,KFREE,LFREE) 64 CALL MEMGET('REAL',KTD1,N2ORBX,WRK,KFREE,LFREE) 65 CALL MEMGET('REAL',KTD2,N2ORBX,WRK,KFREE,LFREE) 66 CALL MEMGET('REAL',KTE ,N2ORBX,WRK,KFREE,LFREE) 67C 68 CALL DZERO(WRK(KTA),N2ORBX) 69 CALL DZERO(WRK(KTB),N2ORBX) 70 CALL DZERO(WRK(KTB1),N2ORBX) 71 CALL DZERO(WRK(KTB2),N2ORBX) 72 CALL DZERO(WRK(KTC1),N2ORBX) 73 CALL DZERO(WRK(KTC2),N2ORBX) 74 CALL DZERO(WRK(KTD1),N2ORBX) 75 CALL DZERO(WRK(KTD2),N2ORBX) 76 CALL DZERO(WRK(KTE),N2ORBX) 77C 78 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,ZYM1,MJWOP) 79 CALL GTZYMT(NSIM,VEC2,KZYV2,ISYMV2,ZYM2,MJWOP) 80C 81C DBGFLG initialization 82C 1 2 3 4 5 6 7 8 9 10 83C DATA DBGFLG/-1,-2,-3,-4,-5,-6,-7,-8,-9,-10/ 84 DATA DBGFLG/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10/ 85C 86C VECA is only available if E3TEST is set 87C 88 VAL = DDOT(KZYVR,VECA,1,ETRS,1) 89C 90 CALL TCASE1(VECA, VEC1, VEC2,WRK(KTA),WRK(KTB), 91 * ETRS,XINDX,ZYM1,ZYM2,ISYRLM, 92 * DEN1,UDV,WRK(KFREE),LFREE, 93 * KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 94 * DBGFLG) 95C 96 CALL TCASE2(VECA,VEC1,VEC2,WRK(KTC1), 97 * ETRS,XINDX,ZYM1,ZYM2,ISYRLM, 98 * DEN1,UDV,WRK(KFREE),LFREE, 99 * KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 100 * DBGFLG) 101C 102 CALL TCASE2(VECA,VEC2,VEC1,WRK(KTC2), 103 * ETRS,XINDX,ZYM2,ZYM1,ISYRLM, 104 * DEN1,UDV,WRK(KFREE),LFREE, 105 * KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1,CMO,MJWOP, 106 * DBGFLG) 107C 108 CALL TCASE3(VECA,VEC1,VEC2,WRK(KTB),WRK(KTB1),WRK(KTC1),WRK(KTD1), 109 * ETRS,XINDX,ZYM1,ZYM2,ISYRLM, 110 * DEN1,UDV,WRK(KFREE),LFREE, 111 * KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 112 * DBGFLG) 113C 114 CALL TCASE3(VECA,VEC2,VEC1,WRK(KTB),WRK(KTB2),WRK(KTC2),WRK(KTD2), 115 * ETRS,XINDX,ZYM2,ZYM1,ISYRLM, 116 * DEN1,UDV,WRK(KFREE),LFREE, 117 * KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1,CMO,MJWOP, 118 * DBGFLG) 119C 120 CALL TCASE4(VECA, VEC1, VEC2,WRK(KTA),WRK(KTB1), 121 * WRK(KTB2),WRK(KTC1),WRK(KTC2),WRK(KTE), 122 * ETRS,XINDX,ZYM1,ZYM2,ISYRLM, 123 * DEN1,UDV,WRK(KFREE),LFREE, 124 * KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 125 * DBGFLG) 126C 127 IF (CRCAL .OR. E3TEST) THEN 128 WRITE (LUPRI,'(A,F20.12)') 'Total contribution from E3SOL:', 129 * DDOT(KZYVR,VECA,1,ETRS,1) - VAL 130 END IF 131C 132 RETURN 133 END 134 SUBROUTINE TOPGET(ZYM1,ZYM2,ZYM3,IKLVL,OVLAP,ISYMDN, 135 * ISYMV1,ISYMV2,ISYMV3,TOP,DEN1,EPS,NUCFLG, 136 * ISYRLM,CMO,WRK,LWRK,IPRLVL,PRSTR,ADDFLG) 137C 138C Output: 139C 140C TOP = sum(l,m) f(l) <L| T(l,m)(k1,k2,..) |R> TE(l,m) 141C 142C Input: 143C 144C OVLAP = <L|R>, DEN1 = <L|...|R> 145C NUCFLG indicates if T(l,m) = TN(l,m) or T(l,m) = TE(l,m) 146C IKLVL is the number of times T(l,m) is to be one-index tranformed 147C 148#include "implicit.h" 149#include "dummy.h" 150C 151#include "maxorb.h" 152#include "priunit.h" 153#include "infdim.h" 154#include "inforb.h" 155#include "infpri.h" 156#include "inftap.h" 157#include "infinp.h" 158#include "infrsp.h" 159C 160 LOGICAL NUCFLG 161 INTEGER ADDFLG 162 CHARACTER*(*)PRSTR 163C 164 DIMENSION DEN1(NASHDI,NASHDI) 165 DIMENSION WRK(*), CMO(*) 166 DIMENSION ZYM1(NORBT,NORBT),ZYM2(NORBT,NORBT),ZYM3(NORBT,NORBT) 167 DIMENSION TOP(N2ORBX) 168 DIMENSION ISYRLM(*) 169C 170 NSIM = 1 171C 172 IF (IKLVL.EQ.0) ISYM = 1 173 IF (IKLVL.EQ.1) ISYM = ISYMV1 174 IF (IKLVL.EQ.2) ISYM = MULD2H(ISYMV1,ISYMV2) 175 IF (IKLVL.EQ.3) ISYM = MULD2H(MULD2H(ISYMV1,ISYMV2),ISYMV3) 176 ISYMT = MULD2H(ISYM,ISYMDN) 177C 178 KFREE = 1 179 CALL MEMGET('REAL',KTNLM ,NLMSOL,WRK,KFREE,LWRK) 180 CALL MEMGET('REAL',KRLMAO,NNBASX,WRK,KFREE,LWRK) 181 CALL MEMGET('REAL',KUCMO ,NORBT*NBAST,WRK,KFREE,LWRK) 182 CALL MEMGET('REAL',KTELM ,N2ORBX,WRK,KFREE,LWRK) 183 CALL MEMGET('REAL',KTLMA ,N2ORBX,WRK,KFREE,LWRK) 184 CALL MEMGET('REAL',KTLMB ,N2ORBX,WRK,KFREE,LWRK) 185Clf 186 CALL MEMGET('REAL',KTLMC ,N2ORBX,WRK,KFREE,LWRK) 187Clf 188 CALL MEMGET('REAL',KFLVEC,NLMSOL,WRK,KFREE,LWRK) 189Clf 190 CALL DZERO(WRK(KTLMC),N2ORBX) 191 192C 193C Unpack symmetry blocked CMO 194C 195 CALL UPKCMO(CMO,WRK(KUCMO)) 196C 197C Calculate f(l) factors. 198C 199 CALL SOLFL(WRK(KFLVEC),EPS,RSOL,LSOLMX) 200C 201C Read nuclear contributions TN(l,m). 202C 203 IF (LUSOL.LE.0) 204 &CALL GPOPEN(LUSOL,FNSOL,'OLD',' ','UNFORMATTED',IDUMMY,.FALSE.) 205 REWIND LUSOL 206 CALL MOLLAB('SOLVRLM ',LUSOL,LUERR) 207 READ (LUSOL) 208 CALL READT(LUSOL,NLMSOL,WRK(KTNLM)) 209C 210C Loop over l,m expansion. 211C 212 LM = 0 213 DO 520 L = 0,LSOLMX 214 READ (LUSOL) L1,(ISYRLM(M),M=1,2*L+1) 215 DO 500 M = -L,L 216 LM = LM + 1 217 IF (ISYRLM(L+M+1) .NE. ISYMT) THEN 218 READ (LUSOL) 219 GO TO 500 220 END IF 221C 222C Read R(l,m) in ao basis, transform to mo basis and unpack. 223C Electronic contribution TE(l,m). 224C 225 CALL READT(LUSOL,NNBASX,WRK(KRLMAO)) 226 CALL UTHU(WRK(KRLMAO),WRK(KTELM),WRK(KUCMO),WRK(KFREE), 227 * NBAST,NORBT) 228 CALL DCOPY(NNORBX,WRK(KTELM),1,WRK(KTLMA),1) 229 CALL DSPTSI(NORBT,WRK(KTLMA),WRK(KTELM)) 230C 231C One-index transform TE(l,m) IKLVL times. 232C The result will be in WRK(KTLMA) and of symmetry ISYM. 233C (ISYM should equal ISYMDN.) 234C 235 CALL DCOPY(N2ORBX,WRK(KTELM),1,WRK(KTLMA),1) 236 ISYM = ISYMT 237 IF (IKLVL.GE.1) THEN 238 CALL DZERO(WRK(KTLMA),N2ORBX) 239 CALL OITH1(ISYMV1,ZYM1,WRK(KTELM),WRK(KTLMA),ISYM) 240 ISYM = MULD2H(ISYM,ISYMV1) 241 END IF 242 IF (IKLVL.GE.2) THEN 243 CALL DZERO(WRK(KTLMB),N2ORBX) 244 CALL OITH1(ISYMV2,ZYM2,WRK(KTLMA),WRK(KTLMB),ISYM) 245 CALL DCOPY(N2ORBX,WRK(KTLMB),1,WRK(KTLMA),1) 246 ISYM = MULD2H(ISYM,ISYMV2) 247 END IF 248 IF (IKLVL.GE.3) THEN 249 CALL DZERO(WRK(KTLMA),N2ORBX) 250 CALL OITH1(ISYMV3,ZYM3,WRK(KTLMB),WRK(KTLMA),ISYM) 251 ISYM = MULD2H(ISYM,ISYMV3) 252 END IF 253C 254C Add the contribution from TE(l,m) or TN(l,m) to the effective operator. 255C 256 IF (NUCFLG) THEN 257 FACT = WRK((KTNLM-1)+LM) 258 ELSE 259 CALL MELONE(WRK(KTLMA),ISYM,DEN1,OVLAP,FACT,200,'TOPGET') 260 END IF 261 FACT = WRK((KFLVEC-1)+LM)*FACT 262 CALL DAXPY(N2ORBX,FACT,WRK(KTELM),1,WRK(KTLMC),1) 263 IF (ADDFLG.GT.0) THEN 264 CALL DAXPY(N2ORBX,FACT,WRK(KTELM),1,TOP,1) 265 END IF 266C 267 500 CONTINUE 268 520 CONTINUE 269C 270 CALL GPCLOSE(LUSOL,'KEEP') 271C 272 IF (IPRRSP.GE.IPRLVL) THEN 273 WRITE(LUPRI,'(/3A,2D22.14)') 'Norm of TOPGET in ', PRSTR, 274 * ' : ',DNRM2(N2ORBX,wrk(ktlmc),1),DNRM2(N2ORBX,top,1) 275 END IF 276C 277 RETURN 278C 279 END 280 SUBROUTINE TCASE1(VECA, VEC1, VEC2,TA,TB, 281 * ETRS,XINDX,ZYM1,ZYM2,ISYRLM, 282 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 283 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 284 * DBGFLG) 285#include "implicit.h" 286C 287 PARAMETER ( D1 = 1.0D0, DM1 = -1.0D0, D2 = 2.0D0, DM2 = -2.0D0 ) 288C 289#include "maxorb.h" 290#include "infdim.h" 291#include "inforb.h" 292#include "wrkrsp.h" 293#include "infrsp.h" 294#include "infpri.h" 295#include "infvar.h" 296#include "qrinf.h" 297#include "infspi.h" 298#include "infden.h" 299#include "infinp.h" 300C 301 DIMENSION ETRS(KZYVR),XINDX(*),ISYRLM(2*LSOLMX+1) 302 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 303 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 304 DIMENSION TA(NORBT,NORBT),TB(NORBT,NORBT) 305 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 306 DIMENSION MJWOP(2,MAXWOP,8) 307 INTEGER DBGFLG(10) 308C 309 LOGICAL TDM, LREF, NORHO2 310C 311C Initialise variables 312C 313 JSPIN = 0 314 TDM = .TRUE. 315 KFREE = 1 316 NORHO2 = .TRUE. 317 NSIM = 1 318C 319C TA = 2*sum(l,m) f(l) <0| TE(l,m) |0> TE(l,m) 320C 321 CALL DZERO(TA,N2ORBX) 322 CALL TOPGET(DUMMY,DUMMY,DUMMY,0,D1,1, 323 * IDUMMY,IDUMMY,IDUMMY,TA,UDV,EPSOL,.FALSE., 324 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TA', 325 $ DBGFLG(1)) 326 CALL DSCAL(N2ORBX,D2,TA,1) 327C 328C TB = 2*sum(l,m) f(l)*( <0|TE(l,m)|0> - Tn(l,m) )*TE(l,m) 329C 330 IF (INERSI) THEN 331 EPS = EPSTAT 332 ELSE 333 EPS = EPSOL 334 END IF 335 CALL DZERO(TB,N2ORBX) 336 CALL TOPGET(DUMMY,DUMMY,DUMMY,0,DUMMY,1, 337 * IDUMMY,IDUMMY,IDUMMY,TB,DUMMY,EPS,.TRUE., 338 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TB', 339 $ DBGFLG(2)) 340clf 341 CALL DSCAL(N2ORBX,DM1,TB,1) 342 CALL TOPGET(DUMMY,DUMMY,DUMMY,0,D1,1, 343 * IDUMMY,IDUMMY,IDUMMY,TB,UDV,EPS,.FALSE., 344 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TB', 345 $ DBGFLG(3)) 346clf 347 CALL DSCAL(N2ORBX,D2,TB,1) 348C 349 IF (MZCONF(ISYMV1) .EQ. 0 .OR. MZCONF(ISYMV2) .EQ. 0) RETURN 350C 351C / <01L| [qj,TB] |02R> + <02L| [qj,TB] |01R> \ 352C | 0 | 353C | <01L| [qj+,TB] |02R> + <02L| [qj+,TB] |01R> | 354C \ 0 / 355C 356C Construct <01L|..|02R> + <02L|..|01R> density 357C 358 ILSYM = MULD2H(IREFSY,ISYMV1) 359 IRSYM = MULD2H(IREFSY,ISYMV2) 360 NCL = MZCONF(ISYMV1) 361 NCR = MZCONF(ISYMV2) 362 KZVARL = MZYVAR(ISYMV1) 363 KZVARR = MZYVAR(ISYMV2) 364 LREF = .FALSE. 365 ISYMDN = MULD2H(ILSYM,IRSYM) 366 CALL DZERO(DEN1,NASHT*NASHT) 367 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 368 * VEC1,VEC2,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM,NORHO2, 369 * XINDX,WRK,KFREE,LFREE,LREF) 370C 371C Make the gradient 372C 373 IF ( MZWOPT(IGRSYM) .GT. 0 ) THEN 374 CALL ORBSX(NSIM,IGRSYM,KZYVR,ETRS,TB,OVLAP,ISYMDN, 375 * DEN1,MJWOP,WRK(KFREE),LFREE) 376 END IF 377C 378 CALL PRIRES(ETRS,VECA,IGRSYM,'TCASE1') 379C 380 RETURN 381 END 382 SUBROUTINE TCASE2(VECA, VEC1, VEC2,TC1, 383 * ETRS,XINDX,ZYM1,ZYM2,ISYRLM, 384 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 385 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 386 * DBGFLG) 387#include "implicit.h" 388C 389 PARAMETER ( D1 = 1.0D0, D2 = 2.0D0 ) 390C 391#include "maxorb.h" 392#include "infdim.h" 393#include "inforb.h" 394#include "wrkrsp.h" 395#include "infrsp.h" 396#include "infpri.h" 397#include "infvar.h" 398#include "qrinf.h" 399#include "infspi.h" 400#include "infden.h" 401#include "infinp.h" 402C 403 DIMENSION ETRS(KZYVR),XINDX(*),ISYRLM(2*LSOLMX+1) 404 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 405 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 406 DIMENSION TC1(NORBT,NORBT) 407 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 408 DIMENSION MJWOP(2,MAXWOP,8) 409 INTEGER DBGFLG(10) 410C 411 LOGICAL TDM, LREF, NORHO2 412C 413C Initialise variables 414C 415 JSPIN = 0 416 TDM = .TRUE. 417 IPRONE = 200 418 KFREE = 1 419 NORHO2 = .TRUE. 420 NSIM = 1 421C 422 CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE) 423 CALL GETREF(WRK(KCREF),MZCONF(1)) 424C 425C TC1 = 2*sum(l,m) f(l) <0| TE(l,m)(k1) |0> TE(l,m) + ... 426C 427 CALL DZERO(TC1,N2ORBX) 428 IF (MZWOPT(ISYMV1).GT.0) THEN 429 CALL TOPGET(ZYM1,DUMMY,DUMMY,1,D1,1, 430 * ISYMV1,IDUMMY,IDUMMY,TC1,UDV,EPSOL,.FALSE., 431 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TC1 cont1', 432 $ DBGFLG(4)) 433 END IF 434C 435C ... + 2*sum(l,m) f(l) ( <01L| TE(l,m) |0> + <0| TE(l,m) |01R> ) TE(l,m) 436C 437 IF (MZCONF(ISYMV1).GT.0) THEN 438C 439C Construct the density matrix <01L|..|0> + <0|..|01R> 440C 441 ILSYM = IREFSY 442 IRSYM = MULD2H(IREFSY,ISYMV1) 443 NCL = MZCONF(1) 444 NCR = MZCONF(ISYMV1) 445 KZVARL = MZCONF(1) 446 KZVARR = MZYVAR(ISYMV1) 447 LREF = .TRUE. 448 ISYMDN = MULD2H(ILSYM,IRSYM) 449 CALL DZERO(DEN1,NASHT*NASHT) 450 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 451 * WRK(KCREF),VEC1,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM, 452 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 453C 454 CALL TOPGET(DUMMY,DUMMY,DUMMY,0,OVLAP,ISYMDN, 455 * IDUMMY,IDUMMY,IDUMMY,TC1,DEN1,EPSOL,.FALSE., 456 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TC1 cont2', 457 $ DBGFLG(5)) 458 END IF 459C 460 CALL DSCAL(N2ORBX,D2,TC1,1) 461C 462 IF (MZCONF(ISYMV2).LE.0) RETURN 463C 464C / 0 \ 465C | Sj(2) | * <0| TC1 |0> 466C | 0 | 467C \ Sj(2)' / 468C 469 IF (IGRSYM.EQ.ISYMV2) THEN 470 OVLAP = D1 471 CALL MELONE(TC1,1,UDV,OVLAP,FACT,IPRONE,'FACT in TCASE2') 472 NZCONF = MZCONF(IGRSYM) 473 NZVAR = MZVAR(IGRSYM) 474 CALL DAXPY(NZCONF,FACT,VEC2,1,ETRS,1) 475 CALL DAXPY(NZCONF,FACT,VEC2(NZVAR+1),1,ETRS(NZVAR+1),1) 476 END IF 477C 478 CALL PRIRES(ETRS,VECA,IGRSYM,'TCASE2') 479C 480 RETURN 481 END 482 SUBROUTINE TCASE3(VECA, VEC1, VEC2,TB,TB1,TC1,TD1, 483 * ETRS,XINDX,ZYM1,ZYM2,ISYRLM, 484 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 485 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 486 * DBGFLG) 487#include "implicit.h" 488C 489 PARAMETER ( D1 = 1.0D0 ) 490C 491#include "maxorb.h" 492#include "infdim.h" 493#include "inforb.h" 494#include "wrkrsp.h" 495#include "infrsp.h" 496#include "infpri.h" 497#include "infvar.h" 498#include "qrinf.h" 499#include "infspi.h" 500#include "infden.h" 501#include "infinp.h" 502C 503 DIMENSION ETRS(KZYVR),XINDX(*),ISYRLM(2*LSOLMX+1) 504 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 505 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 506 DIMENSION TB(NORBT,NORBT),TB1(NORBT,NORBT) 507 DIMENSION TC1(NORBT,NORBT),TD1(NORBT,NORBT) 508 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 509 DIMENSION MJWOP(2,MAXWOP,8) 510 INTEGER DBGFLG(10) 511C 512 LOGICAL LCON, LORB 513 LOGICAL TDM, LREF, NORHO2 514C 515C Initialise variables 516C 517 JSPIN = 0 518 TDM = .TRUE. 519 KFREE = 1 520 IPRONE = -1 521 NORHO2 = .TRUE. 522 NSIM = 1 523C 524C TD1 = TB1 + TC1, TB1 = TB(k1) 525C 526 CALL DZERO(TB1,N2ORBX) 527 CALL DZERO(TD1,N2ORBX) 528 CALL OITH1(ISYMV1,ZYM1,TB,TB1,1) 529 CALL DAXPY(N2ORBX,D1,TB1,1,TD1,1) 530 CALL DAXPY(N2ORBX,D1,TC1,1,TD1,1) 531C 532 IF (MZCONF(ISYMV2).LE.0) RETURN 533C 534 CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE) 535 CALL GETREF(WRK(KCREF),MZCONF(1)) 536C 537C / <0| [qj,TD1] |02R> + <02L| [qj,TD1] |0> \ 538C | <j| TD1 |02R> | 539C | <0| [qj+,TD1] |02R> + <02L| [qj+,TD1] |0> | 540C \ -<02L| TD1 |j> / 541C 542C Construct the density matrix <02L|..|0> + <0|..|02R> 543C 544 ILSYM = IREFSY 545 IRSYM = MULD2H(IREFSY,ISYMV2) 546 NCL = MZCONF(1) 547 NCR = MZCONF(ISYMV2) 548 KZVARL = MZCONF(1) 549 KZVARR = MZYVAR(ISYMV2) 550 LREF = .TRUE. 551 ISYMDN = MULD2H(ILSYM,IRSYM) 552 CALL DZERO(DEN1,NASHT*NASHT) 553 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 554 * WRK(KCREF),VEC2,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM, 555 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 556C 557C Make the gradient 558C 559 ISYMST = MULD2H(IGRSYM,IREFSY) 560 IF ( ISYMST .EQ. IREFSY ) THEN 561 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 562 ELSE 563 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 564 END IF 565 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 566 LREF = .FALSE. 567 NZYVEC = MZYVAR(ISYMV2) 568 NZCVEC = MZCONF(ISYMV2) 569 CALL RSP1GR(NSIM,KZYVR,IDUMMY,JSPIN,IGRSYM,JSPIN,ISYMV2,ETRS, 570 * VEC2,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,TD1, 571 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 572C 573 CALL PRIRES(ETRS,VECA,IGRSYM,'TCASE3') 574C 575 RETURN 576 END 577 SUBROUTINE TCASE4(VECA, VEC1, VEC2,TA,TB1,TB2,TC1,TC2,TE, 578 * ETRS,XINDX,ZYM1,ZYM2,ISYRLM, 579 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 580 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 581 * DBGFLG) 582#include "implicit.h" 583C 584 PARAMETER ( D1 = 1.0D0, D2 = 2.0D0, DH = 0.5D0 ) 585C 586#include "maxorb.h" 587#include "infdim.h" 588#include "inforb.h" 589#include "wrkrsp.h" 590#include "infrsp.h" 591#include "infpri.h" 592#include "infvar.h" 593#include "qrinf.h" 594#include "infspi.h" 595#include "infden.h" 596#include "infinp.h" 597C 598 DIMENSION ETRS(KZYVR),XINDX(*),ISYRLM(2*LSOLMX+1) 599 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 600 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 601 DIMENSION TA(NORBT,NORBT),TE(NORBT,NORBT) 602 DIMENSION TB1(NORBT,NORBT),TB2(NORBT,NORBT) 603 DIMENSION TC1(NORBT,NORBT),TC2(NORBT,NORBT) 604 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 605 DIMENSION MJWOP(2,MAXWOP,8) 606 INTEGER DBGFLG(10) 607C 608 LOGICAL LCON, LORB 609 LOGICAL TDM, LREF, NORHO2 610C 611C Initialise variables 612C 613 JSPIN = 0 614 TDM = .TRUE. 615 KFREE = 1 616 IPRONE = 100 617 NORHO2 = .TRUE. 618 NSIM = 1 619C 620 CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE) 621 CALL GETREF(WRK(KCREF),MZCONF(1)) 622C 623C TE = 1/2 * TB1(k2) + 1/2 * TB2(k1) + TC1(k2) + TC2(k1) + ... 624C 625 CALL DZERO(TE,N2ORBX) 626 CALL OITH1(ISYMV2,ZYM2,TB1,TE,ISYMV1) 627 CALL OITH1(ISYMV1,ZYM1,TB2,TE,ISYMV2) 628 CALL DSCAL(N2ORBX,DH,TE,1) 629 CALL OITH1(ISYMV2,ZYM2,TC1,TE,ISYMV1) 630 CALL OITH1(ISYMV1,ZYM1,TC2,TE,ISYMV2) 631C 632C ... + ( S(1)S(2)' + S(2)S(1)' ) * TA + ... 633C 634 IF ((ISYMV1.EQ.ISYMV2) .AND. (MZCONF(ISYMV1).GT.0)) THEN 635 NZCONF = MZCONF(ISYMV1) 636 NZVAR = MZVAR(ISYMV1) 637 FACT = DDOT(NZCONF,VEC1,1,VEC2(NZVAR+1),1) + 638 * DDOT(NZCONF,VEC2,1,VEC1(NZVAR+1),1) 639 CALL DAXPY(N2ORBX,FACT,TA,1,TE,1) 640 END IF 641C 642C ... + sum(l,m) f(l) <0| TE(l,m)(k1,k2) |0> TE(l,m) 643C + sum(l,m) f(l) <0| TE(l,m)(k2,k1) |0> TE(l,m) + ... 644C 645 IF (MZWOPT(ISYMV1).GT.0 .AND. MZWOPT(ISYMV2).GT.0) THEN 646 CALL TOPGET(ZYM1,ZYM2,DUMMY,2,D1,1, 647 * ISYMV1,ISYMV2,IDUMMY,TE,UDV,EPSOL,.FALSE., 648 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TE cont1a', 649 $ DBGFLG(6)) 650 CALL TOPGET(ZYM2,ZYM1,DUMMY,2,D1,1, 651 * ISYMV2,ISYMV1,IDUMMY,TE,UDV,EPSOL,.FALSE., 652 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TE cont1b', 653 $ DBGFLG(7)) 654 END IF 655C 656C ... + 2*sum(l,m) f(l) ( <01L| TE(l,m)(k2) |0> + 657C <0| TE(l,m)(k2) |01R> ) TE(l,m) + ... 658C 659C Put the factor two into one of the vectors. 660C 661 CALL DSCAL(KZYV1,D2,VEC1,1) 662 CALL DSCAL(NORBT*NORBT,D2,ZYM1,1) 663C 664 IF (MZCONF(ISYMV1).GT.0 .AND. MZWOPT(ISYMV2).GT.0) THEN 665C 666C Construct the density matrix <01L|..|0> + <0|..|01R> 667C 668 ILSYM = IREFSY 669 IRSYM = MULD2H(IREFSY,ISYMV1) 670 NCL = MZCONF(1) 671 NCR = MZCONF(ISYMV1) 672 KZVARL = MZCONF(1) 673 KZVARR = MZYVAR(ISYMV1) 674 LREF = .TRUE. 675 ISYMDN = MULD2H(ILSYM,IRSYM) 676 CALL DZERO(DEN1,NASHT*NASHT) 677 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 678 * WRK(KCREF),VEC1,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM, 679 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 680C 681 CALL TOPGET(ZYM2,DUMMY,DUMMY,1,OVLAP,ISYMDN, 682 * ISYMV2,IDUMMY,IDUMMY,TE,DEN1,EPSOL,.FALSE., 683 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TE cont2a', 684 $ DBGFLG(8)) 685 END IF 686C 687C ... + 2*sum(l,m) f(l) ( <02L| TE(l,m)(k1) |0> + 688C <0| TE(l,m)(k1) |02R> ) TE(l,m) + ... 689C 690C The factor two is already included in one of the vectors. 691C 692 IF (MZCONF(ISYMV2).GT.0 .AND. MZWOPT(ISYMV1).GT.0) THEN 693C 694C Construct the density matrix <02L|..|0> + <0|..|02R> 695C 696 ILSYM = IREFSY 697 IRSYM = MULD2H(IREFSY,ISYMV2) 698 NCL = MZCONF(1) 699 NCR = MZCONF(ISYMV2) 700 KZVARL = MZCONF(1) 701 KZVARR = MZYVAR(ISYMV2) 702 LREF = .TRUE. 703 ISYMDN = MULD2H(ILSYM,IRSYM) 704 CALL DZERO(DEN1,NASHT*NASHT) 705 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 706 * WRK(KCREF),VEC2,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM, 707 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 708C 709 CALL TOPGET(ZYM1,DUMMY,DUMMY,1,OVLAP,ISYMDN, 710 * ISYMV1,IDUMMY,IDUMMY,TE,DEN1,EPSOL,.FALSE., 711 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TE cont2b', 712 $ DBGFLG(9)) 713 END IF 714C 715C ... + 2*sum(l,m) f(l) ( <01L| TE(l,m) |02R> + 716C <02L| TE(l,m) |01R> ) TE(l,m) + ... 717C 718C The factor two is already included in one of the vectors. 719C 720 IF (MZCONF(ISYMV1) .GT. 0 .AND. MZCONF(ISYMV2) .GT. 0) THEN 721C 722C Construct <01L|..|02R> + <02L|..|01R> density 723C 724 ILSYM = MULD2H(IREFSY,ISYMV1) 725 IRSYM = MULD2H(IREFSY,ISYMV2) 726 NCL = MZCONF(ISYMV1) 727 NCR = MZCONF(ISYMV2) 728 KZVARL = MZYVAR(ISYMV1) 729 KZVARR = MZYVAR(ISYMV2) 730 LREF = .FALSE. 731 ISYMDN = MULD2H(ILSYM,IRSYM) 732 CALL DZERO(DEN1,NASHT*NASHT) 733 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 734 * VEC1,VEC2,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM,NORHO2, 735 * XINDX,WRK,KFREE,LFREE,LREF) 736C 737 CALL TOPGET(DUMMY,DUMMY,DUMMY,0,OVLAP,ISYMDN, 738 * IDUMMY,IDUMMY,IDUMMY,TE,DEN1,EPSOL,.FALSE., 739 * ISYRLM,CMO,WRK(KFREE),LFREE,100,'TE cont3', 740 $ DBGFLG(10)) 741 END IF 742C 743C / <0| [qj ,TE] |0> \ 744C | <j| TE |0> | 745C | <0| [qj+,TE] |0> | 746C \ -<0| TE |j> / 747C 748 ISYMDN = 1 749 OVLAP = D1 750 ISYMV = IREFSY 751 ISYMST = MULD2H(IGRSYM,IREFSY) 752 IF ( ISYMST .EQ. IREFSY ) THEN 753 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 754 ELSE 755 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 756 END IF 757 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 758 LREF = .TRUE. 759 NZYVEC = MZCONF(1) 760 NZCVEC = MZCONF(1) 761 CALL RSP1GR(NSIM,KZYVR,IDUMMY,JSPIN,IGRSYM,JSPIN,ISYMV,ETRS, 762 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,TE, 763 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 764C 765C Restore the vector 766C 767 CALL DSCAL(KZYV1,DH,VEC1,1) 768C 769 CALL PRIRES(ETRS,VECA,IGRSYM,'TCASE4') 770C 771 RETURN 772 END 773 SUBROUTINE C3SOL(VECA,VEC1,VEC2,ETRS,XINDX,ZYM1,ZYM2, 774 & UDV,WRK,LWRK,KZYVR,KZYV1,KZYV2, 775 & IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP,ISYRLM) 776C 777C Purpose: 778C Memeory efficient routine for computing solvent contribution 779C to E[3] times two vectors. Replaces E3SOL in SCF calculations. 780C 781#include "implicit.h" 782#include "dummy.h" 783C 784 PARAMETER ( D0=0.0D0, D1=1.0D0 ) 785C 786#include "maxorb.h" 787#include "inforb.h" 788#include "infdim.h" 789#include "infinp.h" 790#include "infvar.h" 791#include "infrsp.h" 792#include "infpri.h" 793#include "rspprp.h" 794#include "infcr.h" 795#include "inftap.h" 796#include "qrinf.h" 797C 798 DIMENSION ETRS(KZYVR),XINDX(*) 799 DIMENSION UDV(NASHDI,NASHDI) 800 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 801 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 802 DIMENSION MJWOP(2,MAXWOP,8),ISYRLM(2*LSOLMX+1) 803C 804 LOGICAL LCON, LORB, LREF 805C 806 NSIM = 1 807 KFREE = 1 808 LFREE = LWRK 809 CALL MEMGET('REAL',KTRES ,N2ORBX,WRK,KFREE,LFREE) 810 CALL MEMGET('REAL',KTELM ,N2ORBX,WRK,KFREE,LFREE) 811 CALL MEMGET('REAL',KTLMA ,N2ORBX,WRK,KFREE,LFREE) 812 CALL MEMGET('REAL',KTLMB ,N2ORBX,WRK,KFREE,LFREE) 813 CALL MEMGET('REAL',KFLST ,NLMSOL,WRK,KFREE,LFREE) 814 CALL MEMGET('REAL',KFLOP ,NLMSOL,WRK,KFREE,LFREE) 815 CALL MEMGET('REAL',KTNLM ,NLMSOL,WRK,KFREE,LFREE) 816 CALL MEMGET('REAL',KRLMAO,NNBASX,WRK,KFREE,LFREE) 817 CALL MEMGET('REAL',KUCMO ,NORBT*NBAST,WRK,KFREE,LFREE) 818 CALL MEMGET('REAL',KCREF ,NCREF,WRK,KFREE,LFREE) 819C 820C Get the reference state 821C 822 CALL GETREF(WRK(KCREF),MZCONF(1)) 823C 824C Zero the final effective operator 825C 826 CALL DZERO(WRK(KTRES),N2ORBX) 827C 828C Unpack the response vectors 829C 830 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,ZYM1,MJWOP) 831 CALL GTZYMT(NSIM,VEC2,KZYV2,ISYMV2,ZYM2,MJWOP) 832C 833C Unpack symmetry blocked CMO 834C 835 CALL UPKCMO(CMO,WRK(KUCMO)) 836C 837C Calculate f(l) factors. 838C 839 CALL SOLFL(WRK(KFLST),EPSTAT,RSOL,LSOLMX) 840 CALL SOLFL(WRK(KFLOP),EPSOL,RSOL,LSOLMX) 841C 842C Read nuclear contributions TN(l,m). 843C 844 IF (LUSOL .LE. 0) CALL GPOPEN(LUSOL,FNSOL,'UNKNOWN',' ', 845 & 'UNFORMATTED',IDUMMY,.FALSE.) 846 REWIND LUSOL 847 CALL MOLLAB('SOLVRLM ',LUSOL,LUERR) 848 READ (LUSOL) 849 CALL READT(LUSOL,NLMSOL,WRK(KTNLM)) 850C 851C Loop over l,m expansion. 852C 853 LM = 0 854 DO 520 L = 0,LSOLMX 855 READ (LUSOL) L1,(ISYRLM(M),M=1,2*L+1) 856 DO 500 M = -L,L 857 LM = LM + 1 858 ISYMT = ISYRLM(L+M+1) 859 IF (ISYMT.NE.1 .AND. ISYMT.NE.ISYMV1 .AND. 860 & ISYMT.NE.ISYMV2 .AND. ISYMT.NE.MULD2H(ISYMV1,ISYMV2)) THEN 861 READ (LUSOL) 862 GO TO 500 863 END IF 864C 865C Read R(l,m) in ao basis, transform to mo basis and unpack. 866C Electronic contribution TE(l,m). 867C 868 CALL READT(LUSOL,NNBASX,WRK(KRLMAO)) 869 CALL UTHU(WRK(KRLMAO),WRK(KTELM),WRK(KUCMO),WRK(KFREE), 870 & NBAST,NORBT) 871 CALL DCOPY(NNORBX,WRK(KTELM),1,WRK(KTLMA),1) 872 CALL DSPTSI(NORBT,WRK(KTLMA),WRK(KTELM)) 873C 874C Create the effective operator: 875C 876C TRES = sum(l,m)[ W(k1,k2) + A1(k2) + A12 ] 877C 878C W(k1,k2) = g(l)*( F1 + F2 )*TELM(k1,k2) 879C A1(k2) = g(l)*F3*TELM(k2) 880C A12 = g(l)*F4*TELM 881C 882C F1 = -TNLM 883C F2 = <0| TELM |0> 884C F3 = 2*<0| TELM(k1) |0> 885C F4 = <0| TELM(k1,k2) |0> 886C 887 F1=D0 888 F2=D0 889 F3=D0 890 F4=D0 891C 892 IF (ISYMT.EQ.1) THEN 893 F1 = -WRK((KTNLM-1)+LM) 894 CALL DCOPY(N2ORBX,WRK(KTELM),1,WRK(KTLMA),1) 895 CALL MELONE(WRK(KTLMA),1,UDV,D1,F2,200,'C3SOL') 896 END IF 897 IF (ISYMT.EQ.ISYMV1) THEN 898 CALL DZERO(WRK(KTLMA),N2ORBX) 899 CALL OITH1(ISYMV1,ZYM1,WRK(KTELM),WRK(KTLMA),ISYMT) 900 CALL MELONE(WRK(KTLMA),1,UDV,D1,F3,200,'C3SOL') 901 F3 = 2*F3 902 END IF 903 IF (ISYMT.EQ.MULD2H(ISYMV1,ISYMV2)) THEN 904 CALL DZERO(WRK(KTLMA),N2ORBX) 905 CALL DZERO(WRK(KTLMB),N2ORBX) 906 CALL OITH1(ISYMV1,ZYM1,WRK(KTELM),WRK(KTLMA),ISYMT) 907 CALL OITH1(ISYMV2,ZYM2,WRK(KTLMA),WRK(KTLMB),ISYMV2) 908 CALL MELONE(WRK(KTLMB),1,UDV,D1,F4,200,'C3SOL') 909 END IF 910C 911 FLST = WRK((KFLST-1)+LM) 912 FLOP = WRK((KFLOP-1)+LM) 913 IF (ISYMT.EQ.MULD2H(ISYMV1,ISYMV2)) THEN 914 FACT = FLOP*F4 915 CALL DAXPY(N2ORBX,FACT,WRK(KTELM),1,WRK(KTRES),1) 916 END IF 917 IF (ISYMT.EQ.ISYMV1) THEN 918 CALL DZERO(WRK(KTLMA),N2ORBX) 919 CALL OITH1(ISYMV2,ZYM2,WRK(KTELM),WRK(KTLMA),ISYMT) 920 FACT = FLOP*F3 921 CALL DAXPY(N2ORBX,FACT,WRK(KTLMA),1,WRK(KTRES),1) 922 END IF 923 IF (ISYMT.EQ.1) THEN 924 IF (INERSI) THEN 925 FACT = FLST*(F1+F2) 926 ELSE 927 FACT = FLOP*(F1+F2) 928 END IF 929 CALL DZERO(WRK(KTLMA),N2ORBX) 930 CALL DZERO(WRK(KTLMB),N2ORBX) 931 CALL OITH1(ISYMV1,ZYM1,WRK(KTELM),WRK(KTLMA),ISYMT) 932 CALL OITH1(ISYMV2,ZYM2,WRK(KTLMA),WRK(KTLMB), 933 & MULD2H(ISYMT,ISYMV1)) 934 CALL DAXPY(N2ORBX,FACT,WRK(KTLMB),1,WRK(KTRES),1) 935 END IF 936 500 CONTINUE 937 520 CONTINUE 938C 939C Make the gradient 940C 941C / <0| [qj ,TRES] |0> \ 942C | 0 | 943C | <0| [qj+,TRES] |0> | 944C \ 0 / 945C 946 ISYMDN = 1 947 OVLAP = D1 948 JSPIN = 0 949 ISYMV = IREFSY 950 ISYMST = MULD2H(IGRSYM,IREFSY) 951 IF ( ISYMST .EQ. IREFSY ) THEN 952 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 953 ELSE 954 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 955 END IF 956 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 957 LREF = .TRUE. 958 NZYVEC = NCREF 959 NZCVEC = NCREF 960 CALL RSP1GR(NSIM,KZYVR,IDUMMY,JSPIN,IGRSYM,JSPIN,ISYMV,ETRS, 961 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,WRK(KTRES), 962 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 963 964C 965 CALL GPCLOSE(LUSOL,'KEEP') 966 RETURN 967 END 968