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 absorp.F: Main author Patrick Norman 20C This implementation is published in: 21C 22C P. Norman, D.M. Bishop, H.J.Aa. Jensen, and J. Oddershede, 23C "Near-resonant absorption in the time-dependent self-consistent 24C field and multiconfigurational self-consistent field approximations", 25C J. Chem. Phys. 115 (2001) 10323-10334. 26C 27C and 28C 29C P. Norman, D.M. Bishop, H.J.Aa. Jensen, and J. Oddershede, 30C "Nonlinear response theory with relaxation: the first 31C hyperpolarizability", 32C J. Chem. Phys. 123 (2005) 194103. 33C ============================================================================ 34C 35 SUBROUTINE ABSORP_INTERFACE 36C 37C Purpose: 38C Read in user settings for imaginary polarizabilities describing 39C absorption in the optical processes. 40C 41#include "implicit.h" 42#include "priunit.h" 43#include "absorp.h" 44#include "abslrs.h" 45#include "codata.h" 46#include "infrsp.h" 47C 48 ABS_ALPHA=ABSLRS_ALPHA 49 ABS_BETA =ABSLRS_BETA 50 ABS_GAMMA=ABSLRS_GAMMA 51 ABS_MCD=ABSLRS_MCD 52 ABS_SHG=ABSLRS_SHG 53c 54 NFREQ_ALPHA=ABS_NFREQ_ALPHA 55 DO I=1,NFREQ_ALPHA 56 FREQ_ALPHA(I)=ABS_FREQ_ALPHA(I) 57 ENDDO 58 NFREQ_BETA_B=ABS_NFREQ_BETA_B 59 NFREQ_BETA_C=ABS_NFREQ_BETA_C 60 DO I=1,NFREQ_BETA_B 61 FREQ_BETA_B(I)=ABS_FREQ_BETA_B(I) 62 ENDDO 63 DO I=1,NFREQ_BETA_C 64 FREQ_BETA_C(I)=ABS_FREQ_BETA_C(I) 65 ENDDO 66 DO I=1,3 67 FREQ_INTERVAL(I)=ABS_FREQ_INTERVAL(I) 68 ENDDO 69c 70 IPRABS=IPRABSLRS 71 DAMPING=ABS_DAMP 72 ABS_ANALYZE=ABSLRS_ANALYZE 73 ABS_INTERVAL=ABSLRS_INTERVAL 74 MAXRM=MAX(MAXRM,ABS_MAXRM) 75 IF (ABS_INTERVAL) ABS_REDUCE = .TRUE. 76 77C 78 IF (ABS_OLDCPP) THEN 79 CALL AROUND('Variable settings for absorption calculation') 80C 81 IF (ABS_ALPHA) THEN 82 WRITE (LUPRI,'(A,L4)') 83 & ' Linear polarizability calculation requested: ABS_ALPHA =', 84 & ABS_ALPHA 85 IF (.NOT.ABS_INTERVAL) WRITE(LUPRI,'(A,5(4F12.8,/,16X))') 86 & ' at frequencies:',(FREQ_ALPHA(I),I=1,NFREQ_ALPHA) 87 END IF 88C 89 IF (ABS_BETA) THEN 90 WRITE (LUPRI,'(A,L4)') 91 & ' Nonlinear polarizability calculation requested: ABS_BETA =', 92 & ABS_BETA 93 WRITE(LUPRI,'(A,5(4F12.8))') 94 & ' B-FREQ:',(FREQ_BETA_B(I),I=1,NFREQ_BETA_B) 95 WRITE(LUPRI,'(A,5(4F12.8))') 96 & ' C-FREQ:',(FREQ_BETA_C(I),I=1,NFREQ_BETA_C) 97 END IF 98C 99 IF (ABS_MCD) THEN 100 WRITE (LUPRI,'(A,L4)') 101 & ' Magnetic circular dichroism requested : ABS_MCD =', 102 & ABS_MCD 103 WRITE(LUPRI,'(A,5(4F12.8))') 104 & ' at frequencies:',(FREQ_BETA_B(I),I=1,NFREQ_BETA_B) 105 END IF 106C 107 WRITE(LUPRI,'(A,L4)') 108 & ' Absorption calc over frequency interval : ABS_INTERVAL =' 109 & , ABS_INTERVAL 110 IF (ABS_INTERVAL) THEN 111 WRITE(LUPRI,'(A,I4)') 112 & ' Number of frequencies per batch : NFREQ_BATCH =' 113 & ,NFREQ_BATCH 114 WRITE(LUPRI,'(3(A,F8.5))') 115 & ' Start:',FREQ_INTERVAL(1),' Stop:',FREQ_INTERVAL(2), 116 & ' Step:',FREQ_INTERVAL(3) 117 END IF 118C 119 IF (NEXCITED_STATES .GT. MXSTATES) THEN 120 WRITE(LUPRI,'(/2A,/)') ' --- Warning the number of excited', 121 & ' exceeds maximun, the maximum value is used ---' 122 NEXCITED_STATES = MXSTATES 123 END IF 124C 125C 126 WRITE(LUPRI,'(A,F9.6,F8.1)') 127 & ' Damping parameter (a.u. and cm-1) : DAMPING =' 128 & ,DAMPING,DAMPING*XTKAYS 129 WRITE(LUPRI,'(A,I4)') 130 & ' Number of states used in start iteration : NEXCIT =' 131 & ,NEXCITED_STATES 132 WRITE(LUPRI,'(A,1P,D8.1)') 133 & ' Threshold of convergence in LR solver : THCLR =' 134 & ,THCLR_ABSORP 135 WRITE(LUPRI,'(A,I4)') 136 & ' Maximum iterations in complex LR solver : MAX_MACRO =' 137 & ,MAX_MACRO 138 WRITE(LUPRI,'(A,I4)') 139 & ' Maximum iterations in real LR solver : MAX_MICRO =' 140 & ,MAX_MICRO 141 WRITE(LUPRI,'(A,I4)') 142 & ' Max iter. in optimal orbital algorithm : MAX_ITORB =' 143 & ,MAX_ITORB 144 WRITE(LUPRI,'(A,I4)') 145 & ' Print level in absorption modules : IPRABS =' 146 & ,IPRABS 147 END IF 148C 149C End of ABSORP_INPUT 150C 151 RETURN 152 END 153 SUBROUTINE ABSCALC(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 154 & XINDX,WRK,LWRK) 155C 156C PURPOSE: 157C Driver routine for the computation of response properties including 158C absorption 159C 160#include "implicit.h" 161#include "priunit.h" 162#include "iratdef.h" 163#include "absorp.h" 164#include "infvar.h" 165#include "infrsp.h" 166 167C 168 DIMENSION CMO(*),UDV(*),PV(*),FOCK(*),FC(*),FV(*),FCAC(*),H2AC(*) 169 DIMENSION XINDX(*),WRK(LWRK) 170C 171 IPRRSP = IPRABS - 1 172 KFREE = 1 173 LFREE = LWRK 174C 175 CALL MEMGET('REAL',KMJWOP,(2*8*MAXWOP + 1)/IRAT,WRK, 176 & KFREE,LFREE) 177C 178 IF (ABS_BETA .OR. ABS_MCD) THEN 179 CALL BETA_SETUP 180 END IF 181C 182C Allocate memory for results of linear absorption calc. 183C 184 IF (ABS_INTERVAL) THEN 185 NFREQ_INTERVAL = 1 + 186 & INT((FREQ_INTERVAL(2)-FREQ_INTERVAL(1))/FREQ_INTERVAL(3)) 187 ELSE 188 NFREQ_INTERVAL = NFREQ_ALPHA 189 END IF 190C 191 CALL MEMGET('REAL',KRES,2*NFREQ_INTERVAL*3*3*2,WRK,KFREE,LFREE) 192C 193C Determine linear response vectors 194C 195 CALL AROUND('Solving Linear Response Equations') 196 CALL ABSVEC1(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 197 & XINDX,WRK(KMJWOP),WRK(KRES),WRK(KFREE),LFREE) 198C 199C Determine nonlinear response functions 200C 201 IF (ABS_BETA.OR.ABS_MCD) THEN 202 CALL AROUND('Evaluate Quadratic Response Functions') 203 CALL ABSQRF(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 204 & XINDX,WRK(KMJWOP),WRK(KFREE),LFREE) 205 END IF 206C 207C Print-out a summary of results 208C 209 CALL ABSRESULT(WRK(KMJWOP),WRK(KRES),WRK(KFREE),LFREE) 210C 211 RETURN 212 END 213 SUBROUTINE ABSVEC1(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 214 & XINDX,MJWOP,RESLRF,WRK,LWRK) 215C 216C PURPOSE: 217C Solve the linear response equations and store response vectors on 218C disk. 219C 220#include "implicit.h" 221#include "priunit.h" 222#include "absorp.h" 223#include "infrsp.h" 224#include "inforb.h" 225#include "wrkrsp.h" 226#include "infvar.h" 227C 228 DIMENSION CMO(*),UDV(*),PV(*),FOCK(*),FC(*),FV(*),FCAC(*),H2AC(*) 229 DIMENSION XINDX(*),MJWOP(2,MAXWOP,8) 230 DIMENSION RESLRF(2,NFREQ_INTERVAL,3,3,2),WRK(LWRK) 231C 232 KFREE = 1 233 LFREE = LWRK 234 CALL MEMGET('REAL',KREDE,MAXRM*MAXRM,WRK,KFREE,LFREE) 235 CALL MEMGET('REAL',KREDS,MAXRM*MAXRM,WRK,KFREE,LFREE) 236 CALL MEMGET('REAL',KREDZ,2*MAXRM*MAXRM,WRK,KFREE,LFREE) 237 CALL MEMGET('INTE',KIBTYP,MAXRM,WRK,KFREE,LFREE) 238 CALL MEMGET('REAL',KEIVAL,MAXRM,WRK,KFREE,LFREE) 239 CALL MEMGET('REAL',KRESID,NEXCITED_STATES,WRK,KFREE,LFREE) 240 CALL MEMGET('REAL',KEIVEC,2*MXFREQ*MAXRM,WRK,KFREE,LFREE) 241 CALL MEMGET('REAL',KREDGD,MAXRM,WRK,KFREE,LFREE) 242 CALL MEMGET('REAL',KREDZGD,2*MAXRM,WRK,KFREE,LFREE) 243C 244 IF (IPRABS.GT.0) CALL TIMER('START ',TIMSTR,TIMEND) 245C 246 DO ISYM=1,NSYM 247 IF (NOPER(ISYM).GT.0) THEN 248 KSYMOP = ISYM 249 CALL RESONANT(WRK(KREDE),WRK(KREDS),WRK(KIBTYP), 250 & WRK(KEIVAL),WRK(KRESID),WRK(KEIVEC),CMO,UDV,PV, 251 & FOCK,FC,FV,FCAC,H2AC,XINDX,MJWOP,WRK(KFREE),LFREE) 252C 253 DO IOPER=1,NOPER(KSYMOP) 254 CALL ABSCTL(IOPER,LABOP(IOPER,KSYMOP), 255 & WRK(KREDE),WRK(KREDS), 256 & WRK(KREDZ),WRK(KREDGD),WRK(KREDZGD), 257 & WRK(KEIVEC),WRK(KIBTYP), 258 & CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC,XINDX, 259 & RESLRF,WRK(KFREE),LFREE) 260 END DO 261 END IF 262 END DO 263 IF (IPRABS.GT.0) CALL TIMER('ABSVEC',TIMSTR,TIMEND) 264 RETURN 265 END 266 SUBROUTINE RESONANT(REDE,REDS,IBTYP,EIGVAL,RESIDUAL,EIGVEC, 267 & CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC,XINDX,MJWOP,WRK,LWRK) 268C 269C PURPOSE: 270C Solve the eigenvalue equation ( E[2] - w*S[2] )*X = 0 for 271C the few lowest excited states to be used as startvectors. 272C 273#include "implicit.h" 274#include "priunit.h" 275#include "pgroup.h" 276#include "absorp.h" 277#include "dummy.h" 278#include "infrsp.h" 279#include "wrkrsp.h" 280#include "infvar.h" 281C 282 DIMENSION REDE(*),REDS(*),IBTYP(*),EIGVAL(*),RESIDUAL(*),EIGVEC(*) 283 DIMENSION CMO(*),UDV(*),PV(*),FOCK(*),FC(*),FV(*),FCAC(*),H2AC(*) 284 DIMENSION XINDX(*),MJWOP(2,MAXWOP,8),WRK(LWRK) 285 CHARACTER*8 BLANK 286C 287 PARAMETER ( BLANK=' ' ) 288C 289 IF (IPRABS.GT.0) THEN 290 WRITE(LUPRI,'(/A/2(A,I4),3A/A)') 291 & ' RESONANT -- Solve the eigenvalue equation'// 292 & ' ( E[2] - w*S[2] )*X = 0', 293 & ' RESONANT -- for the',NEXCITED_STATES, 294 & ' lowest excited states of symmetry',KSYMOP, 295 & ' ( ',REP(KSYMOP-1),')', 296 & ' RESONANT -- to be used as startvectors in LR solver.' 297 END IF 298C 299 KFREE = 1 300 LFREE = LWRK 301C 302C Initialize variables 303C 304 THCRSP = THCPP_ABSORP 305 MAXIT = MAX_MICRO 306 CALL RSPVAR(UDV,FOCK,FC,FV,FCAC,H2AC,XINDX,WRK(KFREE),LFREE) 307 IF (ABS_BETA .OR. ABS_MCD .OR. ABS_ANALYZE) THEN 308 CALL SETZY(MJWOP) 309 END IF 310C 311 IF (NEXCITED_STATES.EQ.0) RETURN 312C 313 KZRED = 0 314 KZYRED = 0 315 KEXCNV = NEXCITED_STATES 316 KEXSTV = KEXCNV 317 KEXSIM = KEXCNV 318 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 319 & .FALSE.,BLANK,BLANK,DUMMY,DUMMY,REDE,REDS, 320 & IBTYP,EIGVAL,RESIDUAL,EIGVEC,XINDX,WRK(KFREE),LFREE) 321C 322 CALL DCOPY(KEXCNV,EIGVAL,1,EXC_ENERGY(1,KSYMOP),1) 323C 324 RETURN 325 END 326 SUBROUTINE ABSREDUC(LABEL,NFREQ_ABS,FREQ_ABS, 327 & REDE,REDS,REDZ,REDGD,REDZGD,IBTYP, 328 & REDVEC,FC,FV,CMO,UDV,PV,XINDX,NBATCH,RESLRF,WRK,LWRK) 329C 330#include "implicit.h" 331#include "priunit.h" 332#include "dummy.h" 333#include "infrsp.h" 334#include "wrkrsp.h" 335#include "absorp.h" 336#include "inftap.h" 337#include "ibndxdef.h" 338#include "qrinf.h" 339C 340C PURPOSE: 341C Create start vectors from the reduced eigenvalue problem 342C ( E[2] - {w+iW}*S[2] )* (NR + iNI) = B[1] 343C 344 CHARACTER*8 LABEL 345 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),XINDX(*),IBTYP(*), 346 & WRK(LWRK),REDE(MAXRM*MAXRM),REDS(MAXRM*MAXRM), 347 & REDZ(2*MAXRM*MAXRM),REDGD(MAXRM),REDZGD(2*MAXRM), 348 & REDVEC(2*MAXRM,MXFREQ),RESLRF(2,NFREQ_INTERVAL,3,3,2), 349 & FREQ_ABS(MXFREQ) 350C 351 KFREE = 1 352 LFREE = LWRK 353 CALL MEMGET('INTE',KIPIV,MAXRM,WRK,KFREE,LFREE) 354 CALL MEMGET('REAL',KBVEC,KZYVAR,WRK,KFREE,LFREE) 355 CALL MEMGET('REAL',KGD,KZYVAR,WRK,KFREE,LFREE) 356 CALL MEMGET('REAL',KREDE,KZYRED*KZYRED,WRK,KFREE,LFREE) 357 CALL MEMGET('REAL',KREDS,KZYRED*KZYRED,WRK,KFREE,LFREE) 358C 359C Construct the reduced gradient with trial vectors from LURSP3 360C 361 CALL GETGPV(LABEL,FC,FV,CMO,UDV,PV,XINDX,ANTSYM, 362 & WRK(KFREE),LFREE) 363 CALL DCOPY(KZYVAR,WRK(KFREE),1,WRK(KGD),1) 364 REWIND (LURSP3) 365 IF (KOFFTY.EQ.1) THEN 366 READ (LURSP3) 367 END IF 368 DO I = 1, KZRED 369 IF (IBTYP(KOFFTY+I).EQ.JBCNDX) THEN 370 CALL READT(LURSP3,KZCONF,WRK(KBVEC)) 371 PZY = DDOT(KZCONF,WRK(KBVEC),1,WRK(KGD),1) 372 PYZ = DDOT(KZCONF,WRK(KBVEC),1,WRK(KGD+KZVAR),1) 373 ELSE 374 CALL READT(LURSP3,KZYWOP,WRK(KBVEC)) 375 PZY = DDOT(KZWOPT,WRK(KBVEC),1,WRK(KGD+KZCONF),1) 376 & + DDOT(KZWOPT,WRK(KBVEC+KZWOPT),1, 377 & WRK(KGD+KZCONF+KZVAR),1) 378 PYZ = DDOT(KZWOPT,WRK(KBVEC),1,WRK(KGD+KZCONF+KZVAR),1) 379 & + DDOT(KZWOPT,WRK(KBVEC+KZWOPT),1,WRK(KGD+KZCONF),1) 380 END IF 381 REDGD(2*I-1) = PZY 382 REDGD(2*I) = PYZ 383 END DO 384C 385C Unpack the triangular packed reduced E[2] and S[2] 386C 387 CALL DSPTGE(KZYRED,REDE,WRK(KREDE)) 388 CALL DSPTGE(KZYRED,REDS,WRK(KREDS)) 389C 390 IF (IPRABS.GE.10) THEN 391 WRITE(LUPRI,'(2(/,5X,A))') ' Reduced B[1] gradient', 392 & '========================' 393 CALL OUTPUT(REDGD,1,2,1,KZRED,2,KZRED,-1,LUPRI) 394 WRITE(LUPRI,'(2(/,5X,A))') ' Reduced E[2] matrix', 395 & '=====================' 396 CALL OUTPAK(REDE,KZYRED,-1,LUPRI) 397 WRITE(LUPRI,'(2(/,5X,A))') ' Reduced S[2] matrix', 398 & '=====================' 399 CALL OUTPAK(REDS,KZYRED,-1,LUPRI) 400 END IF 401C 402 DO IFREQ = 1,NFREQ_ABS 403C 404C Put the reduced gradient B[1] in a complex form 405C (with zero imaginary part) for later call to linear equation 406C solver ZSYSV, and construct the reduced ( E[2] - {w+iW}*S[2] ) 407C matrix. 408C 409 CALL DZERO(REDZGD,2*KZYRED) 410 CALL DCOPY(KZYRED,REDGD,1,REDZGD,2) 411 CALL DZERO(REDZ,2*MAXRM*MAXRM) 412 CALL DCOPY(KZYRED*KZYRED,WRK(KREDE),1,REDZ,2) 413 CALL DAXPY(KZYRED*KZYRED,-FREQ_ABS(IFREQ),WRK(KREDS),1, 414 & REDZ,2) 415 CALL DAXPY(KZYRED*KZYRED,-DAMPING,WRK(KREDS),1,REDZ(2),2) 416C 417 IF (IPRABS.GT.10) THEN 418 WRITE(LUPRI,'(/,5X,A,/,5X,2(A,D11.4),/,5X,A)') 419 & ' Reduced ( E[2] - {w+iW}*S[2] ) matrix', 420 & ' with w =', FREQ_ABS(IFREQ),' and W =', DAMPING, 421 & '========================================' 422 CALL COUTPUT(REDZ,1,KZYRED,1,KZYRED,KZYRED,KZYRED, 423 & -1,LUPRI) 424 END IF 425C 426 CALL ZSYSV('L',KZYRED,1,REDZ,KZYRED,WRK(KIPIV),REDZGD, 427 & KZYRED,WRK(KFREE),LFREE,INFO) 428C 429C Store solution of response vector in reduced space in REDVEC. 430C Real and imaginary parts of each element in the response vector 431C are stored subsequently. 432C 433 CALL DCOPY(2*KZYRED,REDZGD,1,REDVEC(1,IFREQ),1) 434C 435 IF (IPRABS.GE.5) THEN 436 WRITE(LUPRI,'(/,5X,A,/,5X,2(A,D11.4),/,5X,A)') 437 & ' Reduced ( E[2] - {w+iW}*S[2] )-1 * B[1] vector', 438 & ' with w =', FREQ_ABS(IFREQ),' and W =', DAMPING, 439 & '================================================' 440 CALL COUTPUT(REDVEC(1,IFREQ),1,2,1,KZRED,2,KZRED,-1,LUPRI) 441 END IF 442C 443 IF (LABEL(2:8).EQ.'DIPLEN') THEN 444 ITYPE = 1 445 ELSE IF (LABEL(2:8).EQ.'ANGMOM') THEN 446 ITYPE = 2 447 ELSE 448 WRITE(LUPRI,'(A)') ' Warning: Unknown operator!' 449 GOTO 99 450 END IF 451 CALL DIPLAB(LABEL,INDEX) 452 RESLRF(1,IFREQ+NBATCH*NFREQ_BATCH,INDEX,INDEX,ITYPE) = 453 & DDOT(KZYRED,REDGD,1,REDVEC(1,IFREQ),2) 454 RESLRF(2,IFREQ+NBATCH*NFREQ_BATCH,INDEX,INDEX,ITYPE) = 455 & DDOT(KZYRED,REDGD,1,REDVEC(2,IFREQ),2) 456 99 CONTINUE 457C 458C End of loop over frequencies 459C 460 END DO 461C 462 IF (IPRABS.GE.2) THEN 463 WRITE(LUPRI,'(/3A,/A,I2,A,I4,A,/A,I7,A,//6X,A,/6X,A)') 464 & ' --- Value of linear polarizability for operator ', 465 & LABEL, 'of <<<',' >>> symmetry',KSYMOP, 466 & ' in the reduced space of dimension', 467 & KZYRED, '. The <<<', 468 & ' >>> full variational space has dimension', 469 & KZYVAR,'. ---', 470 & ' Frequency Damping Real Imaginary', 471 & '----------------------------------------------' 472 DO IFREQ=1,NFREQ_ABS 473 WRITE(LUPRI,'(6X,F10.4,F11.6,2F12.6)') FREQ_ABS(IFREQ), 474 & DAMPING, (RESLRF(I,IFREQ+NBATCH*NFREQ_BATCH, 475 & INDEX,INDEX,ITYPE),I=1,2) 476 END DO 477 CALL FLSHFO(LUPRI) 478 END IF 479C 480 RETURN 481 END 482 SUBROUTINE FINDMAXN(X,IXLEN,IMAX,NMAX) 483C Put indices to the NMAX elements of X with largest magnitude 484C in IMAX. 485 IMPLICIT NONE 486 DOUBLE PRECISION X 487 INTEGER IXLEN,IMAX,NMAX,I,J,K,N 488 DIMENSION X(IXLEN),IMAX(NMAX) 489 N = 1 490 IMAX(1) = 1 491 DO I=2,IXLEN 492C Check if X(I) is larger than elements already stored 493 J = N 494 DO WHILE (J.GE.1.AND.ABS(X(I)).GT.ABS(X(IMAX(J)))) 495 J=J-1 496 ENDDO 497C now put I after J is there is space 498 IF (J.LT.NMAX) THEN 499 N = MIN(NMAX,N+1) 500 DO K=N,J+2,-1 501 IMAX(K) = IMAX(K-1) 502 ENDDO 503 IMAX(J+1) = I 504 ENDIF 505 ENDDO 506 NMAX = N 507 END 508 SUBROUTINE ABSRESULT(MJWOP,RESLRF,WRK,LWRK) 509C 510#include "implicit.h" 511#include "priunit.h" 512#include "absorp.h" 513#include "inforb.h" 514#include "codata.h" 515#include "qrinf.h" 516#include "infvar.h" 517#include "inftap.h" 518#include "infdim.h" 519C 520 LOGICAL FOUND, CONV 521C 522 DIMENSION MJWOP(2,MAXWOP,8),RESLRF(2,NFREQ_INTERVAL,3,3,2) 523 DIMENSION WRK(LWRK),ISRCORBR(5),ISRCORBI(5) 524C 525 KFREE = 1 526 LFREE = LWRK 527C 528 CALL TITLER('FINAL RESULTS FOR ABSORPTION','*',120) 529 CALL AROUND('Excitation energies for dipole allowed states') 530 WRITE(LUPRI,'(16X,A,/17X,A4,A6,A21,/16X,A)') 531 & '==========================================', 532 & 'Sym.','State','Excitation energy', 533 & '==========================================' 534 DO ISYM=1,NSYM 535 IF (NOPER(ISYM).GT.0) THEN 536 DO ISTATE=1,NEXCITED_STATES 537 WRITE(LUPRI,'(13X,2I6,F14.6,A,F10.4,A)') 538 & ISYM,ISTATE,EXC_ENERGY(ISTATE,ISYM),' a.u.', 539 & EXC_ENERGY(ISTATE,ISYM)*XTEV,' eV' 540 END DO 541 END IF 542 END DO 543C 544 CALL AROUND('Polarizability with damping') 545 CALL PRSYMB(LUPRI,'-',66,1) 546 WRITE(LUPRI,'(A3,2A10,A12,2A16)') 547 & ' No','A-oper','B-oper','Frequency','Real part','Imag part' 548 CALL PRSYMB(LUPRI,'-',66,1) 549 DO IFREQ=1,NFREQ_INTERVAL 550 IF (ABS_INTERVAL) THEN 551 FREQ = FREQ_INTERVAL(1) + (IFREQ-1)*FREQ_INTERVAL(3) 552 ELSE 553 FREQ = FREQ_ALPHA(IFREQ) 554 END IF 555 WRITE(LUPRI,'(3(/,I3,2A10,F12.6,2F16.6))') 556 & 6*IFREQ-5,'XDIPLEN','XDIPLEN',FREQ, 557 & RESLRF(1,IFREQ,1,1,1),RESLRF(2,IFREQ,1,1,1), 558 & 6*IFREQ-4,'YDIPLEN','YDIPLEN',FREQ, 559 & RESLRF(1,IFREQ,2,2,1),RESLRF(2,IFREQ,2,2,1), 560 & 6*IFREQ-3,'ZDIPLEN','ZDIPLEN',FREQ, 561 & RESLRF(1,IFREQ,3,3,1),RESLRF(2,IFREQ,3,3,1) 562 WRITE(LUPRI,'(3(I3,2A10,F12.6,2F16.6,/))') 563 & 6*IFREQ-2,'XDIPLEN','YDIPLEN',FREQ, 564 & RESLRF(1,IFREQ,1,2,1),RESLRF(2,IFREQ,1,2,1), 565 & 6*IFREQ-1,'XDIPLEN','ZDIPLEN',FREQ, 566 & RESLRF(1,IFREQ,1,3,1),RESLRF(2,IFREQ,1,3,1), 567 & 6*IFREQ,'YDIPLEN','ZDIPLEN',FREQ, 568 & RESLRF(1,IFREQ,2,3,1),RESLRF(2,IFREQ,2,3,1) 569 WRITE(LUPRI,'(6X,A14,3X,F12.6,2F16.6)') 570 & 'Averaged value',FREQ, 571 & (RESLRF(1,IFREQ,1,1,1)+RESLRF(1,IFREQ,2,2,1)+ 572 & RESLRF(1,IFREQ,3,3,1))/3, 573 & (RESLRF(2,IFREQ,1,1,1)+RESLRF(2,IFREQ,2,2,1)+ 574 & RESLRF(2,IFREQ,3,3,1))/3 575 END DO 576 WRITE(LUPRI,'(A)')' ' 577 CALL PRSYMB(LUPRI,'-',66,1) 578 WRITE (LUPRI,'(A,F10.6,A,F7.4,A,F7.1,A)') 579 & ' Damping parameter equals', 580 & DAMPING,' au =', 581 & DAMPING*XTEV,' eV =', 582 & DAMPING*XTKAYS,' cm-1' 583C 584 IF (ABS_BETA) THEN 585 CALL AROUND('First-order hyperpolarizability'// 586 & ' with damping') 587 ELSE IF (ABS_MCD) THEN 588 CALL AROUND('Magnetic circular dichroism'// 589 & ' with damping') 590 END IF 591 IF (ABS_BETA.OR.ABS_MCD) THEN 592 CALL PRSYMB(LUPRI,'-',94,1) 593 WRITE(LUPRI,'(A3,6A10,2A16)') 594 & ' No','A-oper','B-oper','C-oper','A-freq', 595 & 'B-freq','C-freq','Real part','Imag part' 596 CALL PRSYMB(LUPRI,'-',94,1) 597 BTMP = 99.9D9 598 CTMP = 99.9D9 599 BTERM = 0.0 600 RMORD = 0.0 601 DO IQRF=1,NQRF 602 IF (BTMP.NE.QRFFRQ(IQRF,1).OR.CTMP.NE.QRFFRQ(IQRF,2)) THEN 603 IF ((ABS_MCD).AND.(.NOT.IQRF.EQ.1)) THEN 604 WRITE(LUPRI,'(5X,A21,7X,3F10.6,2F16.6)') 605 & 'Orientational average',(QRFFRQ(IQRF-1,I),I=1,3), 606 & BTERM,RMORD 607 BTERM = 0.0 608 RMORD = 0.0 609 END IF 610 WRITE(LUPRI,'(A)') ' ' 611 BTMP = QRFFRQ(IQRF,1) 612 CTMP = QRFFRQ(IQRF,2) 613 END IF 614 IF (ABS_BETA) THEN 615 WRITE(LUPRI,'(I4,3A10,3F10.6,2F16.6)') IQRF, 616 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3), 617 & (RES_BETA(IQRF,I), I=1,2) 618 ELSE IF (ABS_MCD) THEN 619 WRITE(LUPRI,'(I4,3A10,3F10.6,2F16.6)') IQRF, 620 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3), 621 & RES_BETA(IQRF,2),RES_BETA(IQRF,1) 622C 623C Add contribution to orientational average 624C 625 IF ((QRFLAB(IQRF,1)(:1).EQ.'X') 626 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Y') 627 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 628 BTERM = BTERM - RES_BETA(IQRF,2) 629 RMORD = RMORD - RES_BETA(IQRF,1) 630 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 631 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Z') 632 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 633 BTERM = BTERM - RES_BETA(IQRF,2) 634 RMORD = RMORD - RES_BETA(IQRF,1) 635 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 636 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'X') 637 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 638 BTERM = BTERM - RES_BETA(IQRF,2) 639 RMORD = RMORD - RES_BETA(IQRF,1) 640 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 641 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Y') 642 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 643 BTERM = BTERM + RES_BETA(IQRF,2) 644 RMORD = RMORD + RES_BETA(IQRF,1) 645 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'X') 646 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Z') 647 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 648 BTERM = BTERM + RES_BETA(IQRF,2) 649 RMORD = RMORD + RES_BETA(IQRF,1) 650 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 651 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'X') 652 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 653 BTERM = BTERM + RES_BETA(IQRF,2) 654 RMORD = RMORD + RES_BETA(IQRF,1) 655 END IF 656 END IF 657 END DO 658 WRITE(LUPRI,'(5X,A21,7X,3F10.6,2F16.6)') 659 & 'Orientational average',(QRFFRQ(NQRF,I), I=1,3), 660 & BTERM,RMORD 661 662 WRITE(LUPRI,'(A)')' ' 663 CALL PRSYMB(LUPRI,'-',94,1) 664 WRITE (LUPRI,'(A,F10.6,A,F7.4,A,F7.1,A)') 665 & ' Damping parameter equals', 666 & DAMPING,' au =', 667 & DAMPING*XTEV,' eV =', 668 & DAMPING*XTKAYS,' cm-1' 669 END IF 670C 671 IF (ABS_ANALYZE) THEN 672 MZYVMX = 2*NVARMA 673 CALL MEMGET('REAL',KVEC,2*MZYVMX,WRK,KFREE,LFREE) 674 CALL AROUND('Analysis of response vectors') 675 DO ISYM=1,NSYM 676 IF (NOPER(ISYM).GT.0) THEN 677 KZYVAR=MZYVAR(ISYM) 678 DO IOPER=1,NOPER(ISYM) 679 DO IFREQ=1,NFREQ_INTERVAL 680 IF (ABS_INTERVAL) THEN 681 FREQ = FREQ_INTERVAL(1) + (IFREQ-1)*FREQ_INTERVAL(3) 682 ELSE 683 FREQ = FREQ_ALPHA(IFREQ) 684 END IF 685 WRITE(LUPRI,'(/A11,A10,2(/,A11,F10.6))') 'Property :', 686 & LABOP(IOPER,ISYM), 687 & 'Frequency:',ABS(FREQ), 688 & 'Damping :',DAMPING 689 CALL REARSP(LURSP,2*KZYVAR,WRK(KVEC), 690 & LABOP(IOPER,ISYM),'COMPLEX ', 691 & ABS(FREQ),0.0D0,ISYM,0,THCLR_ABSORP, 692 & FOUND,CONV,ANTSYM) 693 IF (.NOT. FOUND) THEN 694 WRITE(LUPRI,'(/A)') 695 & ' Response vector not found on file RSPVEC' 696 CALL QUIT('Response vector not found on file RSPVEC') 697 ELSE IF(.NOT. CONV) THEN 698 WRITE (LUPRI,'(/A)') ' @WARNING: '// 699 & 'Response vector not converged on file RSPVEC' 700 END IF 701C 702 DNORM_RE=DNRM2(KZYVAR,WRK(KVEC),1) 703 DNORM_IM=DNRM2(KZYVAR,WRK(KVEC+KZYVAR),1) 704C 705 IF (IPRABS.GT.1) THEN 706 WRITE(LUPRI,'(/A,2F14.8)') 707 & ' Norm of vector (real and imag):', 708 & DNORM_RE, DNORM_IM 709C 710 WRITE(LUPRI,'(/7X,2A5,2(2X,2(5X,A2,5X)))') 711 & 'occ','vir','ZR','YR','ZI','YI' 712 DO I=1,KZYVAR/2 713 WRITE(LUPRI,'(I5,2X,2I5,2(2X,2F12.8))') 714 & I,MJWOP(1,I,ISYM),MJWOP(2,I,ISYM), 715 & WRK(KVEC-1+I),WRK(KVEC+KZYVAR/2-1+I), 716 & WRK(KVEC+KZYVAR-1+I),WRK(KVEC+3*KZYVAR/2-1+I) 717 END DO 718C 719 END IF 720 721 WRITE(LUPRI,'(/A,2F14.8)') 722 & ' Norm of vector (real and imag):', 723 & DNORM_RE, DNORM_IM 724C Sum occupied orbital contributions into aggregates 725 MAXOCC=0 726 DO I=1,KZYVAR/2 727 MAXOCC = MAX(MAXOCC,MJWOP(1,I,ISYM)) 728 ENDDO 729 CALL MEMGET('REAL',KORBWR,MAXOCC,WRK,KFREE,LFREE) 730 CALL MEMGET('REAL',KORBWI,MAXOCC,WRK,KFREE,LFREE) 731c real part 732 DO I=1,MAXOCC 733 WRK(KORBWR+I-1) = 0.0D0 734 ENDDO 735 DO I=1,KZYVAR/2 736 WRK(KORBWR+MJWOP(1,I,ISYM)-1) = 737 & WRK(KORBWR+MJWOP(1,I,ISYM)-1) + 738 & WRK(KVEC-1+I)**2 + WRK(KVEC+KZYVAR/2-1+I)**2 739 ENDDO 740c imag part 741 DO I=1,MAXOCC 742 WRK(KORBWI+I-1) = 0.0D0 743 ENDDO 744 DO I=1,KZYVAR/2 745 WRK(KORBWI+MJWOP(1,I,ISYM)-1) = 746 & WRK(KORBWI+MJWOP(1,I,ISYM)-1) + 747 & WRK(KVEC+KZYVAR-1+I)**2 + 748 & WRK(KVEC+3*KZYVAR/2-1+I)**2 749 ENDDO 750 NSRC=5 751 CALL FINDMAXN(WRK(KORBWR),MAXOCC,ISRCORBR,NSRC) 752 NSRC=5 753 CALL FINDMAXN(WRK(KORBWI),MAXOCC,ISRCORBI,NSRC) 754 WRITE (LUPRI,*) 755 & 'Important occupied orbital contributions (normalized)' 756 WRITE (LUPRI,*) 757 & ' occ real occ imag' 758 DO I=1,NSRC 759 WRITE (LUPRI,'(I5,F7.3,I6,F7.3)') ISRCORBR(I), 760 & WRK(KORBWR+ISRCORBR(I)-1)/DNORM_RE**2, 761 & ISRCORBI(I), 762 & WRK(KORBWI+ISRCORBI(I)-1)/DNORM_IM**2 763 ENDDO 764 END DO 765 END DO 766 END IF 767 END DO 768 END IF 769C 770 RETURN 771 END 772 SUBROUTINE ABSRESID(IOP,LABEL,NFREQ_ABS,FREQ_ABS, 773 & CONVERGED,REDVEC,IBTYP,FC,CMO,UDV,PV,XINDX,WRK,LWRK) 774C 775C PURPOSE: 776C Compute the complex residual vector to the linear response 777C equation. 778C 779C R = B[1] - ( E[2] - {w+iW}*S[2] )*(NR + iNI) 780C 781C or equivalently the real and imaginary parts 782C 783C RR = B[1] - E[2]*NR + S[2]*( w*NR - W*NI ) 784C RI = - E[2]*NI + S[2]*( w*NI + W*NR ) 785C 786C From file we read: 787C LURSP3 contains trial vectors 788C LURSP5 contains sigma vectors 789C 790C / b1 \ / b2 \ 791C NR = sum(i=1,KZRED) | | * REDVEC(4i-3) + | | * REDVEC(4i-1) 792C \ b2 /_i \ b1 /_i 793C 794C / b1 \ / b2 \ 795C NI = sum(i=1,KZRED) | | * REDVEC(4i-2) + | | * REDVEC(4i) 796C \ b2 /_i \ b1 /_i 797C 798#include "implicit.h" 799#include "priunit.h" 800#include "dummy.h" 801#include "infrsp.h" 802#include "wrkrsp.h" 803#include "absorp.h" 804#include "inftap.h" 805#include "ibndxdef.h" 806C 807 LOGICAL CONVERGED 808 CHARACTER*8 LABEL 809 DIMENSION REDVEC(2*MAXRM,MXFREQ),IBTYP(MAXRM),FREQ_ABS(MXFREQ), 810 & CMO(*),UDV(*),PV(*),FC(*),XINDX(*),WRK(LWRK) 811C 812 KFREE = 1 813 LFREE = LWRK 814 CALL MEMGET('REAL',KRR,KZYVAR,WRK,KFREE,LFREE) 815 CALL MEMGET('REAL',KRI,KZYVAR,WRK,KFREE,LFREE) 816 CALL MEMGET('REAL',KBVEC,KZYVAR,WRK,KFREE,LFREE) 817 CALL MEMGET('REAL',KSLI,KZYVAR,WRK,KFREE,LFREE) 818C 819 CONVERGED = .TRUE. 820C 821C Get the gradient 822C 823 DO IFREQ=1,NFREQ_ABS 824C 825C Initialize residual vector 826C 827 CALL GETGPV(LABEL,FC,DUMMY,CMO,UDV,PV,XINDX,ANTSYM, 828 & WRK(KFREE),LFREE) 829 CALL DCOPY(KZYVAR,WRK(KFREE),1,WRK(KRR),1) 830 CALL DZERO(WRK(KRI),KZYVAR) 831C 832 IF (IPRABS.GE.5) THEN 833 WRITE(LUPRI,'(2(/,5X,A))') ' B[1] gradient', 834 & '========================' 835 CALL OUTPUT(WRK(KRR),1,KZVAR,1,2,KZVAR,2,1,LUPRI) 836 END IF 837C 838C Read trial and sigma vectors and add contributions to residual vector 839C 840 REWIND (LURSP3) 841 REWIND (LURSP5) 842 IF (KOFFTY.EQ.1) THEN 843 READ (LURSP3) 844 READ (LURSP5) 845 END IF 846 DO I = 1, KZRED 847C 848 FR1 = REDVEC(4*I-3,IFREQ) 849 FR2 = REDVEC(4*I-1,IFREQ) 850 FI1 = REDVEC(4*I-2,IFREQ) 851 FI2 = REDVEC(4*I,IFREQ) 852C 853 FR3 = (FREQ_ABS(IFREQ)*FR1-DAMPING*FI1) 854 FR4 = (FREQ_ABS(IFREQ)*FR2-DAMPING*FI2) 855 FI3 = (FREQ_ABS(IFREQ)*FI1+DAMPING*FR1) 856 FI4 = (FREQ_ABS(IFREQ)*FI2+DAMPING*FR2) 857C 858C Sigma vectors used that equal E[2]*N 859C 860 CALL READT(LURSP5,KZYVAR,WRK(KBVEC)) 861C 862 CALL DAXPY(KZYVAR,-FR1,WRK(KBVEC),1,WRK(KRR),1) 863 CALL DAXPY(KZVAR,-FR2,WRK(KBVEC+KZVAR),1,WRK(KRR),1) 864 CALL DAXPY(KZVAR,-FR2,WRK(KBVEC),1,WRK(KRR+KZVAR),1) 865 CALL DAXPY(KZYVAR,-FI1,WRK(KBVEC),1,WRK(KRI),1) 866 CALL DAXPY(KZVAR,-FI2,WRK(KBVEC+KZVAR),1,WRK(KRI),1) 867 CALL DAXPY(KZVAR,-FI2,WRK(KBVEC),1,WRK(KRI+KZVAR),1) 868C 869C Trial vectors used to perform S[2]*N 870C 871 IF (IBTYP(KOFFTY+I).EQ.JBCNDX) THEN 872 CALL READT(LURSP3,KZCONF,WRK(KBVEC)) 873 CALL DZERO(WRK(KSLI),KZYVAR) 874 CALL RSPSLI(1,0,WRK(KBVEC),DUMMY, 875 & UDV,WRK(KSLI),XINDX,WRK(KFREE),LFREE) 876 ELSE 877 CALL READT(LURSP3,KZYWOP,WRK(KBVEC)) 878 CALL DZERO(WRK(KSLI),KZYVAR) 879 CALL RSPSLI(0,1,DUMMY,WRK(KBVEC), 880 & UDV,WRK(KSLI),XINDX,WRK(KFREE),LFREE) 881 END IF 882C 883 CALL DAXPY(KZYVAR,FR3,WRK(KSLI),1,WRK(KRR),1) 884 CALL DAXPY(KZVAR,-FR4,WRK(KSLI+KZVAR),1,WRK(KRR),1) 885 CALL DAXPY(KZVAR,-FR4,WRK(KSLI),1,WRK(KRR+KZVAR),1) 886 CALL DAXPY(KZYVAR,FI3,WRK(KSLI),1,WRK(KRI),1) 887 CALL DAXPY(KZVAR,-FI4,WRK(KSLI+KZVAR),1,WRK(KRI),1) 888 CALL DAXPY(KZVAR,-FI4,WRK(KSLI),1,WRK(KRI+KZVAR),1) 889C 890C End loop over trial and sigma vectors 891C 892 END DO 893C 894 DNORM_RR = DNRM2(KZYVAR,WRK(KRR),1) 895 DNORM_RI = DNRM2(KZYVAR,WRK(KRI),1) 896 DNORM_RT = SQRT(DNORM_RR**2 + DNORM_RI**2) 897C 898 DNORM_NR = DNRM2(KZYRED,REDVEC(1,IFREQ),2) 899 DNORM_NI = DNRM2(KZYRED,REDVEC(2,IFREQ),2) 900 DNORM_NT = SQRT( DNORM_NR**2 + DNORM_NI**2 ) 901C 902 RESID(1,IFREQ,IOP,KSYMOP) = DNORM_RR/DNORM_NT 903 RESID(2,IFREQ,IOP,KSYMOP) = DNORM_RI/DNORM_NT 904 RESID(3,IFREQ,IOP,KSYMOP) = DNORM_RT/DNORM_NT 905C 906 IF (DNORM_RT/DNORM_NT.GT.THCLR_ABSORP) CONVERGED=.FALSE. 907C 908C End loop over frequencies 909C 910 END DO 911C 912C Print norm of residual vector 913C 914 IF (IPRABS.GE.2) THEN 915 WRITE(LUPRI,'(/,1X,A,1P,D9.2,/,1X,A,/,1X,A,I4,A)') 916 & 'Convergence of RSP solution vectors, threshold =', 917 & THCLR_ABSORP, 918 & '-------------------------------------------------------------', 919 & '(dimension of reduced space:',KZYRED,')' 920 DO IFREQ = 1,NFREQ_ABS 921 WRITE(LUPRI,'(1X,A,F7.4,3X,A,F10.6,3X,A,1P,3D9.2)') 922 & 'Frequency:',FREQ_ABS(IFREQ), 923 & 'Damping:',DAMPING, 924 & 'Residual:',(RESID(I,IFREQ,IOP,KSYMOP),I=1,3) 925 END DO 926 END IF 927C 928 RETURN 929 END 930 SUBROUTINE ABSLR(IOP,LABEL,NFREQ_ABS,FREQ_ABS, 931 & REDE,REDS,REDZ,REDGD,REDZGD,IBTYP, 932 & REDVEC,CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) 933C 934#include "implicit.h" 935#include "priunit.h" 936#include "dummy.h" 937#include "infrsp.h" 938#include "wrkrsp.h" 939#include "absorp.h" 940#include "inftap.h" 941#include "ibndxdef.h" 942C 943C PURPOSE: 944C Solve the coupled equations in this iteration 945C 946C (i) ( E[2] - w*S[2] ) NR = B[1] - W*S[2]*NI 947C (ii) ( E[2] - w*S[2] ) NI = W*S[2]*NR 948C 949 CHARACTER*8 LABEL,BLANK 950 PARAMETER (BLANK=' ', D0=0.0D0,DSQRT2=1.4142135623731D0) 951 DIMENSION REDE(MAXRM*MAXRM),REDS(MAXRM*MAXRM),REDZ(2*MAXRM*MAXRM), 952 & REDGD(MAXRM),REDZGD(2*MAXRM),REDVEC(2*MAXRM,MXFREQ), 953 & IBTYP(MAXRM),FREQ_ABS(MXFREQ) 954 DIMENSION CMO(*),UDV(*),PV(*),FOCK(*),FC(*),FCAC(*),H2AC(*) 955 DIMENSION XINDX(*),WRK(LWRK) 956C 957 KFREE = 1 958 LFREE = LWRK 959 CALL MEMGET('REAL',KGD1,KZYVAR,WRK,KFREE,LFREE) 960 CALL MEMGET('REAL',KGD2,KZYVAR,WRK,KFREE,LFREE) 961 CALL MEMGET('REAL',KSLI,KZYVAR,WRK,KFREE,LFREE) 962 CALL MEMGET('REAL',KBVEC,KZYVAR,WRK,KFREE,LFREE) 963C 964 DO IFREQ=1,NFREQ_ABS 965C 966 IF (RESID(3,IFREQ,IOP,KSYMOP).LT.THCLR_ABSORP) THEN 967 IF (IPRABS.GE.2) WRITE(LUPRI,'(/,1X,A,F8.4,A)') 968 & '----- Frequency:',FREQ_ABS(IFREQ),' converged -----' 969 GO TO 100 970 END IF 971C 972C Construct the right-hand-side of equations (i) and (ii) 973C 974 CALL GETGPV(LABEL,FC,DUMMY,CMO,UDV,PV,XINDX,ANTSYM, 975 & WRK(KFREE),LFREE) 976 CALL DCOPY(KZYVAR,WRK(KFREE),1,WRK(KGD1),1) 977 CALL DZERO(WRK(KGD2),KZYVAR) 978C 979C Read trial vectors and add contributions to gradient 980C 981 REWIND (LURSP3) 982 IF (KOFFTY.EQ.1) THEN 983 READ (LURSP3) 984 END IF 985 DO I = 1, KZRED 986 FR1 = DAMPING * REDVEC(4*I-3,IFREQ) 987 FR2 = DAMPING * REDVEC(4*I-1,IFREQ) 988 FI1 = DAMPING * REDVEC(4*I-2,IFREQ) 989 FI2 = DAMPING * REDVEC(4*I,IFREQ) 990 IF (IBTYP(KOFFTY+I).EQ.JBCNDX) THEN 991 CALL READT(LURSP3,KZCONF,WRK(KBVEC)) 992 CALL DZERO(WRK(KSLI),KZYVAR) 993 CALL RSPSLI(1,0,WRK(KBVEC),DUMMY, 994 & UDV,WRK(KSLI),XINDX,WRK(KFREE),LFREE) 995 CALL DAXPY(KZCONF,-FI1,WRK(KSLI),1,WRK(KGD1),1) 996 CALL DAXPY(KZCONF,FI2,WRK(KSLI),1,WRK(KGD1+KZVAR),1) 997 CALL DAXPY(KZCONF,FR1,WRK(KSLI),1,WRK(KGD2),1) 998 CALL DAXPY(KZCONF,-FR2,WRK(KSLI),1,WRK(KGD2+KZVAR),1) 999 ELSE 1000 CALL READT(LURSP3,KZYWOP,WRK(KBVEC)) 1001 CALL DZERO(WRK(KSLI),KZYVAR) 1002 CALL RSPSLI(0,1,DUMMY,WRK(KBVEC), 1003 & UDV,WRK(KSLI),XINDX,WRK(KFREE),LFREE) 1004 CALL DAXPY(KZWOPT,-FI1,WRK(KSLI),1,WRK(KGD1+KZCONF),1) 1005 CALL DAXPY(KZWOPT,-FI1,WRK(KSLI+KZWOPT),1, 1006 & WRK(KGD1+KZCONF+KZVAR),1) 1007 CALL DAXPY(KZWOPT,FI2,WRK(KSLI+KZWOPT),1, 1008 & WRK(KGD1+KZCONF),1) 1009 CALL DAXPY(KZWOPT,FI2,WRK(KSLI),1, 1010 & WRK(KGD1+KZCONF+KZVAR),1) 1011C 1012 CALL DAXPY(KZWOPT,FR1,WRK(KSLI),1,WRK(KGD2+KZCONF),1) 1013 CALL DAXPY(KZWOPT,FR1,WRK(KSLI+KZWOPT),1, 1014 & WRK(KGD2+KZCONF+KZVAR),1) 1015 CALL DAXPY(KZWOPT,-FR2,WRK(KSLI+KZWOPT),1, 1016 & WRK(KGD2+KZCONF),1) 1017 CALL DAXPY(KZWOPT,-FR2,WRK(KSLI),1, 1018 & WRK(KGD2+KZCONF+KZVAR),1) 1019 END IF 1020 END DO 1021C 1022C Get norm of solution vectors and gradients 1023C 1024 DNORM_NR = DNRM2(KZYRED,REDVEC(1,IFREQ),2) 1025 DNORM_NI = DNRM2(KZYRED,REDVEC(2,IFREQ),2) 1026 DNORM_NT = SQRT( DNORM_NR**2 + DNORM_NI**2 ) 1027 DNORM_GD1 = DNRM2(KZYVAR,WRK(KGD1),1) 1028 DNORM_GD2 = DNRM2(KZYVAR,WRK(KGD2),1) 1029C 1030 IF (IPRABS.GE.2) WRITE(LUPRI,'(/,1X,A,F7.4,4(/,1X,A,1P,D9.2))') 1031 & 'Frequency = ',FREQ_ABS(IFREQ), 1032 & 'Norm of gradient vector 1 =',DNORM_GD1, 1033 & 'Norm of solution vector 1 =',DNORM_NR, 1034 & 'Norm of gradient vector 2 =',DNORM_GD2, 1035 & 'Norm of solution vector 2 =',DNORM_NI 1036C 1037C Solve Eqs. (i) and (ii) 1038C 1039 MAXIT = MAX_MACRO 1040 KEXCNV = 1 1041 KEXSTV = KEXCNV 1042 KEXSIM = KEXCNV 1043 RESTLR = .TRUE. 1044C 1045 IF (RESID(1,IFREQ,IOP,KSYMOP).GT.RESID(2,IFREQ,IOP,KSYMOP))THEN 1046C 1047 IF (RESID(1,IFREQ,IOP,KSYMOP).LT.THCLR_ABSORP/DSQRT2) THEN 1048 IF (IPRABS.GE.2) WRITE(LUPRI,'(/,1X,A,1P,D9.2,A)') 1049 & '----- Skip real equation because residual is only:', 1050 & RESID(1,IFREQ,IOP,KSYMOP), ' -----' 1051 ELSE 1052 IF (IPRABS.GE.2) WRITE(LUPRI,'(/,/,1X,A,F8.4,A)') 1053 & '----- Solving real equation for frequency:', 1054 & FREQ_ABS(IFREQ),' ------' 1055 THCRSP = 0.1D0*THCLR_ABSORP*DNORM_NT/DNORM_NR/DSQRT2 1056 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 1057 & .TRUE.,BLANK,BLANK,WRK(KGD1),REDGD,REDE,REDS, 1058 & IBTYP,FREQ_ABS(IFREQ),RESIDUAL,REDVEC,XINDX, 1059 & WRK(KFREE),LFREE) 1060 END IF 1061C 1062 ELSE 1063C 1064 IF (RESID(2,IFREQ,IOP,KSYMOP).LT.THCLR_ABSORP/DSQRT2) THEN 1065 IF (IPRABS.GE.2) WRITE(LUPRI,'(/,1X,A,1P,D9.2,A)') 1066 & '----- Skip imaginary equation because residual is only:', 1067 & RESID(2,IFREQ,IOP,KSYMOP), ' -----' 1068 ELSE 1069 IF (IPRABS.GE.2) WRITE(LUPRI,'(/,/,1X,A,F8.4,A)') 1070 & '----- Solving imaginary equation for frequency:', 1071 & FREQ_ABS(IFREQ),' ------' 1072 THCRSP = 0.5D0*THCLR_ABSORP*DNORM_NT/DNORM_NI/DSQRT2 1073 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 1074 & .TRUE.,BLANK,BLANK,WRK(KGD2),REDGD,REDE,REDS, 1075 & IBTYP,FREQ_ABS(IFREQ),RESIDUAL,REDVEC,XINDX, 1076 & WRK(KFREE),LFREE) 1077 END IF 1078C 1079 END IF 1080C 1081C End of loop over frequencies 1082C 1083 100 CONTINUE 1084 END DO 1085C 1086 RETURN 1087 END 1088 SUBROUTINE ABSWRT(LUIN,LUOUT,LABEL,IOP,NFREQ_ABS,FREQ_ABS, 1089 & IBTYP,REDVEC,WRK,LWRK) 1090C 1091C PURPOSE: 1092C 1093C If LUIN=LURSP3, construct and write response vectors N to file 1094C If LUIN=LURSP5, construct and write vectors E[2]*N to file 1095C 1096C From file we read: 1097C LUIN contains trial (LUIN=LURSP3) vectors or 1098C sigma vectors (LUIN=LURSP5). Output vectors are written to LUOUT 1099C which normally is LURSP, apart from when we perform a reduction of 1100C the dimension of the reduced space. 1101C 1102C / b1 \ / b2 \ 1103C NR = sum(i=1,KZRED) | | * REDVEC(4i-3) + | | * REDVEC(4i-1) 1104C \ b2 /_i \ b1 /_i 1105C 1106C / b1 \ / b2 \ 1107C NI = sum(i=1,KZRED) | | * REDVEC(4i-2) + | | * REDVEC(4i) 1108C \ b2 /_i \ b1 /_i 1109C 1110#include "implicit.h" 1111#include "priunit.h" 1112#include "infrsp.h" 1113#include "wrkrsp.h" 1114#include "absorp.h" 1115#include "ibndxdef.h" 1116C 1117 CHARACTER*8 LABEL,BLANK 1118 PARAMETER (D0 = 0.0D0, D1=1.0D0) 1119 DIMENSION REDVEC(2*MAXRM,MXFREQ),IBTYP(MAXRM), 1120 & FREQ_ABS(MXFREQ),WRK(LWRK) 1121C 1122C Do not allocate with MEMGET since that will give some extra 1123C bytes of information in between real and imaginary parts 1124C 1125 KBVEC = 1 1126 KNR = KBVEC + KZYVAR 1127 KNI = KNR + KZYVAR 1128 KFREE = KNI + KZYVAR 1129 LFREE = LWRK - KFREE 1130 IF (LFREE.LT.0) CALL ERRWRK('ABSWRT',KFREE,LWRK) 1131C 1132 DO IFREQ=1,NFREQ_ABS 1133C 1134C Read trial vectors and add contributions to response vector 1135C 1136 REWIND (LUIN) 1137 IF (KOFFTY.EQ.1) THEN 1138 READ (LUIN) 1139 END IF 1140C 1141 CALL DZERO(WRK(KNR),2*KZYVAR) 1142 DO I = 1, KZRED 1143 FR1 = REDVEC(4*I-3,IFREQ) 1144 FR2 = REDVEC(4*I-1,IFREQ) 1145 FI1 = REDVEC(4*I-2,IFREQ) 1146 FI2 = REDVEC(4*I,IFREQ) 1147C 1148 IF (IBTYP(KOFFTY+I).EQ.JBCNDX) THEN 1149 CALL READT(LUIN,KZCONF,WRK(KBVEC)) 1150 CALL DAXPY(KZCONF,FR1,WRK(KBVEC),1,WRK(KNR),1) 1151 CALL DAXPY(KZCONF,FR2,WRK(KBVEC),1,WRK(KNR+KZVAR),1) 1152 CALL DAXPY(KZCONF,FI1,WRK(KBVEC),1,WRK(KNI),1) 1153 CALL DAXPY(KZCONF,FI2,WRK(KBVEC),1,WRK(KNI+KZVAR),1) 1154 ELSE 1155 CALL READT(LUIN,KZYWOP,WRK(KBVEC)) 1156C 1157 CALL DAXPY(KZWOPT,FR1,WRK(KBVEC),1,WRK(KNR+KZCONF),1) 1158 CALL DAXPY(KZWOPT,FR1,WRK(KBVEC+KZWOPT),1, 1159 & WRK(KNR+KZCONF+KZVAR),1) 1160 CALL DAXPY(KZWOPT,FR2,WRK(KBVEC+KZWOPT),1, 1161 & WRK(KNR+KZCONF),1) 1162 CALL DAXPY(KZWOPT,FR2,WRK(KBVEC),1, 1163 & WRK(KNR+KZCONF+KZVAR),1) 1164C 1165 CALL DAXPY(KZWOPT,FI1,WRK(KBVEC),1,WRK(KNI+KZCONF),1) 1166 CALL DAXPY(KZWOPT,FI1,WRK(KBVEC+KZWOPT),1, 1167 & WRK(KNI+KZCONF+KZVAR),1) 1168 CALL DAXPY(KZWOPT,FI2,WRK(KBVEC+KZWOPT),1, 1169 & WRK(KNI+KZCONF),1) 1170 CALL DAXPY(KZWOPT,FI2,WRK(KBVEC),1, 1171 & WRK(KNI+KZCONF+KZVAR),1) 1172 END IF 1173 END DO 1174C 1175C Write response vectors to file 1176C 1177 CALL WRTRSP(LUOUT,2*KZYVAR,WRK(KNR),LABEL,'COMPLEX ', 1178 & FREQ_ABS(IFREQ),D0,KSYMOP,0, 1179 & RESID(3,IFREQ,IOP,KSYMOP),D1) 1180C 1181 IF (IPRABS.GE.10) THEN 1182 WRITE(LUPRI,'(/A,A10,A,I4,/A,F12.8,A,I4)') 1183 & ' Response vector in ABSWRT for operator',LABEL, 1184 & ' of symmetry', 1185 & KSYMOP,' and frequency',FREQ_ABS(IFREQ), 1186 & ' is written to file LU=',LUOUT 1187 CALL PRSYMB(LUPRI,'=',72,1) 1188 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 1189 CALL OUTPUT(WRK(KNR),1,KZYVAR/2,1,4,KZYVAR/2,4,1,LUPRI) 1190 END IF 1191C 1192 END DO 1193C 1194 RETURN 1195 END 1196 SUBROUTINE ABSCTL(IOPER,LABEL, 1197 & REDE,REDS,REDZ,REDGD,REDZGD,REDVEC,IBTYP, 1198 & CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC,XINDX, 1199 & RESLRF,WRK,LWRK) 1200C 1201#include "implicit.h" 1202#include "priunit.h" 1203#include "pgroup.h" 1204#include "infrsp.h" 1205#include "wrkrsp.h" 1206#include "absorp.h" 1207#include "inftap.h" 1208#include "ibndxdef.h" 1209C 1210C PURPOSE: 1211C 1212C Solve the complex LR equation 1213C 1214C ( E[2] - {w+iW}*S[2] )* (NR + iNI) = B[1] 1215C 1216C or equivalently the pair of real equations 1217C 1218C ( E[2] - w*S[2] )* NR = B[1] - W*S[2]*NI 1219C ( E[2] - w*S[2] )* NI = W*S[2]*NR 1220C 1221 LOGICAL CONVERGED 1222 CHARACTER LABEL*8 1223 DIMENSION FREQ_ABS(MXFREQ) 1224 DIMENSION REDE(MAXRM*MAXRM),REDS(MAXRM*MAXRM),REDZ(2*MAXRM*MAXRM), 1225 & REDGD(MAXRM),REDZGD(2*MAXRM),REDVEC(2*MAXRM,MXFREQ), 1226 & IBTYP(MAXRM) 1227 DIMENSION CMO(*),UDV(*),PV(*),FOCK(*),FC(*),FV(*),FCAC(*),H2AC(*) 1228 DIMENSION XINDX(*),RESLRF(2,NFREQ_INTERVAL,3,3,2),WRK(LWRK) 1229C 1230 KFREE = 1 1231 LFREE = LWRK 1232C 1233C Copy the frequencies into local variables to be used in the solver. 1234C One can thereby choose to change the frequencies for certain operators. 1235C 1236 IF (LABEL(2:8).EQ.'ANGMOM') THEN 1237 NFREQ_ABS = 1 1238 FREQ_ABS(1) = 0.0D0 1239 ELSE 1240 NFREQ_ABS = NFREQ_ALPHA 1241 CALL DCOPY(NFREQ_ALPHA,FREQ_ALPHA,1,FREQ_ABS,1) 1242 END IF 1243C 1244C Check if this an absorption calculation over a specified frequency 1245C interval in which case we divide the interval in batches of freqs. 1246C 1247C NFREQ_BATCH = number of frequencies per batch 1248C 1249 IBATCH = 0 1250 10 CONTINUE 1251C 1252 IF (ABS_INTERVAL) THEN 1253 NFREQ_ABS=MIN(NFREQ_BATCH,NFREQ_INTERVAL-IBATCH*NFREQ_BATCH) 1254 DO I=1,NFREQ_ABS 1255 FREQ_ABS(I)=FREQ_INTERVAL(1) + 1256 & (I-1+IBATCH*NFREQ_BATCH)*FREQ_INTERVAL(3) 1257 END DO 1258 END IF 1259C 1260 IF (IPRABS.GE.0) THEN 1261 WRITE(LUPRI,'(/A/3A,I2,3A/A,5F12.8/,(11X,5F12.8))') 1262 & ' ABSVEC1 -- Solving linear response equations'// 1263 & ' ( E[2] - {w-iW}*S[2] ) N = B[1] ', 1264 & ' ABSVEC1 -- for operator ', LABEL, ' of symmetry', 1265 & KSYMOP,' ( ',REP(KSYMOP-1),') at the frequencies:', 1266 & ' ABSVEC1 --',(FREQ_ABS(I),I=1,NFREQ_ABS) 1267 END IF 1268C 1269C Check if response solution vectors already exist on file unit 1270C LURSP with name RSPVEC. 1271C 1272 CALL CHKONFILE(LURSP,CONVERGED,LABEL,KSYMOP, 1273 & NFREQ_ABS,FREQ_ABS,THCLR_ABSORP,WRK(KFREE)) 1274 IF (CONVERGED) THEN 1275 WRITE(LUPRI,'(/A,I3)') 1276 & ' ABSCTL: converged response vectors found on file'// 1277 & ' RSPVEC with unit=',LURSP 1278 CALL GETLRF(LURSP,LABEL,NFREQ_ABS,FREQ_ABS, 1279 & FC,FV,CMO,UDV,PV,XINDX,IBATCH,RESLRF, 1280 & WRK(KFREE),LFREE) 1281 GOTO 990 1282 END IF 1283C 1284C Reduce the dimension of the reduced space 1285C 1286 IF (ABS_REDUCE) THEN 1287 CALL ABSREDUC(LABEL,NFREQ_ABS,FREQ_ABS, 1288 & REDE,REDS,REDZ,REDGD,REDZGD,IBTYP, 1289 & REDVEC,FC,FV,CMO,UDV,PV,XINDX,IBATCH,RESLRF, 1290 & WRK(KFREE),LFREE) 1291 CALL REDSPACE(LABEL,IOPER,NFREQ_ABS,FREQ_ABS, 1292 & REDE,REDS,REDVEC,IBTYP, 1293 & FC,FV,CMO,UDV,PV,XINDX,WRK(KFREE),LFREE) 1294 END IF 1295C 1296 ITER = 0 1297 100 CONTINUE 1298C 1299 IF (IPRABS.GE.2) THEN 1300 WRITE(LUPRI,'(/A,/A,I3,A,/A)') 1301 & ' =======================================', 1302 & ' ---- Macro iteration number',ITER,' ------', 1303 & ' =======================================' 1304C 1305 END IF 1306C 1307C Solve for resonse vector in reduced space 1308C 1309 CALL ABSREDUC(LABEL,NFREQ_ABS,FREQ_ABS, 1310 & REDE,REDS,REDZ,REDGD,REDZGD,IBTYP, 1311 & REDVEC,FC,FV,CMO,UDV,PV,XINDX,IBATCH,RESLRF, 1312 & WRK(KFREE),LFREE) 1313C 1314C Compute the residual in this iteration 1315C 1316 CALL ABSRESID(IOPER,LABEL,NFREQ_ABS,FREQ_ABS, 1317 & CONVERGED,REDVEC,IBTYP,FC,CMO,UDV,PV,XINDX,WRK(KFREE),LFREE) 1318C 1319C Check status 1320C 1321 IF (CONVERGED) GOTO 900 1322 IF (ITER.GE.MAX_MACRO) GOTO 910 1323C 1324C Not converged, expand reduced space by solving coupled LR equations 1325C 1326 CALL ABSLR(IOPER,LABEL,NFREQ_ABS,FREQ_ABS, 1327 & REDE,REDS,REDZ,REDGD,REDZGD, 1328 & IBTYP,REDVEC,CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC,XINDX, 1329 & WRK(KFREE),LFREE) 1330C 1331C Iteration completed 1332C 1333 ITER = ITER + 1 1334 GOTO 100 1335C 1336C Final messages from ABSCTL 1337C 1338 900 CONTINUE 1339 IF (IPRABS.GE.0) THEN 1340 WRITE(LUPRI,'(/,1X,A,I4,A)') 1341 & '*** ABSCTL: THE REQUESTED',NFREQ_ABS, 1342 & ' SOLUTION VECTORS CONVERGED' 1343 END IF 1344 GOTO 920 1345 910 CONTINUE 1346 WRITE(LUPRI,'(/A,I4,A/A)') 1347 & ' *** ABSCTL WARNING: THE REQUESTED',NFREQ_ABS, 1348 & ' SOLUTION VECTORS NOT CONVERGED', 1349 & ' --- MAXIMUM NUMBER OF ITERATIONS REACHED ---' 1350 GOTO 920 1351C 1352C Construct and write response vectors to file 1353C 1354 920 CONTINUE 1355 CALL ABSWRT(LURSP3,LURSP,LABEL,IOPER,NFREQ_ABS,FREQ_ABS, 1356 & IBTYP,REDVEC,WRK(KFREE),LFREE) 1357 CALL GETLRF(LURSP,LABEL,NFREQ_ABS,FREQ_ABS, 1358 & FC,FV,CMO,UDV,PV,XINDX,IBATCH,RESLRF, 1359 & WRK(KFREE),LFREE) 1360C 1361 IF (IPRABS.GE.0) THEN 1362 WRITE(LUPRI,'(2(/,1X,A),1P,D9.2,/,1X,A,/,1X,A,I4,A)') 1363 & 'RSP solution vectors written to file.', 1364 & 'Convergence of RSP solution vectors, threshold =', 1365 & THCLR_ABSORP, 1366 & '-------------------------------------------------------------', 1367 & '(dimension of reduced space:',KZYRED,')' 1368 DO IFREQ = 1,NFREQ_ABS 1369 WRITE(LUPRI,'(1X,A,F7.4,4X,A,F10.6,4X,A,1P,D9.2)') 1370 & 'Frequency:',FREQ_ABS(IFREQ), 1371 & 'Damping:',DAMPING, 1372 & 'Residual:',RESID(3,IFREQ,IOPER,KSYMOP) 1373 END DO 1374 WRITE(LUPRI,'(/2A,3(/A))') 1375 & '@ Value of linear response result for operator: ', 1376 & LABEL, 1377 & '@ --------------------------------------------', 1378 & '@ Frequency Damping Real Imaginary', 1379 & '@ --------------------------------------------' 1380 IF (LABEL(2:8).EQ.'DIPLEN') THEN 1381 ITYPE = 1 1382 ELSE IF (LABEL(2:8).EQ.'ANGMOM') THEN 1383 ITYPE = 2 1384 ELSE 1385 WRITE(LUPRI,'(A)') ' Warning: Unknown operator!' 1386 GOTO 99 1387 END IF 1388 CALL DIPLAB(LABEL,INDEX) 1389 DO IFREQ=1,NFREQ_ABS 1390 WRITE(LUPRI,'(A,F10.4,F11.6,2F12.6)') '@', 1391 & FREQ_ABS(IFREQ), DAMPING, 1392 & (RESLRF(I,IFREQ+IBATCH*NFREQ_BATCH, 1393 & INDEX,INDEX,ITYPE),I=1,2) 1394 END DO 1395 99 CONTINUE 1396 END IF 1397C 1398C If this an absorption calculation over a specified frequency 1399C interval then we need to do another batch of freqs. 1400C 1401 990 CONTINUE 1402 IBATCH = IBATCH + 1 1403 IF (ABS_INTERVAL .AND. 1404 & NFREQ_INTERVAL-IBATCH*NFREQ_BATCH.GT.0) THEN 1405 GOTO 10 1406 END IF 1407C 1408 RETURN 1409 END 1410 SUBROUTINE GETLRF(LU,LABEL,NFREQ_ABS,FREQ_ABS, 1411 & FC,FV,CMO,UDV,PV,XINDX,NBATCH,RESLRF, 1412 & WRK,LWRK) 1413#include "implicit.h" 1414#include "priunit.h" 1415#include "wrkrsp.h" 1416#include "absorp.h" 1417C 1418 CHARACTER*8 LABEL,LABGD 1419 LOGICAL FOUND,CONV 1420 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),XINDX(*), 1421 & WRK(LWRK),RESLRF(2,NFREQ_INTERVAL,3,3,2), 1422 & FREQ_ABS(MXFREQ) 1423C 1424 KFREE = 1 1425 LFREE = LWRK 1426C 1427 CALL DIPLAB(LABEL,INDEX) 1428 IF (LABEL(2:8).EQ.'DIPLEN') THEN 1429 ITYPE = 1 1430 ELSE IF (LABEL(2:8).EQ.'ANGMOM') THEN 1431 ITYPE = 2 1432 ELSE 1433 WRITE(LUPRI,'(A)') ' Warning: Unknown operator!' 1434 RETURN 1435 END IF 1436C 1437 CALL MEMGET('REAL',KVEC,2*KZYVAR,WRK,KFREE,LFREE) 1438 CALL MEMGET('REAL',KGD,KZYVAR,WRK,KFREE,LFREE) 1439C 1440 DO IOPER=1,NOPER(KSYMOP) 1441 LABGD = LABOP(IOPER,KSYMOP) 1442 IF (LABEL(2:8).EQ.LABGD(2:8)) THEN 1443 CALL DIPLAB(LABGD,INXGD) 1444 CALL GETGPV(LABGD,FC,FV,CMO,UDV,PV,XINDX,ANTSYM, 1445 & WRK(KGD),LFREE) 1446 DO IFREQ=1,NFREQ_ABS 1447 CALL REARSP(LU,2*KZYVAR,WRK(KVEC),LABEL,'COMPLEX ', 1448 & ABS(FREQ_ABS(IFREQ)),0.0D0,KSYMOP,0,THCLR_ABSORP, 1449 & FOUND,CONV,ANTSYM) 1450 RESLRF(1,IFREQ+NBATCH*NFREQ_BATCH,INXGD,INDEX,ITYPE) = 1451 & DDOT(KZYVAR,WRK(KGD),1,WRK(KVEC),1) 1452 RESLRF(2,IFREQ+NBATCH*NFREQ_BATCH,INXGD,INDEX,ITYPE) = 1453 & DDOT(KZYVAR,WRK(KGD),1,WRK(KVEC+KZYVAR),1) 1454 END DO 1455 END IF 1456 END DO 1457C 1458 RETURN 1459 END 1460 SUBROUTINE CHKONFILE(LU,FOUND,LABEL,ISYM,NFREQ_ABS,FREQ_ABS, 1461 & THD,FLAGS) 1462#include "priunit.h" 1463C 1464 LOGICAL FOUND,FLAGS(NFREQ_ABS) 1465 CHARACTER*8 LABEL,LAB1,LAB2 1466 INTEGER LU,ISYM,NFREQ_ABS,ISYM1,ISYM2,LEN,NBAS,NORB 1467 REAL*8 FREQ_ABS(NFREQ_ABS),THD,FREQ1,FREQ2,ANTSYM,RSD,EMCSCF 1468C 1469 FOUND = .TRUE. 1470 DO I=1,NFREQ_ABS 1471 FLAGS(I)=.FALSE. 1472 END DO 1473C 1474 REWIND(LU) 1475 READ(LU) 1476 100 READ(LU,END=900,ERR=900) LAB1,LAB2,FREQ1,FREQ2,ISYM1,ISYM2, 1477 & ANTSYM,RSD,LEN,EMCSCF,NBAS,NORB 1478C 1479 IF (LAB1.NE.LABEL .OR. LAB2 .NE.'COMPLEX ' .OR. 1480 & ISYM1.NE.ISYM .OR. ISYM2.NE.0 .OR. 1481 & RSD .GT.THD .OR. FREQ2.NE.0.0D0) GOTO 100 1482C 1483 DO I=1,NFREQ_ABS 1484 IF (FREQ1.EQ.FREQ_ABS(I)) THEN 1485 FLAGS(I)=.TRUE. 1486 END IF 1487 END DO 1488C 1489 READ(LU,END=900,ERR=900) 1490C 1491 GOTO 100 1492C 1493 900 CONTINUE 1494 DO I=1,NFREQ_ABS 1495 FOUND = FOUND .AND. FLAGS(I) 1496 END DO 1497 RETURN 1498 END 1499 SUBROUTINE REDSPACE(LABEL,IOPER,NFREQ_ABS,FREQ_ABS, 1500 & REDE,REDS,REDVEC,IBTYP,FC,FV,CMO,UDV,PV,XINDX,WRK,LWRK) 1501#include "implicit.h" 1502#include "priunit.h" 1503#include "infrsp.h" 1504#include "wrkrsp.h" 1505#include "absorp.h" 1506#include "inftap.h" 1507#include "inforb.h" 1508#include "ibndxdef.h" 1509C 1510 CHARACTER*8 LABEL 1511 LOGICAL FOUND,CONV,READFLAG 1512 DIMENSION REDVEC(2*MAXRM,MXFREQ),IBTYP(MAXRM), 1513 & REDE(MAXRM*MAXRM),REDS(MAXRM*MAXRM), 1514 & FREQ_ABS(MXFREQ),WRK(LWRK) 1515 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),XINDX(*) 1516C 1517 IF (IPRABS.GE.2) WRITE(LUPRI,'(2(/A))') 1518 & ' INFO: lowering dimension of reduced space' 1519 CALL FLSHFO(LUPRI) 1520C 1521C Reduction of dimension in reduced space only implented for SCF 1522C type wave functions. 1523C 1524 IF (KZCONF .GT. 0) RETURN 1525C 1526 KBVEC = 1 1527 KSVEC = KBVEC + KZYVAR 1528 KB = KSVEC + KZYVAR 1529 KS = KB + 2*KZYVAR 1530 KFREE = KS + 2*KZYVAR 1531 LFREE = LWRK - KFREE 1532 IF (LFREE.LT.0) CALL ERRWRK('ABSWRT',KFREE,LWRK) 1533C 1534C Write response vectors N and sigma vectors E[2]*N to respective 1535C temporary files LUBTMP and LUSTMP. Unit numbers are set in GPOPEN. 1536C 1537 LUBTMP=-1 1538 LUSTMP=-1 1539 CALL GPOPEN(LUBTMP,'RSPBTMP','NEW','SEQUENTIAL','UNFORMATTED', 1540 & 0,.FALSE.) 1541 CALL GPOPEN(LUSTMP,'RSPSTMP','NEW','SEQUENTIAL','UNFORMATTED', 1542 & 0,.FALSE.) 1543 WRITE(LUBTMP) NISH,NASH,NORB,NBAS,NSYM 1544 WRITE(LUBTMP) 'EOFLABEL' 1545 WRITE(LUSTMP) NISH,NASH,NORB,NBAS,NSYM 1546 WRITE(LUSTMP) 'EOFLABEL' 1547 CALL ABSWRT(LURSP3,LUBTMP,LABEL,IOPER,NFREQ_ABS,FREQ_ABS, 1548 & IBTYP,REDVEC,WRK(KFREE),LFREE) 1549 CALL ABSWRT(LURSP5,LUSTMP,LABEL,IOPER,NFREQ_ABS,FREQ_ABS, 1550 & IBTYP,REDVEC,WRK(KFREE),LFREE) 1551C 1552C Read response and sigma vectors from temporary files and perform 1553C modified Gram-Schmidt orthogonalization of the respective real and 1554C imaginary parts of the response vectors (these then form the new 1555C set of trial vectors which are written to LURSP3). The 1556C corresponding sigma vectors are contructed and written to LURSP5. 1557C We keep NBOLD of the old trial vectors in order to avoid linear 1558C dependence. 1559C 1560 NBOLD = 0 1561 NLINDEP = 0 1562 IFREQ=0 1563 READFLAG = .TRUE. 1564 DO I=1+NBOLD,2*NFREQ_ABS+NBOLD 1565 IF (READFLAG) THEN 1566 IFREQ=IFREQ+1 1567 CALL REARSP(LUBTMP,LENGTH,WRK(KB),LABEL,'COMPLEX ', 1568 & ABS(FREQ_ABS(IFREQ)),0.0D0, 1569 & KSYMOP,0,THCLR,FOUND,CONV,ANTSYM) 1570 IF (LENGTH.NE.2*KZYVAR) CALL QUIT(' REDSPACE: read error') 1571 CALL REARSP(LUSTMP,LENGTH,WRK(KS),LABEL,'COMPLEX ', 1572 & ABS(FREQ_ABS(IFREQ)),0.0D0, 1573 & KSYMOP,0,THCLR,FOUND,CONV,ANTSYM) 1574 IF (LENGTH.NE.2*KZYVAR) CALL QUIT(' REDSPACE: read error') 1575 CALL DCOPY(KZYVAR,WRK(KB),1,WRK(KBVEC),1) 1576 CALL DCOPY(KZYVAR,WRK(KS),1,WRK(KSVEC),1) 1577 READFLAG = .FALSE. 1578C 1579 IF (IPRABS.GE.2) THEN 1580 CALL GETGPV(LABEL,FC,FV,CMO,UDV,PV,XINDX,ANTSYM, 1581 & WRK(KFREE),LFREE) 1582 VAL1 = DDOT(KZYVAR,WRK(KFREE),1,WRK(KB),1) 1583 VAL2 = DDOT(KZYVAR,WRK(KFREE),1,WRK(KB+KZYVAR),1) 1584 WRITE(LUPRI,'(2A,F10.6,2F16.10)') 1585 & ' Value of linear response for ', 1586 & LABEL,FREQ_ABS(IFREQ),VAL1,VAL2 1587 END IF 1588C 1589 IF (IPRABS.GE.10) THEN 1590 WRITE(LUPRI,'(/A,A10,A,I4,/A,F12.8,A,I4)') 1591 & ' Response vector in REDSPACE for operator',LABEL, 1592 & ' of symmetry', 1593 & KSYMOP,' and frequency',FREQ_ABS(IFREQ), 1594 & ' is read from file LUBTMP=',LUBTMP 1595 CALL PRSYMB(LUPRI,'=',72,1) 1596 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 1597 CALL OUTPUT(WRK(KB),1,KZVAR,1,4,KZVAR,4,1,LUPRI) 1598 WRITE(LUPRI,'(/A,A10,A,I4,/A,F12.8,A,I4)') 1599 & ' Sigma vector in REDSPACE for operator',LABEL, 1600 & ' of symmetry', 1601 & KSYMOP,' and frequency',FREQ_ABS(IFREQ), 1602 & ' is read from file LUSTMP=',LUSTMP 1603 CALL PRSYMB(LUPRI,'=',72,1) 1604 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 1605 CALL OUTPUT(WRK(KS),1,KZVAR,1,4,KZVAR,4,1,LUPRI) 1606 END IF 1607 ELSE 1608 CALL DCOPY(KZYVAR,WRK(KB+KZYVAR),1,WRK(KBVEC),1) 1609 CALL DCOPY(KZYVAR,WRK(KS+KZYVAR),1,WRK(KSVEC),1) 1610 READFLAG = .TRUE. 1611 END IF 1612C 1613C First KZYVAR elements of the space for complex response vectors 1614C are now available as scratch. 1615C 1616 KBTMP=KB 1617 KSTMP=KS 1618C 1619 REWIND (LURSP3) 1620 REWIND (LURSP5) 1621C 1622 DO J=1,I-1-NLINDEP 1623 CALL READT(LURSP3,KZYVAR,WRK(KBTMP)) 1624 CALL READT(LURSP5,KZYVAR,WRK(KSTMP)) 1625C Orthogonalize against trial vectors (Z,Y) 1626 FAC = DDOT(KZYVAR,WRK(KBTMP),1,WRK(KBVEC),1) 1627 CALL DAXPY(KZYVAR,-FAC,WRK(KBTMP),1,WRK(KBVEC),1) 1628 CALL DAXPY(KZYVAR,-FAC,WRK(KSTMP),1,WRK(KSVEC),1) 1629C Orthogonalize against trial vectors (Y,Z) 1630 FAC = DDOT(KZVAR,WRK(KBTMP+KZVAR),1,WRK(KBVEC),1) + 1631 & DDOT(KZVAR,WRK(KBTMP),1,WRK(KBVEC+KZVAR),1) 1632 CALL DAXPY(KZVAR,-FAC,WRK(KBTMP+KZVAR),1,WRK(KBVEC),1) 1633 CALL DAXPY(KZVAR,-FAC,WRK(KBTMP),1,WRK(KBVEC+KZVAR),1) 1634 CALL DAXPY(KZVAR,-FAC,WRK(KSTMP+KZVAR),1,WRK(KSVEC),1) 1635 CALL DAXPY(KZVAR,-FAC,WRK(KSTMP),1,WRK(KSVEC+KZVAR),1) 1636 END DO 1637C 1638C Check linear dependence. 1639C 1640 BNORM = DNRM2(KZYVAR,WRK(KBVEC),1) 1641 IF (BNORM.LT.1.0D-6) THEN 1642 IF (IPRABS.GE.0) WRITE(LUPRI,'(A)') 1643 & ' INFO: vector not used as trial vector'// 1644 & ' due to linear dependence' 1645 NLINDEP = NLINDEP + 1 1646 ELSE 1647 FAC = 1.0D0/BNORM 1648 CALL DSCAL(KZYVAR,FAC,WRK(KBVEC),1) 1649 CALL DSCAL(KZYVAR,FAC,WRK(KSVEC),1) 1650C Orthogonalize trial vector (Z,Y) against (Y,Z) using 1651C symemtric orthogonalization to keep paired structure 1652 OVLP=2.0D0*DDOT(KZVAR,WRK(KBVEC),1,WRK(KBVEC+KZVAR),1) 1653 C1=0.5D0*(1/SQRT(1+OVLP) + 1/SQRT(1-OVLP)) 1654 C2=0.5D0*(1/SQRT(1+OVLP) - 1/SQRT(1-OVLP)) 1655 CALL DCOPY(KZYVAR,WRK(KBVEC),1,WRK(KBTMP),1) 1656 CALL DCOPY(KZYVAR,WRK(KSVEC),1,WRK(KSTMP),1) 1657 CALL DSCAL(KZYVAR,C1,WRK(KBVEC),1) 1658 CALL DSCAL(KZYVAR,C1,WRK(KSVEC),1) 1659 CALL DAXPY(KZVAR,C2,WRK(KBTMP),1,WRK(KBVEC+KZVAR),1) 1660 CALL DAXPY(KZVAR,C2,WRK(KBTMP+KZVAR),1,WRK(KBVEC),1) 1661 CALL DAXPY(KZVAR,C2,WRK(KSTMP),1,WRK(KSVEC+KZVAR),1) 1662 CALL DAXPY(KZVAR,C2,WRK(KSTMP+KZVAR),1,WRK(KSVEC),1) 1663C 1664 CALL WRITT(LURSP3,KZYVAR,WRK(KBVEC)) 1665 CALL WRITT(LURSP5,KZYVAR,WRK(KSVEC)) 1666 IF (IPRABS.GE.10) THEN 1667 WRITE(LUPRI,'(/A,I4,A,I2)') 1668 & ' INFO: Normalized trial vector',I, 1669 & ' is written to file LURSP3=',LURSP3 1670 IF (IPRABS.GE.0) THEN 1671 CALL PRSYMB(LUPRI,'=',72,1) 1672 WRITE(LUPRI,'(4X,2A15)') 'Z','Y' 1673 CALL OUTPUT(WRK(KBVEC),1,KZVAR,1,2, 1674 & KZVAR,2,1,LUPRI) 1675 CALL FLSHFO(LUPRI) 1676 END IF 1677 END IF 1678 END IF 1679C 1680 END DO 1681C 1682 CALL GPCLOSE(LUBTMP,'DELETE') 1683 CALL GPCLOSE(LUSTMP,'DELETE') 1684C 1685C Construct the reduced E[2] and S[2] matrices in triangular packed 1686C format. 1687C 1688 KZRED = 2*NFREQ_ABS+NBOLD-NLINDEP 1689 KZYRED = 2*KZRED 1690C 1691 DO I=1,KZRED 1692 IROW=2*I-1 1693 JOFF1 = IROW*(IROW-1)/2 1694 JOFF2 = JOFF1+IROW 1695 REWIND(LURSP3) 1696 DO IDUM=1,I 1697 CALL READT(LURSP3,KZYVAR,WRK(KBVEC)) 1698 END DO 1699 REWIND(LURSP3) 1700 REWIND(LURSP5) 1701 DO J=1,I 1702 CALL READT(LURSP3,KZYVAR,WRK(KB)) 1703 CALL READT(LURSP5,KZYVAR,WRK(KS)) 1704 ICOL=2*J-1 1705 REDE(JOFF1+ICOL) = 1706 & DDOT(KZYVAR,WRK(KBVEC),1,WRK(KS),1) 1707 IF (I.NE.J) REDE(JOFF1+ICOL+1) = 1708 & DDOT(KZVAR,WRK(KBVEC),1,WRK(KS+KZVAR),1) + 1709 & DDOT(KZVAR,WRK(KBVEC+KZVAR),1,WRK(KS),1) 1710 REDE(JOFF2+ICOL) = 1711 & DDOT(KZVAR,WRK(KBVEC+KZVAR),1,WRK(KS),1) + 1712 & DDOT(KZVAR,WRK(KBVEC),1,WRK(KS+KZVAR),1) 1713 REDE(JOFF2+ICOL+1) = 1714 & DDOT(KZVAR,WRK(KBVEC+KZVAR),1,WRK(KS+KZVAR),1) + 1715 & DDOT(KZVAR,WRK(KBVEC),1,WRK(KS),1) 1716 CALL DSCAL(KZVAR,-1.0D0,WRK(KB+KZVAR),1) 1717 REDS(JOFF1+ICOL) = 1718 & DDOT(KZYVAR,WRK(KBVEC),1,WRK(KB),1) 1719 REDS(JOFF2+ICOL) = 1720 & DDOT(KZVAR,WRK(KBVEC+KZVAR),1,WRK(KB),1) + 1721 & DDOT(KZVAR,WRK(KBVEC),1,WRK(KB+KZVAR),1) 1722 CALL DSCAL(KZYVAR,-1.0D0,WRK(KB),1) 1723 IF (I.NE.J) REDS(JOFF1+ICOL+1) = 1724 & DDOT(KZVAR,WRK(KBVEC),1,WRK(KB+KZVAR),1) + 1725 & DDOT(KZVAR,WRK(KBVEC+KZVAR),1,WRK(KB),1) 1726 REDS(JOFF2+ICOL+1) = 1727 & DDOT(KZVAR,WRK(KBVEC+KZVAR),1,WRK(KB+KZVAR),1) + 1728 & DDOT(KZVAR,WRK(KBVEC),1,WRK(KB),1) 1729 END DO 1730 END DO 1731 CALL DSCAL(KZYRED*(KZYRED+1)/2,2.0D0,REDS,1) 1732C 1733 IF (IPRABS.GE.10) THEN 1734 WRITE(LUPRI,'(/A)') 'Reduced E[2] after reduction = ' 1735 CALL OUTPAK(REDE,KZYRED,1,LUPRI) 1736 WRITE(LUPRI,'(/A)') 'Reduced S[2] after reduction = ' 1737 CALL OUTPAK(REDS,KZYRED,1,LUPRI) 1738 CALL FLSHFO(LUPRI) 1739 END IF 1740C 1741 RETURN 1742 END 1743 SUBROUTINE BETA_SETUP 1744#include "implicit.h" 1745#include "priunit.h" 1746#include "inforb.h" 1747#include "absorp.h" 1748C 1749 PARAMETER (THRZERO = 1.0D-6) 1750C 1751 LOGICAL DOHYP 1752 CHARACTER*8 ALAB,BLAB,CLAB 1753C 1754 NQRF = 0 1755 NFREQ_ALPHA = 0 1756C 1757 DO ICFR=1,NFREQ_BETA_C 1758 DO IBFR=1,NFREQ_BETA_B 1759 FREQB = FREQ_BETA_B(IBFR) 1760 FREQC = FREQ_BETA_C(ICFR) 1761 FREQA = -(FREQB+FREQC) 1762 IF (ABS_SHG .AND. .NOT.FREQB.EQ.FREQC) GOTO 100 1763 DO ISYMA=1,NSYM 1764 DO ISYMB=1,NSYM 1765 ISYMC = MULD2H(ISYMA,ISYMB) 1766 IF (NOPER(ISYMA).GT.0 .AND. NOPER(ISYMB).GT.0 .AND. 1767 & NOPER(ISYMC).GT.0) THEN 1768 DO IAOP=1,NOPER(ISYMA) 1769 DO IBOP=1,NOPER(ISYMB) 1770 DO ICOP=1,NOPER(ISYMC) 1771 ALAB = LABOP(IAOP,ISYMA) 1772 BLAB = LABOP(IBOP,ISYMB) 1773 CLAB = LABOP(ICOP,ISYMC) 1774 CALL QRCHK(DOHYP,ALAB,BLAB,CLAB, 1775 & ISYMA,ISYMB,ISYMC,FREQA,FREQB,FREQC) 1776 IF (DOHYP) THEN 1777 NQRF = NQRF +1 1778 QRFLAB(NQRF,1) = ALAB 1779 QRFLAB(NQRF,2) = BLAB 1780 QRFLAB(NQRF,3) = CLAB 1781 QRFSYM(NQRF,1) = ISYMA 1782 QRFSYM(NQRF,2) = ISYMB 1783 QRFSYM(NQRF,3) = ISYMC 1784 QRFFRQ(NQRF,1) = FREQA 1785 QRFFRQ(NQRF,2) = FREQB 1786 QRFFRQ(NQRF,3) = FREQC 1787 END IF 1788 END DO 1789 END DO 1790 END DO 1791 END IF 1792 END DO 1793 END DO 1794 100 CONTINUE 1795 END DO 1796 END DO 1797C 1798 CALL GPDSRT(NFREQ_ALPHA,FREQ_ALPHA,THRZERO) 1799C 1800 IF (IPRABS.GE.0) THEN 1801 CALL AROUND('Setup of Hyperpolarizability Calculation') 1802 WRITE (LUPRI,'(2(/A),I4,A)') 1803 & ' This calculations requires the solution of linear response', 1804 & ' equations for electric dipole operators at',NFREQ_ALPHA, 1805 & ' frequencies:' 1806 WRITE(LUPRI,'(/A,5(4F12.8,/,9X))') 1807 & ' LR FREQ:',(FREQ_ALPHA(I),I=1,NFREQ_ALPHA) 1808 WRITE (LUPRI,'(/A,I4,A)') 1809 & ' and the evaluation of',NQRF,' quadratic response functions:' 1810 WRITE(LUPRI,'(/2A,/A3,6A12,/2A)') 1811 & '--------------------------------------', 1812 & '-------------------------------------', 1813 & ' No','A-oper','B-oper','C-oper', 1814 & 'A-freq','B-freq','C-freq', 1815 & '--------------------------------------', 1816 & '-------------------------------------' 1817 DO IQRF=1,NQRF 1818 WRITE(LUPRI,'(I4,3A12,3F12.8)') IQRF, 1819 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3) 1820 END DO 1821 WRITE(LUPRI,'(2A)') 1822 & '--------------------------------------', 1823 & '-------------------------------------' 1824 END IF 1825C 1826C End of BETA_SETUP 1827C 1828 RETURN 1829 END 1830 SUBROUTINE GPDSRT(N,X,THD) 1831#include "implicit.h" 1832C 1833C======================================================================== 1834C 1835C Purpose: Sort the elements of a vector X containing N double-precision 1836C real numbers. Remove repeated occurrences of elements which 1837C are seperated by less than THD. 1838C 1839C In: X - vector of double-precision real numbers 1840C N - number of elements in vector 1841C THD - threshold for two elements being considered equal 1842C 1843C Out: X - vector containing unique elements, now sorted 1844C N - number of unique elements 1845C 1846C Author: Patrick Norman, 2003. 1847C 1848C======================================================================== 1849C 1850 DIMENSION X(N) 1851C 1852 IF (N.LE.1) RETURN 1853C 1854 M = N 1855 J = 1 1856C 1857 100 CONTINUE 1858 I = J+1 1859 200 CONTINUE 1860C 1861 IF (ABS(X(J)-X(I)).LT.THD) THEN 1862 X(I) = X(M) 1863 M = M-1 1864 IF (I.LE.M) GOTO 200 1865 END IF 1866 IF (X(I).LT.X(J)) THEN 1867 T = X(J) 1868 X(J) = X(I) 1869 X(I) = T 1870 END IF 1871C 1872 I = I+1 1873 IF (I.LE.M) GOTO 200 1874C 1875 J = J+1 1876 IF (J.LT.M ) GOTO 100 1877C 1878 N = M 1879C 1880 RETURN 1881 END 1882 SUBROUTINE ABSQRF(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 1883 & XINDX,MJWOP,WRK,LWRK) 1884C 1885C PURPOSE: 1886C Compute the quadratic response functions with response vectors 1887C found on file. 1888C 1889#include "implicit.h" 1890#include "priunit.h" 1891#include "absorp.h" 1892#include "abslrs.h" 1893#include "wrkrsp.h" 1894#include "infvar.h" 1895#include "maxaqn.h" 1896#include "maxorb.h" 1897#include "mxcent.h" 1898#include "symmet.h" 1899#include "cbilrs.h" 1900#include "infrsp.h" 1901#include "inforb.h" 1902#include "qrinf.h" 1903#include "infdim.h" 1904C 1905 PARAMETER ( D0=0.0D0 ) 1906 DIMENSION CMO(*),UDV(*),PV(*),FOCK(*),FC(*),FV(*),FCAC(*),H2AC(*) 1907 DIMENSION XINDX(*),MJWOP(2,MAXWOP,8),WRK(LWRK) 1908C 1909 DIMENSION HYPVAL(2),TMPVAL(2),VAL(4) 1910 CHARACTER*8 ALAB,BLAB,CLAB 1911 LOGICAL DOCAL 1912C 1913 KFREE = 1 1914 LFREE = LWRK 1915C 1916 ISPINA = 0 1917 ISPINB = 0 1918 ISPINC = 0 1919C 1920 MZYVMX = 2*NVARMA 1921 CALL MEMGET('REAL',KVECA,2*MZYVMX,WRK,KFREE,LFREE) 1922 CALL MEMGET('REAL',KVECB,2*MZYVMX,WRK,KFREE,LFREE) 1923 CALL MEMGET('REAL',KVECC,2*MZYVMX,WRK,KFREE,LFREE) 1924 CALL MEMGET('REAL',KGRAD,2*MZYVMX,WRK,KFREE,LFREE) 1925C 1926 LURSPRES = -1 1927 CALL GPOPEN(LURSPRES,'RESULTS.RSP','UNKNOWN','SEQUENTIAL', 1928 & 'FORMATTED',0,.FALSE.) 1929 DO IQRF=1,NQRF 1930C 1931 HYPVAL(1) = D0 1932 HYPVAL(2) = D0 1933C 1934 ALAB = QRFLAB(IQRF,1) 1935 BLAB = QRFLAB(IQRF,2) 1936 CLAB = QRFLAB(IQRF,3) 1937 ISYMA = QRFSYM(IQRF,1) 1938 ISYMB = QRFSYM(IQRF,2) 1939 ISYMC = QRFSYM(IQRF,3) 1940 FREQA = QRFFRQ(IQRF,1) 1941 FREQB = QRFFRQ(IQRF,2) 1942 FREQC = QRFFRQ(IQRF,3) 1943C 1944C Check if we perhaps have this result on file already 1945C K.Ruud, back programming again in august 2011 1946C 1947 CALL ABCCHK(DOCAL,LURSPRES,ISYMA,ISYMB,ISYMC,FREQA,FREQB,FREQC, 1948 & ALAB,BLAB,CLAB,HYPVAL) 1949 IF (DOCAL) THEN 1950C 1951 WRITE(LUPRI,'(/,2(A,I5),A,3(/A,A10,I5,F10.6))') 1952 & ' Quadratic response function no',IQRF,' of',NQRF,'.', 1953 & ' A operator, symmetry, frequency: ',ALAB,ISYMA,FREQA, 1954 & ' B operator, symmetry, frequency: ',BLAB,ISYMB,FREQB, 1955 & ' C operator, symmetry, frequency: ',CLAB,ISYMC,FREQC 1956C 1957 IF (ABSLRS.AND.(NCONF.LE.1)) THEN 1958 CALL ABS_QRRDVE(ISYMA,ISYMB,ISYMC,ALAB,BLAB,CLAB, 1959 & FREQA,FREQB,FREQC,KZYVA,KZYVB,KZYVC,WRK(KVECA),WRK(KVECB), 1960 & WRK(KVECC)) 1961 ELSE 1962 CALL QRRDVE(ISYMA,ISYMB,ISYMC,ALAB,BLAB,CLAB, 1963 & FREQA,FREQB,FREQC,KZYVA,KZYVB,KZYVC,WRK(KVECA),WRK(KVECB), 1964 & WRK(KVECC)) 1965 ENDIF 1966C 1967C Calculate Na B[2] Nc type terms (two permutations) 1968C 1969 CALL DZERO(WRK(KGRAD),2*MZYVMX) 1970 CALL X2INIT(1,KZYVA,KZYVC,ISYMA,ISPINA,ISYMC,ISPINC,WRK(KVECA), 1971 & WRK(KVECC),WRK(KGRAD),XINDX,UDV,PV,BLAB,ISYMB,ISPINB, 1972 & CMO,MJWOP,WRK(KFREE),LFREE) 1973 VAL(1) = DDOT(KZYVA,WRK(KGRAD),1,WRK(KVECA),1) 1974 VAL(2) = DDOT(KZYVA,WRK(KGRAD),1,WRK(KVECA+KZYVA),1) 1975 TMPVAL(1) = VAL(1) 1976 TMPVAL(2) = VAL(2) 1977 HYPVAL(1) = HYPVAL(1) + VAL(1) 1978 HYPVAL(2) = HYPVAL(2) + VAL(2) 1979C 1980 CALL DZERO(WRK(KGRAD),2*MZYVMX) 1981 CALL X2INIT(1,KZYVA,KZYVB,ISYMA,ISPINA,ISYMB,ISPINB,WRK(KVECA), 1982 & WRK(KVECB),WRK(KGRAD),XINDX,UDV,PV,CLAB,ISYMC,ISPINC, 1983 & CMO,MJWOP,WRK(KFREE),LFREE) 1984 VAL(1) = DDOT(KZYVA,WRK(KGRAD),1,WRK(KVECA),1) 1985 VAL(2) = DDOT(KZYVA,WRK(KGRAD),1,WRK(KVECA+KZYVA),1) 1986 TMPVAL(1) = TMPVAL(1) + VAL(1) 1987 TMPVAL(2) = TMPVAL(2) + VAL(2) 1988 HYPVAL(1) = HYPVAL(1) + VAL(1) 1989 HYPVAL(2) = HYPVAL(2) + VAL(2) 1990C 1991 IF (DAMPING.GT.D0) THEN 1992 CALL DZERO(WRK(KGRAD),2*MZYVMX) 1993 CALL X2INIT(1,KZYVA,KZYVC,ISYMA,ISPINA,ISYMC,ISPINC, 1994 & WRK(KVECA),WRK(KVECC+KZYVC),WRK(KGRAD),XINDX,UDV,PV, 1995 & BLAB,ISYMB,ISPINB,CMO,MJWOP,WRK(KFREE),LFREE) 1996 VAL(1) = -DDOT(KZYVA,WRK(KGRAD),1,WRK(KVECA+KZYVA),1) 1997 VAL(2) = DDOT(KZYVA,WRK(KGRAD),1,WRK(KVECA),1) 1998 TMPVAL(1) = TMPVAL(1) + VAL(1) 1999 TMPVAL(2) = TMPVAL(2) + VAL(2) 2000 HYPVAL(1) = HYPVAL(1) + VAL(1) 2001 HYPVAL(2) = HYPVAL(2) + VAL(2) 2002C 2003 CALL DZERO(WRK(KGRAD),2*MZYVMX) 2004 CALL X2INIT(1,KZYVA,KZYVB,ISYMA,ISPINA,ISYMB,ISPINB, 2005 & WRK(KVECA),WRK(KVECB+KZYVB),WRK(KGRAD),XINDX,UDV,PV, 2006 & CLAB,ISYMC,ISPINC,CMO,MJWOP,WRK(KFREE),LFREE) 2007 VAL(1) = -DDOT(KZYVA,WRK(KGRAD),1,WRK(KVECA+KZYVA),1) 2008 VAL(2) = DDOT(KZYVA,WRK(KGRAD),1,WRK(KVECA),1) 2009 TMPVAL(1) = TMPVAL(1) + VAL(1) 2010 TMPVAL(2) = TMPVAL(2) + VAL(2) 2011 HYPVAL(1) = HYPVAL(1) + VAL(1) 2012 HYPVAL(2) = HYPVAL(2) + VAL(2) 2013 END IF 2014C 2015 IF (IPRABS.GE.0) WRITE(LUPRI,'(A15,4F15.6)') 2016 & ' Na X[2] Ny ',TMPVAL(1),TMPVAL(2),HYPVAL(1),HYPVAL(2) 2017C 2018C Calculate Nb A[2] Nc type terms (two permutations) 2019C 2020 CALL DZERO(WRK(KGRAD),2*MZYVMX) 2021 CALL A2INIT(1,KZYVB,KZYVC,ISYMB,ISPINB,ISYMC,ISPINC,1, 2022 & WRK(KVECC),WRK(KGRAD),XINDX,UDV,PV,ALAB,ISYMA,ISPINA, 2023 & CMO,MJWOP,WRK(KFREE),LFREE) 2024 VAL(1) = DDOT(KZYVB,WRK(KGRAD),1,WRK(KVECB),1) 2025 VAL(2) = DDOT(KZYVB,WRK(KGRAD),1,WRK(KVECB+KZYVB),1) 2026 TMPVAL(1) = VAL(1) 2027 TMPVAL(2) = VAL(2) 2028 HYPVAL(1) = HYPVAL(1) + VAL(1) 2029 HYPVAL(2) = HYPVAL(2) + VAL(2) 2030C 2031 CALL DZERO(WRK(KGRAD),2*MZYVMX) 2032 CALL A2INIT(1,KZYVC,KZYVB,ISYMC,ISPINC,ISYMB,ISPINB,1, 2033 & WRK(KVECB),WRK(KGRAD),XINDX,UDV,PV,ALAB,ISYMA,ISPINA, 2034 & CMO,MJWOP,WRK(KFREE),LFREE) 2035 VAL(1) = DDOT(KZYVC,WRK(KGRAD),1,WRK(KVECC),1) 2036 VAL(2) = DDOT(KZYVC,WRK(KGRAD),1,WRK(KVECC+KZYVC),1) 2037 TMPVAL(1) = TMPVAL(1) + VAL(1) 2038 TMPVAL(2) = TMPVAL(2) + VAL(2) 2039 HYPVAL(1) = HYPVAL(1) + VAL(1) 2040 HYPVAL(2) = HYPVAL(2) + VAL(2) 2041C 2042 IF (DAMPING.GT.D0) THEN 2043 CALL DZERO(WRK(KGRAD),2*MZYVMX) 2044 CALL A2INIT(1,KZYVB,KZYVC,ISYMB,ISPINB,ISYMC,ISPINC,1, 2045 & WRK(KVECC+KZYVC),WRK(KGRAD),XINDX,UDV,PV,ALAB,ISYMA, 2046 & ISPINA,CMO,MJWOP,WRK(KFREE),LFREE) 2047 VAL(1) = -DDOT(KZYVB,WRK(KGRAD),1,WRK(KVECB+KZYVB),1) 2048 VAL(2) = DDOT(KZYVB,WRK(KGRAD),1,WRK(KVECB),1) 2049 TMPVAL(1) = TMPVAL(1) + VAL(1) 2050 TMPVAL(2) = TMPVAL(2) + VAL(2) 2051 HYPVAL(1) = HYPVAL(1) + VAL(1) 2052 HYPVAL(2) = HYPVAL(2) + VAL(2) 2053C 2054 CALL DZERO(WRK(KGRAD),2*MZYVMX) 2055 CALL A2INIT(1,KZYVC,KZYVB,ISYMC,ISPINC,ISYMB,ISPINB,1, 2056 & WRK(KVECB+KZYVB),WRK(KGRAD),XINDX,UDV,PV,ALAB,ISYMA, 2057 & ISPINA,CMO,MJWOP,WRK(KFREE),LFREE) 2058 VAL(1) = -DDOT(KZYVC,WRK(KGRAD),1,WRK(KVECC+KZYVC),1) 2059 VAL(2) = DDOT(KZYVC,WRK(KGRAD),1,WRK(KVECC),1) 2060 TMPVAL(1) = TMPVAL(1) + VAL(1) 2061 TMPVAL(2) = TMPVAL(2) + VAL(2) 2062 HYPVAL(1) = HYPVAL(1) + VAL(1) 2063 HYPVAL(2) = HYPVAL(2) + VAL(2) 2064 END IF 2065C 2066 IF (IPRABS.GE.0) WRITE(LUPRI,'(A15,4F15.6)') 2067 & ' Nx A[2] Ny ',TMPVAL(1),TMPVAL(2),HYPVAL(1),HYPVAL(2) 2068C 2069C Calculate Na*(E[3]-w*S[3]-i*W*R[3])*Nb*Nc type terms (two permutations) 2070C 2071 CALL T3DRV(1,ISYMA,ISYMB,ISYMC,WRK(KVECB),WRK(KVECC), 2072 & .FALSE.,WRK(KVECA),-FREQB,-FREQC,XINDX,UDV,PV,MJWOP, 2073 & WRK(KFREE),LFREE,CMO,FC,FV) 2074 IF (IPRABS.GE.10) THEN 2075 WRITE(LUPRI,'(/A)') ' E[3] vector ' 2076 CALL PRSYMB(LUPRI,'=',72,1) 2077 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 2078 CALL OUTPUT(WRK(KFREE),1,KZYVA/2,1,4,KZYVA/2,4,1,LUPRI) 2079 END IF 2080 VAL(1) = -DDOT(KZYVA,WRK(KVECA),1,WRK(KFREE),1) 2081 VAL(2) = -DDOT(KZYVA,WRK(KVECA+KZYVA),1,WRK(KFREE),1) 2082 VAL(3) = D0 2083 VAL(4) = D0 2084 CALL R3DRV(ISYMA,ISYMB,ISYMC,KZYVA,KZYVB,KZYVC, 2085 & WRK(KVECA),WRK(KVECB),WRK(KVECC),DAMPING, 2086 & XINDX,UDV,MJWOP,WRK(KFREE),LFREE,VAL(3)) 2087 CALL R3DRV(ISYMA,ISYMC,ISYMB,KZYVA,KZYVC,KZYVB, 2088 & WRK(KVECA),WRK(KVECC),WRK(KVECB),DAMPING, 2089 & XINDX,UDV,MJWOP,WRK(KFREE),LFREE,VAL(3)) 2090 TMPVAL(1) = VAL(1) + VAL(3) 2091 TMPVAL(2) = VAL(2) + VAL(4) 2092 HYPVAL(1) = HYPVAL(1) + VAL(1) + VAL(3) 2093 HYPVAL(2) = HYPVAL(2) + VAL(2) + VAL(4) 2094C 2095 IF (DAMPING.GT.D0) THEN 2096 CALL T3DRV(1,ISYMA,ISYMB,ISYMC,WRK(KVECB+KZYVB),WRK(KVECC), 2097 & .FALSE.,WRK(KVECA),-FREQB,-FREQC,XINDX,UDV,PV,MJWOP, 2098 & WRK(KFREE),LFREE,CMO,FC,FV) 2099 VAL(1) = -DDOT(KZYVA,WRK(KVECA),1,WRK(KFREE),1) 2100 VAL(2) = -DDOT(KZYVA,WRK(KVECA+KZYVA),1,WRK(KFREE),1) 2101 VAL(3) = D0 2102 VAL(4) = D0 2103 CALL R3DRV(ISYMA,ISYMB,ISYMC,KZYVA,KZYVB,KZYVC, 2104 & WRK(KVECA),WRK(KVECB+KZYVB),WRK(KVECC),DAMPING, 2105 & XINDX,UDV,MJWOP,WRK(KFREE),LFREE,VAL(3)) 2106 CALL R3DRV(ISYMA,ISYMC,ISYMB,KZYVA,KZYVC,KZYVB, 2107 & WRK(KVECA),WRK(KVECC),WRK(KVECB+KZYVB),DAMPING, 2108 & XINDX,UDV,MJWOP,WRK(KFREE),LFREE,VAL(3)) 2109 TMPVAL(1) = VAL(1) + VAL(3) 2110 TMPVAL(2) = VAL(2) + VAL(4) 2111 TMPVAL(1) = TMPVAL(1) - VAL(2) - VAL(4) 2112 TMPVAL(2) = TMPVAL(2) + VAL(1) + VAL(3) 2113 HYPVAL(1) = HYPVAL(1) - VAL(2) - VAL(4) 2114 HYPVAL(2) = HYPVAL(2) + VAL(1) + VAL(3) 2115C 2116 CALL T3DRV(1,ISYMA,ISYMB,ISYMC,WRK(KVECB),WRK(KVECC+KZYVC), 2117 & .FALSE.,WRK(KVECA),-FREQB,-FREQC,XINDX,UDV,PV,MJWOP, 2118 & WRK(KFREE),LFREE,CMO,FC,FV) 2119 VAL(1) = -DDOT(KZYVA,WRK(KVECA),1,WRK(KFREE),1) 2120 VAL(2) = -DDOT(KZYVA,WRK(KVECA+KZYVA),1,WRK(KFREE),1) 2121 VAL(3) = D0 2122 VAL(4) = D0 2123 CALL R3DRV(ISYMA,ISYMB,ISYMC,KZYVA,KZYVB,KZYVC, 2124 & WRK(KVECA),WRK(KVECB),WRK(KVECC+KZYVC),DAMPING, 2125 & XINDX,UDV,MJWOP,WRK(KFREE),LFREE,VAL(3)) 2126 CALL R3DRV(ISYMA,ISYMC,ISYMB,KZYVA,KZYVC,KZYVB, 2127 & WRK(KVECA),WRK(KVECC+KZYVC),WRK(KVECB),DAMPING, 2128 & XINDX,UDV,MJWOP,WRK(KFREE),LFREE,VAL(3)) 2129 TMPVAL(1) = TMPVAL(1) - VAL(2) - VAL(4) 2130 TMPVAL(2) = TMPVAL(2) + VAL(1) + VAL(3) 2131 HYPVAL(1) = HYPVAL(1) - VAL(2) - VAL(4) 2132 HYPVAL(2) = HYPVAL(2) + VAL(1) + VAL(3) 2133C 2134 CALL T3DRV(1,ISYMA,ISYMB,ISYMC,WRK(KVECB+KZYVB), 2135 & WRK(KVECC+KZYVC), 2136 & .FALSE.,WRK(KVECA),-FREQB,-FREQC,XINDX,UDV,PV,MJWOP, 2137 & WRK(KFREE),LFREE,CMO,FC,FV) 2138 VAL(1) = -DDOT(KZYVA,WRK(KVECA),1,WRK(KFREE),1) 2139 VAL(2) = -DDOT(KZYVA,WRK(KVECA+KZYVA),1,WRK(KFREE),1) 2140 VAL(3) = D0 2141 VAL(4) = D0 2142 CALL R3DRV(ISYMA,ISYMB,ISYMC,KZYVA,KZYVB,KZYVC, 2143 & WRK(KVECA),WRK(KVECB+KZYVB),WRK(KVECC+KZYVC),DAMPING, 2144 & XINDX,UDV,MJWOP,WRK(KFREE),LFREE,VAL(3)) 2145 CALL R3DRV(ISYMA,ISYMC,ISYMB,KZYVA,KZYVC,KZYVB, 2146 & WRK(KVECA),WRK(KVECC+KZYVC),WRK(KVECB+KZYVB),DAMPING, 2147 & XINDX,UDV,MJWOP,WRK(KFREE),LFREE,VAL(3)) 2148 TMPVAL(1) = TMPVAL(1) - VAL(1) - VAL(3) 2149 TMPVAL(2) = TMPVAL(2) - VAL(2) - VAL(4) 2150 HYPVAL(1) = HYPVAL(1) - VAL(1) - VAL(3) 2151 HYPVAL(2) = HYPVAL(2) - VAL(2) - VAL(4) 2152 END IF 2153C 2154 IF (IPRABS.GE.0) WRITE(LUPRI,'(A15,4F15.6)') 2155 & ' Na T[3] NbNc ', 2156 & TMPVAL(1),TMPVAL(2),HYPVAL(1),HYPVAL(2) 2157C 2158C End of conditional for having to actually calculate the response function 2159C 2160 END IF 2161C 2162C Save value for beta, which equals minus the value of the QRF. 2163C 2164 RES_BETA(IQRF,1) = -HYPVAL(1) 2165 RES_BETA(IQRF,2) = -HYPVAL(2) 2166C 2167C For restart purposes, we write the response function to file 2168C 2169 WRITE(LURSPRES,'(I4,3A10,3F10.6,2F16.6)') IQRF, 2170 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3), 2171 & RES_BETA(IQRF,1),RES_BETA(IQRF,2) 2172C 2173C End loop over QRFs 2174C 2175 END DO 2176 CALL GPCLOSE(LURSPRES,'KEEP') 2177C 2178 RETURN 2179 END 2180 SUBROUTINE QRCHK(DOHYP,ALAB,BLAB,CLAB,ISYMA,ISYMB,ISYMC, 2181 & FREQA,FREQB,FREQC) 2182#include "implicit.h" 2183#include "priunit.h" 2184#include "absorp.h" 2185#include "abslrs.h" 2186C 2187 PARAMETER (THRZERO = 1.0D-6) 2188 LOGICAL DOHYP,NEWFRQ 2189 CHARACTER*8 ALAB,BLAB,CLAB,LAB(3) 2190 DIMENSION FREQ(3),ISYM(3) 2191C 2192 DOHYP = .TRUE. 2193C 2194C Operator requirement for magnetic circular dichroism 2195C 2196 IF (ABS_MCD .OR. ABSLRS_MCD) THEN 2197 IF (ALAB(2:8).NE.'DIPLEN' .OR. BLAB(2:8).NE.'DIPLEN' .OR. 2198 & CLAB(2:8).NE.'ANGMOM') DOHYP = .FALSE. 2199 END IF 2200C 2201C Check if equivalent QRF is in the list already. 2202C 2203 IF (DAMPING .EQ. 0.0D0) THEN 2204C Overall permutational symmetry. 2205 JCTR=3 2206 KCTR=1 2207 LCTR=1 2208 ELSE 2209C Only intrinsic permutational symmetry, so A operator has to match 2210C first operator in list of response functions. 2211 JCTR=1 2212 KCTR=2 2213 LCTR=2 2214 END IF 2215 DO IQRF = 1,NQRF 2216 DO J = 1,JCTR 2217 DO K = KCTR,3 2218 IF (K.NE.J) THEN 2219 DO L = LCTR,3 2220 IF (L.NE.K .AND. L.NE.J) THEN 2221C 2222 IF ( ALAB.EQ.QRFLAB(IQRF,J) .AND. 2223 & BLAB.EQ.QRFLAB(IQRF,K) .AND. 2224 & CLAB.EQ.QRFLAB(IQRF,L) .AND. 2225 & ISYMA.EQ.QRFSYM(IQRF,J) .AND. 2226 & ISYMB.EQ.QRFSYM(IQRF,K) .AND. 2227 & ISYMC.EQ.QRFSYM(IQRF,L) .AND. 2228 & ABS( FREQA-QRFFRQ(IQRF,J)).LT.THRZERO .AND. 2229 & ABS( FREQB-QRFFRQ(IQRF,K)).LT.THRZERO .AND. 2230 & ABS( FREQC-QRFFRQ(IQRF,L)).LT.THRZERO ) THEN 2231 DOHYP = .FALSE. 2232 END IF 2233C 2234 END IF 2235 END DO 2236 END IF 2237 END DO 2238 END DO 2239 END DO 2240C 2241 IF (DOHYP) THEN 2242C 2243C Check if this QRF will inflict new LR solver frequencies. 2244C 2245 FREQ(1) = FREQA 2246 FREQ(2) = FREQB 2247 FREQ(3) = FREQC 2248 DO I=1,3 2249 NEWFRQ = .TRUE. 2250 DO IFR=1,NFREQ_ALPHA 2251 IF (ABS(FREQ_ALPHA(IFR)-ABS(FREQ(I))).LT.THRZERO) THEN 2252 NEWFRQ = .FALSE. 2253 END IF 2254 END DO 2255 IF (NEWFRQ) THEN 2256 IF (NFREQ_ALPHA.GE.MXFREQ) THEN 2257 WRITE(LUPRI,'(2(/A),I4,A,/A)') 2258 & ' The specified calculation requires more than the allowed', 2259 & ' number of frequencies in the LR solver MXFREQ=',MXFREQ,'.', 2260 & ' The program will stop.' 2261 CALL QUIT('Too many frequencies in LR solver.') 2262 END IF 2263 NFREQ_ALPHA = NFREQ_ALPHA + 1 2264 FREQ_ALPHA(NFREQ_ALPHA) = ABS(FREQ(I)) 2265 END IF 2266 END DO 2267C 2268 END IF 2269C 2270 IF (NQRF.GE.MXQRF .AND. DOHYP) THEN 2271 WRITE(LUPRI,'(2(/A),I4,A,/A)') 2272 & ' The specified calculation requires more than the allowed', 2273 & ' number of quadratic response functions MXQRF=',MXQRF,'.', 2274 & ' The program will stop' 2275 CALL QUIT('Too many quadratic response functions specified.') 2276 END IF 2277C 2278 RETURN 2279 END 2280 SUBROUTINE QRRDVE(ISYMA,ISYMB,ISYMC,ALAB,BLAB,CLAB, 2281 & FREQA,FREQB,FREQC,KZYVA,KZYVB,KZYVC,VECA,VECB,VECC) 2282C 2283C PURPOSE: Read in response vectors for quadratic response. 2284C 2285#include "implicit.h" 2286#include "priunit.h" 2287#include "absorp.h" 2288#include "inftap.h" 2289#include "rspprp.h" 2290#include "inflr.h" 2291#include "qrinf.h" 2292C 2293 LOGICAL FOUND, CONV 2294 CHARACTER*8 ALAB,BLAB,CLAB,BLANK 2295 PARAMETER ( D0=0.0D0, DM1=-1.0D0 ) 2296 DIMENSION VECA(*),VECB(*),VECC(*) 2297C 2298 KZYVA = MZYVAR(ISYMA) 2299 KZYVB = MZYVAR(ISYMB) 2300 KZYVC = MZYVAR(ISYMC) 2301C 2302 IF (IPRABS.GE.2) THEN 2303 WRITE(LUPRI,'(2(/A),2(/A,3(I10)),/A,3A10,/A,3F10.6)') 2304 & ' Variables in QRRDVE', 2305 & ' ==================================================== ', 2306 & ' KZYVA,KZYVB,KZYVC: ',KZYVA,KZYVB,KZYVC, 2307 & ' ISYMA,ISYMB,ISYMC: ',ISYMA,ISYMB,ISYMC, 2308 & ' ALAB,BLAB,CLAB : ',ALAB,BLAB,CLAB, 2309 & ' FREQA,FREQB,FREQC: ',FREQA,FREQB,FREQC 2310 END IF 2311C 2312C Read in Na 2313C 2314 CALL REARSP(LURSP,2*KZYVA,VECA,ALAB,'COMPLEX ',ABS(FREQA),D0, 2315 & ISYMA,0,THCLR,FOUND,CONV,ANTSYM) 2316 IF (.NOT. (FOUND .AND. CONV)) THEN 2317 IF (.NOT. FOUND) THEN 2318 WRITE (LUPRI,'(/3A,F8.5,A,I3,/A)') ' Response label ',ALAB, 2319 & ' with frequency ',FREQA, ' and symmetry', 2320 & ISYMA,' not found on file RSPVEC' 2321 CALL QUIT('Response vector not found on file RSPVEC') 2322 ELSE 2323 WRITE (LUPRI,'(/3A,F8.5,/A,I3,A)') ' @WARNING----'// 2324 & ' Response label ',ALAB, 2325 & ' with frequency ',FREQA, ' and symmetry', 2326 & ISYMA,' not converged on file RSPVEC' 2327 END IF 2328 END IF 2329 IF (FREQA .GT. D0) THEN 2330 CALL DSWAP(KZYVA/2,VECA,1,VECA(1+KZYVA/2),1) 2331 CALL DSWAP(KZYVA/2,VECA(1+KZYVA),1,VECA(1+KZYVA+KZYVA/2),1) 2332 CALL DSCAL(KZYVA,DM1,VECA,1) 2333 END IF 2334C 2335 IF (IPRABS.GE.10) THEN 2336 WRITE(LUPRI,'(/A)') ' Na vector in ABS_READVEC ' 2337 CALL PRSYMB(LUPRI,'=',72,1) 2338 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 2339 CALL OUTPUT(VECA,1,KZYVA/2,1,4,KZYVA/2,4,1,LUPRI) 2340 END IF 2341C 2342C Read in Nb 2343C 2344 CALL REARSP(LURSP,2*KZYVB,VECB,BLAB,'COMPLEX ',ABS(FREQB),D0, 2345 & ISYMB,0,THCLR,FOUND,CONV,ANTSYM) 2346 IF (.NOT. (FOUND .AND. CONV)) THEN 2347 IF (.NOT. FOUND) THEN 2348 WRITE (LUPRI,'(/3A,F8.5,A,I3,/A)') ' Response label ',BLAB, 2349 & ' with frequency ',FREQB, ' and symmetry', 2350 & ISYMB,' not found on file RSPVEC' 2351 CALL QUIT('Response vector not found on file') 2352 ELSE 2353 WRITE (LUPRI,'(/3A,F8.5,/A,I3,A)') ' @WARNING----'// 2354 & ' Response label ',BLAB, 2355 & ' with frequency ',FREQB, ' and symmetry', 2356 & ISYMB,' not converged on file RSPVEC' 2357 END IF 2358 END IF 2359 IF (FREQB .LT. D0) THEN 2360 CALL DSWAP(KZYVB/2,VECB,1,VECB(1+KZYVB/2),1) 2361 CALL DSWAP(KZYVB/2,VECB(1+KZYVB),1,VECB(1+KZYVB+KZYVB/2),1) 2362 CALL DSCAL(KZYVB,DM1,VECB,1) 2363 END IF 2364C 2365 IF (IPRABS.GE.10) THEN 2366 WRITE(LUPRI,'(/A)') ' Nb vector in ABS_READVEC ' 2367 CALL PRSYMB(LUPRI,'=',72,1) 2368 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 2369 CALL OUTPUT(VECB,1,KZYVB/2,1,4,KZYVB/2,4,1,LUPRI) 2370 END IF 2371C 2372C Read in Nc 2373C 2374 CALL REARSP(LURSP,2*KZYVC,VECC,CLAB,'COMPLEX ',ABS(FREQC),D0, 2375 & ISYMC,0,THCLR,FOUND,CONV,ANTSYM) 2376 IF (.NOT. (FOUND .AND. CONV)) THEN 2377 IF (.NOT. FOUND) THEN 2378 WRITE (LUPRI,'(/3A,F8.5,A,I3,/A)') ' Response label ',CLAB, 2379 & ' with frequency ',FREQC, ' and symmetry', 2380 & ISYMC,' not found on file RSPVEC' 2381 CALL QUIT('Response vector not found on file') 2382 ELSE 2383 WRITE (LUPRI,'(/3A,F8.5,/A,I3,A)') ' @WARNING----'// 2384 & ' Response label ',CLAB, 2385 & ' with frequency ',FREQC, ' and symmetry', 2386 & ISYMC,' not converged on file RSPVEC' 2387 END IF 2388 END IF 2389 IF (FREQC .LT. D0) THEN 2390 CALL DSWAP(KZYVC/2,VECC,1,VECC(1+KZYVC/2),1) 2391 CALL DSWAP(KZYVC/2,VECC(1+KZYVC),1,VECC(1+KZYVC+KZYVC/2),1) 2392 CALL DSCAL(KZYVC,DM1,VECC,1) 2393 END IF 2394C 2395 IF (IPRABS.GE.10) THEN 2396 WRITE(LUPRI,'(/A)') ' Nc vector in ABS_READVEC ' 2397 CALL PRSYMB(LUPRI,'=',72,1) 2398 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 2399 CALL OUTPUT(VECC,1,KZYVC/2,1,4,KZYVC/2,4,1,LUPRI) 2400 END IF 2401C 2402 RETURN 2403 END 2404 SUBROUTINE R3DRV(ISYMA,ISYMB,ISYMC,KZYVA,KZYVB,KZYVC, 2405 & VECA,VECB,VECC,FACTOR, 2406 & XINDX,UDV,MJWOP,WRK,LWRK,RESULT) 2407C 2408C Compute the third-order relaxation contribution times -i*FACTOR (where 2409C FACTOR is equal to the common damping parameter). Results for real and 2410C imaginary parts are added to RESULT(1) and RESULT(2), respectively. 2411C 2412C Na*R[3]*Nb*Nc 2413C 2414#include "implicit.h" 2415#include "priunit.h" 2416#include "infvar.h" 2417#include "wrkrsp.h" 2418#include "infrsp.h" 2419#include "inforb.h" 2420#include "infdim.h" 2421#include "qrinf.h" 2422C 2423 LOGICAL TDM, NORHO2 2424c 2425 PARAMETER ( D0=0.0D0, DH=0.5D0, D1=1.0D0, NTERMS=10 ) 2426 DIMENSION WRK(LWRK),XINDX(*),MJWOP(2,MAXWOP,8),UDV(NASHDI,NASHDI) 2427 DIMENSION VECA(KZYVA),VECB(KZYVB),VECC(KZYVC),RESULT(2) 2428 DIMENSION TMP(2,NTERMS) 2429C 2430 KFREE=1 2431 LFREE=LWRK 2432C 2433C There is no third-order relaxation contribution when the 2434C damping parameter is equal to zero. There is also no contribution 2435C in Hartree-Fock since the average value of triple exciteations must 2436C vanish 2437C 2438 IF (FACTOR.EQ.D0 .OR. TDHF) RETURN 2439C 2440 CALL DZERO(TMP,NTERMS) 2441 CALL MEMGET('REAL',KZYMAT,NORBT*NORBT,WRK,KFREE,LFREE) 2442 CALL MEMGET('REAL',KDEN1,NASHT*NASHT,WRK,KFREE,LFREE) 2443 CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE) 2444C 2445 CALL GETREF(WRK(KCREF),MZCONF(1)) 2446C 2447 IPRTMP=IPRRSP 2448 IPRRSP=201 2449 IPRONE = 0 2450 CALL DSWAP(KZYVA/2,VECA(1),1,VECA(1+KZYVA/2),1) 2451 CALL DSWAP(KZYVA/2,VECA(1+KZYVA),1,VECA(1+KZYVA+KZYVA/2),1) 2452C 2453C 1/2*<0|[nC,[nB,nA]]|0> 2454C 2455 IF (MZWOPT(ISYMA).GT.0 .AND. MZWOPT(ISYMB).GT.0 .AND. 2456 & MZWOPT(ISYMC).GT.0) THEN 2457 CALL TRZYM2(VECA,VECB,VECC,KZYVA,KZYVB,KZYVC, 2458 & ISYMA,ISYMB,ISYMC,WRK(KZYMAT),MJWOP,WRK(KFREE),LFREE) 2459 CALL MELONE(WRK(KZYMAT),1,UDV,D1,TMP(1,1),IPRONE,'R3DRV') 2460 TMP(1,1) = TMP(1,1)*DH 2461 CALL TRZYM2(VECA(1+KZYVA),VECB,VECC,KZYVA,KZYVB,KZYVC, 2462 & ISYMA,ISYMB,ISYMC,WRK(KZYMAT),MJWOP,WRK(KFREE),LFREE) 2463 CALL MELONE(WRK(KZYMAT),1,UDV,D1,TMP(2,1),IPRONE,'R3DRV') 2464 TMP(2,1) = TMP(2,1)*DH 2465 END IF 2466C 2467C <0|[nA,nC] - [nA,nC](+)|0B>, where (+) is the dagger 2468C 2469 IF (MZWOPT(ISYMA).GT.0 .AND. MZCONF(ISYMB).GT.0 .AND. 2470 & MZWOPT(ISYMC).GT.0) THEN 2471 ILSYM = IREFSY 2472 IRSYM = MULD2H(IREFSY,ISYMB) 2473 NCL = MZCONF(1) 2474 NCR = MZCONF(ISYMB) 2475 TDM=.TRUE. 2476 NORHO2=.TRUE. 2477 CALL DZERO(WRK(KDEN1),NASHT*NASHT) 2478 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,WRK(KCREF),VECB, 2479 & WRK(KDEN1),DUMMY,0,0,TDM,NORHO2,XINDX,WRK(KFREE), 2480 & 1,LFREE) 2481 OVLAP = D0 2482 IF (ILSYM.EQ.IRSYM) 2483 & OVLAP = DDOT(NCL,WRK(KCREF),1,VECB,1) 2484C 2485 CALL MEMGET('REAL',KTMP,NORBT*NORBT,WRK,KFREE,LFREE) 2486 CALL TRZYMT(1,VECA,VECC,KZYVA,KZYVC,ISYMA,ISYMC,WRK(KTMP), 2487 & MJWOP,WRK,LWRK) 2488 CALL DCOPY(NORBT*NORBT,WRK(KTMP),1,WRK(KZYMAT),1) 2489 CALL DGEMM('T','N',NORBT,NORBT,NORBT,-1.0D0, 2490 & WRK(KTMP),NORBT,1.0D0,1,1.0D0,WRK(KZYMAT),NORBT) 2491 CALL MELONE(ZYMAT,ISYMV1,DEN1,OVLAP,TMP(1,2),IPRONE,'R3DRV') 2492C 2493 CALL TRZYMT(1,VECA(1+KZYVA),VECC,KZYVA,KZYVC,ISYMA,ISYMC, 2494 & WRK(KTMP),MJWOP,WRK,LWRK) 2495 CALL DCOPY(NORBT*NORBT,WRK(KTMP),1,WRK(KZYMAT),1) 2496 CALL DGEMM('T','N',NORBT,NORBT,NORBT,-1.0D0, 2497 & WRK(KTMP),NORBT,1.0D0,1,1.0D0,WRK(KZYMAT),NORBT) 2498 CALL MELONE(ZYMAT,ISYMV1,WRK(KDEN1),OVLAP,TMP(2,2), 2499 & IPRONE,'R3DRV') 2500 END IF 2501C 2502C 1/2*<0|nCnB + nC(+)nB(+)|0A>, where (+) is the dagger 2503C 2504 IF (MZCONF(ISYMA).GT.0 .AND. MZWOPT(ISYMB).GT.0 .AND. 2505 & MZWOPT(ISYMC).GT.0) THEN 2506 CALL MEMGET('REAL',KDEN2,N2ASHX*N2ASHX,WRK,KFREE,LFREE) 2507 ILSYM = IREFSY 2508 IRSYM = MULD2H(IREFSY,ISYMA) 2509 NCL = MZCONF(1) 2510 NCR = MZCONF(ISYMA) 2511 TDM=.TRUE. 2512 NORHO2=.FALSE. 2513 CALL DZERO(WRK(KDEN1),NASHT*NASHT) 2514 CALL DZERO(WRK(KDEN2),N2ASHX*N2ASHX) 2515 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,WRK(KCREF),VECA, 2516 & WRK(KDEN1),WRK(KDEN2),0,0,TDM,NORHO2,XINDX, 2517 & WRK(KFREE),1,LFREE) 2518 OVLAP = D0 2519 IF (ILSYM.EQ.IRSYM) 2520 & OVLAP = DDOT(NCL,WRK(KCREF),1,VECA,1) 2521 CALL TWOMOM(ISYMB,ISYMC,KZYVB,KZYVC,VECB,VECC,MJWOP, 2522 & WRK(KDEN1),WRK(KDEN2),WRK(KFREE),LFREE, 2523 & TMP(1,3)) 2524 END IF 2525C 2526 IPRRSP=IPRTMP 2527C 2528C Multiply by -i*W (minus i times the gamma factor). 2529C 2530 RESULT(1) = FACTOR*TMP(1,2) 2531 RESULT(2) = -FACTOR*TMP(1,1) 2532C 2533 RETURN 2534 END 2535 SUBROUTINE TWOMOM(ISYMB,ISYMC,KZYVB,KZYVC,VECB,VECC,MJWOP, 2536 & DEN1,DEN2,WRK,LWRK,RESULT) 2537C 2538C Calculate <0|nCnB + nC(+)nB(+)|0A> given the one- and two-electron 2539C transition densities. 2540C 2541#include "implicit.h" 2542#include "priunit.h" 2543#include "infvar.h" 2544#include "wrkrsp.h" 2545#include "infrsp.h" 2546#include "inforb.h" 2547#include "infdim.h" 2548#include "qrinf.h" 2549C 2550 DIMENSION WRK(LWRK),MJWOP(2,MAXWOP,8) 2551 DIMENSION DEN1(NASHDI,NASHDI),DEN2(N2ASHX,N2ASHX) 2552 DIMENSION VECB(KZYVB),VECC(KZYVC),RESULT(2) 2553C 2554 KFREE=1 2555 LFREE=LWRK 2556 CALL QUIT('TWOMOM: to be implemented....') 2557C 2558 RETURN 2559 END 2560 SUBROUTINE ABCCHK(DOCAL,LURSPRES,ISYMA,ISYMB,ISYMC,FREQA,FREQB, 2561 & FREQC,ALAB,BLAB,CLAB,HYPVAL) 2562C 2563#include "implicit.h" 2564#include "iratdef.h" 2565#include "absorp.h" 2566C 2567C PURPOSE: 2568C Check if this quadratic response calculation with absorption may 2569C already exist. To allow for flexible restart of long MCD jobs 2570C K.Ruud, August 2011. Based on BCCHK 2571C 2572 PARAMETER (THD = 1.0D-8) 2573C 2574 LOGICAL DOCAL 2575 CHARACTER*8 LABEL(3), ALAB, BLAB, CLAB 2576 DIMENSION FREQ(3),VALUE(2), HYPVAL(2) 2577 SAVE N 2578 DATA N/0/ 2579#include "priunit.h" 2580C 2581 DOCAL = .TRUE. 2582 REWIND (LURSPRES) 2583 346 READ (LURSPRES,'(I4,3A10,3F10.6,2F16.6)',END=347) 2584 & IQRF, (LABEL(I),I=1,3), (FREQ(I),I=1,3),VALUE(1), 2585 & VALUE(2) 2586 IF ((LABEL(1) .EQ. ALAB) .AND. (LABEL(2) .EQ. BLAB) .AND. 2587 & (LABEL(3) .EQ. CLAB) .AND. 2588 & ABS(FREQB-FREQ(2)).LT.THD 2589 & .AND. ABS(FREQC-FREQ(3)).LT.THD) THEN 2590 WRITE(LUPRI,'(/A)') ' The following quadratic '// 2591 & 'response function has already been calculated' 2592 WRITE (LUPRI,'(3A10,3F10.6)') 2593 & (LABEL(I),I=1,3), (FREQ(I),I=1,3) 2594 DOCAL = .FALSE. 2595 HYPVAL(1) = -VALUE(1) 2596 HYPVAL(2) = -VALUE(2) 2597 END IF 2598 CALL FLSHFO(LUPRI) 2599 GOTO 346 2600 347 CONTINUE 2601#if defined (VAR_MFDS) 2602 BACKSPACE (LURSPRES) 2603#endif 2604 IF (DOCAL) THEN 2605 N = N + 1 2606 IF (N .GT. MXQRF) THEN 2607 WRITE (LUPRI,'(/A,I3)') 2608 & 'ABCCHK: # unique calculations .gt. MXCALC =',MXQRF 2609 CALL QUIT('ABCCHK: # unique calculations .gt. MXQRF') 2610 END IF 2611 END IF 2612C 2613 RETURN 2614 END 2615