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 /* Deck abarsp */ 20 SUBROUTINE ABARSP(CICLC,HFCLC,TRPCLC,OOTV,IOPSYM,EXCLC, 21 * EXVAL,NEXVAL,NABATY,NABAOP,LABEL,LUGDVE, 22 * LUSOVE,LUREVE,THRCLC,MAXCLC,IPRCLC,MXRM, 23 * MXPHP,WRK,LWRK) 24C 25C CICLC - true for CI calculations 26C HFCLC - true for RHF - closed shell or one electron in one 27C active orbital 28C TRPCLC - true for triplet perturbation operators 29C OOTV - true for optimal orbital trial vectors 30C IOPSYM - symmetry of perturbation operators 31C EXCLC - true for excitation energy calculations, 32C false for linear response equations 33C EXVAL - calculated excitation energies (output) 34C - frequency for linear response (input) 35C NEXVAL - number of excitation energies/frequencies 36C NABATY - 1 for real operator, -1 for imaginary operator (for each 37C operator) 38C NABAOP - number of right-hand sides 39C LUGDVE - unit number for right-hand sides 40C LUSOVE - unit number for solutions 41C LUREVE - unit number for residuals 42C THRCLC - threshold for convergence 43C MAXCLC - maximum number of iterations 44C IPRCLC - print level 45C MXRM - maximum size of reduced space 46C MXPHP - maximum size for explicit subblock of configuration 47C Hessian 48C In common: 49C 50C NEWCMO - true : transform integrals to MO basis because possibly new CMO coefficients, 51C normally true first time routine is called at a new geometry. 52C 53#include "implicit.h" 54#include "iratdef.h" 55#include "dummy.h" 56#include "maxorb.h" 57C 58 DIMENSION WRK(LWRK),EXVAL(NEXVAL),NABATY(NABAOP) 59C 60#include "mxcent.h" 61 PARAMETER (D0 = 0.0D0, D1 = 1.0D0) 62 CHARACTER*8 LABEL(NABAOP) 63C 64C Used from common blocks: 65C ABAINF : DODRCT, TDA_SINGLET, TDA_TRIPLET, ? 66C exeinf.h : FTRCTL, NEWCMO, ITRLVL_LAST, LVLDRC_LAST 67C INFLIN : NCONRF 68C INFINP : DIRFCK,HSORHF,FLAG(:),... 69C INFORB : NASHT,... 70C infvar.h : NVAR, ... (reset of FLAG(27) 71C INFTRA : USEDRC 72C infrsp.h : TDA, ... 73C 74#include "abainf.h" 75#include "exeinf.h" 76#include "inflin.h" 77#include "infinp.h" 78#include "cbisol.h" 79#include "inforb.h" 80#include "infvar.h" 81#include "infpri.h" 82#include "wrkrsp.h" 83#include "infrsp.h" 84#include "infdim.h" 85#include "inftra.h" 86#include "inftap.h" 87#include "rspprp.h" 88#include "inflr.h" 89#include "infpp.h" 90#include "infhyp.h" 91#include "priunit.h" 92#include "infsop.h" 93#include "abslrs.h" 94#include "qm3.h" 95C 96 LOGICAL CICLC,HFCLC,TRPCLC,EXCLC,OOTV,FLAGSV,EX, USEDRC_SAVE 97 LOGICAL FLAG27_save, TDA_save 98C 99 CALL QENTER('ABARSP') 100 101 USEDRC_SAVE = USEDRC 102 FLAG27_SAVE = FLAG(27) ! needs to be restored after triplet response 103 ! otherwise susequent call of ABARSPN for singlet response may fail /Feb.2018 - hjaaj 104 ! (SIRIFC is not read again, so it is not restored this way) 105 TDA_save = TDA 106 107 IF (MMPCM) CALL MMPCMINIT() 108 TIMRSP = SECOND() 109 IPRRSP = IPRCLC - 2 110 IF (IPRRSP .GE. 0) THEN 111 WRITE (LUPRI,'(//A/A)') 112 * ' THIS IS OUTPUT FROM THE MCSCF AND SOPPA RESPONSE SOLVER', 113 * ' ---------------------------------------------------------' 114 IF (EXCLC .AND. TRPCLC) THEN 115 WRITE(LUPRI,'(2(/A,I6),/A,1P,D13.2)') 116 * ' Symmetry of triplet excitation operator',IOPSYM, 117 * ' Number of triplet excitations:',NEXVAL, 118 * ' Convergence threshold:',THRCLC 119 IF (TDA_TRIPLET) WRITE(LUPRI,'(A)') ' * TDA is used.' 120 ELSE IF (EXCLC) THEN 121 WRITE(LUPRI,'(2(/A,I6),/A,1P,D13.2)') 122 * ' Symmetry of singlet excitation operator',IOPSYM, 123 * ' Number of singlet excitations:',NEXVAL, 124 * ' Convergence threshold:',THRCLC 125 IF (TDA_SINGLET) WRITE(LUPRI,'(A)') ' * TDA is used.' 126 ELSE IF (TRPCLC) THEN 127 WRITE(LUPRI,'(/3(/A,I6),/A,1P,D13.2)') 128 * ' Symmetry of triplet property operator', IOPSYM, 129 * ' Number of operators for triplet linear response equations:' 130 * ,NABAOP, ' Number of response frequencies:',NEXVAL, 131 * ' Convergence threshold:',THRCLC 132 IF (TDA_TRIPLET) WRITE(LUPRI,'(A)') ' * TDA is used.' 133 ELSE 134 WRITE(LUPRI,'(/3(/A,I6),/A,1P,D13.2)') 135 * ' Symmetry of singlet property operator', IOPSYM, 136 * ' Number of operators for singlet linear response equations:' 137 * ,NABAOP, ' Number of response frequencies:',NEXVAL, 138 * ' Convergence threshold:',THRCLC 139 IF (TDA_SINGLET) WRITE(LUPRI,'(A)') ' * TDA is used.' 140 ENDIF 141 END IF 142 CALL FLSHFO(LUPRI) 143 144C 145 REWIND LUGDVE 146 REWIND LUSOVE 147 REWIND LUREVE 148C 149C CONSTRUCTION OF ORBITAL DIAGONAL WITH AVDIA 150C 151 AVDIA = .TRUE. 152C 153C ************************* 154C ***** INPUT SECTION ***** 155C ************************* 156C 157C DEFINE CALCULATION BY FLAGS AND LOGICALS 158C 159 IF (TRPCLC) THEN 160 TRPLET = .TRUE. 161 TRPFLG = .TRUE. 162 TDA = TDA_TRIPLET 163 ELSE 164 TRPLET = .FALSE. 165 TRPFLG = .FALSE. 166 TDA = TDA_SINGLET 167 ENDIF 168C 169 IF (CICLC) THEN 170 RSPCI = .TRUE. 171 ELSE 172 RSPCI = .FALSE. 173 ENDIF 174C 175 IF (CICLC.AND.HFCLC) THEN 176 WRITE (LUPRI,'(//,A,/,A,/)') 177 * ' ERROR: THE CALCULATION IS SPECIFIED AS BOTH CI AND HF' 178 * ,' CONSEQUENTLY THE CALCULATION MUST STOP' 179 CALL QUIT('ABARSP: INCORRECT SPECIFICATION OF CICLC AND HFCLC') 180 ENDIF 181C 182 KSYMOP = IOPSYM 183 OPTORB = OOTV 184 MAXRM = MXRM 185 MAXPHP = MXPHP 186C 187 IF (EXCLC) THEN 188 THCPP = THRCLC 189 MAXITP = MAXCLC 190 IPRPP = IPRCLC - 2 191 NPPCNV(KSYMOP) = NEXVAL 192 NPPSIM(KSYMOP) = NEXVAL 193 NPPSTV(KSYMOP) = NEXVAL 194 ELSE 195 THCLR = THRCLC 196 MAXITL = MAXCLC 197 IPRLR = IPRCLC - 2 198 NFREQ = NEXVAL 199 DO 120 I = 1,NEXVAL 200 FREQ(I) = EXVAL(I) 201 120 CONTINUE 202 ENDIF 203C 204C 205 IF (NEXVAL.GT.0) CALL RSPSET 206C 207 CALL FLSHFO(LUPRI) 208C 209C PERFORM INTEGRAL TRANSFORMATION, if needed 210C 211 IF (NEWCMO .OR. FTRCTL) THEN 212 IF ( (NASHT .GT. 1 .AND. .NOT.HSROHF) .OR. .NOT.DIRFCK ) THEN ! no mo integrals needed for ao-direct HF or DFT 213 KCMO = 1 214 KWTRA = KCMO + NCMOT 215 LWTRA = LWRK - KWTRA 216 IF (LWTRA.LT.0) CALL ERRWRK('ABARSP 1',-(KWTRA-1),LWRK) 217 REWIND LUSIFC 218 IF (RSPCI) THEN 219 CALL MOLLAB('CIRESPON',LUSIFC,LUERR) 220 ELSE 221 CALL MOLLAB(LBSIFC,LUSIFC,LUERR) 222 END IF 223 READ (LUSIFC) 224 READ (LUSIFC) 225 CALL READT (LUSIFC,NCMOT,WRK(KCMO)) 226 227 USEDRC = .TRUE. 228 IF (RSPCI) THEN 229 USEDRC = .FALSE. 230 JTRLVL = 0 231 ELSE IF (SOPPA .OR. NEXVAL.EQ.0) THEN ! SOPPA or quadratic response 232 JTRLVL = -10 233 ELSE 234 JTRLVL = -4 235 ENDIF 236 237 IF (NEWCMO .OR. FTRCTL .OR. 238 & (ABS(JTRLVL) .GT. ITRLVL_LAST) .OR. 239 & (USEDRC .AND. LVLDRC_LAST .LT. 0) ) THEN 240 FLAGSV = DORSP 241 DORSP = .TRUE. 242 CALL SIR_INTOPEN 243 DORSP = FLAGSV 244 CALL TRACTL(JTRLVL,WRK(KCMO),WRK(KWTRA),LWTRA) 245 END IF 246 NEWCMO = .FALSE. 247 END IF 248 END IF 249C 250C ORGANIZE CALCULATION FOR EACH PERTURBATION OPERATOR 251C 252 IF (NEXVAL.EQ.0) THEN 253 KWSYM = 1 254 LWSYM = LWRK - KWSYM 255 IF (LWSYM.LT.0) CALL ERRWRK('ABARSP ',KWSYM,LWRK) 256 CALL RSPSYM(WRK(KWSYM),LWSYM) 257 END IF 258C 259C ALLOCATE WORK SPACE FOR MATRICES THAT WILL BE KEPT DURING THE 260C WHOLE RESPONSE CALCULATION AND READ IN THE MATRICES 261C 262 KFREE = 1 263 LFREE = LWRK 264 CALL MEMGET2('REAL','INDX',KINDX,LCINDX,WRK,KFREE,LFREE) 265 CALL MEMGET2('REAL','CMO', KCMO ,NCMOT ,WRK,KFREE,LFREE) 266 LUDV = NACTT*NACTT 267 CALL MEMGET('REAL',KUDV ,LUDV ,WRK,KFREE,LFREE) 268 IF (RSPCI) THEN 269 CALL MEMGET2('REAL','PVX', KPVX ,0,WRK,KFREE,LFREE) 270 CALL MEMGET2('REAL','FOCK',KFOCK,0,WRK,KFREE,LFREE) 271 CALL MEMGET2('REAL','FC', KFC ,0,WRK,KFREE,LFREE) 272 CALL MEMGET2('REAL','FV', KFV ,0,WRK,KFREE,LFREE) 273 ELSE 274 IF (TRPFLG) THEN 275C NEED BOTH TRIPLET AND SINGLET TWO ELECTRON DENSITY MATRICES 276 LPVX = 2*LPVMAT 277 ELSE 278C NEED ONLY SINGLET TWO ELECTRON DENSITY MATRIX 279 LPVX = LPVMAT 280 END IF 281 CALL MEMGET2('REAL','PVX', KPVX ,LPVX ,WRK,KFREE,LFREE) 282 CALL MEMGET2('REAL','FOCK',KFOCK,N2ORBT,WRK,KFREE,LFREE) 283 CALL MEMGET2('REAL','FC', KFC ,NNORBT,WRK,KFREE,LFREE) 284 CALL MEMGET2('REAL','FV', KFV ,NNORBT,WRK,KFREE,LFREE) 285 END IF 286 CALL MEMGET2('REAL','FCAC',KFCAC,NNASHX,WRK,KFREE,LFREE) 287 CALL MEMGET2('REAL','H2AC',KH2AC,NNASHX*NNASHX,WRK,KFREE,LFREE) 288 KTOT = KFREE 289 KWRK1 = KFREE 290 LWRK1 = LFREE 291C 292 IPRRSP = IPRCLC - 5 293C 294C For SOPPA : 295C 296 IF (SOPPA) THEN 297C 298C Initialize XINDX 299C 300 A2EXIST=.FALSE. 301 CALL DZERO(WRK(KINDX),LCINDX) 302C 303C Find address array's for SOPPA calculation 304C 305 CALL SET2SOPPA(WRK(KINDX+KABSAD-1),WRK(KINDX+KABTAD-1), 306 * WRK(KINDX+KIJSAD-1),WRK(KINDX+KIJTAD-1), 307 * WRK(KINDX+KIJ1AD-1),WRK(KINDX+KIJ2AD-1), 308 * WRK(KINDX+KIJ3AD-1),WRK(KINDX+KIADR1-1)) 309C 310 ENDIF 311C 312C ... suppress printing of SCF/MCSCF energy for IPRRSP .gt. 0 313C for normal print levels (.le. 5) /960905-hjaaj 314 CALL RSPMC(WRK(KCMO),WRK(KUDV),WRK(KPVX),WRK(KFOCK),WRK(KFC), 315 * WRK(KFV),WRK(KFCAC),WRK(KH2AC),WRK(KINDX),WRK(KWRK1), 316 * LWRK1) 317C CALL RSPMC(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC,IADR1,WRK,LWRK) 318 IPRRSP = IPRCLC - 2 319C 320 IF (TDHF .NEQV. HFCLC) THEN 321 WRITE(LUPRI,'(//A,L10/A,L10)') 322 & ' ABARSP WARNING: "HFCLC" in parameter list is',HFCLC, 323 & ' but "TDHF" after RSPMC is',TDHF 324 NWARN = NWARN + 1 325 END IF 326C 327 IF ( .NOT.TDHF ) THEN 328 CALL GETCIX(WRK(KINDX),IREFSY,IREFSY,WRK(KTOT),LFREE,0) 329 END IF 330C 331 IF ( (.NOT.TDHF) .AND. (.NOT.RSPCI) ) THEN 332C 333C ... CALCULATE ONE- AND TWO- BODY DENSITY MATRICES 334C ( THE SYMMETRIC MC TWO BODY DENSITY MATRIX CANNOT BE USED IN 335C RESPONSE CALCULATION ) 336C 337 KCREF = KWRK1 338 KTOT = KCREF + NCREF 339 LFREE = LWRK - KTOT 340 IF (LFREE.LT.0) CALL ERRWRK('ABARSP 2',-(KTOT-1),LWRK) 341 KFREE = 1 342C 343 CALL GETREF(WRK(KCREF),NCREF) 344 IF ( IPRRSP.GT.110 ) THEN 345 WRITE(LUPRI,'(/A)')' ***** ONE BODY DENSITY MATRIX '// 346 & 'from SIRIFC file' 347 CALL OUTPUT(WRK(KUDV),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 348 END IF 349 ISPIN1 = 0 350 ISPIN2 = 0 351 CALL RSPDM(IREFSY,IREFSY,NCREF,NCREF,WRK(KCREF),WRK(KCREF), 352 * WRK(KUDV),WRK(KPVX), 353 * ISPIN1,ISPIN2,.FALSE.,.FALSE., 354 * WRK(KINDX),WRK(KTOT),KFREE,LFREE) 355C CALL RSPDM(ILSYM,IRSYM,NCLDIM,NCRDIM,CL,CR, RHO1,RHO2, 356C * ISPIN1,ISPIN2,TDM,NORHO2,XNDXCI,WORK, 357C * KFREE,LFREE) 358C 359 IF (TRPLET) THEN 360 ISPIN1 = 1 361 ISPIN2 = 1 362 CALL RSPDM(IREFSY,IREFSY,NCREF,NCREF,WRK(KCREF),WRK(KCREF), 363 * WRK(KUDV),WRK(KPVX+LPVMAT), 364 * ISPIN1,ISPIN2,.FALSE.,.FALSE., 365 * WRK(KINDX),WRK(KTOT),KFREE,LFREE) 366 END IF 367C 368 IF ( IPRRSP.GT.110 ) THEN 369 WRITE(LUPRI,'(/A)') ' ** ONE BODY DENSITY MATRIX FROM RSPDM' 370 CALL OUTPUT(WRK(KUDV),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 371 END IF 372C 373C For a single open shell, one can use the TDHF code, because for 374C the open shell the one electron density matrix is 1 and the two 375C electron density matrix is 0. 376C 377 ELSE IF (TDHF .AND. NASHT .EQ. 1) THEN 378 WRK(KUDV) = D1 379 WRK(KPVX) = D0 380 IF (TRPLET) WRK(KPVX+1) = D0 381 END IF 382C 383 CALL FLSHFO(LUPRI) 384C 385C Open RSPVEC for response vectors, existing vectors may exist for restart 386C 387 CALL GPINQ('RSPVEC','EXIST',EX) 388 IF (.NOT.EX) THEN 389 CALL GPOPEN(LURSP,'RSPVEC','NEW',' ','UNFORMATTED',IDUMMY, 390 & .FALSE.) 391 WRITE (LURSP) NISH,NASH,NORB,NBAS,NSYM 392 WRITE (LURSP) 'EOFLABEL' 393 ELSE 394 CALL GPOPEN(LURSP,'RSPVEC','OLD',' ','UNFORMATTED',IDUMMY, 395 & .FALSE.) 396 END IF 397C 398 IF (ABSLRS) THEN 399 LUABSVECS = -1 400 CALL GPINQ('ABSVECS','EXIST',EX) 401 IF (.NOT.EX) THEN 402 CALL GPOPEN(LUABSVECS,'ABSVECS','NEW',' ',' ',IDUMMY,.FALSE.) 403 WRITE(LUABSVECS) 'EOFLABEL' 404 ELSE 405 CALL GPOPEN(LUABSVECS,'ABSVECS','OLD',' ',' ',IDUMMY,.FALSE.) 406 END IF 407 ENDIF 408C 409C Open files for trial and sigma vectors 410C 411 CALL GPOPEN(LURSP3,' ','UNKNOWN',' ','UNFORMATTED',IDUMMY,.FALSE.) 412 CALL GPOPEN(LURSP4,' ','UNKNOWN',' ','UNFORMATTED',IDUMMY,.FALSE.) 413 CALL GPOPEN(LURSP5,' ','UNKNOWN',' ','UNFORMATTED',IDUMMY,.FALSE.) 414C 415cs 416C IF (NEXVAL.EQ.0) THEN 417Ckr+hjaaj-nov 96: are we sure that we may not have a linear 418C response where accidentally NEXVAL .eq. 0 419C e.g. because of symmetry ? 420cs 421cs SUBROUTINE QRCALC(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 422cs * XINDX,WRK,LWRK) 423cs praticamente chiama la QRVEC di rspvec.F 424cs supponiamo che parta a calcolare 425cs 426cs CALL QRCALC(WRK(KCMO),WRK(KUDV),WRK(KPVX),WRK(KFOCK), 427cs * WRK(KFC),WRK(KFV),WRK(KFCAC),WRK(KH2AC), 428cs * WRK(KINDX),WRK(KWRK1),LWRK1) 429C ELSE 430 IF (IPRRSP.GE.0) THEN 431 WRITE(LUPRI,'(//A,I5)') 432 * ' ---------- Symmetry of excitation/property operator', 433 * KSYMOP 434 IF (EXCLC) THEN 435 WRITE(LUPRI,'(/A,I5)') ' Number of excitations: ',NEXVAL 436 ELSE 437 WRITE(LUPRI,'(2(/A,I5))') 438 * ' Number of operators for linear response equations: ' 439 * ,NABAOP, ' Number of response frequencies: ',NEXVAL 440 ENDIF 441 END IF 442C 443C DEFINE VARIABLES THAT DEPEND ON SYMMETRY 444C 445 CALL RSPVAR(WRK(KUDV),WRK(KFOCK),WRK(KFC),WRK(KFV), 446 * WRK(KFCAC),WRK(KH2AC),WRK(KINDX),WRK(KWRK1),LWRK1) 447C CALL RSPVAR(UDV,FOCK,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) 448 IF ( KZVAR.EQ.0) THEN 449 NWARN = NWARN + 1 450 WRITE(LUPRI,'(/2A)')' ****ABARSP WARNING******', 451 * ' NUMBER OF VARIABLES IN THIS SYMMETRY IS ZERO' 452 GO TO 100 453 END IF 454C 455 IF (EXCLC) THEN 456C 457C ... FIND EXCITATION ENERGIES AND TRANSITION MOMENTS 458C 459 IF (IPRRSP .GT. 0) TIMPP = SECOND() 460 CALL RSPLEX(LUSOVE,EXVAL, 461 * WRK(KCMO),WRK(KUDV),WRK(KPVX), 462 * WRK(KFC),WRK(KFV),WRK(KFCAC),WRK(KH2AC), 463 * WRK(KINDX),WRK(KWRK1),LWRK1) 464 IF (IPRRSP .GT. 0) THEN 465 TIMPP = SECOND() - TIMPP 466 WRITE (LUPRI,'(//A,F10.2,A,I2)') 467 * 'Time used in polarization propagator calculation is', 468 * TIMPP,' CPU seconds for symmetry',KSYMOP 469 END IF 470C 471 ELSE 472 IF (IPRRSP .GT. 0) TIMLR = SECOND() 473 CALL RSPLLE(LUGDVE,LUSOVE,LUREVE,NABAOP,NABATY,LABEL, 474 * WRK(KCMO),WRK(KUDV),WRK(KPVX),WRK(KFC),WRK(KFV), 475 * WRK(KFCAC),WRK(KH2AC),WRK(KINDX),WRK(KFOCK), 476 * WRK(KWRK1),LWRK1) 477 IF (IPRRSP .GT. 0) THEN 478 TIMLR = SECOND() - TIMLR 479 WRITE (LUPRI,'(//A,F10.2,A,I2)') 480 * ' Time used in linear response calculation is', 481 * TIMLR,' CPU seconds for symmetry',KSYMOP 482 END IF 483 END IF 484C 485C ENDIF 486C ... end if for QRCALC test above 487C 488 CALL FLSHFO(LUPRI) 489C 490 100 CONTINUE 491C 492 CALL GPCLOSE(LURSP3,'KEEP') 493 CALL GPCLOSE(LURSP4,'DELETE') 494 CALL GPCLOSE(LURSP5,'KEEP') 495 CALL GPCLOSE(LURSP ,'KEEP') 496 IF (LUSOL.GT.0) CALL GPCLOSE(LUSOL ,'KEEP') 497 IF (ABSLRS) CALL GPCLOSE(LUABSVECS,'KEEP') 498C 499 TIMRSP = SECOND() - TIMRSP 500 IF (IPRRSP .GT. 0) WRITE (LUPRI,'(/A,F13.2,A)') 501 * ' Total time used in response solver is',TIMRSP, 502 * ' CPU seconds.' 503C 504C ******************************************* 505C 506 USEDRC = USEDRC_SAVE 507 TDA = TDA_save 508 IF (FLAG(27) .NEQV. FLAG27_SAVE) THEN 509C ... call SETCI with new FLAG(27) and 510C reset Sirius CI common blocks for CSFs/determinants 511 FLAG(27) = FLAG27_SAVE 512 CALL SETCI(NCONF,NCDETS,LSYM,WRK,LWRK,0) 513C 514 NVAR = NCONF + NWOPT 515 NVARH = NCONF + NWOPH 516 NVARMA = NCONMA + NWOPMA 517 NCONDI = MAX(1,NCONF) 518 END IF 519 CALL QEXIT('ABARSP') 520 RETURN 521 END 522C /* Deck rsplex */ 523 SUBROUTINE RSPLEX(LUSOVE,EXVAL,CMO,UDV,PV,FC,FV,FCAC, 524 * H2AC,XINDX,WRK,LWRK) 525C 526C Purpose: 527C CONTROL CALCULATION OF EXCITATION ENERGIES AND WRITE SOLUTION 528C VECTORS ON DISK (LUSOVE) 529C 530#include "implicit.h" 531#include "dummy.h" 532 CHARACTER*8 BLANK 533 DIMENSION EXVAL(*) 534 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*) 535 DIMENSION XINDX(*),WRK(LWRK) 536C 537 PARAMETER ( D100 = 100.0D0, D0 = 0.0D0, BLANK = ' ' ) 538#include "codata.h" 539#include "priunit.h" 540C 541C 542C Used from common blocks: 543C /INFRSP/ : most items (/INFRSP/ gives control information for 544C the response calculation(s) ) 545C /WRKRSP/ : 546C 547#include "infrsp.h" 548#include "inftap.h" 549#include "wrkrsp.h" 550#include "rspprp.h" 551#include "infpp.h" 552#include "inforb.h" 553#include "infpri.h" 554#include "qm3.h" 555C 556 CALL QENTER('RSPLEX') 557C space allocation for reduced E(2) and reduced S(2) 558 KREDE = 1 559 KREDS = KREDE + MAXRM*MAXRM 560 KIBTYP = KREDS + MAXRM*MAXRM 561 KEIVAL = KIBTYP + MAXRM 562 KRESID = KEIVAL + MAXRM 563 KEIVEC = KRESID + MAXRM 564 KWRK1 = KEIVEC + MAXRM*MAXRM 565 LWRK1 = LWRK + 1 - KWRK1 566 IF (IPRPP .GT. 2) THEN 567 WRITE(LUPRI,*)' IN RSPLEX: MAXRM ',MAXRM 568 WRITE(LUPRI,*)' IN RSPLEX: LWRK ,LWRK1',LWRK,LWRK1 569 WRITE(LUPRI,*)' IN RSPLEX: THCPP ',THCPP 570 END IF 571C 572C 971201-hjaaj: changed test from 3*KZYVAR to 2*KZYVAR; 573C i.e. we let RSPCTL do the testing but we are sure of 574C enough memory for RSPEVE below. 575C 576 IF (LWRK1 .LT. 2*KZYVAR) THEN 577 WRITE (LUPRI,9000) LWRK1,2*KZYVAR 578 CALL QTRACE(LUPRI) 579 CALL QUIT('ERROR, RSPLEX: INSUFFICIENT WORK SPACE') 580 ENDIF 581 9000 FORMAT(/' RSPLEX, work space too small for 2 (z,y)-vectors', 582 * /' had:',I10,', need more than:',I10) 583C 584 KZRED = 0 585 KZYRED = 0 586 THCRSP = THCPP 587 IPRRSP = IPRPP 588 MAXIT = MAXITP 589C 590C Call RSPCTL to solve propagator eigen problem 591C 592 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 593 * .FALSE.,BLANK,BLANK,DUMMY,DUMMY,WRK(KREDE),WRK(KREDS), 594 * WRK(KIBTYP),WRK(KEIVAL),WRK(KRESID),WRK(KEIVEC), 595 * XINDX,WRK(KWRK1),LWRK1) 596C CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 597C * LINEQ,GD,REDGD,REDE,REDS, 598C * IBTYP,EIVAL,EIVEC,XINDX,WRK,LWRK) 599C 600C CALCULATE EIGENVECTORS AND store on LUSOVE for later 601C eigenvalues are stored in EXVAL 602C 603C 1) MAXIMUM NUMBER OF TRIAL VECTORS 604C 2) ALLOCATE WORK SPACE 605C 606 NSIM = MIN(KEXCNV, (LWRK1-KZYVAR)/KZYVAR) 607 KBVECS = KWRK1 608 KWRK2 = KBVECS + NSIM*KZYVAR 609 LWRK2 = LWRK - KWRK2 610C 611 DO 500 ISIM = 1,KEXCNV,NSIM 612 NBX = MIN( NSIM,(KEXCNV+1-ISIM) ) 613 CALL RSPEVE(WRK(KIBTYP),WRK(KEIVAL),WRK(KEIVEC),WRK(KBVECS), 614 * WRK(KWRK2),NBX,(ISIM-1)) 615C CALL RSPEVE(IBTYP,EIVAL,EIVEC,BVECS,WRK,NBX,IBOFF) 616 DO 550 INUM = 1,NBX 617 EXVAL(ISIM-1+INUM) = WRK(KEIVAL-1+ISIM-1+INUM) 618 CALL WRITT(LUSOVE,KZYVAR,WRK(KBVECS+(INUM-1)*KZYVAR)) 619 IF (IPRPP .GT. 0) THEN 620 WRITE(LUPRI,'(/A,I5,//A,1P,G16.8,A,3(/20X,G16.8,A),/)') 621 * ' STATE NO:',(ISIM-1+INUM), 622 * ' EXCITATION ENERGY :',WRK(KEIVAL-1+ISIM-1+INUM), 623 * ' au', 624 * WRK(KEIVAL-1+ISIM-1+INUM)*XTEV, ' eV', 625 * WRK(KEIVAL-1+ISIM-1+INUM)*XTKAYS,' cm-1', 626 * WRK(KEIVAL-1+ISIM-1+INUM)*XKJMOL,' kj / mole' 627 IF (SOPPA) THEN 628 TXNRM = DNRM2(KZCONF,WRK(KBVECS+(INUM-1)*KZYVAR),1) 629 TYNRM = DNRM2(KZCONF,WRK(KBVECS+(INUM-1) 630 & *KZYVAR+KZVAR),1) 631 T2P2HN = D100*(TXNRM*TXNRM + TYNRM*TYNRM) 632 TPHNRM = D100 - T2P2HN 633 WRITE(LUPRI,'(2(/A,F8.2,A))') 634 & ' SOPPA p-h weight in excitation operator:', 635 & TPHNRM,' %', 636 & ' SOPPA 2p-2h weight in excitation operator:', 637 & T2P2HN,' %' 638 END IF 639 END IF 640 IF (IPRPP .GT. 2) THEN 641 WRITE (LUPRI,'(/A,I3)') 642 & ' EIGENVECTOR FOR STATE NO.',(ISIM-1+INUM) 643 CALL RSPPRO(WRK(KBVECS+(INUM-1)*KZYVAR+KZCONF),KZVAR, 644 * UDV,LUPRI) 645C CALL RSPPRC(WRK(KBVECS+(INUM-1)*KZYVAR),KZCONF,KZVAR,LUPRI) 646 CALL RSPANC(WRK(KBVECS+(INUM-1)*KZYVAR),KZCONF,KZVAR, 647 * MULD2H(KSYMOP,IREFSY),XINDX,MULD2H,LUPRI) 648 END IF 649 CALL WRTRSP(LURSP,KZYVAR,WRK(KBVECS + (INUM-1)*KZYVAR), 650 & 'EXCITLAB',BLANK,WRK(KEIVAL-1+ISIM-1+INUM),D0, 651 & KSYMOP,0,WRK(KRESID-1+ISIM-1+INUM),D0) 652 550 CONTINUE 653 500 CONTINUE 654C 655C *** END OF RSPLEX 656C 657 CALL QEXIT('RSPLEX') 658 RETURN 659 END 660C /* Deck rsplle */ 661 SUBROUTINE RSPLLE(LUGDVE,LUSOVE,LUREVE,NABAOP,NABATY,LABEL,CMO, 662 * UDV,PV,FC,FV,FCAC,H2AC,XINDX,FOCK,WRK,LWRK) 663C 664#include "implicit.h" 665#include "iratdef.h" 666#include "dummy.h" 667#include "mxcent.h" 668C 669 CHARACTER*8 BLANK, LABEL(NABAOP) 670 PARAMETER (BLANK = ' ') 671 LOGICAL RESFLG(8) 672 DIMENSION NABATY(*) 673 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*) 674 DIMENSION XINDX(*),FOCK(*),WRK(LWRK) 675C 676#include "priunit.h" 677#include "gdvec.h" 678#include "infpri.h" 679#include "infrsp.h" 680#include "wrkrsp.h" 681#include "rspprp.h" 682#include "inflr.h" 683#include "inftap.h" 684#include "infdim.h" 685#include "inforb.h" 686#include "abainf.h" 687#include "absorp.h" 688#include "abslrs.h" 689#include "infvar.h" 690#include "dftcom.h" 691C 692 PARAMETER ( DM1 = -1.0D0 ) 693C 694C Save flag to control that RESPONANT is only called once 695C for each symmetry 696C 697 SAVE RESFLG 698 DATA RESFLG/.FALSE.,.FALSE.,.FALSE.,.FALSE., 699 & .FALSE.,.FALSE.,.FALSE.,.FALSE./ 700C 701C DETERMINE SECOND ORDER MOLECULAR PROPERTIES 702C 703 CALL QENTER('RSPLLE') 704C 705 KFREE = 1 706 LFREE = LWRK 707 IF (ABSLRS .AND. NCONF.LE.1) THEN 708 CALL MEMGET('REAL',KRES,2*3*3*NFREQ_INTERVAL,WRK,KFREE, 709 & LFREE) 710 CALL MEMGET('REAL',KGD,KZYVAR,WRK,KFREE,LFREE) 711 CALL MEMGET('REAL',KXSOL,2*NFREQ_ALPHA*KZYVAR, 712 & WRK,KFREE,LFREE) 713 CALL MEMGET('REAL',KMJWOP,(2*8*MAXWOP + 1)/IRAT,WRK, 714 & KFREE,LFREE) 715 CALL ABS2INTER(WRK(KMJWOP),KZVAR) 716 CALL GETGPV(LABEL,FC,FV,CMO,UDV,PV,XINDX,ANTSYM, 717 & WRK(KFREE),LFREE) 718 CALL DCOPY(KZYVAR,WRK(KFREE),1,WRK(KGD),1) 719 GDNRM=DNRM2(2*KZVAR,WRK(KGD),1) 720 IF (GDNRM .LE. 1.0d-8) THEN 721 WRK(KXSOL:4*KZVAR-1)=0.0d0 722 RES0=0.0d0 723 DO K=1,NFREQ_ALPHA 724C CALL WRITE_XVEC(LUABSVECS,4*KZVAR,WRK(KXSOL),LABEL, 725C & ABS_FREQ_ALPHA(K),RES0) 726 CALL WRITE_XVEC2(LUABSVECS,4*KZVAR,WRK(KXSOL),LABEL, 727 & ' ',ABS_FREQ_ALPHA(K),0.0D0,RES0) 728 ENDDO 729 WRITE(LUPRI,*) 'Vectors equal to zero' 730 ELSE 731 NGD=1 732 CALL ABS_CTL(LABEL,KZVAR,WRK(KGD),NGD,WRK(KXSOL), 733 & ABS_FREQ_ALPHA,NFREQ_ALPHA,LUABSVECS, 734 & WRK(KMJWOP),CMO,UDV,FC,FCAC,FV,PV,XINDX, 735 & WRK(KRES),WRK(KFREE),LFREE) 736 ENDIF 737 ELSE 738 CALL MEMGET('REAL',KREDE,MAXRM*MAXRM,WRK,KFREE,LFREE) 739 CALL MEMGET('REAL',KREDS,MAXRM*MAXRM,WRK,KFREE,LFREE) 740 CALL MEMGET('REAL',KREDZ,2*MAXRM*MAXRM,WRK,KFREE,LFREE) 741 CALL MEMGET('INTE',KIBTYP,MAXRM,WRK,KFREE,LFREE) 742 CALL MEMGET('REAL',KEIVAL,MAXRM,WRK,KFREE,LFREE) 743 CALL MEMGET('REAL',KRESID,MAXRM,WRK,KFREE,LFREE) 744 CALL MEMGET('REAL',KEIVEC,2*MXFREQ*MAXRM,WRK,KFREE,LFREE) 745 CALL MEMGET('REAL',KREDGD,MAXRM,WRK,KFREE,LFREE) 746 CALL MEMGET('REAL',KREDZGD,2*MAXRM,WRK,KFREE,LFREE) 747 CALL MEMGET('REAL',KGD,KZYVAR,WRK,KFREE,LFREE) 748C 749 IF (ABSORP .OR. ABSLRS) THEN 750C 751C Nonzero damping parameter corresponds to complex response including 752C absorption in the system 753C 754 IF (.NOT.RESFLG(KSYMOP)) THEN 755 CALL RESONANT(WRK(KREDE),WRK(KREDS),WRK(KIBTYP), 756 & WRK(KEIVAL),WRK(KRESID),WRK(KEIVEC),CMO,UDV,PV, 757 & FOCK,FC,FV,FCAC,H2AC,XINDX,DUMMY,WRK(KFREE),LFREE) 758 RESFLG(KSYMOP) = .TRUE. 759 END IF 760C 761 CALL MEMGET('REAL',KRES,2*3*3*NFREQ_INTERVAL,WRK,KFREE,LFREE) 762 CALL ABSCTL(KOPER,LABEL,WRK(KREDE),WRK(KREDS), 763 & WRK(KREDZ),WRK(KREDGD),WRK(KREDZGD), 764 & WRK(KEIVEC),WRK(KIBTYP), 765 & CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC,XINDX, 766 & WRK(KRES),WRK(KFREE),LFREE) 767C 768 ELSE 769C 770C Infinite lifetime of the excited states, 771C corresponds to the "regular" real response 772C 773 IF (LFREE.LT.3*KZYVAR) THEN 774 WRITE (LUPRI,9100) LFREE,3*KZYVAR 775 CALL QTRACE(LUPRI) 776 CALL QUIT('RSPLLE, INSUFFICIENT SPACE TO SOLVE '// 777 & 'LINEAR EQUATIONS') 778 ENDIF 779 9100 FORMAT(/' RSPLLE, work space too small for 3 (z,y)-vectors', 780 & /' had:',I10,', need more than:',I10) 781C WORK SPACE FOR RSPEVE AND RSPRES 782C WORK SPACE EACH TIME RSPRES IS CALLED 783 MWRK = KZYVAR 784 IF (.NOT. RSPCI ) MWRK = MWRK + NCREF + LACIMX 785C Work space for each residual vector 786 MRES = MAX(NORBT*NORBT,2*N2ASHX) + KZYVAR 787C NUMBER OF SIMULTANEOUS VECTORS 788 MXSIM = (LFREE-MWRK)/MRES 789 MXSIM = MIN(NFREQ,MXSIM) 790 IF (MXSIM .LE. 0) THEN 791 WRITE(LUPRI,*) 'RSPLLE: insufficient memory for RSPRES' 792 WRITE(LUPRI,*) ' had:',LFREE,' need more than:',MWRK+MRES 793 CALL QUIT('RSPLLE: insufficient memory for RSPRES') 794 END IF 795C 796 KWRKE = KFREE 797 KBVECS = KWRKE + KZYVAR 798C 799 KZRED = 0 800 KZYRED = 0 801 THCRSP = THCLR 802 IPRRSP = IPRLR 803 MAXIT = MAXITL 804C 805 DO IOP = 1,NABAOP 806 IF (IPRLR .GE. 0) 807 & WRITE (LUPRI,'(//A,I3,/A,I3,A,I3,2A/A,(T25,5F10.6))') 808 & ' RSPLLE -- linear response calculation for symmetry', 809 & KSYMOP, 810 & ' RSPLLE -- operator no.',IOP,' out of',NABAOP,': ', 811 & LABEL(IOP), 812 & ' RSPLLE -- frequencies :',(FREQ(I),I=1,NFREQ) 813 IF (ABASOP .AND. IPRLR .GE. 0) 814 & WRITE (LUPRI,'(/A/A,I8,/A,I8,/A,I8)') 815 & ' SOPPA calculation:' 816 & ,' p-h variables. - KZWOPT :',KZWOPT 817 & ,' 2p-2h variables. - KZCONF :',KZCONF 818 & ,' Total number of variables - KZVAR :',KZVAR 819C 820C READ IN GD VECTOR FROM LUGDVE 821C 822C added by Bin Gao, June 4, 2009 823C for general GD vectors 824 if ( NABATY(IOP) .le. -2 ) then 825 call READT( LUGDVE, KZYVAR, WRK(KGD) ) 826 DFTIMG = .true. 827 else if ( NABATY(IOP) .ge. 2 ) then 828 call READT( LUGDVE, KZYVAR, WRK(KGD) ) 829 DFTIMG = .false. 830 else 831 CALL READT(LUGDVE,KZVAR,WRK(KGD)) 832C 833 IF (NABATY(IOP).EQ.-1) THEN 834 DFTIMG = .TRUE. 835 ELSE 836 DFTIMG = .FALSE. 837 END IF 838 IF (NABATY(IOP).EQ.1) THEN 839 CALL DCOPY(KZVAR,WRK(KGD),1,WRK(KGD+KZVAR),1) 840 CALL DSCAL(KZVAR,DM1,WRK(KGD+KZVAR),1) 841 ELSE IF (NABATY(IOP).EQ.-1) THEN 842 CALL DCOPY(KZVAR,WRK(KGD),1,WRK(KGD+KZVAR),1) 843 ELSE 844 WRITE(LUPRI,'(/2A,I5,A,I5)') 845 & ' NABATY(IOP) NOT AN ALLOWED VALUE' 846 & ,' IOP=',IOP,' NABATY(IOP)=',NABATY(IOP) 847 CALL QUIT(' RSPLLE: INCORRECT NABATY(IOP)') 848 END IF 849 end if 850 851 DNORM_GD = DNRM2(KZYVAR,WRK(KGD),1) 852 IF (SOPPA) THEN ! for small test cases you may have that 2p-2h non-zero, but p-h zero 853 KGDO_X = KGD + KZCONF 854 KGDO_Y = KGDO_X + KZVAR 855 DNORM_GDORB = DDOT(KZWOPT,WRK(KGDO_X),1,WRK(KGDO_X),1) 856 & + DDOT(KZWOPT,WRK(KGDO_Y),1,WRK(KGDO_Y),1) 857 DNORM_GDORB = SQRT(DNORM_GDORB) 858 IF (DNORM_GDORB .LT. 1.0D-9) THEN 859 WRITE (LUPRI,'(//A/3A,1P,D10.2/A,D10.2)') 860 & ' Solving SOPPA linear response skipped because norm of', 861 & ' p-h property vector '//LABEL//' is only',DNORM_GDORB, 862 & ' although 2p-2h property vecor has norm ',DNORM_GD 863 DNORM_GD = 0.0D0 864 END IF 865 END IF 866C 867 DO I = 1,NFREQ 868 WRK(KEIVAL-1+I) = FREQ(I) 869 END DO 870 KEXSIM = NFREQ 871 KEXCNV = NFREQ 872C 873C Call RSPCTL to solve linear set of response equations 874C 875 IF (DNORM_GD .LT. 1.0D-9) THEN 876 877 WRITE (LUPRI,'(//A/3A,1P,D10.2)') 878 & ' Solving linear response skipped because norm of', 879 & ' property vector '//LABEL//' is only',DNORM_GD 880 881 ELSE ! else DNORM_GD is .ge. 1.0D-9, and we solve. 882 883 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 884 & .TRUE.,LABEL,BLANK,WRK(KGD),WRK(KREDGD), 885 & WRK(KREDE),WRK(KREDS),WRK(KIBTYP),WRK(KEIVAL), 886 & WRK(KRESID),WRK(KEIVEC),XINDX,WRK(KFREE),LFREE) 887 888 END IF ! IF (DNORM_GD .LT. 1.0D-9) THEN .. ELSE .. 889 890 DFTIMG = .FALSE. 891C 892 DO IFREQ = 1,NFREQ,MXSIM 893 NBX = MIN(MXSIM,(NFREQ+1-IFREQ)) 894 IBOFF = IFREQ - 1 895 IF (DNORM_GD .LT. 1.0D-9) THEN 896 CALL DZERO(WRK(KBVECS),NBX*KZYVAR) 897 ELSE 898 CALL RSPEVE(WRK(KIBTYP),WRK(KEIVAL),WRK(KEIVEC), 899 & WRK(KBVECS),WRK(KWRKE),NBX,IBOFF) 900 END IF 901 DO IVEC = 1,NBX 902 IBV = (IVEC-1)*KZYVAR + KBVECS 903 CALL WRITT(LUSOVE,KZYVAR,WRK(IBV)) 904 IF (IPRLR .GE. 2) THEN 905 WRITE (LUPRI,'(/A,I5)') ' EIGENVECTOR NUMBER ', 906 & IFREQ-1+IVEC 907 CALL RSPPRO(WRK(IBV+KZCONF),KZVAR,UDV,LUPRI) 908 CALL RSPANC(WRK(IBV),KZCONF,KZVAR, 909 & MULD2H(KSYMOP,IREFSY),XINDX,MULD2H,LUPRI) 910 END IF 911C added by Bin Gao, June 5, 2009 912C for general case 913 ANTSYM = sign( DM1, dble( NABATY(IOP) ) ) 914C NATSYM = 1.0D0 915 CALL WRTRSP(LURSP,KZYVAR,WRK(IBV),LABEL, 916 & BLANK,WRK(KEIVAL-1+IVEC),D0,KSYMOP,0, 917 & WRK(KRESID-1+IVEC),ANTSYM) 918 END DO ! IVEC 919 CALL RSPRES(WRK(KIBTYP),WRK(KEIVAL),WRK(KEIVEC), 920 * WRK(KGD),WRK(KFREE),XINDX,LFREE,UDV,NBX,IBOFF) 921 DO IVEC = 1,NBX 922 IBV = (IVEC-1)*KZYVAR + KFREE 923 CALL DSCAL(KZYVAR,DM1,WRK(IBV),1) 924 CALL WRITT(LUREVE,KZYVAR,WRK(IBV)) 925 END DO ! IVEC 926 END DO ! IFREQ 927C 928C End loop over B operators 929C 930 END DO 931 END IF 932 END IF 933C 934 CALL QEXIT('RSPLLE') 935 RETURN 936 END 937C /* Deck prpaba */ 938 SUBROUTINE PRPABA(WORD,GPLON,LWRK) 939C 940C GET RIGHT HAND SIDE LONDON VECTOR in GPLON(1) 941C for RESPONSE module 942C 943#include "implicit.h" 944#include "priunit.h" 945#include "iratdef.h" 946C 947 CHARACTER*(*)WORD 948 DIMENSION GPLON(LWRK) 949C 950C Used from common blocks: 951C INFDIM : NVARMA 952C INFLIN : NVARPT 953C 954#include "mxcent.h" 955#include "maxmom.h" 956#include "maxorb.h" 957#include "maxaqn.h" 958C 959#include "inftap.h" 960#include "nuclei.h" 961#include "symmet.h" 962C 963#include "infdim.h" 964#include "inflin.h" 965C 966 LOGICAL OLDDX 967C 968 CALL GPOPEN(LUGDI,ABAGDI,'OLD','DIRECT',' ',IRAT*NVARMA,OLDDX) 969C 970 IF (WORD(1:7).EQ.'XLONMAG') THEN 971 IREC = 3*NUCDEP + IPTAX(1,2) 972 ELSE IF (WORD(1:7).EQ.'YLONMAG') THEN 973 IREC = 3*NUCDEP + IPTAX(2,2) 974 ELSE IF (WORD(1:7).EQ.'ZLONMAG') THEN 975 IREC = 3*NUCDEP + IPTAX(3,2) 976 ELSE 977 WRITE(LUPRI,'(2A)') 'Wrong label in PRPABA : ',WORD 978 CALL QUIT('Wrong Label in PRPABA') 979 ENDIF 980C 981 CALL READDX(LUGDI,IREC,IRAT*NVARPT,GPLON) 982C imaginary operator, i.e. Y = Z: 983 CALL DCOPY(NVARPT,GPLON(1),1,GPLON(1+NVARPT),1) 984C 985 CALL GPCLOSE(LUGDI,'KEEP') 986C 987C END OF PRPABA 988C 989 RETURN 990 END 991