1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck lnrinp */ 20 SUBROUTINE LNRINP(WORD) 21C 22C Nov. and Dec. 93 23C Written by K.L.Bak and P.Joergensen using EXCITA as template. 24C Purpose: To enable calculations of frequency dependent second 25C order response properties in ABACUS. In particular polariza- 26C bilities and vibrational raman optical activity (VROA). 27C 28C Oct - Dev. 99: Extended by K.L.Bak for calculation of 29C vibrational g-factors. 30C 31#include "implicit.h" 32#include "priunit.h" 33#include "mxcent.h" 34#include "maxorb.h" 35#include "cbiexc.h" 36#include "gnrinf.h" 37#include "cbilnr.h" 38#include "abainf.h" 39#include "anrinf.h" 40#include "dorps.h" 41#include "nuclei.h" 42#include "absorp.h" 43#include "abslrs.h" 44#include "maxaqn.h" 45#include "symmet.h" 46#include "codata.h" 47#include "pcmlog.h" 48C 49 PARAMETER (NTABLE = 18) 50 LOGICAL NEWDEF, LOCFIN 51 CHARACTER*8 DIPLEN(3),DIPLOC(3), ANGMOM(3), LONMAG(3) 52 CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 53C 54 DATA DIPLEN/'XDIPLEN','YDIPLEN','ZDIPLEN'/ 55 DATA DIPLOC/'XDIPLOC','YDIPLOC','ZDIPLOC'/ 56 DATA ANGMOM/'XANGMOM','YANGMOM','ZANGMOM'/ 57 DATA LONMAG/'XLONMAG','YLONMAG','ZLONMAG'/ 58C 59 DATA TABLE /'.SKIP ', '.WAVELE','.MAX IT','.THRESH', 60 & '.MAXRED', '.MAXPHP','.STATIC','.VIBGRE', 61 & '.OPTORB', '.FREQUE','.BATCH ','.PRINT ', 62 & '.DAMPIN', '.STOP ','.FREQ I','.REDUCE', 63 & '.NOREBD','.OLDCPP'/ 64C 65 NEWDEF = (WORD .EQ. '*ABALNR') 66 ICHANG = 0 67 NFRVLT = 0 68 IF (NEWDEF) THEN 69 LOCFIN = .FALSE. 70 WORD1 = WORD 71 100 CONTINUE 72 READ (LUCMD, '(A7)') WORD 73 CALL UPCASE(WORD) 74 PROMPT = WORD(1:1) 75 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 76 GO TO 100 77 ELSE IF (PROMPT .EQ. '.') THEN 78 ICHANG = ICHANG + 1 79 DO 200 I = 1, NTABLE 80 IF (TABLE(I) .EQ. WORD) THEN 81 GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17, 82 & 18),I 83 END IF 84 200 CONTINUE 85 IF (WORD .EQ. '.OPTION') THEN 86 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 87 GO TO 100 88 END IF 89 WRITE (LUPRI,'(/3A/)') ' Keyword "',WORD, 90 * '" not recognized in LNRINP.' 91 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 92 CALL QUIT('Illegal keyword under *ABALNR.') 93 1 CONTINUE 94 SKIP = .TRUE. 95 ICHANG = ICHANG + 1 96 GO TO 100 97 2 CONTINUE ! .WAVELE 98 READ(LUCMD,*) NWAVEL 99 IF (NWAVEL .LT. 0) THEN 100 CALL QUIT('# wavelengths negative under *ABALNR') 101 ELSE IF (NWAVEL .EQ. 0) THEN 102 GO TO 100 103 ELSE 104 IF (.NOT. LOCFIN) THEN 105 NFRVAL = 0 106 LOCFIN = .TRUE. 107 END IF 108 END IF 109 NFRVLT = NFRVLT + NWAVEL 110 NFTOT = NWAVEL + NFRVAL 111 IF (NFTOT .GT. MXFR) THEN 112 NWAVEL = MXFR - NFRVAL 113 END IF 114 READ(LUCMD,*) (FRVAL(NFRVAL+I),I=1,NWAVEL) 115 DO I = NFRVAL+1,NFRVAL+NWAVEL 116 IF (ABS(FRVAL(I)) .LT. 1.0D-12) THEN 117 WRITE(LUPRI,*) 118 & I,'th .WAVELE wavelength too small:',FRVAL(I) 119 WRITE(LUPRI,*) 'Use frequency input instead!' 120 CALL QUIT('Wavelength too small under *ABALNR') 121 END IF 122 FRVAL(I) = XTNM/FRVAL(I) 123 END DO 124 NFRVAL = NFRVAL + NWAVEL 125 ICHANG = ICHANG + 1 126 GO TO 100 127 3 CONTINUE 128 READ (LUCMD,*) INPVAL 129 ICHANG = ICHANG + 1 130 IF (INPVAL .EQ. MAXITE) THEN 131 ICHANG = ICHANG - 1 132 ELSE 133 MAXITE = INPVAL 134 END IF 135 ABS_MAXITER=INPVAL 136 GO TO 100 137C .THRESH 138 4 CONTINUE 139 READ (LUCMD,*) DTHCLN 140 ICHANG = ICHANG + 1 141 IF (DTHCLN .EQ. THCLNR) THEN 142 ICHANG = ICHANG - 1 143 ELSE 144 THCLNR = DTHCLN 145 END IF 146 GO TO 100 147 5 CONTINUE 148 READ (LUCMD,*) IMXRM 149 ICHANG = ICHANG + 1 150 IF (IMXRM .EQ. MXRM) THEN 151 ICHANG = ICHANG - 1 152 ELSE 153 MXRM = IMXRM 154 END IF 155 ABS_MAXRM=IMXRM 156 GO TO 100 157 6 CONTINUE 158 READ (LUCMD,*) IMXPHP 159 ICHANG = ICHANG + 1 160 IF (IMXPHP .EQ. MXPHP) THEN 161 ICHANG = ICHANG - 1 162 ELSE 163 MXPHP = IMXPHP 164 END IF 165 GO TO 100 166 7 CONTINUE 167 STATIC = .TRUE. 168 ICHANG = ICHANG + 1 169 GO TO 100 170 8 CONTINUE ! .VIBGREP 171CSPAS:20/3-2011: allowing for perturbations in only one irrep 172 VIBGIR = .TRUE. 173 READ(LUCMD,*) IRVIBG 174 GO TO 100 175 9 CONTINUE 176 OOTV = .TRUE. 177 ICHANG = ICHANG + 1 178 GO TO 100 179 10 CONTINUE ! .FREQUE 180 READ(LUCMD,*) NRDFR 181 IF (NRDFR .LT. 0) THEN 182 CALL QUIT('# frequencies negative under *ABALNR') 183 ELSE 184 IF (.NOT. LOCFIN) THEN 185 NFRVAL = 0 186 LOCFIN = .TRUE. 187 END IF 188 IF (NRDFR .EQ. 0) GO TO 100 189 END IF 190 NFRVLT = NFRVLT + NRDFR 191 NFTOT = NRDFR + NFRVAL 192 IF (NFTOT .GT. MXFR) THEN 193 NRDFR = MXFR - NFRVAL 194 END IF 195 READ(LUCMD,*) (FRVAL(NFRVAL+I),I=1,NRDFR) 196 NFRVAL = NFRVAL + NRDFR 197 ICHANG = ICHANG + 1 198 GO TO 100 199 11 CONTINUE 200 READ(LUCMD,*) NFREQ_BATCH 201 IF (NFREQ_BATCH.GT.MXFR) THEN 202 NFREQ_BATCH = MXFR 203 END IF 204 GO TO 100 205 12 CONTINUE 206 READ (LUCMD,*) IPRLNR 207 ICHANG = ICHANG + 1 208 IF (IPRLNR .EQ. IPRDEF) ICHANG = ICHANG - 1 209 GO TO 100 210 13 CONTINUE 211 ABSORP=.TRUE. 212 ICHANG = ICHANG + 1 213 READ (LUCMD,*) DAMPING 214 ABS_DAMP=DAMPING 215 IF (.NOT. ABS_OLDCPP) THEN 216 ABSLRS=.TRUE. 217 ENDIF 218 GO TO 100 219 14 CONTINUE 220 CUT = .TRUE. 221 ICHANG = ICHANG + 1 222 GO TO 100 223 15 CONTINUE 224 ABS_INTERVAL = .TRUE. 225 ABS_REDUCE = .TRUE. 226 READ(LUCMD,*) (FREQ_INTERVAL(I), I=1,3) 227 GO TO 100 228 16 CONTINUE 229 ABS_REDUCE=.TRUE. 230 GO TO 100 231 17 CONTINUE 232 ABS_REBD=.FALSE. 233 GO TO 100 234 18 CONTINUE 235 ABS_OLDCPP = .TRUE. 236 ABSLRS =.FALSE. 237 GO TO 100 238 ELSE IF (PROMPT .EQ. '*') THEN 239 GO TO 300 240 ELSE 241 WRITE (LUPRI,'(/,3A,/)') ' Prompt "',WORD, 242 * '" not recognized in LNRINP.' 243 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 244 CALL QUIT('Illegal prompt under *ABALNR.') 245 END IF 246 END IF 247 300 CONTINUE 248C 249 IF (THR_REDFAC .GT. 0.0D0) THEN 250 ICHANG = ICHANG + 1 251 WRITE (LUPRI,'(3A,1P,D10.2)') '@ INFO ',WORD1, 252 & ' thresholds multiplied with general factor',THR_REDFAC 253 DTHCLN = DTHCLN*THR_REDFAC 254 END IF 255 CALL IZERO(NOPER,8) 256 DO ILABEL=1,3 257 ISYM = ISYMAX(ILABEL,1) + 1 258 NOPER(ISYM) = NOPER(ISYM) + 1 259 LABOP(NOPER(ISYM),ISYM) = DIPLEN(ILABEL) 260 IF(PCM.AND.LOCFLD) THEN 261 LABOP(NOPER(ISYM),ISYM) = DIPLOC(ILABEL) 262 END IF 263 END DO 264 IF (STATIC) THEN 265 DO ILABEL=1,3 266 ISYM = ISYMAX(ILABEL,2) + 1 267 NOPER(ISYM) = NOPER(ISYM) + 1 268 LABOP(NOPER(ISYM),ISYM) = ANGMOM(ILABEL) 269 END DO 270 DO ILABEL=1,3 271 ISYM = ISYMAX(ILABEL,2) + 1 272 NOPER(ISYM) = NOPER(ISYM) + 1 273 LABOP(NOPER(ISYM),ISYM) = LONMAG(ILABEL) 274 END DO 275 END IF 276C 277 IF ((VROA .OR. RAMAN .OR. OPTROT) .AND. NODIFC) THEN 278 NWARN = NWARN + 1 279 WRITE (LUPRI,'(/2A/A/A/)') 280 & ' WARNING: Raman properties calculation is NOT ', 281 & ' implemented with the .NODIFC keyword active.', 282 & ' WARNING: Raman properties calculation ignored', 283 & ' WARNING: Try calculation again without specifying .NODIFC.' 284 ROAA = .FALSE. 285 ROAG = .FALSE. 286 ELSEIF (VROA) THEN 287 ROAA = .TRUE. 288 ROAG = .TRUE. 289 ELSE IF (OPTROT) THEN 290 ROAA = .FALSE. 291 ROAG = .TRUE. 292 ELSE IF (RAMAN) THEN 293 ROAA = .FALSE. 294 ROAG = .FALSE. 295 ALFA = .TRUE. 296 ELSE 297 ROAA = .FALSE. 298 ROAG = .FALSE. 299 END IF 300C 301 IF (VIB_G .AND. NODIFC) THEN 302 WRITE(LUPRI,'(/3A/)') 303 & ' WARNING: Vibrational g-factors are NOT', 304 & ' implemented with the .NODIFC keyword active.', 305 & ' Try calculation again without specifying .NODIFC.' 306 VIB_G = .FALSE. 307 END IF 308C 309 IF (ICHANG .GT. 0) THEN 310 CALL HEADER('Changes of defaults for *ABALNR:',0) 311 IF (SKIP) THEN 312 WRITE (LUPRI,'(A)') ' LNRABA skipped in this run.' 313 ELSE 314 IF (ABSORP) THEN 315 IF (ABS_INTERVAL) THEN 316 IF (FREQ_INTERVAL(1).GT.FREQ_INTERVAL(2) .OR. 317 & (FREQ_INTERVAL(2)-FREQ_INTERVAL(1)).LT. 318 & FREQ_INTERVAL(3).OR. FREQ_INTERVAL(3).LE.0.0D0) 319 & THEN 320 CALL QUIT('.FREQ_INTERVAL not specify correctly') 321 END IF 322 DSMALL=1.0000001 323 NFREQ_INTERVAL = 1 + 324 & INT(DSMALL*(FREQ_INTERVAL(2)-FREQ_INTERVAL(1))/ 325 & FREQ_INTERVAL(3)) 326 NFRVAL = NFREQ_INTERVAL 327 DO I=1,NFRVAL 328 FRVAL(I)=FREQ_INTERVAL(1) + (I-1)*FREQ_INTERVAL(3) 329 END DO 330 ENDIF 331 332 ABS_MAXRM = MAX(MXRM,ABS_MAXRM) 333 MXRM = MAX(MXRM,ABS_MAXRM) 334 335 NFREQ_ALPHA = NFRVAL 336 NFREQ_INTERVAL = NFRVAL 337 ABS_NFREQ_ALPHA = NFRVAL 338 ABS_NFREQ_INTERVAL = NFRVAL 339 DO I = 1,ABS_NFREQ_ALPHA 340 FREQ_ALPHA(I) = FRVAL(I) 341 ABS_FREQ_ALPHA(I) = FRVAL(I) 342 END DO 343 THCLR_ABSORP = THCLNR 344 IPRABS = IPRLNR 345 MAX_MACRO = MAXITE 346 END IF 347 IF (NFRVAL .GT. 0) THEN 348 IF (NFRVLT .EQ. 0) NFRVLT = NFRVAL 349 J = MIN(MXFR,NFRVAL) 350 WRITE (LUPRI,'(A,I4/A,5F10.6:/,(33X,5F10.6))') 351 & ' Number of frequencies :',NFRVLT, 352 & ' Frequencies :', 353 & (FRVAL(I), I=1,J) 354 IF (NFRVLT .GT. MXFR) THEN 355 WRITE(LUPRI,'(/A/A,I5/A)') 356 & 'You have asked for too many frequencies under *ABALNR', 357 & 'Current maximum is MXFR =',MXFR, 358 & 'Reduce number of frequencies or increase'// 359 & ' MXFR in cbilnr.h and recompile.' 360 CALL QUIT( 361 & 'You have asked for too many frequencies under *ABALNR') 362 END IF 363 END IF 364 IF (STATIC) WRITE(LUPRI,'(A)') 365 & ' (Also) using static approximation for'// 366 & ' optical activity G tensor.' 367 IF (ABSORP) THEN 368 WRITE(LUPRI,'(A,F10.6,A)') 369 & ' Damping parameter equals :',DAMPING, ' a.u.' 370 END IF 371 WRITE (LUPRI,'(A,I4)') 372 & ' Print level :',IPRLNR 373 WRITE (LUPRI,'(A,1P,D10.2)') 374 & ' Lin Resp convergence threshold :',THCLNR 375 WRITE (LUPRI,'(A,I4)') 376 & ' Max. Lin Resp iterations :',MAXITE 377 IF (OOTV) WRITE (LUPRI,'(A)') 378 & ' Optimal orbital microiterations used.' 379 IF (CUT) THEN 380 WRITE (LUPRI,'(/A)') ' Program is stopped after LNRABA.' 381 END IF 382 END IF 383 END IF 384C 385 RETURN 386 END 387C /* Deck lnrini */ 388 SUBROUTINE LNRINI 389C 390C Initialize with default values for LNRABA module. 391C Values can be changed by user in input in LNRINP. 392C 393#include "implicit.h" 394#include "mxcent.h" 395#include "cbilnr.h" 396#include "cbiexc.h" 397#include "abainf.h" 398C 399! in cbilnr.h : 400 THCLNR = 5.D-05 401 FRVAL(1) = 0.0D0 402 NFRVAL = 1 403 IPRLNR = IPRDEF ! IPRDEF from abainf.h 404 ALFA = ABA_ALPHA ! ABA_ALPHA from abainf.h 405 STATIC = .FALSE. 406! in cbiexc.h : 407 MAXITE = 60 408 MXRM = 400 409 MXPHP = 0 410 IPR1IN = IPRDEF ! IPRDEF from abainf.h 411 SKIP = .FALSE. 412 CUT = .FALSE. 413 OOTV = .FALSE. 414! in abainf.h : 415 VIBGIR = .FALSE. ! SPAS:20/3-2011: allowing for perturbations in only one irrep 416 RETURN 417 END 418C /* Deck lnraba */ 419 SUBROUTINE LNRABA(POLDD,POLDQ,POLDL,POLDA,POLVL,POLVV,FOVIBG, 420 & WORK,LWORK,PASS) 421C 422C Frequency-dependent Linear Response module in ABACUS 423C 424C Uses the general linear response solver in the RESPONS program unit 425C 426#include "implicit.h" 427#include "dummy.h" 428#include "mxcent.h" 429#include "maxorb.h" 430#include "iratdef.h" 431#include "priunit.h" 432#include "cbilnr.h" 433 LOGICAL PASS 434 LOGICAL CICLC, HFCLC, TRIPLE, EXECLC, FOUND, CONV 435 DIMENSION WORK(LWORK),SNDPRP(2) 436 DIMENSION POLDD(2,3,3,MXFR), POLDQ(2,3,3,3,MXFR) 437 DIMENSION POLDL(2,3,3,MXFR), POLDA(2,3,3,MXFR) 438 DIMENSION POLVL(2,3,3,MXFR), POLVV(2,3,3,MXFR) 439 DIMENSION POLDLS(3,3), POLDAS(3,3) 440CSPAS:6/10-08: trying to add mass independent vibrational g-factor 441C DIMENSION FOVIBG(3*NATOMS,3*NATOMS,NSYM,MXFR) 442 DIMENSION FOVIBG(3*NATOMS+3,3*NATOMS+3,NSYM,MXFR) 443CKeinSPASmehr 444 CHARACTER*8 LABEL1, LABEL2, LABINT(3*MXCOOR), BLANK 445 PARAMETER (DP5=0.5D0) 446 PARAMETER (LOCDBG = 3) 447C 448#include "orgcom.h" 449#include "cbiexc.h" 450#include "inflin.h" 451#include "infvar.h" 452#include "infdim.h" 453#include "inforb.h" 454#include "nuclei.h" 455#include "inftap.h" 456#include "infrsp.h" 457#include "wrkrsp.h" 458#include "maxmom.h" 459#include "maxaqn.h" 460#include "symmet.h" 461#include "abainf.h" 462#include "gnrinf.h" 463#include "infsop.h" 464#include "absorp.h" 465#include "abslrs.h" 466#include "pcmlog.h" 467C 468 IF (SKIP) RETURN 469 CALL QENTER('LNRABA') 470 CALL TIMER('START ',TIMEIN,TIMOUT) 471C 472 IF (IPRLNR .GE. 0) 473 & CALL TITLER('Solving Linear Response Equations','*',118) 474C 475 IPRRSP = IPRLNR 476C 477C Get reference state 478C =================== 479C 480C 1. Work Allocations: 481C 482 IF (ABASOP) THEN 483 LUDV = NORBT * NORBT 484 LPVX = LPVMAT 485 ELSE 486 LUDV = N2ASHX 487 LPVX = 0 488 ENDIF 489 KFREE = 1 490 LFREE = LWORK 491C 492 CALL MEMGET('REAL',KCMO ,NCMOT ,WORK,KFREE,LFREE) 493 CALL MEMGET('REAL',KUDV ,LUDV ,WORK,KFREE,LFREE) 494 CALL MEMGET('REAL',KPVX ,LPVX ,WORK,KFREE,LFREE) 495 CALL MEMGET('REAL',KXINDX,LCINDX,WORK,KFREE,LFREE) 496C 497 KWORK1 = KFREE 498 LWORK1 = LFREE 499C 500 CALL RD_SIRIFC('CMO',FOUND,WORK(KCMO)) 501 IF (.NOT.FOUND) CALL QUIT('LNRABA error: CMO not found on SIRIFC') 502 IF (NASHT .GT. 0) THEN 503 CALL RD_SIRIFC('DV',FOUND,WORK(KWORK1)) 504 IF (.NOT.FOUND) 505 & CALL QUIT('LNRABA error: DV not found on SIRIFC') 506 CALL DSPTSI(NASHT,WORK(KWORK1),WORK(KUDV)) 507 END IF 508C 509 ISYM = 1 510 CALL LNRVAR(ISYM,IPRLNR,WORK(KWORK1),LWORK1) 511C 512 CALL GETCIX(WORK(KXINDX),IREFSY,IREFSY,WORK(KWORK1),LWORK1,0) 513C 514C 515C SOPPA : 516C 517 IF (ABASOP) THEN 518C 519C Initialize XINDX 520C 521 CALL DZERO(WORK(KXINDX),LCINDX) 522C 523C Find address array's for SOPPA calculation 524C 525 CALL SET2SOPPA(WORK(KXINDX+KABSAD-1),WORK(KXINDX+KABTAD-1), 526 * WORK(KXINDX+KIJSAD-1),WORK(KXINDX+KIJTAD-1), 527 * WORK(KXINDX+KIJ1AD-1),WORK(KXINDX+KIJ2AD-1), 528 * WORK(KXINDX+KIJ3AD-1),WORK(KXINDX+KIADR1-1)) 529C 530C 531 REWIND (LUSIFC) 532 IF (CCPPA) THEN 533 CALL MOLLAB('CCSDINFO',LUSIFC,LUPRI) 534 ELSE 535 CALL MOLLAB('MP2INFO ',LUSIFC,LUPRI) 536 ENDIF 537C 538C reads the MP2 or CCSD correlation coefficients into PV 539C 540 CALL READT (LUSIFC,LPVMAT,WORK(KPVX)) 541C 542 IF (IPRLNR.GT.10) THEN 543 IF (CCPPA) THEN 544 WRITE(LUPRI,'(/A)')' EXCIT1 : CCSD correlation ', 545 & 'coefficients' 546 ELSE 547 WRITE(LUPRI,'(/A,A)')' EXCIT1 :', 548 & ' MP2 correlation coefficients' 549 ENDIF 550 CALL OUTPUT(WORK(KPVX),1,LPVMAT,1,1,LPVMAT,1,1,LUPRI) 551 END IF 552C 553C reads the MP2 or CCSD second order one particle density matrix 554C 555 CALL READT (LUSIFC,NORBT*NORBT,WORK(KUDV)) 556C 557C UDV contains the MP2 one-density. Remove the diagonal 558C contribution from the zeroth order. (Added in MP2FAC) 559C 560 IF (IPRLNR.GT.10) THEN 561 IF (CCPPA) THEN 562 WRITE(LUPRI,'(/A)')' RSPMC : CCSD density' 563 ELSE 564 WRITE(LUPRI,'(/A)')' RSPMC : MP2 density' 565 END IF 566 CALL OUTPUT(WORK(KUDV),1,NORBT*NORBT,1,1,NORBT*NORBT,1,1, 567 & LUPRI) 568 END IF 569C 570 CALL SOPUDV(WORK(KUDV)) 571 END IF 572C 573C Construct property-integrals and write to LUPROP 574C ================================================ 575C 576C 2. Work Allocations: 577C 578 KIDSYM = KWORK1 579 KIDADR = KIDSYM + 9*MXCENT 580 KWORK2 = KIDADR + 9*MXCENT 581 LWORK2 = LWORK - KWORK2 582C 583 NLBTOT = 0 584C 585 IPR1IN=MAX(0,IPRLNR-5) 586 IF (ALFA .OR. ROAA .OR. ROAG) THEN 587 NCOMP = 0 588 NPATOM = 0 589Clf 590 IF(PCM.AND.LOCFLD) THEN 591 CALL GET1IN(DUMMY,'DIPLOC ',NCOMP,WORK(KWORK2),LWORK2, 592 & LABINT,WORK(KIDSYM),WORK(KIDADR), 593 & IDUMMY,.TRUE.,NPATOM,.TRUE.,DUMMY,.FALSE.,DUMMY,IPR1IN) 594 NLAB = 3 595 CALL LABCOP(NLAB,NLBTOT,LABINT,WORK(KIDSYM),LABAPP,LABSYM) 596 ELSE 597 CALL GET1IN(DUMMY,'DIPLEN ',NCOMP,WORK(KWORK2),LWORK2, 598 & LABINT,WORK(KIDSYM),WORK(KIDADR), 599 & IDUMMY,.TRUE.,NPATOM,.TRUE.,DUMMY,.FALSE.,DUMMY,IPR1IN) 600 NLAB = 3 601 CALL LABCOP(NLAB,NLBTOT,LABINT,WORK(KIDSYM),LABAPP,LABSYM) 602 END IF 603 ENDIF 604 IF (ROAA) THEN 605 NCOMP = 0 606 NPATOM = 0 607 DIPTMX = DIPORG(1) 608 DIPTMY = DIPORG(2) 609 DIPTMZ = DIPORG(3) 610 DIPORG(1) = 0.0D0 611 DIPORG(2) = 0.0D0 612 DIPORG(3) = 0.0D0 613 CALL GET1IN(DUMMY,'THETA ',NCOMP,WORK(KWORK2),LWORK2, 614 & LABINT,WORK(KIDSYM),WORK(KIDADR), 615 & IDUMMY,.TRUE.,NPATOM,.TRUE.,DUMMY,.FALSE.,DUMMY, 616 & IPR1IN) 617 DIPORG(1) = DIPTMX 618 DIPORG(2) = DIPTMY 619 DIPORG(3) = DIPTMZ 620 NLAB = 6 621 CALL LABCOP(NLAB,NLBTOT,LABINT,WORK(KIDSYM),LABAPP,LABSYM) 622 END IF 623C 624 IF (ROAG) THEN 625C 626 CALL LABCOP(1,NLBTOT,'XLONMAG ',ISYMAX(1,2),LABAPP,LABSYM) 627 CALL LABCOP(1,NLBTOT,'YLONMAG ',ISYMAX(2,2),LABAPP,LABSYM) 628 CALL LABCOP(1,NLBTOT,'ZLONMAG ',ISYMAX(3,2),LABAPP,LABSYM) 629C 630 NCOMP = 0 631 NPATOM = 0 632 CALL GET1IN(DUMMY,'ANGMOM ',NCOMP,WORK(KWORK2),LWORK2, 633 & LABINT,WORK(KIDSYM),WORK(KIDADR), 634 & IDUMMY,.TRUE.,NPATOM,.TRUE.,DUMMY,.FALSE.,DUMMY, 635 & IPR1IN) 636 NLAB = 3 637 CALL LABCOP(NLAB,NLBTOT,LABINT,WORK(KIDSYM),LABAPP,LABSYM) 638C 639 IF (MVEOR) THEN 640 NCOMP = 0 641 NPATOM = 0 642 CALL GET1IN(DUMMY,'DIPVEL ',NCOMP,WORK(KWORK2),LWORK2, 643 & LABINT,WORK(KIDSYM),WORK(KIDADR), 644 & IDUMMY,.TRUE.,NPATOM,.TRUE.,DUMMY,.FALSE.,DUMMY, 645 & IPR1IN) 646 NLAB = 3 647 CALL LABCOP(NLAB,NLBTOT,LABINT,WORK(KIDSYM),LABAPP,LABSYM) 648 END IF 649C 650 END IF 651C 652CSPAS:6/10-08: trying to add mass independent vibrational g-factor 653C 654 IF (VIB_G) THEN 655 NCOMP = 0 656 NPATOM = 0 657 CALL GET1IN(DUMMY,'DIPVEL ',NCOMP,WORK(KWORK2),LWORK2, 658 & LABINT,WORK(KIDSYM),WORK(KIDADR), 659 & IDUMMY,.TRUE.,NPATOM,.TRUE.,DUMMY,.FALSE.,DUMMY, 660 & IPR1IN) 661 NLAB = 3 662 CALL LABCOP(NLAB,NLBTOT,LABINT,WORK(KIDSYM),LABAPP,LABSYM) 663 END IF 664C 665CKeinSPASmehr 666C 667C Start calculation of polarizabilities and raman optical activity. 668C ================================================================= 669C 670 IF (ALFA .OR. ROAA .OR. ROAG) THEN 671C 672C Set variables and logicals 673C 674 CICLC = .FALSE. 675 HFCLC = NASHT .LE. 1 676 TRIPLE = .FALSE. 677 EXECLC = .FALSE. 678 NABATY = 1 679 TRPLET = .FALSE. 680 NABAOP = 1 681 BLANK = ' ' 682C 683 IF (NFRVAL.LE.0) GOTO 999 684C 685C Zero the property complex tensors 686C 687 CALL DZERO(SNDPRP,2) 688 CALL DZERO(POLDD,2*9*MXFR) 689 CALL DZERO(POLDQ,2*27*MXFR) 690 CALL DZERO(POLDL,2*9*MXFR) 691 CALL DZERO(POLDA,2*9*MXFR) 692 CALL DZERO(POLVL,2*9*MXFR) 693 CALL DZERO(POLVV,2*9*MXFR) 694 CALL DZERO(POLDLS,9) 695 CALL DZERO(POLDAS,9) 696C 697C Start the calculations 698C 699C tbp, jan. 2004: pass twice to get DIPVEL as first operator 700C for velocity gauge optical rotation adjusted 701C to ensure vanishing rotation at zero frequency 702C (the modified velocity gauge). 703C 704 IFRZER = 0 705 NFRSAV = NFRVAL 706 IF (ROAG .AND. MVEOR .AND. NFRVAL.GT.0) THEN 707 NUMRUN = 2 708 IFRZER = 0 709 XMINFR = 1.0D15 710 DO IFRVAL = 1,NFRVAL 711 TSTFR = ABS(FRVAL(IFRVAL)) 712 IF (TSTFR.LT.1.0D-10 .AND. TSTFR.LT.XMINFR) THEN 713 XMINFR = TSTFR 714 IFRZER = IFRVAL 715 END IF 716 END DO 717 ELSE 718 NUMRUN = 1 719 END IF 720C 721 DO IRUN = 1,NUMRUN 722C 723 LUSOVE = -1 724 LUGDVE = -1 725 LUREVE = -1 726 CALL GPOPEN(LUSOVE,' ','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 727 CALL GPOPEN(LUGDVE,' ','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 728 CALL GPOPEN(LUREVE,' ','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 729 IF (IRUN .EQ. 2) THEN 730 IF (IFRZER .EQ. 0) THEN 731 IF (NFRVAL.GE.MXFR) THEN 732 CALL QUIT('Too many frequencies in LNRABA') 733 ELSE 734 NFRVAL = NFRVAL + 1 735 NFREQ_ALPHA = NFRVAL 736 IFRZER = NFRVAL 737 FRVAL(IFRZER) = 0.0D0 738 FREQ_ALPHA(IFRZER) = FRVAL(IFRZER) 739 END IF 740 END IF 741 END IF 742C 743 DO ISYM=1,NSYM 744CSPAS:20/3-2011: allowing for VIB_G perturbations in only one irrep 745 IF (VIB_G .AND. VIBGIR .AND. (ISYM .NE. IRVIBG)) CYCLE 746CKeinSPASmehr 747 KSYMOP = ISYM 748 IF (ABSLRS) THEN 749 ABS_KLRED(1)=0 750 ABS_KLRED(2)=0 751C 752 LUSB = -1 753 LUAB = -1 754 LUSS = -1 755 LUAS = -1 756C 757 CALL GPOPEN(LUSB,'ABS_SB','NEW',' ',' ',IDUMMY,.FALSE.) 758 CALL GPOPEN(LUAB,'ABS_AB','NEW',' ',' ',IDUMMY,.FALSE.) 759 CALL GPOPEN(LUSS,'ABS_SS','NEW',' ',' ',IDUMMY,.FALSE.) 760 CALL GPOPEN(LUAS,'ABS_AS','NEW',' ',' ',IDUMMY,.FALSE.) 761C 762 LUE1RED = -1 763 LUE2RED = -1 764 LUSRED = -1 765 CALL GPOPEN(LUE1RED,'ABS_E1RED','NEW', 766 & ' ',' ',IDUMMY,.FALSE.) 767 CALL GPOPEN(LUE2RED,'ABS_E2RED','NEW', 768 & ' ',' ',IDUMMY,.FALSE.) 769 CALL GPOPEN(LUSRED ,'ABS_SRED' ,'NEW', 770 & ' ',' ',IDUMMY,.FALSE.) 771 ENDIF 772C 773 DO IOPER=1,NOPER(KSYMOP) 774 KOPER=IOPER 775 LABEL1 = LABOP(IOPER,KSYMOP) 776 IDIP = 0 777 IF ((LABEL1 .EQ.'XDIPLEN').OR.(LABEL1 .EQ.'XDIPLOC')) 778 $ IDIP = 1 779 IF ((LABEL1 .EQ.'YDIPLEN').OR.(LABEL1 .EQ.'YDIPLOC')) 780 $ IDIP = 2 781 IF ((LABEL1 .EQ.'ZDIPLEN').OR.(LABEL1 .EQ.'ZDIPLOC')) 782 $ IDIP = 3 783 IF ((LABEL1(2:7) .EQ.'DIPLEN') .OR. (LABEL1(2:7) .EQ. 784 & 'DIPLOC')) THEN 785 NABATY = 1 786 ELSE IF ((LABEL1(2:7) .EQ. 'ANGMOM') .OR. 787 & (LABEL1(2:7) .EQ. 'LONMAG')) THEN 788 NABATY = -1 789 ELSE 790 WRITE (LUPRI,'(3A)') 'ERROR Operator label "',LABEL1, 791 & '" not recognised' 792 CALL QUIT('Unrecognised operator in LNRABA') 793 END IF 794 IF (IRUN .EQ. 2) THEN 795 IF (IDIP .EQ. 1) THEN 796 LABEL1 = 'XDIPVEL ' 797 ELSE IF (IDIP .EQ. 2) THEN 798 LABEL1 = 'YDIPVEL ' 799 ELSE IF (IDIP .EQ. 3) THEN 800 LABEL1 = 'ZDIPVEL ' 801 ELSE 802 GO TO 300 ! cycle operator loop 803 END IF 804 NABATY = -1 ! DIPVEL is anti-Hermitian 805 END IF 806C 807 CALL LNRVAR(ISYM,IPRLNR,WORK(KWORK2),LWORK2) 808C 809 KGD1 = KWORK1 810 KWRKG1 = KGD1 811 LWRKG1 = LWORK - KWRKG1 812 KSLV = KGD1 + 2*NVARPT 813 KLAST = KSLV + 4*NVARPT 814 IF (KLAST.GT.LWORK) CALL STOPIT('LNRABA',' ',KLAST,LWORK) 815 KWRK = KLAST 816 LWRK = LWORK - KLAST + 1 817C 818C Find right hand side (gradient) of first operator and write to file 819C =================================================================== 820C 821 CALL GETGPV(LABEL1,DUMMY,DUMMY,WORK(KCMO),WORK(KUDV), 822 & WORK(KPVX),WORK(KXINDX),ANTSYM,WORK(KWRKG1),LWRKG1) 823 REWIND LUGDVE 824 CALL WRITT(LUGDVE,2*NVARPT,WORK(KWRKG1)) 825 IF (IPRLNR.GE.LOCDBG) THEN 826 WRITE (LUPRI,'(/3A,1P,D15.6)') 827 & 'GP Vector, label: ',LABEL1,' Norm: ', 828 & DSQRT(DDOT(2*NVARPT,WORK(KWRKG1),1,WORK(KWRKG1),1)) 829 IF (IPRLNR.GT.10) THEN 830 CALL OUTPUT(WORK(KGD1),1,NVARPT,1,2,NVARPT, 831 & 2,1,LUPRI) 832 END IF 833 END IF 834C 835C Calculate linear response vector and write to file 836C ================================================== 837C 838 CALL ABARSP(CICLC,HFCLC,TRIPLE,OOTV,ISYM,EXECLC, 839 & FRVAL,NFRVAL,NABATY,NABAOP,LABEL1,LUGDVE,LUSOVE, 840 & LUREVE,THCLNR,MAXITE,IPRRSP,MXRM,MXPHP, 841 & WORK(KWRK),LWRK) 842C 843C Loop over the second property operators 844C ======================================= 845C 846 DO 200 IPRLBL = 1, NLBTOT 847C 848C Find label and symmetry of second operator 849C ========================================== 850C 851 LABEL2 = LABAPP(IPRLBL) 852 KSYM = LABSYM(IPRLBL) 853C 854C Only ANGMOM in second run (for vel. and mod. vel. opt. rot.). 855C ============================================================= 856C 857 IF (IRUN.EQ.2 .AND. LABEL2(2:7).NE.'ANGMOM') 858 & GO TO 200 859C 860C If symmetry of first operator equals symmetry of 861C second operator, that is if ISYM = KSYM, then 862C ================================================ 863C 864 IF (KSYM.EQ.ISYM) THEN 865C 866 IF (IPRLNR.GE.LOCDBG) THEN 867 WRITE(LUPRI,'(/,A,A,A,A,A,/)') 868 & '***** RESPONSE CALCULATION FOR ',LABEL1,' ', 869 & LABEL2,' OPERATOR DOUBLE *****' 870 CALL FLSHFO(LUPRI) 871 END IF 872C 873C Find right hand side (gradient) for second operator 874C =================================================== 875C 876 CALL GETGPV(LABEL2,DUMMY,DUMMY,WORK(KCMO), 877 & WORK(KUDV),WORK(KPVX),WORK(KXINDX),ANTSYM, 878 & WORK(KWRKG1),LWRKG1) 879C 880 IF (IPRLNR.GE.LOCDBG) THEN 881 WRITE (LUPRI,'(/3A,1P,D15.6)') 882 & 'GP Vector, label: ',LABEL2,' Norm: ', 883 & DSQRT(DDOT(2*NVARPT,WORK(KGD1),1,WORK(KGD1),1)) 884 IF (IPRLNR.GT.10) THEN 885 CALL OUTPUT(WORK(KGD1),1,NVARPT,1,2,NVARPT, 886 & 2,1,LUPRI) 887 END IF 888 ENDIF 889C 890C Form second order properties SNDPRP 891C =================================== 892C 893 IF (.NOT.ABSORP) REWIND LUSOVE 894 DO 100 IFRVAL = 1,NFRVAL 895 IF (ABSORP) THEN 896 IF (ABSLRS) THEN 897 LUABSVECS = -1 898 CALL GPOPEN(LUABSVECS,'ABSVECS','OLD',' ', 899 & ' ',IDUMMY,.FALSE.) 900C CALL READ_XVEC(LUABSVECS,4*NVARPT, 901C & WORK(KSLV),LABEL1,ISYM, 902C & FRVAL(IFRVAL),ABS_THCLR,FOUND,CONV) 903 CALL READ_XVEC2(LUABSVECS,4*NVARPT, 904 & WORK(KSLV),LABEL1,' ',ISYM, 905 & FRVAL(IFRVAL),0.0D0,ABS_THCLR,FOUND,CONV) 906 CALL GPCLOSE(LUABSVECS,'KEEP') 907 ENDIF 908 CALL GPOPEN(LURSP,'RSPVEC','OLD',' ', 909 & ' ',IDUMMY,.FALSE.) 910 CALL REARSP(LURSP,4*NVARPT,WORK(KSLV), 911 & LABEL1,'COMPLEX ',FRVAL(IFRVAL),0.0D0, 912 & ISYM,0,THCLR_ABSORP,FOUND,CONV,ANTSYM) 913 CALL GPCLOSE(LURSP,'KEEP') 914c ENDIF 915 SNDPRP(1)=DDOT(2*NVARPT,WORK(KSLV),1, 916 & WORK(KGD1),1) 917 SNDPRP(2)=DDOT(2*NVARPT,WORK(KSLV+2*NVARPT),1, 918 & WORK(KGD1),1) 919 ELSE 920 CALL READT(LUSOVE,2*NVARPT,WORK(KSLV)) 921 SNDPRP(1)=DDOT(2*NVARPT,WORK(KSLV),1, 922 & WORK(KGD1),1) 923 END IF 924C 925 IF (IPRLNR.GE.LOCDBG) THEN 926 WRITE (LUPRI,'(A,I4,A,1P,D15.6)') 927 & 'Solution Vector no. ',IFRVAL,' Norm: ', 928 & DSQRT(DDOT(2*NVARPT,WORK(KSLV),1, 929 & WORK(KSLV),1)) 930 IF (IPRLNR.GT.10) THEN 931 CALL OUTPUT(WORK(KSLV),1,NVARPT,1,2, 932 & NVARPT,2,1,LUPRI) 933 END IF 934 ENDIF 935 IF (IPRLNR.GE.(LOCDBG-1) .OR. IDIP.EQ.0) THEN 936 WRITE (LUPRI,'(/,A,F15.8)') 937 & ' Frequency = ',FRVAL(IFRVAL) 938 WRITE (LUPRI,'(4A,F15.8)') 939 & ' Second order property for ', 940 & LABEL1,LABEL2,' = ',SNDPRP(1) 941 ENDIF 942C 943C Write properties into the various property matrices 944C =================================================== 945C Skip this section if first operator is not a dipole 946C operator (tbp, 2004). 947C 948 IF (IDIP.GT.0 .AND. IDIP.LT.4) THEN 949 IF (LABEL2(2:7).EQ.'DIPLEN') THEN 950 IF (LABEL2(1:1).EQ.'X') THEN 951 POLDD(1,IDIP,1,IFRVAL) = SNDPRP(1) 952 POLDD(2,IDIP,1,IFRVAL) = SNDPRP(2) 953 END IF 954 IF (LABEL2(1:1).EQ.'Y') THEN 955 POLDD(1,IDIP,2,IFRVAL) = SNDPRP(1) 956 POLDD(2,IDIP,2,IFRVAL) = SNDPRP(2) 957 END IF 958 IF (LABEL2(1:1).EQ.'Z') THEN 959 POLDD(1,IDIP,3,IFRVAL) = SNDPRP(1) 960 POLDD(2,IDIP,3,IFRVAL) = SNDPRP(2) 961 END IF 962Clf 963 ELSE IF (LABEL2(2:7).EQ.'DIPLOC') THEN 964 IF (LABEL2(1:1).EQ.'X') THEN 965 POLDD(1,IDIP,1,IFRVAL) = SNDPRP(1) 966 POLDD(2,IDIP,1,IFRVAL) = SNDPRP(2) 967 END IF 968 IF (LABEL2(1:1).EQ.'Y') THEN 969 POLDD(1,IDIP,2,IFRVAL) = SNDPRP(1) 970 POLDD(2,IDIP,2,IFRVAL) = SNDPRP(2) 971 END IF 972 IF (LABEL2(1:1).EQ.'Z') THEN 973 POLDD(1,IDIP,3,IFRVAL) = SNDPRP(1) 974 POLDD(2,IDIP,3,IFRVAL) = SNDPRP(2) 975 END IF 976Clf 977 ELSE IF (LABEL2(3:8).EQ.'THETA ') THEN 978 IF (LABEL2(1:2).EQ.'XX') THEN 979 POLDQ(1,IDIP,1,1,IFRVAL) = SNDPRP(1) 980 POLDQ(2,IDIP,1,1,IFRVAL) = SNDPRP(2) 981 END IF 982 IF (LABEL2(1:2).EQ.'XY') THEN 983 POLDQ(1,IDIP,1,2,IFRVAL) = SNDPRP(1) 984 POLDQ(2,IDIP,1,2,IFRVAL) = SNDPRP(2) 985 END IF 986 IF (LABEL2(1:2).EQ.'XY') THEN 987 POLDQ(1,IDIP,2,1,IFRVAL) = SNDPRP(1) 988 POLDQ(2,IDIP,2,1,IFRVAL) = SNDPRP(2) 989 END IF 990 IF (LABEL2(1:2).EQ.'XZ') THEN 991 POLDQ(1,IDIP,1,3,IFRVAL) = SNDPRP(1) 992 POLDQ(2,IDIP,1,3,IFRVAL) = SNDPRP(2) 993 END IF 994 IF (LABEL2(1:2).EQ.'XZ') THEN 995 POLDQ(1,IDIP,3,1,IFRVAL) = SNDPRP(1) 996 POLDQ(2,IDIP,3,1,IFRVAL) = SNDPRP(2) 997 END IF 998 IF (LABEL2(1:2).EQ.'YY') THEN 999 POLDQ(1,IDIP,2,2,IFRVAL) = SNDPRP(1) 1000 POLDQ(2,IDIP,2,2,IFRVAL) = SNDPRP(2) 1001 END IF 1002 IF (LABEL2(1:2).EQ.'YZ') THEN 1003 POLDQ(1,IDIP,2,3,IFRVAL) = SNDPRP(1) 1004 POLDQ(2,IDIP,2,3,IFRVAL) = SNDPRP(2) 1005 END IF 1006 IF (LABEL2(1:2).EQ.'YZ') THEN 1007 POLDQ(1,IDIP,3,2,IFRVAL) = SNDPRP(1) 1008 POLDQ(2,IDIP,3,2,IFRVAL) = SNDPRP(2) 1009 END IF 1010 IF (LABEL2(1:2).EQ.'ZZ') THEN 1011 POLDQ(1,IDIP,3,3,IFRVAL) = SNDPRP(1) 1012 POLDQ(2,IDIP,3,3,IFRVAL) = SNDPRP(2) 1013 END IF 1014 ELSE IF (LABEL2(2:7).EQ.'LONMAG') THEN 1015C 1016C Multiply with minus the Bohr-magneton (-0.5) to create the magnetic 1017C dipole operator from the angular momentum operator 1018C 1019 IF (LABEL2(1:1).EQ.'X') THEN 1020 POLDL(1,IDIP,1,IFRVAL) = -DP5*SNDPRP(1) 1021 POLDL(2,IDIP,1,IFRVAL) = -DP5*SNDPRP(2) 1022 END IF 1023 IF (LABEL2(1:1).EQ.'Y') THEN 1024 POLDL(1,IDIP,2,IFRVAL) = -DP5*SNDPRP(1) 1025 POLDL(2,IDIP,2,IFRVAL) = -DP5*SNDPRP(2) 1026 END IF 1027 IF (LABEL2(1:1).EQ.'Z') THEN 1028 POLDL(1,IDIP,3,IFRVAL) = -DP5*SNDPRP(1) 1029 POLDL(2,IDIP,3,IFRVAL) = -DP5*SNDPRP(2) 1030 END IF 1031 ELSE IF (LABEL2(2:7).EQ.'ANGMOM') THEN 1032 IF (IRUN .EQ. 2) THEN 1033 IF (IPRLNR.GT.(LOCDBG+10)) THEN 1034 WRITE(LUPRI,*) 1035 & '...property goes into POLVL' 1036 END IF 1037 IF (LABEL2(1:1).EQ.'X') THEN 1038 POLVL(1,IDIP,1,IFRVAL) = 1039 & DP5*SNDPRP(1) 1040 POLVL(2,IDIP,1,IFRVAL) = 1041 & DP5*SNDPRP(2) 1042 END IF 1043 IF (LABEL2(1:1).EQ.'Y') THEN 1044 POLVL(1,IDIP,2,IFRVAL) = 1045 & DP5*SNDPRP(1) 1046 POLVL(2,IDIP,2,IFRVAL) = 1047 & DP5*SNDPRP(2) 1048 END IF 1049 IF (LABEL2(1:1).EQ.'Z') THEN 1050 POLVL(1,IDIP,3,IFRVAL) = 1051 & DP5*SNDPRP(1) 1052 POLVL(2,IDIP,3,IFRVAL) = 1053 & DP5*SNDPRP(2) 1054 END IF 1055 ELSE 1056 IF (LABEL2(1:1).EQ.'X') THEN 1057 POLDA(1,IDIP,1,IFRVAL) = 1058 & -DP5*SNDPRP(1) 1059 POLDA(2,IDIP,1,IFRVAL) = 1060 & -DP5*SNDPRP(2) 1061 END IF 1062 IF (LABEL2(1:1).EQ.'Y') THEN 1063 POLDA(1,IDIP,2,IFRVAL) = 1064 & -DP5*SNDPRP(1) 1065 POLDA(2,IDIP,2,IFRVAL) = 1066 & -DP5*SNDPRP(2) 1067 END IF 1068 IF (LABEL2(1:1).EQ.'Z') THEN 1069 POLDA(1,IDIP,3,IFRVAL) = 1070 & -DP5*SNDPRP(1) 1071 POLDA(2,IDIP,3,IFRVAL) = 1072 & -DP5*SNDPRP(2) 1073 END IF 1074 END IF 1075C-tbp ELSE IF (IPRLNR.LE.2) THEN 1076C ... hjaaj mar 2004: this property is not 1077C predefined, should never happen, but 1078C we make sure the value is printed. 1079C (the value was printed above if IPRLNR.gt.2) 1080C-tbp WRITE (LUPRI,'(/A,F15.8,A/4A,F15.8)') 1081C-tbp& ' Frequency = ',FRVAL(IFRVAL),' au', 1082C-tbp& ' Second order property for ', 1083C-tbp& LABEL1,LABEL2,' = ',SNDPRP(1) 1084 END IF 1085 END IF 1086 100 CONTINUE 1087 ENDIF 1088 200 CONTINUE 1089C 1090 300 CONTINUE 1091C 1092C End loop over operator in this symmetry 1093C 1094 END DO 1095C 1096 IF (ABSLRS) THEN 1097 CALL GPCLOSE(LUSB,'DELETE') 1098 CALL GPCLOSE(LUAB,'DELETE') 1099 CALL GPCLOSE(LUSS,'DELETE') 1100 CALL GPCLOSE(LUAS,'DELETE') 1101C 1102 CALL GPCLOSE(LUE1RED,'DELETE') 1103 CALL GPCLOSE(LUE2RED,'DELETE') 1104 CALL GPCLOSE(LUSRED, 'DELETE') 1105 ENDIF 1106C 1107C End loop over symmetries 1108C 1109 END DO 1110C 1111C Close and delete files. 1112C 1113 IF (LUSOVE.GT.0) CALL GPCLOSE(LUSOVE,'DELETE') 1114 IF (LUGDVE.GT.0) CALL GPCLOSE(LUGDVE,'DELETE') 1115 IF (LUREVE.GT.0) CALL GPCLOSE(LUREVE,'DELETE') 1116C 1117C End repetition loop (IRUN) 1118C 1119 END DO 1120C 1121 NFRVAL = NFRSAV 1122 NFREQ_ALPHA = NFRVAL 1123 IF (ROAG .AND. MVEOR) THEN 1124 IF (IFRZER.LT.1 .OR. IFRZER.GT.MXFR) THEN 1125 CALL QUIT('Logical error in LNRABA') 1126 END IF 1127 DO IFRVAL = 1,NFRVAL 1128 DO IANG = 1,3 1129 DO IDIP = 1,3 1130 DO IR = 1,2 1131 POLVV(IR,IDIP,IANG,IFRVAL) = 1132 & POLVL(IR,IDIP,IANG,IFRVAL) 1133 & - POLVL(IR,IDIP,IANG,IFRZER) 1134 END DO 1135 END DO 1136 END DO 1137 END DO 1138 IF ((IPRLNR.GE.LOCDBG) .AND. (NFRVAL.GT.0)) THEN 1139 DO IFRVAL = 1,NFRVAL 1140 WRITE(LUPRI,'(/,A,1P,D15.6)') 'Frequency: ',FRVAL(IFRVAL) 1141 WRITE(LUPRI,'(A)') 'POLDL: London' 1142 WRITE(LUPRI,'(A)') 'POLDA: No-London' 1143 WRITE(LUPRI,'(A)') 'POLVL: Velocity' 1144 WRITE(LUPRI,'(A)') 'POLVV: Modified velocity' 1145 WRITE(LUPRI,'(A,A,/,A,A)') 1146 & ' Ei Bj POLDL POLDA POLVL', 1147 & ' POLVV', 1148 & '-------------------------------------------------', 1149 & '---------------------' 1150 DO IANG = 1,3 1151 DO IDIP = 1,3 1152 WRITE(LUPRI,'(1X,I2,1X,I2,1X,F15.7,1X,F15.7,1X,F15.7,1X,F15.7)') 1153 & IDIP,IANG,POLDL(1,IDIP,IANG,IFRVAL), 1154 & POLDA(1,IDIP,IANG,IFRVAL), 1155 & POLVL(1,IDIP,IANG,IFRVAL), 1156 & POLVV(1,IDIP,IANG,IFRVAL) 1157 IF (ABSORP) THEN 1158 WRITE(LUPRI,'(7X,F15.7,1X,F15.7,1X,F15.7,1X,F15.7)') 1159 & POLDL(2,IDIP,IANG,IFRVAL), 1160 & POLDA(2,IDIP,IANG,IFRVAL), 1161 & POLVL(2,IDIP,IANG,IFRVAL), 1162 & POLVV(2,IDIP,IANG,IFRVAL) 1163 END IF 1164 END DO 1165 END DO 1166 WRITE(LUPRI,'(A,A)') 1167 & '-------------------------------------------------', 1168 & '---------------------' 1169 END DO 1170 WRITE(LUPRI,'(/,A,1P,D15.6)') 'Frequency: ',FRVAL(IFRZER) 1171 WRITE(LUPRI,'(A)') 'POLVL: Velocity' 1172 WRITE(LUPRI,'(A,/,A)') 1173 & ' Ei Bj POLVL', 1174 & '----------------------' 1175 DO IANG = 1,3 1176 DO IDIP = 1,3 1177 WRITE(LUPRI,'(1X,I2,1X,I2,1X,F15.7)') 1178 & IDIP,IANG,POLVL(1,IDIP,IANG,IFRVAL) 1179 IF (ABSORP) THEN 1180 WRITE(LUPRI,'(7X,F15.7)') 1181 & POLVL(2,IDIP,IANG,IFRVAL) 1182 END IF 1183 END DO 1184 END DO 1185 WRITE(LUPRI,'(A)') 1186 & '----------------------' 1187 END IF 1188 END IF 1189C 1190C Skip point if no frequencies specified 1191C STATIC = static approximation to optical rotation G tensor 1192 IF (STATIC) 1193 & CALL LNRAST(POLDLS,POLDAS,WORK(KCMO),WORK(KWORK1),LWORK1) 1194C 1195 999 CONTINUE 1196C 1197C-tbp CALL GPCLOSE(LUSOVE,'DELETE') 1198C-tbp CALL GPCLOSE(LUGDVE,'DELETE') 1199C-tbp CALL GPCLOSE(LUREVE,'DELETE') 1200C 1201 END IF ! for IF (ALFA .OR. ROAA .OR. ROAG) THEN 1202C 1203C Start calculation of vibrational g-factors. 1204C =========================================== 1205C 1206 IF (VIB_G) THEN 1207 CALL HEADER('Calculating terms for vibrational g-factor',-1) 1208C 1209C Set variables and logicals 1210C 1211 CICLC = .FALSE. 1212 HFCLC = NASHT .LE. 1 1213 TRIPLE = .FALSE. 1214 EXECLC = .FALSE. 1215 NABATY = -1 1216 LUSOVE = -1 1217 LUGDVE = -1 1218 LUREVE = -1 1219 NABAOP = 1 1220 BLANK = ' ' 1221 IFRVAL = 1 ! SPAS:8/10-08 there is no frequency for VIB_G 1222C 1223 CALL GPOPEN(LUSOVE,' ','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 1224 CALL GPOPEN(LUGDVE,' ','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 1225 CALL GPOPEN(LUREVE,' ','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 1226C 1227CSPAS:6/10-08: trying to add mass independent vibrational g-factor 1228C CALL DZERO(FOVIBG,3*NATOMS*3*NATOMS*NSYM*NFRVAL) 1229 CALL DZERO(FOVIBG,(3*NATOMS+3)*(3*NATOMS+3)*NSYM*NFRVAL) 1230CKeinSPASmehr 1231C 1232C Loop over the first operators in the linear response function. 1233C These are derivatives with respect to the nuclear cartesian 1234C coordinates. 1235C ============================================================== 1236C 1237 NBR = 0 1238C 1239 DO ISYM = 1,NSYM 1240 IF (VIBGIR .AND. (ISYM .NE. IRVIBG)) CYCLE 1241C 1242 CALL LNRVAR(ISYM,IPRLNR,WORK(KWORK1),LWORK1) 1243C 1244C 3. Work Allocations: 1245C 1246 KGD1 = KWORK1 1247 KWRKG1 = KGD1 1248 LWRKG1 = LWORK - KWRKG1 1249 KSLV = KGD1 + 2*NVARPT 1250 KLAST = KSLV + 2*NVARPT 1251 IF (KLAST.GT.LWORK) CALL STOPIT('LNRABA',' ',KLAST,LWORK) 1252 KWRK = KLAST 1253 LWRK = LWORK - KLAST + 1 1254C 1255 DO IATOM = 1,NUCIND 1256 DO ICOOR = 1,3 1257C 1258C Determine the displacement coordinate number ISCOOR 1259C of symmetry ISYM. If the displacement coordinate 1260C does not have symmetry ISYM, ISCOOR is set to zero. 1261C =================================================== 1262C 1263 ISCOOR = IPTCNT(3*(IATOM - 1) + ICOOR,ISYM-1,1) 1264C 1265 IF (ISCOOR.GT.0) THEN 1266C 1267 NBR = NBR + 1 1268C 1269 IF (NBR .NE. ISCOOR) THEN 1270 WRITE(LUPRI,'(A/A,2I8)') 1271 & 'Keld, why are NBR and ISCOOR different?', 1272 & 'IATOM,NBR,ISCOOR',IATOM,NBR,ISCOOR 1273 END IF 1274C 1275C Find right hand side for first operator 1276C and write to file 1277C ======================================= 1278C 1279 CALL RDS2X(NBR,WORK(KUDV),WORK(KXINDX), 1280 & WORK(KWRKG1),LWRKG1) 1281 REWIND LUGDVE 1282 CALL WRITT(LUGDVE,2*NVARPT,WORK(KWRKG1)) 1283C 1284 IF (IPRLNR.GT.3) THEN 1285 WRITE (LUPRI,'(/A)') 'GP Vector, label: d/dR' 1286 CALL OUTPUT(WORK(KWRKG1),1,NVARPT,1,2, 1287 & NVARPT,2,1,LUPRI) 1288 ENDIF 1289C 1290C Calculate eigenvector for first operator 1291C (i.e. derivatives with respect to the 1292C nuclear cartesian coordinates) and write to file. 1293C ================================================= 1294C 1295 CALL ABARSP(CICLC,HFCLC,TRIPLE,OOTV,ISYM,EXECLC, 1296 & FRVAL,NFRVAL,NABATY,NABAOP,'d/dR ', 1297 & LUGDVE,LUSOVE,LUREVE,THCLNR,MAXITE, 1298 & IPRRSP,MXRM,MXPHP,WORK(KWRK),LWRK) 1299C 1300C Loop over the second operators in the linear 1301C response function. 1302C These are derivatives with respect to the nuclear 1303C cartesian coordinates. 1304C ================================================= 1305C 1306 DO JATOM = 1,NUCIND 1307 DO JCOOR = 1,3 1308C 1309C Determine the displacement coordinate number 1310C JSCOOR of symmetry ISYM. If the displacement 1311C coordinate does not have symmetry ISYM, 1312C JSCOOR is set to zero. 1313C ============================================ 1314C 1315 JSCOOR = IPTCNT(3*(JATOM-1)+JCOOR,ISYM-1,1) 1316C 1317 IF (JSCOOR.GT.0) THEN 1318C 1319C Find right hand side for second operator 1320C ======================================== 1321C 1322 CALL RDS2X(JSCOOR,WORK(KUDV),WORK(KXINDX), 1323 & WORK(KWRKG1),LWRKG1) 1324C 1325 IF (IPRLNR.GT.13) THEN 1326 WRITE (LUPRI,'(/A)') 1327 & 'GP Vector, label: d/dR' 1328 CALL OUTPUT(WORK(KWRKG1),1,2*NVARPT,1, 1329 & 1,2*NVARPT,1,1,LUPRI) 1330 ENDIF 1331 1332C Form the second order property: 1333C vibrational g-factor. 1334C =============================== 1335C 1336 REWIND LUSOVE 1337C 1338CSPAS:8/10-08 there is no frequency for VIB_G 1339C DO IFRVAL = 1,NFRVAL 1340CKeinSPASmehr 1341C 1342 CALL READT(LUSOVE,2*NVARPT,WORK(KSLV)) 1343C 1344 FOVIBG(NBR,JSCOOR,ISYM,IFRVAL) 1345 & = DDOT(2*NVARPT,WORK(KSLV),1, 1346 & WORK(KGD1),1) 1347C 1348 IF (IPRLNR.GT.13) THEN 1349 WRITE (LUPRI,'(2A)') 1350 & 'Solution Vector, label: ','d/dR' 1351 CALL OUTPUT(WORK(KSLV),1,2*NVARPT,1, 1352 & 1,2*NVARPT,1,1,LUPRI) 1353 ENDIF 1354C 1355 IF (IPRLNR.GT.2) THEN 1356CSPAS:8/10-08 there is no frequency for VIB_G 1357C WRITE (LUPRI,'(A,F15.8)') 1358C & 'Frequency = ',FRVAL(IFRVAL) 1359CKeinSPASmehr 1360 WRITE (LUPRI,'(4A,F15.8)') 1361 & 'Second order property for ', 1362 & 'd/dR','d/dR',' = ', 1363 & FOVIBG(NBR,JSCOOR,ISYM,IFRVAL) 1364 ENDIF 1365C 1366CSPAS:8/10-08 there is no frequency for VIB_G 1367C END DO 1368CKeinSPASmehr 1369C 1370 END IF 1371C 1372 END DO 1373 END DO 1374C 1375CSPAS:6/10-08: trying to add mass independent vibrational g-factor 1376C 1377C Loop over dipole velocity operators 1378C =================================== 1379C 1380 DO JPRLBL = 1, NLBTOT 1381C 1382C Find label and symmetry of operator 1383C =================================== 1384C 1385 LABEL2 = LABAPP(JPRLBL) 1386 KSYM = LABSYM(JPRLBL) 1387C 1388C Check for dipole velocity operator 1389C ================================== 1390C 1391 IF (LABEL2(2:7).EQ.'DIPVEL') THEN 1392C 1393 IF (LABEL2(1:1).EQ.'X') JDIPV = 1 1394 IF (LABEL2(1:1).EQ.'Y') JDIPV = 2 1395 IF (LABEL2(1:1).EQ.'Z') JDIPV = 3 1396C 1397C If symmetry of first operator equals symmetry of 1398C second operator, that is if ISYM = KSYM, then 1399C ================================================ 1400C 1401 IF (KSYM.EQ.ISYM) THEN 1402C 1403C Find right hand side (gradient) 1404C for second operator 1405C =============================== 1406C 1407 CALL GETGPV(LABEL2,DUMMY,DUMMY,WORK(KCMO), 1408 & WORK(KUDV),WORK(KPVX), 1409 & WORK(KXINDX),ANTSYM, 1410 & WORK(KWRKG1),LWRKG1) 1411C 1412 IF (IPRLNR.GT.3) THEN 1413 WRITE (LUPRI,'(/2A)') 1414 & 'GP Vector, label: ',LABEL2 1415 CALL OUTPUT(WORK(KWRKG1),1,2*NVARPT,1, 1416 & 1,2*NVARPT,1,1,LUPRI) 1417 ENDIF 1418C 1419C Form the second order property: 1420C << X/Y/ZDIPVEL ; d/dR >> 1421C =============================== 1422C 1423 REWIND LUSOVE 1424C 1425CSPAS:8/10-08 there is no frequency for VIB_G 1426C DO IFRVAL = 1,NFRVAL 1427CKeinSPASmehr 1428C 1429 CALL READT(LUSOVE,2*NVARPT,WORK(KSLV)) 1430C 1431 FOVIBG(NBR,3*NATOMS+JDIPV,ISYM,IFRVAL) 1432 & = DDOT(2*NVARPT,WORK(KSLV),1, 1433 & WORK(KGD1),1) 1434C 1435 IF (IPRLNR.GT.3) THEN 1436 WRITE (LUPRI,'(2A)') 1437 & 'Solution Vector, label: ','d/dR' 1438 CALL OUTPUT(WORK(KSLV),1,2*NVARPT,1, 1439 & 1,2*NVARPT,1,1,LUPRI) 1440 ENDIF 1441C 1442 IF (IPRLNR.GT.2) THEN 1443CSPAS:8/10-08 there is no frequency for VIB_G 1444C WRITE (LUPRI,'(A,F15.8)') 1445C & 'Frequency = ',FRVAL(IFRVAL) 1446CKeinSPASmehr 1447 WRITE (LUPRI,'(4A,F15.8)') 1448 & 'Second order property for ', 1449 & LABEL2,'d/dR',' = ', 1450 & FOVIBG(NBR,3*NATOMS+JDIPV,ISYM, 1451 & IFRVAL) 1452 ENDIF 1453C 1454CSPAS:8/10-08 there is no frequency for VIB_G 1455C END DO 1456CKeinSPASmehr 1457C 1458 ENDIF 1459C 1460 ENDIF 1461C 1462 END DO 1463C 1464CKeinSPASmehr 1465C 1466C 1467 END IF 1468C 1469 END DO 1470 END DO 1471C 1472CSPAS:6/10-08: trying to add mass independent vibrational g-factor 1473C 1474C Loop over dipole velocity operators 1475C =================================== 1476C 1477 DO IPRLBL = 1, NLBTOT 1478C 1479C Find label and symmetry of operator 1480C =================================== 1481C 1482 LABEL1 = LABAPP(IPRLBL) 1483 KSYMOP = LABSYM(IPRLBL) 1484C 1485C Check for dipole velocity operator 1486C ================================== 1487C 1488 IF (LABEL1(2:7).EQ.'DIPVEL') THEN 1489C 1490 IF (LABEL1(1:1).EQ.'X') IDIPV = 1 1491 IF (LABEL1(1:1).EQ.'Y') IDIPV = 2 1492 IF (LABEL1(1:1).EQ.'Z') IDIPV = 3 1493C 1494C If symmetry of first operator is correct 1495C that is if ISYM = KSYMOP, then 1496C ======================================== 1497C 1498 IF (KSYMOP.EQ.ISYM) THEN 1499C 1500C Find right hand side (gradient) for first operator 1501C ================================================== 1502C 1503 CALL GETGPV(LABEL1,DUMMY,DUMMY,WORK(KCMO), 1504 & WORK(KUDV),WORK(KPVX),WORK(KXINDX), 1505 & ANTSYM,WORK(KWRKG1),LWRKG1) 1506 REWIND LUGDVE 1507 CALL WRITT(LUGDVE,2*NVARPT,WORK(KWRKG1)) 1508C 1509 IF (IPRLNR.GT.3) THEN 1510 WRITE (LUPRI,'(/2A)') 1511 & 'GP Vector, label: ',LABEL1 1512 CALL OUTPUT(WORK(KWRKG1),1,2*NVARPT,1, 1513 & 1,2*NVARPT,1,1,LUPRI) 1514 ENDIF 1515C 1516C Calculate linear response vector and write to file 1517C ================================================== 1518C 1519 CALL ABARSP(CICLC,HFCLC,TRIPLE,OOTV,ISYM,EXECLC, 1520 & FRVAL,NFRVAL,NABATY,NABAOP,LABEL1,LUGDVE, 1521 & LUSOVE,LUREVE,THCLNR,MAXITE,IPRRSP,MXRM,MXPHP, 1522 & WORK(KWRK),LWRK) 1523C 1524C Loop over the second operators in the linear 1525C response function. 1526C These are derivatives with respect to the nuclear 1527C cartesian coordinates. 1528C ================================================= 1529C 1530 DO JATOM = 1,NUCIND 1531 DO JCOOR = 1,3 1532C 1533C Determine the displacement coordinate number 1534C JSCOOR of symmetry ISYM. If the displacement 1535C coordinate does not have symmetry ISYM, 1536C JSCOOR is set to zero. 1537C ============================================ 1538C 1539 JSCOOR = IPTCNT(3*(JATOM-1)+JCOOR,ISYM-1,1) 1540C 1541 IF (JSCOOR.GT.0) THEN 1542C 1543C Find right hand side for second operator 1544C ======================================== 1545C 1546 CALL RDS2X(JSCOOR,WORK(KUDV),WORK(KXINDX), 1547 & WORK(KWRKG1),LWRKG1) 1548C 1549 IF (IPRLNR.GT.3) THEN 1550 WRITE (LUPRI,'(/A)') 1551 & 'GP Vector, label: d/dR' 1552 CALL OUTPUT(WORK(KWRKG1),1,2*NVARPT,1, 1553 & 1,2*NVARPT,1,1,LUPRI) 1554 ENDIF 1555 1556C Form the second order property: 1557C << d/dR ; X/Y/ZDIPVEL >> 1558C =============================== 1559C 1560 REWIND LUSOVE 1561C 1562 CALL READT(LUSOVE,2*NVARPT,WORK(KSLV)) 1563C 1564 FOVIBG(3*NATOMS+IDIPV,JSCOOR,ISYM,IFRVAL) 1565 & = DDOT(2*NVARPT,WORK(KSLV),1, 1566 & WORK(KGD1),1) 1567C 1568 IF (IPRLNR.GT.3) THEN 1569 WRITE (LUPRI,'(2A)') 1570 & 'Solution Vector, label: ',LABEL1 1571 CALL OUTPUT(WORK(KSLV),1,2*NVARPT,1, 1572 & 1,2*NVARPT,1,1,LUPRI) 1573 ENDIF 1574C 1575 IF (IPRLNR.GT.2) THEN 1576 WRITE (LUPRI,'(4A,F15.8)') 1577 & 'Second order property for ', 1578 & 'd/dR',LABEL1,' = ', 1579 & FOVIBG(3*NATOMS+IDIPV,JSCOOR,ISYM, 1580 & IFRVAL) 1581 ENDIF 1582C 1583 END IF 1584C 1585 END DO 1586 END DO 1587C 1588C Loop over dipole velocity operators 1589C =================================== 1590C 1591 DO JPRLBL = 1, NLBTOT 1592C 1593C Find label and symmetry of operator 1594C =================================== 1595C 1596 LABEL2 = LABAPP(JPRLBL) 1597 KSYM = LABSYM(JPRLBL) 1598C 1599C Check for dipole velocity operator 1600C ================================== 1601C 1602 IF (LABEL2(2:7).EQ.'DIPVEL') THEN 1603C 1604 IF (LABEL2(1:1).EQ.'X') JDIPV = 1 1605 IF (LABEL2(1:1).EQ.'Y') JDIPV = 2 1606 IF (LABEL2(1:1).EQ.'Z') JDIPV = 3 1607C 1608C If symmetry of first operator equals symmetry of 1609C second operator, that is if ISYM = KSYM, then 1610C ================================================ 1611C 1612 IF (KSYM.EQ.ISYM) THEN 1613C 1614C Find right hand side (gradient) 1615C for second operator 1616C =============================== 1617C 1618 CALL GETGPV(LABEL2,DUMMY,DUMMY,WORK(KCMO), 1619 & WORK(KUDV),WORK(KPVX), 1620 & WORK(KXINDX),ANTSYM, 1621 & WORK(KWRKG1),LWRKG1) 1622C 1623 IF (IPRLNR.GT.3) THEN 1624 WRITE (LUPRI,'(/2A)') 1625 & 'GP Vector, label: ',LABEL2 1626 CALL OUTPUT(WORK(KWRKG1),1,2*NVARPT,1, 1627 & 1,2*NVARPT,1,1,LUPRI) 1628 ENDIF 1629C 1630C Form the second order property: 1631C << X/Y/ZDIPVEL ; X/Y/ZDIPVEL >> 1632C =============================== 1633C 1634 REWIND LUSOVE 1635C 1636 CALL READT(LUSOVE,2*NVARPT,WORK(KSLV)) 1637C 1638 FOVIBG(3*NATOMS+IDIPV,3*NATOMS+JDIPV,ISYM, 1639 & IFRVAL) 1640 & = DDOT(2*NVARPT,WORK(KSLV),1, 1641 & WORK(KGD1),1) 1642C 1643 IF (IPRLNR.GT.3) THEN 1644 WRITE (LUPRI,'(2A)') 1645 & 'Solution Vector, label: ',LABEL1 1646 CALL OUTPUT(WORK(KSLV),1,2*NVARPT,1, 1647 & 1,2*NVARPT,1,1,LUPRI) 1648 ENDIF 1649C 1650 IF (IPRLNR.GT.2) THEN 1651 WRITE (LUPRI,'(4A,F15.8)') 1652 & 'Second order property for ', 1653 & LABEL2,LABEL1,' = ', 1654 & FOVIBG(3*NATOMS+IDIPV, 1655 & 3*NATOMS+JDIPV,ISYM, 1656 & IFRVAL) 1657 ENDIF 1658C 1659 ENDIF 1660C 1661 ENDIF 1662C 1663 END DO 1664C 1665 ENDIF 1666C 1667 ENDIF 1668C 1669 END DO 1670C 1671CKeinSPASmehr 1672C 1673 END DO 1674C 1675 CALL GPCLOSE(LUSOVE,'DELETE') 1676 CALL GPCLOSE(LUGDVE,'DELETE') 1677 CALL GPCLOSE(LUREVE,'DELETE') 1678C 1679 END IF 1680C 1681 CALL TIMER ('LNRABA',TIMEIN,TIMOUT) 1682 PASS = .TRUE. 1683 IF (CUT) THEN 1684 WRITE (LUPRI,'(/A/A)') 1685 & ' Program stopped after LNRABA as requested.', 1686 & ' No restart file has been written.' 1687 CALL QUIT(' ***** End of ABACUS (.STOP in *ABALNR) *****') 1688 END IF 1689C 1690 CALL QEXIT('LNRABA') 1691 RETURN 1692 END 1693C 1694C /* Deck lnrast */ 1695 SUBROUTINE LNRAST(POLDLS,POLDAS,CMO,WORK,LWORK) 1696C 1697C _ST_atic approximation to optical rotation G tensors. 1698C 1699#include "implicit.h" 1700#include "dummy.h" 1701#include "mxcent.h" 1702#include "maxorb.h" 1703#include "iratdef.h" 1704#include "priunit.h" 1705#include "cbilnr.h" 1706 LOGICAL FOUND, CONV 1707 DIMENSION WORK(LWORK), SNDPRP(2), CMO(*) 1708 DIMENSION POLDLS(3,3), POLDAS(3,3) 1709 DIMENSION DTOT(N2ORBX), DTOTSV(N2ORBX) 1710 DIMENSION DIPSLV(N2ORBX), ANGSLV(N2ORBX) 1711 CHARACTER*8 LABEL1, LABEL2, BLANK 1712 PARAMETER (DP5=0.5D0) 1713#include "cbiexc.h" 1714#include "inflin.h" 1715#include "infvar.h" 1716#include "infdim.h" 1717#include "inforb.h" 1718#include "nuclei.h" 1719#include "inftap.h" 1720#include "infrsp.h" 1721#include "wrkrsp.h" 1722#include "maxmom.h" 1723#include "maxaqn.h" 1724#include "symmet.h" 1725#include "abainf.h" 1726#include "gnrinf.h" 1727#include "infsop.h" 1728#include "absorp.h" 1729 BLANK = ' ' 1730C 1731C Construct one-electron density matrix 1732C 1733 CALL MKDENS(DTOT,WORK,LWORK,IPRINT) 1734C 1735 DO ISYM=1,NSYM 1736 KSYMOP = ISYM 1737 DO IOPER=1,NOPER(KSYMOP) 1738 KSLV1 = 1 1739 KSLV2 = KSLV1 + 2*NVARPT 1740 LABEL1 = LABOP(IOPER,KSYMOP) 1741 IF (LABEL1(2:7).EQ.'DIPLEN') THEN 1742 IF (LABEL1 .EQ.'XDIPLEN') IDIP = 1 1743 IF (LABEL1 .EQ.'YDIPLEN') IDIP = 2 1744 IF (LABEL1 .EQ.'ZDIPLEN') IDIP = 3 1745 CALL GPOPEN(LURSP,'RSPVEC','OLD',' ', 1746 & ' ',IDUMMY,.FALSE.) 1747 CALL REARSP(LURSP,2*NVARPT,WORK(KSLV1),LABEL1, 1748 & BLANK,0.0D0,0.0D0,ISYM,0, 1749 & THCLNR,FOUND,CONV,ANTSYM) 1750 CALL GPCLOSE(LURSP,'KEEP') 1751 IF (.NOT.FOUND) THEN 1752 WRITE (LUPRI,'(2A)') 'Solution Vector, label: ' 1753 & ,LABEL1,' not found on RSPVEC' 1754 CALL QUIT('Solution vector 1 not found on RSPVEC') 1755 ELSE 1756 IF (IPRLNR.GT.3) THEN 1757 WRITE(LUPRI,'(2A)') 'Solution Vector, '// 1758 & 'label: ',LABEL1 1759 CALL OUTPUT(WORK(KSLV1),1,NVARPT,1,2,NVARPT, 1760 & 2,1,LUPRI) 1761 END IF 1762 CALL UPWOP(WORK(KSLV1),DIPSLV,.TRUE.) 1763 CALL DGEMM('N','T',NORBT,NORBT,NORBT,1.D0, 1764 & DTOT,NORBT,DIPSLV,NORBT,0.D0,DTOTSV,NORBT) 1765 DO IOPER2=1,NOPER(KSYMOP) 1766 LABEL2 = LABOP(IOPER2,KSYMOP) 1767 IF (LABEL2(2:7).EQ.'ANGMOM') THEN 1768 CALL GPOPEN(LURSP,'RSPVEC','OLD',' ', 1769 & ' ',IDUMMY,.FALSE.) 1770 CALL REARSP(LURSP,2*NVARPT,WORK(KSLV2), 1771 & LABEL2,BLANK,0.0D0,0.0D0,ISYM,0, 1772 & THCLNR,FOUND,CONV,ANTSYM) 1773 CALL GPCLOSE(LURSP,'KEEP') 1774 IF (.NOT.FOUND) THEN 1775 WRITE (LUPRI,'(2A)') 'Solution Vector, '// 1776 & 'label: ',LABEL2,' not found on RSPVEC' 1777 CALL QUIT('Solution vector 2 not found '// 1778 & 'on RSPVEC') 1779 ELSE 1780 IF (IPRLNR.GT.3) THEN 1781 WRITE(LUPRI,'(2A)') 'Solution '// 1782 & 'Vector, label: ', LABEL2 1783 CALL OUTPUT(WORK(KSLV2),1,NVARPT,1,2, 1784 & NVARPT,2,1,LUPRI) 1785 END IF 1786 CALL UPWOP(WORK(KSLV2),ANGSLV,.FALSE.) 1787 SNDPRP(1)=DDOT(N2ORBX,DTOTSV,1, 1788 & ANGSLV,1) 1789 IF (LABEL2(1:1).EQ.'X') 1790 & POLDAS(IDIP,1) = SNDPRP(1) 1791 IF (LABEL2(1:1).EQ.'Y') 1792 & POLDAS(IDIP,2) = SNDPRP(1) 1793 IF (LABEL2(1:1).EQ.'Z') 1794 & POLDAS(IDIP,3) = SNDPRP(1) 1795 END IF 1796 ELSE IF (LABEL2(2:7).EQ.'LONMAG') THEN 1797 CALL GPOPEN(LURSP,'RSPVEC','OLD',' ', 1798 & ' ',IDUMMY,.FALSE.) 1799 CALL REARSP(LURSP,2*NVARPT,WORK(KSLV2), 1800 & LABEL2,BLANK,0.0D0,0.0D0,ISYM,0, 1801 & THCLNR,FOUND,CONV,ANTSYM) 1802 CALL GPCLOSE(LURSP,'KEEP') 1803 IF (.NOT.FOUND) THEN 1804 WRITE (LUPRI,'(2A)') 'Solution Vector, '// 1805 & 'label: ',LABEL2,' not found on RSPVEC' 1806 CALL QUIT('Solution vector 2 not found '// 1807 & 'on RSPVEC') 1808 ELSE 1809 IF (IPRLNR.GT.3) THEN 1810 WRITE(LUPRI,'(2A)') 'Solution '// 1811 & 'Vector, label: ',LABEL2 1812 CALL OUTPUT(WORK(KSLV2),1,NVARPT,1,2, 1813 & NVARPT,2,1,LUPRI) 1814 END IF 1815 CALL UPWOP(WORK(KSLV2),ANGSLV,.FALSE.) 1816 SNDPRP(1)=DDOT(N2ORBX,DTOTSV,1, 1817 & ANGSLV,1) 1818 IF (LABEL2(1:1).EQ.'X') 1819 & POLDLS(IDIP,1) = SNDPRP(1) 1820 IF (LABEL2(1:1).EQ.'Y') 1821 & POLDLS(IDIP,2) = SNDPRP(1) 1822 IF (LABEL2(1:1).EQ.'Z') 1823 & POLDLS(IDIP,3) = SNDPRP(1) 1824 END IF 1825 END IF 1826 END DO 1827 END IF 1828 END IF 1829 END DO 1830 END DO 1831 WRITE(LUPRI,9003) 1832 WRITE(LUPRI,9007) 1833 WRITE(LUPRI,9009) 'Bx','By','Bz' 1834 WRITE(LUPRI,9008) 'Ex',(POLDLS(1,JDIP), JDIP=1,3) 1835 WRITE(LUPRI,9008) 'Ey',(POLDLS(2,JDIP), JDIP=1,3) 1836 WRITE(LUPRI,9008) 'Ez',(POLDLS(3,JDIP), JDIP=1,3) 1837C 1838 WRITE(LUPRI,9004) 1839 WRITE(LUPRI,9007) 1840 WRITE(LUPRI,9009) 'Bx','By','Bz' 1841 WRITE(LUPRI,9008) 'Ex',(POLDAS(1,JDIP), JDIP=1,3) 1842 WRITE(LUPRI,9008) 'Ey',(POLDAS(2,JDIP), JDIP=1,3) 1843 WRITE(LUPRI,9008) 'Ez',(POLDAS(3,JDIP), JDIP=1,3) 1844 DO IFRVAL=1,NFRVAL 1845 WRITE(LUPRI,9005) FRVAL(IFRVAL) 1846 WRITE(LUPRI,9007) 1847 WRITE(LUPRI,9009) 'Bx','By','Bz' 1848 WRITE(LUPRI,9008) 'Ex',(POLDLS(1,JDIP)*FRVAL(IFRVAL), 1849 & JDIP=1,3) 1850 WRITE(LUPRI,9008) 'Ey',(POLDLS(2,JDIP)*FRVAL(IFRVAL), 1851 & JDIP=1,3) 1852 WRITE(LUPRI,9008) 'Ez',(POLDLS(3,JDIP)*FRVAL(IFRVAL), 1853 & JDIP=1,3) 1854 1855 WRITE(LUPRI,9006) FRVAL(IFRVAL) 1856 WRITE(LUPRI,9007) 1857 WRITE(LUPRI,9009) 'Bx','By','Bz' 1858 WRITE(LUPRI,9008) 'Ex',(POLDAS(1,JDIP)*FRVAL(IFRVAL), 1859 & JDIP=1,3) 1860 WRITE(LUPRI,9008) 'Ey',(POLDAS(2,JDIP)*FRVAL(IFRVAL), 1861 & JDIP=1,3) 1862 WRITE(LUPRI,9008) 'Ez',(POLDAS(3,JDIP)*FRVAL(IFRVAL), 1863 & JDIP=1,3) 1864 END DO 1865 9003 FORMAT (//10X,'London G tensor using static approximation') 1866 9004 FORMAT (//10X,'No-London G tensor using static approximation') 1867 9005 FORMAT (//10X,'London G tensor stat. appr. for freq ',F12.6, 1868 & ' au') 1869 9006 FORMAT (/10X,'No-London G tensor stat. appr. for freq ',F12.6, 1870 & ' au') 1871 9007 FORMAT (10X,'------------------------------------------------') 1872 9008 FORMAT (16X,A,2X,3F12.7) 1873 9009 FORMAT (27X,3(A,10X)/) 1874C 1875 RETURN 1876 END 1877C 1878C /* Deck lnrvar */ 1879 SUBROUTINE LNRVAR(ISYM,IPRINT,WORK,LWORK) 1880C 1881#include "implicit.h" 1882#include "mxcent.h" 1883#include "maxorb.h" 1884 DIMENSION WORK(LWORK) 1885#include "inflin.h" 1886#include "infvar.h" 1887#include "infrsp.h" 1888#include "inforb.h" 1889#include "wrkrsp.h" 1890#include "abainf.h" 1891C 1892 LOGICAL TRIPLE 1893C 1894C Set variables for response modules 1895C ================================== 1896C 1897 TRIPLE = .FALSE. 1898 CALL ABAVAR(ISYM,TRIPLE,IPRINT,WORK,LWORK) 1899 IREFSY = LSYMRF 1900 NCREF = NCONRF 1901 KSYMST = LSYMST 1902 KSYMOP = LSYMPT 1903 KZWOPT = NWOPT 1904 IF ((NASHT .EQ. 0).AND.(.NOT.ABASOP)) THEN 1905 KZCONF = 0 1906 ELSE IF (NASHT .EQ. 1) THEN 1907 KZCONF = 0 ! also zero for ROHF 1908 ELSE 1909 KZCONF = NCONST 1910 END IF 1911 KZVAR = KZWOPT + KZCONF 1912 KZYWOP = 2*KZWOPT 1913 KZYCON = 2*KZCONF 1914 KZYVAR = 2*KZVAR 1915C 1916 RETURN 1917 END 1918C /* Deck betagm */ 1919 FUNCTION BETAGM(ALFA,GM) 1920C 1921C Calculate Beta(G')**2 = 1922C BETAGM = 0.5*(3*ALFA(I,J)*GM(I,J) - ALFA(I,I)*GM(J,J)) 1923C 1924#include "implicit.h" 1925 PARAMETER ( D0 = 0.0D0 , DP5 = 0.5D0 , D3 = 3.0D0 ) 1926 DIMENSION ALFA(3,3),GM(3,3) 1927C 1928 BETAGM = D0 1929 DO 100 I = 1,3 1930 DO 200 J = 1,3 1931 BETAGM = BETAGM + 1932 & DP5*(D3*ALFA(I,J)*GM(I,J)-ALFA(I,I)*GM(J,J)) 1933 200 CONTINUE 1934 100 CONTINUE 1935C 1936 RETURN 1937 END 1938C /* Deck betaal */ 1939 FUNCTION BETAAL(ALFA) 1940C 1941C Calculate Beta(alfa)**2 = 1942C BETAAL = 0.5 * (3*ALFA(I,J)*ALFA(I,J) - ALFA(I,I)*ALFA(J,J)) 1943C 1944#include "implicit.h" 1945 PARAMETER ( D0 = 0.0D0 , DP5 = 0.5D0 , D3 = 3.0D0 ) 1946 DIMENSION ALFA(3,3) 1947C 1948 BETAAL = D0 1949 DO 100 I = 1,3 1950 DO 200 J = 1,3 1951 BETAAL = BETAAL + DP5*(D3*ALFA(I,J)*ALFA(I,J) 1952 * -ALFA(I,I)*ALFA(J,J)) 1953 200 CONTINUE 1954 100 CONTINUE 1955C 1956 RETURN 1957 END 1958C /* Deck betaa */ 1959 FUNCTION BETAA(ALFA,A,OMEGA) 1960C 1961C Calculate Beta(A)**2 = 1962C BETAA = 0.5*OMEGA*ALFA(I,J)*EPSILON(I,K,L)*A(K,L,J) 1963C 1964#include "implicit.h" 1965 PARAMETER ( D0 = 0.0D0 , DP5 = 0.5D0 , D1 = 1.0D0 , DM1 = -1.0D0) 1966 DIMENSION ALFA(3,3), A(3,3,3) 1967 DIMENSION EPSILO(3,3,3) 1968 DATA EPSILO / D0, D0, D0, D0, D0, DM1, D0, D1, D0, 1969 & D0, D0, D1, D0, D0, D0, DM1, D0, D0, 1970 & D0, DM1, D0, D1, D0, D0, D0, D0, D0/ 1971C 1972 BETAA = D0 1973 DO 100 I = 1,3 1974 DO 200 J = 1,3 1975 DO 300 K = 1,3 1976 DO 400 L = 1,3 1977 BETAA = BETAA + DP5*OMEGA*ALFA(I,J)*EPSILO(I,K,L)* 1978 * A(K,L,J) 1979 400 CONTINUE 1980 300 CONTINUE 1981 200 CONTINUE 1982 100 CONTINUE 1983C 1984 RETURN 1985 END 1986C /* Deck alfmn */ 1987 FUNCTION ALFMN(ALFA) 1988C 1989C Calculate AlfaMean = ALFMN = ( (1/3) * ALFA(I,I) ) 1990C 1991#include "implicit.h" 1992 PARAMETER ( D0 = 0.0D0, D3 = 3.0D0 ) 1993 DIMENSION ALFA(3,3) 1994C 1995 ALFMN = D0 1996 DO 100 I = 1,3 1997 ALFMN = ALFMN + ALFA(I,I) 1998 100 CONTINUE 1999C 2000 ALFMN = ALFMN / D3 2001C 2002 RETURN 2003 END 2004C /* Deck gmmn */ 2005 FUNCTION GMMN(GM) 2006C 2007C Calculate G' = GMMN = ( (1/3) * GM(I,I) ) 2008C 2009#include "implicit.h" 2010 PARAMETER ( D0 = 0.0D0, D3 = 3.0D0 ) 2011 DIMENSION GM(3,3) 2012C 2013 GMMN = D0 2014 DO 100 I = 1,3 2015 GMMN = GMMN + GM(I,I) 2016 100 CONTINUE 2017C 2018 GMMN = GMMN / D3 2019C 2020 RETURN 2021 END 2022C /* Deck raminl */ 2023 FUNCTION RAMINL(ALFMN,BETAAL) 2024C 2025C Calculate RamanIntensity = RAMIN = (45 * ALFMN**2) + (4 * BETA**2) 2026C for linear pol. incident radiation perpendicular to scattering plane 2027C 2028#include "implicit.h" 2029 PARAMETER ( D0 = 0.0D0, D4 = 4.0D0, D45 = 45.0D0 ) 2030C 2031 RAMINL = (D45 * ALFMN**2) + (D4 * BETAAL) 2032C 2033 RETURN 2034 END 2035C /* Deck depoll */ 2036 FUNCTION DEPOLL(ALFMN,BETAAL) 2037C 2038C Calculate the Depolarization Ratio for rightangle scattering, 2039C linear polarized incident light, perpendicular (=parall/perpend) 2040C plane = DEPOLR = (3 * BETA**2) / (45 * ALFMN**2 + 4 * BETA**2) 2041C 2042#include "implicit.h" 2043 PARAMETER ( D0 = 0.0D0, D3 = 3.0D0, D4 = 4.0D0, D45 = 45.0D0 ) 2044C 2045 A = D3 * BETAAL 2046 B = (D45 * ALFMN**2) + (D4 * BETAAL) 2047C 2048 IF ( B .GT. 0.D-6 ) THEN 2049 DEPOLL = A / B 2050 ELSE 2051 DEPOLL = D0 2052 ENDIF 2053C 2054 RETURN 2055 END 2056C /* Deck raminn */ 2057 FUNCTION RAMINN(ALFMN,BETAAL) 2058C 2059C Calculate RamanIntensity = RAMIN = (45 * ALFMN**2) + (7 * BETA**2) 2060C for natural & circular pol. incident radiation 2061C 2062#include "implicit.h" 2063 PARAMETER ( D0 = 0.0D0, D7 = 7.0D0, D45 = 45.0D0 ) 2064C 2065 RAMINN = (D45 * ALFMN**2) + (D7 * BETAAL) 2066C 2067 RETURN 2068 END 2069C /* Deck depoln */ 2070 FUNCTION DEPOLN(ALFMN,BETAAL) 2071C 2072C Calculate the Depolarization Ratio for rightangle scattering, 2073C naturel or circ. pol. incident light (=parall/perpend). 2074C plane = DEPOLR = (6 * BETA**2) / (45 * ALFMN**2 + 7 * BETA**2) 2075C 2076#include "implicit.h" 2077 PARAMETER ( D0 = 0.0D0, D6 = 6.0D0, D7 = 7.0D0, D45 = 45.0D0 ) 2078C 2079 A = D6 * BETAAL 2080 B = (D45 * ALFMN**2) + (D7 * BETAAL) 2081C 2082 IF ( B .GT. 0.D-6 ) THEN 2083 DEPOLN = A / B 2084 ELSE 2085 DEPOLN = D0 2086 ENDIF 2087C 2088 RETURN 2089 END 2090C /* Deck deltaz */ 2091 FUNCTION DELTAZ(BETAGM,BETAA) 2092C 2093C Calculate the difference differential scattering cross section 2094C in depolarized rightangle scattering = 2095C DELTAZ = 4/c ( 6*BETAGM - 2*BETAA) 2096C 2097#include "implicit.h" 2098#include "codata.h" 2099 PARAMETER ( D2 = 2.0D0, D4 = 4.0D0, D6 = 6.0D0 ) 2100C 2101 DELTAZ = D4 / CVEL * (D6 * BETAGM - D2 * BETAA) 2102C 2103 RETURN 2104 END 2105C /* Deck deltax */ 2106 FUNCTION DELTAX(BETAGM,BETAA,ALFAMN,GMMN) 2107C 2108C Calculate the difference differential scattering cross section 2109C in polarized rightangle scattering = 2110C DELTAX = 4/c ( 45*ALFAMN*GMMN+7*BETAGM+BETAA) 2111C 2112#include "implicit.h" 2113#include "codata.h" 2114 PARAMETER ( D4 = 4.0D0, D7 = 7.0D0, D45 = 45.0D0) 2115C 2116 DELTAX = D4 / CVEL * (D45*ALFAMN*GMMN+D7*BETAGM+BETAA) 2117C 2118 RETURN 2119 END 2120C /* Deck delta0 */ 2121 FUNCTION DELTA0(BETAGM,BETAA,ALFAMN,GMMN) 2122C 2123C Calculate the difference differential scattering cross section 2124C in forward scattering with no analyzer = 2125C DELTA0 = 4/c (180*ALFAMN*GMMN+4*BETAGM-4*BETAA) 2126C 2127#include "implicit.h" 2128#include "codata.h" 2129 PARAMETER ( D4 = 4.0D0, D180 = 18.0D1 ) 2130C 2131 DELTA0 = D4 / CVEL * (D180*ALFAMN*GMMN+D4*BETAGM-D4*BETAA) 2132C 2133 RETURN 2134 END 2135C /* Deck deltab */ 2136 FUNCTION DELTAB(BETAGM,BETAA) 2137C 2138C Calculate the difference differential scattering cross section 2139C in backward scattering with no analyzer = 2140C DELTAB = 4/c ( 24*BETAGM+8*BETAA) 2141C 2142#include "implicit.h" 2143#include "codata.h" 2144 PARAMETER ( D4 = 4.0D0, D8 = 8.0D0, D24 = 24.0D0 ) 2145C 2146 DELTAB = D4 / CVEL * (D24*BETAGM + D8*BETAA) 2147C 2148 RETURN 2149 END 2150C /* Deck cid */ 2151 FUNCTION CID(DELTA,RMINN) 2152C 2153C Calculate the Circular Intensity Difference (CID) for all 2154C arragements. CID = - DELTA / (2*RMINN) 2155C RMIN(0) and (180) with no analyzer = 2*RMIN. 2156C Rmin (depolarized) = DEPLN * RMINN 2157#include "implicit.h" 2158C 2159 PARAMETER ( D0 = 0.0D0, DP5 = 0.5D0 ) 2160C 2161 IF ( RMINN .GT. 0.D-6 ) THEN 2162 CID = DELTA / RMINN * DP5 2163 ELSE 2164 CID = D0 2165 ENDIF 2166C 2167 RETURN 2168 END 2169C /* Deck boltzk */ 2170 FUNCTION BOLTZK(FREQ) 2171C 2172C Calculate the Boltzmannfactor for Spectra simulation 2173C BOLTZK = 1/(1-EXP(-h*freq/k/T)), FREQ is in Hartree = 2Pi*freq 2174C 2175#include "implicit.h" 2176#include "codata.h" 2177 PARAMETER ( TSTAND = 298.15D0 ) 2178C 2179 BOLTZK = 1-EXP(-(FREQ*AUTK)/TSTAND) 2180 RETURN 2181 END 2182C /* Deck crossk */ 2183 FUNCTION CROSSK(FRVAL,FREQ) 2184C 2185C Calculate the Constant for the differential scattering cross section 2186C K = 16/90*Pi**4*((freq0-freq)/C)**3*freq0/C 2187C Combined with the freq. factor out of PLACZEK approximation 2188C B = SQRT(1/(4Pi*freq)) 2189C CROSSK = K*B**2 2190C ****************** FREQ is in Hartree = 2Pi*freq 2191C 2192#include "implicit.h" 2193#include "codata.h" 2194 PARAMETER ( D180 = 180.D0) 2195C 2196 CROSSK = (FRVAL-FREQ)**3*FRVAL/FREQ/D180/CVEL**3 2197 RETURN 2198 END 2199C /* Deck lnrout */ 2200 SUBROUTINE LNROUT(POLDD,POLDL,POLDA,POLVL,POLVV,FOVIBG,IPRINT, 2201 & WORK,LWORK) 2202C 2203#include "implicit.h" 2204#include "priunit.h" 2205#include "absorp.h" 2206#include "abslrs.h" 2207#include "inforb.h" 2208#include "mxcent.h" 2209#include "maxaqn.h" 2210#include "maxorb.h" 2211#include "codata.h" 2212 PARAMETER (FACTOR = (288.0D-30)*(PI**2)*XFMOL*(XTANG**4), 2213 & D3 = 3.0D0) 2214#include "nuclei.h" 2215#include "symmet.h" 2216#include "abainf.h" 2217#include "cbilnr.h" 2218 DIMENSION POLDD(2,3,3,NFRVAL), POLDL(2,3,3,NFRVAL), 2219 & POLDA(2,3,3,NFRVAL), 2220 & POLVL(2,3,3,NFRVAL), POLVV(2,3,3,NFRVAL), 2221 & BETANL(2), BETAL(2), ALPHNL(2), ALPHAL(2), 2222 & BETAVL(2), BETAVV(2), ALPHVL(2), ALPHVV(2) 2223cmbh: helper strings 2224 character mlab1*8,mlab2*8,mlab3*8 2225cmbh end 2226 CHARACTER COORDI*21,COORDJ*21 2227CSPAS:6/10-08: trying to add mass independent vibrational g-factor 2228C DIMENSION FOVIBG(3*NATOMS,3*NATOMS,NSYM,MXFR) 2229 DIMENSION FOVIBG(3*NATOMS+3,3*NATOMS+3,NSYM,MXFR) 2230CKeinSPASmehr 2231 DIMENSION WORK(LWORK) 2232 2233C 2234 IF (ALFA .AND. IPRINT .GE. 0) THEN 2235 CALL TITLER('Frequency dependent polarizabilities','+',118) 2236 DO 100 IFRVAL = 1,NFRVAL 2237 WRITE(LUPRI,9001) FRVAL(IFRVAL) 2238 WRITE(LUPRI,9002) 2239 IF (ABSORP) WRITE(LUPRI,'(/,10X,A)') 'Real part:' 2240 WRITE(LUPRI,9003) 'Ex','Ey','Ez' 2241 WRITE(LUPRI,9004) 'Ex',(POLDD(1,1,JDIP,IFRVAL), JDIP=1,3) 2242 WRITE(LUPRI,9004) 'Ey',(POLDD(1,2,JDIP,IFRVAL), JDIP=1,3) 2243 WRITE(LUPRI,9004) 'Ez',(POLDD(1,3,JDIP,IFRVAL), JDIP=1,3) 2244cmbh: print out polarizabilities for MidasCpp 2245 mlab1='XDIPLEN' 2246 mlab2='YDIPLEN' 2247 mlab3='ZDIPLEN' 2248 call wripro(POLDD(1,1,1,IFRVAL),"SCF/DFT ",2, 2249 * mlab1,mlab1,mlab1,mlab1, 2250 * FRVAL(IFRVAL),FRVAL(IFRVAL),FRVAL(IFRVAL),1,0, 2251 * 0,0) 2252 call wripro(POLDD(1,1,2,IFRVAL),"SCF/DFT ",2, 2253 * mlab1,mlab2,mlab1,mlab2, 2254 * FRVAL(IFRVAL),FRVAL(IFRVAL),FRVAL(IFRVAL),1,0, 2255 * 0,0) 2256 call wripro(POLDD(1,1,3,IFRVAL),"SCF/DFT ",2, 2257 * mlab1,mlab3,mlab1,mlab3, 2258 * FRVAL(IFRVAL),FRVAL(IFRVAL),FRVAL(IFRVAL),1,0, 2259 * 0,0) 2260 call wripro(POLDD(1,2,2,IFRVAL),"SCF/DFT ",2, 2261 * mlab2,mlab2,mlab2,mlab2, 2262 * FRVAL(IFRVAL),FRVAL(IFRVAL),FRVAL(IFRVAL),1,0, 2263 * 0,0) 2264 call wripro(POLDD(1,2,3,IFRVAL),"SCF/DFT ",2, 2265 * mlab2,mlab3,mlab2,mlab3, 2266 * FRVAL(IFRVAL),FRVAL(IFRVAL),FRVAL(IFRVAL),1,0, 2267 * 0,0) 2268 call wripro(POLDD(1,3,3,IFRVAL),"SCF/DFT ",2, 2269 * mlab3,mlab3,mlab3,mlab3, 2270 * FRVAL(IFRVAL),FRVAL(IFRVAL),FRVAL(IFRVAL),1,0, 2271 * 0,0) 2272 2273cmbh end 2274 IF (ABSORP) THEN 2275 WRITE(LUPRI,'(/,10X,A)') 'Imaginary part:' 2276 WRITE(LUPRI,9003) 'Ex','Ey','Ez' 2277 WRITE(LUPRI,9004) 'Ex',(POLDD(2,1,JDIP,IFRVAL), JDIP=1,3) 2278 WRITE(LUPRI,9004) 'Ey',(POLDD(2,2,JDIP,IFRVAL), JDIP=1,3) 2279 WRITE(LUPRI,9004) 'Ez',(POLDD(2,3,JDIP,IFRVAL), JDIP=1,3) 2280 END IF 2281 POLISO = (POLDD(1,1,1,IFRVAL) + 2282 & POLDD(1,2,2,IFRVAL) + 2283 & POLDD(1,3,3,IFRVAL)) / D3 2284 WRITE (LUPRI,'(/A,1P,G14.7)') 2285 & 'Isotropic polarizability: ',POLISO 2286 100 CONTINUE 2287 ENDIF 2288 IF (ROAG) THEN 2289 IF (ABSORP) THEN 2290 CALL TITLER('Optical Rotation and Electronic Circular '// 2291 & 'Dichroism','+',102) 2292 ELSE 2293 CALL TITLER('Optical rotation','+',118) 2294 END IF 2295 IF (IPRINT .GE. 0) THEN 2296 DO IFRVAL = 1,NFRVAL 2297 WRITE(LUPRI,9005) FRVAL(IFRVAL) 2298 WRITE(LUPRI,9007) 2299 IF (ABSORP) WRITE(LUPRI,'(/10X,A)') 'Real part:' 2300 WRITE(LUPRI,9009) 'Bx','By','Bz' 2301 WRITE(LUPRI,9008) 'Ex',(POLDL(1,1,JDIP,IFRVAL), JDIP=1,3) 2302 WRITE(LUPRI,9008) 'Ey',(POLDL(1,2,JDIP,IFRVAL), JDIP=1,3) 2303 WRITE(LUPRI,9008) 'Ez',(POLDL(1,3,JDIP,IFRVAL), JDIP=1,3) 2304 IF (ABSORP) THEN 2305 WRITE(LUPRI,'(/10X,A)') 'Imaginary part:' 2306 WRITE(LUPRI,9009) 'Bx','By','Bz' 2307 WRITE(LUPRI,9008) 'Ex',(POLDL(2,1,JDIP,IFRVAL), 2308 & JDIP=1,3) 2309 WRITE(LUPRI,9008) 'Ey',(POLDL(2,2,JDIP,IFRVAL), 2310 & JDIP=1,3) 2311 WRITE(LUPRI,9008) 'Ez',(POLDL(2,3,JDIP,IFRVAL), 2312 & JDIP=1,3) 2313 END IF 2314 END DO 2315 DO IFRVAL = 1,NFRVAL 2316 WRITE(LUPRI,9006) FRVAL(IFRVAL) 2317 WRITE(LUPRI,9007) 2318 IF (ABSORP) WRITE(LUPRI,'(/10X,A)') 'Real part:' 2319 WRITE(LUPRI,9009) 'Bx','By','Bz' 2320 WRITE(LUPRI,9008) 'Ex',(POLDA(1,1,JDIP,IFRVAL), JDIP=1,3) 2321 WRITE(LUPRI,9008) 'Ey',(POLDA(1,2,JDIP,IFRVAL), JDIP=1,3) 2322 WRITE(LUPRI,9008) 'Ez',(POLDA(1,3,JDIP,IFRVAL), JDIP=1,3) 2323 IF (ABSORP) THEN 2324 WRITE(LUPRI,'(/10X,A)') 'Imaginary part:' 2325 WRITE(LUPRI,9009) 'Bx','By','Bz' 2326 WRITE(LUPRI,9008) 'Ex',(POLDA(2,1,JDIP,IFRVAL), 2327 & JDIP=1,3) 2328 WRITE(LUPRI,9008) 'Ey',(POLDA(2,2,JDIP,IFRVAL), 2329 & JDIP=1,3) 2330 WRITE(LUPRI,9008) 'Ez',(POLDA(2,3,JDIP,IFRVAL), 2331 & JDIP=1,3) 2332 END IF 2333 END DO 2334 IF (MVEOR) THEN 2335 DO IFRVAL = 1,NFRVAL 2336 TSTFR = ABS(FRVAL(IFRVAL)) 2337 IF (TSTFR .LT. 1.0D-8) THEN 2338 WRITE(LUPRI,9010) FRVAL(IFRVAL) 2339 WRITE(LUPRI,9007) 2340 WRITE(LUPRI,'(10X,A)') '--- undefined ---' 2341 ELSE 2342 FRQINV = 1.0D0/FRVAL(IFRVAL) 2343 CALL DSCAL(18,FRQINV,POLVL(1,1,1,IFRVAL),1) 2344 WRITE(LUPRI,9010) FRVAL(IFRVAL) 2345 WRITE(LUPRI,9007) 2346 IF (ABSORP) WRITE(LUPRI,'(/10X,A)') 'Real part:' 2347 WRITE(LUPRI,9009) 'Bx','By','Bz' 2348 WRITE(LUPRI,9008) 'Ex',(POLVL(1,1,JDIP,IFRVAL), 2349 & JDIP=1,3) 2350 WRITE(LUPRI,9008) 'Ey',(POLVL(1,2,JDIP,IFRVAL), 2351 & JDIP=1,3) 2352 WRITE(LUPRI,9008) 'Ez',(POLVL(1,3,JDIP,IFRVAL), 2353 & JDIP=1,3) 2354 IF (ABSORP) THEN 2355 WRITE(LUPRI,'(/10X,A)') 'Imaginary part:' 2356 WRITE(LUPRI,9009) 'Bx','By','Bz' 2357 WRITE(LUPRI,9008) 'Ex',(POLVL(2,1,JDIP,IFRVAL), 2358 & JDIP=1,3) 2359 WRITE(LUPRI,9008) 'Ey',(POLVL(2,2,JDIP,IFRVAL), 2360 & JDIP=1,3) 2361 WRITE(LUPRI,9008) 'Ez',(POLVL(2,3,JDIP,IFRVAL), 2362 & JDIP=1,3) 2363 END IF 2364 END IF 2365 END DO 2366 DO IFRVAL = 1,NFRVAL 2367 TSTFR = ABS(FRVAL(IFRVAL)) 2368 IF (TSTFR .LT. 1.0D-8) THEN 2369 WRITE(LUPRI,9011) FRVAL(IFRVAL) 2370 WRITE(LUPRI,9007) 2371 WRITE(LUPRI,'(10X,A)') '--- undefined ---' 2372 ELSE 2373 FRQINV = 1.0D0/FRVAL(IFRVAL) 2374 CALL DSCAL(18,FRQINV,POLVV(1,1,1,IFRVAL),1) 2375 WRITE(LUPRI,9011) FRVAL(IFRVAL) 2376 WRITE(LUPRI,9007) 2377 IF (ABSORP) WRITE(LUPRI,'(/10X,A)') 'Real part:' 2378 WRITE(LUPRI,9009) 'Bx','By','Bz' 2379 WRITE(LUPRI,9008) 'Ex',(POLVV(1,1,JDIP,IFRVAL), 2380 & JDIP=1,3) 2381 WRITE(LUPRI,9008) 'Ey',(POLVV(1,2,JDIP,IFRVAL), 2382 & JDIP=1,3) 2383 WRITE(LUPRI,9008) 'Ez',(POLVV(1,3,JDIP,IFRVAL), 2384 & JDIP=1,3) 2385 IF (ABSORP) THEN 2386 WRITE(LUPRI,'(/10X,A)') 'Imaginary part:' 2387 WRITE(LUPRI,9009) 'Bx','By','Bz' 2388 WRITE(LUPRI,9008) 'Ex',(POLVV(2,1,JDIP,IFRVAL), 2389 & JDIP=1,3) 2390 WRITE(LUPRI,9008) 'Ey',(POLVV(2,2,JDIP,IFRVAL), 2391 & JDIP=1,3) 2392 WRITE(LUPRI,9008) 'Ez',(POLVV(2,3,JDIP,IFRVAL), 2393 & JDIP=1,3) 2394 END IF 2395 END IF 2396 END DO 2397 END IF 2398 END IF 2399 TMASS = 0.0D0 2400 DO IATOM = 1, NUCIND 2401 DO ISYMOP = 0, MAXOPR 2402 IF (IAND(ISYMOP,ISTBNU(IATOM)) .EQ. 0) THEN 2403 NATTYP = IZATOM(IATOM) 2404 IF (NATTYP .NE. 0 .AND. .NOT. NOORBT(IATOM)) THEN 2405 AMASS = DISOTP(NATTYP,ISOTOP(IATOM),'MASS') 2406 TMASS = TMASS + AMASS 2407 END IF 2408 END IF 2409 END DO 2410 END DO 2411 FACTOT = FACTOR*XTKAYS*XTKAYS 2412 IF (ABSORP) THEN 2413 WRITE(LUPRI,'(2(/,10X,A))') 2414 & 'Optical rotation and Electronic Circular Dicroism summary', 2415 & '---------------------------------------------------------' 2416 WRITE(LUPRI,'(/,10X,A)') 2417 & 'ORD is given in units:'// 2418 & ' degrees/(dm g/cm^3)' 2419 WRITE(LUPRI,'(10X,A)') 2420 & 'CD extinction coefficient is given in units:'// 2421 & ' L/(mol cm)' 2422 ELSE 2423 WRITE(LUPRI,'(2(/,10X,A))') 2424 & 'Optical rotation summary', 2425 & '------------------------' 2426 WRITE(LUPRI,'(/,10X,A)') 2427 & 'ORD is given in units:'// 2428 & ' degrees/(dm g/cm^3)' 2429 END IF 2430 DO IFRVAL = 1, NFRVAL 2431 TSTFR = ABS(FRVAL(IFRVAL)) 2432 IF (TSTFR .LT. 1.0D-14) THEN 2433 WRITE (LUPRI,'(/10X,A,F12.6,A)') 2434 & 'Frequency: ',FRVAL(IFRVAL),' au' 2435 BETAL(1) = 0.0D0 2436 BETANL(1) = 0.0D0 2437 BETAL(2) = 0.0D0 2438 BETANL(2) = 0.0D0 2439 BETAVL(1) = 0.0D0 2440 BETAVV(1) = 0.0D0 2441 BETAVL(2) = 0.0D0 2442 BETAVV(2) = 0.0D0 2443 ELSE 2444 WRITE (LUPRI,'(/10X,A,F12.6,A,F12.6,A)') 2445 & 'Frequency: ',FRVAL(IFRVAL),' au = ', 2446 & XTNM/FRVAL(IFRVAL), ' nm' 2447 XXXDIV = 1.0D0/(D3*FRVAL(IFRVAL)) 2448 BETAL(1) =-(POLDL(1,1,1,IFRVAL) + POLDL(1,2,2,IFRVAL) 2449 & + POLDL(1,3,3,IFRVAL))*XXXDIV 2450 BETANL(1)=-(POLDA(1,1,1,IFRVAL) + POLDA(1,2,2,IFRVAL) 2451 & + POLDA(1,3,3,IFRVAL))*XXXDIV 2452 BETAL(2) =-(POLDL(2,1,1,IFRVAL) + POLDL(2,2,2,IFRVAL) 2453 & + POLDL(2,3,3,IFRVAL))*XXXDIV 2454 BETANL(2)=-(POLDA(2,1,1,IFRVAL) + POLDA(2,2,2,IFRVAL) 2455 & + POLDA(2,3,3,IFRVAL))*XXXDIV 2456 BETAVL(1)=-(POLVL(1,1,1,IFRVAL) + POLVL(1,2,2,IFRVAL) 2457 & + POLVL(1,3,3,IFRVAL))*XXXDIV 2458 BETAVV(1)=-(POLVV(1,1,1,IFRVAL) + POLVV(1,2,2,IFRVAL) 2459 & + POLVV(1,3,3,IFRVAL))*XXXDIV 2460 BETAVL(2)=-(POLVL(2,1,1,IFRVAL) + POLVL(2,2,2,IFRVAL) 2461 & + POLVL(2,3,3,IFRVAL))*XXXDIV 2462 BETAVV(2)=-(POLVV(2,1,1,IFRVAL) + POLVV(2,2,2,IFRVAL) 2463 & + POLVV(2,3,3,IFRVAL))*XXXDIV 2464 END IF 2465 XXXCON = FACTOT*FRVAL(IFRVAL)*FRVAL(IFRVAL)/TMASS 2466 ALPHAL(1)=XXXCON*BETAL(1) 2467 ALPHNL(1)=XXXCON*BETANL(1) 2468 ALPHAL(2)=XXXCON*BETAL(2) 2469 ALPHNL(2)=XXXCON*BETANL(2) 2470 ALPHVL(1)=XXXCON*BETAVL(1) 2471 ALPHVV(1)=XXXCON*BETAVV(1) 2472 ALPHVL(2)=XXXCON*BETAVL(2) 2473 ALPHVV(2)=XXXCON*BETAVV(2) 2474 IF (ABSORP) THEN 2475 WRITE(LUPRI,'(10X,A,F12.6,A,F12.2,A)') 2476 & 'Damping : ',DAMPING,' au = ', 2477 & DAMPING*XTKAYS,' cm-1' 2478 IF (MVEOR) THEN 2479 WRITE(LUPRI,'(8(/10X,2(A,F15.6)))') 2480 & 'Beta (London) :', 2481 & BETAL(1), ' + i ',BETAL(2), 2482 & 'Beta (No-London) :', 2483 & BETANL(1),' + i ',BETANL(2), 2484 & 'Beta (Velocity) :', 2485 & BETAVL(1), ' + i ',BETAVL(2), 2486 & 'Beta (Modified Velocity) :', 2487 & BETAVV(1),' + i ',BETAVV(2), 2488 & 'Optical rotation (London) :', 2489 & ALPHAL(1),' + i ',ALPHAL(2), 2490 & 'Optical rotation (No-London) :', 2491 & ALPHNL(1),' + i ',ALPHNL(2), 2492 & 'Optical rotation (Velocity) :', 2493 & ALPHVL(1),' + i ',ALPHVL(2), 2494 & 'Optical rotation (Mod. Vel.) :', 2495 & ALPHVV(1),' + i ',ALPHVV(2) 2496 ELSE 2497 WRITE(LUPRI,'(2(/10X,2(A,F12.6)))') 2498 & 'Beta (London) :', 2499 & BETAL(1), ' + i ',BETAL(2), 2500 & 'Beta (No-London) :', 2501 & BETANL(1),' + i ',BETANL(2) 2502 WRITE(LUPRI,'(4(/10X,A,F16.6))') 2503 & 'Optical rotation (London) :', 2504 & ALPHAL(1), 2505 & 'Optical rotation (No-London) :', 2506 & ALPHNL(1), 2507 & 'Circular dichroism (London) :', 2508 & ALPHAL(2)*TMASS/3.2988D5, 2509 & 'Circular dichroism (No-London) :', 2510 & ALPHNL(2)*TMASS/3.2988D5 2511 END IF 2512 ELSE 2513 IF (MVEOR) THEN 2514 WRITE(LUPRI,'(8(/10X,A,F16.6))') 2515 & 'Beta (London) :',BETAL(1), 2516 & 'Beta (No-London) :',BETANL(1), 2517 & 'Beta (Velocity) :',BETAVL(1), 2518 & 'Beta (Modified Velocity) :',BETAVV(1), 2519 & 'Optical rotation (London) :',ALPHAL(1), 2520 & 'Optical rotation (No-London) :',ALPHNL(1), 2521 & 'Optical rotation (Velocity) :',ALPHVL(1), 2522 & 'Optical rotation (Mod. Vel.) :',ALPHVV(1) 2523 ELSE 2524 WRITE(LUPRI,'(4(/10X,A,F16.6))') 2525 & 'Beta (London) :',BETAL(1), 2526 & 'Beta (No-London) :',BETANL(1), 2527 & 'Optical rotation (London) :',ALPHAL(1), 2528 & 'Optical rotation (No-London) :',ALPHNL(1) 2529 END IF 2530 END IF 2531 IF (IPRINT .GE. 1) THEN 2532 WRITE(LUPRI,'(/A,1P,D16.9)') 2533 & 'Conversion, beta --> alpha: ',XXXCON 2534 END IF 2535 END DO 2536 ENDIF 2537C 2538 IF (ABSORP .AND. IPRINT.GE.0 .AND. (.NOT. ABSLRS)) THEN 2539 WRITE(LUPRI,'(/10X,A4,A6,A21,/10X,A)') 2540 & 'Sym.','State','Excitation energy', 2541 & '---------------------------------------------------' 2542 DO ISYM=1,NSYM 2543 IF (NOPER(ISYM).GT.0) THEN 2544 DO ISTATE=1,NEXCITED_STATES 2545 WRITE(LUPRI,'(6X,2I6,F12.4,A,F9.3,A,F10.2,A)') 2546 & ISYM,ISTATE,EXC_ENERGY(ISTATE,ISYM),' a.u.', 2547 & EXC_ENERGY(ISTATE,ISYM)*XTEV,' eV', 2548 & XTNM/EXC_ENERGY(ISTATE,ISYM),' nm' 2549 END DO 2550 END IF 2551 END DO 2552 END IF 2553C 2554 IF (VIB_G) THEN 2555 NCOOR = 3*NUCDEP 2556 KCSTRA = 1 2557 KSCTRA = KCSTRA + NCOOR*NCOOR 2558 KLAST = KSCTRA + NCOOR*NCOOR 2559 IF (KLAST .GT. LWORK) CALL STOPIT('LNROUT',' ',KLAST,LWORK) 2560 CALL TRACOR(WORK(KCSTRA),WORK(KSCTRA),1,NCOOR,0) 2561C 2562 CALL HEADER('Vibrational g-factor',-1) 2563CSPAS:8/10-08 there is no frequency for VIB_G 2564C DO IFRVAL = 1,NFRVAL 2565 IFRVAL = 1 2566CSPAS WRITE(LUPRI,9110) FRVAL(IFRVAL) 2567CKeinSPASmehr 2568 WRITE(LUPRI,9111) 2569 WRITE(LUPRI,'(16X,A,21X,A,5X,A,8X,A)') 2570 & 'Coordinates','Sym','(au)','(cm-1 u)' 2571 WRITE(LUPRI,'(4X,2A)') 2572 & ('--------------------------------------',I=1,2) 2573 DO ISYM = 1,NSYM 2574CSPAS:20/3-2011: allowing for VIB_G perturbations in only one irrep 2575 IF (.NOT.VIBGIR .OR. (ISYM .EQ. IRVIBG)) THEN 2576CKeinSPASmehr 2577 DO ISCOOR = 1,3*NATOMS 2578 CALL GENCOR(WORK(KCSTRA),NCOOR,ISCOOR,COORDI,ICORSY) 2579 IF (ICORSY .EQ. ISYM) THEN 2580 DO JSCOOR = 1,3*NATOMS 2581 CALL GENCOR(WORK(KCSTRA),NCOOR,JSCOOR,COORDJ, 2582 & JCORSY) 2583 IF (JCORSY .EQ. ISYM) THEN 2584 VAL = FOVIBG(ISCOOR,JSCOOR,ISYM,IFRVAL) 2585 WRITE(LUPRI,'(5X,2A,I3,F12.6,F14.3)') 2586 & COORDI, COORDJ, ISYM, 2587 & VAL,VAL*XTKAYS/XFAMU 2588 END IF 2589 END DO 2590CSPAS:6/10-08: trying to add mass independent vibrational g-factor 2591 DO JDIPV = 1,3 2592 IF (JDIPV.EQ.1) THEN 2593 COORDJ = "XDIPVEL" 2594 JDIPSY = ISYMAX(1,1) + 1 2595 ELSE IF (JDIPV.EQ.2) THEN 2596 COORDJ = "YDIPVEL" 2597 JDIPSY = ISYMAX(2,1) + 1 2598 ELSE IF (JDIPV.EQ.3) THEN 2599 COORDJ = "ZDIPVEL" 2600 JDIPSY = ISYMAX(3,1) + 1 2601 END IF 2602 IF (JDIPSY .EQ. ISYM) THEN 2603 WRITE(LUPRI,'(5X,2A,I3,F12.6)') 2604 & COORDI, COORDJ, ISYM, 2605 & FOVIBG(ISCOOR,3*NATOMS+JDIPV,ISYM,IFRVAL) 2606 END IF 2607 END DO 2608CKeinSPASmehr 2609 END IF 2610 END DO 2611CSPAS:6/10-08: trying to add mass independent vibrational g-factor 2612 DO IDIPV = 1,3 2613 IF (IDIPV.EQ.1) THEN 2614 COORDI = "XDIPVEL" 2615 IDIPSY = ISYMAX(1,1) + 1 2616 ELSE IF (IDIPV.EQ.2) THEN 2617 COORDI = "YDIPVEL" 2618 IDIPSY = ISYMAX(2,1) + 1 2619 ELSE IF (IDIPV.EQ.3) THEN 2620 COORDI = "ZDIPVEL" 2621 IDIPSY = ISYMAX(3,1) + 1 2622 END IF 2623 IF (IDIPSY .EQ. ISYM) THEN 2624 DO JSCOOR = 1,3*NATOMS 2625 CALL GENCOR(WORK(KCSTRA),NCOOR,JSCOOR,COORDJ, 2626 & JCORSY) 2627 IF (JCORSY .EQ. ISYM) THEN 2628 WRITE(LUPRI,'(5X,2A,I3,F12.6)') 2629 & COORDI, COORDJ, ISYM, 2630 & FOVIBG(3*NATOMS+IDIPV,JSCOOR,ISYM,IFRVAL) 2631 END IF 2632 END DO 2633 DO JDIPV = 1,3 2634 IF (JDIPV.EQ.1) THEN 2635 COORDJ = "XDIPVEL" 2636 JDIPSY = ISYMAX(1,1) + 1 2637 ELSE IF (JDIPV.EQ.2) THEN 2638 COORDJ = "YDIPVEL" 2639 JDIPSY = ISYMAX(2,1) + 1 2640 ELSE IF (JDIPV.EQ.3) THEN 2641 COORDJ = "ZDIPVEL" 2642 JDIPSY = ISYMAX(3,1) + 1 2643 END IF 2644 IF (JDIPSY .EQ. ISYM) THEN 2645 WRITE(LUPRI,'(5X,2A,I3,F12.6)') 2646 & COORDI, COORDJ, ISYM, 2647 & FOVIBG(3*NATOMS+IDIPV,3*NATOMS+JDIPV,ISYM,IFRVAL) 2648 END IF 2649 END DO 2650 END IF 2651 END DO 2652CKeinSPASmehr 2653CSPAS:20/3-2011: allowing for VIB_G perturbations in only one irrep 2654 END IF 2655CKeinSPASmehr 2656 END DO 2657 WRITE(LUPRI,9111) 2658CSPAS:8/10-08 there is no frequency for VIB_G 2659C END DO 2660CKeinSPASmehr 2661 END IF 2662C 2663 9001 FORMAT (/10X,'Polarizability tensor for frequency ',F12.6,' au') 2664 9002 FORMAT (10X,'----------------------------------------------------' 2665 & ) 2666 9003 FORMAT (16X,3(A14,1X),/) 2667 9004 FORMAT (16X,A,2X,1P,3(G14.7,1X)) 2668 9005 FORMAT (/10X,'London G tensor for frequency ',F12.6,' au') 2669 9006 FORMAT (/10X,'No-London G tensor for frequency ',F12.6,' au') 2670 9007 FORMAT (10X,'-------------------------------------------------') 2671 9008 FORMAT (16X,A,2X,3F12.7) 2672 9009 FORMAT (27X,3(A,10X)/) 2673 9010 FORMAT (/10X,'Velocity G tensor for frequency ',F12.6,' au') 2674 9011 FORMAT (/10X,'Mod. Vel. G tensor for frequency ',F12.6,' au') 2675 9110 FORMAT (14X,'Vibrational g-factors for frequency ',F12.6,' au',/) 2676 9111 FORMAT (4X,76('=')) 2677C 2678 RETURN 2679 END 2680