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 newabs.F: Main author Joanna Kauczor 20C This implementation is published in: 21C 22C Reference C: 23C J. Kauczor, P. Jorgensen, and P. Norman, 24C "On the efficiency of algorithms for solving Hartree-Fock and Kohn-Sham 25C response equations", 26C J. Chem. Theory Comput. 7 (2011) 1610. 27C ============================================================================ 28C 29 SUBROUTINE ABSORP_INPUT(WORD) 30C 31C Purpose: 32C Read in user settings for imaginary polarizabilities describing 33C absorption in the optical processes. 34C 35C Used from 36C codata.h : XTKAYS 37C gnrinf.h : THR_REDFAC 38C symmet.h : ISYMAX 39C infvar.h : NCONF 40C inftab.h : LUPROP 41#include "implicit.h" 42#include "priunit.h" 43#include "codata.h" 44#include "gnrinf.h" 45#include "maxorb.h" 46C MXCORB for symmet 47#include "maxaqn.h" 48C MXQN for symmet 49#include "mxcent.h" 50C MXCENT for symmet 51#include "symmet.h" 52C ISYMAX 53#include "infvar.h" 54C NCONF 55#include "abslrs.h" 56#include "absorp.h" 57#include "abscrs.h" 58#include "inforb.h" 59#include "inftap.h" 60#include "infinp.h" 61#include "infrsp.h" 62C 63 LOGICAL NEWDEF, ALLCOMP, NXTLAB 64 PARAMETER ( NTABLE = 40 ) 65 PARAMETER ( THRFRQ=1.0D-14, THD=1.0D-6 ) 66 CHARACTER PROMPT*1, WORD*7, WORD1*7, TABLE(NTABLE)*7 67 CHARACTER*8 DIPLEN(3),ANGMOM(3),THETA(6),RTNLBL(2),LABFND 68 CHARACTER*8 DIPVEL(3) 69 INTEGER NUMBER_ORBS(8), NUMBER_EXORBS(8) 70C 71 logical, external :: fun_is_ready_for_qr 72c 73 DATA DIPLEN/'XDIPLEN','YDIPLEN','ZDIPLEN'/ 74 DATA DIPVEL/'XDIPVEL','YDIPVEL','ZDIPVEL'/ 75 DATA ANGMOM/'XANGMOM','YANGMOM','ZANGMOM'/ 76 DATA THETA/'XXTHETA','XYTHETA','XZTHETA','YYTHETA','YZTHETA', 77 & 'ZZTHETA'/ 78 DATA TABLE /'.ALPHA ','.FREQUE','.THCLR ','.MAXIT ','.MAXRM ', 79 & '.PRINT ','.DAMPIN','.ANALYZ','.FREQ I','.MCD ', 80 & '.MCHD ','.NSCD ','.BETA ','.BFREQ ','.BFREQI', 81 & '.CFREQ ','.CFREQI','.SHG ','.XXCOMP','.YYCOMP', 82 & '.ZZCOMP','.IMAG F','.NBATCH','.NOREBD','.OLDCPP', 83 & '.MAX MI','.MAXITO','.EXCITA','.MAX MA','.THCPP ', 84 & '.REDUCE','.BATCH ','.IDRI ','.GAMMA ','.DFREQ ', 85 & '.DFREQI','.ELEMNT','.OTOFGA','.ALLELE','.VELOCI'/ 86C 87 NEWDEF = (WORD .EQ. '*ABSORP') 88 ICHANG = 0 89C 90C Set default values 91 ALLCOMP = .TRUE. !All components of requested tensor 92 ICOMP = 0 !ICHOMP=1,2,3 --> XX,YY,ZZCOMP 93 ABSLRS = .TRUE. !New linear response solver 94 ABSLRS_MCD = .FALSE. !Magnetic circular dichroism 95 ABSLRS_MCHD = .FALSE. !Magneto-chiral dichroism 96 ABSLRS_NSCD = .FALSE. !Nuclear-spin circular dichroism 97 ABSLRS_IDRI = .FALSE. !Intensity dependent refractive index 98 ABS_GAMMA_ELEMENT = .FALSE. !Specific elements for GAMMA 99 ABS_GAMMA_ALLELE = .FALSE. !All elements for GAMMA 100 ABS_OTOF_GAMMA = .FALSE. !One-to-one frequencies for GAMMA 101C 102 IF (NEWDEF) THEN 103C Common default value of the damping parameter is set to 104C be 1000 cm-1 = 4.556333D-3 a.u. 105 ABS_DAMP = 1000/XTKAYS 106 100 CONTINUE 107C 108C Read in input 109C 110 READ (LUCMD, '(A7)') WORD 111 CALL UPCASE(WORD) 112 PROMPT = WORD(1:1) 113 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') GOTO 100 114 IF (PROMPT .EQ. '.') THEN 115 ICHANG = ICHANG + 1 116 DO I = 1, NTABLE 117 IF (TABLE(I) .EQ. WORD) 118 & GOTO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15, 119 & 16,17,18,19,20,21,22,23,24,25,26,27,28,29,30, 120 & 31,32,33,34,35,36,37,38,39,40),I 121 END DO 122 IF (WORD .EQ. '.OPTION') THEN 123 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 124 GOTO 100 125 END IF 126 WRITE (LUPRI,'(/,3A,/)') ' KEYWORD "',WORD, 127 * '" NOT RECOGNIZED IN ABSORP_INPUT.' 128 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 129 CALL QUIT(' ILLEGAL KEYWORD IN ABSORP_INPUT ') 130C 131C .ALPHA 132 1 CONTINUE 133 ABSLRS_ALPHA=.TRUE. 134 ABSLRS_BETA=.FALSE. 135 ABSLRS_GAMMA=.FALSE. 136 GOTO 100 137C 138C .FREQUE 139 2 CONTINUE 140 IF (ABSLRS_ALPHA .OR. ABS_VELOCI) THEN 141 READ (LUCMD,*) ABS_NFREQ_ALPHA 142 IF (ABS_NFREQ_ALPHA.LE.NMXFREQ) THEN 143 READ (LUCMD,*) (ABS_FREQ_ALPHA(J), 144 & J=1,ABS_NFREQ_ALPHA) 145 ELSE 146 WRITE (LUPRI,'(3(/,A,I5),/)') 147 & ' NUMBER OF FREQUENCIES SPECIFIED :', 148 & ABS_NFREQ_ALPHA, 149 & ' IS GREATER THAN THE ALLOWED NUMBER:',NMXFREQ, 150 & ' THE NUMBER IS RESET TO THE MAXIMUM:',NMXFREQ 151 READ (LUCMD,*) (ABS_FREQ_ALPHA(J),J=1,NMXFREQ), 152 & (FFFF,J=NMXFREQ+1,ABS_NFREQ_ALPHA) 153 ABS_NFREQ_ALPHA = NMXFREQ 154 END IF 155 CALL GPDSRT(ABS_NFREQ_ALPHA,ABS_FREQ_ALPHA,THD) 156 ELSE 157 WRITE(LUPRI,1000) WORD 158 READ(LUCMD,*) 159 READ(LUCMD,*) 160 END IF 161 GOTO 100 162C 163C .THCLR 164 3 CONTINUE 165 READ (LUCMD,*) ABS_THCLR 166 IF (ABS_OLDCPP) THCLR_ABSORP=ABS_THCLR 167 GOTO 100 168C 169C .MAXIT 170 4 CONTINUE 171 READ (LUCMD,*) ABS_MAXITER 172 GOTO 100 173C 174C .MAXRM 175 5 CONTINUE 176 READ (LUCMD,*) ABS_MAXRM 177 GOTO 100 178C 179C .PRINT 180 6 CONTINUE 181 READ (LUCMD,*) IPRABSLRS 182 GOTO 100 183C 184C .DAMPIN 185 7 CONTINUE 186 READ (LUCMD,*) ABS_DAMP 187 GOTO 100 188C 189C .ANALYZ 190 8 CONTINUE 191 ABSLRS_ANALYZE=.TRUE. 192 GOTO 100 193C 194C .FREQ I 195 9 CONTINUE 196 ABSLRS_INTERVAL = .TRUE. 197c ABS_REDUCE = .TRUE. 198 READ(LUCMD,*) (ABS_FREQ_INTERVAL(I), I=1,3) 199 GOTO 100 200C 201C .MCD 202 10 CONTINUE 203 ABSLRS_ALPHA=.FALSE. 204 ABSLRS_MCD =.TRUE. 205 GOTO 100 206C 207C .MCHD 208 11 CONTINUE 209 ABSLRS_ALPHA=.FALSE. 210 ABSLRS_MCHD =.TRUE. 211 GOTO 100 212C 213C .NSCD 214 12 CONTINUE 215 ABSLRS_ALPHA=.FALSE. 216 ABSLRS_NSCD =.TRUE. 217 GOTO 100 218C 219C .BETA 220 13 CONTINUE 221 ABSLRS_ALPHA=.FALSE. 222 ABSLRS_BETA =.TRUE. 223 ABSLRS_GAMMA =.FALSE. 224 GOTO 100 225C 226C .BFREQ 227 14 CONTINUE 228 IF (ABSLRS_BETA .OR. ABSLRS_MCD .OR. ABSLRS_MCHD .OR. 229 & ABSLRS_NSCD) THEN 230 READ (LUCMD,*) ABS_NFREQ_BETA_B 231 IF (ABS_NFREQ_BETA_B.LE.NMXFREQ) THEN 232 READ (LUCMD,*) (ABS_FREQ_BETA_B(J), 233 & J=1,ABS_NFREQ_BETA_B) 234 ELSE 235 WRITE (LUPRI,'(3(/,A,I5),/)') 236 & ' NUMBER OF FREQUENCIES SPECIFIED :', 237 & ABS_NFREQ_BETA_B, 238 & ' IS GREATER THAN THE ALLOWED NUMBER:',NMXFREQ, 239 & ' THE NUMBER IS RESET TO THE MAXIMUM:',NMXFREQ 240 READ (LUCMD,*) (ABS_FREQ_BETA_B(J),J=1,NMXFREQ), 241 & (FFFF,J=NMXFREQ+1,ABS_NFREQ_BETA_B) 242 ABS_NFREQ_BETA_B = NMXFREQ 243 END IF 244 ELSE IF (ABSLRS_GAMMA .OR. ABSLRS_IDRI) THEN 245 READ (LUCMD,*) ABS_NFREQ_GAMMA_B 246 IF (ABS_NFREQ_GAMMA_B .LE. NMXFREQ) THEN 247 READ (LUCMD,*) (ABS_FREQ_GAMMA_B(J), 248 & J = 1,ABS_NFREQ_GAMMA_B) 249 IF (ABSLRS_IDRI) THEN 250 ABS_NFREQ_GAMMA_C = ABS_NFREQ_GAMMA_B 251 ABS_NFREQ_GAMMA_D = ABS_NFREQ_GAMMA_B 252 DO J= 1,ABS_NFREQ_GAMMA_B 253 ABS_FREQ_GAMMA_C(J) = -ABS_FREQ_GAMMA_B(J) 254 ABS_FREQ_GAMMA_D(J) = ABS_FREQ_GAMMA_B(J) 255 END DO 256 END IF 257 ELSE 258 WRITE (LUPRI,'(3(/,A,I5),/)') 259 & ' NUMBER OF FREQUENCIES SPECIFIED :', 260 & ABS_NFREQ_GAMMA_B, 261 & ' IS GREATER THAN THE ALLOWED NUMBER:',NMXFREQ, 262 & ' THE NUMBER IS RESET TO THE MAXIMUM:',NMXFREQ 263 READ (LUCMD,*) (ABS_FREQ_GAMMA_B(J),J=1,NMXFREQ), 264 & (FFFF,J=NMXFREQ+1,ABS_NFREQ_GAMMA_B) 265 ABS_NFREQ_GAMMA_B = NMXFREQ 266 IF (ABSLRS_IDRI) THEN 267 DO J=1,ABS_NFREQ_GAMMA_B 268 ABS_FREQ_GAMMA_C(J)= -ABS_FREQ_GAMMA_B(J) 269 ABS_FREQ_GAMMA_D(J)= ABS_FREQ_GAMMA_B(J) 270 END DO 271 END IF 272 END IF 273 ELSE 274 WRITE(LUPRI,1010) WORD 275 READ(LUCMD,*) 276 READ(LUCMD,*) 277 END IF 278 GOTO 100 279C 280C .BFREQI 281 15 CONTINUE 282 ABSQRS_INTERVAL = .TRUE. 283c ABS_REDUCE = .TRUE. 284 READ(LUCMD,*) (ABS_BFREQ_INTERVAL(I), I=1,3) 285 GOTO 100 286C 287C .CFREQ 288 16 CONTINUE 289 IF (ABSLRS_BETA) THEN 290 READ (LUCMD,*) ABS_NFREQ_BETA_C 291 IF (ABS_NFREQ_BETA_C.LE.NMXFREQ) THEN 292 READ (LUCMD,*) (ABS_FREQ_BETA_C(J), 293 & J=1,ABS_NFREQ_BETA_C) 294 ELSE 295 WRITE (LUPRI,'(3(/,A,I5),/)') 296 & ' NUMBER OF FREQUENCIES SPECIFIED :',ABS_NFREQ_BETA_C, 297 & ' IS GREATER THAN THE ALLOWED NUMBER:',NMXFREQ, 298 & ' THE NUMBER IS RESET TO THE MAXIMUM:',NMXFREQ 299 READ (LUCMD,*) (ABS_FREQ_BETA_C(J),J=1,NMXFREQ), 300 & (FFFF,J=NMXFREQ+1,ABS_NFREQ_BETA_C) 301 ABS_NFREQ_BETA_C = NMXFREQ 302 END IF 303 ELSE IF (ABSLRS_GAMMA) THEN 304 READ (LUCMD,*) ABS_NFREQ_GAMMA_C 305 IF (ABS_NFREQ_GAMMA_C .LE. NMXFREQ) THEN 306 READ (LUCMD,*) (ABS_FREQ_GAMMA_C(J), 307 & J = 1,ABS_NFREQ_GAMMA_C) 308 ELSE 309 WRITE (LUPRI,'(3(/,A,I5),/)') 310 & ' NUMBER OF FREQUENCIES SPECIFIED :', 311 & ABS_NFREQ_GAMMA_C, 312 & ' IS GREATER THAN THE ALLOWED NUMBER:',NMXFREQ, 313 & ' THE NUMBER IS RESET TO THE MAXIMUM:',NMXFREQ 314 READ (LUCMD,*) (ABS_FREQ_GAMMA_C(J),J=1,NMXFREQ), 315 & (FFFF,J=NMXFREQ+1,ABS_NFREQ_GAMMA_C) 316 ABS_NFREQ_GAMMA_C = NMXFREQ 317 END IF 318 ELSE 319 WRITE(LUPRI,1010) WORD 320 READ(LUCMD,*) 321 READ(LUCMD,*) 322 END IF 323 GOTO 100 324C 325C .CFREQI 326 17 CONTINUE 327 ABSQRS_CINTERVAL = .TRUE. 328c ABS_REDUCE = .TRUE. 329 READ(LUCMD,*) (ABS_CFREQ_INTERVAL(I), I=1,3) 330 GOTO 100 331C 332C .SHG 333 18 CONTINUE 334 ABSLRS_ALPHA=.FALSE. 335 ABSLRS_BETA =.TRUE. 336 ABSLRS_GAMMA=.FALSE. 337 ABSLRS_SHG =.TRUE. 338 GOTO 100 339C 340C .XXCOMP 341 19 CONTINUE 342 ALLCOMP=.FALSE. 343 ICOMP=1 344 GOTO 100 345C 346C .YYCOMP 347 20 CONTINUE 348 ALLCOMP=.FALSE. 349 ICOMP=2 350 GOTO 100 351C 352C .ZZCOMP 353 21 CONTINUE 354 ALLCOMP=.FALSE. 355 ICOMP=3 356 GOTO 100 357C 358C .IMAG F 359 22 CONTINUE 360 ABS_IMFREQ=.TRUE. 361 ABS_DAMP=0.0d0 362 GOTO 100 363C 364C .NBATCH 365 23 CONTINUE 366 ABS_BATCH=.TRUE. 367 READ(LUCMD,*) ABS_NBATCHMAX 368 GOTO 100 369C 370C .NOREBD 371 24 CONTINUE 372 ABS_REBD = .FALSE. 373 GOTO 100 374C 375C .OLDCPP 376 25 CONTINUE 377 ABS_OLDCPP = .TRUE. 378 GOTO 100 379C 380C .MAX MI 381 26 CONTINUE 382 IF (ABS_OLDCPP) THEN 383 READ (LUCMD,*) MAX_MICRO 384 ELSE 385 WRITE(LUPRI,1020) WORD 386 READ(LUCMD,*) 387 END IF 388 GOTO 100 389C 390C .MAXITO 391 27 CONTINUE 392 IF (ABS_OLDCPP) THEN 393 READ (LUCMD,*) MAX_ITORB 394 IF (MAX_ITORB .GT. 0) OPTORB = .TRUE. 395 ELSE 396 WRITE(LUPRI,1020) WORD 397 READ(LUCMD,*) 398 END IF 399 GOTO 100 400C 401C .EXCITA 402 28 CONTINUE 403 IF (ABS_OLDCPP) THEN 404 READ (LUCMD,*) NEXCITED_STATES 405 ELSE 406 WRITE(LUPRI,1020) WORD 407 READ(LUCMD,*) 408 END IF 409 GOTO 100 410C 411C .MAX MA 412 29 CONTINUE 413 IF (ABS_OLDCPP) THEN 414 READ (LUCMD,*) MAX_MACRO 415 ELSE 416 WRITE(LUPRI,1020) WORD 417 READ(LUCMD,*) 418 END IF 419 GOTO 100 420C 421C .THCPP 422 30 CONTINUE 423 IF (ABS_OLDCPP) THEN 424 READ (LUCMD,*) THCPP_ABSORP 425 ABS_THCLR=THCPP_ABSORP 426 ELSE 427 WRITE(LUPRI,1020) WORD 428 READ(LUCMD,*) 429 END IF 430 GOTO 100 431C 432C .REDUCE 433 31 CONTINUE 434 IF (ABS_OLDCPP) THEN 435 ABS_REDUCE=.TRUE. 436 ELSE 437 WRITE(LUPRI,1020) WORD 438 END IF 439 GOTO 100 440C 441C .BATCH 442 32 CONTINUE 443 READ(LUCMD,*) NFREQ_BATCH 444 IF (NFREQ_BATCH.GT.MXFREQ) THEN 445 NFREQ_BATCH = MXFREQ 446 END IF 447 GOTO 100 448C 449C .IDRI 450 33 CONTINUE 451 ABSLRS_ALPHA = .FALSE. 452 ABSLRS_BETA = .FALSE. 453 ABSLRS_GAMMA = .FALSE. 454 ABSLRS_IDRI = .TRUE. 455 GOTO 100 456C 457C .GAMMA 458 34 CONTINUE 459 ABSLRS_ALPHA = .FALSE. 460 ABSLRS_BETA = .FALSE. 461 ABSLRS_GAMMA = .TRUE. 462 GOTO 100 463C 464C .DFREQ 465 35 CONTINUE 466 IF (ABSLRS_GAMMA) THEN 467 READ (LUCMD,*) ABS_NFREQ_GAMMA_D 468 IF (ABS_NFREQ_GAMMA_D .LE. NMXFREQ) THEN 469 READ (LUCMD,*) (ABS_FREQ_GAMMA_D(J), 470 & J=1,ABS_NFREQ_GAMMA_D) 471 ELSE 472 WRITE (LUPRI,'(3(/,A,I5),/)') 473 & ' NUMBER OF FREQUENCIES SPECIFIED :',ABS_NFREQ_GAMMA_D, 474 & ' IS GREATER THAN THE ALLOWED NUMBER:',NMXFREQ, 475 & ' THE NUMBER IS RESET TO THE MAXIMUM:',NMXFREQ 476 READ (LUCMD,*) (ABS_FREQ_GAMMA_D(J), 477 & J=1,NMXFREQ), 478 & (FFFF,J=NMXFREQ+1,ABS_NFREQ_GAMMA_D) 479 ABS_NFREQ_GAMMA_D = NMXFREQ 480 END IF 481 ELSE 482 WRITE(LUPRI,1010) WORD 483 READ(LUCMD,*) 484 READ(LUCMD,*) 485 END IF 486 GOTO 100 487C 488C .DFREQI 489 36 CONTINUE 490 GOTO 100 491C 492C .ELEMNT 493 37 CONTINUE 494 IF (ABSLRS_GAMMA .AND. .NOT. ABS_GAMMA_ALLELE) THEN 495 ABS_GAMMA_ELEMENT = .TRUE. 496 READ (LUCMD,*) ABS_NELEMENTS_GAMMA 497 IF (ABS_NELEMENTS_GAMMA .GT. 81) THEN 498 CALL QUIT(' Too many elements specified for'// 499 & ' gamma tensor: ',ABS_NELEMENTS_GAMMA, 500 & ', maximum = 81) ') 501 END IF 502 READ (LUCMD,*) 503 & (ABS_ELEMENTS_GAMMA(J),J=1,ABS_NELEMENTS_GAMMA) 504 END IF 505 GOTO 100 506C 507C .OTOFGA 508 38 CONTINUE 509 IF (ABSLRS_GAMMA .OR. ABSLRS_IDRI) THEN 510 ABS_OTOF_GAMMA = .TRUE. 511 END IF 512 GOTO 100 513C 514C .ALLELE 515 39 CONTINUE 516 IF (ABSLRS_GAMMA .AND. .NOT. ABS_GAMMA_ELEMENT) THEN 517 ABS_GAMMA_ALLELE = .TRUE. 518 ENDIF 519 GOTO 100 520 40 CONTINUE 521 ABS_VELOCI=.TRUE. 522 GOTO 100 523C 524 ELSE IF (PROMPT .EQ. '*') THEN 525 GOTO 300 526 ELSE 527 WRITE (LUPRI,'(/,3A,/)') ' PROMPT "',WORD, 528 * '" NOT RECOGNIZED IN ABSORPTION INPUT.' 529 CALL QUIT(' ILLEGAL PROMPT IN ABSORPTION INPUT ') 530 END IF 531 GOTO 100 532 END IF 533 300 CONTINUE 534 IF (THR_REDFAC .GT. 0.0D0) THEN 535 ICHANG = ICHANG + 1 536 WRITE (LUPRI,'(3A,1P,D10.2)') '@ INFO ',WORD1, 537 & ' thresholds multiplied with general factor',THR_REDFAC 538 ABS_THCLR = ABS_THCLR*THR_REDFAC 539 THCPP_ABSORP=ABS_THCLR 540 END IF 541C 542 ABS_MAXRM=MAX(ABS_MAXRM,MAXRM) 543C 544C Process user input 545C 546#ifndef DISABLE_XC_RESPONSE_SANITY_CHECK 547! sanity check for incomplete/untested xc functional implementations 548 IF ((ABSLRS_BETA .OR. ABSLRS_GAMMA .OR. ABSLRS_MCD .OR. 549 & ABSLRS_MCHD .OR. ABSLRS_NSCD) .AND. DODFT) THEN 550 IF (.NOT. fun_is_ready_for_qr()) THEN 551 WRITE(LUPRI, *) 'ERROR: functional not fully implemented' 552 WRITE(LUPRI, *) ' or tested for QR' 553 WRITE(LUPRI, *) 'to disable this stop recompile' 554 WRITE(LUPRI, *) 'with -DDISABLE_XC_RESPONSE_SANITY_CHECK' 555 WRITE(LUPRI, *) 'note that GGAKEY functionals are always' 556 WRITE(LUPRI, *) 'stopped although they could be correct' 557 CALL QUIT('functional not fully implemented/tested for QR') 558 END IF 559 END IF 560#endif 561 562C 563 ABSLRS = ABSLRS_ALPHA .OR. ABSLRS_BETA .OR. ABSLRS_GAMMA .OR. 564 & ABSLRS_MCD .OR. ABSLRS_MCHD .OR. ABSLRS_NSCD .OR. 565 & ABSLRS_GAMMA .OR. ABSLRS_IDRI .OR. ABS_VELOCI 566 IF (.NOT.ABSLRS) THEN 567 WRITE(LUPRI,'(/A)') ' Absorption input ignored because:' 568 WRITE(LUPRI,'(A,L2)')' No process requested:', 569 & ' ABS_ALPHA = ', ABSLRS 570 ELSE 571C 572C Put operators in lists 573C 574 CALL IZERO(ABS_NOPER,8) 575 CALL IZERO(NOPER,8) 576 DO ILABEL=1,3 577 IF (ALLCOMP .OR. ICOMP.EQ.ILABEL) THEN 578 ISYM = ISYMAX(ILABEL,1)+1 579 ABS_NOPER(ISYM) = ABS_NOPER(ISYM) + 1 580 NOPER(ISYM)=NOPER(ISYM)+1 581 IF (.NOT. ABS_VELOCI) THEN 582 ABS_LABOP(ABS_NOPER(ISYM),ISYM) = DIPLEN(ILABEL) 583 LABOP(NOPER(ISYM),ISYM) = DIPLEN(ILABEL) 584 ELSE 585 ABS_LABOP(ABS_NOPER(ISYM),ISYM) = DIPVEL(ILABEL) 586 LABOP(NOPER(ISYM),ISYM) = DIPVEL(ILABEL) 587 END IF 588 IF (ABSLRS_MCD) THEN 589 ISYM = ISYMAX(ILABEL,2)+1 590 ABS_NOPER(ISYM) = ABS_NOPER(ISYM) + 1 591 NOPER(ISYM) = NOPER(ISYM) + 1 592 ABS_LABOP(ABS_NOPER(ISYM),ISYM) = ANGMOM(ILABEL) 593 LABOP(NOPER(ISYM),ISYM) = ANGMOM(ILABEL) 594 ELSE IF (ABSLRS_MCHD) THEN 595 ISYM2 = ISYMAX(ILABEL,2)+1 596 ABS_NOPER(ISYM2) = ABS_NOPER(ISYM2) + 1 597 ABS_LABOP(ABS_NOPER(ISYM2),ISYM2) = ANGMOM(ILABEL) 598 NOPER(ISYM2) = NOPER(ISYM2) + 1 599 LABOP(NOPER(ISYM2),ISYM2) = ANGMOM(ILABEL) 600 DO JLABEL=ILABEL,3 601 ISYM3 = ISYMAX(JLABEL,1)+1 602 ISYM4 = MULD2H(ISYM,ISYM3) 603 ABS_NOPER(ISYM4) = ABS_NOPER(ISYM4) + 1 604 NOPER(ISYM4) = NOPER(ISYM4) + 1 605 IF (ILABEL .EQ.1) THEN 606 KLABEL=JLABEL 607 ELSE 608 KLABEL=JLABEL+ILABEL 609 ENDIF 610 ABS_LABOP(ABS_NOPER(ISYM4),ISYM4) = THETA(KLABEL) 611 LABOP(NOPER(ISYM4),ISYM4) = THETA(KLABEL) 612 ENDDO 613 END IF 614 END IF 615 END DO 616 IF (ABSLRS_NSCD) THEN 617 !stop if we are using symmetry. 618 if (NSYM.gt.1) then 619 WRITE (LUPRI,*) 620 & 'NSCD not yet working symmetry: remove sym and resubmit' 621 CALL QUIT('NSCD not fully working with symmetry') 622 end if 623 CALL GPOPEN(LUPROP,'AOPROPER','UNKNOWN',' ', 624 & 'UNFORMATTED',IDUMMY,.FALSE.) 625 333 CONTINUE 626 IF (NXTLAB(LABFND,RTNLBL,LUPROP)) THEN 627 IF (LABFND(1:3) .EQ. 'PSO') THEN 628 !IND is symmetry of operator, resumed from file. Sonia 629 IND = (ICHAR(RTNLBL(1)(1:1))-ICHAR('0')) 630 ABS_NOPER(IND) = ABS_NOPER(IND) + 1 631 NOPER(IND) = NOPER(IND) + 1 632 ABS_LABOP(ABS_NOPER(IND),IND) = LABFND 633 LABOP(NOPER(IND),IND) = LABFND 634 ENDIF 635 GO TO 333 636 END IF 637 CALL GPCLOSE(LUPROP,'KEEP') 638 END IF 639C 640 IF (ABSLRS_INTERVAL) THEN 641 IF( ABS_FREQ_INTERVAL(1).GT.ABS_FREQ_INTERVAL(2) .OR. 642 & (ABS_FREQ_INTERVAL(2)-ABS_FREQ_INTERVAL(1)).LT. 643 & ABS_FREQ_INTERVAL(3) 644 & .OR. ABS_FREQ_INTERVAL(3).LE.0.0D0 ) THEN 645 CALL QUIT('.FREQ_INTERVAL not specify correctly') 646 END IF 647 END IF 648 IF (ABSQRS_INTERVAL) THEN 649 IF( ABS_BFREQ_INTERVAL(1).GT.ABS_BFREQ_INTERVAL(2) .OR. 650 & (ABS_BFREQ_INTERVAL(2)-ABS_BFREQ_INTERVAL(1)).LT. 651 & ABS_BFREQ_INTERVAL(3) 652 & .OR. ABS_BFREQ_INTERVAL(3).LE.0.0D0 ) THEN 653 CALL QUIT('.BFREQI not specified correctly') 654 END IF 655 DSMALL=1.0000001 656 ABS_NFREQ_BETA_B = 1 + 657 & INT(DSMALL*(ABS_BFREQ_INTERVAL(2)-ABS_BFREQ_INTERVAL(1))/ 658 & ABS_BFREQ_INTERVAL(3)) 659 DO I=1,ABS_NFREQ_BETA_B 660 ABS_FREQ_BETA_B(I)=ABS_BFREQ_INTERVAL(1) + 661 & (I-1)*ABS_BFREQ_INTERVAL(3) 662 END DO 663 ENDIF 664 IF (ABSQRS_CINTERVAL) THEN 665 IF( ABS_CFREQ_INTERVAL(1).GT.ABS_CFREQ_INTERVAL(2) .OR. 666 & (ABS_CFREQ_INTERVAL(2)-ABS_CFREQ_INTERVAL(1)).LT. 667 & ABS_CFREQ_INTERVAL(3) 668 & .OR. ABS_CFREQ_INTERVAL(3).LE.0.0D0 ) THEN 669 CALL QUIT('.CFREQI not specified correctly') 670 END IF 671 DSMALL=1.0000001 672 ABS_NFREQ_BETA_C = 1 + 673 & INT(DSMALL*(ABS_CFREQ_INTERVAL(2)-ABS_CFREQ_INTERVAL(1))/ 674 & ABS_CFREQ_INTERVAL(3)) 675 DO I=1,ABS_NFREQ_BETA_C 676 ABS_FREQ_BETA_C(I)=ABS_CFREQ_INTERVAL(1) + 677 & (I-1)*ABS_CFREQ_INTERVAL(3) 678 END DO 679 ENDIF 680C 681 IF (ABS_OLDCPP) THEN 682 ABSORP=.TRUE. 683 ABSLRS=.FALSE. 684 WRITE(LUPRI,'(/A)')'.OLDCPP specify' 685 WRITE(LUPRI,'(/A)')'the old CPP solver is used' 686 ELSE 687 CALL AROUND('Variable settings for absorption calculation') 688C 689 IF (ABSLRS_ALPHA) THEN 690 WRITE (LUPRI,'(A,L4)') 691 & ' Linear polarizability calculation requested: 692 & ABSLRS_ALPHA =',ABSLRS_ALPHA 693 IF (.NOT.ABSLRS_INTERVAL) WRITE(LUPRI,'(A,5(4F12.8,/,16X))') 694 & ' at frequencies:',(ABS_FREQ_ALPHA(I),I=1,ABS_NFREQ_ALPHA) 695 END IF 696C 697 IF (ABSLRS_BETA) THEN 698 WRITE (LUPRI,'(A,L4)') 699 & ' Nonlinear polarizability calculation requested:ABSLRS_BETA =', 700 & ABSLRS_BETA 701 WRITE(LUPRI,'(A,5(4F12.8))') 702 & ' B-FREQ:',(ABS_FREQ_BETA_B(I),I=1,ABS_NFREQ_BETA_B) 703 WRITE(LUPRI,'(A,5(4F12.8))') 704 & ' C-FREQ:',(ABS_FREQ_BETA_C(I),I=1,ABS_NFREQ_BETA_C) 705 END IF 706C 707 IF (ABSLRS_GAMMA .OR. ABSLRS_IDRI) THEN 708 IF (ABSLRS_GAMMA) THEN 709 WRITE (LUPRI,'(2A,L4)') 710 & ' Nonlinear polarizability calculation requested:ABSLRS_GAMMA ', 711 & '=',ABSLRS_GAMMA 712 ELSE IF (ABSLRS_IDRI) THEN 713 WRITE (LUPRI,'(2A,L4)') 714 & ' Nonlinear polarizability calculation requested:ABSLRS_IDRI ', 715 & '=',ABSLRS_IDRI 716 END IF 717 WRITE(LUPRI,'(A,5(4F12.8))') 718 & ' B-FREQ:',(ABS_FREQ_GAMMA_B(I),I=1,ABS_NFREQ_GAMMA_B) 719 WRITE(LUPRI,'(A,5(4F12.8))') 720 & ' C-FREQ:',(ABS_FREQ_GAMMA_C(I),I=1,ABS_NFREQ_GAMMA_C) 721 WRITE(LUPRI,'(A,5(4F12.8))') 722 & ' D-FREQ:',(ABS_FREQ_GAMMA_D(I),I=1,ABS_NFREQ_GAMMA_D) 723 END IF 724cC 725 IF (ABSLRS_MCD) THEN 726 WRITE (LUPRI,'(A,L4)') 727 & ' Magnetic circular dichroism requested : ABSLRS_MCD = ', 728 & ABSLRS_MCD 729 WRITE(LUPRI,'(A,5(4F12.8))') 730 & ' at frequencies:',(ABS_FREQ_BETA_B(I), 731 & I=1,ABS_NFREQ_BETA_B) 732 END IF 733 IF (ABSLRS_MCHD) THEN 734 WRITE (LUPRI,'(A,L4)') 735 & ' Magneto-chiral circular dichroism requested: ABSLRS_MCHD=', 736 & ABSLRS_MCHD 737 WRITE(LUPRI,'(A,5(4F12.8))') 738 & ' at frequencies:',(ABS_FREQ_BETA_B(I), 739 & I=1,ABS_NFREQ_BETA_B) 740 END IF 741cC 742 IF (ABSLRS_NSCD) THEN 743 WRITE (LUPRI,'(A,L4)') 744 & ' Nuclear-spin circular dichroism/OR requested: ABSLRS_NSCD=', 745 & ABSLRS_NSCD 746 WRITE(LUPRI,'(A,5(4F12.8))') 747 & ' at frequencies:',(ABS_FREQ_BETA_B(I), 748 & I=1,ABS_NFREQ_BETA_B) 749 END IF 750! 751 WRITE(LUPRI,'(A,L4)') 752 & ' Absorption calc over frequency interval : 753 & ABSLRS_INTERVAL =', ABSLRS_INTERVAL 754 IF (ABSLRS_INTERVAL) THEN 755c IF (NCONF.GT.1 ) THEN 756c CALL QUIT('Error occured! Not implemented!') 757c ELSE 758 WRITE(LUPRI,'(A,I4)') 759 & ' Number of frequencies : ABS_NFREQ_INTERVAL=' 760 & ,ABS_NFREQ_INTERVAL 761 ENDIF 762 WRITE(LUPRI,'(3(A,F8.5))') 763 & ' Start:',ABS_FREQ_INTERVAL(1),'Stop:',ABS_FREQ_INTERVAL(2), 764 & ' Step:',ABS_FREQ_INTERVAL(3) 765c END IF 766C 767 WRITE(LUPRI,'(A,F9.6,F8.1)') 768 & ' Damping parameter (a.u. and cm-1) : DAMPING =' 769 & ,ABS_DAMP,ABS_DAMP*XTKAYS 770 WRITE(LUPRI,'(A,1P,D8.1)') 771 & ' Threshold of convergence in LR solver : THCLR =' 772 & ,ABS_THCLR 773 WRITE(LUPRI,'(A,I4)') 774 & ' Maximum iterations in complex LR solver : MAXIT =' 775 & ,ABS_MAXITER 776 WRITE(LUPRI,'(A,I4)') 777 & ' Maximum size of reduced space : MAXRM =' 778 & ,ABS_MAXRM 779 WRITE(LUPRI,'(A,I4)') 780 & ' Print level in absorption modules : PRINT =' 781 & ,IPRABSLRS 782 END IF 783 ENDIF 784 CALL ABSORP_INTERFACE 785C 786C Messages 787C 788 1000 FORMAT(/, 789 & ' Keyword:',A8,' ignored since an ALPHA calc not specified.', 790 & /,' If asked for, the .ALPHA keyword should be placed first', 791 & /,' of *ABSORPTION input keys.', 792 & /,' The program will continue.') 793 794 1010 FORMAT(/, 795 & ' Keyword:',A8,' ignored since a BETA or GAMMA calc not ' 796 & 'specified.', 797 & /,' If asked for, the .BETA or .GAMMA keyword should be placed ' 798 & 'first of *ABSORPTION input keys.', 799 & /,' The program will continue.') 800 801 1020 FORMAT(/, 802 & ' Keyword:',A8,' ignored,calc with old cpp not specified.', 803 & /,' If asked for, the .OLDCPP keyword should be placed first', 804 & /,' of *ABSORPTION input keys.', 805 & /,' The program will continue.') 806c 807C 808C End of ABSORP_INPUT 809C 810 RETURN 811 END 812 SUBROUTINE ABSLRSCALC(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 813 & XINDX,WRK,LWRK) 814C 815C PURPOSE: 816C Driver routine for the computation of response properties including 817C absorption 818C 819C Used from: 820C infvar.h : MAXWOP 821C inforb.h : NSYM 822C iradef.h : IRAT 823#include "implicit.h" 824#include "dummy.h" 825#include "priunit.h" 826#include "iratdef.h" 827#include "abslrs.h" 828#include "abscrs.h" 829#include "infvar.h" 830#include "inforb.h" 831C 832 DIMENSION CMO(*),UDV(*),PV(*),FOCK(*),FC(*),FV(*),FCAC(*),H2AC(*) 833 DIMENSION XINDX(*),WRK(LWRK) 834 LOGICAL EX 835C 836 KFREE = 1 837 LFREE = LWRK 838C 839 LUABSVECS = -1 840 CALL GPINQ('ABSVECS','EXIST',EX) 841 IF (.NOT.EX) THEN 842 CALL GPOPEN(LUABSVECS,'ABSVECS','NEW',' ',' ',IDUMMY,.FALSE.) 843 WRITE(LUABSVECS) 'EOFLABEL' 844 ELSE 845 CALL GPOPEN(LUABSVECS,'ABSVECS','OLD',' ',' ',IDUMMY,.FALSE.) 846 END IF 847C 848 CALL MEMGET2('REAL','MJWOP',KMJWOP,(2*8*MAXWOP + 1)/IRAT,WRK, 849 & KFREE,LFREE) 850C 851C Determine the QR that is to be calculated 852C ========================================= 853C 854 IF (ABSLRS_BETA .OR. ABSLRS_MCD .OR. ABSLRS_MCHD .OR. ABSLRS_NSCD) 855 & THEN 856 CALL ABS_BETA_SETUP 857 END IF 858 859C 860C Determine the CR that is to be calculated 861C ========================================= 862C 863 IF (ABSLRS_GAMMA .OR. ABSLRS_IDRI) THEN 864 CALL ABS_GAMMA_SETUP 865 END IF 866C 867C Allocate memory for results of linear absorption calc 868C ===================================================== 869C 870 DSMALL=1.0000001 871 IF (ABSLRS_INTERVAL) THEN 872 ABS_NFREQ_INTERVAL = 1 + 873 & INT(DSMALL*(ABS_FREQ_INTERVAL(2)-ABS_FREQ_INTERVAL(1))/ 874 & ABS_FREQ_INTERVAL(3)) 875 ELSE 876 ABS_NFREQ_INTERVAL = ABS_NFREQ_ALPHA 877 END IF 878C 879 CALL MEMGET('REAL',KRES,2*ABS_NFREQ_INTERVAL*3*3*4,WRK, 880 & KFREE,LFREE) 881 CALL DZERO(WRK(KRES),2*ABS_NFREQ_INTERVAL*3*3*4) 882C 883C Determine linear response vectors 884C ================================= 885C 886 CALL AROUND('Solving Linear Response Equations') 887C default: Use the new SCF solver (see Reference C) 888 CALL ABS_LRS(NSYM,CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 889 & XINDX,WRK(KMJWOP),WRK(KRES),WRK(KFREE),LFREE) 890C 891C Determine double-indexed response vectors for cubic response 892C ============================================================ 893C 894 IF (ABSLRS_GAMMA .OR. ABSLRS_IDRI) THEN 895 CALL ABS_NXY(CMO, UDV, PV, XINDX, WRK(KMJWOP), FC, FCAC, FV, 896 & H2AC, FOCK, WRK(KRES), WRK(KFREE), LFREE) 897 END IF 898C 899C Evaluation of QRF 900C ================= 901C 902 IF (ABSLRS_BETA .OR. ABSLRS_MCD .OR. ABSLRS_MCHD .OR. ABSLRS_NSCD) 903 & THEN 904 CALL AROUND('Evaluate Quadratic Response Functions') 905 CALL ABSQRF(CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 906 & XINDX,WRK(KMJWOP),WRK(KFREE),LFREE) 907 END IF 908C 909C Evaluation of CRF 910C ================= 911C 912 IF (ABSLRS_GAMMA .OR. ABSLRS_IDRI) THEN 913 CALL AROUND('Evaluate Cubic Response Functions') 914 CALL ABS_CRF(CMO, UDV, PV, FOCK, FC, FV, FCAC, H2AC, 915 & XINDX, WRK(KMJWOP), WRK(KFREE), LFREE) 916 END IF 917C 918C Print obtained response results 919C =============================== 920C 921 CALL ABSLRSRESULT(WRK(KMJWOP),WRK(KRES),WRK(KFREE),LFREE) 922c 923 CALL GPCLOSE(LUABSVECS,'KEEP') 924C 925 RETURN 926 END 927 928 SUBROUTINE ABS_LRS(NSYM,CMO,UDV,PV,FOCK,FC,FV,FCAC,H2AC, 929 & XINDX,MJWOP,RESLRF,WRK,LWRK) 930C 931C PURPOSE: 932C Solve the linear response equations and store response vectors on 933C disk. 934C 935#include "implicit.h" 936#include "dummy.h" 937#include "priunit.h" 938#include "abslrs.h" 939c 940 LOGICAL CONVERGED 941 CHARACTER*8 LABEL, BLANK 942 PARAMETER (BLANK=' ') 943 DIMENSION CMO(*),UDV(*),PV(*),FOCK(*),FC(*),FV(*),FCAC(*),H2AC(*) 944c DIMENSION XINDX(*),MJWOP(2,MAXWOP,8) 945 DIMENSION XINDX(*),MJWOP(*) 946 DIMENSION RESLRF(2,ABS_NFREQ_INTERVAL,3,3,4),WRK(LWRK) 947C local 948 REAL*8 thd 949 logical negf 950C 951C 952 IF (ABS_NGD) THEN 953 NGD=ABS_NFREQ_INTERVAL 954 ELSE 955 NGD=1 956 ENDIF 957 IF (IPRABSLRS.GT.0) CALL TIMER('START ',TIMSTR,TIMEND) 958C 959 negf=.false. 960 DO ISYM=1,NSYM 961 IF (ABS_NOPER(ISYM).GT.0) THEN 962C 963 LFREE=LWRK 964 KFREE=1 965C 966 CALL ABSINTER(ISYM,UDV,FOCK,FC,FV,FCAC,H2AC,XINDX,MJWOP, 967 & KZVAR,WRK(KFREE),LFREE) 968C 969 KGD = 1 970 kfreqs = KGD + 4*NGD*KZVAR 971 KXSOL = kfreqs + abs_nfreq_interval 972 KFREE = KXSOL + 4*ABS_NFREQ_INTERVAL*KZVAR 973 LFREE = LWRK - KFREE 974C 975 ABS_KLRED(1)=0 976 ABS_KLRED(2)=0 977C 978 LUSB = -1 979 LUAB = -1 980 LUSS = -1 981 LUAS = -1 982C 983 CALL GPOPEN(LUSB,'ABS_SB','NEW',' ',' ',IDUMMY,.FALSE.) 984 CALL GPOPEN(LUAB,'ABS_AB','NEW',' ',' ',IDUMMY,.FALSE.) 985 CALL GPOPEN(LUSS,'ABS_SS','NEW',' ',' ',IDUMMY,.FALSE.) 986 CALL GPOPEN(LUAS,'ABS_AS','NEW',' ',' ',IDUMMY,.FALSE.) 987C 988 LUE1RED = -1 989 LUE2RED = -1 990 LUSRED = -1 991 CALL GPOPEN(LUE1RED,'ABS_E1RED','NEW', 992 & ' ',' ',IDUMMY,.FALSE.) 993 CALL GPOPEN(LUE2RED,'ABS_E2RED','NEW', 994 & ' ',' ',IDUMMY,.FALSE.) 995 CALL GPOPEN(LUSRED ,'ABS_SRED' ,'NEW', 996 & ' ',' ',IDUMMY,.FALSE.) 997C 998 DO IOPER=1,ABS_NOPER(ISYM) 999 LABEL=ABS_LABOP(IOPER,ISYM) 1000 IF (ABSLRS_INTERVAL) THEN 1001 DO I=1,ABS_NFREQ_INTERVAL 1002 ABS_FREQ_ALPHA(I)=ABS_FREQ_INTERVAL(1) + 1003 & (I-1)*ABS_FREQ_INTERVAL(3) 1004 END DO 1005 END IF 1006 IF (IPRABSLRS.GE.0) THEN 1007 WRITE(LUPRI,'(/2A,/3A,I2,A/,A,5(4F12.8,/,11X))') 1008 & ' ABSLRS -- Solving linear response equations', 1009 & '( E[2] - {w+iW}*S[2] ) N = B[1] ', 1010 & ' ABSLRS -- for operator ', LABEL, 1011 & ' of symmetry',ISYM,' at frequencies:', 1012 & ' ABSLRS --',(ABS_FREQ_ALPHA(I), 1013 & I=1,ABS_NFREQ_INTERVAL) 1014 END IF 1015C 1016 CALL DZERO(WRK(KGD),4*KZVAR*NGD) 1017 CALL DZERO(WRK(KXSOL),4*KZVAR*ABS_NFREQ_INTERVAL) 1018 CALL GETGPV(LABEL,FC,FV,CMO,UDV,PV,XINDX,ANTSYM, 1019 & WRK(KFREE),LFREE) 1020 CALL DCOPY(2*KZVAR,WRK(KFREE),1,WRK(KGD),1) 1021 GDNRM=DNRM2(2*KZVAR,WRK(KGD),1) 1022 RES0=0.0d0 1023 IF (GDNRM .LE. 1.0d-6) THEN 1024 WRK(KXSOL:4*KZVAR-1)=0.0d0 1025 DO K=1,ABS_NFREQ_INTERVAL 1026C CALL WRITE_XVEC(LUABSVECS,4*KZVAR,WRK(KXSOL),LABEL, 1027C & ABS_FREQ_ALPHA(K),RES0) 1028 CALL WRITE_XVEC2(LUABSVECS,4*KZVAR,WRK(KXSOL),LABEL, 1029 & BLANK,ABS_FREQ_ALPHA(K),0.0d0,RES0) 1030 ENDDO 1031 WRITE(LUPRI,*) 'Vectors equal to zero' 1032 ELSE 1033 nfreqs=abs_nfreq_interval 1034 do i=1,abs_nfreq_interval 1035 if (abs_freq_alpha(i) .lt. 0.d0) negf=.true. 1036 wrk(kfreqs-1+i)=abs(abs_freq_alpha(i)) 1037 enddo 1038 if (negf) then 1039 call gpdsrt(nfreqs,wrk(kfreqs),thd) 1040 endif 1041 CALL ABS_CTL(LABEL,KZVAR,WRK(KGD),NGD,WRK(KXSOL), 1042 & wrk(kfreqs),nfreqs,LUABSVECS, 1043 & MJWOP,CMO,UDV,FC,FCAC,FV,PV,XINDX,RESLRF, 1044 & WRK(KFREE),LFREE) 1045 ENDIF 1046C 1047 CALL GET_LRF(LUABSVECS,LABEL,ABS_NFREQ_INTERVAL, 1048 & ABS_FREQ_ALPHA,FC,FV,CMO,UDV,PV,XINDX, 1049 & RESLRF,KZVAR,WRK(KFREE),LFREE) 1050C 1051 END DO 1052C 1053 CALL GPCLOSE(LUSB,'DELETE') 1054 CALL GPCLOSE(LUAB,'DELETE') 1055 CALL GPCLOSE(LUSS,'DELETE') 1056 CALL GPCLOSE(LUAS,'DELETE') 1057C 1058 CALL GPCLOSE(LUE1RED,'DELETE') 1059 CALL GPCLOSE(LUE2RED,'DELETE') 1060 CALL GPCLOSE(LUSRED, 'DELETE') 1061C 1062 END IF 1063 END DO 1064 IF (IPRABSLRS.GT.0) CALL TIMER('ABSVEC',TIMSTR,TIMEND) 1065C 1066 RETURN 1067 END 1068 SUBROUTINE ABSINTER(ISYM,UDV,FOCK,FC,FV,FCAC,H2AC,XINDX,MJWOP, 1069 & ABS_KZVAR,WRK,LWRK) 1070C 1071#include "implicit.h" 1072#include "priunit.h" 1073C 1074C Used from: 1075C infvar.h : MAXWOP 1076C inforb.h : NSYM 1077#include "abslrs.h" 1078#include "absorp.h" 1079!MAXWOP 1080#include "infvar.h" 1081!NSYM 1082#include "inforb.h" 1083c 1084#include "infrsp.h" 1085#include "wrkrsp.h" 1086#include "infdim.h" 1087#include "qrinf.h" 1088#include "maxorb.h" 1089#include "maxash.h" 1090#include "infind.h" 1091#include "infpri.h" 1092c 1093 DIMENSION UDV(*),FOCK(*),FC(*),FV(*),FCAC(*),H2AC(*) 1094 DIMENSION XINDX(*),MJWOP(2,MAXWOP,8),WRK(LWRK) 1095 INTEGER ABS_KZVAR,ISYM 1096C local 1097 INTEGER KFREE,LFREE 1098c 1099 KFREE=1 1100 LFREE=LWRK 1101C 1102 KSYMOP = ISYM 1103 ABS_GRADSYM = ISYM 1104 ABS_NSYM = NSYM 1105 DO I = 1,NSYM 1106 ABS_NASH(I) = NASH(I) 1107 ABS_NISH(I) = NISH(I) 1108 ABS_NORB(I) = NORB(I) 1109 DO J = 1,NSYM 1110 ABS_MULD2H(I,J) = MULD2H(I,J) 1111 ENDDO 1112 ENDDO 1113 LUABSPRI=LUPRI 1114 CALL RSPVAR(UDV,FOCK,FC,FV,FCAC,H2AC,XINDX,WRK(KFREE),LFREE) 1115 CALL SETZY(MJWOP) 1116 ABS_KZVAR=KZVAR 1117C 1118 DAMPING=ABS_DAMP 1119 THCLR_ABSORP=ABS_THCLR 1120c 1121 MAXRM=ABS_MAXRM 1122 NFREQ_INTERVAL=ABS_NFREQ_INTERVAL 1123 IPRABS=IPRABSLRS 1124C 1125 RETURN 1126 END 1127 SUBROUTINE ABS2INTER(MJWOP,ABS_KZVAR) 1128C 1129C Used from: 1130C infvar.h: MAXWOP 1131C inforb.h: NSYM 1132C 1133#include "implicit.h" 1134#include "priunit.h" 1135#include "abslrs.h" 1136#include "absorp.h" 1137#include "infvar.h" 1138#include "inforb.h" 1139c 1140#include "infrsp.h" 1141#include "wrkrsp.h" 1142#include "infdim.h" 1143#include "qrinf.h" 1144#include "maxorb.h" 1145#include "maxash.h" 1146#include "infind.h" 1147#include "infpri.h" 1148c 1149 DIMENSION MJWOP(2,MAXWOP,8) 1150 INTEGER ABS_KZVAR,ISYM 1151c 1152 ABS_GRADSYM = KSYMOP 1153 ABS_NSYM = NSYM 1154 DO I = 1,NSYM 1155 ABS_NASH(I) = NASH(I) 1156 ABS_NISH(I) = NISH(I) 1157 ABS_NORB(I) = NORB(I) 1158 DO J = 1,NSYM 1159 ABS_MULD2H(I,J) = MULD2H(I,J) 1160 ENDDO 1161 ENDDO 1162 LUABSPRI=LUPRI 1163 CALL SETZY(MJWOP) 1164 ABS_KZVAR=KZVAR 1165C 1166 ABS_DAMP=DAMPING 1167 ABS_THCLR=THCLR_ABSORP 1168c 1169 ABS_MAXRM=MAXRM 1170 ABS_NFREQ_INTERVAL=NFREQ_INTERVAL 1171 ABS_NFREQ_ALPHA=NFREQ_ALPHA 1172 IPRABSLRS=IPRABS 1173 DO I=1,ABS_NFREQ_INTERVAL 1174 ABS_FREQ_ALPHA(I)=FREQ_ALPHA(I) 1175 ENDDO 1176C 1177 RETURN 1178 END 1179 SUBROUTINE ABS_LINE2(KZVAR,XTMP,ZYMB,ISYMDM,ISYMB,NNMAX, 1180 & CMO,FC,FCAC,FV,UDV,XINDX,MJWOP,WRK,LWRK) 1181C 1182c 1183C PURPOSE: 1184C 1185C 1186c 1187#include "implicit.h" 1188#include "thrzer.h" 1189#include "dummy.h" 1190#include "priunit.h" 1191C 1192C Used from: 1193C inforb.h: NORB,NORBT,NASHT,NISHT,N2BASX,NBAST,NBAS,NSYM,MULD2H,ICMO 1194C infinp.h: DODFT,HSROHF,DIRFCK 1195C dftcom.h: DFT 1196C thrzer.h: THZER 1197C maxorb.h: MXCORB (needed for inforb) 1198#include "abslrs.h" 1199#include "maxorb.h" 1200#include "inforb.h" 1201#include "infinp.h" 1202#include "dftcom.h" 1203#include "pcmlog.h" 1204#include "gnrinf.h" 1205C 1206 INTEGER NNMAX,KZVAR 1207 DIMENSION XTMP(2*KZVAR,NNMAX) 1208 DIMENSION ZYMB(NORBT,NORBT,NNMAX) 1209 DIMENSION CMO(*),FC(*),FCAC(*),UDV(*),FV(*),WRK(LWRK) 1210 DIMENSION ISYMDM(2*NNMAX),IFCTYP(2*NNMAX),KDMAO(NNMAX), 1211 & KFMAO(NNMAX), 1212 & KFMMO(NNMAX),KDVAO(NNMAX),KFVAO(NNMAX),KFVMO(NNMAX) 1213C 1214 DIMENSION MJWOP(*),XINDX(*) 1215 PARAMETER (D2 = 2.0D0) 1216 PARAMETER (D1 = 1.0D0) 1217C local 1218 INTEGER KZYVAR 1219C 1220 CALL QENTER('ABS_LINE2') 1221C 1222C Allocate work space for density and Fock matrices 1223C 1224 KFREE=1 1225 LFREE=LWRK 1226 KZYVAR=2*KZVAR 1227 IF (NASHT.GT.0) THEN 1228 NTMAX=2*NNMAX 1229 ELSE 1230 NTMAX=NNMAX 1231 ENDIF 1232c 1233 CALL MEMGET('REAL',KDENS,NTMAX*N2BASX,WRK,KFREE,LFREE) 1234 CALL MEMGET('REAL',KFOCK,NTMAX*N2BASX,WRK,KFREE,LFREE) 1235 1236c 1237 DO I=1,NNMAX 1238 KDMAO(I) = KDENS +(I-1)*N2BASX 1239 KFMAO(I) = KFOCK +(I-1)*N2BASX 1240 ENDDO 1241 IF (NASHT.GT.0) THEN 1242 DO I=1,NNMAX 1243 KDVAO(I)=KDENS+(NNMAX+I-1)*N2BASX 1244 KFVAO(I)=KFOCK+(NNMAX+I-1)*N2BASX 1245 ENDDO 1246 ENDIF 1247C 1248 CALL DZERO(WRK(KDMAO(1)),NTMAX*N2BASX) 1249c 1250 NKD=0 1251C 1252 DO I=1,NNMAX 1253 CALL GTZYMT(1,XTMP(1,I),KZYVAR,ISYMB,ZYMB(1,1,I),MJWOP) 1254C 1255C Construct modified density matrices in AO basis 1256C 1257 IF (NASHT.GE.0) THEN 1258 CALL ABS_ZYTRA(ZYMB(1,1,I),CMO,UDV,WRK(KDMAO(I)), 1259 & WRK(KDVAO(I)),WRK(KFREE),LFREE) 1260 ELSE 1261 CALL ABS_ZYTRA(ZYMB(1,1,I),CMO,UDV,WRK(KDMAO(I)), 1262 & VDUMMY,WRK(KFREE),LFREE) 1263 ENDIF 1264 END DO 1265C 1266 IF (NISHT .GT. 0) THEN 1267 CALL DSCAL(NNMAX*N2BASX,D2,WRK(KDMAO(1)),1) 1268 ENDIF 1269 IF ((NASHT.GT.0 .AND. DODFT).OR.HSROHF) THEN 1270 CALL DAXPY(NNMAX*N2BASX,D1,WRK(KDVAO(1)),1,WRK(KDMAO(1)),1) 1271 CALL DSCAL(NNMAX*N2BASX,-D1,WRK(KDVAO(1)),1) 1272 END IF 1273 IF (NASHT .GT. 0) THEN 1274 DO I=1,NNMAX 1275 ISYMDM(NNMAX+I)=ISYMB 1276 ENDDO 1277 ENDIF 1278C 1279C Construct Fock matrices in AO basis 1280C 1281C IFCTYP = XY 1282C X indicates symmetry about diagonal 1283C X = 0 No symmetry 1284C X = 1 Symmetric 1285C X = 2 Anti-symmetric 1286C Y indicates contributions 1287C Y = 0 no contribution ! 1288C Y = 1 Coulomb 1289C Y = 2 Exchange 1290C Y = 3 Coulomb + Exchange 1291C 1292C Check if density matrix is unsymmetric (IX=0), 1293C symmetric (IX=10), antisymmetric (IX=20), or zero matrix (IX=30) 1294C to threshold THRZER 1295C 1296 DO I = 1,NTMAX 1297 IX = 10 * MATSYM(NBAST,NBAST,WRK(KDENS+(I-1)*N2BASX),THRZER) 1298 IF (IX .EQ. 30) THEN 1299C if zero density matrix, do nothing 1300 IFCTYP(I) = 0 1301 ELSE IF (IX .EQ. 10) THEN 1302C if symmetric, Coulomb+exchange 1303 IFCTYP(I) = IX + 3 1304 ELSE IF (IX .EQ. 20) THEN 1305C if antisymmetric density matrix, only exchange 1306 IFCTYP(I) = IX + 2 1307 ELSE IF (IX .EQ. 0) THEN 1308C if general, Coulomb+exchange 1309 IFCTYP(I) = IX + 3 1310 ELSE 1311 WRITE(LUPRI,'(/A,I3)') 'ABSLIN: error'// 1312 & ' density matrix symmetry:',I 1313 CALL QUIT('ABSLIN: error density matrix symmetry') 1314 END IF 1315 END DO 1316C 1317 IF ((NASHT.GT.0.AND.DODFT).OR.HSROHF) THEN 1318C Change IFCTYP for all Fa matrices 1319 DO I=1,NNMAX 1320c IF (TRPLET) THEN 1321c IFCTYP(I)=10*(IFCTYP(I)/10) + 2 1322c IFCTYP(NNMAX+I)=10*(IFCTYP(NNMAX+I)/10) + 3 1323c ELSE 1324 IFCTYP(I)=10*(IFCTYP(I)/10) + 3 1325 IFCTYP(NNMAX+I)=10*(IFCTYP(NNMAX+I)/10) + 2 1326c END IF 1327 END DO 1328 END IF 1329C 1330 CALL DZERO(WRK(KFMAO(1)),NTMAX*N2BASX) 1331 CALL SIRFCK(WRK(KFMAO(1)),WRK(KDMAO(1)),NTMAX,ISYMDM,IFCTYP, 1332 & DIRFCK,WRK(KFREE),LFREE) 1333C 1334 IF ((NASHT.GT.0 .AND. DODFT).OR.HSROHF) THEN 1335C 1336C Calculated matrices are (FC+FV,-FV(exch)=FV-Q) 1337C Input fock matrices (form SIRIFC) are of the form 1338C (FC+Q,FV-Q), where Q=FV(coul) + 2*FV(exch) 1339C adapt to this form (which works in RSPORB) 1340C 1341C FC+Q 1342 CALL DAXPY(NNMAX*N2BASX,-D1,WRK(KFVAO(1)),1,WRK(KFMAO(1)),1) 1343 END IF 1344C Transform Fock matrices to MO basis 1345C Reuse memory used for density matrices 1346C 1347 DO I=1,NNMAX 1348 KFMMO(I) = KDENS +(I-1)*NORBT*NORBT 1349 ENDDO 1350 IF (NASHT.GT.0) THEN 1351 DO I=1,NNMAX 1352 KFVMO(I) = KDENS +(NNMAX+I-1)*NORBT*NORBT 1353 ENDDO 1354 ENDIF 1355C 1356 NKD=0 1357 NASIM=4 1358 KINT=INT(NNMAX/NASIM) 1359 KMOD=MOD(NNMAX,NASIM) 1360 KTOT=KINT 1361 IF (KMOD.NE.0) THEN 1362 KTOT=KINT+1 1363 ENDIF 1364 DO K=1,KTOT 1365 IF (K.GT.KINT) NASIM=KMOD 1366 DO I=1,NASIM 1367c DO I=1,NNMAX 1368 NKD=NKD+1 1369 CALL FAOMO(ISYMB,WRK(KFMAO(NKD)),WRK(KFMMO(NKD)),CMO, 1370 & WRK(KFREE),LFREE) 1371 IF (NASHT .GT.0) THEN 1372 CALL DZERO(WRK(KFVMO(NKD)),NORBT*NORBT) 1373 DO 300 ISYM1=1,NSYM 1374 ISYM2 = MULD2H(ISYM1,ABS_GRADSYM) 1375 NORB1 = NORB(ISYM1) 1376 NORB2 = NORB(ISYM2) 1377 IF (NORB1.EQ.0 .OR. NORB2.EQ.0) GO TO 300 1378 CALL AUTPV(ISYM1,ISYM2,CMO(ICMO(ISYM1)+1), 1379 & CMO(ICMO(ISYM2)+1),WRK(KFVAO(NKD)),NBAS,NBAST, 1380 & WRK(KFVMO(NKD)),NORB,NORBT,WRK(KFREE),LFREE) 1381 300 CONTINUE 1382 END IF 1383 ENDDO 1384c add DFT contribution 1385 IF (DODFT) THEN 1386C Special case of alpha density only (old code) 1387 IF((NISHT.EQ.0) .AND. (NASHT.GT.0)) THEN 1388 DO IOSIM = 1, NASIM 1389 call dft_lin_respab(WRK(KFMMO(NKD-NASIM+IOSIM)), 1390 & WRK(KFVMO(NKD-NASIM+IOSIM)),CMO, 1391 & ZYMB(1,1,NKD-NASIM+IOSIM),.FALSE., 1392 & ISYMB,WRK(KFREE),LFREE,IPRDFT) 1393 END DO 1394 ENDIF 1395C Proceed normal way otherwise 1396 IF ((NASHT.GT.0) .AND. (NISHT.GT.0)) THEN 1397 call dft_lin_respab_b(NASIM,WRK(KFMMO(NKD-NASIM+1)), 1398 & WRK(KFVMO(NKD-NASIM+1)),CMO, 1399 & ZYMB(1,1,NKD-NASIM+1),.FALSE., 1400 & ISYMB,WRK(KFREE),LFREE,IPRDFT) 1401 ELSE 1402c call dft_lin_respf(NNMAX,WRK(KFMMO(1)),CMO,ZYMB, 1403c & .FALSE.,ISYMB,WRK(KFREE),LFREE,IPRDFT) 1404 call dft_lin_respf(NASIM,WRK(KFMMO(NKD-NASIM+1)),CMO, 1405 & ZYMB(1,1,NKD-NASIM+1), 1406 & .FALSE.,ISYMB,WRK(KFREE),LFREE,IPRDFT) 1407 END IF 1408 END IF 1409 END DO 1410 CALL DZERO(XTMP,KZYVAR*NNMAX) 1411 IF (NASHT.GT.0) THEN 1412 CALL MEMGET('REAL',KDUM,N2ORBX,WRK,KFREE,LFREE) 1413 call DZERO(WRK(KDUM),N2ORBX) 1414 CALL FCKOIN(NNMAX,FC,FV,ZYMB,WRK(KFMMO(1)),WRK(KFVMO(1))) 1415 CALL RSPORB(.TRUE.,NNMAX,FC,WRK(KFMMO(1)),WRK(KFVMO(1)), 1416 & WRK(KDUM),WRK(KDUM),UDV, 1417 & XTMP,ZYMB) 1418 ELSE 1419 CALL FCKOIN(NNMAX,FC,DUMMY,ZYMB,WRK(KFMMO(1)),VDUMMY) 1420 CALL RSPORB(.TRUE.,NNMAX,FC,WRK(KFMMO(1)),VDUMMY, 1421 & VDUMMY,VDUMMY,UDV,XTMP,ZYMB) 1422c CALL ORBEX(ISYMB,ISYMB,KZYVAR,WRK(KFMMO(I)),VDUMMY, 1423c & VDUMMY,VDUMMY,XTMP(1,I),1.0D0,UDV,MJWOP) 1424 ENDIF 1425C 1426 CALL QEXIT('ABS_LINE2') 1427 RETURN 1428 END 1429 SUBROUTINE ABS_LINE2_INTER(KZVAR,VECS,VECB,XTMP,XTEMP,ISYMB,NLT, 1430 & NNMAX,KNVEC,CMO,FC,FCAC,FV,UDV,XINDX,MJWOP, 1431 & WRK,LWRK) 1432C 1433C PURPOSE: 1434C 1435 use pelib_interface, only: use_pelib 1436#include "implicit.h" 1437#include "abslrs.h" 1438#include "thrzer.h" 1439#include "dummy.h" 1440#include "priunit.h" 1441C Used from 1442C inforb.h: NORB,NORBT,NASHT,NISHT,N2BASX,NBAST,NBAS,NSYM,MULD2H,ICMO 1443C infinp.h: DODFT,HSROHF,DIRFCK 1444C dftcom.h: DFT 1445C thrzer.h: THZER 1446C maxorb.h: MXCORB (needed for inforb) 1447#include "inforb.h" 1448#include "maxorb.h" 1449#include "infinp.h" 1450#include "dftcom.h" 1451c 1452#include "pcmlog.h" 1453#include "gnrinf.h" 1454C 1455 INTEGER NLT,NNMAX,KZVAR 1456 INTEGER KNVEC(2) 1457 DIMENSION VECS(KZVAR,NLT) 1458 DIMENSION VECB(KZVAR,NLT),XTMP(2*KZVAR,NNMAX),XTEMP(2*KZVAR,NNMAX) 1459 DIMENSION CMO(*),FC(*),FCAC(*),UDV(*),FV(*),WRK(LWRK) 1460 DIMENSION ISYMDM(2*NNMAX) 1461C 1462c DIMENSION MJWOP(2,MAXWOP,8) 1463 DIMENSION MJWOP(*),XINDX(*) 1464 PARAMETER (D2 = 2.0D0) 1465 PARAMETER (D1 = 1.0D0) 1466C local 1467 INTEGER KZYVAR 1468C 1469 CALL QENTER('ABS_LINE2') 1470C 1471C Allocate work space for density and Fock matrices 1472C 1473 KFREE=1 1474 LFREE=LWRK 1475 KZYVAR=2*KZVAR 1476c 1477 NNS=KNVEC(1) 1478 NNA=KNVEC(2) 1479C 1480 NKD=0 1481C 1482C Unpack the trial vector of specific parity 1483c 1484 DO I=1,MIN(NNS,NNA) 1485 NKD=NKD+1 1486 ISYMDM(NKD)=ISYMB 1487 CALL DCOPY(KZVAR,VECB(1,I),1,XTMP(1,NKD),1) 1488 CALL DCOPY(KZVAR,VECB(1,I),1,XTMP(KZVAR+1,NKD),1) 1489 CALL DAXPY(KZVAR,1.0d0,VECB(1,NNS+I),1,XTMP(1,NKD),1) 1490 CALL DAXPY(KZVAR,-1.0d0,VECB(1,NNS+I),1,XTMP(KZVAR+1,NKD),1) 1491 END DO 1492 IF (NNS.NE.NNA) THEN 1493 IF (NNS.GE.NNA) THEN 1494 F1=1.0d0 1495 NNT=0 1496 ELSE 1497 F1=-1.0d0 1498 NNT=NNS 1499 ENDIF 1500 DO I=MIN(NNA,NNS)+1,NNMAX 1501 NKD=NKD+1 1502 ISYMDM(NKD)=ISYMB 1503 CALL DCOPY(KZVAR,VECB(1,NNT+I),1,XTMP(1,NKD),1) 1504 CALL DCOPY(KZVAR,VECB(1,NNT+I),1,XTMP(KZVAR+1,NKD),1) 1505 CALL DSCAL(KZVAR,F1,XTMP(KZVAR+1,NKD),1) 1506 ENDDO 1507 ENDIF 1508 CALL DCOPY(2*KZVAR*NNMAX,XTMP,1,XTEMP,1) 1509c 1510 IF (ABS_BATCH) THEN 1511 if (nnmax .gt. abs_nbatchmax) then 1512 KINT=INT(NNMAX/ABS_NBATCHMAX) 1513 KMOD=MOD(NNMAX,ABS_NBATCHMAX) 1514 KTOT=KINT 1515 KNMAX=ABS_NBATCHMAX 1516 IF (KMOD.NE.0) THEN 1517 KTOT=KINT+1 1518 ENDIF 1519 else 1520 kint = 1 1521 ktot = 1 1522 knmax = nnmax 1523 endif 1524 ELSE 1525 KINT=1 1526 KTOT=1 1527 KNMAX=NNMAX 1528 ENDIF 1529 KNTOT=KNMAX 1530 CALL MEMGET('REAL',KZYMB,KNMAX*NORBT*NORBT,WRK,KFREE,LFREE) 1531 DO K=1,KTOT 1532 IF (K.GT.KINT) KNTOT=KMOD 1533 IF (IPRABSLRS.GE.1) THEN 1534 WRITE(LUABSPRI,'(A,I4,A,I4)') 1535 & 'Number of linear transformations',KNTOT,' in batch ',K 1536 END IF 1537 CALL ABS_LINE2(KZVAR,XTMP(1,(K-1)*KNMAX+1),WRK(KZYMB),ISYMDM, 1538 & ISYMB,KNTOT,CMO,FC,FCAC,FV, 1539 & UDV,XINDX,MJWOP,WRK(KFREE),LFREE) 1540 ENDDO 1541C 1542C IF (FLAG(16) .OR. PCM .OR. QM3 .OR. QMNPMM) THEN 1543 IF (EMBEDDING .OR. USE_PELIB()) THEN 1544 CALL RSPSLV(0,NNMAX,VDUMMY,XTEMP,CMO,XINDX,UDV, 1545 & XTMP,WRK(KFREE),LFREE) 1546 ENDIF 1547 IF (QMMM) THEN 1548 CALL QUIT('ERROR: CPP cannot be used with QMMM') 1549 END IF 1550C 1551 CALL DZERO(VECS,NLT*KZVAR) 1552 DO I=1,MIN(NNS,NNA) 1553 CALL DCOPY(KZVAR,XTMP(1,I),1,VECS(1,I),1) 1554 CALL DAXPY(KZVAR,1.0d0,XTMP(KZVAR+1,I),1,VECS(1,I),1) 1555 CALL DCOPY(KZVAR,XTMP(1,I),1,VECS(1,NNS+I),1) 1556 CALL DAXPY(KZVAR,-1.0d0,XTMP(KZVAR+1,I),1,VECS(1,NNS+I),1) 1557 CALL DSCAL(KZVAR,0.5d0,VECS(1,I),1) 1558 CALL DSCAL(KZVAR,0.5d0,VECS(1,NNS+I),1) 1559 ENDDO 1560 IF (MIN(NNS,NNA).NE.NNMAX) THEN 1561 IF (NNS.GE.NNA) THEN 1562 F1=1.0d0 1563 NNT=0 1564 ELSE 1565 F1=-1.0d0 1566 NNT=NNS 1567 ENDIF 1568 DO I=MIN(NNS,NNA)+1,NNMAX 1569 CALL DCOPY(KZVAR,XTMP(1,I),1,VECS(1,NNT+I),1) 1570 CALL DAXPY(KZVAR,F1,XTMP(KZVAR+1,I),1,VECS(1,NNT+I),1) 1571 CALL DSCAL(KZVAR,0.5d0,VECS(1,NNT+I),1) 1572 ENDDO 1573 ENDIF 1574 1575 IF (IPRABSLRS.GE.200) THEN 1576 WRITE(LUPRI,'(//A)') ' SIGMA VECTORS (IN ORDER REAL-SYM, 1577 & REAL-ANTYSYM., IMG-ANTYSYM., IMG-SYM.' 1578 CALL OUTPUT(VECS,1,KZVAR,1,NLT,KZVAR,NLT,1,LUPRI) 1579 ENDIF 1580C 1581 CALL QEXIT('ABS_LINE2') 1582 RETURN 1583 END 1584 SUBROUTINE ABS_ZYTRA(ZYM,CMO,DV,DXCAO,DXVAO,WRK,LWRK) 1585C 1586C 1587C Purpose: Computes DAO = CMO * ZYM[0,+;-,0] * transp(CMO) 1588C = CMO * transp( CMO * transp(ZYM[]) ) 1589C 1590C ZYM[0,+;-,0](M,N) = +ZYM(M,N), M virtual, N occupied 1591C = -ZYM(M,N), N virtual, M occupied 1592C = 0, otherwhise 1593C 1594C Input: 1595C ZYM(*) kappa matrix 1596C CMO(*) MO coefficients 1597C TMP(*) temporary matrix 1598C 1599#include "implicit.h" 1600C 1601C Used from 1602C maxorb.h: MXCORB (needed for inforb) 1603C infinp.h: NORBT,NASHT,NISHT,NBAST,N2BASX,NSYM,MULD2H,ICMO,NORB,NBAS,NASH,NISH,IORB,IBAS,ICMO,IASH 1604#include "priunit.h" 1605#include "abslrs.h" 1606#include "inforb.h" 1607 DIMENSION ZYM(NORBT,*),WRK(LWRK) 1608 DIMENSION CMO(*),DV(NASHT,*),DXCAO(*),DXVAO(*) 1609 PARAMETER (HALF = 0.5D0, D2 = 2.0D0, D4 = 4.0D0) 1610 PARAMETER (HALFM = -0.5D0, D1 = 1.D0, DM1 = -1.D0) 1611C 1612 CALL QENTER('ABS_ZYTRA') 1613C 1614C loop over symmetries 1615C 1616 KFREE=1 1617 LFREE=LWRK 1618 LDXAO1 = MAX(NASHT,NISHT)*NBAST 1619 CALL MEMGET('REAL',KDXAO1,LDXAO1,WRK,KFREE,LFREE) 1620 CALL MEMGET('REAL',KDXAO2,N2BASX,WRK,KFREE,LFREE) 1621C 1622C First DXVAO, if nasht.gt.0: 1623C 1624 IF (NASHT .GT. 0) THEN 1625C 1626C ************************************************* 1627C D-I and D-II: D = D-I - D-II Active matrices 1628C ************************************************* 1629C DXVAO-I: loop over symmetries then put result in DXVAO 1630C 1631 CALL MEMGET('REAL',KDXV ,NASHT*NORBT,WRK,KFREE,LFREE) 1632 CALL DZERO(WRK(KDXAO2),N2BASX) 1633 DO 2000 ISYM1 = 1,NSYM 1634 NASH1 = NASH(ISYM1) 1635 NISH1 = NISH(ISYM1) 1636 ISYM2 = MULD2H(ISYM1,ABS_GRADSYM) 1637 IORB1 = IORB(ISYM1) 1638 IORB2 = IORB(ISYM2) 1639 NBAS1 = NBAS(ISYM1) 1640 NBAS2 = NBAS(ISYM2) 1641 NORB1 = NORB(ISYM1) 1642 NORB2 = NORB(ISYM2) 1643 ICMO1 = ICMO(ISYM1) 1644 ICMO2 = ICMO(ISYM2) 1645C ** step 1 of active dm: 1646C ***** first calculate one-index transformation of second index 1647C ***** of DV(uv), the active density matrix. 1648C DXV(p,u) = sum(v) Bo(v,p) DV(v,u) 1649 IASH1 = IASH(ISYM1) 1650 IF (NASH1 .EQ. 0 .OR. NORB2 .EQ. 0) GO TO 2000 1651 CALL DGEMM('T','N',NORB2,NASH1,NASH1,1.D0, 1652 & ZYM(IORB1+NISH1+1,IORB2+1),NORBT, 1653 & DV(IASH1+1,IASH1+1),NASHT,0.D0, 1654 & WRK(KDXV),NORB2) 1655C ** step 2 of active dm: 1656 CALL DGEMM('N','N',NBAS2,NASH1,NORB2,1.D0,CMO(ICMO2+1), 1657 & NBAS2,WRK(KDXV),NORB2,0.D0,WRK(KDXAO1),NBAS2) 1658 IOFMOV = ICMO1 + 1 + NISH1*NBAS1 1659 IDXAO2 = KDXAO2 + IBAS(ISYM2)*NBAST + IBAS(ISYM1) 1660 CALL DGEMM('N','T',NBAS1,NBAS2,NASH1,1.D0,CMO(IOFMOV),NBAS1, 1661 & WRK(KDXAO1),NBAS2,0.D0,WRK(IDXAO2),NBAST) 1662C ** this symmetry block finished 1663 2000 CONTINUE 1664C 1665 CALL DCOPY(N2BASX,WRK(KDXAO2),1,DXVAO,1) 1666C 1667C DXVAO-II: loop over symmetries then add results to DXVAO 1668C only chnage from previous loop is that first 1669C MPATB is here MPAB 1670C 1671 CALL DZERO(WRK(KDXAO2),N2BASX) 1672 DO 3000 ISYM1 = 1,NSYM 1673 ISYM2 = MULD2H(ISYM1,ABS_GRADSYM) 1674 NASH1 = NASH(ISYM1) 1675 NISH1 = NISH(ISYM1) 1676 NASH2 = NASH(ISYM2) 1677 NISH2 = NISH(ISYM2) 1678 NBAS1 = NBAS(ISYM1) 1679 NBAS2 = NBAS(ISYM2) 1680 IORB1 = IORB(ISYM1) 1681 IORB2 = IORB(ISYM2) 1682 NORB1 = NORB(ISYM1) 1683 NORB2 = NORB(ISYM2) 1684 ICMO1 = ICMO(ISYM1) 1685 ICMO2 = ICMO(ISYM2) 1686C ** step 1 of active dm: 1687C ***** first calculate one-index transformation of second index 1688C ***** of DV(uv), the active density matrix. 1689C DXV(p,u) = sum(v) Bo(v,p) DV(v,u) 1690 IASH1 = IASH(ISYM1) 1691 IASH2 = IASH(ISYM2) 1692 IF( NASH2 .EQ. 0 .OR. NORB1 .EQ. 0) GO TO 3000 1693 CALL DGEMM('N','N',NORB1,NASH2,NASH2,1.D0, 1694 & ZYM(IORB1+1,IORB2+NISH2+1),NORBT, 1695 & DV(IASH2+1,IASH2+1),NASHT,0.D0,WRK(KDXV),NORB1) 1696C ** step 2 of active dm: 1697 CALL DGEMM('N','N',NBAS1,NASH2,NORB1,1.D0,CMO(ICMO1+1), 1698 & NBAS1,WRK(KDXV),NORB1,0.D0,WRK(KDXAO1),NBAS1) 1699 IOFMOV = ICMO2 + 1 + NISH2*NBAS2 1700 IDXAO2 = KDXAO2 + IBAS(ISYM2)*NBAST + IBAS(ISYM1) 1701 CALL DGEMM('N','T',NBAS1,NBAS2,NASH2,1.D0,WRK(KDXAO1),NBAS1, 1702 & CMO(IOFMOV),NBAS2,0.D0,WRK(IDXAO2),NBAST) 1703C ** this symmetry block finished 1704 3000 CONTINUE 1705C 1706C subtract D_II from D-I 1707 CALL DAXPY(N2BASX,DM1,WRK(KDXAO2),1,DXVAO,1) 1708 END IF 1709C 1710C DXVAO finished, now DXCAO, if nisht.gt.0: 1711 IF (NISHT .GT. 0) THEN 1712 CALL DZERO(WRK(KDXAO2),N2BASX) 1713 DO 4000 ISYM1 = 1,ABS_NSYM 1714 NISH1 = NISH(ISYM1) 1715 IF (NISH1 .EQ. 0) GO TO 4000 1716 ISYM2 = MULD2H(ISYM1,ABS_GRADSYM) 1717 NBAS1 = NBAS(ISYM1) 1718 NBAS2 = NBAS(ISYM2) 1719 IORB1 = IORB(ISYM1) 1720 IORB2 = IORB(ISYM2) 1721 NORB1 = NORB(ISYM1) 1722 NORB2 = NORB(ISYM2) 1723 ICMO1 = ICMO(ISYM1) 1724 ICMO2 = ICMO(ISYM2) 1725C 1726 IF (NORB2 .NE. 0) THEN 1727 CALL DGEMM('N','T',NBAS2,NISH1,NORB2,1.D0,CMO(ICMO2+1), 1728 & NBAS2,ZYM(IORB1+1,IORB2+1),NORBT,0.D0, 1729 & WRK(KDXAO1),NBAS2) 1730 IDXAO2 = KDXAO2 + IBAS(ISYM2)*NBAST + IBAS(ISYM1) 1731 CALL DGEMM('N','T',NBAS1,NBAS2,NISH1,1.D0,CMO(ICMO1+1), 1732 & NBAS1,WRK(KDXAO1),NBAS2,0.D0,WRK(IDXAO2),NBAST) 1733 END IF 1734 4000 CONTINUE 1735 CALL DCOPY(N2BASX,WRK(KDXAO2),1,DXCAO,1) 1736c 1737 CALL DZERO(WRK(KDXAO2),N2BASX) 1738 DO 5000 ISYM1 = 1,NSYM 1739 NISH1 = NISH(ISYM1) 1740 IF (NISH1 .EQ. 0) GO TO 5000 1741 ISYM2 = MULD2H(ISYM1,ABS_GRADSYM) 1742 NBAS1 = NBAS(ISYM1) 1743 NBAS2 = NBAS(ISYM2) 1744 IORB1 = IORB(ISYM1) 1745 IORB2 = IORB(ISYM2) 1746 NORB1 = NORB(ISYM1) 1747 NORB2 = NORB(ISYM2) 1748 ICMO1 = ICMO(ISYM1) 1749 ICMO2 = ICMO(ISYM2) 1750 IF (NORB2 .NE. 0) THEN 1751 CALL DGEMM('N','N',NBAS2,NISH1,NORB2,1.D0,CMO(ICMO2+1), 1752 & NBAS2,ZYM(IORB2+1,IORB1+1),NORBT,0.D0, 1753 & WRK(KDXAO1),NBAS2) 1754 IDXAO2 = KDXAO2 + IBAS(ISYM1)*NBAST + IBAS(ISYM2) 1755 CALL DGEMM('N','T',NBAS2,NBAS1,NISH1,1.D0,WRK(KDXAO1),NBAS2, 1756 & CMO(ICMO1+1),NBAS1,0.D0,WRK(IDXAO2),NBAST) 1757 END IF 1758C 1759 5000 CONTINUE 1760C 1761 CALL DAXPY(N2BASX,-1.0d0,WRK(KDXAO2),1,DXCAO,1) 1762 ENDIF 1763C 1764 IF ( IPRABSLRS.GE.200 ) THEN 1765 WRITE (LUPRI,'(/A)') 1766 & 'Final result in ABS_ZYTRA' 1767 CALL OUTPUT(DXCAO,1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 1768 END IF 1769C 1770 CALL QEXIT('ABS_ZYTRA') 1771C 1772 RETURN 1773 END 1774 SUBROUTINE ABSLRSRESULT(MJWOP,RESLRF,WRK,LWRK) 1775C 1776#include "implicit.h" 1777#include "priunit.h" 1778C 1779C Used from 1780C inforb.h: NSYM 1781C codata.h: XTEV 1782C qrinf.h: MZYVAR (kzyvar=mzyvar(isym)) 1783C infvar.h: MAXWOP 1784C infdim.h: NVARMA (when ABS_ANALYZE) 1785#include "abslrs.h" 1786#include "inforb.h" 1787#include "codata.h" 1788#include "qrinf.h" 1789#include "infvar.h" 1790#include "infdim.h" 1791#include "absorp.h" 1792!Information on nuclei/nuclear moments required for NSCD/NSOR. Sonia 1793#include "mxcent.h" 1794#include "maxorb.h" 1795#include "maxaqn.h" 1796#include "chrnos.h" 1797#include "nuclei.h" 1798#include "symmet.h" 1799#include "cbiher.h" 1800#include "spnout.h" 1801 1802C 1803 LOGICAL FOUND, CONV 1804C 1805 DIMENSION MJWOP(2,MAXWOP,8),RESLRF(2,ABS_NFREQ_INTERVAL,3,3,4) 1806 DIMENSION WRK(LWRK),ISRCORBR(5),ISRCORBI(5) 1807 !Magneto-chiral dichroism averages. Janusz 1808 DIMENSION AREAABB(NFREQ_BETA_B),AREABBA(NFREQ_BETA_B) 1809 DIMENSION AIMAABB(NFREQ_BETA_B),AIMABBA(NFREQ_BETA_B) 1810 DIMENSION GREABG(NFREQ_BETA_B),GIMABG(NFREQ_BETA_B) 1811 DIMENSION BIREFR(NFREQ_BETA_B),DICHRO(NFREQ_BETA_B) 1812 !Nuclear spin CD/OR averages. Sonia 1813 CHARACTER*8 LABINT, LABX, LABY, LABZ, BLANK 1814 PARAMETER (BLANK=' ') 1815 DOUBLE PRECISION XNSCD(NPATOM), XNSOR(NPATOM) 1816C 1817 KFREE = 1 1818 LFREE = LWRK 1819C 1820 CALL TITLER('FINAL RESULTS FOR ABSORPTION','*',120) 1821C 1822 CALL AROUND('Polarizability with damping') 1823 CALL PRSYMB(LUPRI,'-',66,1) 1824 WRITE(LUPRI,'(A3,2A10,A12,2A16)') 1825 & ' No','A-oper','B-oper','Frequency','Real part','Imag part' 1826 CALL PRSYMB(LUPRI,'-',66,1) 1827 DO IFREQ=1,ABS_NFREQ_INTERVAL 1828 IF (ABSLRS_INTERVAL) THEN 1829 FREQ = ABS_FREQ_INTERVAL(1) + (IFREQ-1)*ABS_FREQ_INTERVAL(3) 1830 ELSE 1831 FREQ = ABS_FREQ_ALPHA(IFREQ) 1832 END IF 1833 IF (.NOT. ABS_VELOCI) THEN 1834 WRITE(LUPRI,'(3(/,I3,2A10,F12.6,2F16.6))') 1835 & 6*IFREQ-5,'XDIPLEN','XDIPLEN',FREQ, 1836 & RESLRF(1,IFREQ,1,1,1),RESLRF(2,IFREQ,1,1,1), 1837 & 6*IFREQ-4,'YDIPLEN','YDIPLEN',FREQ, 1838 & RESLRF(1,IFREQ,2,2,1),RESLRF(2,IFREQ,2,2,1), 1839 & 6*IFREQ-3,'ZDIPLEN','ZDIPLEN',FREQ, 1840 & RESLRF(1,IFREQ,3,3,1),RESLRF(2,IFREQ,3,3,1) 1841 WRITE(LUPRI,'(3(I3,2A10,F12.6,2F16.6,/))') 1842 & 6*IFREQ-2,'XDIPLEN','YDIPLEN',FREQ, 1843 & RESLRF(1,IFREQ,1,2,1),RESLRF(2,IFREQ,1,2,1), 1844 & 6*IFREQ-1,'XDIPLEN','ZDIPLEN',FREQ, 1845 & RESLRF(1,IFREQ,1,3,1),RESLRF(2,IFREQ,1,3,1), 1846 & 6*IFREQ,'YDIPLEN','ZDIPLEN',FREQ, 1847 & RESLRF(1,IFREQ,2,3,1),RESLRF(2,IFREQ,2,3,1) 1848 WRITE(LUPRI,'(6X,A14,3X,F12.6,2F16.6)') 1849 & 'Averaged value',FREQ, 1850 & (RESLRF(1,IFREQ,1,1,1)+RESLRF(1,IFREQ,2,2,1)+ 1851 & RESLRF(1,IFREQ,3,3,1))/3, 1852 & (RESLRF(2,IFREQ,1,1,1)+RESLRF(2,IFREQ,2,2,1)+ 1853 & RESLRF(2,IFREQ,3,3,1))/3 1854 ELSE IF (ABS_VELOCI) THEN 1855 WRITE(LUPRI,'(3(/,I3,2A10,F12.6,2F16.6))') 1856 & 6*IFREQ-5,'XDIPVEL','XDIPVEL',FREQ, 1857 & RESLRF(1,IFREQ,1,1,4),RESLRF(2,IFREQ,1,1,4), 1858 & 6*IFREQ-4,'YDIPVEL','YDIPVEL',FREQ, 1859 & RESLRF(1,IFREQ,2,2,4),RESLRF(2,IFREQ,2,2,4), 1860 & 6*IFREQ-3,'ZDIPVEL','ZDIPVEL',FREQ, 1861 & RESLRF(1,IFREQ,3,3,4),RESLRF(2,IFREQ,3,3,4) 1862 WRITE(LUPRI,'(3(I3,2A10,F12.6,2F16.6,/))') 1863 & 6*IFREQ-2,'XDIPVEL','YDIPVEL',FREQ, 1864 & RESLRF(1,IFREQ,1,2,4),RESLRF(2,IFREQ,1,2,4), 1865 & 6*IFREQ-1,'XDIPVEL','ZDIPVEL',FREQ, 1866 & RESLRF(1,IFREQ,1,3,4),RESLRF(2,IFREQ,1,3,4), 1867 & 6*IFREQ,'YDIPVEL','ZDIPVEL',FREQ, 1868 & RESLRF(1,IFREQ,2,3,4),RESLRF(2,IFREQ,2,3,4) 1869 WRITE(LUPRI,'(6X,A14,3X,F12.6,2F16.6)') 1870 & 'Averaged value',FREQ, 1871 & (RESLRF(1,IFREQ,1,1,4)+RESLRF(1,IFREQ,2,2,4)+ 1872 & RESLRF(1,IFREQ,3,3,4))/3, 1873 & (RESLRF(2,IFREQ,1,1,4)+RESLRF(2,IFREQ,2,2,4)+ 1874 & RESLRF(2,IFREQ,3,3,4))/3 1875 END IF 1876 END DO 1877 WRITE(LUPRI,'(A)')' ' 1878 CALL PRSYMB(LUPRI,'-',66,1) 1879 WRITE (LUPRI,'(A,F10.6,A,F7.4,A,F7.1,A)') 1880 & ' Damping parameter equals', 1881 & ABS_DAMP,' au =', 1882 & ABS_DAMP*XTEV,' eV =', 1883 & ABS_DAMP*XTKAYS,' cm-1' 1884c 1885 IF (ABSLRS_BETA) THEN 1886 CALL AROUND('First-order hyperpolarizability'// 1887 & ' with damping') 1888 ELSE IF (ABSLRS_MCD) THEN 1889 CALL AROUND('Magnetic circular dichroism'// 1890 & ' with damping') 1891 END IF 1892c 1893!-------------------------------------------------------------------! 1894! Magneto Optical Activity (CD and OR) ! 1895!-------------------------------------------------------------------! 1896c 1897 IF (ABSLRS_BETA .OR. ABSLRS_MCD) THEN 1898 CALL PRSYMB(LUPRI,'-',94,1) 1899 WRITE(LUPRI,'(A3,6A10,2A16)') 1900 & ' No','A-oper','B-oper','C-oper','A-freq', 1901 & 'B-freq','C-freq','Real part','Imag part' 1902 CALL PRSYMB(LUPRI,'-',94,1) 1903 BTMP = 99.9D9 1904 CTMP = 99.9D9 1905 BTERM = 0.0 1906 RMORD = 0.0 1907 DO IQRF=1,NQRF 1908 IF (BTMP.NE.QRFFRQ(IQRF,1).OR.CTMP.NE.QRFFRQ(IQRF,2)) THEN 1909 IF ((ABSLRS_MCD).AND.(.NOT.IQRF.EQ.1)) THEN 1910 WRITE(LUPRI,'(5X,A21,7X,3F10.6,2F16.6)') 1911 & 'Orientational average',(QRFFRQ(IQRF-1,I),I=1,3), 1912 & BTERM,RMORD 1913 BTERM = 0.0 1914 RMORD = 0.0 1915 END IF 1916 WRITE(LUPRI,'(A)') ' ' 1917 BTMP = QRFFRQ(IQRF,1) 1918 CTMP = QRFFRQ(IQRF,2) 1919 END IF 1920 IF (ABSLRS_BETA) THEN 1921 WRITE(LUPRI,'(I5,3A10,3F10.6,2F16.6)') IQRF, 1922 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3), 1923 & (RES_BETA(IQRF,I), I=1,2) 1924 ELSE IF (ABSLRS_MCD) THEN 1925 WRITE(LUPRI,'(I3,3A10,3F10.6,2F16.6)') IQRF, 1926 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3), 1927 & RES_BETA(IQRF,2),RES_BETA(IQRF,1) 1928C 1929C Add contribution to orientational average 1930C 1931 IF ((QRFLAB(IQRF,1)(:1).EQ.'X') 1932 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Y') 1933 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 1934 BTERM = BTERM - RES_BETA(IQRF,2) 1935 RMORD = RMORD - RES_BETA(IQRF,1) 1936 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 1937 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Z') 1938 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 1939 BTERM = BTERM - RES_BETA(IQRF,2) 1940 RMORD = RMORD - RES_BETA(IQRF,1) 1941 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 1942 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'X') 1943 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 1944 BTERM = BTERM - RES_BETA(IQRF,2) 1945 RMORD = RMORD - RES_BETA(IQRF,1) 1946 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 1947 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Y') 1948 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 1949 BTERM = BTERM + RES_BETA(IQRF,2) 1950 RMORD = RMORD + RES_BETA(IQRF,1) 1951 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'X') 1952 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Z') 1953 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 1954 BTERM = BTERM + RES_BETA(IQRF,2) 1955 RMORD = RMORD + RES_BETA(IQRF,1) 1956 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 1957 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'X') 1958 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 1959 BTERM = BTERM + RES_BETA(IQRF,2) 1960 RMORD = RMORD + RES_BETA(IQRF,1) 1961 END IF 1962 END IF 1963 END DO 1964 WRITE(LUPRI,'(5X,A21,7X,3F10.6,2F16.6)') 1965 & 'Orientational average',(QRFFRQ(NQRF,I), I=1,3), 1966 & BTERM,RMORD 1967 WRITE(LUPRI,'(A)')' ' 1968 CALL PRSYMB(LUPRI,'-',94,1) 1969 WRITE (LUPRI,'(A,F10.6,A,F7.4,A,F7.1,A)') 1970 & ' Damping parameter equals', 1971 & DAMPING,' au =', 1972 & DAMPING*XTEV,' eV =', 1973 & DAMPING*XTKAYS,' cm-1' 1974 END IF 1975!-------------------------------------------------------------------! 1976! Nuclear Spin Optical Activity (CD and OR) - S. Coriani, 2013 ! 1977!-------------------------------------------------------------------! 1978 IF (ABSLRS_NSCD) THEN 1979 1980 write(lupri,*)'NPATOM=', NPATOM 1981 write(lupri,*) (IPATOM(I),I = 1, NPATOM) 1982 1983 CALL AROUND('Nuclear Spin optical rotation/circular dichroism'// 1984 & ' with damping') 1985 CALL PRSYMB(LUPRI,'-',95,1) 1986 WRITE(LUPRI,'(A3,6A10,2A16)') 1987 & ' No','A-oper','B-oper','C-oper','A-freq', 1988 & 'B-freq','C-freq','Real part','Imag part' 1989 CALL PRSYMB(LUPRI,'-',95,1) 1990 BTMP = 99.9D9 1991 CTMP = 99.9D9 1992 !NPATOM is the number of nuclei requested in input for which 1993 !PSO has been computed 1994 CALL DZERO(XNSCD,NPATOM) 1995 CALL DZERO(XNSOR,NPATOM) 1996 !Sonia: note that RES_BETA has opposite sign than std 1997 !code since scaled by -1. I change the sign again 1998 !here to adapt to old output 1999 DO IQRF=1,NQRF 2000 WRITE(LUPRI,'(I3,3A10,3F10.6,2F16.6)') IQRF, 2001 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3), 2002 & -RES_BETA(IQRF,2),-RES_BETA(IQRF,1) 2003 END DO 2004 2005 CALL PRSYMB(LUPRI,'=',97,1) 2006 WRITE(LUPRI,'(A10,2X,3A10,2X,2A14,2X,2A17)') 2007 & ' Averages ',' A-freq ',' B-freq ',' C-freq ', 2008 & ' 6*Re<<A,B,C>> ',' 6*Im<<A,B,C>> ',' CD(mu rad/M cm) ', 2009 & ' OR(mu rad/M cm) ' 2010 CALL PRSYMB(LUPRI,'=',97,1) 2011 2012 DO IDX=1,NPATOM 2013 IWHICH=IPATOM(IDX) 2014 !quick fix to resume the info on Gval of atom 2015 do ii=1,6 2016 xmmom = DISOTP(IZATOM(IWHICH),ii,'MMOM') 2017 xspin = DISOTP(IZATOM(IWHICH),ii,'SPIN') 2018 if (ABS(xmmom).gt.1.0D-4) exit 2019 end do 2020 Unit_Factor = xmmom*(1.8687D0) !(micro rad / (M cm)) 2021 !--------------------------------------------------- 2022 2023 DO IQRF=1,NQRF 2024 IF (BTMP.NE.QRFFRQ(IQRF,1).OR.CTMP.NE.QRFFRQ(IQRF,2)) THEN 2025 IF (.NOT.IQRF.EQ.1) THEN 2026 WRITE(LUPRI,'(1X,A5,A2,A1,1X,3F10.6,4F16.6)') 2027 & 'Atom(',NAMN(IWHICH),')', 2028 & (QRFFRQ(IQRF-1,I),I=1,3), 2029 & XNSCD(IDX),XNSOR(IDX), 2030 & XNSCD(IDX)*Unit_Factor*QRFFRQ(IQRF-1,2), 2031 & XNSOR(IDX)*Unit_Factor*QRFFRQ(IQRF-1,2) 2032 XNSCD(IDX) = 0.0D0 2033 XNSOR(IDX) = 0.0D0 2034 END IF 2035 WRITE(LUPRI,'(A)') ' ' 2036 BTMP = QRFFRQ(IQRF,1) 2037 CTMP = QRFFRQ(IQRF,2) 2038 END IF 2039C Add contribution to orientational average 2040 !write (lupri,*)'Name of nucleus is ', NAMN(IDX) 2041 !WARNING: FOLLOWING ONLY TESTED FOR NO SYMMETRY CASE 2042 KSYMOP = 1 2043 !ISCORX = IPTCNT(3*(IDX-1)+1,KSYMOP-1,2) 2044 !ISCORY = IPTCNT(3*(IDX-1)+2,KSYMOP-1,2) 2045 !ISCORZ = IPTCNT(3*(IDX-1)+3,KSYMOP-1,2) 2046 ISCORX = IPTCNT(3*(IWHICH-1)+1,KSYMOP-1,2) 2047 ISCORY = IPTCNT(3*(IWHICH-1)+2,KSYMOP-1,2) 2048 ISCORZ = IPTCNT(3*(IWHICH-1)+3,KSYMOP-1,2) 2049 IF (ISCORX .GT. 0) THEN 2050 IFIRST = ISCORX/100 2051 ISECND = MOD(ISCORX,100)/10 2052 ITHIRD = MOD(MOD(ISCORX,100),10) 2053 LABX = 'PSO '//CHRNOS(IFIRST)// 2054 & CHRNOS(ISECND)// 2055 & CHRNOS(ITHIRD)//' ' 2056! write(lupri,*)'Test label PSO X:', LABX 2057 END IF 2058 IF (ISCORY .GT. 0) THEN 2059 IFIRST = ISCORY/100 2060 ISECND = MOD(ISCORY,100)/10 2061 ITHIRD = MOD(MOD(ISCORY,100),10) 2062 LABY = 'PSO '//CHRNOS(IFIRST)// 2063 & CHRNOS(ISECND)// 2064 & CHRNOS(ITHIRD)//' ' 2065 !write(lupri,*)'Test label PSO Y:', LABY 2066 END IF 2067 IF (ISCORZ .GT. 0) THEN 2068 IFIRST = ISCORZ/100 2069 ISECND = MOD(ISCORZ,100)/10 2070 ITHIRD = MOD(MOD(ISCORZ,100),10) 2071 LABZ = 'PSO '//CHRNOS(IFIRST)// 2072 & CHRNOS(ISECND)// 2073 & CHRNOS(ITHIRD)//' ' 2074 !write(lupri,*)'Test label PSO Z:', LABZ 2075 END IF 2076 2077! write(lupri,*)' QRFLAB(IQRF,3)(1:8)', 2078! & QRFLAB(IQRF,3)(1:8) 2079 2080 IF ((QRFLAB(IQRF,1)(:1).EQ.'X') 2081 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Y') 2082 & .AND.(QRFLAB(IQRF,3)(1:8).EQ.LABZ)) THEN 2083 2084! write(lupri,*)' XYZ ' 2085 XNSCD(IDX) = XNSCD(IDX) - RES_BETA(IQRF,2) 2086 XNSOR(IDX) = XNSOR(IDX) - RES_BETA(IQRF,1) 2087 2088 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2089 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Z') 2090 & .AND.(QRFLAB(IQRF,3)(1:8).EQ.LABX)) THEN 2091! write(lupri,*)' YZX ' 2092 XNSCD(IDX) = XNSCD(IDX) - RES_BETA(IQRF,2) 2093 XNSOR(IDX) = XNSOR(IDX) - RES_BETA(IQRF,1) 2094 2095 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2096 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'X') 2097 & .AND.(QRFLAB(IQRF,3)(1:8).EQ.LABY)) THEN 2098! write(lupri,*)' ZXY ' 2099 XNSCD(IDX) = XNSCD(IDX) - RES_BETA(IQRF,2) 2100 XNSOR(IDX) = XNSOR(IDX) - RES_BETA(IQRF,1) 2101 2102 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2103 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Y') 2104 & .AND.(QRFLAB(IQRF,3)(1:8).EQ.LABX)) THEN 2105! write(lupri,*)' ZYX ' 2106 XNSCD(IDX) = XNSCD(IDX) + RES_BETA(IQRF,2) 2107 XNSOR(IDX) = XNSOR(IDX) + RES_BETA(IQRF,1) 2108 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'X') 2109 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'Z') 2110 & .AND.(QRFLAB(IQRF,3)(1:8).EQ.LABY)) THEN 2111! write(lupri,*)' XZY ' 2112 XNSCD(IDX) = XNSCD(IDX) + RES_BETA(IQRF,2) 2113 XNSOR(IDX) = XNSOR(IDX) + RES_BETA(IQRF,1) 2114 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2115 & .AND.(QRFLAB(IQRF,2)(:1).EQ.'X') 2116 & .AND.(QRFLAB(IQRF,3)(1:8).EQ.LABZ)) THEN 2117! write(lupri,*)' YXZ ' 2118 XNSCD(IDX) = XNSCD(IDX) + RES_BETA(IQRF,2) 2119 XNSOR(IDX) = XNSOR(IDX) + RES_BETA(IQRF,1) 2120 END IF 2121 END DO 2122 WRITE(LUPRI,'(1X,A5,A2,A1,1X,3F10.6,4F16.6)') 2123 & 'Atom(',NAMN(IWHICH),')', 2124 & (QRFFRQ(NQRF,I), I=1,3), 2125 & XNSCD(IDX),XNSOR(IDX), 2126 & XNSCD(IDX)*Unit_Factor*QRFFRQ(NQRF,2), 2127 & XNSOR(IDX)*Unit_Factor*QRFFRQ(NQRF,2) 2128 END DO 2129 WRITE(LUPRI,'(A)')' ' 2130 CALL PRSYMB(LUPRI,'-',94,1) 2131 WRITE (LUPRI,'(A,F10.6,A,F7.4,A,F7.1,A)') 2132 & ' Damping parameter equals', 2133 & DAMPING,' au =', 2134 & DAMPING*XTEV,' eV =', 2135 & DAMPING*XTKAYS,' cm-1' 2136 END IF 2137c ------------------------------------------------------------------- 2138c Magneto-Chiral Dichroism and Birefringence (J. Cukras, 2013) 2139c ------------------------------------------------------------------- 2140 2141 IF ( ABSLRS_MCHD ) THEN 2142 LUFU=-1 2143 CALL GPOPEN(LUFU,'MCHD_QRF.list','NEW',' ','FORMATTED', 2144 & IDUMMY,.FALSE.) 2145 CALL AROUND('Magneto-chiral circular dichroism'// 2146 & ' with damping') 2147 CALL PRSYMB(LUFU,'-',94,1) 2148 WRITE(LUFU,'(A3,6A10,2A16)') 2149 & ' No','A-oper','B-oper','C-oper','A-freq', 2150 & 'B-freq','C-freq','G Real part','G Imag part' 2151 CALL PRSYMB(LUFU,'-',94,1) 2152 BTMP = 99.9D9 2153 CTMP = 99.9D9 2154 GTERMRE = 0.0 2155 GTERMIM = 0.0 2156 A1TERMRE = 0.0 2157 A1TERMIM = 0.0 2158 A2TERMRE = 0.0 2159 A2TERMIM = 0.0 2160 IFREQ=0 2161C Loop for G 2162 DO IQRF=1,NQRF 2163 IF (BTMP.NE.QRFFRQ(IQRF,1).OR.CTMP.NE.QRFFRQ(IQRF,2)) THEN 2164 IF (.NOT.IQRF.EQ.1) THEN 2165 WRITE(LUFU,'(5X,A21,7X,3F10.6,2F16.6)') 2166 & 'G_abg orientational average',(QRFFRQ(IQRF-1,I),I=1,3), 2167 & GTERMRE,GTERMIM 2168 IFREQ=IFREQ+1 2169 GREABG(IFREQ) = GTERMRE 2170 GIMABG(IFREQ) = GTERMIM 2171 GTERMRE = 0.0 2172 GTERMIM = 0.0 2173 END IF 2174 WRITE(LUFU,'(A)') ' ' 2175 BTMP = QRFFRQ(IQRF,1) 2176 CTMP = QRFFRQ(IQRF,2) 2177 END IF 2178 WRITE(LUFU,'(I5,3A10,3F10.6,2F16.6)') IQRF, 2179 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3), 2180 & RES_BETA(IQRF,1),RES_BETA(IQRF,2) 2181 2182C Contribution to orientational avarage for 2183C G_{\alpha,\beta,\gamma} 2184 2185 IF ((QRFLAB(IQRF,1)(:1).EQ.'X') 2186 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'YA') 2187 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 2188 GTERMRE = GTERMRE - RES_BETA(IQRF,1) 2189 GTERMIM = GTERMIM - RES_BETA(IQRF,2) 2190 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2191 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'ZA') 2192 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 2193 GTERMRE = GTERMRE - RES_BETA(IQRF,1) 2194 GTERMIM = GTERMIM - RES_BETA(IQRF,2) 2195 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2196 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XA') 2197 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 2198 GTERMRE = GTERMRE - RES_BETA(IQRF,1) 2199 GTERMIM = GTERMIM - RES_BETA(IQRF,2) 2200 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2201 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'YA') 2202 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 2203 GTERMRE = GTERMRE + RES_BETA(IQRF,1) 2204 GTERMIM = GTERMIM + RES_BETA(IQRF,2) 2205 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'X') 2206 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'ZA') 2207 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 2208 GTERMRE = GTERMRE + RES_BETA(IQRF,1) 2209 GTERMIM = GTERMIM + RES_BETA(IQRF,2) 2210 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2211 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XA') 2212 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 2213 GTERMRE = GTERMRE + RES_BETA(IQRF,1) 2214 GTERMIM = GTERMIM + RES_BETA(IQRF,2) 2215 END IF 2216 2217 END DO 2218 WRITE(LUFU,'(5X,A21,7X,3F10.6,2F16.6)') 2219 & 'G_abg orientational average',(QRFFRQ(NQRF,I), I=1,3), 2220 & GTERMRE,GTERMIM 2221 IFREQ=IFREQ+1 2222 GREABG(IFREQ) = GTERMRE 2223 GIMABG(IFREQ) = GTERMIM 2224 2225! CALL PRSYMB(LUPRI,'-',94,1) 2226! WRITE(LUPRI,'(A3,6A10,2A16)') 2227! & ' No','A-oper','B-oper','C-oper','A-freq', 2228! & 'B-freq','C-freq','A Real part','A Imag part' 2229! CALL PRSYMB(LUPRI,'-',94,1) 2230 BTMP = 99.9D9 2231 CTMP = 99.9D9 2232 IFREQ=0 2233C Loop for A 2234 DO IQRF=1,NQRF 2235 IF (BTMP.NE.QRFFRQ(IQRF,1).OR.CTMP.NE.QRFFRQ(IQRF,2)) THEN 2236 IF (.NOT.IQRF.EQ.1) THEN 2237 WRITE(LUFU,'(5X,A21,7X,3F10.6,2F16.6)') 2238 & 'A_aabb orientational average',(QRFFRQ(IQRF-1,I),I=1,3), 2239 & A1TERMRE,A1TERMIM 2240 WRITE(LUFU,'(5X,A21,7X,3F10.6,2F16.6)') 2241 & 'A_abba orientational average',(QRFFRQ(IQRF-1,I),I=1,3), 2242 & A2TERMRE,A2TERMIM 2243 IFREQ=IFREQ+1 2244 AREAABB(IFREQ) = A1TERMRE 2245 AIMAABB(IFREQ) = A1TERMIM 2246 AREABBA(IFREQ) = A2TERMRE 2247 AIMABBA(IFREQ) = A2TERMIM 2248 A1TERMRE = 0.0 2249 A1TERMIM = 0.0 2250 A2TERMRE = 0.0 2251 A2TERMIM = 0.0 2252 END IF 2253c WRITE(LUPRI,'(A)') ' ' 2254 BTMP = QRFFRQ(IQRF,1) 2255 CTMP = QRFFRQ(IQRF,2) 2256 END IF 2257 WRITE(LUFU,'(I5,3A10,3F10.6,2F16.6)') IQRF, 2258 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3), 2259 & RES_BETA(IQRF,2),RES_BETA(IQRF,1) 2260c Orientational avarage form 2261c A_{\alpha,\alpha\beta,\beta} 2262 IF ((QRFLAB(IQRF,1)(:1).EQ.'X') 2263 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XX') 2264 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 2265 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2266 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2267C Part of the abba 2268 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2269 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2270 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2271 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'YY') 2272 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 2273 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2274 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2275C Part of the abba 2276 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2277 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2278 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2279 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'ZZ') 2280 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 2281 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2282 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2283C Part of the abba 2284 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2285 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2286 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'X') 2287 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XY') 2288 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 2289 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2290 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2291 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2292 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'YZ') 2293 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 2294 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2295 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2296 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'X') 2297 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XZ') 2298 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 2299 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2300 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2301 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2302 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XY') 2303 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 2304 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2305 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2306 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2307 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'YZ') 2308 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 2309 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2310 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2311 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2312 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XZ') 2313 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 2314 A1TERMRE = A1TERMRE - RES_BETA(IQRF,2) 2315 A1TERMIM = A1TERMIM - RES_BETA(IQRF,1) 2316C Orientational avarage form 2317C A_{\alpha,\beta\beta,\alpha} 2318C ELSE IF ((QRFLAB(IQRF,1)(:1).EQ.'X') 2319C & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XX') 2320C & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 2321C A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2322C A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2323C ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2324C & .AND.(QRFLAB(IQRF,2)(:2).EQ.'YY') 2325C & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 2326C A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2327C A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2328C ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2329C & .AND.(QRFLAB(IQRF,2)(:2).EQ.'ZZ') 2330C & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 2331C A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2332C A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2333 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'X') 2334 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'YY') 2335 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 2336 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2337 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2338 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2339 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'ZZ') 2340 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 2341 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2342 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2343 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'X') 2344 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'ZZ') 2345 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'X')) THEN 2346 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2347 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2348 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Y') 2349 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XX') 2350 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Y')) THEN 2351 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2352 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2353 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2354 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'YY') 2355 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 2356 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2357 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2358 ELSE IF((QRFLAB(IQRF,1)(:1).EQ.'Z') 2359 & .AND.(QRFLAB(IQRF,2)(:2).EQ.'XX') 2360 & .AND.(QRFLAB(IQRF,3)(:1).EQ.'Z')) THEN 2361 A2TERMRE = A2TERMRE - RES_BETA(IQRF,2) 2362 A2TERMIM = A2TERMIM - RES_BETA(IQRF,1) 2363 END IF 2364 END DO 2365 2366 WRITE(LUFU,'(5X,A21,7X,3F10.6,2F16.6)') 2367 & 'A_aabb orientational average',(QRFFRQ(NQRF,I),I=1,3), 2368 & A1TERMRE,A1TERMIM 2369 WRITE(LUFU,'(5X,A21,7X,3F10.6,2F16.6)') 2370 & 'A_abba orientational average',(QRFFRQ(NQRF,I),I=1,3), 2371 & A2TERMRE,A2TERMIM 2372 IFREQ=IFREQ+1 2373 AREAABB(IFREQ) = A1TERMRE 2374 AIMAABB(IFREQ) = A1TERMIM 2375 AREABBA(IFREQ) = A2TERMRE 2376 AIMABBA(IFREQ) = A2TERMIM 2377 2378 WRITE(LUFU,'(A)')' ' 2379 CALL PRSYMB(LUFU,'-',94,1) 2380 WRITE (LUPRI,'(A,F10.6,A,F7.4,A,F7.1,A)') 2381 & ' Damping parameter equals', 2382 & DAMPING,' au =', 2383 & DAMPING*XTEV,' eV =', 2384 & DAMPING*XTKAYS,' cm-1' 2385 WRITE(LUPRI,'(T5,A)') '','NOTE:', 2386 &'The real and imaginary part of A_abba should be zero.', 2387 &'The quadratic response functions are collected in a separate' 2388 &//' file called MCHD_QRF.list.', 2389 &'You can copy it from the tmp dir using the -get flag of the' 2390 &//' dalton run script.', 2391 &'The precision of the property depends linearly on the' 2392 &//' parameter .THCLR set in the input.', 2393 &'All values are given in a.u.' 2394 WRITE(LUPRI,'(A)') '' 2395 CALL PRSYMB(LUPRI,'-',122,0) 2396 WRITE(LUPRI,'(A6,A6,A14,A14,A14,A14,A11,A11,A16,A16)') '', 2397 & 'FREQ', 2398 & 'Re(G_abg)','Im(G_abg)', 2399 & 'Re(A_aabb)','Im(A_aabb)', 2400 & 'Re(A_abba)','Im(A_abba)', 2401 & "Birefringence",'Dichroism' 2402 CALL PRSYMB(LUPRI,'-',122,0) 2403 DO IFREQ=1,NFREQ_BETA_B 2404 BIREFR(IFREQ)=(FREQ_BETA_B(IFREQ)/5.0D0)*AIMAABB(IFREQ)* 2405 & 0.5D0+ 2406 & GREABG(IFREQ)*0.25D0 2407 DICHRO(IFREQ)=(FREQ_BETA_B(IFREQ)/5.0D0)*AREAABB(IFREQ)*0.5D0+ 2408 & GIMABG(IFREQ)*0.25D0 2409 WRITE(LUPRI,'(A6,F6.3,F14.6,F14.6, 2410 & F14.6,F14.6,F11.6,F11.6,F16.8,F16.8)') 2411 & "#MCHD#",FREQ_BETA_B(IFREQ),0.25D0*GREABG(IFREQ), 2412 & 0.25D0*GIMABG(IFREQ),0.5D0*AREAABB(IFREQ),0.5D0*AIMAABB(IFREQ), 2413 & 0.5D0*AREABBA(IFREQ),0.5D0*AIMABBA(IFREQ), 2414 & BIREFR(IFREQ), 2415 & DICHRO(IFREQ) 2416 END DO 2417 CALL GPCLOSE(LUFU,'KEEP') 2418 END IF 2419 2420C ------------------------------------------------------------------ 2421 2422 IF (ABSLRS_ANALYZE) THEN 2423 MZYVMX = 2*NVARMA 2424 CALL MEMGET('REAL',KVEC,2*MZYVMX,WRK,KFREE,LFREE) 2425 CALL AROUND('Analysis of response vectors') 2426 DO ISYM=1,NSYM 2427 IF (ABS_NOPER(ISYM).GT.0) THEN 2428 KZYVAR=MZYVAR(ISYM) 2429 DO IOPER=1,ABS_NOPER(ISYM) 2430 DO IFREQ=1,ABS_NFREQ_INTERVAL 2431 IF (ABSLRS_INTERVAL) THEN 2432 FREQ = ABS_FREQ_INTERVAL(1) + 2433 & (IFREQ-1)*ABS_FREQ_INTERVAL(3) 2434 ELSE 2435 FREQ = ABS_FREQ_ALPHA(IFREQ) 2436 END IF 2437 WRITE(LUPRI,'(/A11,A10,2(/,A11,F10.6))') 'Property :', 2438 & ABS_LABOP(IOPER,ISYM), 2439 & 'Frequency:',ABS(FREQ), 2440 & 'Damping :',ABS_DAMP 2441C CALL READ_XVEC(LUABSVECS,2*KZYVAR,WRK(KVEC), 2442C & ABS_LABOP(IOPER,ISYM),ISYM, 2443C & ABS(FREQ),ABS_THCLR,FOUND,CONV) 2444 CALL READ_XVEC2(LUABSVECS,2*KZYVAR,WRK(KVEC), 2445 & ABS_LABOP(IOPER,ISYM),BLANK,ISYM, 2446 & ABS(FREQ),0.0d0,ABS_THCLR,FOUND,CONV) 2447 IF (.NOT. FOUND) THEN 2448 WRITE(LUPRI,'(/A)') 2449 & ' Response vector not found on file LUABSVECS' 2450 CALL QUIT('Response vector not found on file') 2451 ELSE IF(.NOT. CONV) THEN 2452 WRITE (LUPRI,'(/A)') ' @WARNING: '// 2453 & 'Response vector not converged on file LUABSVECS' 2454 END IF 2455C 2456 DNORM_RE=DNRM2(KZYVAR,WRK(KVEC),1) 2457 DNORM_IM=DNRM2(KZYVAR,WRK(KVEC+KZYVAR),1) 2458C 2459 IF (IPRABSLRS.GT.1) THEN 2460 WRITE(LUPRI,'(/A,2F14.8)') 2461 & ' Norm of vector (real and imag):', 2462 & DNORM_RE, DNORM_IM 2463C 2464 WRITE(LUPRI,'(/7X,2A5,2(2X,2(5X,A2,5X)))') 2465 & 'occ','vir','ZR','YR','ZI','YI' 2466 DO I=1,KZYVAR/2 2467 WRITE(LUPRI,'(I5,2X,2I5,2(2X,2F12.8))') 2468 & I,MJWOP(1,I,ISYM),MJWOP(2,I,ISYM), 2469 & WRK(KVEC-1+I),WRK(KVEC+KZYVAR/2-1+I), 2470 & WRK(KVEC+KZYVAR-1+I),WRK(KVEC+3*KZYVAR/2-1+I) 2471 END DO 2472C 2473 END IF 2474 2475 WRITE(LUPRI,'(/A,2F14.8)') 2476 & ' Norm of vector (real and imag):', 2477 & DNORM_RE, DNORM_IM 2478C 2479C Sum occupied orbital contributions into aggregates 2480C Sonia 2 Joanna: I changed MAXOCC to MAXOCCU because MAXOCC is 2481C a constant in one of the new common blocks I had to include for NSCD 2482C 2483 MAXOCCU=0 2484 DO I=1,KZYVAR/2 2485 MAXOCCU = MAX(MAXOCCU,MJWOP(1,I,ISYM)) 2486 ENDDO 2487 CALL MEMGET('REAL',KORBWR,MAXOCCU,WRK,KFREE,LFREE) 2488 CALL MEMGET('REAL',KORBWI,MAXOCCU,WRK,KFREE,LFREE) 2489c real part 2490 DO I=1,MAXOCCU 2491 WRK(KORBWR+I-1) = 0.0D0 2492 ENDDO 2493 DO I=1,KZYVAR/2 2494 WRK(KORBWR+MJWOP(1,I,ISYM)-1) = 2495 & WRK(KORBWR+MJWOP(1,I,ISYM)-1) + 2496 & WRK(KVEC-1+I)**2 + WRK(KVEC+KZYVAR/2-1+I)**2 2497 ENDDO 2498c imag part 2499 DO I=1,MAXOCCU 2500 WRK(KORBWI+I-1) = 0.0D0 2501 ENDDO 2502 DO I=1,KZYVAR/2 2503 WRK(KORBWI+MJWOP(1,I,ISYM)-1) = 2504 & WRK(KORBWI+MJWOP(1,I,ISYM)-1) + 2505 & WRK(KVEC+KZYVAR-1+I)**2 + 2506 & WRK(KVEC+3*KZYVAR/2-1+I)**2 2507 ENDDO 2508 NSRC=5 2509 CALL FINDMAXN(WRK(KORBWR),MAXOCCU,ISRCORBR,NSRC) 2510 NSRC=5 2511 CALL FINDMAXN(WRK(KORBWI),MAXOCCU,ISRCORBI,NSRC) 2512 WRITE (LUPRI,*) 2513 & 'Important occupied orbital contributions (normalized)' 2514 WRITE (LUPRI,*) 2515 & ' occ real occ imag' 2516 DO I=1,NSRC 2517 WRITE (LUPRI,'(I5,F7.3,I6,F7.3)') ISRCORBR(I), 2518 & WRK(KORBWR+ISRCORBR(I)-1)/DNORM_RE**2, 2519 & ISRCORBI(I), 2520 & WRK(KORBWI+ISRCORBI(I)-1)/DNORM_IM**2 2521 ENDDO 2522 END DO 2523 END DO 2524 END IF 2525 END DO 2526 END IF 2527C 2528 RETURN 2529 END 2530 SUBROUTINE GET_LRF(LU,LABEL,NFREQ_ABS,FREQ_ABS, 2531 & FC,FV,CMO,UDV,PV,XINDX,RESLRF,KZVAR, 2532 & WRK,LWRK) 2533#include "implicit.h" 2534#include "priunit.h" 2535#include "abslrs.h" 2536C 2537 CHARACTER*8 LABEL,LABGD, BLANK 2538 PARAMETER (BLANK=' ') 2539 LOGICAL FOUND,CONV 2540 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),XINDX(*), 2541 & WRK(LWRK),RESLRF(2,ABS_NFREQ_INTERVAL,3,3,4), 2542 & FREQ_ABS(NMXFREQ) 2543 parameter (d0=-1.D0-16) 2544C 2545 KFREE = 1 2546 LFREE = LWRK 2547C 2548 IF (LABEL(2:8).EQ.'DIPLEN') THEN 2549 ITYPE = 1 2550 ELSE IF (LABEL(2:8).EQ.'ANGMOM') THEN 2551 ITYPE = 2 2552 ELSE IF (LABEL(3:7).EQ.'THETA') THEN 2553 ITYPE = 3 2554 ELSE IF (LABEL(2:8).EQ.'DIPVEL') THEN 2555 ITYPE = 4 2556 ELSE 2557 !WRITE(LUPRI,'(A)') ' Warning: Unknown operator!' 2558 ! - could e.g. be PSO for .NSCD 2559 RETURN 2560 END IF 2561 CALL DIPLAB(LABEL,INDEX) 2562C 2563 2564 KZYVAR=2*KZVAR 2565 CALL MEMGET('REAL',KVEC,2*KZYVAR,WRK,KFREE,LFREE) 2566 CALL MEMGET('REAL',KGD,2*KZYVAR,WRK,KFREE,LFREE) 2567 2568 IF (.NOT.(ABS_GDCOMPLEX.OR.ABS_NGD)) THEN 2569C 2570 DO IOPER=1,ABS_NOPER(ABS_GRADSYM) 2571 LABGD = ABS_LABOP(IOPER,ABS_GRADSYM) 2572 IF (LABEL(2:8).EQ.LABGD(2:8)) THEN 2573 CALL DIPLAB(LABGD,INXGD) 2574 CALL GETGPV(LABGD,FC,FV,CMO,UDV,PV,XINDX,ANTSYM, 2575 & WRK(KGD),LFREE) 2576 DO IFREQ=1,ABS_NFREQ_INTERVAL 2577C CALL READ_XVEC(LU,2*KZYVAR,WRK(KVEC), 2578C & LABEL,ABS_GRADSYM, 2579C & abs(FREQ_ABS(IFREQ)),ABS_THCLR,FOUND,CONV) 2580c & ABS(FREQ_ABS(IFREQ),ABS_THCLR,FOUND,CONV) 2581 CALL READ_XVEC2(LU,2*KZYVAR,WRK(KVEC), 2582 & LABEL,BLANK,ABS_GRADSYM, 2583 & abs(FREQ_ABS(IFREQ)),0.0d0,ABS_THCLR,FOUND,CONV) 2584 if (freq_abs(ifreq) .lt. d0) then 2585 call abs_vec_swap(2*kzyvar,wrk(kvec)) 2586 endif 2587 RESLRF(1,IFREQ,INXGD,INDEX,ITYPE) = 2588 & DDOT(KZYVAR,WRK(KGD),1,WRK(KVEC),1) 2589 RESLRF(2,IFREQ,INXGD,INDEX,ITYPE) = 2590 & DDOT(KZYVAR,WRK(KGD), 2591 & 1,WRK(KVEC+KZYVAR),1) 2592 END DO 2593 END IF 2594 END DO 2595 ELSE ! Complex and/or frequency-dependent rhs 2596 CALL QUIT('Collection of terms have not yet been implemented') 2597 END IF 2598C 2599 RETURN 2600 END 2601 SUBROUTINE ABS_QRRDVE(ISYMA,ISYMB,ISYMC,ALAB,BLAB,CLAB, 2602 & FREQA,FREQB,FREQC,KZYVA,KZYVB,KZYVC,VECA,VECB,VECC) 2603C 2604C PURPOSE: Read in response vectors for quadratic response. 2605C 2606#include "implicit.h" 2607#include "priunit.h" 2608#include "absorp.h" 2609#include "abslrs.h" 2610#include "inftap.h" 2611#include "rspprp.h" 2612#include "inflr.h" 2613#include "qrinf.h" 2614C 2615 LOGICAL FOUND, CONV 2616 CHARACTER*8 ALAB,BLAB,CLAB,BLANK 2617 PARAMETER (BLANK=' ') 2618 PARAMETER ( D0=0.0D0, DM1=-1.0D0) 2619 DIMENSION VECA(*),VECB(*),VECC(*) 2620C 2621 KZYVA = MZYVAR(ISYMA) 2622 KZYVB = MZYVAR(ISYMB) 2623 KZYVC = MZYVAR(ISYMC) 2624C 2625 IF (IPRABS.GE.2) THEN 2626 WRITE(LUPRI,'(2(/A),2(/A,3(I10)),/A,3A10,/A,3F10.6)') 2627 & ' Variables in QRRDVE', 2628 & ' ==================================================== ', 2629 & ' KZYVA,KZYVB,KZYVC: ',KZYVA,KZYVB,KZYVC, 2630 & ' ISYMA,ISYMB,ISYMC: ',ISYMA,ISYMB,ISYMC, 2631 & ' ALAB,BLAB,CLAB : ',ALAB,BLAB,CLAB, 2632 & ' FREQA,FREQB,FREQC: ',FREQA,FREQB,FREQC 2633 END IF 2634C 2635C Read in Na 2636C 2637C CALL READ_XVEC(LUABSVECS,2*KZYVA,VECA,ALAB,ISYMA, 2638C & ABS(FREQA),ABS_THCLR,FOUND,CONV) 2639 CALL READ_XVEC2(LUABSVECS,2*KZYVA,VECA,ALAB,BLANK,ISYMA, 2640 & ABS(FREQA),0.0d0,ABS_THCLR,FOUND,CONV) 2641 IF (.NOT. (FOUND .AND. CONV)) THEN 2642 IF (.NOT. FOUND) THEN 2643 WRITE (LUPRI,'(/3A,F8.5,A,I3,/A)') ' Response label ',ALAB, 2644 & ' with frequency ',FREQA, ' and symmetry', 2645 & ISYMA,' not found on file LUABSVECS' 2646 CALL QUIT('Response vector not found on file') 2647 ELSE 2648 WRITE (LUPRI,'(/3A,F8.5,/A,I3,A)') ' @WARNING----'// 2649 & ' Response label ',ALAB, 2650 & ' with frequency ',FREQA, ' and symmetry', 2651 & ISYMA,' not converged on file LUABSVECS' 2652 END IF 2653 END IF 2654 IF (FREQA .GT. D0) THEN 2655 CALL DSWAP(KZYVA/2,VECA,1,VECA(1+KZYVA/2),1) 2656 CALL DSWAP(KZYVA/2,VECA(1+KZYVA),1,VECA(1+KZYVA+KZYVA/2),1) 2657 CALL DSCAL(KZYVA,DM1,VECA,1) 2658 END IF 2659C 2660 IF (IPRABS.GE.10) THEN 2661 WRITE(LUPRI,'(/A)') ' Na vector in ABS_READVEC ' 2662 CALL PRSYMB(LUPRI,'=',72,1) 2663 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 2664 CALL OUTPUT(VECA,1,KZYVA/2,1,4,KZYVA/2,4,1,LUPRI) 2665 END IF 2666C 2667C Read in Nb 2668C 2669C CALL READ_XVEC(LUABSVECS,2*KZYVB,VECB,BLAB,ISYMB, 2670C & ABS(FREQB),ABS_THCLR,FOUND,CONV) 2671 CALL READ_XVEC2(LUABSVECS,2*KZYVB,VECB,BLAB,BLANK,ISYMB, 2672 & ABS(FREQB),0.0d0,ABS_THCLR,FOUND,CONV) 2673 IF (.NOT. (FOUND .AND. CONV)) THEN 2674 IF (.NOT. FOUND) THEN 2675 WRITE (LUPRI,'(/3A,F8.5,A,I3,/A)') ' Response label ',BLAB, 2676 & ' with frequency ',FREQB, ' and symmetry', 2677 & ISYMB,' not found on file LUABSVECS' 2678 CALL QUIT('Response vector not found on file') 2679 ELSE 2680 WRITE (LUPRI,'(/3A,F8.5,/A,I3,A)') ' @WARNING----'// 2681 & ' Response label ',BLAB, 2682 & ' with frequency ',FREQB, ' and symmetry', 2683 & ISYMB,' not converged on file LUABSVECS' 2684 END IF 2685 END IF 2686 IF (FREQB .LT. D0) THEN 2687 CALL DSWAP(KZYVB/2,VECB,1,VECB(1+KZYVB/2),1) 2688 CALL DSWAP(KZYVB/2,VECB(1+KZYVB),1,VECB(1+KZYVB+KZYVB/2),1) 2689 CALL DSCAL(KZYVB,DM1,VECB,1) 2690 END IF 2691C 2692 IF (IPRABS.GE.10) THEN 2693 WRITE(LUPRI,'(/A)') ' Nb vector in ABS_READVEC ' 2694 CALL PRSYMB(LUPRI,'=',72,1) 2695 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 2696 CALL OUTPUT(VECB,1,KZYVB/2,1,4,KZYVB/2,4,1,LUPRI) 2697 END IF 2698C 2699C Read in Nc 2700C 2701C CALL READ_XVEC(LUABSVECS,2*KZYVC,VECC,CLAB,ISYMC, 2702C & ABS(FREQC),ABS_THCLR,FOUND,CONV) 2703 CALL READ_XVEC2(LUABSVECS,2*KZYVC,VECC,CLAB,BLANK,ISYMC, 2704 & ABS(FREQC),0.0d0,ABS_THCLR,FOUND,CONV) 2705 IF (.NOT. (FOUND .AND. CONV)) THEN 2706 IF (.NOT. FOUND) THEN 2707 WRITE (LUPRI,'(/3A,F8.5,A,I3,/A)') ' Response label ',CLAB, 2708 & ' with frequency ',FREQC, ' and symmetry', 2709 & ISYMC,' not found on file LUABSVECS' 2710 CALL QUIT('Response vector not found on file') 2711 ELSE 2712 WRITE (LUPRI,'(/3A,F8.5,/A,I3,A)') ' @WARNING----'// 2713 & ' Response label ',CLAB, 2714 & ' with frequency ',FREQC, ' and symmetry', 2715 & ISYMC,' not converged on file LUABSVECS' 2716 END IF 2717 END IF 2718 IF (FREQC .LT. D0) THEN 2719 CALL DSWAP(KZYVC/2,VECC,1,VECC(1+KZYVC/2),1) 2720 CALL DSWAP(KZYVC/2,VECC(1+KZYVC),1,VECC(1+KZYVC+KZYVC/2),1) 2721 CALL DSCAL(KZYVC,DM1,VECC,1) 2722 END IF 2723C 2724 IF (IPRABS.GE.10) THEN 2725 WRITE(LUPRI,'(/A)') ' Nc vector in ABS_READVEC ' 2726 CALL PRSYMB(LUPRI,'=',72,1) 2727 WRITE(LUPRI,'(4X,4A15)') 'ZR','YR','ZI','YI' 2728 CALL OUTPUT(VECC,1,KZYVC/2,1,4,KZYVC/2,4,1,LUPRI) 2729 END IF 2730C 2731 RETURN 2732 END 2733 SUBROUTINE ABS_BETA_SETUP 2734#include "implicit.h" 2735#include "priunit.h" 2736#include "inforb.h" 2737#include "absorp.h" 2738#include "abslrs.h" 2739C 2740 PARAMETER (THRZERO = 1.0D-6) 2741C 2742 LOGICAL DOHYP 2743 CHARACTER*8 ALAB,BLAB,CLAB 2744C 2745 NQRF = 0 2746 NFREQ_ALPHA = 0 2747 NFREQ_BETA_B=ABS_NFREQ_BETA_B 2748 NFREQ_BETA_C=ABS_NFREQ_BETA_C 2749 DO I=1,ABS_NFREQ_BETA_B 2750 FREQ_BETA_B(I)=ABS_FREQ_BETA_B(I) 2751 ENDDO 2752 DO I=1,ABS_NFREQ_BETA_C 2753 FREQ_BETA_C(I)=ABS_FREQ_BETA_C(I) 2754 ENDDO 2755 ABS_ALPHA=ABSLRS_ALPHA 2756 ABS_BETA=ABSLRS_BETA 2757 ABS_MCD=ABSLRS_MCD 2758 DAMPING=ABS_DAMP 2759C 2760 DO ICFR=1,ABS_NFREQ_BETA_C 2761 DO IBFR=1,ABS_NFREQ_BETA_B 2762 FREQB = ABS_FREQ_BETA_B(IBFR) 2763 FREQC = ABS_FREQ_BETA_C(ICFR) 2764 FREQA = -(FREQB+FREQC) 2765 IF (ABSLRS_SHG .AND. .NOT.FREQB.EQ.FREQC) GOTO 100 2766 DO ISYMA=1,NSYM 2767 DO ISYMB=1,NSYM 2768 ISYMC = MULD2H(ISYMA,ISYMB) 2769 IF (ABS_NOPER(ISYMA).GT.0 .AND. ABS_NOPER(ISYMB).GT.0 2770 & .AND. ABS_NOPER(ISYMC).GT.0) THEN 2771 DO IAOP=1,ABS_NOPER(ISYMA) 2772 DO IBOP=1,ABS_NOPER(ISYMB) 2773 DO ICOP=1,ABS_NOPER(ISYMC) 2774 ALAB = ABS_LABOP(IAOP,ISYMA) 2775 BLAB = ABS_LABOP(IBOP,ISYMB) 2776 CLAB = ABS_LABOP(ICOP,ISYMC) 2777 CALL ABS_QRCHK(DOHYP,ALAB,BLAB,CLAB, 2778 & ISYMA,ISYMB,ISYMC,FREQA,FREQB,FREQC) 2779 IF (DOHYP) THEN 2780 NQRF = NQRF +1 2781 QRFLAB(NQRF,1) = ALAB 2782 QRFLAB(NQRF,2) = BLAB 2783 QRFLAB(NQRF,3) = CLAB 2784 QRFSYM(NQRF,1) = ISYMA 2785 QRFSYM(NQRF,2) = ISYMB 2786 QRFSYM(NQRF,3) = ISYMC 2787 QRFFRQ(NQRF,1) = FREQA 2788 QRFFRQ(NQRF,2) = FREQB 2789 QRFFRQ(NQRF,3) = FREQC 2790 END IF 2791 END DO 2792 END DO 2793 END DO 2794 END IF 2795 END DO 2796 END DO 2797 100 CONTINUE 2798 END DO 2799 END DO 2800C 2801 CALL GPDSRT(NFREQ_ALPHA,FREQ_ALPHA,THRZERO) 2802c ABS_NFREQ_ALPHA = NFREQ_ALPHA 2803c DO I=1,ABS_NFREQ_ALPHA 2804c ABS_FREQ_ALPHA(I)=FREQ_ALPHA(I) 2805c ENDDO 2806! IF (IPRABS.GE.0) THEN 2807 CALL AROUND('Setup of Hyperpolarizability Calculation') 2808 WRITE (LUPRI,'(2(/A),I4,A)') 2809 & ' This calculations requires the solution of linear response', 2810 & ' equations for electric dipole operators at',ABS_NFREQ_ALPHA, 2811 & ' frequencies:' 2812 WRITE(LUPRI,'(/A,5(4F12.8,/,9X))') 2813 & ' LR FREQ:',(ABS_FREQ_ALPHA(I),I=1,ABS_NFREQ_ALPHA) 2814 WRITE (LUPRI,'(/A,I5,A)') 2815 & ' and the evaluation of',NQRF,' quadratic response functions:' 2816 WRITE(LUPRI,'(/2A,/A3,6A12,/2A)') 2817 & '--------------------------------------', 2818 & '-------------------------------------', 2819 & ' No','A-oper','B-oper','C-oper', 2820 & 'A-freq','B-freq','C-freq', 2821 & '--------------------------------------', 2822 & '-------------------------------------' 2823 DO IQRF=1,NQRF 2824 WRITE(LUPRI,'(I4,3A12,3F12.8)') IQRF, 2825 & (QRFLAB(IQRF,I), I=1,3),(QRFFRQ(IQRF,I), I=1,3) 2826 END DO 2827 WRITE(LUPRI,'(2A)') 2828 & '--------------------------------------', 2829 & '-------------------------------------' 2830! END IF 2831C 2832C End of BETA_SETUP 2833C 2834 RETURN 2835 END 2836 SUBROUTINE ABS_QRCHK(DOHYP,ALAB,BLAB,CLAB,ISYMA,ISYMB,ISYMC, 2837 & FREQA,FREQB,FREQC) 2838#include "implicit.h" 2839#include "priunit.h" 2840#include "absorp.h" 2841#include "abslrs.h" 2842C 2843 PARAMETER (THRZERO = 1.0D-6) 2844 LOGICAL DOHYP,NEWFRQ 2845 CHARACTER*8 ALAB,BLAB,CLAB,LAB(3) 2846 DIMENSION FREQ(3),ISYM(3) 2847C 2848 DOHYP = .TRUE. 2849C 2850C Operator requirement for magnetic circular dichroism 2851C 2852 IF (ABS_MCD .OR. ABSLRS_MCD) THEN 2853 IF (ALAB(2:8).NE.'DIPLEN' .OR. BLAB(2:8).NE.'DIPLEN' .OR. 2854 & CLAB(2:8).NE.'ANGMOM') DOHYP = .FALSE. 2855 END IF 2856C 2857C Operator requirement for MCHD 2858C 2859 IF (ABSLRS_MCHD) THEN 2860 IF (ALAB(2:8).NE.'DIPLEN' .OR. 2861 & (BLAB(2:8).NE.('ANGMOM') .AND. BLAB(3:8).NE.'THETA') .OR. 2862 & CLAB(2:8).NE.'ANGMOM') DOHYP = .FALSE. 2863 END IF 2864C 2865C Operator requirement for NSCD 2866C 2867 IF (ABSLRS_NSCD) THEN 2868 IF (ALAB(2:8).NE.'DIPLEN' .OR. BLAB(2:8).NE.('DIPLEN') .OR. 2869 & CLAB(1:3).NE.'PSO') DOHYP = .FALSE. 2870 END IF 2871c 2872C Check if equivalent QRF is in the list already. 2873c 2874 IF (DAMPING .EQ. 0.0D0) THEN 2875C Overall permutational symmetry. 2876 JCTR=3 2877 KCTR=1 2878 LCTR=1 2879 ELSE 2880C Only intrinsic permutational symmetry, so A operator has to match 2881C first operator in list of response functions. 2882 JCTR=1 2883 KCTR=2 2884 LCTR=2 2885 END IF 2886 DO IQRF = 1,NQRF 2887 DO J = 1,JCTR 2888 DO K = KCTR,3 2889 IF (K.NE.J) THEN 2890 DO L = LCTR,3 2891 IF (L.NE.K .AND. L.NE.J) THEN 2892C 2893 IF ( ALAB.EQ.QRFLAB(IQRF,J) .AND. 2894 & BLAB.EQ.QRFLAB(IQRF,K) .AND. 2895 & CLAB.EQ.QRFLAB(IQRF,L) .AND. 2896 & ISYMA.EQ.QRFSYM(IQRF,J) .AND. 2897 & ISYMB.EQ.QRFSYM(IQRF,K) .AND. 2898 & ISYMC.EQ.QRFSYM(IQRF,L) .AND. 2899 & ABS( FREQA-QRFFRQ(IQRF,J)).LT.THRZERO .AND. 2900 & ABS( FREQB-QRFFRQ(IQRF,K)).LT.THRZERO .AND. 2901 & ABS( FREQC-QRFFRQ(IQRF,L)).LT.THRZERO ) THEN 2902 DOHYP = .FALSE. 2903 END IF 2904C 2905 END IF 2906 END DO 2907 END IF 2908 END DO 2909 END DO 2910 END DO 2911C 2912 IF (DOHYP) THEN 2913C 2914C Check if this QRF will inflict new LR solver frequencies. 2915C 2916 FREQ(1) = FREQA 2917 FREQ(2) = FREQB 2918 FREQ(3) = FREQC 2919 DO I=1,3 2920 NEWFRQ = .TRUE. 2921 DO IFR=1,ABS_NFREQ_ALPHA 2922 IF (ABS(ABS_FREQ_ALPHA(IFR)-ABS(FREQ(I))).LT.THRZERO) THEN 2923 NEWFRQ = .FALSE. 2924 END IF 2925 END DO 2926 IF (NEWFRQ) THEN 2927 IF (ABS_NFREQ_ALPHA.GE.NMXFREQ) THEN 2928 WRITE(LUPRI,'(2(/A),I4,A,/A)') 2929 & ' The specified calculation requires more than the allowed', 2930 & ' number of frequencies in the LR solver NMXFREQ=',NMXFREQ,'.', 2931 & ' The program will stop.' 2932 CALL QUIT('Too many frequencies in LR solver.') 2933 END IF 2934 ABS_NFREQ_ALPHA = ABS_NFREQ_ALPHA + 1 2935 ABS_FREQ_ALPHA(ABS_NFREQ_ALPHA) = ABS(FREQ(I)) 2936 END IF 2937 END DO 2938C 2939 END IF 2940C 2941 IF (NQRF.GE.MXQRF .AND. DOHYP) THEN 2942 WRITE(LUPRI,'(2(/A),I4,A,/A)') 2943 & ' The specified calculation requires more than the allowed', 2944 & ' number of quadratic response functions MXQRF=',MXQRF,'.', 2945 & ' The program will stop' 2946 CALL QUIT('Too many quadratic response functions specified.') 2947 END IF 2948C 2949 RETURN 2950 END 2951 SUBROUTINE ABS_VEC_SWAP(NDIM,XVEC) 2952C 2953 DIMENSION XVEC(NDIM,4) 2954 REAL YVEC(NDIM,4) 2955 parameter (dm1=-1.0d0) 2956 2957 call dzero(yvec,4*ndim) 2958 call dcopy(ndim,xvec(1,1),1,yvec(1,2),1) 2959 call dcopy(ndim,xvec(1,2),1,yvec(1,1),1) 2960 call dscal(2*ndim,dm1,yvec(1,1),1) 2961 call dcopy(ndim,xvec(1,3),1,yvec(1,4),1) 2962 call dcopy(ndim,xvec(1,4),1,yvec(1,3),1) 2963 call dcopy(4*ndim,yvec,1,xvec,1) 2964 2965 RETURN 2966 END 2967! end of DALTON/rsp/abscomplex.F 2968