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 herctl */ 20#ifndef PRG_DIRAC 21 SUBROUTINE HERCTL(WORK,LWORK) 22C 23C The code for calculation of one- and two-electron integrals 24C was written by T. Helgaker in 1984 at the University of Oslo. 25C 26C General contractions were implemented by T. Helgaker at the 27C University of Aarhus in Feb-Mar 1988. 28C 29C Symmetry was implemented by P. R. Taylor and T. Helgaker at 30C NASA Ames in Apr 1988. 31C 32C The supermatrix code is written by O. Kvalheim at the University 33C of Bergen. 34C 35C Spin-orbit integrals by O. Vahtras and T. Helgaker at the 36C University of Oslo, Nov 1989. 37C 38C Cartesian and spherical moments integrals by T. Helgaker, 39C University of Oslo, Sep and Oct 1990 40C 41C Integrals for indirect spin-spin coupling tensors by T. Helgaker 42C and O. Vahtras at the University of Oslo, Feb 1991. 43C 44C Half-derivative overlap integrals for NACMES by T. Helgaker, 45C at University of Aarhus, Jun 1991 46C 47C Overlap matrix differentiated with respect to external magetic 48C field, K. Ruud & T. Helgaker, Oct 1991 49C 50C Electronic angular momentum around fixed center or nuclei, 51C K. Ruud & T. Helgaker, Oct 1991 52C 53C One-electron contribution to the magnetic moment of molecules, 54C K. Ruud & T. Helgaker, Nov. 1991 55C 56C Kinetic energy integrals, K. Ruud, Nov. 1991 57C 58C Cosine and sine integrals, T. Helgaker, Jun. 1993 59C 60C Mass-velocity and Darwin integrals 61C S. Kirpekar & H.J.Aa. Jensen, Jul. 1993 62C 63C Magnetic field derivative integrals of dipole length 64C K.Ruud, Aug.-93 65C 66 67#include "implicit.h" 68#include "priunit.h" 69#include "dummy.h" 70#include "mxcent.h" 71#include "qm3.h" 72#include "maxorb.h" 73C 74#include "cbiher.h" 75#include "cbihr1.h" 76#include "cbihr2.h" 77#include "cbihrs.h" 78#include "cbieri.h" 79#include "gnrinf.h" 80#include "exeinf.h" 81#include "inftap.h" 82#include "huckel.h" 83#include "r12int.h" 84#include "pcmlog.h" 85C 86Cholesky 87#include "ccdeco.h" 88Cholesky 89C 90C inftra : USE_INTSORT 91C 92#include "inftra.h" 93C 94 95 LOGICAL SET, FEXIST 96C 97 DIMENSION WORK(LWORK) 98C 99 CALL QENTER('HERCTL') 100 CALL GETTIM(TIMHER,WALHER) 101C 102C ************************* 103C ***** Input Section ***** 104C ************************* 105C 106 TIMSTR = SECOND() 107 CALL HERINP(WORK,LWORK) 108 IF (TSTINP) THEN 109 WRITE (LUPRI,'(/15X,A)') 110 & '*** End of input test for HERMIT ***' 111 CALL QUIT('*** End of input test for HERMIT ***') 112 END IF 113 TIMINP = SECOND() - TIMSTR 114 CALL FLSHFO(LUPRI) 115C 116C ************************************ 117C ***** Calculation of Integrals ***** 118C ************************************ 119C 120 121 CALL GPINQ('AOPROPER','EXIST',FEXIST) 122 IF (FEXIST) THEN 123 IF (LUPROP .LE. 0) CALL GPOPEN(LUPROP,'AOPROPER','OLD',' ',' ', 124 & IDUMMY,.FALSE.) 125 CALL GPCLOSE(LUPROP,'DELETE') 126 END IF 127C 128 IF (QM3) THEN 129 QM3LO1 = .TRUE. 130 QM3LO2 = .TRUE. 131 LONUPO = .FALSE. 132 LOELFD = .FALSE. 133 CALL GPOPEN(LUQM3E,'ELFDMM','UNKNOWN',' ',' ', 134 & IDUMMY,.FALSE.) 135 CALL GPOPEN(LUQM3P,'POTMM','UNKNOWN',' ',' ', 136 & IDUMMY,.FALSE.) 137 FORQM3 = .TRUE. 138C make sure QM3 integrals are evaluated for relevant properties; 139C is reset to .false. before return (is used in other routines 140C later to signal when integrals are for QM3 and when not). 141 ELSE 142 FORQM3 = .FALSE. 143 END IF 144C 145 RUNQM3 = .FALSE. 146C 147 IF (DORLM) CALL GPOPEN(LUSOL, 148 & 'AOSOLINT','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 149 150C Delete old MO integrals integrals 151 IF (USE_INTSORT) THEN 152 IF (LUORDA .LE. 0) CALL DAOPEN(LUORDA,'AOORDINT') 153 CALL DARMOV(LUORDA) 154 IF (LUMINT .LE. 0) CALL DAOPEN(LUMINT,'MOTWOINT') 155 CALL DARMOV(LUMINT) 156 ELSE 157 CALL GPINQ('AOSRTINT.DA','EXIST',FEXIST) 158 IF (FEXIST) THEN 159 CALL SYSTEM('rm -f AOSRTINT.DA*') 160 END IF 161 CALL GPINQ('MOTWOINT','EXIST',FEXIST) 162 IF (FEXIST) THEN 163 CALL GPOPEN(LUINTM,'MOTWOINT', 164 & 'UNKNOWN',' ',' ',IDUMMY,.FALSE.) 165 CALL GPCLOSE(LUINTM,'DELETE') 166 END IF 167 END IF 168 169C Delete any existing old 2-el. AO integrals 170C (they may have been generated in previous run 171C even if RUNTWO is .false. now ! This will for example be 172C the case if DIRCAL and DOMP2)) 173 174C Let us give a chance to the user who knows what she/he does: 175C if RMAOTWO false then it is user's responsibility that 176C all the needed AOTWOINT/AOSR2IT/AUSUPINT files are available 177C 178 IF (RMAOTWO) THEN 179 IF (LUINTA .LE. 0) CALL GPOPEN(LUINTA,'AOTWOINT', 180 & 'UNKNOWN',' ',' ',IDUMMY,.FALSE.) 181 CALL GPCLOSE(LUINTA,'DELETE') 182C 183 IF (LUSRINT .LE. 0) CALL GPOPEN(LUSRINT,'AOSR2INT', 184 & 'UNKNOWN',' ',' ',IDUMMY,.FALSE.) 185 CALL GPCLOSE(LUSRINT,'DELETE') 186 IF (LUSUPM .LE. 0) CALL GPOPEN(LUSUPM,'AOSUPINT', 187 & 'UNKNOWN',' ',' ',IDUMMY,.FALSE.) 188 CALL GPCLOSE(LUSUPM,'DELETE') 189 END IF ! RMAOTWO 190 191 ! Postpone calculation of two-electron integrals; 192 ! 2-el integrals will be calculated with MAKE_AOTWOINT 193 ! calls when needed. 194 ! Except for all R12-type integrals. 195 196 IF (RUNTWO) THEN 197 IF (.NOT.(R12CAL .OR. DCCR12 .OR. R12EIN 198 & .OR. R12INT .OR. U12INT .OR. U21INT)) THEN 199 WRITE(LUPRI,'(/A)') 'Calculation of 2-electron integrals'// 200 & ' are postponed until they are needed.' 201 RUNTWO = .FALSE. 202 END IF 203 END IF 204 205 206C IF (RUNTWO .or. RUNERI) : 207C Open LUINTA and LUSRINT for new 2-el. integrals 208C to be calculated in the CALL HERINT call below 209 210 IF ((RUNTWO .OR. RUNERI) .AND. .NOT. DCCR12) THEN 211 CALL GPOPEN(LUINTA,'AOTWOINT','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 212 IF (DOSRIN) CALL 213 & GPOPEN(LUSRINT,'AOSR2INT','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 214 END IF 215C 216 CALL GPOPEN(LUONEL,'AOONEINT','OLD',' ',' ',IDUMMY,.FALSE.) 217 CALL MOLLAB('INTEGRAL',LUONEL,LUPRI) 218 CALL HERINT(WORK,LWORK) 219 220 IF (PCM) THEN 221 CALL IEFMAT(WORK,LWORK) 222 END IF 223C 224 IF ( (QM3) .AND. (.NOT.INTDIR) ) THEN 225 DO 752 KI = 1, ISYTP 226 IF (MDLWRD(KI)(1:3) .EQ. 'SPC') THEN 227 IF ( .NOT. (LONUPO) ) THEN 228 CALL QUIT('Error, NUCPOT integrals not computed') 229 END IF 230 END IF 231 752 CONTINUE 232C 233 DO 753 KI = 1, ISYTP 234 IF (MDLWRD(KI)(1:5) .EQ. 'SPC_E') THEN 235 IF ( .NOT. (LOELFD) ) THEN 236 CALL QUIT('Error, NELFLD integrals not computed') 237 END IF 238 END IF 239 753 CONTINUE 240 END IF 241 242#if defined(BUILD_GEN1INT) 243C test suite of Gen1Int interface, added by Bin Gao, Oct. 7, 2011 244 if (TEST_GEN1INT) 245 & call gen1int_host_test(LWORK, WORK, LUPRI, IPRDEF) 246#endif 247 248C 249C ********************************** 250C ***** P Supermatrix Elements ***** 251C ********************************** 252C 253C Integrals on AOSUPINT 254C HJAaJ Oct 2003: CALL FORMSUP has been moved to SIRIUS and 255C is only done when relevant ... 256C 257C ********************************** 258C ***** Integral sorting ***** 259C ********************************** 260C 261C Integrals on AOORDINT 262C 263 IF (USE_INTSORT .AND. CHOINT) THEN 264 WRITE(LUPRI,'(/A)') 265 & ' * INFO: .SORT INT ignored when .CHOLESKY' 266 USE_INTSORT = .FALSE. 267 END IF 268 IF (USE_INTSORT) THEN 269 CALL GPINQ('AOTWOINT','EXIST',FEXIST) 270 IF (FEXIST) THEN 271 CALL TIMER('START ',TIMSTR,TIMEND) 272 CALL SORTAO(WORK,LWORK) 273 CALL TIMER('SORTAO',TIMSTR,TIMEND) 274 CALL FLSHFO(LUPRI) 275 ELSE 276 NWARN = NWARN + 1 277 WRITE(LUPRI,'(/A)') 278 & ' * WARNING: .SORT INT requested but no integrals to sort!' 279 END IF 280 END IF 281C 282C ****************************** 283C ***** End of Calculation ***** 284C ****************************** 285C 286 CALL GPCLOSE(LUONEL,'KEEP') 287 IF (LUPROP .GT. 0) CALL GPCLOSE(LUPROP,'KEEP') 288 IF (LUQM3E .GT. 0) CALL GPCLOSE(LUQM3E,'KEEP') 289 IF (LUQM3P .GT. 0) CALL GPCLOSE(LUQM3P,'KEEP') 290 IF (LUSOL .GT. 0) CALL GPCLOSE(LUSOL ,'KEEP') 291C 292 IF (LUINTA .GT. 0) CALL GPCLOSE(LUINTA ,'KEEP') 293 IF (LUSRINT .GT. 0) CALL GPCLOSE(LUSRINT,'KEEP') 294C 295 IF (RUNTWO .OR. RUNERI) FTRCTL = .TRUE. 296C ..transformation control, old AO/MO files cannot be used any more. 297C 298 CALL GETTIM(TEND,WEND) 299 TIMHER = TEND - TIMHER 300 WALHER = WEND - WALHER 301 CALL TIMTXT(' Total CPU time used in HERMIT:',TIMHER,LUPRI) 302 CALL TIMTXT(' Total wall time used in HERMIT:',WALHER,LUPRI) 303 CALL QEXIT('HERCTL') 304 RETURN 305 END 306#endif 307C /* Deck herinp */ 308 SUBROUTINE HERINP(WORK,LWORK) 309C 310C --- General Input for HERMIT --- 311C 312C Trygve Helgaker, extensively extended by others 313C 314C 950602-vebjorn: 315C Added flag HRINPC to ensure that HERMIT input processing is done only once. 316C 317#include "implicit.h" 318#include "priunit.h" 319#include "iratdef.h" 320#include "mxcent.h" 321#include "maxorb.h" 322#include "maxaqn.h" 323C 324C NTABLE has been increased 118 (WK/UniKA/07-25-2005). 325 PARAMETER (NDIR = 9, NTABLE = 146) 326 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, D4 = 4.0D0, D10 = 1.0D1) 327 CHARACTER WORD*7, PROMPT*1, TABDIR(NDIR)*7, TABLE(NTABLE)*7, 328 & WORD1*7 329 DIMENSION WORK(LWORK) 330C 331#ifdef PRG_DIRAC 332#include "dcbgen.h" 333#include "hrunit.h" 334#include "cbieri.h" 335#else 336C 337C inftra : USE_INTSORT 338C 339#include "dummy.h" 340#include "gnrinf.h" 341#include "cbihrs.h" 342#include "inftap.h" 343#include "infinp.h" 344#include "inftra.h" 345#endif 346C 347#include "cbiher.h" 348#include "cbihr1.h" 349#include "nuclei.h" 350#include "ccom.h" 351#include "orgcom.h" 352 353#include "efield.h" 354#include "cbisol.h" 355#include "qm3.h" 356#include "qmmm.h" 357#include "huckel.h" 358#include "drw2el.h" 359#include "elweak.h" 360#include "r12int.h" 361#include "codata.h" 362#include "denfit.h" 363#include "infpar.h" 364 365C Keywords for R12 have been added (WK/UniKA/04-11-2002). 366cLig <> added RANGMO and RPSO 367C Bin Gao, December 17, 2009; adds integrals S2MBRA, S2MKET, and S2MMIX 368 DATA TABDIR /'*END OF', '*READIN', '*ONEINT', '*TWOINT', 369 & '*SUPINT', '*ER2INT', '*SORINT', '*DENFIT', 370 & '*QM3INP'/ 371C adds test suite of Gen1Int interface as term 239 372C Bin Gao, Oct. 3, 2011 373! 101 102 103 104 374 DATA TABLE /'.PRINT ', '.INPTES', '.NOSUP ', '.SPIN-O', 375! 105 106 107 108 376 & '.DIPLEN', '.NO HAM', '.SOTEST', '.DIPVEL', 377! 109 110 111 112 378 & '.QUADRU', '.PHASEO', '.SECMOM', '.SUPONL', 379! 113 114 115 116 380 & '.CARMOM', '.SPHMOM', '.FC ', '.PSO ', 381! 117 118 119 120 382 & '.SD ', '.DSO ', '.POINTS', '.SELECT', 383! 121 122 123 124 384 & '.QUASUM', '.SD+FC ', '.PROPRI', '.HDO ', 385! 125 126 127 128 386 & '.S1MAG ', '.S2MAG ', '.ANGMOM', '.ANGLON', 387! 129 130 131 132 388 & '.LONMOM', '.MAGMOM', '.S1MAGT', '.MGMOMT', 389! 133 134 135 136 390 & '.KINENE', '.S2MAGT', '.DSUSNL', '.DSUSLL', 391! 137 138 139 140 392 & '.DSUSLH', '.DIASUS', '.DSUTST', '.NSTNOL', 393! 141 142 143 144 394 & '.NSTLON', '.NST ', '.NSNLTS', '.NSLTST', 395! 145 146 147 148 396 & '.NELFLD', '.NSTTST', '.EFGCAR', '.EFGSPH', 397! 149 150 151 152 398 & '.S1MAGL', '.S1MAGR', '.HDOBR ', '.S1MLT ', 399! 153 154 155 156 400 & '.HDOBRT', '.S1MRT ', '.NUCPOT', '.NPOTST', 401! 157 158 159 160 402 & '.MGMO2T', '.MGMTHR', '.HBDO ', '.SUSCGO', 403! 161 162 163 164 404 & '.NSTCGO', '.EXPIKR', '.MASSVE', '.DARWIN', 405! 165 166 167 168 406 & '.CM-1 ', '.CM-2 ', '.SQHDOL', '.SQHDOR', 407! 169 170 171 172 408 & '.NOTWO ', '.GAUGEO', '.DIPORG', '.NO2SO ', 409! 173 174 175 176 410 & '.S1ELE ', '.S1ELB ', '.ONEELD', '.THETA ', 411! 177 178 179 180 412 & '.NUCMOD', '.SORT I', '.DIPGRA', '.QUAGRA', 413! 181 182 183 184 414 & '.OCTGRA', '.ROTSTR', '.THIRDM', '.SOFIEL', 415! 185 186 187 188 416 & '.SOMAGM', '.DEROVL', '.DERHAM', '.ELGDIA', 417! 189 190 191 192 418 & '.ELGDIL', '.MNF-SO', '.DPTOVL', '.DPTPOT', 419! 193 194 195 196 420 & '.XDDXR3', '.AD2DAR', '.PVIOLA', '.WEINBG', 421! 197 198 199 200 422 & '.FINDPT', '.FNPROP', '.PVP ', '.1ELPOT', 423! 201 202 203 204 424 & '.QDBINT', '.QDBTST', '.R12 ', '.R12GTG', 425! 205 206 207 208 426 & '.AUXBAS', '.NOTV12', '.R12INT', '.U12INT', 427! 209 210 211 212 428 & '.U21INT', '.DCCR12', '.RANGMO', '.RPSO ', 429! 213 214 215 216 430 & '.DPTPXP', '.NOPICH', '.OZ-KE ', '.PSO-KE', 431! 217 218 219 220 432 & '.DNS-KE', '.SD-KE ', '.FC-KE ', '.DSO-KE', 433! 221 222 223 224 434 & '.PSO-OZ', '.SQHD2O', '.R12STG', '.R12RSG', 435! 225 226 227 228 436 & '.R12RGG', '.R12ERF', '.R12RRF', '.DIPLOC', 437! 229 230 231 232 438 & '.2-EL D', '.EFBDER', '.EFBDR2', '.MAGQDP', 439! 233 234 235 236 440 & '.MQDPTS', '.ORB-OR', '.DERAM ', '.DIPANH', 441! 237 238 239 240 442 & '.NEWTRA', '.KINADI', '.GENINT', '.S2MBRA', 443! 241 242 243 244 444 & '.S2MKET', '.S2MMIX', '.LRINTS', '.NORMAO', 445! 245 246 446 & '.TWOBYT', '.LFDIPL'/ 447C 448C 449 CALL QENTER('HERINP') 450C 451 RMAOTWO = .TRUE. 452C 453 IPRDEF = IPRUSR + 1 454 CALL GPOPEN(LUCMD,'DALTON.INP','OLD',' ','FORMATTED',IDUMMY, 455 & .FALSE.) 456C 457C Check if input has been processed earlier. 458C HRINPC variable makes sure HERMIT input is only read once 459C 460 IF (HRINPC) GOTO 1000 461 HRINPC = .TRUE. 462C 463C *** Initialize *** 464C 465 ! CALL HERINI - Apr. 2011: HERINI now called in EXEHER /hjaaj 466 ! so Hermit variables are also initialized if only 467 ! ABACUS 468 469 !Initialize defaults for density fitting 470 CALL DENSFIT_INI 471C 472C If WESTA - calculate integrals for WESTA program 473C 474 IF (WESTA) THEN 475 SQHDOR = .TRUE. 476 SQHD2O = .TRUE. 477 DIPLEN = .TRUE. 478 DIPVEL = .TRUE. 479 ANGMOM = .TRUE. 480 ONEPRP = .TRUE. 481 END IF 482C 483 CALL TITLER('Output from **INTEGRALS input processing (HERMIT)', 484 & '*',118) 485C 486C **** Find Hermit input ***** 487C 488 REWIND (LUCMD,IOSTAT=IOS) 489C ... IOSTAT to avoid program abort on some systems 490C if reading input from a terminal 491 900 READ (LUCMD,'(A7)',ERR=910,END=920) WORD1 492 CALL UPCASE(WORD1) 493 IF ((WORD1 .EQ. '**HERMI') .OR. (WORD1 .EQ. '*HERMIT') .OR. 494 & (WORD1 .EQ. '**INTEG')) THEN 495 GO TO 930 496 ELSE 497 GO TO 900 498 END IF 499 910 CONTINUE 500 NWARN = NWARN + 1 501 WRITE (LUPRI,'(//A)') 502 & 'WARNING **INTEGRALS input: error reading Dalton input file' 503 920 CONTINUE 504 WORD = '**END O' 505 WORD1 = '**END O' 506 WRITE (LUPRI,'(/A)') 507 & ' - Using defaults, no **INTEGRALS input found' 508 GOTO 300 509 930 CONTINUE 510C 511 CALL TITLER('Output from HERMIT input processing','*',118) 512C 513C ***** Process input for COMMON /CBIHER/ ***** 514C 515 100 READ (LUCMD, '(A7)') WORD 516 CALL UPCASE(WORD) 517 PROMPT = WORD(1:1) 518 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 519 GO TO 100 520 ELSE IF (PROMPT .EQ. '.') THEN 521 DO I = 1, NTABLE 522 IF (TABLE(I) .EQ. WORD) THEN 523 GO TO (101,102,103,104,105,106,107,108,109,110, 524 & 111,112,113,114,115,116,117,118,119,120, 525 & 121,122,123,124,125,126,127,128,129,130, 526 & 131,132,133,134,135,136,137,138,139,140, 527 & 141,142,143,144,145,146,147,148,149,150, 528 & 151,152,153,154,155,156,157,158,159,160, 529 & 161,162,163,164,165,166,167,168,169,170, 530 & 171,172,173,174,175,176,177,178,179,180, 531 & 181,182,183,184,185,186,187,188,189,190, 532 & 191,192,193,194,195,196,197,198,199,200, 533 & 201,202,203,204,205,206,207,208,209,210, 534 & 211,212,213,214,215,216,217,218,219,220, 535 & 221,222,223,224,225,226,227,228,229,230, 536 & 231,232,233,234,235,236,237,238,239,240, 537 & 241,242,243,244,245,246), I 538 END IF 539 END DO 540 IF (WORD .EQ. '.OPTION') THEN 541 CALL PRTAB(NDIR,TABDIR, WORD1//' input keywords',LUPRI) 542 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 543 GO TO 100 544 END IF 545 WRITE (LUPRI,'(/4A/)') 546 & 'ERROR: Keyword ',WORD,' not recognized under ',WORD1 547 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 548 CALL QUIT('Illegal keyword under '//WORD1) 549 101 CONTINUE 550 READ (LUCMD,*) IPRDEF 551 GO TO 100 552 102 CONTINUE 553 TSTINP = .TRUE. 554 GO TO 100 555 103 CONTINUE 556 SUPMAT = .FALSE. 557 GO TO 100 558 104 CONTINUE 559 SPNORB = .TRUE. 560 ONEPRP = .TRUE. 561 GO TO 100 562 105 CONTINUE 563 DIPLEN = .TRUE. 564 ONEPRP = .TRUE. 565 GO TO 100 566 106 CONTINUE 567 HAMILT = .FALSE. 568 SUPMAT = .FALSE. 569 GO TO 100 570 107 CONTINUE 571 SOTEST = .TRUE. 572 GO TO 100 573 108 CONTINUE 574 DIPVEL = .TRUE. 575 ONEPRP = .TRUE. 576 GO TO 100 577 109 CONTINUE 578 QUADRU = .TRUE. 579 ONEPRP = .TRUE. 580 GO TO 100 581 110 CONTINUE 582 READ (LUCMD, *,IOSTAT=IOS) (ORIGIN(I),I=1,3) 583 IF (IOS.NE.0) THEN 584 WRITE(LUPRI,*) 'Error in reading (ORIGIN(I),I=1,3)!' 585 WRITE(LUPRI,*) '(ORIGIN(I),I=1,3):', (ORIGIN(I),I=1,3) 586 CALL QUIT('Error in reading (ORIGIN(I),I=1,3)!') 587 ENDIF 588 GO TO 100 589 111 CONTINUE 590 SECMOM = .TRUE. 591 ONEPRP = .TRUE. 592 GO TO 100 593 112 CONTINUE 594 SUPMAT = .TRUE. 595 HAMILT = .FALSE. 596 NOTWO = .TRUE. 597 ONEPRP = .FALSE. 598 GO TO 100 599 113 CONTINUE 600 CARMOM = .TRUE. 601 READ (LUCMD,*) IORCAR_input 602 IORCAR = MAX(IORCAR,IORCAR_input) 603 ONEPRP = .TRUE. 604 GO TO 100 605 114 CONTINUE 606 SPHMOM = .TRUE. 607 READ (LUCMD,*) IORSPH 608 ONEPRP = .TRUE. 609 GO TO 100 610 115 CONTINUE 611 FERMI = .TRUE. 612 ONEPRP = .TRUE. 613 GO TO 100 614 116 CONTINUE 615 PSO = .TRUE. 616 ONEPRP = .TRUE. 617 GO TO 100 618 117 CONTINUE 619 SPIDIP = .TRUE. 620 ONEPRP = .TRUE. 621 GO TO 100 622 118 CONTINUE 623 DSO = .TRUE. 624 ONEPRP = .TRUE. 625 GO TO 100 626 119 CONTINUE 627 READ (LUCMD,*) NPQUAD 628 GO TO 100 629 120 CONTINUE 630 READ (LUCMD, *) NPATOM 631 IF (NPATOM .GT. MXCENT) THEN 632 WRITE (LUPRI,'(//2A/A,I8/A,I6)') 633 & 'ERROR: Too many atoms selected under ',WORD, 634 & 'ERROR: Number of atoms selected:',NPATOM, 635 & 'ERROR: Number of atoms allowed: ',MXCENT 636 CALL QUIT('Error in '//WORD1) 637 END IF 638 READ (LUCMD, *) (IPATOM(I),I=1,NPATOM) 639 ALLATM = .FALSE. 640 GO TO 100 641 121 CONTINUE 642 TRIANG = .FALSE. 643 GO TO 100 644 122 CONTINUE 645 SDFC = .TRUE. 646 ONEPRP = .TRUE. 647 GO TO 100 648 123 CONTINUE 649 PROPRI = .TRUE. 650 GO TO 100 651 124 CONTINUE 652 HDO = .TRUE. 653 ONEPRP = .TRUE. 654 GO TO 100 655 125 CONTINUE 656 S1MAG = .TRUE. 657 ONEPRP = .TRUE. 658 GO TO 100 659 126 CONTINUE 660 S2MAG = .TRUE. 661 ONEPRP = .TRUE. 662 GOTO 100 663 127 CONTINUE 664 ANGMOM = .TRUE. 665 ONEPRP = .TRUE. 666 GO TO 100 667 128 CONTINUE 668 ANGLON = .TRUE. 669 ONEPRP = .TRUE. 670 GO TO 100 671 129 CONTINUE 672 LONMOM = .TRUE. 673 ONEPRP = .TRUE. 674 GO TO 100 675 130 CONTINUE 676 MAGMOM = .TRUE. 677 ONEPRP = .TRUE. 678 GO TO 100 679 131 CONTINUE 680 S1MAG = .TRUE. 681 S1MAGT = .TRUE. 682 ONEPRP = .TRUE. 683 GOTO 100 684 132 CONTINUE 685 MGMOMT = .TRUE. 686 LONMOM = .TRUE. 687 HAMILT = .TRUE. 688 ONEPRP = .TRUE. 689 GOTO 100 690 133 CONTINUE 691 KINENE = .TRUE. 692 ONEPRP = .TRUE. 693 GOTO 100 694 134 CONTINUE 695 S2MAG = .TRUE. 696 S2MAGT = .TRUE. 697 ONEPRP = .TRUE. 698 GOTO 100 699 135 CONTINUE 700 DSUSNL = .TRUE. 701 ONEPRP = .TRUE. 702 GOTO 100 703 136 CONTINUE 704 DSUSLL = .TRUE. 705 ONEPRP = .TRUE. 706 GOTO 100 707 137 CONTINUE 708 DSUSLH = .TRUE. 709 ONEPRP = .TRUE. 710 GOTO 100 711 138 CONTINUE 712 DIASUS = .TRUE. 713 ONEPRP = .TRUE. 714 GOTO 100 715 139 CONTINUE 716 DSUTST = .TRUE. 717 DSUSLL = .TRUE. 718 ANGLON = .TRUE. 719 ONEPRP = .TRUE. 720 GOTO 100 721 140 CONTINUE 722 NUCSNL = .TRUE. 723 ONEPRP = .TRUE. 724 GOTO 100 725 141 CONTINUE 726 NUCSLO = .TRUE. 727 ONEPRP = .TRUE. 728 GOTO 100 729 142 CONTINUE 730 NUCSHI = .TRUE. 731 ONEPRP = .TRUE. 732 GOTO 100 733 143 CONTINUE 734 NSNLTS = .TRUE. 735 NUCSLO = .TRUE. 736 PSO = .TRUE. 737 ONEPRP = .TRUE. 738 GOTO 100 739 144 CONTINUE 740 NSLTST = .TRUE. 741 NELFLD = .TRUE. 742 NUCSNL = .TRUE. 743 ONEPRP = .TRUE. 744 GOTO 100 745 145 CONTINUE 746 NELFLD = .TRUE. 747 ONEPRP = .TRUE. 748 GOTO 100 749 146 CONTINUE 750 NSTTST = .TRUE. 751 NUCSLO = .TRUE. 752 NUCSNL = .TRUE. 753 NUCSHI = .TRUE. 754 ONEPRP = .TRUE. 755 GOTO 100 756 147 CONTINUE 757 EFGCAR = .TRUE. 758 ONEPRP = .TRUE. 759 GOTO 100 760 148 CONTINUE 761 EFGSPH = .TRUE. 762 ONEPRP = .TRUE. 763 GOTO 100 764 149 CONTINUE 765 S1MAGL = .TRUE. 766 ONEPRP = .TRUE. 767 GOTO 100 768 150 CONTINUE 769 S1MAGR = .TRUE. 770 ONEPRP = .TRUE. 771 GOTO 100 772 151 CONTINUE 773 HDOBR = .TRUE. 774 ONEPRP = .TRUE. 775 GOTO 100 776 152 CONTINUE 777 S1MLT = .TRUE. 778 S1MAGL = .TRUE. 779 ONEPRP = .TRUE. 780 GOTO 100 781 153 CONTINUE 782 HDOBR = .TRUE. 783 HDOBRT = .TRUE. 784 DIPVEL = .TRUE. 785 ONEPRP = .TRUE. 786 GOTO 100 787 154 CONTINUE 788 S1MRT = .TRUE. 789 S1MAGR = .TRUE. 790 ONEPRP = .TRUE. 791 GOTO 100 792 155 CONTINUE 793 NUCPOT = .TRUE. 794 ONEPRP = .TRUE. 795 GOTO 100 796 156 CONTINUE 797 NUCPOT = .TRUE. 798 HAMILT = .TRUE. 799 NPOTST = .TRUE. 800 ONEPRP = .TRUE. 801 GOTO 100 802 157 CONTINUE 803 MGMO2T = .TRUE. 804 ONEPRP = .TRUE. 805 GOTO 100 806 158 CONTINUE 807 READ (LUCMD,*) PRTHRS 808 GO TO 100 809 159 CONTINUE 810 HBDO = .TRUE. 811 ONEPRP = .TRUE. 812 GOTO 100 813 160 CONTINUE 814 SUSCGO = .TRUE. 815 ONEPRP = .TRUE. 816 GOTO 100 817 161 CONTINUE 818 NSTCGO = .TRUE. 819 ONEPRP = .TRUE. 820 GOTO 100 821 162 CONTINUE 822 EXPIKR = .TRUE. 823 ONEPRP = .TRUE. 824 READ (LUCMD,*) (EXPKR(I),I=1,3) 825 GOTO 100 826 163 CONTINUE 827 MASSVL = .TRUE. 828 ONEPRP = .TRUE. 829 GOTO 100 830 164 CONTINUE 831 DARWIN = .TRUE. 832 ONEPRP = .TRUE. 833 GOTO 100 834 165 CONTINUE 835 READ (LUCMD,'(A7)') FIELD1 836 IF (.NOT. ((FIELD1 .EQ. 'X-FIELD') 837 & .OR. (FIELD1 .EQ. 'Y-FIELD') 838 & .OR. (FIELD1 .EQ. 'Z-FIELD') 839 & .OR. (FIELD1 .EQ. 'XYZ-ALL'))) THEN 840 WRITE (LUPRI,'(/,3A,/)') ' Field direction "',FIELD1, 841 & '" illegal' 842 CALL QUIT('Illegal field directions for CM-1 integrals') 843 END IF 844 CM1 = .TRUE. 845 ONEPRP = .TRUE. 846 GOTO 100 847 166 CONTINUE 848 READ (LUCMD,'(A7)') FIELD2 849C Modified by Bin Gao, December 14, 2009 850C adding XYZ-ALL 851 IF (.NOT. ((FIELD2 .EQ. 'X-FIELD') 852 & .OR. (FIELD2 .EQ. 'Y-FIELD') 853 & .OR. (FIELD2 .EQ. 'Z-FIELD') 854 & .OR. (FIELD2 .EQ. 'XYZ-ALL'))) THEN 855 WRITE (LUPRI,'(/,3A,/)') ' Field direction "',FIELD2, 856 & '" illegal' 857 CALL QUIT('Illegal field directions for CM-2 integrals') 858 END IF 859 CM2 = .TRUE. 860 ONEPRP = .TRUE. 861 GOTO 100 862 167 CONTINUE ! .SQHDOL - Bra differentiated half-derivative overlap matrix 863 SQHDOL = .TRUE. 864 ONEPRP = .TRUE. 865 GOTO 100 866 168 CONTINUE ! .SQHDOR - Ket differentiated half-derivative overlap matrix 867 SQHDOR = .TRUE. 868 ONEPRP = .TRUE. 869 GOTO 100 870 169 CONTINUE 871 NOTWO = .TRUE. 872 SUPMAT = .FALSE. 873 GO TO 100 874 170 CONTINUE 875 READ (LUCMD, *, IOSTAT=IOS) (GAGORG(I),I=1,3) 876 IF (IOS.NE.0) THEN 877 WRITE(LUPRI,*) 'Error in reading (GAGORG(I),I=1,3)!' 878 WRITE(LUPRI,*) '(GAGORG(I),I=1,3):', (GAGORG(I),I=1,3) 879 CALL QUIT('Error in reading (GAGORG(I),I=1,3)!') 880 ENDIF 881 GO TO 100 882 171 CONTINUE 883 READ (LUCMD, *, IOSTAT=IOS) (DIPORG(I),I=1,3) 884 IF (IOS.NE.0) THEN 885 WRITE(LUPRI,*) 'Error in reading (DIPORG(I),I=1,3)!' 886 WRITE(LUPRI,*) '(DIPORG(I),I=1,3):', (DIPORG(I),I=1,3) 887 CALL QUIT('Error in reading (DIPORG(I),I=1,3)!') 888 ENDIF 889 GO TO 100 890 172 CONTINUE 891c ach 892 no2so=.true. 893 GOTO 100 894 173 CONTINUE 895 S1ELE = .TRUE. 896 ONEPRP = .TRUE. 897 GO TO 100 898 174 CONTINUE 899 S1ELB = .TRUE. 900 ONEPRP = .TRUE. 901 GO TO 100 902 175 CONTINUE 903 ONEELD = .TRUE. 904 ONEPRP = .TRUE. 905 GO TO 100 906 176 CONTINUE 907 THETA = .TRUE. 908 ONEPRP = .TRUE. 909 GO TO 100 910 177 CONTINUE ! .NUCMOD 911 CALL QUIT( 912 & 'ERROR, .NUCMOD has been moved to *MOLBAS input section') 913 GO TO 100 914! .SORT Integrals 915 178 CONTINUE 916 SUPMAT = .FALSE. 917 USE_INTSORT = .TRUE. 918 NEWTRA = .TRUE. 919 GO TO 100 920 179 CONTINUE 921 DPLGRA = .TRUE. 922 ONEPRP = .TRUE. 923 GO TO 100 924 180 CONTINUE 925 QUAGRA = .TRUE. 926 ONEPRP = .TRUE. 927 GO TO 100 928 181 CONTINUE 929 OCTGRA = .TRUE. 930 ONEPRP = .TRUE. 931 GO TO 100 932 182 CONTINUE 933 ROTSTR = .TRUE. 934 ONEPRP = .TRUE. 935 GO TO 100 936 183 CONTINUE 937 THRMOM = .TRUE. 938 ONEPRP = .TRUE. 939 GO TO 100 940 184 CONTINUE 941 SOFLD = .TRUE. 942 ONEPRP = .TRUE. 943 GO TO 100 944 185 CONTINUE 945 SOMM = .TRUE. 946 ONEPRP = .TRUE. 947 GO TO 100 948 186 CONTINUE 949 DEROVL = .TRUE. 950 ONEPRP = .TRUE. 951 GO TO 100 952 187 CONTINUE 953 DERHAM = .TRUE. 954 ONEPRP = .TRUE. 955 GO TO 100 956 188 CONTINUE 957 ELGDIA = .TRUE. 958 ONEPRP = .TRUE. 959 GO TO 100 960 189 CONTINUE 961 ELGDIL = .TRUE. 962 ONEPRP = .TRUE. 963 GO TO 100 964 190 CONTINUE 965 MNF_SO = .TRUE. 966 ONEPRP = .TRUE. 967 GO TO 100 968 191 CONTINUE 969 DPTOVL = .TRUE. 970 ONEPRP = .TRUE. 971 GO TO 100 972 192 CONTINUE 973 DPTPOT = .TRUE. 974 ONEPRP = .TRUE. 975 GO TO 100 976 193 CONTINUE 977 XDDXR3 = .TRUE. 978 ONEPRP = .TRUE. 979 GO TO 100 980 194 CONTINUE 981 AD2DAR = .TRUE. 982 READ (LUCMD, *) DARFAC 983 GO TO 100 984 195 CONTINUE 985 PVIOLA = .TRUE. 986 SPNORB = .TRUE. 987 ONEPRP = .TRUE. 988 GO TO 100 989 196 CONTINUE 990 READ(LUCMD,*) BGWEIN 991 BGWEIL = .TRUE. 992 BGWEIN = D1 - D4*BGWEIN 993 GOTO 100 994 197 CONTINUE 995 ONEPRP = .TRUE. 996 DPTPOT = .TRUE. 997 FINDPT = .TRUE. 998 READ (LUCMD, *) DPTFAC 999 GOTO 100 1000 198 CONTINUE 1001 READ (LUCMD, *) X10FNP 1002 PVFINN = D10 ** X10FNP 1003 GOTO 100 1004 199 CONTINUE 1005 PVPINT = .TRUE. 1006 ONEPRP = .TRUE. 1007 GO TO 100 1008 200 CONTINUE 1009 POTENE = .TRUE. 1010 ONEPRP = .TRUE. 1011 GO TO 100 1012 201 CONTINUE 1013C Modified by Bin Gao, December 14, 2009: adds XYZ-ALL 1014 READ (LUCMD,'(A7)') FIELD3 1015 IF (.NOT. ((FIELD3 .EQ. 'XX-FGRD') 1016 & .OR. (FIELD3 .EQ. 'XY-FGRD') 1017 & .OR. (FIELD3 .EQ. 'XZ-FGRD') 1018 & .OR. (FIELD3 .EQ. 'YY-FGRD') 1019 & .OR. (FIELD3 .EQ. 'YZ-FGRD') 1020 & .OR. (FIELD3 .EQ. 'ZZ-FGRD') 1021 & .OR. (FIELD3 .EQ. 'XYZ-ALL'))) THEN 1022 WRITE (LUPRI,'(/4A/)') 'ERROR: Field direction "',FIELD3, 1023 & '" illegal for input ',WORD 1024 CALL QUIT('Illegal field direction for QDBINT integrals') 1025 END IF 1026 QDBINT = .TRUE. 1027 ONEPRP = .TRUE. 1028 GOTO 100 1029 202 CONTINUE 1030 QDBTST = .TRUE. 1031 THETA = .TRUE. 1032 ONEPRP = .TRUE. 1033 GOTO 100 1034 203 CONTINUE ! .R12 1035C R12 calculation requested with linear r12 terms (WK/UniKA/04-11-2002). 1036 R12CAL = .TRUE. 1037 CARMOM = .TRUE. 1038 IORCAR = MAX(IORCAR,2) 1039 ONEPRP = .TRUE. 1040 GOTO 100 1041 204 CONTINUE 1042C R12 calculation requested with Gaussian-type geminals (WK/UniKA/03-24-2005). 1043 R12CAL = .TRUE. 1044 CARMOM = .TRUE. 1045 IORCAR = MAX(IORCAR,2) 1046 ONEPRP = .TRUE. 1047 R12EOR = .TRUE. 1048C SLATER = .TRUE. 1049 READ(LUCMD,*) NTOGAM 1050 IF (NTOGAM .LE. 0 .OR. NTOGAM .GT. MAXGAM) 1051 & CALL QUIT('Invalid number of Gaussian-damped terms.') 1052 DO NTGAM = 1, NTOGAM 1053 READ(LUCMD,*) GAMAB(NTGAM), GAMAX(NTGAM) 1054 END DO 1055 GOTO 100 1056 205 CONTINUE 1057C Calculation with multiple basis sets requested (WK/UniKA/04-11-2002). 1058 CALL QUIT('ERROR, .AUXBAS has been renamed to .R12AUX'// 1059 & ' and moved to *MOLBAS') 1060 GOTO 100 1061 206 CONTINUE 1062C No calculation of 1/r12 integrals (WK/UniKA/04-11-2002). 1063 NOTV12 = .TRUE. 1064 V12INT = .FALSE. 1065 GOTO 100 1066 207 CONTINUE 1067C Calculation of integrals of the operator r12 (WK/UniKA/04-11-2002). 1068 R12INT = .TRUE. 1069 GOTO 100 1070 208 CONTINUE 1071C Calculation of integrals of the operator U12 (WK/UniKA/04-11-2002). 1072 U12INT = .TRUE. 1073 GOTO 100 1074 209 CONTINUE 1075C Calculation of integrals of the operator U21 (WK/UniKA/04-11-2002). 1076 U21INT = .TRUE. 1077 GOTO 100 1078 210 CONTINUE 1079C Integrals are computed for DIRCCR12 program (WK/UniKA/25-11-2002). 1080 DCCR12 = .TRUE. 1081 SUPMAT = .FALSE. 1082 GOTO 100 1083cLig added what to do for RANGMO and RPSO options 1084 211 CONTINUE 1085 RANGMO = .TRUE. 1086 ONEPRP = .TRUE. 1087 GO TO 100 1088 212 CONTINUE 1089 RPSO = .TRUE. 1090 ONEPRP = .TRUE. 1091 GO TO 100 1092 213 CONTINUE 1093 PXPINT = .TRUE. 1094 ONEPRP = .TRUE. 1095 GO TO 100 1096 214 CONTINUE 1097 NOPICH = .TRUE. 1098 GO TO 100 1099 215 CONTINUE 1100 OZKE = .TRUE. 1101 ONEPRP = .TRUE. 1102 GO TO 100 1103 216 CONTINUE 1104 PSOKE = .TRUE. 1105 ONEPRP = .TRUE. 1106 GO TO 100 1107 217 CONTINUE 1108 DNSKE = .TRUE. 1109 ONEPRP = .TRUE. 1110 GO TO 100 1111 1112 218 CONTINUE 1113 SDKE = .TRUE. 1114 ONEPRP = .TRUE. 1115 GO TO 100 1116 219 CONTINUE 1117 FCKE = .TRUE. 1118 ONEPRP = .TRUE. 1119 GO TO 100 1120 220 CONTINUE 1121 DSOKE = .TRUE. 1122 ONEPRP = .TRUE. 1123 GO TO 100 1124 221 CONTINUE 1125 PSOOZ = .TRUE. 1126 ONEPRP = .TRUE. 1127 GO TO 100 1128 222 CONTINUE ! .SQH2DO - Ket second differentiated half-derivative overlap matrix 1129 SQHD2O = .TRUE. 1130 ONEPRP = .TRUE. 1131 GO TO 100 1132 223 CONTINUE 1133C R12 calculation requested with fitted Slater-type geminal (WK/UniKA/03-24-2005). 1134 R12CAL = .TRUE. 1135 CARMOM = .TRUE. 1136 IORCAR = MAX(IORCAR,2) 1137 ONEPRP = .TRUE. 1138 R12EOR = .TRUE. 1139 SLATER = .TRUE. 1140 READ(LUCMD,*) ALPSTG 1141 READ(LUCMD,*) NTOGAM 1142 IF (NTOGAM .LE. 0 .OR. NTOGAM .GT. MAXGAM) 1143 & CALL QUIT('Invalid number of Gaussians.') 1144 CALL STGFIT(ALPSTG,NTOGAM,GAMAB,GAMAX) 1145 GOTO 100 1146 224 CONTINUE 1147C R12 calculation requested with fitted "r12-times-Slater" geminal (WK/UniKA/03-24-2005). 1148 R12CAL = .TRUE. 1149 CARMOM = .TRUE. 1150 IORCAR = MAX(IORCAR,2) 1151 ONEPRP = .TRUE. 1152 R12EOR = .TRUE. 1153 READ(LUCMD,*) ALPSTG 1154 READ(LUCMD,*) NTOGAM 1155 IF (NTOGAM .LE. 0 .OR. NTOGAM .GT. MAXGAM) 1156 & CALL QUIT('Invalid number of Gaussians.') 1157 CALL RSTGFT(ALPSTG,NTOGAM,GAMAB,GAMAX) 1158 GOTO 100 1159 225 CONTINUE 1160C R12 calculation requested with Gaussian-damped r12 integrals (WK/UniKA/03-24-2005). 1161 R12CAL = .TRUE. 1162 CARMOM = .TRUE. 1163 IORCAR = MAX(IORCAR,2) 1164 ONEPRP = .TRUE. 1165 R12EOR = .TRUE. 1166 READ(LUCMD,*) NTOGAM 1167 IF (NTOGAM .LE. 0 .OR. NTOGAM .GT. MAXGAM) 1168 & CALL QUIT('Invalid number of Gaussians.') 1169 DO NTGAM = 1, NTOGAM 1170 READ(LUCMD,*) GAMAB(NTGAM), GAMAX(NTGAM) 1171 END DO 1172 GOTO 100 1173 226 CONTINUE 1174C R12 calculation requested with fitted ERFC-type geminal (WK/UniKA/03-29-2005). 1175 R12CAL = .TRUE. 1176 CARMOM = .TRUE. 1177 IORCAR = MAX(IORCAR,2) 1178 ONEPRP = .TRUE. 1179 R12EOR = .TRUE. 1180 SLATER = .TRUE. 1181 READ(LUCMD,*) ALPSTG 1182 READ(LUCMD,*) NTOGAM 1183 IF (NTOGAM .LE. 0 .OR. NTOGAM .GT. MAXGAM) 1184 & CALL QUIT('Invalid number of Gaussians.') 1185 CALL ERFFIT(ALPSTG,NTOGAM,GAMAB,GAMAX) 1186 GOTO 100 1187 227 CONTINUE 1188C R12 calculation requested with fitted "r12-times-ERFC" geminal (WK/UniKA/03-29-2005). 1189 R12CAL = .TRUE. 1190 CARMOM = .TRUE. 1191 IORCAR = MAX(IORCAR,2) 1192 ONEPRP = .TRUE. 1193 R12EOR = .TRUE. 1194 READ(LUCMD,*) ALPSTG 1195 READ(LUCMD,*) NTOGAM 1196 IF (NTOGAM .LE. 0 .OR. NTOGAM .GT. MAXGAM) 1197 & CALL QUIT('Invalid number of Gaussians.') 1198 CALL RRFFIT(ALPSTG,NTOGAM,GAMAB,GAMAX) 1199 GOTO 100 1200 228 CONTINUE 1201 DIPLEN=.TRUE. 1202 ONEPRP=.TRUE. 1203 GO TO 100 1204 229 CONTINUE 1205 DAR2EL = .TRUE. 1206 GO TO 100 1207 230 CONTINUE 1208 EFBDER = .TRUE. 1209 ONEPRP = .TRUE. 1210 GO TO 100 1211 231 CONTINUE 1212 EFB2DR = .TRUE. 1213 ONEPRP = .TRUE. 1214 GO TO 100 1215 232 CONTINUE 1216 MAGQDP = .TRUE. 1217 ONEPRP = .TRUE. 1218 GO TO 100 1219 233 CONTINUE 1220 MAGQDP = .TRUE. 1221 ANGMOM = .TRUE. 1222 ONEPRP = .TRUE. 1223 MQDPTS = .TRUE. 1224 GOTO 100 1225 234 CONTINUE 1226 ORBORB = .TRUE. 1227 GOTO 100 1228 235 CONTINUE 1229 DERAM = .TRUE. 1230 ONEPRP = .TRUE. 1231 GO TO 100 1232 236 CONTINUE 1233 DIPANH = .TRUE. 1234 ONEPRP = .TRUE. 1235 GO TO 100 1236! .NEWTRA ("new transformation" from 1988 (!) 1237! this option is without calling SORTAO 1238 237 CONTINUE 1239 NEWTRA = .TRUE. 1240 GO TO 100 1241 238 CONTINUE 1242 KINADI = .TRUE. 1243 ONEPRP = .TRUE. 1244 GOTO 100 1245 239 CONTINUE 1246#if defined(BUILD_GEN1INT) 1247 TEST_GEN1INT = .true. 1248#endif 1249 goto 100 1250C copied from linsca abacus/herdrv.F, Bin Gao, December 17, 2009 1251 240 CONTINUE 1252 S2MBRA = .TRUE. 1253 ONEPRP = .TRUE. 1254 GO TO 100 1255 241 CONTINUE 1256 S2MKET = .TRUE. 1257 ONEPRP = .TRUE. 1258 GO TO 100 1259 242 CONTINUE 1260 S2MMIX = .TRUE. 1261 ONEPRP = .TRUE. 1262 GO TO 100 1263 243 CONTINUE 1264c jim-dbg : setting all integrals for LRESC module 1265 ONEPRP = .TRUE. 1266 DOLRINTS=.TRUE. 1267 FERMI = .TRUE. ! 115.Fermi 1268 KINENE = .TRUE. ! 133.Kinenerg 1269 SOFLD = .TRUE. ! 184.somf 1270 DPTOVL = .TRUE. ! 191. 1271 SPIDIP = .TRUE. ! 117.Sd 1272 ANGMOM = .TRUE. ! 127. 1273 DARWIN = .TRUE. ! 164 1274 MASSVL = .TRUE. ! 163 1275 NSTCGO = .TRUE. ! 161.Dia 1276 PSO = .TRUE. ! 116.Pso 1277 DIPVEL = .TRUE. ! 108.DipVel 1278 OZKE = .TRUE. ! 215.Lkin 1279 PSOKE = .TRUE. ! 216.PsoKin 1280 DNSKE = .TRUE. ! 217.DiaKin 1281 PSOOZ = .TRUE. ! 221. PSO-OZ for AngPso 1282Cxx el problema es que el origen de gauge debe estar en el atomo q elegimos, LRATOM. 1283c 1 - Leer INP file y buscar LRATOM 1284c 2 - hacer que GAGORG = CORD(LRATOM) 1285c Por ahora GAGORG = (0.0 0.0 0.0) 1286c ACLARARLO EN OUT FILE !! 1287cx write(lupri,*)' At INTEGRAL SECT to do GAUGEO on LRATOM' 1288cx write(lupri,*)' Calling lrscinp for SELECT' 1289cx CALL LRSCINP('.SELECT') 1290cx write(lupri,*)' Despues. LRATOM ', LRATOM 1291cx GAGORG(1) = 0.000000 1292cx GAGORG(2) = 0.000000 1293cx GAGORG(3) = 0.000000 1294c ================================================= 1295c ------------------------------------------------- 1296css x READ (LUCMD,*) (LRGAUG(IS), IS = 1, 3) 1297css x DO ICENT = 1, NUCIND 1298c x NAME = NAMEX(3*ICENT)(1:4) 1299css x WRITE (LUPRI,'(2X,A,3X," : ",3(A1,2X,A,F15.10))') 1300css & x NAME, '1' , 'x' , CORD(1,ICENT), 1301css & x '2' , 'y' , CORD(2,ICENT), 1302css & x '3' , 'z' , CORD(3,ICENT) 1303css x ENDDO 1304cs READ (LUCMD, *, IOSTAT=IOS) (GAGORG(I),I=1,3) 1305cs IF (IOS.NE.0) THEN 1306cs WRITE(LUPRI,*) 'Error in reading (GAGORG(I),I=1,3)!' 1307cs WRITE(LUPRI,*) '(GAGORG(I),I=1,3):', (GAGORG(I),I=1,3) 1308cs CALL QUIT('Error in reading (GAGORG(I),I=1,3)!') 1309cs ENDIF 1310 GO TO 100 1311 244 CONTINUE 1312 RMAOTWO = .FALSE. 1313 GO TO 100 1314 245 CONTINUE ! .TWOBYTEPACKING 1315 TWOBYTEPACKING = .TRUE. 1316 GO TO 100 1317 246 CONTINUE 1318 LFDIPLN = .TRUE. 1319 ONEPRP = .TRUE. 1320 GO TO 100 1321 ELSE IF (PROMPT .EQ. '*') THEN 1322 GO TO 300 1323 ELSE 1324 WRITE (LUPRI,'(/3A/)') ' Prompter "',PROMPT,'" illegal' 1325 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 1326 CALL QUIT('Illegal prompt in '//WORD1) 1327 END IF 1328 300 CONTINUE 1329 1330 IF (.NOT. BGWEIL) BGWEIN = WEINBG 1331cLig <> 1332 NMRISS = FERMI .OR. PSO .OR. SPIDIP .OR. DSO .OR. SDFC 1333 & .OR. RPSO 1334 IF (TSTINP) WRITE (LUPRI,'(/A/)') ' @@ Input test run only !!!' 1335 WRITE (LUPRI,'(/A,I5)') ' Default print level: ',IPRDEF 1336 IF (SPNORB .AND. SOTEST) THEN 1337 WRITE (LUPRI,'(/A)') 1338 * ' FATAL ERROR in input: .SPIN-ORBIT and .SOTEST cannot '/ 1339 * /'both be specified.' 1340 CALL QUIT('Error in **INTEGRALS input.') 1341 END IF 1342 IF (HAMILT) THEN 1343 IF (NOTWO) THEN 1344 WRITE (LUPRI,'(/A)') 1345 & ' Calculation of one-electron Hamiltonian integrals.' 1346 ELSE 1347 WRITE (LUPRI,'(/A)') 1348 & ' Calculation of one- and two-electron Hamiltonian integrals.' 1349 END IF 1350 ELSE 1351 IF (NOTWO) WRITE (LUPRI,'(/A)') 1352 & ' Two-electron integrals not calculated.' 1353 END IF 1354 IF (R12INT) WRITE (LUPRI,'(A)') 1355 & ' Calculation of two-electron integrals over r12.' 1356 IF (U12INT) WRITE (LUPRI,'(A)') 1357 & ' Calculation of two-electron integrals over [T1,r12].' 1358 IF (U21INT) WRITE (LUPRI,'(A)') 1359 & ' Calculation of two-electron integrals over [T2,r12].' 1360 IF (DAR2EL) WRITE (LUPRI,'(/A)') 1361 & ' Calculate two-electron Darwin integrals' 1362 IF (AD2DAR) WRITE (LUPRI,'(/A/A,F18.14/)') 1363 & ' Two-electron Darwin integrals are added to standard', 1364 & ' repulsion integrals with perturbation parameter:',DARFAC 1365 IF (SPNORB .AND. .NOT.no2so) WRITE (LUPRI,'(/A)') 1366 & ' Calculate two-electron spin-orbit integrals' 1367 IF (ORBORB) WRITE (LUPRI,'(/A)') 1368 & ' Calculate two-electron orbit-orbit integrals' 1369 IF (TWOBYTEPACKING .AND. NBASIS .LE. 255) WRITE (LUPRI,'(/A)') 1370 & ' (Two byte packing of indices on files requested.)' 1371 IF (FINDPT) THEN 1372 WRITE (LUPRI,'(/A/A,F18.14/A,F18.14,A/)') 1373 & ' A direct relativistic perturbation is added to', 1374 & ' Hamiltonian and metric with perturbation parameter:',DPTFAC, 1375 & ' (perturbation parameter * c^{-2} = ',DPTFAC*ALPHAC**2,')' 1376 IF (NODTOT .GT. 0) CALL PARQUIT('FINDPT') 1377 IF (DIRCAL) CALL QUIT( 1378 & 'Direct calculation (.DIRECT) is not compatible with .FINDPT') 1379 END IF 1380 IF (ONEPRP) THEN 1381 WRITE (LUPRI,'(/A)') 1382 & ' The following one-electron property integrals are' 1383 & //' calculated as requested:' 1384 WRITE (LUPRI,'(10X,A)') '- overlap integrals' 1385 IF (DIPLEN) WRITE (LUPRI,'(10X,A)') '- dipole length integrals' 1386 IF (DIPVEL) WRITE (LUPRI,'(10X,A)') 1387 & '- dipole velocity integrals' 1388C 1389 IF (QUADRU) WRITE (LUPRI,'(10X,A)') 1390 & '- quadrupole moment integrals' 1391 IF (THETA) WRITE (LUPRI,'(10X,A)') 1392 & '- traceless quadrupole moment integrals' 1393 IF (SECMOM) WRITE (LUPRI,'(10X,A)')'- second moment integrals' 1394C 1395 IF (THRMOM) WRITE (LUPRI,'(10X,A)')'- third moment integrals' 1396C 1397 IF (CARMOM) THEN 1398 IF (IORCAR .GT. 0) THEN 1399 WRITE (LUPRI,'(10X,A,I2,A)') 1400 & '- Cartesian multipole moment integrals of orders', 1401 & ABS(IORCAR),' and lower' 1402 ELSE 1403 WRITE (LUPRI,'(10X,A,I2)') 1404 & '- Cartesian multipole moment integrals of order', 1405 & ABS(IORCAR) 1406 END IF 1407 END IF 1408 IF (SPHMOM) THEN 1409 IF (IORSPH .GT. 0) THEN 1410 WRITE (LUPRI,'(10X,A,I2,A)') 1411 & '- Spherical multipole moment integrals of orders', 1412 & ABS(IORSPH),' and lower' 1413 ELSE 1414 WRITE (LUPRI,'(10X,A,I2)') 1415 & '- Spherical multipole moment integrals of order', 1416 & ABS(IORSPH) 1417 END IF 1418 END IF 1419C 1420 IF (KINENE) WRITE (LUPRI,'(10X,A)') 1421 & '- electronic kinetic energy' 1422 IF (KINADI) WRITE (LUPRI,'(10X,A)') 1423 & '- adiabatic kinetic energy from nuclear derivatives' 1424 IF (MASSVL) WRITE (LUPRI,'(10X,A)') 1425 & '- mass velocity integrals' 1426 IF (DARWIN) WRITE (LUPRI,'(10X,A)') 1427 & '- 1-electron Darwin integrals' 1428 IF (SPNORB) WRITE (LUPRI,'(10X,A)') 1429 & '- spatial spin-orbit integrals' 1430c ach 1431 IF (no2so) WRITE (LUPRI,'(10X,A)') 1432 & 'but no two-electron spin-orbit integrals' 1433 IF (NMRISS) THEN 1434 IF (FERMI) THEN 1435 WRITE (LUPRI,'(10X,A)')'- Fermi contact integrals' 1436 WRITE (LUPRI,'(10X,A)') 1437 & ' (Dirac delta function integrals)' 1438 END IF 1439 IF (PSO) THEN 1440 WRITE (LUPRI,'(10X,A)') 1441 & '- paramagnetic spin-orbit integrals' 1442 WRITE (LUPRI,'(10X,A)') 1443 & ' (nuclear moment - electron orbit coupling)' 1444 END IF 1445 IF (SPIDIP) THEN 1446 WRITE (LUPRI,'(10X,A)')'- spin-dipole integrals' 1447 WRITE (LUPRI,'(10X,A)') 1448 & ' (electron spin - nuclear moment coupling)' 1449 END IF 1450 IF (DSO) THEN 1451 WRITE (LUPRI,'(10X,A)') 1452 & '- diamagnetic spin-orbit integrals' 1453 WRITE (LUPRI,'(10X,A)') 1454 & ' (indirect nuclear dipole - dipole coupling)' 1455 END IF 1456 IF (SDFC) THEN 1457 WRITE (LUPRI,'(10X,A)') 1458 & '- spin-dipole + Fermi contact integrals' 1459 WRITE (LUPRI,'(10X,A)') 1460 & ' (electron spin - nuclear magnetic field coupling)' 1461 END IF 1462cLig 1463 IF (RPSO) THEN 1464 WRITE (LUPRI,'(10X,A)') 1465 & '- CTOCD-DZ diamagnetic shielding integrals' 1466 END IF 1467cLig 1468 END IF 1469 IF (HDO) WRITE (LUPRI,'(10X,A)') 1470 & '- antisymmetrized half-derivative overlap integrals' 1471 IF (S1MAG) WRITE (LUPRI,'(10X,A)') 1472 & '- first magnetic derivatives of overlap integrals' 1473 IF (S1MAGT) WRITE (LUPRI,'(10X,A)') 1474 & '- test of first magnetic derivative of overlap integrals' 1475 IF (S2MAG) WRITE (LUPRI,'(10X,A)') 1476 & '- second magnetic derivatives of overlap integrals' 1477 IF (S2MAGT) WRITE (LUPRI,'(10X,A)') 1478 & '- test of second magnetic derivatives of overlap integrals' 1479 IF (ANGMOM) WRITE (LUPRI,'(10X,A)') 1480 & '- electronic angular momentum around the origin' 1481C copied from linsca abacus/herdrv.F, Bin Gao, December 17, 2009 1482 IF (S2MBRA) WRITE (LUPRI,'(10X,A)') 1483 & '- second magnetic derivatives of overlap integrals, '// 1484 & 'bra part' 1485 IF (S2MKET) WRITE (LUPRI,'(10X,A)') 1486 & '- second magnetic derivatives of overlap integrals, '// 1487 & 'ket part' 1488 IF (S2MMIX) WRITE (LUPRI,'(10X,A)') 1489 & '- second magnetic derivatives of overlap integrals, '// 1490 & 'mixed bra and ket part' 1491 IF (ANGLON) WRITE (LUPRI,'(10X,A)') 1492 & '- electronic angular momentum around the nuclei' 1493 IF (LONMOM) WRITE (LUPRI,'(10X,A)') 1494 & '- London orbital contribution to angular momentum' 1495 IF (MAGMOM) WRITE (LUPRI,'(10X,A)') 1496 & '- one-electron contribution to magnetic moment' 1497 IF (MGMOMT) WRITE (LUPRI,'(10X,A)') 1498 & '- test of London contribution to angular momentum' 1499 IF (DSUSNL) WRITE (LUPRI,'(10X,A)') 1500 & '- Magnetic susceptibility without London orbital contribution' 1501 IF (DSUSLL) WRITE (LUPRI,'(10X,A)') 1502 & '- Angular London orbital contribution to magnetic susc.' 1503 IF (DSUSLH) WRITE (LUPRI,'(10X,A)') 1504 & '- London orbital contribution to magnetic susceptibility' 1505 IF (DIASUS) WRITE (LUPRI,'(10X,A)') 1506 & '- Magnetic susceptibility integrals' 1507 IF (DSUTST) WRITE (LUPRI,'(10X,A)') 1508 & '- Test of London orbital contr. to magnetic susc. integrals' 1509 IF (NUCSNL) WRITE (LUPRI,'(10X,A)') 1510 & '- Nuclear shieldings without London orbital contribution' 1511 IF (NUCSLO) WRITE (LUPRI,'(10X,A)') 1512 & '- London orbital contribution to nuclear shieldings' 1513 IF (NUCSHI) WRITE (LUPRI,'(10X,A)') 1514 & '- Nuclear shielding tensor integrals' 1515 IF (NSNLTS) WRITE (LUPRI,'(10X,A)') 1516 & '- Test of London orbital contribution to nuclear shieldings' 1517 IF (ELGDIA) WRITE (LUPRI,'(10X,A)') 1518 & '- Diamagnetic one-electron spin-orbit (no-London)' 1519 IF (ELGDIL) WRITE (LUPRI,'(10X,A)') 1520 & '- Diamagnetic one-electron spin-orbit (London)' 1521 IF (MNF_SO) WRITE (LUPRI,'(10X,A)') 1522 & '- Mean field spin-orbit integrals (AMFI)' 1523 IF (NELFLD) WRITE (LUPRI,'(10X,A)') 1524 & '- Electric field at the nuclei' 1525 IF (NSNLTS) WRITE(LUPRI,'(10X,A)') 1526 & '- Test of non-London orbital contr. to nuclear shieldings' 1527 IF (NSTTST) WRITE (LUPRI,'(10X,A)') 1528 & '- Test of nuclear shielding tensor integrals' 1529 IF (EFGCAR) WRITE (LUPRI,'(10X,A)') 1530 & '- Cartesien electric field gradient integrals' 1531 IF (EFGSPH) WRITE (LUPRI,'(10X,A)') 1532 & '- Spherical electric field gradient integrals' 1533 IF (S1MAGL) WRITE (LUPRI,'(10X,A)') 1534 & '- Bra-differentiated overlap matrix with respect to B' 1535 IF (S1MAGR) WRITE (LUPRI,'(10X,A)') 1536 & '- Ket-differentiated overlap matrix with respect to B' 1537 IF (HBDO) WRITE (LUPRI,'(10X,A)') 1538 & '-Half B-differentiated overlap matrix' 1539 IF (HDOBR) WRITE (LUPRI,'(10X,A)') 1540 & '- Ket-differentiated hdo-integrals with respect to B' 1541 IF (S1MLT) WRITE (LUPRI,'(10X,A)') 1542 & '- Test of bra-diff. overlap matrix with respect to B' 1543 IF (S1MRT) WRITE (LUPRI,'(10X,A)') 1544 & '- Test of ket-diff. overlap matrix with respect to B' 1545 IF (HDOBRT) WRITE (LUPRI,'(10X,A)') 1546 & '- Test of ket-diff. hdo-integrals with respect to B' 1547 IF (SQHDOL) WRITE (LUPRI,'(10X,A)') 1548 & '- Bra differentiated half-derivative overlap matrix' 1549 IF (SQHDOR) WRITE (LUPRI,'(10X,A)') 1550 & '- Ket differentiated half-derivative overlap matrix' 1551 IF (SQHD2O) WRITE (LUPRI,'(10X,A)') 1552 & '- Ket second differentiated half-derivative overlap matrix' 1553 IF (NUCPOT) WRITE (LUPRI,'(10X,A)') 1554 & '- Potential energy at the nuclei' 1555 IF (NPOTST) WRITE (LUPRI,'(10X,A)') 1556 & '- Test of nuclear potential energy' 1557 IF (MGMO2T) WRITE (LUPRI,'(10X,A)') 1558 & '- Test of two-electron part of magnetic moment' 1559 IF (SUSCGO) WRITE (LUPRI,'(10X,A)') 1560 & '- Diamagnetic magnetizability using common gauge origin' 1561 IF (NSTCGO) WRITE (LUPRI,'(10X,A)') 1562 & '- Diamagnetic shielding tensor using common gauge origin' 1563 IF (EXPIKR) WRITE (LUPRI,'(10X,A)') 1564 & '- Cosine and sine integals' 1565 IF (DPLGRA) WRITE (LUPRI,'(10X,A)') 1566 & '- Dipole gradient integrals' 1567 IF (DIPANH) WRITE (LUPRI,'(10X,A)') 1568 & '- Anharmonic corrections to dipole gradients' 1569 IF (QUAGRA) WRITE (LUPRI,'(10X,A)') 1570 & '- Quadrupole gradient integrals' 1571 IF (OCTGRA) WRITE (LUPRI,'(10X,A)') 1572 & '- Octupole gradient integrals' 1573 IF (ROTSTR) WRITE (LUPRI,'(10X,A)') 1574 & '- Rotational strength integrals' 1575 IF (SOFLD) WRITE (LUPRI,'(10X,A)') 1576 & '- Magnetic-field correction to one-electron SO integrals' 1577 IF (SOMM) WRITE (LUPRI,'(10X,A)') 1578 & '- Magnetic-moment correction to one-electron SO integrals' 1579 IF (POTENE) WRITE (LUPRI,'(10X,A)') 1580 & '- Potential energy Hamiltonian integrals' 1581 IF (PVPINT) WRITE (LUPRI,'(10X,A)') 1582 & '- pVp integrals of the Douglas-Kroll Hamiltonian' 1583 IF (PXPINT) WRITE (LUPRI,'(10X,A)') 1584 & '- small component dipole length integrals for DPT' 1585 IF (DKTRAN) WRITE (LUPRI,'(10X,A)') 1586 & '- The second order Douglas-Kroll-Hess transformation'// 1587 & ' will be applied (DKH2)' 1588 IF (CM1) THEN 1589 WRITE (LUPRI,'(10X,A)') 1590 & '- First order magnetic derivative of electric field' 1591 WRITE (LUPRI,'(12X,A,A1,A)') 1592 & 'Electric field applied in ',FIELD1(1:1),'-direction' 1593 END IF 1594 IF (CM2) THEN 1595 WRITE (LUPRI,'(10X,A)') 1596 & '- Second order magnetic derivative of electric field' 1597 WRITE (LUPRI,'(12X,A,A1,A)') 1598 & 'Electric field applied in ',FIELD2(1:1),'-direction' 1599 END IF 1600 IF (QDBINT) THEN 1601 WRITE (LUPRI,'(10X,A)') 1602 & '- First order magnetic derivative of electric '// 1603 & 'field gradient' 1604Cajt qdb WRITE (LUPRI,'(12X,A,A1,A)') 1605Cajt qdb & 'Electric field gradient applied in ',FIELD3(1:2), 1606 WRITE (LUPRI,'(12X,A,A7,A)') 1607 & 'Electric field gradient applied in ',FIELD3(1:7), 1608Cajt qdb 1609 & '-direction' 1610 END IF 1611 IF (EFBDER) WRITE (LUPRI,'(10X,A)') 1612 & '- 1.order London orbital correction to electric field at' 1613 & //' a position in space' 1614 IF (EFB2DR) WRITE (LUPRI,'(10X,A)') 1615 & '- 2.order London orbital correction to electric field at' 1616 & //' a position in space' 1617 IF (QDBTST) WRITE (LUPRI,'(10X,A)') 1618 & 'Magnetic-field derivative integrals of electric '// 1619 & 'field gradients will be tested' 1620 IF (DEROVL) WRITE (LUPRI,'(10X,A)') 1621 & '- Geometrical derivatives of overlap integrals' 1622 IF (DERAM) WRITE (LUPRI,'(10X,A)') 1623 & '- Geometrical derivatives of angular momentum integrals' 1624 IF (DERHAM) WRITE (LUPRI,'(10X,A)') 1625 & '- Geometrical derivatives of one-electron Hamiltonian '// 1626 & 'integrals' 1627 IF (MAGQDP) WRITE (LUPRI,'(10X,A)') 1628 & '- Magnetic quadrupole integrals calculated' 1629 IF (MQDPTS) WRITE (LUPRI,'(10X,A)') 1630 & '- Test of magnetic quadrupole integrals calculated' 1631 IF (PSOKE) WRITE (LUPRI,'(10X,A)') 1632 & '- kinetic energy correction to paramagnetic spin-orbit' 1633 IF (DNSKE) WRITE (LUPRI,'(10X,A)') 1634 & '- kinetic energy correction to diamagnetic shielding' 1635 IF (OZKE) WRITE (LUPRI,'(10X,A)') 1636 & '- kinetic energy correction to orbital Zeeman' 1637 IF (SDKE) WRITE (LUPRI,'(10X,A)') 1638 & '- kinetic energy correction to spin-dipole operator' 1639 IF (FCKE) WRITE (LUPRI,'(10X,A)') 1640 & '- kinetic energy correction to Fermi Contact operator' 1641 IF (DSOKE) WRITE (LUPRI,'(10X,A)') 1642 & '- kinetic energy correction to diamagnetic spin-orbit '// 1643 & 'operator' 1644 IF (PSOOZ) WRITE (LUPRI,'(10X,A)') 1645 & '- orbital Zeeman correction to paramagnetic spin-orbit'// 1646 & ' integral' 1647 IF (DPTPOT) WRITE (LUPRI,'(10X,A,A)') 1648 & '- small component potential energy for DPT' 1649 IF (DPTOVL) WRITE (LUPRI,'(10X,A)') 1650 & '- small component overlap integrals for DPT' 1651 IF (XDDXR3) WRITE (LUPRI,'(10X,A)') 1652 & '-d/dB d/dK < 1/Rk^3 > type of integrals. ' 1653 IF (PROPRI) WRITE (LUPRI,'(/A)') 1654 & ' All one-electron property integrals are printed.' 1655 IF (S1ELE) WRITE (LUPRI,'(10X,A)') 1656 & '- first electric derivatives of overlap integrals,'// 1657 & ' Type A' 1658 IF (S1ELB) WRITE (LUPRI,'(10X,A)') 1659 & '- first electric derivatives of overlap integrals,'// 1660 & ' Type B' 1661 IF (ONEELD) WRITE (LUPRI,'(10X,A)') 1662 & '- first electric derivatives of one-electron'// 1663 & ' Hamiltonian integrals' 1664 IF (PVIOLA) WRITE (LUPRI,'(10X,A)') 1665 & '- parity violation electroweak interaction' 1666cLig added message to display for RPSO and RAGNMO optiion 1667 IF (RANGMO) WRITE (LUPRI,'(10X,A)') 1668 & '- Product between (r - r_go) and l_go for the'// 1669 & ' computation of CTOCD-DZ magnetizability in an'// 1670 & ' analytical way' 1671cLig 1672 IF (LFDIPLN) WRITE (LUPRI,'(10X,A)') 1673 & '- Effective dipole operator - only possible with PE library' 1674 END IF 1675 1676 WRITE (LUPRI,'(4(/A,3F20.12))') 1677 & ' Center of mass (bohr):', (CMXYZ(I),I=1,3), 1678 & ' Operator center (bohr):', (ORIGIN(I),I=1,3), 1679 & ' Gauge origin (bohr):', (GAGORG(I),I=1,3), 1680 & ' Dipole origin (bohr):', (DIPORG(I),I=1,3) 1681 1682 IF (EXPIKR) WRITE (LUPRI,'(/A,3F20.15)') 1683 & ' Wave numbers for exp(ikr):', (EXPKR(I),I=1,3) 1684 IF (SOTEST) WRITE (LUPRI,'(/A)') 1685 * ' Test of spatial spin-orbit integrals.' 1686 IF (.NOT.HAMILT) WRITE (LUPRI,'(/A)') 1687 * ' Ordinary (field-free non-relativistic) Hamiltonian '/ 1688 * /'integrals not calculated.' 1689 IF (MGMO2T) WRITE (LUPRI,'(/A,1P,D10.2)') 1690 & ' Threshold for testing two-electron part of magnetic moment:' 1691 & ,PRTHRS 1692 IF (NMRISS) THEN 1693 IF (ALLATM) THEN 1694 WRITE (LUPRI,'(/2A)') 1695 & ' Integrals for all indirect spin-spin', 1696 & ' coupling and/or shielding tensors are calculated.' 1697 ELSE 1698 WRITE (LUPRI,'(/2A/)') 1699 & ' Indirect spin-spin integrals involving the following', 1700 & ' nuclei are calculated:' 1701 WRITE (LUPRI,'(10X,20I3)') (IPATOM(I),I = 1, NPATOM) 1702 END IF 1703 IF (DSO) THEN 1704 WRITE (LUPRI,'(/2A,I3)') 1705 & ' Number of integration points for diamagnetic', 1706 & ' spin-orbit integrals: ',NPQUAD 1707 IF (.NOT.TRIANG) WRITE (LUPRI,'(A)') 1708 & ' Integrals for symmetry related coupling tensors' 1709 & //' JAB and JBA calculated.' 1710 END IF 1711 END IF 1712C 1713C **** Process input for various program sections ***** 1714C 1715 CALL ER2INI 1716 301 PROMPT = WORD(1:1) 1717 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 1718 READ (LUCMD, '(A7)') WORD 1719 GO TO 301 1720 ELSE IF (PROMPT .EQ. '*') THEN 1721 IF (WORD(1:2) .EQ. '**') GO TO 399 1722 CALL UPCASE(WORD) 1723 DO I = 1, NDIR 1724 IF (WORD .EQ. TABDIR(I)) THEN 1725 GO TO (399,302,303,304,305,306,307,308,309), I 1726 END IF 1727 END DO 1728 WRITE (LUPRI,'(/3A/)') ' Input module ',WORD,' nonexistent.' 1729 CALL PRTAB(NDIR,TABDIR,WORD1//' input modules',LUPRI) 1730 CALL QUIT('Illegal input directive in '//WORD1) 1731 ELSE 1732 WRITE (LUPRI,'(/3A/)') ' Prompter "',PROMPT, 1733 & '" illegal or out of order.' 1734 CALL PRTAB(NDIR,TABDIR,WORD1//' input modules',LUPRI) 1735 CALL QUIT('Error in prompt in '//WORD1//', input section.') 1736 END IF 1737 1738! DATA TABDIR /'*END OF', '*READIN', '*ONEINT', '*TWOINT', ! TABDIR(1:4) 1739! & '*SUPINT', '*ER2INT', '*SORINT', '*DENFIT', ! TABDIR(5:8) 1740! & '*QM3INP'/ ! TABDIR(9) 1741 1742C 302 CALL REAINP(WORD,RELCAL) 1743 302 CALL QUIT('Input ERROR, *READIN has been moved to **DALTON'// 1744 & ' section as *MOLBAS') 1745 GO TO 301 1746 303 CALL HR1INP(WORD) 1747 GO TO 301 1748 304 CALL HR2INP(WORD) 1749 GO TO 301 1750 305 CALL HRSINP(WORD) 1751 GO TO 301 1752 306 CALL ER2INP(WORD) 1753 GO TO 301 1754 307 CALL SRTINP(WORD) 1755 GO TO 301 1756 !Density-fitting input keywords 1757 308 CALL DENSFIT_INP(WORD) 1758 GO TO 301 1759 309 CALL QUIT('ERROR, *QM3 has been moved to **DALTON section') 1760 GO TO 301 1761C 1762 399 CONTINUE 1763C 1764 CALL HR1INP(WORD) 1765 CALL HR2INP(WORD) 1766 CALL HRSINP(WORD) 1767 CALL ER2INP(WORD) 1768 CALL SRTINP(WORD) 1769 1000 CONTINUE 1770C 1771 IF (QMMM) CALL READMM(WORK,LWORK) 1772 1773 CALL SETDCH 1774 CALL GPCLOSE(LUCMD,'KEEP') 1775 IF (THRSUP .LT. 0.0D0) THRSUP = THRS 1776C 1777 IF (DORLM .AND. .NOT. CAVUSR) 1778 & WRITE(LUPRI,'(/A,3F12.6)') ' Cavity center (center of mass):', 1779 & (CAVORG(I), I = 1, 3) 1780C 1781 IF (TESTIN) THEN 1782 CALL QUIT('*** End of input test for Molecule specification') 1783 END IF 1784C 1785 CALL QEXIT('HERINP') 1786 RETURN 1787C 1788 END 1789C /* Deck hr1inp */ 1790 SUBROUTINE HR1INP(WORD) 1791#include "implicit.h" 1792#include "priunit.h" 1793#include "mxcent.h" 1794 PARAMETER (NTABLE = 10) 1795C 1796 LOGICAL SET, NEWDEF 1797 CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 1798C 1799#include "orgcom.h" 1800#include "cbiher.h" 1801#include "cbihr1.h" 1802C 1803 SAVE SET 1804 DATA TABLE /'.SKIP ', '.PRINT ', '.SOLVEN', '.NOT AL', '.CAVORG', 1805 & '.FNMC ', 'XXXXXXX', 'XXXXXXX', 'XXXXXXX', 'XXXXXXX'/ 1806 DATA SET/.FALSE./ 1807C 1808 CALL QENTER('HR1INP') 1809 IF (SET) THEN 1810 IF (WORD .NE. '*END OF' .AND. WORD(1:2) .NE. '**') THEN 1811 969 READ (LUCMD, '(A7)') WORD 1812 CALL UPCASE(WORD) 1813 PROMPT = WORD(1:1) 1814 IF (PROMPT .NE. '*') GO TO 969 1815 END IF 1816 GOTO 199 1817 END IF 1818C 1819 SET = .TRUE. 1820C 1821C Initialize /CBIHR1/ 1822C 1823 RUNONE = .TRUE. 1824 IPRONE = IPRDEF 1825 DORLM = .FALSE. 1826 ALLRLM = .TRUE. 1827 CAVUSR = .FALSE. 1828 FNMC = .FALSE. 1829C 1830 NEWDEF = WORD .EQ. '*ONEINT' 1831 ICHANG = 0 1832 IF (NEWDEF) THEN 1833 WORD1 = WORD 1834 100 CONTINUE 1835 READ (LUCMD, '(A7)') WORD 1836 CALL UPCASE(WORD) 1837 PROMPT = WORD(1:1) 1838 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 1839 GO TO 100 1840 ELSE IF (PROMPT .EQ. '.') THEN 1841 ICHANG = ICHANG + 1 1842 DO 200 I = 1, NTABLE 1843 IF (TABLE(I) .EQ. WORD) THEN 1844 GO TO (1,2,3,4,5,6,7,8,9,10), I 1845 END IF 1846 200 CONTINUE 1847 IF (WORD .EQ. '.OPTION') THEN 1848 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 1849 GO TO 100 1850 END IF 1851 WRITE (LUPRI,'(/,3A,/)') ' Keyword "',WORD, 1852 * '" not recognized in ONEINP.' 1853 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 1854 CALL QUIT('Illegal keyword in ONEINP.') 1855 1 CONTINUE 1856 RUNONE = .FALSE. 1857 GO TO 100 1858 2 CONTINUE 1859 READ (LUCMD,*) IPRONE 1860 IF (IPRONE .EQ. IPRDEF) ICHANG = ICHANG - 1 1861 GO TO 100 1862 3 CONTINUE 1863 DORLM = .TRUE. 1864 READ (LUCMD,*) LMAX 1865 GO TO 100 1866 4 CONTINUE 1867 ALLRLM = .FALSE. 1868 GO TO 100 1869 5 CONTINUE 1870 READ (LUCMD,*) (CAVORG(I),I = 1, 3) 1871 CAVUSR = .TRUE. 1872 GO TO 100 1873 6 CONTINUE 1874 FNMC = .TRUE. 1875 GO TO 100 1876 7 CONTINUE 1877 GO TO 100 1878 8 CONTINUE 1879 GO TO 100 1880 9 CONTINUE 1881 GO TO 100 1882 10 CONTINUE 1883 GO TO 100 1884 ELSE IF (PROMPT .EQ. '*') THEN 1885 GO TO 300 1886 ELSE 1887 WRITE (LUPRI,'(/,3A,/)') ' Prompt "',WORD, 1888 * '" not recognized in ONEINP.' 1889 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 1890 CALL QUIT('Illegal prompt in ONEINP.') 1891 END IF 1892 END IF 1893 300 CONTINUE 1894 IF (ICHANG .EQ. 0) GOTO 199 1895 IF (NEWDEF) THEN 1896 CALL HEADER('Changes of defaults for ONEINP:',1) 1897 IF (.NOT.RUNONE) THEN 1898 WRITE (LUPRI,'(A)') ' No one-electron integrals calculated.' 1899 ELSE 1900 IF (IPRONE .NE. IPRDEF) WRITE (LUPRI,'(A,I5)') 1901 & ' Print level in ONEINT:',IPRONE 1902 END IF 1903 IF (DORLM) THEN 1904 WRITE (LUPRI,'(A/A,I2)') 1905 & ' One-electron RLM integrals calculated.', 1906 & ' Maximum L quantum number: ', LMAX 1907 IF (ALLRLM) THEN 1908 WRITE (LUPRI,'(A)') ' All symmetries saved on file.' 1909 ELSE 1910 WRITE (LUPRI,'(A)') 1911 & ' Only totally symmetric integrals saved on file.' 1912 END IF 1913 IF (CAVUSR) WRITE(LUPRI,'(A,3F15.10)') 1914 & ' User supplied cavity center',(CAVORG(I),I=1,3) 1915 END IF 1916 IF (FNMC) WRITE (LUPRI,'(/,2A)') ' Finite nuclear mass', 1917 & ' correction by modifying the kinetic energy integrals.' 1918 END IF 1919 199 CALL QEXIT('HR1INP') 1920 RETURN 1921 END 1922C /* Deck hr2inp */ 1923 SUBROUTINE HR2INP(WORD) 1924 1925#include "implicit.h" 1926#include "priunit.h" 1927#include "mxcent.h" 1928 PARAMETER (D1 = 1.0D0,D0 = 0.0D0) 1929 PARAMETER (NTABLE = 16) 1930C 1931 LOGICAL SET, NEWDEF 1932 CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 1933C 1934#include "cbiher.h" 1935#include "cbihr2.h" 1936#ifdef PRG_DIRAC 1937#include "dcbgen.h" 1938#else 1939#include "gnrinf.h" 1940#endif 1941 SAVE SET 1942 DATA TABLE /'.SKIP ', '.PRINT ', '.PANAS ', '.RETURN', '.SOFOCK', 1943 & '.TIME ', '.ICEDIF', '.IFTHRS', '.THRFAC', '.ERF ', 1944 & '.ERFEXP', '.DOSRIN', '.KSPT2M', '.ERFGAU', '.COMLAM', 1945 & 'xXXXXXX'/ 1946 DATA SET/.FALSE./ 1947C 1948 CALL QENTER('HR2INP') 1949 IF (SET) THEN 1950 IF (WORD .NE. '*END OF' .AND. WORD(1:2) .NE. '**') THEN 1951 969 READ (LUCMD, '(A7)') WORD 1952 CALL UPCASE(WORD) 1953 PROMPT = WORD(1:1) 1954 IF (PROMPT .NE. '*') GO TO 969 1955 END IF 1956 GOTO 199 1957 END IF 1958C 1959 SET = .TRUE. 1960C 1961C Initialize /CBIHR2/ 1962C 1963 RUNTWO = .NOT.NOTWO 1964 IPRTWO = IPRDEF 1965 IPRNTA = 0 1966 IPRNTB = 0 1967 IPRNTC = 0 1968 IPRNTD = 0 1969 RTNTWO = .FALSE. 1970 TKTIME = .FALSE. 1971 SOFOCK = .FALSE. 1972 ICDIFF = 1 1973 IEDIFF = 1 1974 IFTHRS = 20 1975 USRSCR = .FALSE. 1976 THRFAC(1) = D1 1977 THRFAC(2) = D1 1978C 1979C Put ERFEXP in /GNRINF/ to awoid TKTIME conflict with twosta.h 1980C 1981C ERFEXP = .FALSE. 1982C 1983 NEWDEF = WORD .EQ. '*TWOINT' 1984 ICHANG = 0 1985 IPRSUM = 0 1986 IF (NEWDEF) THEN 1987 WORD1 = WORD 1988 100 CONTINUE 1989 READ (LUCMD, '(A7)') WORD 1990 CALL UPCASE(WORD) 1991 PROMPT = WORD(1:1) 1992 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 1993 GO TO 100 1994 ELSE IF (PROMPT .EQ. '.') THEN 1995 ICHANG = ICHANG + 1 1996 DO 200 I = 1, NTABLE 1997 IF (TABLE(I) .EQ. WORD) THEN 1998 GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15), I 1999 END IF 2000 200 CONTINUE 2001 IF (WORD .EQ. '.OPTION') THEN 2002 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 2003 GO TO 100 2004 END IF 2005 WRITE (LUPRI,'(/,3A,/)') ' Keyword "',WORD, 2006 * '" not recognized in HR2INP.' 2007 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 2008 CALL QUIT('Illegal keyword in HR2INP.') 2009 1 CONTINUE 2010 RUNTWO = .FALSE. 2011 GO TO 100 2012 2 CONTINUE 2013 READ (LUCMD, *) IPRTWO, 2014 & IPRNTA, IPRNTB, IPRNTC, IPRNTD 2015 IPRSUM = IPRNTA + IPRNTB + IPRNTC + IPRNTD 2016 IF (IPRTWO .EQ. IPRDEF .AND. IPRSUM .EQ. 0) THEN 2017 ICHANG = ICHANG - 1 2018 END IF 2019 GO TO 100 2020 3 CONTINUE 2021 READ (LUCMD,*,ERR=35) PANAS 2022 GOTO 36 2023 35 PANAS = 0.25D0 2024 BACKSPACE LUCMD 2025 36 CONTINUE 2026C 2027C We cannot use new integral code for Panas correction 2028C 2029 SEGBAS = .FALSE. 2030 GO TO 100 2031 4 CONTINUE 2032 RTNTWO = .TRUE. 2033 GO TO 100 2034 5 CONTINUE 2035C&&&& SOFOCK - construction of Fock matrices in SO-basis 2036 SOFOCK = .TRUE. 2037 GO TO 100 2038 6 CONTINUE 2039 TKTIME = .TRUE. 2040 GO TO 100 2041 7 CONTINUE 2042C&&&& ICEDIF Separate screening of Coulomb and exchange contributions 2043C&&&& in direct SCF 2044 READ (LUCMD,*) ICDIFF,IEDIFF 2045 GO TO 100 2046 8 CONTINUE 2047C&&&& IFTHRS Screening threshold in direct construction of Fock matrices 2048 READ (LUCMD,*) IFTHRS 2049 USRSCR = .TRUE. 2050 GO TO 100 2051 9 CONTINUE 2052C&&& THRFAC: Factors to multiply LL-integral threshold for SL- and SS - integrals 2053C&&& This option only used in DIRAC 2054 READ(LUCMD,*) THRFAC(1),THRFAC(2) 2055 GO TO 100 2056 10 CONTINUE 2057C&&& ERF : Value of \chi in separation of two-electron operator 2058C&&& Vee=Vsr+Vlr ; Vlr=erf(\chi*r_12)/(r_12) 2059 READ (LUCMD,*,ERR=37) CHIVAL 2060 DOSRIN = .TRUE. 2061 GO TO 100 2062 37 CALL QUIT('Error reading CHIVAL for *TWOINT/.ERF') 2063 11 CONTINUE 2064C&&& ERFEXP: Value of \chi in separation of two-electron operator 2065C&&& NEW VERSION March 2010 mnp+hjaj 2066 ERFEXP(0) = .TRUE. 2067 ERFEXP(2) = .TRUE. 2068 READ (LUCMD,*,ERR=39) CHIVAL 2069 DOSRIN = .TRUE. 2070 GO TO 100 2071 39 CALL QUIT('Error reading CHIVAL for *TWOINT/.ERFEXP') 2072 12 CONTINUE 2073C&&& DOSRIN : Calculate short-range 2-elec. integrals needed for 2074C&&& computing the short-range Hartree term (U_sr) in MC-srDFT 2075 DOSRIN = .TRUE. 2076 WRITE(LUPRI,*) 'INFO: .DOSRIN option is obsolete' 2077 ! is set by .ERF* options, which always are needed 2078 ! for DOSRIN = .TRUE. 2079 GO TO 100 2080C&&& KSPT2M 2081 13 CONTINUE 2082 READ (LUCMD,*,ERR=38) CHIVAL 2083 CHI1ST = .TRUE. 2084 RUNTWO = .TRUE. 2085 GO TO 100 2086 38 CALL QUIT('Error reading CHIVAL for *TWOINT/.MU ') 2087 14 CONTINUE 2088C&&& ERFGAU: Value of \chi in separation of two-electron operator 2089C&&& Vee=Vsr+Vlr ; Vlr=erf(\chi*r_12)+N*exp(\chi)/(r_12) 2090 ERFEXP(0) = .TRUE. 2091 ERFEXP(1) = .TRUE. 2092 READ (LUCMD,*,ERR=141) CHIVAL 2093 DOSRIN = .TRUE. 2094 GO TO 100 2095 141 CALL QUIT('Error reading CHIVAL for *TWOINT/.ERFGAU') 2096 15 CONTINUE 2097 2098C K. Sharkas and J. Toulouse beg 2099C COMLAM : COMplement LAMbda: value of lambda for linear separation of the electron-electron interaction: Vee= lambda/r_12 + (1-lambda)/r_12 2100C K. Sharkas, A. Savin, H. J. Aa. Jensen, J. Toulouse, J. Chem. Phys. 137, 044104 (2012) 2101C For now, range separation must also be activated (with very large mu) when using the linear separation. 2102C Example of input: 2103C**INTEGRALS 2104C*TWOINT 2105C.DOSRIN 2106C.COMLAM 2107C 0.25 (value of lambda) 2108C.ERF 2109C 10000000000 (large value of mu) 2110C 2111C See in srdft.F for the associated lambda-dependent density functionals. 2112 2113 COMLAM = .TRUE. 2114 READ (LUCMD,*,ERR=151) VLAMBDA 2115 GO TO 100 2116 151 CALL QUIT('Error reading Lambda for *TWOINT/.COMLAM') 2117C K. Sharkas and J. Toulouse end 2118 16 CONTINUE 2119C&&& not in use 2120 GO TO 100 2121 ELSE IF (PROMPT .EQ. '*') THEN 2122 GO TO 300 2123 ELSE 2124 WRITE (LUPRI,'(/,3A,/)') ' Prompt "',WORD, 2125 & '" not recognized in HR2INP.' 2126 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 2127 CALL QUIT('Illegal prompt in HR2INP.') 2128 END IF 2129 END IF 2130 300 CONTINUE 2131 ICEDIF = ICDIFF + 2*IEDIFF 2132C IF (ICHANG .EQ. 0.AND..NOT.(DIRCAL.OR.CHI1ST)) GOTO 199 2133 IF (NEWDEF) THEN 2134 CALL HEADER('Set-up from HR2INP:',1) 2135 IF (.NOT.(RUNTWO.OR.DIRCAL)) THEN 2136 WRITE (LUPRI,'(A)') ' No two-electron integrals calculated.' 2137 ELSE 2138 WRITE (LUPRI,'(A,I5)') ' Print level in TWOINT:',IPRTWO 2139 IF (IPRSUM .GT. 0) THEN 2140 WRITE (LUPRI,'(A,4I3)') 2141 & ' Extra output for the following shells:', 2142 & IPRNTA, IPRNTB, IPRNTC, IPRNTD 2143 IF (RTNTWO) WRITE (LUPRI,'(A)') 2144 & ' Program will exit TWOINT after these shells.' 2145 END IF 2146 IF (TKTIME) WRITE (LUPRI,'(/,2A)') ' Detailed timing for', 2147 & ' integral calculation will be provided.' 2148 IF (PANAS .NE. 0.0D0) WRITE (LUPRI,'(/,A,F10.5)') 2149 & ' Coulomb integrals screened with a factor of',PANAS 2150 IF (CHIVAL .NE. -1.0D0) THEN 2151 IF (CHI1ST) THEN 2152 DIRCAL = .TRUE. 2153 WRITE (LUPRI,'(2A,E15.5)') 2154 & ' Two-electron integrals calculated', 2155 & ' for KS-PT2 with mu = ',CHIVAL 2156 WRITE (LUPRI,'(A)') 2157 & ' DFT optimization is run .DIRECT' 2158 ELSE 2159 IF (ERFEXP(1)) THEN 2160 WRITE (LUPRI,410) 2161 ELSE IF (ERFEXP(2)) THEN 2162 WRITE (LUPRI,420) 2163 ELSE 2164 WRITE (LUPRI,400) 2165 ENDIF 2166 IF (CHIVAL .GT. 0.1d0) then 2167 WRITE (LUPRI,'(A,F10.5)') 2168 & ' with the coupling parameter : ',CHIVAL 2169 ELSE 2170 WRITE (LUPRI,'(A,1P,G12.5)') 2171 & ' with the coupling parameter : ',CHIVAL 2172 END IF 2173 IF (COMLAM) THEN 2174 WRITE (LUPRI,'(2A,E15.5)') 2175 & ' Two-electron integrals calculated', 2176 & ' for linear coupling with Lambda = ',VLAMBDA 2177 ENDIF 2178 ENDIF 2179 END IF 2180 IF (DIRCAL) THEN 2181 IF(SOFOCK) THEN 2182 WRITE(LUPRI,'(A)') 2183 & ' * Direct calculation of Fock matrices in SO-basis.' 2184 ELSE 2185 WRITE(LUPRI,'(A)') 2186 & ' * Direct calculation of Fock matrices in AO-basis.' 2187 ENDIF 2188 2189 IF (USRSCR) THEN 2190 FCKTHR = -IFTHRS 2191 FCKTHR = 10.0D0**FCKTHR 2192 WRITE(LUPRI,'(2A,1P,D9.2)') 2193 & ' * Screening threshold in direct Fock ', 2194 & 'matrix construction: ',FCKTHR 2195 IF(IFTHRS.GE.16) WRITE(LUPRI,'(4X,A)') 2196 & '---> WARNING : Screening turned off !' 2197 ELSE 2198 WRITE(LUPRI,'(A)') 2199 & ' * Program controlled screening thresholds used for this.' 2200 END IF 2201 IF(ICDIFF.EQ.1) WRITE(LUPRI,'(A)') 2202 & ' * Separate density screening of Coulomb integral batches' 2203 IF(IEDIFF.EQ.1) WRITE(LUPRI,'(A)') 2204 & ' * Separate density screening of exchange integral batches' 2205 ENDIF 2206 IF (RELCAL .AND. (THRFAC(1).NE.D1.OR.THRFAC(2).NE.D1)) THEN 2207 WRITE(LUPRI,'(A,2(/3X,A,1P,D10.3))') 2208 + ' * Threshold factors for omitting integrals:', 2209 + 'SL-integrals: ',THRFAC(1), 2210 + 'SS-integrals: ',THRFAC(2) 2211 END IF 2212 END IF 2213 END IF ! NEWDEF 2214 400 FORMAT(/' srDFT-hybrid : Using an Erf type two-elec. operator') 2215 410 FORMAT(/' srDFT-hybrid : Using an Erf+Gau type two-elec. operator' 2216 & ) 2217 420 FORMAT(/' srDFT-hybrid : Using an Erf+Exp type two-elec. operator' 2218 & ) 2219 199 CALL QEXIT('HR2INP') 2220 RETURN 2221 END 2222#ifndef PRG_DIRAC 2223C /* Deck hrsinp */ 2224 SUBROUTINE HRSINP(WORD) 2225C 2226C Input for FRMSUP routine (FoRM SUPermatrix 2-el integrals) 2227C 2228#include "implicit.h" 2229#include "priunit.h" 2230#include "mxcent.h" 2231 PARAMETER (NTABLE = 10) 2232C 2233 LOGICAL SET, NEWDEF 2234 CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 2235C 2236#include "gnrinf.h" 2237#include "cbiher.h" 2238#include "cbihrs.h" 2239 SAVE SET 2240 DATA TABLE /'.SKIP ', '.PRINT ', '.NOSYMM', '.ONLY J', '.THRESH', 2241 & 'XXXXXXX', 'XXXXXXX', 'XXXXXXX', 'XXXXXXX', 'XXXXXXX'/ 2242 DATA SET/.FALSE./ 2243C 2244 CALL QENTER('HRSINP') 2245 IF (SET) THEN 2246 IF (WORD .NE. '*END OF' .AND. WORD(1:2) .NE. '**') THEN 2247 969 READ (LUCMD, '(A7)') WORD 2248 CALL UPCASE(WORD) 2249 PROMPT = WORD(1:1) 2250 IF (PROMPT .NE. '*') GO TO 969 2251 END IF 2252 GOTO 199 2253 END IF 2254C 2255 SET = .TRUE. 2256C 2257C Initialize /CBIHRS/ 2258C 2259 RUNSUP = SUPMAT 2260 IPRSUP = IPRDEF 2261 NOSSUP = .FALSE. 2262 ONLY_J = .FALSE. 2263 THRSUP = -1.0D0 2264C 2265 NEWDEF = WORD .EQ. '*SUPINT' 2266 ICHANG = 0 2267 IF (NEWDEF) THEN 2268 WORD1 = WORD 2269 100 CONTINUE 2270 READ (LUCMD, '(A7)') WORD 2271 CALL UPCASE(WORD) 2272 PROMPT = WORD(1:1) 2273 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 2274 GO TO 100 2275 ELSE IF (PROMPT .EQ. '.') THEN 2276 ICHANG = ICHANG + 1 2277 DO 200 I = 1, NTABLE 2278 IF (TABLE(I) .EQ. WORD) THEN 2279 GO TO (1,2,3,4,5,6,7,8,9,10), I 2280 END IF 2281 200 CONTINUE 2282 IF (WORD .EQ. '.OPTION') THEN 2283 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 2284 GO TO 100 2285 END IF 2286 WRITE (LUPRI,'(/,3A,/)') ' Keyword "',WORD, 2287 * '" not recognized in SUPINP.' 2288 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 2289 CALL QUIT('Illegal keyword in SUPINP.') 2290 1 CONTINUE 2291 RUNSUP = .FALSE. 2292 SUPMAT = .FALSE. 2293 GO TO 100 2294 2 CONTINUE 2295 READ (LUCMD,*) IPRSUP 2296 IF (IPRSUP .EQ. IPRDEF) ICHANG = ICHANG - 1 2297 GO TO 100 2298 3 CONTINUE 2299 NOSSUP = .TRUE. 2300 GO TO 100 2301 4 CONTINUE 2302 ONLY_J = .TRUE. 2303 GO TO 100 2304 5 CONTINUE 2305 READ (LUCMD,*) THRSUP 2306 IF (THR_REDFAC .GT. 0.0D0) THEN 2307 ICHANG = ICHANG + 1 2308 WRITE (LUPRI,'(3A,1P,D10.2)') '@ INFO ',WORD1, 2309 & ' thresholds multiplied with general factor',THR_REDFAC 2310 THRSUP = THRSUP*THR_REDFAC 2311 END IF 2312 GO TO 100 2313 6 CONTINUE 2314 GO TO 100 2315 7 CONTINUE 2316 GO TO 100 2317 8 CONTINUE 2318 GO TO 100 2319 9 CONTINUE 2320 GO TO 100 2321 10 CONTINUE 2322 GO TO 100 2323 ELSE IF (PROMPT .EQ. '*') THEN 2324 GO TO 300 2325 ELSE 2326 WRITE (LUPRI,'(/,3A,/)') ' Prompt "',WORD, 2327 * '" not recognized in SUPINP.' 2328 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 2329 CALL QUIT('Illegal prompt in SUPINP.') 2330 END IF 2331 END IF 2332 300 CONTINUE 2333 IF (ICHANG .EQ. 0) GOTO 199 2334 IF (NEWDEF) THEN 2335 CALL HEADER('Changes of defaults for *SUPINT:',1) 2336 IF (IPRSUP .NE. IPRDEF) WRITE (LUPRI,'(A,I5)') 2337 & ' Print level :',IPRSUP 2338 IF (ONLY_J) THEN 2339 WRITE(LUPRI,'(A)') ' Only Coulomb (no exchange)' 2340 ELSE 2341 WRITE(LUPRI,'(A)') ' P-supermatrix integrals,'// 2342 & ' (if requested for DFT: with scaled exchange)' 2343 END IF 2344 IF (THRSUP .NE. -1.0D0) WRITE (LUPRI,'(A,D12.2)') 2345 & ' Threshold for supermatrix integrals:',THRSUP 2346 IF (NOSSUP) THEN 2347 WRITE (LUPRI,'(A)') ' No symmetry used in SUPINT.' 2348 ELSE 2349 WRITE (LUPRI,'(A)') ' Symmetry used in SUPINT.' 2350 END IF 2351 END IF 2352 199 CALL QEXIT('HRSINP') 2353 RETURN 2354 END 2355C /* Deck herint */ 2356 SUBROUTINE HERINT(WORK,LWORK) 2357#include "implicit.h" 2358#include "priunit.h" 2359#include "iratdef.h" 2360#include "mxcent.h" 2361#include "maxorb.h" 2362#include "aovec.h" 2363#include "maxaqn.h" 2364#include "dummy.h" 2365 PARAMETER (D0 = 0.0D0) 2366C 2367#include "gnrinf.h" 2368#include "cbiher.h" 2369#include "cbihr1.h" 2370#include "cbihr2.h" 2371#include "cbihrs.h" 2372#include "cbieri.h" 2373#include "orgcom.h" 2374#include "nuclei.h" 2375#include "huckel.h" 2376#include "inftap.h" 2377#include "efield.h" 2378#include "r12int.h" 2379#include "drw2el.h" 2380#include "qm3.h" 2381#include "pcmlog.h" 2382C 2383 DIMENSION RLMORI(3) 2384 CHARACTER*8 LABELT(3), LABELS(6) 2385C 2386 DIMENSION WORK(LWORK) 2387#include "memint.h" 2388C 2389C Control routine for calculation of undifferentiated one- and 2390C two-electron Hamiltonian integrals and transformation of 2391C two-electron integrals to P supermatrix elements. 2392C 2393 IF (SKIP) RETURN 2394 CALL QENTER('HERINT') 2395C 2396 IF (IPRDEF .EQ. 1) CALL TITLER('Output from HERINT','*',126) 2397C 2398 I2TYP = 0 2399C 2400C ********************************** 2401C ***** One-Electron Integrals ***** 2402C ********************************** 2403C 2404C Both Hamiltionian and property one-electron integrals. 2405C 2406 IF (RUNONE) THEN 2407 IF (IPRDEF .GT. 1) CALL TITLER('Output from HERONE','*',126) 2408 2409! Open standard property integral file AOPROPER if LUPROP .le. 0 2410! (if calculating integrals for DKH, then LUPROP is already opened 2411! as AODKHINT) 2412 IF (LUPROP .LE. 0) 2413 & CALL GPOPEN(LUPROP,'AOPROPER','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 2414C 2415C ***** These integrals are used to modify the Hamiltonian **** 2416C ***** thus calculated first of all if requested **** 2417C ***** DPT potential energy integrals **** 2418C 2419 IF (ONEPRP .AND. DPTPOT) THEN 2420 CALL TIMER('START ',TIMSTR,TIMEND) 2421 CALL PR1INT('DPTPO1 ',WORK,LWORK,IDUMMY, 2422 & IDUMMY,TRIANG,PROPRI,IPRONE) 2423 CALL FLSHFO(LUPRI) 2424 CALL PR1INT('DPTPO2 ',WORK,LWORK,IDUMMY, 2425 & IDUMMY,TRIANG,PROPRI,IPRONE) 2426 CALL TIMER('DPTPOT',TIMSTR,TIMEND) 2427 CALL FLSHFO(LUPRI) 2428 END IF 2429C 2430C *********************************************** 2431C ***** Overlap, kinetic energy and nuclear ***** 2432C ***** attraction integrals ***** 2433C *********************************************** 2434C 2435C Contract any primitive Douglas-Kroll integrals if necessary 2436C 2437 IF (DKTRAN .AND. .NOT. PVPINT) THEN 2438 CALL DKH_CONT(PROPRI,IPRONE,WORK,LWORK) 2439 END IF 2440C 2441C ***************************************** 2442C ***** R(l,m) integrals for solvent ***** 2443C ***************************************** 2444C 2445 IF (DORLM) THEN 2446 CALL TIMER('START ',TIMSTR,TIMEND) 2447 CALL RLMSET_CBISOL(LMAX) 2448 CALL RLMNUC(WORK,LWORK,IPRONE) 2449 CALL TIMER('RLMNUC',TIMSTR,TIMEND) 2450 CALL DCOPY(3,DIPORG,1,RLMORI,1) 2451 CALL DCOPY(3,CAVORG,1,DIPORG,1) 2452 DO 100 IORDER = 0, LMAX 2453 CALL PR1INT('SOLVENT',WORK,LWORK,IORDER, 2454 & IDUMMY,TRIANG,PROPRI,IPRONE) 2455 100 CONTINUE 2456 CALL DCOPY(3,RLMORI,1,DIPORG,1) 2457 CALL TIMER('RLMINT',TIMSTR,TIMEND) 2458 CALL FLSHFO(LUPRI) 2459 ELSE 2460! some /cbisol/ variables are used for allocation, also when DORLM false /hjaaj 2461 CALL RLMSET_CBISOL(0) 2462 END IF 2463C 2464C ***************************************** 2465C Integrals written on LUONEL - 2466C Some previous information has been written in READIN 2467C NOTE: must be called after CALL RLMSET_CBISOL 2468C because /CBISOL/ is used for allocation in ONEDRV 2469C ***************************************** 2470C 2471 IF (.NOT. ONLYOV) THEN 2472 IF (HAMILT) THEN 2473 CALL TIMER('START ',TIMSTR,TIMEND) 2474 CALL ONEDRV(WORK,LWORK,IPRONE,.FALSE.,0,.FALSE.,.TRUE., 2475 & .TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,PCM) 2476 CALL TIMER('ONEDRV',TIMSTR,TIMEND) 2477 CALL FLSHFO(LUPRI) 2478 END IF 2479 END IF 2480C 2481C ******************************************* 2482C ***** One-electron property integrals ***** 2483C ******************************************* 2484C 2485 IF (DOHUCKEL) THEN 2486C 2487C ***** Overlap integrals for Huckel ***** 2488C 2489 CALL TIMER('START ',TIMSTR,TIMEND) 2490 CALL PR1INT('HUCKEL ',WORK,LWORK,IDUMMY, 2491 & IDUMMY,TRIANG,PROPRI,IPRONE) 2492 CALL TIMER('HUCKEL',TIMSTR,TIMEND) 2493 CALL FLSHFO(LUPRI) 2494 END IF 2495 END IF ! RUNONE 2496 IF (RUNONE .OR. ONEPRP) THEN 2497C 2498C ***** Overlap integrals ***** 2499C 2500! HJAaJ Jan 2011: Always calculate overlap integrals 2501 CALL TIMER('START ',TIMSTR,TIMEND) 2502 CALL PR1INT('OVERLAP ',WORK,LWORK,IDUMMY, 2503 & IDUMMY,TRIANG,PROPRI,IPRONE) 2504 CALL TIMER('OVERLAP',TIMSTR,TIMEND) 2505 CALL FLSHFO(LUPRI) 2506C 2507 IF (ONLYOV) THEN 2508 CALL QUIT('Overlap integrals calculated') 2509 END IF 2510C 2511C ***** Modified overlap integrals ***** 2512C for population analysis /hjaaj Mar 2006 2513C 2514 CALL TIMER('START ',TIMSTR,TIMEND) 2515 CALL PR1INT('POPOVLP',WORK,LWORK,IDUMMY, 2516 & IDUMMY,TRIANG,PROPRI,IPRONE) 2517 CALL TIMER('POPOVLP',TIMSTR,TIMEND) 2518 CALL FLSHFO(LUPRI) 2519C 2520C R**2 and R**4 integrals /hjaaj Nov 2017 2521C 2522 CALL TIMER('START ',TIMSTR,TIMEND) 2523 IORDER = 2 2524 CALL PR1INT('R2 ',WORK,LWORK,IORDER, 2525 & IDUMMY,TRIANG,PROPRI,IPRONE) 2526 CALL TIMER('R2 ',TIMSTR,TIMEND) 2527 CALL FLSHFO(LUPRI) 2528 CALL TIMER('START ',TIMSTR,TIMEND) 2529 IORDER = 4 2530 CALL PR1INT('R4 ',WORK,LWORK,IORDER, 2531 & IDUMMY,TRIANG,PROPRI,IPRONE) 2532 CALL TIMER('R4 ',TIMSTR,TIMEND) 2533 CALL FLSHFO(LUPRI) 2534C 2535C ***** Dipole length integrals ***** 2536C 2537! HJAaJ Jan 2011: Always calculate dipole and quadrupole integrals 2538! and kinetic energy integrals 2539! (For wave function analysis in sirpop.F and RESPONS) 2540! IF (DIPLEN) THEN 2541 CALL TIMER('START ',TIMSTR,TIMEND) 2542 CALL PR1INT('DIPLEN ',WORK,LWORK,IDUMMY, 2543 & IDUMMY,TRIANG,PROPRI,IPRONE) 2544 CALL TIMER('DIPLEN',TIMSTR,TIMEND) 2545 CALL FLSHFO(LUPRI) 2546! END IF 2547C 2548C ***** Quadrupole integrals ***** 2549C 2550! IF (QUADRU) THEN 2551 CALL TIMER('START ',TIMSTR,TIMEND) 2552 CALL PR1INT('QUADRUP',WORK,LWORK,IDUMMY, 2553 & IDUMMY,TRIANG,PROPRI,IPRONE) 2554 CALL TIMER('QUADRUP',TIMSTR,TIMEND) 2555 CALL FLSHFO(LUPRI) 2556! END IF 2557C 2558C ***** Kinetic energy integrals ***** 2559C 2560! IF (KINENE) THEN 2561 CALL TIMER('START ',TIMSTR,TIMEND) 2562 CALL PR1INT('KINENER',WORK,LWORK,IDUMMY, 2563 & IDUMMY,TRIANG,PROPRI,IPRONE) 2564 CALL TIMER('KINENE',TIMSTR,TIMEND) 2565 CALL FLSHFO(LUPRI) 2566! END IF 2567 2568 IF (KINADI) THEN ! reduced mass diagonal correction 2569 CALL TIMER('START ',TIMSTR,TIMEND) 2570 CALL PR1INT('KINADIA',WORK,LWORK,IDUMMY, 2571 & IDUMMY,TRIANG,PROPRI,IPRONE) 2572 CALL TIMER('KINADI',TIMSTR,TIMEND) 2573 CALL FLSHFO(LUPRI) 2574 END IF 2575 2576C ------------------------------------------------------ 2577 2578 IF (ONEPRP) THEN 2579 2580C ------------------------------------------------------ 2581 2582C 2583C ***** Dipole velocity integrals ***** 2584C 2585 IF (DIPVEL) THEN 2586 CALL TIMER('START ',TIMSTR,TIMEND) 2587 CALL PR1INT('DIPVEL ',WORK,LWORK,IDUMMY, 2588 & IDUMMY,TRIANG,PROPRI,IPRONE) 2589 CALL TIMER('DIPVEL',TIMSTR,TIMEND) 2590 CALL FLSHFO(LUPRI) 2591 END IF 2592C 2593C ***** Effective dipole integrals ***** 2594C - only when polarizable environment turned on 2595C using PE library - otherwise they are just 2596C standard dipole integrals 2597 2598 IF (LFDIPLN) THEN 2599 CALL TIMER('START ',TIMSTR,TIMEND) 2600 CALL PR1INT('LFDIPLN',WORK,LWORK,IDUMMY, 2601 & IDUMMY,TRIANG,PROPRI,IPRONE) 2602 CALL TIMER('LFDIPL',TIMSTR,TIMEND) 2603 CALL FLSHFO(LUPRI) 2604 END IF 2605C 2606C ***** Quadrupole integrals ***** 2607C (Other than 'QUADRUP' 2608C 2609 IF (THETA) THEN 2610 CALL TIMER('START ',TIMSTR,TIMEND) 2611 CALL PR1INT('THETA ',WORK,LWORK,IDUMMY, 2612 & IDUMMY,TRIANG,PROPRI,IPRONE) 2613 CALL TIMER('THETA ',TIMSTR,TIMEND) 2614 CALL FLSHFO(LUPRI) 2615 END IF 2616C 2617C ***** Second moments integrals ***** 2618C 2619 IF (SECMOM) THEN 2620 CALL TIMER('START ',TIMSTR,TIMEND) 2621 CALL PR1INT('SECMOM ',WORK,LWORK,IDUMMY, 2622 & IDUMMY,TRIANG,PROPRI,IPRONE) 2623 CALL TIMER('SECMOM',TIMSTR,TIMEND) 2624 CALL FLSHFO(LUPRI) 2625 END IF 2626C 2627C ***** Third moments integrals ***** 2628C 2629 IF (THRMOM) THEN 2630 CALL TIMER('START ',TIMSTR,TIMEND) 2631 CALL PR1INT('THRMOM ',WORK,LWORK,IDUMMY, 2632 & IDUMMY,TRIANG,PROPRI,IPRONE) 2633 CALL TIMER('THRMOM',TIMSTR,TIMEND) 2634 CALL FLSHFO(LUPRI) 2635 END IF 2636C 2637C ***** Potential energy integrals ***** 2638C 2639 IF (POTENE) THEN 2640 CALL TIMER('START ',TIMSTR,TIMEND) 2641 CALL PR1INT('POTENER',WORK,LWORK,IDUMMY, 2642 & IDUMMY,TRIANG,PROPRI,IPRONE) 2643 CALL TIMER('POTENE',TIMSTR,TIMEND) 2644 CALL FLSHFO(LUPRI) 2645 END IF 2646C 2647C ***** Douglas-Kroll-Hess pVp integrals ***** 2648C 2649 IF (PVPINT) THEN 2650 CALL TIMER('START ',TIMSTR,TIMEND) 2651 CALL PR1INT('PVPINT ',WORK,LWORK,IDUMMY, 2652 & IDUMMY,TRIANG,PROPRI,IPRONE) 2653 CALL TIMER('PVPINT',TIMSTR,TIMEND) 2654 CALL FLSHFO(LUPRI) 2655 END IF 2656C 2657C ***** Mass-velocity integrals ***** 2658C 2659 IF (MASSVL) THEN 2660 CALL TIMER('START ',TIMSTR,TIMEND) 2661 CALL PR1INT('MASSVEL',WORK,LWORK,IDUMMY, 2662 & IDUMMY,TRIANG,PROPRI,IPRONE) 2663 CALL TIMER('MASSVL',TIMSTR,TIMEND) 2664 CALL FLSHFO(LUPRI) 2665 END IF 2666C 2667C ***** Darwin integrals ***** 2668C 2669 IF (DARWIN) THEN 2670 CALL TIMER('START ',TIMSTR,TIMEND) 2671 CALL PR1INT('DARWIN ',WORK,LWORK,IDUMMY, 2672 & IDUMMY,TRIANG,PROPRI,IPRONE) 2673 CALL TIMER('DARWIN',TIMSTR,TIMEND) 2674 CALL FLSHFO(LUPRI) 2675 END IF 2676C 2677C ***** Spatial one-electron spin-orbit integrals ***** 2678C 2679 IF (SPNORB) THEN 2680 CALL TIMER('START ',TIMSTR,TIMEND) 2681 CALL PR1INT('SPNORB ',WORK,LWORK,IDUMMY, 2682 & IDUMMY,TRIANG,PROPRI,IPRONE) 2683 CALL TIMER('SPNORB',TIMSTR,TIMEND) 2684 CALL FLSHFO(LUPRI) 2685 END IF 2686C 2687C ***** Cartesian moments integrals ***** 2688C 2689 IF (CARMOM) THEN 2690 ISTR = 0 2691 IF (IORCAR .LT. 0) ISTR = ABS(IORCAR) 2692 DO 200 IORDER = ISTR, ABS(IORCAR) 2693 CALL TIMER('START ',TIMSTR,TIMEND) 2694 CALL PR1INT('CARMOM ',WORK,LWORK,IORDER, 2695 & IDUMMY,TRIANG,PROPRI,IPRONE) 2696 CALL TIMER('CARMOM',TIMSTR,TIMEND) 2697 CALL FLSHFO(LUPRI) 2698 200 CONTINUE 2699 END IF 2700C 2701C ***** Spherical moments integrals ***** 2702C 2703 IF (SPHMOM) THEN 2704 ISTR = 0 2705 IF (IORSPH .LT. 0) ISTR = ABS(IORSPH) 2706 DO 300 IORDER = ISTR, ABS(IORSPH) 2707 CALL TIMER('START ',TIMSTR,TIMEND) 2708 CALL PR1INT('SPHMOM ',WORK,LWORK,IORDER, 2709 & IDUMMY,TRIANG,PROPRI,IPRONE) 2710 CALL TIMER('SPHMOM',TIMSTR,TIMEND) 2711 CALL FLSHFO(LUPRI) 2712 300 CONTINUE 2713 END IF 2714C 2715C ***** Fermi contact integrals ***** 2716C 2717 IF (FERMI) THEN 2718 CALL TIMER('START ',TIMSTR,TIMEND) 2719 CALL PR1INT('FERMI C',WORK,LWORK,IDUMMY, 2720 & IDUMMY,TRIANG,PROPRI,IPRONE) 2721 CALL TIMER('FERMI ',TIMSTR,TIMEND) 2722 CALL FLSHFO(LUPRI) 2723 END IF 2724C 2725C ***** Paramagnetic spin-orbit integrals ***** 2726C 2727 IF (PSO) THEN 2728 CALL TIMER('START ',TIMSTR,TIMEND) 2729 CALL PR1INT('PSO ',WORK,LWORK,IDUMMY, 2730 & IDUMMY,TRIANG,PROPRI,IPRONE) 2731 CALL TIMER('PSO ',TIMSTR,TIMEND) 2732 CALL FLSHFO(LUPRI) 2733 END IF 2734C 2735C ***** spin-dipole integrals ***** 2736C 2737 IF (SPIDIP) THEN 2738 CALL TIMER('START ',TIMSTR,TIMEND) 2739 CALL PR1INT('SPIN-DI',WORK,LWORK,IDUMMY, 2740 & IDUMMY,TRIANG,PROPRI,IPRONE) 2741 CALL TIMER('SPIN-D',TIMSTR,TIMEND) 2742 CALL FLSHFO(LUPRI) 2743 END IF 2744C 2745C ***** spin-dipole + Fermi contact integrals ***** 2746C 2747 IF (SDFC) THEN 2748 CALL TIMER('START ',TIMSTR,TIMEND) 2749 CALL PR1INT('SDFC ',WORK,LWORK,IDUMMY, 2750 & IDUMMY,TRIANG,PROPRI,IPRONE) 2751 CALL TIMER('SDFC ',TIMSTR,TIMEND) 2752 CALL FLSHFO(LUPRI) 2753 END IF 2754C 2755C ***** diamagnetic spin-orbit integrals ***** 2756C 2757 IF (DSO) THEN 2758 CALL TIMER('START ',TIMSTR,TIMEND) 2759 CALL PR1INT('DSO ',WORK,LWORK,IDUMMY, 2760 & NPQUAD,TRIANG,PROPRI,IPRONE) 2761 CALL TIMER('DSO ',TIMSTR,TIMEND) 2762 CALL FLSHFO(LUPRI) 2763 END IF 2764C 2765C ***** Half-derivative overlap integrals ***** 2766C 2767 IF (HDO) THEN 2768 CALL TIMER('START ',TIMSTR,TIMEND) 2769 CALL PR1INT('HDO ',WORK,LWORK,IDUMMY, 2770 & IDUMMY,TRIANG,PROPRI,IPRONE) 2771 CALL TIMER('HDO ',TIMSTR,TIMEND) 2772 CALL FLSHFO(LUPRI) 2773 END IF 2774 IF (SQHD2O) THEN ! second derivative 2775 CALL TIMER('START ',TIMSTR,TIMEND) 2776 CALL PR1INT('SQHD2OR',WORK,LWORK,IDUMMY, 2777 & IDUMMY,TRIANG,PROPRI,IPRONE) 2778 CALL TIMER('SQHD2OR',TIMSTR,TIMEND) 2779 CALL FLSHFO(LUPRI) 2780 END IF 2781C 2782C ***** Bra-half derivative overlap integrals **** 2783C 2784 IF (SQHDOL) THEN 2785 CALL TIMER('START ',TIMSTR,TIMEND) 2786 CALL PR1INT('SQHDOL ',WORK,LWORK,IDUMMY, 2787 & IDUMMY,TRIANG,PROPRI,IPRONE) 2788 CALL TIMER('SQHDOL',TIMSTR,TIMEND) 2789 CALL FLSHFO(LUPRI) 2790 END IF 2791C 2792C ***** Ket-half derivative overlap integrals ***** 2793C 2794 IF (SQHDOR) THEN 2795 CALL TIMER('START ',TIMSTR,TIMEND) 2796 CALL PR1INT('SQHDOR ',WORK,LWORK,IDUMMY, 2797 & IDUMMY,TRIANG,PROPRI,IPRONE) 2798 CALL TIMER('SQHDOR',TIMSTR,TIMEND) 2799 CALL FLSHFO(LUPRI) 2800 END IF 2801C 2802C ***** First magnetic derivative of overlap integrals ***** 2803C 2804 IF (S1MAG) THEN 2805 CALL TIMER('START ',TIMSTR,TIMEND) 2806 CALL PR1INT('S1MAG ',WORK,LWORK,IDUMMY, 2807 & IDUMMY,TRIANG,PROPRI,IPRONE) 2808 CALL TIMER('S1MAG ',TIMSTR,TIMEND) 2809 CALL FLSHFO(LUPRI) 2810 END IF 2811C 2812C **** Bra-differentiated overlap integrals (with respect to B) **** 2813C 2814 IF (S1MAGL) THEN 2815 CALL TIMER('START ',TIMSTR,TIMEND) 2816 CALL PR1INT('S1MAGL ',WORK,LWORK,IDUMMY, 2817 & IDUMMY,TRIANG,PROPRI,IPRONE) 2818 CALL TIMER('S1MAGL',TIMSTR,TIMEND) 2819 CALL FLSHFO(LUPRI) 2820 END IF 2821C 2822C **** Ket-differentiated overlap integrals (with respect to B) **** 2823C 2824 IF (S1MAGR) THEN 2825 CALL TIMER('START ',TIMSTR,TIMEND) 2826 CALL PR1INT('S1MAGR ',WORK,LWORK,IDUMMY, 2827 & IDUMMY,TRIANG,PROPRI,IPRONE) 2828 CALL TIMER('S1MAGR',TIMSTR,TIMEND) 2829 CALL FLSHFO(LUPRI) 2830 END IF 2831C 2832C ***** Half B-differentiated overlap matrix ***** 2833C 2834 IF (HBDO) THEN 2835 CALL TIMER('START ',TIMSTR,TIMEND) 2836 CALL PR1INT('HBDO ',WORK,LWORK,IDUMMY, 2837 & IDUMMY,TRIANG,PROPRI,IPRONE) 2838 CALL TIMER('HBDO ',TIMSTR,TIMEND) 2839 CALL FLSHFO(LUPRI) 2840 END IF 2841C 2842C ***** Ket-differentiated hdo-integrals with respect to B ***** 2843C 2844 IF (HDOBR) THEN 2845 CALL TIMER('START ',TIMSTR,TIMEND) 2846 CALL PR1INT('HDOBR ',WORK,LWORK,IDUMMY, 2847 & IDUMMY,TRIANG,PROPRI,IPRONE) 2848 CALL TIMER('HDOBR ',TIMSTR,TIMEND) 2849 CALL FLSHFO(LUPRI) 2850 END IF 2851C 2852C *** Test of ket-differentiated hdo-integrals with respect to B *** 2853C 2854 IF (HDOBRT) THEN 2855 LABELT(1) = 'XDIPVEL ' 2856 LABELT(2) = 'YDIPVEL ' 2857 LABELT(3) = 'ZDIPVEL ' 2858 CALL HDBTST(WORK,LWORK,IPRINT,LABELT,NATOM,ORIGIN) 2859 END IF 2860C 2861C ***** Test of first magnetic derivative of overlap integrals ***** 2862C 2863 IF (S1MAGT) THEN 2864 LABELT(1) = 'dS/dBX ' 2865 LABELT(2) = 'dS/dBY ' 2866 LABELT(3) = 'dS/dBZ ' 2867 CALL MG1TST(WORK,LWORK,IPRONE,'OVERLAP ',LABELT,ORIGIN, 2868 & 'SM1 ') 2869 CALL FLSHFO(LUPRI) 2870 END IF 2871C 2872C ***** Test of bra-diff. overlap matrix with respect to B ***** 2873C 2874 IF (S1MLT) THEN 2875 LABELT(1) = 'd<S|/dBX' 2876 LABELT(2) = 'd<S|/dBY' 2877 LABELT(3) = 'd<S|/dBZ' 2878 CALL MG1TST(WORK,LWORK,IPRONE,'OVERLAP ',LABELT,ORIGIN, 2879 & 'S1ML') 2880 CALL FLSHFO(LUPRI) 2881 END IF 2882C 2883C ***** Test of ket-diff. overlap matrix with respect to B ***** 2884C Need to be rewritten for both operator center and gauge origin 2885C K.Ruud, jul.-94 2886C 2887 IF (S1MRT) THEN 2888 WRITE (LUPRI,'(A)') 'This option is currently disabled'// 2889 & ' (conflicting origin definitions)' 2890 CALL QUIT('Unimplemented module S1MRT') 2891 LABELT(1) = 'd|S>/dBX' 2892 LABELT(2) = 'd|S>/dBY' 2893 LABELT(3) = 'd|S>/dBZ' 2894 CALL MG1TST(WORK,LWORK,IPRONE,'OVERLAP ',LABELT,ORIGIN, 2895 & 'S1MR') 2896 CALL FLSHFO(LUPRI) 2897 END IF 2898C 2899C ***** Second magnetic derivatives of overlap integrals ***** 2900C 2901 IF (S2MAG) THEN 2902 CALL TIMER('START ',TIMSTR,TIMEND) 2903 CALL PR1INT('S2MAG ',WORK,LWORK,IDUMMY, 2904 & IDUMMY,TRIANG,PROPRI,IPRONE) 2905 CALL TIMER('S2MAG ',TIMSTR,TIMEND) 2906 CALL FLSHFO(LUPRI) 2907 END IF 2908C 2909C ***** Test of second magnetic derivatives of overlap integrals **** 2910C 2911 IF (S2MAGT) THEN 2912 LABELS(1) = 'dS/dB2XX' 2913 LABELS(2) = 'dS/dB2XY' 2914 LABELS(3) = 'dS/dB2XZ' 2915 LABELS(4) = 'dS/dB2YY' 2916 LABELS(5) = 'dS/dB2YZ' 2917 LABELS(6) = 'dS/dB2ZZ' 2918 CALL MGTST2(WORK,LWORK,IPRONE,'OVERLAP ',LABELS) 2919 CALL FLSHFO(LUPRI) 2920 END IF 2921C 2922C ***** Second magnetic derivatives of overlap integrals ***** 2923Cdj bra-diff. part 2924C 2925 IF (S2MBRA) THEN 2926 CALL TIMER('START ',TIMSTR,TIMEND) 2927 CALL PR1INT('S2MBRA ',WORK,LWORK,IDUMMY, 2928 & IDUMMY,TRIANG,PROPRI,IPRONE) 2929 CALL TIMER('S2MBRA',TIMSTR,TIMEND) 2930 CALL FLSHFO(LUPRI) 2931 END IF 2932C 2933C ***** Second magnetic derivatives of overlap integrals ***** 2934Cdj ket-diff. part 2935C 2936 IF (S2MKET) THEN 2937 CALL TIMER('START ',TIMSTR,TIMEND) 2938 CALL PR1INT('S2MKET ',WORK,LWORK,IDUMMY, 2939 & IDUMMY,TRIANG,PROPRI,IPRONE) 2940 CALL TIMER('S2MKET',TIMSTR,TIMEND) 2941 CALL FLSHFO(LUPRI) 2942 END IF 2943C 2944C ***** Second magnetic derivatives of overlap integrals ***** 2945Cdj mixed bra-diff. ket-diff. part 2946C 2947 IF (S2MMIX) THEN 2948 CALL TIMER('START ',TIMSTR,TIMEND) 2949 CALL PR1INT('S2MMIX ',WORK,LWORK,IDUMMY, 2950 & IDUMMY,TRIANG,PROPRI,IPRONE) 2951 CALL TIMER('S2MMIX',TIMSTR,TIMEND) 2952 CALL FLSHFO(LUPRI) 2953 END IF 2954C 2955C ***** Electronic angular momentum around the origin ***** 2956C 2957 IF (ANGMOM) THEN 2958 CALL TIMER('START ',TIMSTR,TIMEND) 2959 CALL PR1INT('ANGMOM ',WORK,LWORK,IDUMMY, 2960 & IDUMMY,TRIANG,PROPRI,IPRONE) 2961 CALL TIMER('ANGMOM',TIMSTR,TIMEND) 2962 CALL FLSHFO(LUPRI) 2963 END IF 2964C 2965C ***** Electronic angular momentum around the nuclei ***** 2966C 2967 IF (ANGLON) THEN 2968 CALL TIMER('START ',TIMSTR,TIMEND) 2969 CALL PR1INT('ANGLON ',WORK,LWORK,IDUMMY, 2970 & IDUMMY,TRIANG,PROPRI,IPRONE) 2971 CALL TIMER('ANGLON',TIMSTR,TIMEND) 2972 CALL FLSHFO(LUPRI) 2973 END IF 2974C 2975C ***** London orbital contribution to angular momentum ***** 2976C 2977 IF (LONMOM) THEN 2978 CALL TIMER('START ',TIMSTR,TIMEND) 2979 CALL PR1INT('LONMOM ',WORK,LWORK,IDUMMY, 2980 & IDUMMY,TRIANG,PROPRI,IPRONE) 2981 CALL TIMER('LONMOM',TIMSTR,TIMEND) 2982 CALL FLSHFO(LUPRI) 2983 END IF 2984C 2985C ***** One-electron contribution to magnetic moment ***** 2986C 2987 IF (MAGMOM) THEN 2988 CALL TIMER('START ',TIMSTR,TIMEND) 2989 CALL PR1INT('MAGMOM ',WORK,LWORK,IDUMMY, 2990 & IDUMMY,TRIANG,PROPRI,IPRONE) 2991 CALL TIMER('MAGMOM',TIMSTR,TIMEND) 2992 CALL FLSHFO(LUPRI) 2993 END IF 2994C 2995C ***** Test of London contribution to angular momentum ***** 2996C 2997 IF (MGMOMT) THEN 2998 LABELT(1) = 'XLONMOM ' 2999 LABELT(2) = 'YLONMOM ' 3000 LABELT(3) = 'ZLONMOM ' 3001 CALL MG1TST(WORK,LWORK,IPRONE,'ONEHAMIL',LABELT,ORIGIN, 3002 & 'LN ') 3003 CALL FLSHFO(LUPRI) 3004 END IF 3005C 3006C ***** Diamagnetic magnetizability using common gauge origin ***** 3007C 3008 IF (SUSCGO) THEN 3009 CALL TIMER('START ',TIMSTR,TIMEND) 3010 CALL PR1INT('DSUSCGO',WORK,LWORK,IDUMMY, 3011 & IDUMMY,TRIANG,PROPRI,IPRONE) 3012 CALL TIMER('SUSCGO',TIMSTR,TIMEND) 3013 CALL FLSHFO(LUPRI) 3014 END IF 3015C 3016C **** Diamagnetic shielding tensor using common gauge origin **** 3017C 3018 IF (NSTCGO) THEN 3019 CALL TIMER('START ',TIMSTR,TIMEND) 3020 CALL PR1INT('NSTCGO ',WORK,LWORK,IDUMMY, 3021 & IDUMMY,TRIANG,PROPRI,IPRONE) 3022 CALL TIMER('NSTCGO',TIMSTR,TIMEND) 3023 CALL FLSHFO(LUPRI) 3024 END IF 3025C 3026C **** Diamagnetic shielding tensor using common gauge origin **** 3027C 3028 IF (ELGDIA) THEN 3029 CALL TIMER('START ',TIMSTR,TIMEND) 3030 CALL PR1INT('ELGDIAN',WORK,LWORK,IDUMMY, 3031 & IDUMMY,TRIANG,PROPRI,IPRONE) 3032 CALL TIMER('ELGDIA',TIMSTR,TIMEND) 3033 CALL FLSHFO(LUPRI) 3034 END IF 3035C 3036C **** Diamagnetic shielding tensor using common gauge origin **** 3037C 3038 IF (ELGDIL) THEN 3039 CALL TIMER('START ',TIMSTR,TIMEND) 3040 CALL PR1INT('ELGDIAL',WORK,LWORK,IDUMMY, 3041 & IDUMMY,TRIANG,PROPRI,IPRONE) 3042 CALL TIMER('ELGDIL',TIMSTR,TIMEND) 3043 CALL FLSHFO(LUPRI) 3044 END IF 3045C 3046C ** Diamagnetic contribution to magnetizability, no london contribution ** 3047C 3048 IF (DSUSNL) THEN 3049 CALL TIMER('START ',TIMSTR,TIMEND) 3050 CALL PR1INT('DSUSNOL',WORK,LWORK,IDUMMY, 3051 & IDUMMY,TRIANG,PROPRI,IPRONE) 3052 CALL TIMER('DSUSNL',TIMSTR,TIMEND) 3053 CALL FLSHFO(LUPRI) 3054 END IF 3055C 3056C ** Magnetizability, angular part of London orbital contribution ** 3057C 3058 IF (DSUSLL) THEN 3059 CALL TIMER('START ',TIMSTR,TIMEND) 3060 CALL PR1INT('DSUSLAN',WORK,LWORK,IDUMMY, 3061 & IDUMMY,TRIANG,PROPRI,IPRONE) 3062 CALL TIMER('DSUSLL',TIMSTR,TIMEND) 3063 CALL FLSHFO(LUPRI) 3064 END IF 3065C 3066C ***** Magnetizability, London orbital contribution ***** 3067C 3068 IF (DSUSLH) THEN 3069 CALL TIMER('START ',TIMSTR,TIMEND) 3070 CALL PR1INT('DSUSLH ',WORK,LWORK,IDUMMY, 3071 & IDUMMY,TRIANG,PROPRI,IPRONE) 3072 CALL TIMER('DSUSLH',TIMSTR,TIMEND) 3073 CALL FLSHFO(LUPRI) 3074 END IF 3075C 3076C ***** Diamagnetic part of magnetizability ***** 3077C 3078 IF (DIASUS) THEN 3079 CALL TIMER('START ',TIMSTR,TIMEND) 3080 CALL PR1INT('DIASUS ',WORK,LWORK,IDUMMY, 3081 & IDUMMY,TRIANG,PROPRI,IPRONE) 3082 CALL TIMER('DIASUS',TIMSTR,TIMEND) 3083 CALL FLSHFO(LUPRI) 3084 END IF 3085C 3086C **** Test of London orbital contribution to magnetizability **** 3087C 3088 IF (DSUTST) THEN 3089 LABELS(1) = 'XXDSUSLL' 3090 LABELS(2) = 'XYDSUSLL' 3091 LABELS(3) = 'XZDSUSLL' 3092 LABELS(4) = 'YYDSUSLL' 3093 LABELS(5) = 'YZDSUSLL' 3094 LABELS(6) = 'ZZDSUSLL' 3095 LABELT(1) = 'XANGLON ' 3096 LABELT(2) = 'YANGLON ' 3097 LABELT(3) = 'ZANGLON ' 3098 CALL SUSTST(WORK,LWORK,IPRONE,LABELT,LABELS) 3099 END IF 3100C 3101C ***** Kinetic energy correction to orbital Zeeman ***** 3102C 3103 IF (OZKE) THEN 3104 CALL TIMER('START ',TIMSTR,TIMEND) 3105 CALL PR1INT('OZKE ',WORK,LWORK,IDUMMY, 3106 & IDUMMY,TRIANG,PROPRI,IPRONE) 3107 CALL TIMER('OZKE ',TIMSTR,TIMEND) 3108 CALL FLSHFO(LUPRI) 3109 END IF 3110C 3111C ***** Kinetic energy correction to paramagnetic SO integrals ***** 3112C 3113 IF (PSOKE) THEN 3114 CALL TIMER('START ',TIMSTR,TIMEND) 3115 CALL PR1INT('PSOKE ',WORK,LWORK,IDUMMY, 3116 & IDUMMY,TRIANG,PROPRI,IPRONE) 3117 CALL TIMER('PSOKE ',TIMSTR,TIMEND) 3118 CALL FLSHFO(LUPRI) 3119 END IF 3120C 3121C ***** Kinetic energy correction to diamagnetic shielding ***** 3122C 3123 IF (DNSKE) THEN 3124 CALL TIMER('START ',TIMSTR,TIMEND) 3125 CALL PR1INT('DNSKE ',WORK,LWORK,IDUMMY, 3126 & IDUMMY,TRIANG,PROPRI,IPRONE) 3127 CALL TIMER('DNSKE ',TIMSTR,TIMEND) 3128 CALL FLSHFO(LUPRI) 3129 END IF 3130C 3131C ***** Kinetic energy correction to spin-dipole integrals ***** 3132C 3133 IF (SDKE) THEN 3134 CALL TIMER('START ',TIMSTR,TIMEND) 3135 CALL PR1INT('SDKE ',WORK,LWORK,IDUMMY, 3136 & IDUMMY,TRIANG,PROPRI,IPRONE) 3137 CALL TIMER('SDKE ',TIMSTR,TIMEND) 3138 CALL FLSHFO(LUPRI) 3139 END IF 3140C 3141C ***** Kinetic energy correction to Fermi contact integrals ***** 3142C 3143 IF (FCKE) THEN 3144 CALL TIMER('START ',TIMSTR,TIMEND) 3145 CALL PR1INT('FCKE ',WORK,LWORK,IDUMMY, 3146 & IDUMMY,TRIANG,PROPRI,IPRONE) 3147 CALL TIMER('FCKE ',TIMSTR,TIMEND) 3148 CALL FLSHFO(LUPRI) 3149 END IF 3150C 3151C ***** Kinetic energy correction to diamagnetic spin-orbit ***** 3152C 3153 IF (DSOKE) THEN 3154 CALL TIMER('START ',TIMSTR,TIMEND) 3155 CALL PR1INT('DSOKE ',WORK,LWORK,IDUMMY, 3156 & NPQUAD,TRIANG,PROPRI,IPRONE) 3157 CALL TIMER('DSOKE ',TIMSTR,TIMEND) 3158 CALL FLSHFO(LUPRI) 3159 END IF 3160C 3161C **** Magnetic-field corrected spin-orbit **** 3162C 3163 IF (PSOOZ) THEN 3164 CALL TIMER('START ',TIMSTR,TIMEND) 3165 CALL PR1INT('PSOOZ ',WORK,LWORK,IDUMMY, 3166 & IDUMMY,TRIANG,PROPRI,IPRONE) 3167 CALL TIMER('PSO-OZ ',TIMSTR,TIMEND) 3168 END IF 3169C 3170C ***** Potential energy at the individual nuclei ***** 3171C 3172 IF (NUCPOT) THEN 3173 CALL TIMER('START ',TIMSTR,TIMEND) 3174 CALL PR1INT('NUCPOT ',WORK,LWORK,IDUMMY, 3175 & IDUMMY,TRIANG,PROPRI,IPRONE) 3176 CALL TIMER('NUCPOT',TIMSTR,TIMEND) 3177 CALL FLSHFO(LUPRI) 3178 END IF 3179C 3180C ***** Test of potential energy **** 3181C 3182 IF (NPOTST) THEN 3183 CALL NPTST(WORK,LWORK,NATOM) 3184 END IF 3185C 3186C ***** Electric field at the individual nuclei ***** 3187C 3188 IF (NELFLD) THEN 3189 CALL TIMER('START ',TIMSTR,TIMEND) 3190 CALL PR1INT('NEFIELD',WORK,LWORK,IDUMMY, 3191 & IDUMMY,TRIANG,PROPRI,IPRONE) 3192 CALL TIMER('NELFLD',TIMSTR,TIMEND) 3193 CALL FLSHFO(LUPRI) 3194 END IF 3195C 3196C ***** Electric field gradient at the individual nuclei, cartesian ***** 3197C 3198 IF (EFGCAR) THEN 3199 CALL TIMER('START ',TIMSTR,TIMEND) 3200 CALL PR1INT('ELFGRDC',WORK,LWORK,IDUMMY, 3201 & IDUMMY,TRIANG,PROPRI,IPRONE) 3202 CALL TIMER('ELFGRD',TIMSTR,TIMEND) 3203 CALL FLSHFO(LUPRI) 3204 END IF 3205C 3206C ***** Electric field gradient at the individual nuclei, spherical ***** 3207C 3208 IF (EFGSPH) THEN 3209 CALL TIMER('START ',TIMSTR,TIMEND) 3210 CALL PR1INT('ELFGRDS',WORK,LWORK,IDUMMY, 3211 & IDUMMY,TRIANG,PROPRI,IPRONE) 3212 CALL TIMER('ELFGRD',TIMSTR,TIMEND) 3213 CALL FLSHFO(LUPRI) 3214 END IF 3215C 3216C ***** Nuclear shielding without London orbital contribution ***** 3217C 3218 IF (NUCSNL) THEN 3219 CALL TIMER('START ',TIMSTR,TIMEND) 3220 CALL PR1INT('NUCSNLO',WORK,LWORK,IDUMMY, 3221 & IDUMMY,TRIANG,PROPRI,IPRONE) 3222 CALL TIMER('NSTNLO',TIMSTR,TIMEND) 3223 CALL FLSHFO(LUPRI) 3224 END IF 3225C 3226C ***** London orbital contribution to nuclear shielding integrals ***** 3227C 3228 IF (NUCSLO) THEN 3229 CALL TIMER('START ',TIMSTR,TIMEND) 3230 CALL PR1INT('NUCSLO ',WORK,LWORK,IDUMMY, 3231 & IDUMMY,TRIANG,PROPRI,IPRONE) 3232 CALL TIMER('NUCSLO',TIMSTR,TIMEND) 3233 CALL FLSHFO(LUPRI) 3234 END IF 3235C 3236C ***** Nuclear shielding tensor integrals ***** 3237C 3238 IF (NUCSHI) THEN 3239 CALL TIMER('START ',TIMSTR,TIMEND) 3240 CALL PR1INT('NUCSHI ',WORK,LWORK,IDUMMY, 3241 & IDUMMY,TRIANG,PROPRI,IPRONE) 3242 CALL TIMER('NUCSHI',TIMSTR,TIMEND) 3243 CALL FLSHFO(LUPRI) 3244 END IF 3245C 3246C ***** Test of London orbital contribution to nuclear shieldings ***** 3247C 3248 IF (NSNLTS) THEN 3249 CALL NS1TST(WORK,LWORK,IPRINT,DOATOM,NATOM) 3250 END IF 3251C 3252C ***** Test of non-London orbital contribution to nuclear shieldings ***** 3253C 3254 IF (NSLTST) THEN 3255 CALL NSTST2(WORK,LWORK,IPRINT,DOATOM,NATOM) 3256 END IF 3257C 3258C ***** Test of nuclear shielding tensor integrals ***** 3259C 3260 IF (NSTTST) THEN 3261 CALL NSTST3(WORK,LWORK,IPRINT,DOATOM,NATOM) 3262 END IF 3263C 3264C ***** Cosine and sine integrals ***** 3265C 3266 IF (EXPIKR) THEN 3267 CALL TIMER('START ',TIMSTR,TIMEND) 3268 CALL DCOPY(3,ORIGIN,1,RLMORI,1) 3269 CALL DCOPY(3,EXPKR,1,ORIGIN,1) 3270 CALL PR1INT('EXPIKR ',WORK,LWORK,IDUMMY, 3271 & IDUMMY,TRIANG,PROPRI,IPRONE) 3272 CALL DCOPY(3,RLMORI,1,ORIGIN,1) 3273 CALL TIMER('EXPIKR',TIMSTR,TIMEND) 3274 CALL FLSHFO(LUPRI) 3275 END IF 3276C 3277C **** First order magnetic derivative of electric field **** 3278C 3279 IF (CM1) THEN 3280 CALL TIMER('START ',TIMSTR,TIMEND) 3281 IF (FIELD1.EQ.'XYZ-ALL') THEN 3282 FIELD1 = 'X-FIELD' 3283 CALL PR1INT('CM1 ',WORK,LWORK,IDUMMY, 3284 & IDUMMY,TRIANG,PROPRI,IPRONE) 3285 FIELD1 = 'Y-FIELD' 3286 CALL PR1INT('CM1 ',WORK,LWORK,IDUMMY, 3287 & IDUMMY,TRIANG,PROPRI,IPRONE) 3288 FIELD1 = 'Z-FIELD' 3289 CALL PR1INT('CM1 ',WORK,LWORK,IDUMMY, 3290 & IDUMMY,TRIANG,PROPRI,IPRONE) 3291 ELSE 3292 CALL PR1INT('CM1 ',WORK,LWORK,IDUMMY, 3293 & IDUMMY,TRIANG,PROPRI,IPRONE) 3294 END IF 3295 CALL TIMER('CM1 ',TIMSTR,TIMEND) 3296 END IF 3297C 3298C **** 1.order B derivative of electric field at a point **** 3299C 3300 IF (EFBDER) THEN 3301 CALL TIMER('START',TIMSTR,TIMEND) 3302 CALL PR1INT('EFIELB1',WORK,LWORK,IDUMMY, 3303 & IDUMMY,TRIANG,PROPRI,IPRONE) 3304 CALL TIMER('EFBDER',TIMSTR,TIMEND) 3305 END IF 3306C 3307C **** 2.order B derivative of electric field at a point **** 3308C 3309 IF (EFB2DR) THEN 3310 CALL TIMER('START',TIMSTR,TIMEND) 3311 CALL PR1INT('EFIELB2',WORK,LWORK,IDUMMY, 3312 & IDUMMY,TRIANG,PROPRI,IPRONE) 3313 CALL TIMER('EFBDR2',TIMSTR,TIMEND) 3314 END IF 3315C 3316C **** Second order magnetic derivative of electric field **** 3317C 3318 IF (CM2) THEN 3319 CALL TIMER('START ',TIMSTR,TIMEND) 3320C Modified by Bin Gao, December 14, 2009 3321C copied from linsca 3322 IF (FIELD2 .EQ. 'XYZ-ALL') THEN 3323 FIELD2 = 'X-FIELD' 3324 CALL PR1INT('CM2 ',WORK,LWORK,IDUMMY, 3325 & IDUMMY,TRIANG,PROPRI,IPRONE) 3326 FIELD2 = 'Y-FIELD' 3327 CALL PR1INT('CM2 ',WORK,LWORK,IDUMMY, 3328 & IDUMMY,TRIANG,PROPRI,IPRONE) 3329 FIELD2 = 'Z-FIELD' 3330 CALL PR1INT('CM2 ',WORK,LWORK,IDUMMY, 3331 & IDUMMY,TRIANG,PROPRI,IPRONE) 3332 ELSE 3333 CALL PR1INT('CM2 ',WORK,LWORK,IDUMMY, 3334 & IDUMMY,TRIANG,PROPRI,IPRONE) 3335 END IF 3336 CALL TIMER('CM2 ',TIMSTR,TIMEND) 3337 END IF 3338C 3339C **** Second order magnetic derivative of electric field **** 3340C 3341 IF (QDBINT) THEN 3342 CALL TIMER('START ',TIMSTR,TIMEND) 3343Cajt qdb CALL PR1INT('QDBINT ',WORK,LWORK,IDUMMY, 3344Cajt qdb & IDUMMY,TRIANG,PROPRI,IPRONE) 3345 IF (FIELD3 .EQ. 'XYZ-ALL') THEN 3346 FIELD3 = 'XX-FGRD' 3347 CALL PR1INT('QDBINT ',WORK,LWORK,IDUMMY, 3348 & IDUMMY,TRIANG,PROPRI,IPRONE) 3349 FIELD3 = 'XY-FGRD' 3350 CALL PR1INT('QDBINT ',WORK,LWORK,IDUMMY, 3351 & IDUMMY,TRIANG,PROPRI,IPRONE) 3352 FIELD3 = 'YY-FGRD' 3353 CALL PR1INT('QDBINT ',WORK,LWORK,IDUMMY, 3354 & IDUMMY,TRIANG,PROPRI,IPRONE) 3355 FIELD3 = 'XZ-FGRD' 3356 CALL PR1INT('QDBINT ',WORK,LWORK,IDUMMY, 3357 & IDUMMY,TRIANG,PROPRI,IPRONE) 3358 FIELD3 = 'YZ-FGRD' 3359 CALL PR1INT('QDBINT ',WORK,LWORK,IDUMMY, 3360 & IDUMMY,TRIANG,PROPRI,IPRONE) 3361 FIELD3 = 'ZZ-FGRD' 3362 CALL PR1INT('QDBINT ',WORK,LWORK,IDUMMY, 3363 & IDUMMY,TRIANG,PROPRI,IPRONE) 3364 ELSE 3365 CALL PR1INT('QDBINT ',WORK,LWORK,IDUMMY, 3366 & IDUMMY,TRIANG,PROPRI,IPRONE) 3367 END IF 3368Cajt qdb 3369 CALL TIMER('QDBINT',TIMSTR,TIMEND) 3370 END IF 3371C 3372C ***** Test of London contribution to angular momentum ***** 3373C 3374 IF (QDBTST) THEN 3375 LABELT(1) = FIELD3(1:2)//'-QDB X' 3376 LABELT(2) = FIELD3(1:2)//'-QDB Y' 3377 LABELT(3) = FIELD3(1:2)//'-QDB Z' 3378 LABELS(1) = FIELD3(1:2)//'THETA ' 3379 CALL MG1TST(WORK,LWORK,IPRONE,LABELS,LABELT,ORIGIN, 3380 & ' ') 3381 CALL FLSHFO(LUPRI) 3382 END IF 3383C 3384C **** First geometrical derivative of dipole integrals **** 3385C 3386 IF (DPLGRA) THEN 3387 CALL TIMER('START ',TIMSTR,TIMEND) 3388 CALL PR1INT('DPLGRA ',WORK,LWORK,IDUMMY, 3389 & IDUMMY,TRIANG,PROPRI,IPRONE) 3390 CALL TIMER('DPLGRA',TIMSTR,TIMEND) 3391 END IF 3392C 3393C **** Second geometrical derivative of dipole integrals **** 3394C 3395 IF (DIPANH) THEN 3396 CALL TIMER('START ',TIMSTR,TIMEND) 3397 CALL PR1INT('DIPANH ',WORK,LWORK,IDUMMY, 3398 & IDUMMY,TRIANG,PROPRI,IPRONE) 3399 CALL TIMER('DIPANH',TIMSTR,TIMEND) 3400 END IF 3401C 3402C **** First geometrical derivative of quadrupole integrals **** 3403C 3404 IF (QUAGRA) THEN 3405 CALL TIMER('START ',TIMSTR,TIMEND) 3406 CALL PR1INT('QUAGRA ',WORK,LWORK,IDUMMY, 3407 & IDUMMY,TRIANG,PROPRI,IPRONE) 3408 CALL TIMER('QUAGRA',TIMSTR,TIMEND) 3409 END IF 3410C 3411C **** First geometrical derivative of octupole integrals **** 3412C 3413 IF (OCTGRA) THEN 3414 CALL TIMER('START ',TIMSTR,TIMEND) 3415 CALL PR1INT('OCTGRA ',WORK,LWORK,IDUMMY, 3416 & IDUMMY,TRIANG,PROPRI,IPRONE) 3417 CALL TIMER('OCTGRA ',TIMSTR,TIMEND) 3418 END IF 3419C 3420C **** Magnetic quadrupole integrals **** 3421C 3422 IF (MAGQDP) THEN 3423 CALL TIMER('START ',TIMSTR,TIMEND) 3424 CALL PR1INT('MAGQDP ',WORK,LWORK,IDUMMY, 3425 & IDUMMY,TRIANG,PROPRI,IPRONE) 3426 CALL TIMER('MAGQDP',TIMSTR,TIMEND) 3427 CALL FLSHFO(LUPRI) 3428 END IF 3429C 3430C ***** Test of magnetic quadrupole integrals ***** 3431C 3432 IF (MQDPTS) THEN 3433 LABELT(1) = 'XXMAGQDP' 3434 LABELT(2) = 'YXMAGQDP' 3435 LABELT(3) = 'ZXMAGQDP' 3436 LABELS(1) = 'XANGMOM ' 3437 CALL MG1TST(WORK,LWORK,IPRONE,LABELS,LABELT,ORIGIN, 3438 & 'MQDP') 3439 LABELT(1) = 'XYMAGQDP' 3440 LABELT(2) = 'YYMAGQDP' 3441 LABELT(3) = 'ZYMAGQDP' 3442 LABELS(1) = 'YANGMOM ' 3443 CALL MG1TST(WORK,LWORK,IPRONE,LABELS,LABELT,ORIGIN, 3444 & 'MQDP') 3445 LABELT(1) = 'XZMAGQDP' 3446 LABELT(2) = 'YZMAGQDP' 3447 LABELT(3) = 'ZZMAGQDP' 3448 LABELS(1) = 'ZANGMOM ' 3449 CALL MG1TST(WORK,LWORK,IPRONE,LABELS,LABELT,ORIGIN, 3450 & 'MQDP') 3451 CALL FLSHFO(LUPRI) 3452 END IF 3453C 3454C **** Rotational strength integrals **** 3455C 3456 IF (ROTSTR) THEN 3457 CALL TIMER('START ',TIMSTR,TIMEND) 3458 CALL PR1INT('ROTSTR ',WORK,LWORK,IDUMMY, 3459 & IDUMMY,TRIANG,PROPRI,IPRONE) 3460 CALL TIMER('ROTSTR ',TIMSTR,TIMEND) 3461 END IF 3462C 3463C **** Magnetic-field corrected spin-orbit **** 3464C 3465 IF (SOFLD) THEN 3466 CALL TIMER('START ',TIMSTR,TIMEND) 3467 CALL PR1INT('SOFIELD',WORK,LWORK,IDUMMY, 3468 & IDUMMY,TRIANG,PROPRI,IPRONE) 3469 CALL TIMER('SOFIELD',TIMSTR,TIMEND) 3470 END IF 3471C 3472C **** ??? integrals **** 3473C 3474 IF (SOMM) THEN 3475 CALL TIMER('START ',TIMSTR,TIMEND) 3476 CALL PR1INT('SOMAGMO',WORK,LWORK,IDUMMY, 3477 & NPQUAD,TRIANG,PROPRI,IPRONE) 3478 CALL TIMER('SOMAGMO',TIMSTR,TIMEND) 3479 END IF 3480C 3481C **** First electric derivative of overlap integrals **** 3482C **** Type A **** 3483 IF (S1ELE) THEN 3484 CALL PR1INT('S1ELE ',WORK,LWORK,IDUMMY, 3485 & IDUMMY,TRIANG,PROPRI,IPRONE) 3486 END IF 3487C 3488C ****Mean-field spin-orbit integrals **** 3489C 3490 IF (MNF_SO) THEN 3491 CALL AMFIINTEGRALS(PROPRI,IPRONE,WORK,LWORK) 3492 END IF 3493C 3494C **** First electric derivative of overlap integrals **** 3495C **** Type B **** 3496 IF (S1ELB) THEN 3497 CALL PR1INT('S1ELB ',WORK,LWORK,IDUMMY, 3498 & IDUMMY,TRIANG,PROPRI,IPRONE) 3499 END IF 3500C 3501C ****First electric deriv. of one-electron Ham. integrals**** 3502C 3503 IF (ONEELD) THEN 3504 CALL PR1INT('ONEELD ',WORK,LWORK,IDUMMY, 3505 & IDUMMY,TRIANG,PROPRI,IPRONE) 3506 END IF 3507C 3508C **** First geometrical derivatives of integrals *** 3509C 3510 IF (DEROVL) THEN 3511 CALL TIMER('START ',TIMSTR,TIMEND) 3512 CALL PR1INT('DEROVLP',WORK,LWORK,IDUMMY,IDUMMY, 3513 & TRIANG,PROPRI,IPRONE) 3514 CALL TIMER('DEROVL',TIMSTR,TIMEND) 3515 END IF 3516C 3517C **** First geometrical derivatives of integrals *** 3518C 3519 IF (DERAM) THEN 3520 CALL TIMER('START ',TIMSTR,TIMEND) 3521 CALL PR1INT('DERANGM',WORK,LWORK,IDUMMY,IDUMMY, 3522 & TRIANG,PROPRI,IPRONE) 3523 CALL TIMER('DERAM ',TIMSTR,TIMEND) 3524 END IF 3525C 3526C **** First geometr. deriv. of 1el hamiltonian integrals *** 3527C 3528 IF (DERHAM) THEN 3529 CALL TIMER('START ',TIMSTR,TIMEND) 3530 CALL PR1INT('DERHAMI',WORK,LWORK,IDUMMY,IDUMMY, 3531 & TRIANG,PROPRI,IPRONE) 3532 CALL TIMER('DERHAM',TIMSTR,TIMEND) 3533 END IF 3534C 3535C ***** DPT overlap integrals ******************************** 3536C 3537 IF (DPTOVL) THEN 3538 CALL TIMER('START ',TIMSTR,TIMEND) 3539 CALL PR1INT('DPTOVL ',WORK,LWORK,IDUMMY, 3540 & IDUMMY,TRIANG,PROPRI,IPRONE) 3541 CALL TIMER('DPTOVL',TIMSTR,TIMEND) 3542 CALL FLSHFO(LUPRI) 3543 END IF 3544C 3545C ***** Parity violating electroweak interaction ***** 3546C 3547 IF (PVIOLA) THEN 3548 CALL TIMER('START ',TIMSTR,TIMEND) 3549 CALL PR1INT('PVIOLA ',WORK,LWORK,IDUMMY, 3550 & IDUMMY,TRIANG,PROPRI,IPRONE) 3551 CALL TIMER('PVIOLA',TIMSTR,TIMEND) 3552 CALL FLSHFO(LUPRI) 3553 END IF 3554C 3555C ***** DPT pso-lookalikes SQUARED ***** 3556C 3557 IF (XDDXR3) THEN 3558 CALL TIMER('START ',TIMSTR,TIMEND) 3559 CALL PR1INT('XDDXR3 ',WORK,LWORK,IDUMMY, 3560 & IDUMMY,TRIANG,PROPRI,IPRONE) 3561 CALL TIMER('XDDXR3',TIMSTR,TIMEND) 3562 CALL FLSHFO(LUPRI) 3563 END IF 3564cLig call to PR1INT form RANGMO and RPSO options 3565C 3566C ***** (r-r')l' integrals, for calculation of 3567C magetizability ***** 3568C 3569 IF (RANGMO) THEN 3570 CALL TIMER('START ',TIMSTR,TIMEND) 3571 CALL PR1INT('RANGMO ',WORK,LWORK,IDUMMY, 3572 & IDUMMY,TRIANG,PROPRI,IPRONE) 3573 CALL TIMER('RANGMO ',TIMSTR,TIMEND) 3574 CALL FLSHFO(LUPRI) 3575 END IF 3576C 3577C ***** (r-r')l'/|r-R_I|**3 integrals, for shieldings ***** 3578C 3579 IF (RPSO) THEN 3580 CALL TIMER('START ',TIMSTR,TIMEND) 3581 CALL PR1INT('RPSO ',WORK,LWORK,IDUMMY, 3582 & IDUMMY,TRIANG,PROPRI,IPRONE) 3583 CALL TIMER('RPSO ',TIMSTR,TIMEND) 3584 CALL FLSHFO(LUPRI) 3585 END IF 3586C 3587C ***** pXp integrals, for DPT correction to dipole ***** 3588C 3589 IF (PXPINT) THEN 3590 CALL TIMER('START ',TIMSTR,TIMEND) 3591 CALL PR1INT('PXPINT ',WORK,LWORK,IDUMMY, 3592 & IDUMMY,TRIANG,PROPRI,IPRONE) 3593 CALL TIMER('PXPINT',TIMSTR,TIMEND) 3594 CALL FLSHFO(LUPRI) 3595 CALL TIMER('START ',TIMSTR,TIMEND) 3596 CALL PR1INT('PRPINT ',WORK,LWORK,IDUMMY, 3597 & IDUMMY,TRIANG,PROPRI,IPRONE) 3598 CALL TIMER('PRPINT',TIMSTR,TIMEND) 3599 CALL FLSHFO(LUPRI) 3600 END IF 3601cLig 3602 END IF ! (ONEPRP) 3603 END IF ! (RUNONE .OR. ONEPRP) 3604 RUNONE = .TRUE. ! enable one-electron integrals for next iteration if e.g. geometry optimization 3605C 3606C Determine number of different integrals in R12 method (WK/UniKA/19-11-2002). 3607 NOPP12 = 0 3608 IF (R12EIN) THEN 3609 NOPP12 = NOPP12 + 1 3610 ELSE 3611 IF (V12INT) NOPP12 = NOPP12 + 1 3612 IF (R12INT) NOPP12 = NOPP12 + 1 3613 IF (U12INT) NOPP12 = NOPP12 + 1 3614 IF (U21INT) NOPP12 = NOPP12 + 1 3615 END IF 3616 IF (NOPP12 .EQ. 0) 3617 &CALL QUIT('Wrong number of integral types in HERINT') 3618C 3619C Determine pointers to various integrals (WK/UniKA/26-11-2002). 3620 IADV12 = 0 3621 IADR12 = 0 3622 IADU12 = 0 3623 IADU21 = 0 3624 IF (R12EIN) THEN 3625 IADV12 = 1 3626 ELSE 3627 IF (V12INT) IADV12 = 1 3628 IF (R12INT) IADR12 = 1 3629 IF (U12INT) IADU12 = 1 3630 IF (U21INT) IADU21 = 1 3631 END IF 3632 IADR12 = IADV12 + IADR12 3633 IADU12 = IADR12 + IADU12 3634 IADU21 = IADU12 + IADU21 3635C 3636C ********************************** 3637C ***** Two-Electron Integrals ***** 3638C ********************************** 3639C 3640C Cauchy-Scwartz integrals for screening 3641C 3642 IF(DIRCAL) THEN 3643 CALL TIMER('START ',TIMSTR,TIMEND) 3644 CALL GABGEN(0,0,0,0,0,.FALSE.,IPRDEF,WORK,LWORK) 3645C CALL GABGEN(IJOB,ITYPE,IGTYP,MAXDIF,JATOM,NOCONT,WORK,LWORK) 3646 CALL TIMER('GABGEN',TIMSTR,TIMEND) 3647 ENDIF 3648C 3649C Integrals on LUINTA 3650C 3651 IF (RUNTWO) THEN 3652 I2TYP = 0 3653C 3654C ******************************************** 3655C ***** Two-electron repulsion integrals ***** 3656C ******************************************** 3657C 3658 IF (HAMILT) THEN 3659 IF (RUNERI) THEN 3660 IF (IPRDEF .GT. 1) CALL TITLER('Output from ERI 2INT', 3661 & '*',126) 3662 IF (CHIVAL .gt. 0.0D0) THEN 3663 CALL QUIT 3664 & ('ERROR: range-separation is not implemented in ERI') 3665 END IF 3666 CALL TIMER('START ',TIMSTR,TIMEND) 3667 CALL MEMGET('REAL',KCCFBT,MXCONT*MXPRIM,WORK,KFREE,LFREE) 3668 CALL MEMGET('INTE',KINDXB,8*MXSHEL*MXCONT,WORK,KFREE, 3669 & LFREE) 3670 CALL ER2INT(WORK(KCCFBT),WORK(KINDXB),WORK(KFREE),LFREE) 3671 CALL MEMREL('HERINT.ER2INT',WORK,KWORK,KWORK,KFREE,LFREE) 3672 CALL TIMER('ERI 2INT',TIMSTR,TIMEND) 3673 CALL FLSHFO(LUPRI) 3674 ELSE 3675 IF (IPRDEF .GT. 1) CALL TITLER('Output from TWOINT', 3676 & '*',126) 3677 NPAO = MXSHEL*MXAOVC 3678 CALL MEMGET('INTE',KJSTRS,NPAO*2,WORK,KFREE,LFREE) 3679 CALL MEMGET('INTE',KNPRIM,NPAO*2,WORK,KFREE,LFREE) 3680 CALL MEMGET('INTE',KNCONT,NPAO*2,WORK,KFREE,LFREE) 3681 CALL MEMGET('INTE',KIORBS,NPAO ,WORK,KFREE,LFREE) 3682 CALL MEMGET('INTE',KJORBS,NPAO ,WORK,KFREE,LFREE) 3683 CALL MEMGET('INTE',KKORBS,NPAO ,WORK,KFREE,LFREE) 3684 CALL PAOVEC(WORK(KJSTRS),WORK(KNPRIM),WORK(KNCONT), 3685 & WORK(KIORBS),WORK(KJORBS),WORK(KKORBS),0, 3686 & .FALSE.,IPRDEF) 3687 CALL MEMREL('HERINT.PAOVEC',WORK,KWORK,KJORBS,KFREE, 3688 & LFREE) 3689 CALL TIMER('START ',TIMSTR,TIMEND) 3690 IRNTYP = 0 3691 CALL TWOINT(WORK(KFREE),LFREE,DUMMY,DUMMY,DUMMY,1,IDUMMY, 3692 & IDUMMY,DUMMY,IDUMMY,NUMDIS,1,IRNTYP,0,0, 3693 & .TRUE.,.TRUE.,.FALSE.,TKTIME, 3694 & IPRTWO,IPRNTA,IPRNTB,IPRNTC,IPRNTD, 3695 & RTNTWO,IDUMMY,I2TYP,WORK(KJSTRS), 3696 & WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 3697 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 3698 & .FALSE.,.false.) 3699 CALL TIMER('TWOINT',TIMSTR,TIMEND) 3700 CALL FLSHFO(LUPRI) 3701 CALL MEMREL('HERINT.TWOINT',WORK,KWORK,KWORK,KFREE,LFREE) 3702 IF (CHI1ST) CHIVAL = -1.0D0 3703 END IF 3704 END IF 3705C 3706C ****************************************************** 3707C ***** Short-range two-electron integrals (MCDFT) ***** 3708C ****************************************************** 3709C 3710 IF(DOSRIN) THEN 3711 IF (PANAS .EQ. D0 .AND. CHIVAL .EQ. -1.d0) 3712 & CALL QUIT('FATAL ERROR : Short-range integrals '// 3713 & 'specified with the unmodified 2e-repulsion!') 3714 SRINTS = .TRUE. 3715 NPAO = MXSHEL*MXAOVC 3716 CALL MEMGET('INTE',KJSTRS,NPAO*2,WORK,KFREE,LFREE) 3717 CALL MEMGET('INTE',KNPRIM,NPAO*2,WORK,KFREE,LFREE) 3718 CALL MEMGET('INTE',KNCONT,NPAO*2,WORK,KFREE,LFREE) 3719 CALL MEMGET('INTE',KIORBS,NPAO ,WORK,KFREE,LFREE) 3720 CALL MEMGET('INTE',KJORBS,NPAO ,WORK,KFREE,LFREE) 3721 CALL MEMGET('INTE',KKORBS,NPAO ,WORK,KFREE,LFREE) 3722 CALL PAOVEC(WORK(KJSTRS),WORK(KNPRIM),WORK(KNCONT), 3723 & WORK(KIORBS),WORK(KJORBS),WORK(KKORBS),0, 3724 & .FALSE.,IPRDEF) 3725 CALL MEMREL('HERINT.PAOVEC',WORK,KWORK,KJORBS,KFREE, 3726 & LFREE) 3727 CALL TIMER('START ',TIMSTR,TIMEND) 3728 IRNTYP = 0 3729 CALL TWOINT(WORK(KFREE),LFREE,DUMMY,DUMMY,DUMMY,1,IDUMMY, 3730 & IDUMMY,DUMMY,IDUMMY,NUMDIS,1,IRNTYP,0,0, 3731 & .TRUE.,.TRUE.,.FALSE.,TKTIME, 3732 & IPRTWO,IPRNTA,IPRNTB,IPRNTC,IPRNTD, 3733 & RTNTWO,IDUMMY,I2TYP,WORK(KJSTRS), 3734 & WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 3735 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 3736 & .FALSE.,.false.) 3737 CALL TIMER('TWOINT short-range',TIMSTR,TIMEND) 3738 CALL FLSHFO(LUPRI) 3739 CALL MEMREL('HERINT.TWOINT',WORK,KWORK,KWORK,KFREE,LFREE) 3740 SRINTS = .FALSE. 3741 ENDIF 3742C 3743C ***************************************************** 3744C ***** Spatial two-electron spin-orbit integrals ***** 3745C ***************************************************** 3746C 3747 IF (SPNORB.and..not.no2so) THEN 3748 NPAO = MXSHEL*MXAOVC 3749 CALL MEMGET('INTE',KJSTRS,NPAO*2,WORK,KFREE,LFREE) 3750 CALL MEMGET('INTE',KNPRIM,NPAO*2,WORK,KFREE,LFREE) 3751 CALL MEMGET('INTE',KNCONT,NPAO*2,WORK,KFREE,LFREE) 3752 CALL MEMGET('INTE',KIORBS,NPAO ,WORK,KFREE,LFREE) 3753 CALL MEMGET('INTE',KJORBS,NPAO ,WORK,KFREE,LFREE) 3754 CALL MEMGET('INTE',KKORBS,NPAO ,WORK,KFREE,LFREE) 3755 CALL PAOVEC(WORK(KJSTRS),WORK(KNPRIM),WORK(KNCONT), 3756 & WORK(KIORBS),WORK(KJORBS),WORK(KKORBS),0, 3757 & .FALSE.,IPRDEF) 3758 CALL MEMREL('HERINT.PAOVEC',WORK,KWORK,KJORBS,KFREE,LFREE) 3759 CALL TIMER('START ',TIMSTR,TIMEND) 3760 CALL GETTIM(SO1CPU,SO1WAL) 3761 IRNTYP = - 2 3762 CALL TWOINT(WORK(KFREE),LFREE,DUMMY,DUMMY,DUMMY,1,IDUMMY, 3763 & IDUMMY, 3764 & DUMMY,IDUMMY,NUMDIS,1,IRNTYP,0,0,.TRUE.,.TRUE., 3765 & .FALSE.,TKTIME,IPRTWO,IPRNTA,IPRNTB,IPRNTC, 3766 & IPRNTD,RTNTWO,IDUMMY,I2TYP,WORK(KJSTRS), 3767 & WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 3768 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 3769 & .FALSE.,.false.) 3770 CALL MEMREL('HERINT.TWOINT',WORK,KWORK,KWORK,KFREE,LFREE) 3771 CALL GETTIM(SO2CPU,SO2WAL) 3772 SOCPU=SO2CPU-SO1CPU 3773 SOWAL=SO2WAL-SO1WAL 3774 CALL TIMER('TWOINT 2-el spin-orbit',TIMSTR,TIMEND) 3775 WRITE(LUPRI,'(2(/A),2(/A,F7.0,A))') 3776 & ' Two-electron spin-orbit integrals', 3777 & ' =================================', 3778 & ' Spin-orbit 2-electron CPU time ',SOCPU,' seconds', 3779 & ' Spin-orbit 2-electron wall time ',SOWAL,' seconds' 3780 CALL FLSHFO(LUPRI) 3781 END IF 3782C 3783C ***************************************** 3784C ***** Two-electron Darwin integrals ***** 3785C ***************************************** 3786C 3787 IF (DAR2EL) THEN 3788 NPAO = MXSHEL*MXAOVC 3789 CALL MEMGET('INTE',KJSTRS,NPAO*2,WORK,KFREE,LFREE) 3790 CALL MEMGET('INTE',KNPRIM,NPAO*2,WORK,KFREE,LFREE) 3791 CALL MEMGET('INTE',KNCONT,NPAO*2,WORK,KFREE,LFREE) 3792 CALL MEMGET('INTE',KIORBS,NPAO ,WORK,KFREE,LFREE) 3793 CALL MEMGET('INTE',KJORBS,NPAO ,WORK,KFREE,LFREE) 3794 CALL MEMGET('INTE',KKORBS,NPAO ,WORK,KFREE,LFREE) 3795 CALL PAOVEC(WORK(KJSTRS),WORK(KNPRIM),WORK(KNCONT), 3796 & WORK(KIORBS),WORK(KJORBS),WORK(KKORBS),0, 3797 & .FALSE.,IPRDEF) 3798 CALL MEMREL('HERINT.PAOVEC',WORK,KWORK,KJORBS,KFREE,LFREE) 3799 CALL TIMER('START ',TIMSTR,TIMEND) 3800 CALL GETTIM(DA1CPU,DA1WAL) 3801 IRNTYP = 0 3802 DO2DAR = .TRUE. 3803 IF (LUINTA .GT. 0) CALL GPCLOSE(LUINTA,'KEEP') 3804 CALL GPOPEN(LUINTA,'AO2ELDAR','UNKNOWN',' ',' ', 3805 & IDUMMY,.FALSE.) 3806 CALL TWOINT(WORK(KFREE),LFREE,DUMMY,DUMMY,DUMMY,1,IDUMMY, 3807 & IDUMMY, 3808 & DUMMY,IDUMMY,NUMDIS,1,IRNTYP,0,0,.TRUE.,.TRUE., 3809 & .FALSE.,TKTIME,IPRTWO,IPRNTA,IPRNTB,IPRNTC, 3810 & IPRNTD,RTNTWO,IDUMMY,I2TYP,WORK(KJSTRS), 3811 & WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 3812 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 3813 & .FALSE.,.false.) 3814 DO2DAR = .FALSE. 3815 IF (LUINTA .GT. 0) CALL GPCLOSE(LUINTA,'KEEP') 3816 CALL MEMREL('HERINT.TWOINT',WORK,KWORK,KWORK,KFREE,LFREE) 3817 CALL GETTIM(DA2CPU,DA2WAL) 3818 DACPU=DA2CPU-DA1CPU 3819 DAWAL=DA2WAL-DA1WAL 3820 CALL TIMER('TWOINT 2-el Darwin',TIMSTR,TIMEND) 3821 WRITE(LUPRI,'(2(/A),2(/A,F5.0,A))') 3822 & ' Two-electron Darwin integrals', 3823 & ' =================================', 3824 & ' Spin-orbit Darwin CPU time ',DACPU,' seconds', 3825 & ' Spin-orbit Darwin wall time ',DAWAL,' seconds' 3826 CALL FLSHFO(LUPRI) 3827 END IF 3828C 3829C ********************************************** 3830C ***** Two-electron orbit-orbit integrals ***** 3831C ********************************************** 3832C 3833 IF (ORBORB) THEN 3834 NPAO = MXSHEL*MXAOVC 3835 CALL MEMGET('INTE',KJSTRS,NPAO*2,WORK,KFREE,LFREE) 3836 CALL MEMGET('INTE',KNPRIM,NPAO*2,WORK,KFREE,LFREE) 3837 CALL MEMGET('INTE',KNCONT,NPAO*2,WORK,KFREE,LFREE) 3838 CALL MEMGET('INTE',KIORBS,NPAO ,WORK,KFREE,LFREE) 3839 CALL MEMGET('INTE',KJORBS,NPAO ,WORK,KFREE,LFREE) 3840 CALL MEMGET('INTE',KKORBS,NPAO ,WORK,KFREE,LFREE) 3841 CALL PAOVEC(WORK(KJSTRS),WORK(KNPRIM),WORK(KNCONT), 3842 & WORK(KIORBS),WORK(KJORBS),WORK(KKORBS),0, 3843 & .FALSE.,IPRDEF) 3844 CALL MEMREL('HERINT.PAOVEC',WORK,KWORK,KJORBS,KFREE,LFREE) 3845 CALL TIMER('START ',TIMSTR,TIMEND) 3846 CALL GETTIM(DA1CPU,DA1WAL) 3847 IRNTYP = 0 3848 BPH2OO = .TRUE. 3849 IF (LUINTA .GT. 0) CALL GPCLOSE(LUINTA,'KEEP') 3850 CALL GPOPEN(LUINTA,'2-ORBORB','UNKNOWN',' ',' ', 3851 & IDUMMY,.FALSE.) 3852 CALL TWOINT(WORK(KFREE),LFREE,DUMMY,DUMMY,DUMMY,1,IDUMMY, 3853 & IDUMMY, 3854 & DUMMY,IDUMMY,NUMDIS,1,IRNTYP,0,0,.TRUE.,.TRUE., 3855 & .FALSE.,TKTIME,IPRTWO,IPRNTA,IPRNTB,IPRNTC, 3856 & IPRNTD,RTNTWO,IDUMMY,I2TYP,WORK(KJSTRS), 3857 & WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 3858 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 3859 & .FALSE.,.false.) 3860 BPH2OO = .FALSE. 3861 IF (LUINTA .GT. 0) CALL GPCLOSE(LUINTA,'KEEP') 3862 CALL MEMREL('HERINT.TWOINT',WORK,KWORK,KWORK,KFREE,LFREE) 3863 CALL GETTIM(DA2CPU,DA2WAL) 3864 DACPU=DA2CPU-DA1CPU 3865 DAWAL=DA2WAL-DA1WAL 3866 CALL TIMER('TWOINT',TIMSTR,TIMEND) 3867 WRITE(LUPRI,'(2(/A),2(/A,F5.0,A))') 3868 & ' Two-electron Darwin integrals', 3869 & ' =================================', 3870 & ' Orbit-orbit CPU time ',DACPU,' seconds', 3871 & ' Orbit-orbit wall time ',DAWAL,' seconds' 3872 CALL FLSHFO(LUPRI) 3873 END IF 3874C 3875C ************************************* 3876C ***** Test spin-orbit integrals ***** 3877C ************************************* 3878C 3879 IF (SOTEST) THEN 3880 CALL SOCHK2(WORK,LWORK,IPRDEF) 3881 CALL FLSHFO(LUPRI) 3882 END IF 3883C 3884C ******************************************************** 3885C ***** Test of two-electron part of magnetic moment ***** 3886C ******************************************************** 3887C 3888 IF (MGMO2T) THEN 3889 CALL MM2TST(WORK,LWORK,PRTHRS,IPRINT) 3890 END IF 3891 END IF 3892C 3893 CALL QEXIT('HERINT') 3894 RETURN 3895 END 3896C 3897 SUBROUTINE MAKE_AOTWOINT(WORK,LWORK) 3898C 3899C 11-jan-2007 Hans Joergen Aa. Jensen 3900C 3901C Calculate two-electron integrals 3902C 3903#include "implicit.h" 3904 DIMENSION WORK(LWORK) 3905#include "priunit.h" 3906C 3907C Used from common blocks: 3908C GNRINF: DIRCAL, DOSRIN 3909C CBIHER: ONEPRP 3910C CBIHR1: RUNONE 3911C CBIHR2: RUNTWO 3912C INFTAP: LUINTA, LUSRINT 3913#include "mxcent.h" 3914#include "gnrinf.h" 3915#include "cbiher.h" 3916#include "cbihr1.h" 3917#include "cbihr2.h" 3918#include "inftap.h" 3919C 3920 LOGICAL AOEXIST, SREXIST 3921 LOGICAL DIRCAL_SAV, RUNONE_SAV, ONEPRP_SAV, RUNTWO_SAV 3922C 3923 CALL QENTER('MAKE_AOTWOINT') 3924 CALL GPINQ('AOTWOINT','EXIST',AOEXIST) 3925 IF (DOSRIN) THEN 3926 CALL GPINQ('AOSR2INT','EXIST',SREXIST) 3927 ELSE 3928 SREXIST = .TRUE. ! to signal not needed 3929 END IF 3930 IF (.NOT. AOEXIST .OR. .NOT. SREXIST) THEN 3931 RUNONE_SAV = RUNONE 3932 ONEPRP_SAV = ONEPRP 3933 RUNTWO_SAV = RUNTWO 3934 DIRCAL_SAV = DIRCAL 3935 RUNONE = .FALSE. 3936 ONEPRP = .FALSE. 3937 RUNTWO = .TRUE. 3938 DIRCAL = .FALSE. 3939C Open LUINTA/LUSRINT for new integrals 3940 WRITE (LUPRI,'(/A)') '---> (Re)generating AOTWOINT' 3941 CALL GPOPEN(LUINTA,'AOTWOINT','UNKNOWN',' ',' ',IDUMMY,.FALSE.) 3942 IF (DOSRIN) THEN 3943 WRITE (LUPRI,'(A)') '---> and (re)generating AOSR2INT' 3944 CALL GPOPEN(LUSRINT,'AOSR2INT','UNKNOWN',' ',' ', 3945 & IDUMMY,.FALSE.) 3946 END IF 3947 IF (LUPROP .LT. 0) 3948 & CALL GPOPEN(LUPROP,'AOPROPER','OLD',' ',' ',IDUMMY,.FALSE.) 3949 CALL HERINT(WORK,LWORK) 3950 RUNONE = RUNONE_SAV 3951 ONEPRP = ONEPRP_SAV 3952 RUNTWO = RUNTWO_SAV 3953 DIRCAL = DIRCAL_SAV 3954 END IF 3955 3956 CALL QEXIT('MAKE_AOTWOINT') 3957 RETURN 3958 END 3959#endif /* ifndef PRG_DIRAC */ 3960C /* Deck setdch */ 3961 SUBROUTINE SETDCH 3962C 3963#include "implicit.h" 3964#include "mxcent.h" 3965#include "maxorb.h" 3966#include "maxaqn.h" 3967C 3968#include "symmet.h" 3969#include "dorps.h" 3970C 3971 DO 100 IREP = 0, MAXREP 3972 DOREPS(IREP) = .TRUE. 3973 100 CONTINUE 3974 RETURN 3975 END 3976C /* Deck srtinp */ 3977 SUBROUTINE SRTINP(WORD) 3978#include "implicit.h" 3979#include "priunit.h" 3980C 3981C inftra : USE_INTSORT 3982C 3983#include "cbisor.h" 3984#include "inftra.h" 3985 PARAMETER (NTABLE = 10) 3986C 3987 LOGICAL SET, NEWDEF 3988 CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 3989C 3990 SAVE SET 3991 DATA TABLE /'xxxxxxx', 'xxxxxx ', '.THRQ ', '.PRINT ', '.IO PRI', 3992 & '.INTSYM', '.DELAO ', '.KEEP ', 'XXXXXXX', 'XXXXXXX'/ 3993 DATA SET/.FALSE./ 3994C 3995 IF (SET) THEN 3996 IF (WORD .NE. '*END OF' .AND. WORD(1:2) .NE. '**') THEN 3997 969 READ (LUCMD, '(A7)') WORD 3998 CALL UPCASE(WORD) 3999 PROMPT = WORD(1:1) 4000 IF (PROMPT .NE. '*') GO TO 969 4001 END IF 4002 RETURN 4003 END IF 4004C 4005 SET = .TRUE. 4006C 4007C Initialize /INTSRT/ 4008C 4009 THRQ2 = 1.0D-15 4010 ISPRINT = 0 4011 ISPRFIO = 0 4012 ISNTSYM = 1 4013 DELAO = .FALSE. 4014 DO I=1,8 4015 ISKEEP(I) = 0 4016 END DO 4017C 4018 NEWDEF = WORD .EQ. '*SORINT' 4019 ICHANG = 0 4020 IF (NEWDEF) THEN 4021 WORD1 = WORD 4022 100 CONTINUE 4023 READ (LUCMD, '(A7)') WORD 4024 CALL UPCASE(WORD) 4025 PROMPT = WORD(1:1) 4026 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 4027 GO TO 100 4028 ELSE IF (PROMPT .EQ. '.') THEN 4029 ICHANG = ICHANG + 1 4030 DO 200 I = 1, NTABLE 4031 IF (TABLE(I) .EQ. WORD) THEN 4032 GO TO (1,2,3,4,5,6,7,8,9,10), I 4033 END IF 4034 200 CONTINUE 4035 IF (WORD .EQ. '.OPTION') THEN 4036 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 4037 GO TO 100 4038 END IF 4039 WRITE (LUPRI,'(/,3A,/)') ' Keyword "',WORD, 4040 * '" not recognized in *SORINT.' 4041 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 4042 CALL QUIT('Illegal keyword in *SORINT.') 4043 1 CONTINUE 4044 GO TO 100 4045 2 CONTINUE 4046 GO TO 100 4047 3 CONTINUE 4048 READ (LUCMD, *) THRQ2 4049 IF (THRQ2.NE.1D-15) ICHANG=ICHANG+1 4050 GO TO 100 4051 4 CONTINUE 4052 READ (LUCMD, *) ISPRINT 4053 IF (ISPRINT.NE.0) ICHANG=ICHANG+1 4054 GO TO 100 4055 5 CONTINUE 4056 READ (LUCMD, *) ISPRFIO 4057 IF (ISPRFIO.NE.0) ICHANG=ICHANG+1 4058 GO TO 100 4059 6 CONTINUE 4060 READ (LUCMD, *) ISNTSYM 4061 IF (ISNTSYM .NE. 1) ICHANG = ICHANG + 1 4062 7 CONTINUE 4063 DELAO = .TRUE. 4064 ICHANG=ICHANG+1 4065 GO TO 100 4066 8 CONTINUE 4067 READ (LUCMD,*) (ISKEEP(I), I=1,8) 4068 GO TO 100 4069 9 CONTINUE 4070 GO TO 100 4071 10 CONTINUE 4072 GO TO 100 4073 ELSE IF (PROMPT .EQ. '*') THEN 4074 GO TO 300 4075 ELSE 4076 WRITE (LUPRI,'(/,3A,/)') ' Prompt "',WORD, 4077 * '" not recognized in *SORINT.' 4078 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 4079 CALL QUIT('Illegal prompt in *SORINT.') 4080 END IF 4081 END IF 4082 300 CONTINUE 4083 IF (ICHANG .EQ. 0) RETURN 4084 IF (USE_INTSORT .AND. NEWDEF) THEN 4085 CALL HEADER('Changes of defaults for SORTAO:',1) 4086 IF (THRQ2 .NE. 1.0D-15) WRITE (LUPRI,'(A,D12.2)') 4087 & ' Threshold for sorted integrals:', THRQ2 4088 IF (ISPRINT .NE. 0) WRITE(LUPRI,'(A,I5)') 4089 & ' Print level in SORTAO',ISPRINT 4090 IF (ISPRFIO .NE. 0) WRITE(LUPRI,'(A,I5)') 4091 & ' Print level in fast I/O routines',ISPRFIO 4092 IF (ISNTSYM.NE.1) WRITE(LUPRI,'(A,I5)') 4093 & ' Integral symmetry ',ISNTSYM 4094 IF (DELAO) WRITE(LUPRI,'(A)') 4095 & ' Delete basic file of two-electron integrals' 4096 IF (ISUM(8,ISKEEP,1) .NE.0) WRITE (LUPRI,'(A,8I2)') 4097 & ' Integrals to be kept in symmetries with "0":', 4098 & (ISKEEP(I),I=1,8) 4099 END IF 4100 RETURN 4101 END 4102C/* Deck herini */ 4103 SUBROUTINE HERINI 4104#include "implicit.h" 4105#include "priunit.h" 4106#include "mxcent.h" 4107#include "maxorb.h" 4108#include "nuclei.h" 4109#include "gnrinf.h" 4110#include "inftra.h" 4111#include "cbirea.h" 4112#include "cbisol.h" 4113#include "cbiher.h" 4114#include "drw2el.h" 4115#include "orgcom.h" 4116#include "elweak.h" 4117#include "r12int.h" 4118#include "infinp.h" 4119#include "qm3.h" 4120C 4121Cholesky 4122#include "ccdeco.h" 4123Cholesky 4124 LOGICAL FIRST_CALL 4125 SAVE FIRST_CALL 4126 DATA FIRST_CALL /.TRUE./ 4127C 4128 IF (.NOT. FIRST_CALL) RETURN 4129 FIRST_CALL = .FALSE. 4130C 4131C Initialize /CBIHER/ 4132C 4133 USE_INTSORT = .FALSE. ! in inftra.h 4134 IF(RELCAL) THEN 4135 HAMILT = .FALSE. 4136 ONEPRP = .FALSE. 4137 DIPVEL = .FALSE. 4138 DIRAC_HER = .TRUE. 4139 SUPMAT = .FALSE. 4140 ELSE 4141 HAMILT = .TRUE. 4142 ONEPRP = .FALSE. 4143 DIPVEL = .FALSE. 4144 DIRAC_HER = .FALSE. 4145 SUPMAT = .NOT. DIRCAL .AND. .NOT. USE_INTSORT 4146 ENDIF 4147 SKIP = .FALSE. 4148 TSTINP = .FALSE. 4149 ONEPRP = .FALSE. 4150 HAMILT = .TRUE. 4151 DIPLEN = .FALSE. 4152 QUADRU = .FALSE. 4153 SPNORB = .FALSE. 4154 no2so = .false. 4155 SOTEST = .FALSE. 4156 NOTWO = DIRCAL 4157 TWOBYTEPACKING = .FALSE. 4158 SECMOM = .FALSE. 4159 CARMOM = .FALSE. 4160 SPHMOM = .FALSE. 4161 FERMI = .FALSE. 4162 PSO = .FALSE. 4163 SPIDIP = .FALSE. 4164 DSO = .FALSE. 4165 SDFC = .FALSE. 4166 PROPRI = .FALSE. 4167 HDO = .FALSE. 4168 S1MAG = .FALSE. 4169 S2MAG = .FALSE. 4170 ANGMOM = .FALSE. 4171 ANGLON = .FALSE. 4172 LONMOM = .FALSE. 4173 MAGMOM = .FALSE. 4174 S1MAGT = .FALSE. 4175 MGMOMT = .FALSE. 4176 KINENE = .FALSE. 4177 S2MAGT = .FALSE. 4178 DSUSNL = .FALSE. 4179 DSUSLL = .FALSE. 4180 DSUSLH = .FALSE. 4181 DIASUS = .FALSE. 4182 DSUTST = .FALSE. 4183 NUCSNL = .FALSE. 4184 NUCSLO = .FALSE. 4185 NUCSHI = .FALSE. 4186 NSTTST = .FALSE. 4187 NSLTST = .FALSE. 4188 NELFLD = .FALSE. 4189 NSNLTS = .FALSE. 4190 EFGCAR = .FALSE. 4191 EFGSPH = .FALSE. 4192 S1MAGL = .FALSE. 4193 S1MAGR = .FALSE. 4194 HDOBR = .FALSE. 4195 S1MLT = .FALSE. 4196 S1MRT = .FALSE. 4197 HDOBRT = .FALSE. 4198 NUCPOT = .FALSE. 4199 NPOTST = .FALSE. 4200 MGMO2T = .FALSE. 4201 HBDO = .FALSE. 4202 SUSCGO = .FALSE. 4203 NSTCGO = .FALSE. 4204 MASSVL = .FALSE. 4205 DARWIN = .FALSE. 4206 CM1 = .FALSE. 4207 CM2 = .FALSE. 4208 QDBINT = .FALSE. 4209 SQHDOL = .FALSE. 4210 SQHDOR = .FALSE. 4211 S1ELE = .FALSE. 4212 S1ELB = .FALSE. 4213 ONEELD = .FALSE. 4214 THETA = .FALSE. 4215 DPLGRA = .FALSE. 4216 QUAGRA = .FALSE. 4217 OCTGRA = .FALSE. 4218 ROTSTR = .FALSE. 4219 THRMOM = .FALSE. 4220 SOFLD = .FALSE. 4221 SOMM = .FALSE. 4222 DEROVL = .FALSE. 4223 DERAM = .FALSE. 4224 DERHAM = .FALSE. 4225 ELGDIA = .FALSE. 4226 ELGDIL = .FALSE. 4227 MNF_SO = .FALSE. 4228 PVPINT = .FALSE. 4229 POTENE = .FALSE. 4230 PXPINT = .FALSE. 4231 NOPICH = .FALSE. 4232 PRTHRS = 1.0D-10 4233 NPQUAD = 40 4234 ALLATM = .TRUE. 4235 TRIANG = .TRUE. 4236 EXPIKR = .FALSE. 4237 DO2DAR = .FALSE. 4238 DAR2EL = .FALSE. 4239 AD2DAR = .FALSE. 4240 ORBORB = .FALSE. 4241 FINDPT = .FALSE. 4242 NFIELD = 0 4243 DPTINT = .FALSE. 4244 NO2DPT = .FALSE. 4245 DPTFAC = 0.D0 4246 DARFAC = 0.D0 4247 S4CENT = .FALSE. 4248 DPTOVL = .FALSE. 4249 DPTPOT = .FALSE. 4250 XDDXR3 = .FALSE. 4251 LFDIPLN = .FALSE. 4252C Initialize R12 logical variables R12CAL (WK/UniKA/04-11-2002). 4253 R12DIA = .TRUE. 4254 R12SVD = .FALSE. 4255 NOAUXB = .FALSE. 4256 R12PRP = .FALSE. 4257 NOTONE = .FALSE. 4258 NOTTWO = .FALSE. 4259 NOTTRE = .FALSE. 4260 IANCC2 = 0 4261 IAPCC2 = 0 4262 R12LEV = 0.D0 4263 R12RST = .FALSE. 4264 R12CAL = .FALSE. 4265 R12XXL = .FALSE. 4266 R12CBS = .FALSE. 4267 R12OLD = .FALSE. 4268 R12SQR = .FALSE. 4269 R12TRA = .FALSE. 4270 R12INT = .FALSE. 4271 R12EIN = .FALSE. 4272 R12ECO = .FALSE. 4273 R12EOR = .FALSE. 4274 SLATER = .FALSE. 4275 R12NOA = .FALSE. 4276 R12NOB = .FALSE. 4277 R12NOP = .FALSE. 4278 R12HYB = .TRUE. 4279 NORXR = .FALSE. 4280 MAGQDP = .FALSE. 4281 MQDPTS = .FALSE. 4282 U12INT = .FALSE. 4283 U21INT = U12INT 4284 V12INT = .TRUE. 4285 NOTV12 = .NOT. V12INT 4286 NOTR12 = .NOT. R12INT 4287 NOTU12 = .NOT. U12INT 4288 ANTICO = .FALSE. 4289 DIRSCF = .FALSE. 4290 U12DIR = .FALSE. 4291 R12DIR = .FALSE. 4292 V12DIR = .FALSE. 4293 DCCR12 = .FALSE. 4294 STEPIV = .FALSE. 4295 CUSP12 = .FALSE. 4296 DUMR12 = .FALSE. 4297 ONEAUX = .FALSE. 4298 INTGAC = 0 4299 INTGAD = 0 4300 COMBSS = .FALSE. 4301 LAUXBS = .FALSE. 4302 LOOPDP = .FALSE. 4303 VCLTHR = 0D0 4304 SVDTHR = 1D-12 4305 DO I = 1,8 4306 NRXR12(I) = 0 4307 LOCFRO(I) = 0 4308 END DO 4309cLig setting RANGMO and RPSO false 4310 RANGMO = .FALSE. 4311 RPSO = .FALSE. 4312cLig 4313 DOLRINTS=.FALSE. ! jim-dgb DOLRINT init to FALSE 4314C 4315 CALL DZERO(ORIGIN,3) 4316chj may 09: GAGORG and DIPORG now initialized in READIN 4317chj CALL DZERO(GAGORG,3) 4318chj CALL DZERO(DIPORG,3) 4319 CALL DZERO(EXPKR ,3) 4320C 4321Cholesky 4322 IF (CHOINT) NOTWO = .TRUE. 4323Cholesky 4324C 4325C Initialize /CBISOL/ (10-Dec-92 th+hjaaj) 4326C -- not used in Hermit; SOLVNT must be false 4327C in order to skip solvent modules in ONEDRV 4328C and cavity center in VIBMAS 4329C 4330 SOLVNT = .FALSE. 4331C 4332 PVFINN = 1.D0 4333 BGWEIL = .FALSE. 4334 PVIOLA = .FALSE. 4335C 4336C For .R12AUX the cartesian moments CMxxxxxx up to order two 4337C are needed. 4338 IF (LMULBS) CALL SET_CARMOM_2 4339C 4340 RETURN 4341 END 4342C/* Deck amfiin */ 4343 SUBROUTINE AMFIINTEGRALS(PROPRI,IPRINT,WORK,LWORK) 4344C 4345C Interface to B.Schimmelpfennig's atomic mean-field integral code 4346C Written by K.Ruud in northern Norway in January 2001. 4347C 4348#include "implicit.h" 4349#include "dummy.h" 4350#include "mxcent.h" 4351C 4352 LOGICAL BREIT, PROPRI 4353 DIMENSION WORK(LWORK) 4354#include "gnrinf.h" 4355#include "inftap.h" 4356#include "nuclei.h" 4357#include "cbirea.h" 4358C 4359C Open Mean-field input file generate in READIN 4360C 4361 LUAMFI_INP = -1 4362 LUMNFPRP = -1 4363 CALL GPOPEN(LUAMFI_INP,'AMFI_MNF.INP','OLD',' ','FORMATTED', 4364 & IDUMMY,.FALSE.) 4365 CALL GPOPEN(LUMNFPRP,'AOPROPER.MNF','UNKNOWN',' ','UNFORMATTED', 4366 & IDUMMY,.FALSE.) 4367C 4368C For the time being, no Breit-integrals (but soon....) 4369C 4370 BREIT = .NOT.DKTRAN 4371CBS SO FAR ONLY BREIT-PAULI !!! Not no-pair 4372C 4373C We loop over symmetry-independent centers 4374C 4375 REWIND (LUMNFPRP) 4376 DO IATOM = 1, NUCIND 4377 CALL AMFI(LUAMFI_INP,LUMNFPRP,BREIT,GAUNUC,GNUEXP(IATOM), 4378 & WORK,LWORK) 4379 END DO 4380C 4381C We should now have all the atomic mean-field SO integrals, 4382C now generate them in symmetry-orbital basis 4383C 4384 REWIND (LUMNFPRP) 4385 CALL SYMTRAFO(LUMNFPRP,PROPRI,IPRINT,WORK,LWORK) 4386C 4387 CALL GPCLOSE(LUAMFI_INP,'KEEP') 4388 CALL GPCLOSE(LUMNFPRP,'DELETE') 4389 RETURN 4390 END 4391C/* Deck rstgft */ 4392 subroutine rstgft(alpha,i,ggan,ggas) 4393C 4394C R12 calculation requested with fitted "r12-times-Slater" geminal (WK/UniKA/03-24-2005). 4395C 4396#include "implicit.h" 4397 dimension ggax(6,6),ggac(6,6),ggas(6),ggan(6) 4398 character*5 fn 4399 data ggax / 4400 &0.5230D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4401 &0.3383D+00,0.2652D+01,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4402 &0.2670D+00,0.1503D+01,0.7928D+01,0.0000D+00,0.0000D+00,0.0000D+00, 4403 &0.2272D+00,0.1077D+01,0.4270D+01,0.1885D+02,0.0000D+00,0.0000D+00, 4404 &0.2011D+00,0.8520D+00,0.2941D+01,0.9710D+01,0.3980D+02,0.0000D+00, 4405 &0.1824D+00,0.7118D+00,0.2252D+01,0.6474D+01,0.1966D+02,0.7792D+02 4406 & / 4407 data ggac / 4408 &0.6498D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4409 &0.4806D+00,0.3307D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4410 &0.3880D+00,0.3184D+00,0.1768D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4411 &0.3259D+00,0.3108D+00,0.1755D+00,0.1102D+00,0.0000D+00,0.0000D+00, 4412 &0.2804D+00,0.3026D+00,0.1785D+00,0.1100D+00,0.7450D-01,0.0000D+00, 4413 &0.2454D+00,0.2938D+00,0.1815D+00,0.1128D+00,0.7502D-01,0.5280D-01 4414 & / 4415c-----data ggax /---old fit without weight function exp(-2x)------ 4416c 1 0.1760d0, 0.0000d0, 0.0000d0, 0.0000d0, 0.0000d0, 0.0000d0, 4417c 2 0.1069d0, 0.4323d0, 0.0000d0, 0.0000d0, 0.0000d0, 0.0000d0, 4418c 3 0.0801d0, 0.2359d0, 0.9192d0, 0.0000d0, 0.0000d0, 0.0000d0, 4419c 4 0.0654d0, 0.1644d0, 0.4662d0, 1.7980d0, 0.0000d0, 0.0000d0, 4420c 5 0.0561d0, 0.1273d0, 0.3080d0, 0.8643d0, 3.3200d0, 0.0000d0, 4421c 6 0.0495d0, 0.1047d0, 0.2289d0, 0.5475d0, 1.5300d0, 5.8660d0 4422c & / 4423c data ggac / 4424c 1 0.2811d0, 0.0000d0, 0.0000d0, 0.0000d0, 0.0000d0, 0.0000d0, 4425c 2 0.1036d0, 0.3999d0, 0.0000d0, 0.0000d0, 0.0000d0, 0.0000d0, 4426c 3 0.0454d0, 0.2352d0, 0.3692d0, 0.0000d0, 0.0000d0, 0.0000d0, 4427c 4 0.0220d0, 0.1459d0, 0.2781d0, 0.3006d0, 0.0000d0, 0.0000d0, 4428c 5 0.0114d0, 0.0929d0, 0.2126d0, 0.2601d0, 0.2354d0, 0.0000d0, 4429c 6 0.0062d0, 0.0603d0, 0.1620d0, 0.2259d0, 0.2212d0, 0.1828d0 4430c----& / 4431 if (i .ge. 1 .and. i .le. 6) then 4432 do k=1,i 4433 ggas(k)=ggax(k,i)*alpha**2 4434 enddo 4435 write(fn,'(a4,i1)') 'sin.',i 4436 LUCFN=-1 4437 CALL GPOPEN(LUCFN,FN,'NEW',' ','FORMATTED',IDUMMY,.FALSE.) 4438 write(LUCFN,'(a,i1,a)') 'set output ''rsg',i,'.ps''' 4439 write(LUCFN,'(a)') 'set term postscript' 4440 write(LUCFN,'(a,f25.20,a)') 'rsg(x)=x*exp(-x*',alpha,')' 4441 do k=1,i 4442 ggan(k)=ggac(k,i) 4443 enddo 4444 do k=1,i 4445 if (k.eq.1) then 4446 if (k.lt.i) then 4447 write(LUCFN,'(a,f25.20,a,f25.20,a)') 4448 & 'fit(x)=',ggan(k),'*x*exp(-x**2*',ggas(k),') + \\' 4449 else 4450 write(LUCFN,'(a,f25.20,a,f25.20,a)') 4451 & 'fit(x)=',ggan(k),'*x*exp(-x**2*',ggas(k),')' 4452 endif 4453 else 4454 if (k.lt.i) then 4455 write(LUCFN,'(f25.20,a,f25.20,a)') 4456 & ggan(k),'*x*exp(-x**2*',ggas(k),') + \\' 4457 else 4458 write(LUCFN,'(f25.20,a,f25.20,a)') 4459 & ggan(k),'*x*exp(-x**2*',ggas(k),')' 4460 endif 4461 endif 4462 enddo 4463 write(LUCFN,'(a,i10,a)') 'set xrange [0:',nint(5d0/alpha),']' 4464 write(LUCFN,'(a,F10.1,a)') 4465 & 'set yrange [0:',nint(5d0/alpha)/10.0d0,']' 4466 write(LUCFN,'(a)') 'plot fit(x), rsg(x)' 4467 call gpclose(LUCFN,'KEEP') 4468 else 4469 call quit('This STG fit is not possible. Sorry.') 4470 end if 4471 end 4472C/* Deck stgfit */ 4473 subroutine stgfit(alpha,i,ggan,ggas) 4474C 4475C R12 calculation requested with fitted Slater-type geminal (WK/UniKA/03-24-2005). 4476C 4477#include "implicit.h" 4478 dimension ggax(6,6),ggac(6,6),ggas(6),ggan(6) 4479 character*5 fn 4480 data ggax / 4481 &0.6853D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4482 &0.4254D+00,0.4520D+01,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4483 &0.3303D+00,0.2321D+01,0.1628D+02,0.0000D+00,0.0000D+00,0.0000D+00, 4484 &0.2783D+00,0.1591D+01,0.7637D+01,0.4574D+02,0.0000D+00,0.0000D+00, 4485 &0.2447D+00,0.1225D+01,0.4924D+01,0.1988D+02,0.1127D+03,0.0000D+00, 4486 &0.2209D+00,0.1004D+01,0.3622D+01,0.1216D+02,0.4587D+02,0.2544D+03 4487 &/ 4488 data ggac / 4489 &0.7354D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4490 &0.5640D+00,0.3102D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4491 &0.4683D+00,0.3087D+00,0.1529D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4492 &0.4025D+00,0.3090D+00,0.1570D+00,0.8898D-01,0.0000D+00,0.0000D+00, 4493 &0.3532D+00,0.3072D+00,0.1629D+00,0.9321D-01,0.5619D-01,0.0000D+00, 4494 &0.3144D+00,0.3037D+00,0.1681D+00,0.9811D-01,0.6024D-01,0.3726D-01 4495 &/ 4496 if (i .ge. 1 .and. i .le. 6) then 4497 do k=1,i 4498 ggas(k)=ggax(k,i)*alpha**2 4499 enddo 4500 write(fn,'(a4,i1)') 'gin.',i 4501 LUCFN=-1 4502 CALL GPOPEN(LUCFN,FN,'NEW',' ','FORMATTED',IDUMMY,.FALSE.) 4503 write(LUCFN,'(a,i1,a)') 'set output ''stg',i,'.ps''' 4504 write(LUCFN,'(a)') 'set term postscript' 4505 write(LUCFN,'(a,f25.20,a)') 'stg(x)=exp(-x*',alpha,')' 4506 do k=1,i 4507 ggan(k)=ggac(k,i) 4508 enddo 4509 do k=1,i 4510 if (k.eq.1) then 4511 if (k.lt.i) then 4512 write(LUCFN,'(a,f25.20,a,f25.20,a)') 4513 & 'fit(x)=',ggan(k),'*exp(-x**2*',ggas(k),') + \\' 4514 else 4515 write(LUCFN,'(a,f25.20,a,f25.20,a)') 4516 & 'fit(x)=',ggan(k),'*exp(-x**2*',ggas(k),')' 4517 endif 4518 else 4519 if (k.lt.i) then 4520 write(LUCFN,'(f25.20,a,f25.20,a)') 4521 & ggan(k),'*exp(-x**2*',ggas(k),') + \\' 4522 else 4523 write(LUCFN,'(f25.20,a,f25.20,a)') 4524 & ggan(k),'*exp(-x**2*',ggas(k),')' 4525 endif 4526 endif 4527 enddo 4528 write(LUCFN,'(a,i10,a)') 'set xrange [0:',nint(3d0/alpha),']' 4529 write(LUCFN,'(a,i10,a)') 'set yrange [0:1]' 4530 write(LUCFN,'(a)') 'plot fit(x), stg(x)' 4531 call gpclose(LUCFN,'KEEP') 4532 else 4533 call quit('This STG fit is not possible. Sorry.') 4534 end if 4535 end 4536C/* Deck erffit */ 4537 subroutine erffit(alpha,i,ggan,ggas) 4538C 4539C R12 calculation requested with fitted ERFC-type geminal (WK/UniKA/03-29-2005). 4540C 4541#include "implicit.h" 4542 dimension ggax(6,6),ggac(6,6),ggas(6),ggan(6) 4543 character*5 fn 4544 data ggax / 4545 &0.1701D+01,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4546 &0.1370D+01,0.8312D+01,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4547 &0.1256D+01,0.4702D+01,0.2877D+02,0.0000D+00,0.0000D+00,0.0000D+00, 4548 &0.1198D+01,0.3495D+01,0.1406D+02,0.8056D+02,0.0000D+00,0.0000D+00, 4549 &0.1161D+01,0.2888D+01,0.9417D+01,0.3575D+02,0.1991D+03,0.0000D+00, 4550 &0.1136D+01,0.2521D+01,0.7177D+01,0.2232D+02,0.8204D+02,0.4515D+03 4551 &/ 4552 data ggac / 4553 &0.7658D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4554 &0.6232D+00,0.2676D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4555 &0.5466D+00,0.2623D+00,0.1308D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4556 &0.4948D+00,0.2604D+00,0.1327D+00,0.7585D-01,0.0000D+00,0.0000D+00, 4557 &0.4563D+00,0.2577D+00,0.1363D+00,0.7886D-01,0.4774D-01,0.0000D+00, 4558 &0.4260D+00,0.2543D+00,0.1393D+00,0.8244D-01,0.5094D-01,0.3157D-01 4559 &/ 4560 if (i .ge. 1 .and. i .le. 6) then 4561 do k=1,i 4562 ggas(k)=ggax(k,i)*alpha**2 4563 enddo 4564 write(fn,'(a4,i1)') 'cin.',i 4565 LUCFN=-1 4566 CALL GPOPEN(LUCFN,FN,'NEW',' ','FORMATTED',IDUMMY,.FALSE.) 4567 write(LUCFN,'(a,i1,a)') 'set output ''etg',i,'.ps''' 4568 write(LUCFN,'(a)') 'set term postscript' 4569 write(LUCFN,'(a,f25.20,a)') 'etg(x)=1.0-erf(x*',alpha,')' 4570 do k=1,i 4571 ggan(k)=ggac(k,i) 4572 enddo 4573 do k=1,i 4574 if (k.eq.1) then 4575 if (k.lt.i) then 4576 write(LUCFN,'(a,f25.20,a,f25.20,a)') 4577 & 'fit(x)=',ggan(k),'*exp(-x**2*',ggas(k),') + \\' 4578 else 4579 write(LUCFN,'(a,f25.20,a,f25.20,a)') 4580 & 'fit(x)=',ggan(k),'*exp(-x**2*',ggas(k),')' 4581 endif 4582 else 4583 if (k.lt.i) then 4584 write(LUCFN,'(f25.20,a,f25.20,a)') 4585 & ggan(k),'*exp(-x**2*',ggas(k),') + \\' 4586 else 4587 write(LUCFN,'(f25.20,a,f25.20,a)') 4588 & ggan(k),'*exp(-x**2*',ggas(k),')' 4589 endif 4590 endif 4591 enddo 4592 write(LUCFN,'(a,i10,a)') 'set xrange [0:',nint(3d0/alpha),']' 4593 write(LUCFN,'(a,i10,a)') 'set yrange [0:1]' 4594 write(LUCFN,'(a)') 'plot fit(x), etg(x)' 4595 call gpclose(LUCFN,'KEEP') 4596 else 4597 call quit('This STG fit is not possible. Sorry.') 4598 end if 4599 end 4600C/* Deck rrffit */ 4601 subroutine rrffit(alpha,i,ggan,ggas) 4602C 4603C R12 calculation requested with fitted "r12-times-ERFC" geminal (WK/UniKA/03-29-2005). 4604C 4605#include "implicit.h" 4606 dimension ggax(6,6),ggac(6,6),ggas(6),ggan(6) 4607 character*5 fn 4608 data ggax / 4609 &0.1492D+01,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4610 &0.1266D+01,0.5257D+01,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4611 &0.1185D+01,0.3351D+01,0.1461D+02,0.0000D+00,0.0000D+00,0.0000D+00, 4612 &0.1143D+01,0.2643D+01,0.8308D+01,0.3411D+02,0.0000D+00,0.0000D+00, 4613 &0.1117D+01,0.2269D+01,0.6008D+01,0.1808D+02,0.7171D+02,0.0000D+00, 4614 &0.1099D+01,0.2037D+01,0.4815D+01,0.1239D+02,0.3603D+02,0.1405D+03 4615 &/ 4616 data ggac / 4617 &0.6939D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4618 &0.5564D+00,0.2815D+00,0.0000D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4619 &0.4835D+00,0.2677D+00,0.1493D+00,0.0000D+00,0.0000D+00,0.0000D+00, 4620 &0.4350D+00,0.2602D+00,0.1460D+00,0.9297D-01,0.0000D+00,0.0000D+00, 4621 &0.3993D+00,0.2535D+00,0.1468D+00,0.9191D-01,0.6278D-01,0.0000D+00, 4622 &0.3716D+00,0.2472D+00,0.1478D+00,0.9346D-01,0.6282D-01,0.4442D-01 4623 &/ 4624 if (i .ge. 1 .and. i .le. 6) then 4625 do k=1,i 4626 ggas(k)=ggax(k,i)*alpha**2 4627 enddo 4628 write(fn,'(a4,i1)') 'din.',i 4629 LUCFN=-1 4630 CALL GPOPEN(LUCFN,FN,'NEW',' ','FORMATTED',IDUMMY,.FALSE.) 4631 write(LUCFN,'(a,i1,a)') 'set output ''reg',i,'.ps''' 4632 write(LUCFN,'(a)') 'set term postscript' 4633 write(LUCFN,'(a,f25.20,a)') 'reg(x)=x-x*erf(x*',alpha,')' 4634 do k=1,i 4635 ggan(k)=ggac(k,i) 4636 enddo 4637 do k=1,i 4638 if (k.eq.1) then 4639 if (k.lt.i) then 4640 write(LUCFN,'(a,f25.20,a,f25.20,a)') 4641 & 'fit(x)=',ggan(k),'*x*exp(-x**2*',ggas(k),') + \\' 4642 else 4643 write(LUCFN,'(a,f25.20,a,f25.20,a)') 4644 & 'fit(x)=',ggan(k),'*x*exp(-x**2*',ggas(k),')' 4645 endif 4646 else 4647 if (k.lt.i) then 4648 write(LUCFN,'(f25.20,a,f25.20,a)') 4649 & ggan(k),'*x*exp(-x**2*',ggas(k),') + \\' 4650 else 4651 write(LUCFN,'(f25.20,a,f25.20,a)') 4652 & ggan(k),'*x*exp(-x**2*',ggas(k),')' 4653 endif 4654 endif 4655 enddo 4656 write(LUCFN,'(a,i10,a)') 'set xrange [0:',nint(30/alpha),']' 4657 write(LUCFN,'(a,F10.1,a)') 4658 & 'set yrange [0:',nint(5d0/alpha)/10.0d0,']' 4659 write(LUCFN,'(a)') 'plot fit(x), reg(x)' 4660 call gpclose(LUCFN,'KEEP') 4661 else 4662 call quit('This STG fit is not possible. Sorry.') 4663 end if 4664 end 4665C -- end of herdrv.F -- 4666