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 19*********************************************************************** 20* * 21* LUCITA, by Jeppe Olsen, DIRAC adaptation by Timo Fleig, 2001 * 22* Dalton adaptation by Jeppe Olsen, Hans Joergen Aa. Jensen and * 23* Stefan Knecht May 08 - Dec 10. * 24* * 25*********************************************************************** 26 27 SUBROUTINE DMPINT(LUINT) 28* 29* Dump integrals in WORK(KINT1),WORK(KINT2) on file LUINT 30* 31 use lucita_energy_types 32 IMPLICIT REAL*8(A-H,O-Z) 33#include "mxpdim.inc" 34#include "glbbas.inc" 35#include "wrkspc.inc" 36* 37 COMMON/CINTFO/I12S,I34S,I1234S,NINT1,NINT2,NBINT1,NBINT2 38* 39 CALL REWINO(LUINT) 40*.1 : One-electron integrals 41 WRITE(LUINT,'(E22.15)') 42 & (WORK(KINT1-1+INT1),INT1=1,NINT1) 43*.2 : Two-electron integrals 44 WRITE(LUINT,'(E22.15)') 45 & (WORK(KINT2-1+INT2),INT2=1,NINT2) 46*.3 : Core energy 47 WRITE(LUINT,'(E22.15)') ECORE_ORIG 48* 49 RETURN 50 END 51*********************************************************************** 52 SUBROUTINE GET_CMOAO(CMO) 53* 54* Obtain AO-MO transformation matrix 55* 56* Jeppe Olsen, November 1997 57* 58 IMPLICIT REAL*8(A-H,O-Z) 59* 60#include "priunit.h" 61* 62#include "mxpdim.inc" 63#include "crun.inc" 64#include "cgas.inc" 65#include "lucinp.inc" 66#include "clunit.inc" 67#include "orbinp.inc" 68*. Output 69 DIMENSION CMO(*) 70 71! write(lupri,*) 'GET_CMOAO: ENVIRO = ',ENVIRO !hjaaj DEBUG 72 IF(ENVIRO(1:6).EQ.'DALTON') THEN 73 CALL GET_CMOAO_DALTON(CMO,NMOS_ENV(1),NAOS_ENV(1),NSMOB) 74 ELSE IF(ENVIRO(1:6).EQ.'MOLCAS') THEN 75*. 76 CALL QUIT('MOLCAS environment not available in this version') 77 ELSE IF(ENVIRO(1:5).EQ.'LUCIA' ) THEN 78*. Read in from LUCIA 1e file : unit 91 79 LU91 = 91 80 CALL GET_CMOAO_LUCIA(CMO,NMOS_ENV,NAOS_ENV,LU91) 81 ELSE IF(ENVIRO(1:4).EQ.'NONE') THEN 82 WRITE(lupri,*) ' GET_CMOAO, Warning : Called with ENVIRO = NONE' 83 WRITE(lupri,*) ' No coefficients read in ' 84 END IF 85* 86 RETURN 87 END 88*********************************************************************** 89 SUBROUTINE GET_CMOAO_DALTON(CMO,NBAS,NMO,NSM) 90* 91* Obtain MO-AO expansion matrix from SIRIUS/DALTON file SIRGEOM 92* 93* Jeppe Olsen, June 1997 94* 95 IMPLICIT REAL*8(A-H,O-Z) 96* 97#include "priunit.h" 98* 99*. Input 100 INTEGER NBAS(*), NMO(*) 101*. Output 102 DIMENSION CMO(*) 103* 104 105! write(lupri,*) 'GET_CMOAO_DALTON: NSM = ',NSM !hjaaj DEBUG 106! write(lupri,*) 'GET_CMOAO_DALTON: NBAS = ',NBAS(1:NSM) !hjaaj DEBUG 107! write(lupri,*) 'GET_CMOAO_DALTON: NMO = ',NMO (1:NSM) !hjaaj DEBUG 108 NTEST = 0 109 ITAP30 = 16 110 OPEN(ITAP30,STATUS='OLD',FORM='UNFORMATTED',FILE='SIRIFC') 111 REWIND ITAP30 112 CALL MOLLAB('TRCCINT ',ITAP30,6) 113*. Skip record containing dimensions of orbitals 114 READ (ITAP30) NSYM 115*. And skip record containing eigenvalues etc 116 READ (ITAP30) 117C READ (ITAP30) NSYMHF,NORBT,NBAST,NCMOT,(NOCC(I),I=1,NSYMHF), 118C * (NORB(I),I=1,NSYMHF),(NBAS(I),I=1,NSYMHF), 119C * POTNUC,EMCSCF 120C 121C 122C READ (ITAP30) (WRK(KEIGVL+I-1),I=1,NORBT), 123C * (IWRK(KEIGSY+I-1),I=1,NORBT) 124*. And then the MO-AO expansion matrix 125 NCOEF = 0 126 DO ISM = 1, NSM 127 NCOEF = NCOEF + NMO(ISM)*NBAS(ISM) 128 END DO 129 READ (ITAP30) (CMO(I),I=1,NCOEF) 130 CLOSE(ITAP30,STATUS='KEEP') 131* 132 NTEST = 000 ! 133! NTEST = 100 ! hjaaj DEBUG 134 IF(NTEST.GE.100) THEN 135 WRITE(lupri,*) ' MO - AO expansion matrix ' 136 WRITE(lupri,*) '=============================' 137 WRITE(lupri,*) 138 CALL APRBLM2(CMO,NBAS,NMO,NSM,0) 139 END IF 140* 141 RETURN 142 END 143*********************************************************************** 144 SUBROUTINE GET_CMOAO_LUCIA(CMO,NMOS,NAOS,LUH) 145* 146* Obtain CMOAO expansion matrix from LUCIA formatted file LUH 147* 148 IMPLICIT REAL*8(A-H,O-Z) 149* 150#include "priunit.h" 151* 152#include "mxpdim.inc" 153#include "lucinp.inc" 154#include "orbinp.inc" 155#include "crun.inc" 156*. Input 157 INTEGER NMOS(*),NAOS(*) 158*. Output 159 DIMENSION CMO(*) 160* 161* Structure of file 162* 1 : Number of syms 163* 2 : NMO's per sym 164* 3 : NAO's per SYM 165* 4 : Number of elements in CMOAO 166* Note : CMOAO and property integrals written in form 167* given by ONEEL_MAT_DISC 168* 169* Jeppe Olsen, Feb. 98 170* 171 WRITE(lupri,*) ' GET_CMOAO_LUCIA, LUH = ', LUH 172 CALL REWINO(LUH) 173*. skip Number of orbital symmetries 174 READ(LUH,*) 175*. skip Number of MO's per symmetry 176 READ(LUH,*) 177*. skip Number of AO's per symmetry 178 READ(LUH,*) 179*. skip read Length of CMO-AO expansion 180 READ(LUH,*) 181*. read CMO-AO expansion matrix 182 CALL ONEEL_MAT_DISC(CMO,1,NSMOB,NAOS,NMOS,LUH,1) 183C ONEEL_MAT_DISC(H,IHSM,NSM,NRPSM,NCPSM,LUH,IFT) 184* 185 NTEST = 000 186 IF(NTEST.GE.100) THEN 187 WRITE(lupri,*) ' MO-AO transformation read in ' 188 CALL PRHONE(CMO,NMOS,1,NSMOB,0) 189C PRHONE(C,NFUNC,M,NSM,IPACK) 190 END IF 191* 192 RETURN 193 END 194*********************************************************************** 195 SUBROUTINE GET_CMOMO(CMOMO) 196* 197* Obtain MO-MO transformation matrix CMOMO for transforming to 198* final set of orbitals 199* 200* Output matrix CMOMO is returned in symmetry packed form 201* 202*. Density matrix is assumed in place 203* 204* Type of final orbitals is provided by the keyword 205* keywords ITRACI_CR, ITRACI_CN 206* 207* ITRACI_CR : COMP => Rotate all orbitals 208* REST => Rotalte only inside orbital subspaces 209* 210* ITRACI_CN : NATU => Transform to natural orbitals 211* ITRACI_CR : CANO => Transform to canonical orbitals 212* 213* Jeppe Olsen, February 1998 ( from FINMO) 214* 215 IMPLICIT REAL*8(A-H,O-Z) 216* 217#include "priunit.h" 218* 219 INTEGER*8 KMAT1, KMAT2, KMAT3, KMAT4 220 ! for addressing of WORK 221#include "mxpdim.inc" 222#include "wrkspc.inc" 223#include "crun.inc" 224#include "glbbas.inc" 225#include "orbinp.inc" 226#include "lucinp.inc" 227#include "cgas.inc" 228 COMMON/CINTFO/I12S,I34S,I1234S,NINT1,NINT2,NBINT1,NBINT2 229*. Output 230 DIMENSION CMOMO(*) 231* 232 NTEST = 100 233 IF(NTEST.GE.1) THEN 234 WRITE(lupri,*) 235 WRITE(lupri,*) ' ====================' 236 WRITE(lupri,*) ' GET_CMOMO in action' 237 WRITE(lupri,*) ' ====================' 238 WRITE(lupri,*) 239 END IF 240 241 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'GETCMO') 242 CALL MEMMAN(KMAT1,NTOOB**2,'ADDL ',2,'MAT1 ') 243 CALL MEMMAN(KMAT2,NTOOB**2,'ADDL ',2,'MAT2 ') 244 CALL MEMMAN(KMAT3,NTOOB**2,'ADDL ',2,'MAT3 ') 245 CALL MEMMAN(KMAT4,2*NTOOB**2,'ADDL ',2,'MAT4 ') 246* 247*. Matrix defining final orbitals 248* 249 IF(ITRACI_CN(1:4).EQ.'CANO' ) THEN 250*. Construct FI+FA in WORK(KMAT1) 251 CALL COPVEC(WORK(KINT1O),WORK(KMAT1),NINT1) 252 CALL FIFAM(WORK(KMAT1)) 253 ELSE IF(ITRACI_CN(1:4).EQ.'NATU' ) THEN 254*. Symmetry order density matrix 255 CALL TYPE_TO_SYM_REO_MAT(WORK(KRHO1),WORK(KMAT2)) 256*. Pack to triangular form 257 CALL TRIPAK_BLKM(WORK(KMAT2),WORK(KMAT1),1,NTOOBS,NSMOB) 258*. multiply by minus one to get natural orbitals 259*. with largest occupations first 260 ONEM = -1.0D0 261 LDIM = 0 262 DO ISM = 1, NSMOB 263 LDIM = LDIM + NTOOBS(ISM)*(NTOOBS(ISM)+1)/2 264 END DO 265 CALL SCALVE(WORK(KMAT1),ONEM,LDIM) 266 IF(NTEST.GE.100) THEN 267 WRITE(lupri,*) ' Packed density matrix ( times - 1 )' 268 CALL APRBLM2(WORK(KMAT1),NACOBS,NACOBS,NSMOB,1) 269 END IF 270 END IF 271* 272* Diagonalize 273* 274 IF(ITRACI_CR(1:4).EQ.'REST') THEN 275*. Diagonalize symmetry-type blocks 276 CALL DIAG_BLKS(WORK(KMAT1),CMOMO,NGSOB,NTOOBS,MXPOBS, 277 & NSMOB,NGAS,WORK(KMAT3),WORK(KMAT4)) 278*. Reorder to assure max diag dominance 279 IREO = 1 280 IF(IREO.NE.0) THEN 281 WRITE(lupri,*) ' CMOMO reordered to assure max. diag. dom.' 282 DO ISM = 1, NSMOB 283 IF(ISM.EQ.1) THEN 284 IOFF = 1 285 ELSE 286 IOFF = IOFF + NTOOBS(ISM-1)**2 287 END IF 288 L = NTOOBS(ISM) 289 CALL GET_DIAG_DOM(CMOMO(IOFF),WORK(KMAT1),L,WORK(KMAT2)) 290 CALL COPVEC(WORK(KMAT1),CMOMO(IOFF),L*L) 291 END DO 292 END IF 293 ELSE IF (ITRACI_CR(1:4).EQ.'COMP') THEN 294*. Diagonalize symmetry blocks 295 CALL DIAG_BLKS(WORK(KMAT1),CMOMO,NACOBS,NTOOBS,MXPOBS, 296 & NSMOB,1,WORK(KMAT3),WORK(KMAT4)) 297*. Reorder to assure max diag dominance 298 IREO = 1 299 IF(IREO.NE.0) THEN 300 WRITE(lupri,*) ' CMOMO reordered to assure max. diag. dom.' 301 DO ISM = 1, NSMOB 302 IF(ISM.EQ.1) THEN 303 IOFF = 1 304 ELSE 305 IOFF = IOFF + NTOOBS(ISM-1)**2 306 END IF 307 L = NTOOBS(ISM) 308 CALL GET_DIAG_DOM(CMOMO(IOFF),WORK(KMAT1),L,WORK(KMAT2)) 309 CALL COPVEC(WORK(KMAT1),CMOMO(IOFF),L*L) 310 END DO 311 END IF 312 END IF 313* 314 IF(NTEST.GE.100) THEN 315 WRITE(lupri,*) ' Output set of MO''s ' 316 CALL APRBLM2(CMOMO,NTOOBS,NTOOBS,NSMOB,0) 317 END IF 318* 319 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'GETCMO') 320 RETURN 321 END 322*********************************************************************** 323* 324* Obtain property integrals with LABEL LABEL from LU91, 325* LUCIA format 326* 327* Jeppe Olsen, Feb.98 328 329 SUBROUTINE GET_H1AO(LABEL,H1AO,IHSM,NBAS) 330* 331* Obtain 1 electron integrals with label LABEL 332* 333* Jeppe Olsen, Feb.98 334 IMPLICIT REAL*8(A-H,O-Z) 335 REAL*8 H1AO(*) 336 INTEGER IHSM, NBAS(*) 337* 338#include "priunit.h" 339* 340 INTEGER*8 KLSCR 341 ! for addressing of WORK 342* 343#include "mxpdim.inc" 344#include "crun.inc" 345#include "orbinp.inc" 346#include "wrkspc.inc" 347#include "lucinp.inc" 348* 349 CHARACTER*8 LABEL 350* 351 IDUM = 0 352 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'GT_H1A') 353* 354 write(lupri,*) 'GET_H1AO: ENVIRO = ',ENVIRO !hjaaj DEBUG 355 IF(ENVIRO(1:6).EQ.'DALTON') THEN 356 LSCR = NTOOB**2 357 CALL MEMMAN(KLSCR,LSCR,'ADDL ',2,'GTH1SC') 358 CALL GET_H1AO_DALTON(LABEL,H1AO,IHSM,WORK(KLSCR),NBAS,NSMOB) 359C GET_H1AO_DALTON(LABEL,H1AO,IHSM,SCR,NBAS,NSM) 360 ELSE IF (ENVIRO(1:5).EQ.'LUCIA') THEN 361 LU91 = 91 362 CALL GET_H1AO_LUCIA(LABEL,H1AO,LU91) 363 END IF 364* 365 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'GT_H1A') 366* 367 RETURN 368 END 369*********************************************************************** 370 SUBROUTINE GET_H1AO_DALTON(LABEL,H1AO,IHSM,SCR,NBAS,NSM) 371* 372*. Obtain one-electron integrals in ao basis from dalton 373* 374* Label of integrals LABEL from FILE AORPROPER 375* 376* Jeppe Olsen 377* 378 IMPLICIT REAL*8(A-H,O-Z) 379* 380#include "priunit.h" 381* 382*. Input 383 CHARACTER*8 LABEL 384 INTEGER NBAS(*) 385#include "multd2h.inc" 386*. output 387 DIMENSION H1AO(*) 388*. Scratch 389 DIMENSION SCR(*) 390* 391 LOGICAL FNDLAB 392* 393 write(lupri,*) 'GET_H1AO_DALTON: LABEL = ',LABEL !hjaaj DEBUG 394 NTEST = 02 395 IF(NTEST.GE.2) THEN 396 WRITE(lupri,*) ' Fetching one-electron integrals with Label ', 397 & LABEL 398 WRITE(lupri,*) ' IHSM NSM', IHSM,NSM 399 END IF 400* 401*. Number of elements in : Complete lower half array 402* Symmetry restricted complete matrix 403* Symmetry restricted lower half matrix 404*-- I am not completely sure about the input format of the integrals 405 NBAST = 0 406 DO ISM = 1, NSM 407 NBAST = NBAST + NBAS(ISM) 408 END DO 409 NINT01 = NBAST*(NBAST+1)/2 410C write(lupri,*) ' IHSM = ', IHSM 411* 412 NINT10 = 0 413 DO IRSM = 1, NSM 414 ICSM = MULTD2H(IHSM,IRSM) 415 NINT10 = NINT10 + NBAS(IRSM)*NBAS(ICSM) 416 END DO 417* 418 NINT11 = 0 419 DO IRSM = 1, NSM 420 ICSM = MULTD2H(IHSM,IRSM) 421 IF(IRSM.GT.ICSM) THEN 422 NINT11 = NINT11 + NBAS(IRSM)*NBAS(ICSM) 423 ELSE IF(IRSM.EQ.ICSM) THEN 424 NINT11 = NINT11 + NBAS(IRSM)*(NBAS(IRSM)+1)/2 425 END IF 426 END DO 427* 428*. Read in integrals, assumed in complete lower half format 429* 430 LUPRP = 15 431 OPEN (LUPRP,STATUS='OLD',FORM='UNFORMATTED',FILE='AOPROPER') 432 REWIND (LUPRP) 433 IF (FNDLAB(LABEL,LUPRP)) THEN 434C write(lupri,*) ' Label obtained' 435 READ(LUPRP) (SCR(I),I=1,NINT01) 436C write(lupri,*) 'integrals readin' 437C call prsym(scr,NBAST) 438C CALL READT(LUPRP,NBAST*(NBAST+1)/2,WRK(KSCR2)) 439 ELSE 440 WRITE(lupri,*)'Property label "',LABEL,'" not found on file' 441 Call Abend2( 'Wrong input or integrals not generated' ) 442 ENDIF 443 CLOSE(LUPRP,STATUS='KEEP') 444* 445C WRITE(lupri,*) ' Number of symmetry apdapted integrals',NINT10 446* 447*. Transfer integrals to symmetry adapted form, complete form 448* 449 IBINT = 1 450*. Loop over symmetry blocks 451 DO IRSM = 1, NSM 452 ICSM = MULTD2H(IHSM,IRSM) 453 NR = NBAS(IRSM) 454 NC = NBAS(ICSM) 455*. Offsets 456 IBR = 1 457 DO ISM = 1, IRSM - 1 458 IBR = IBR + NBAS(ISM) 459 END DO 460 IBC = 1 461 DO ISM = 1, ICSM - 1 462 IBC = IBC + NBAS(ISM) 463 END DO 464*. Complete block, stored in usual column wise fashion 465 DO ICORB = 1, NC 466 DO IRORB = 1, NR 467 ICABS = IBC + ICORB -1 468 IRABS = IBR + IRORB -1 469 ICRMX = MAX(ICABS,IRABS) 470 ICRMN = MIN(ICABS,IRABS) 471 H1AO(IBINT-1 + (ICORB-1)*NR+IRORB) = 472 & SCR(ICRMX*(ICRMX-1)/2+ICRMN) 473 END DO 474 END DO 475 IBINT = IBINT + NR*NC 476 END DO 477* 478 IF(NTEST.GE.100) THEN 479 WRITE(lupri,*) ' One-electron integrals obtained from AOPROPER' 480 CALL PRSYM(SCR,NBAST) 481* 482 WRITE(lupri,*) ' One-electron integrals in packed form' 483 CALL PRHONE(H1AO,NBAS,IHSM,NSM,0) 484C PRHONE(H,NFUNC,IHSM,NSM,IPACK) 485 END IF 486* 487 RETURN 488 END 489*********************************************************************** 490 SUBROUTINE GET_H1AO_LUCIA(LABEL,H1,LUH) 491* 492* 493* Obtain property integrals with LABEL LABEL from LU91, 494* LUCIA format 495* 496* Jeppe Olsen, Feb.98 497* 498 IMPLICIT REAL*8(A-H,O-Z) 499* 500#include "priunit.h" 501* 502#include "mxpdim.inc" 503#include "lucinp.inc" 504#include "orbinp.inc" 505#include "crun.inc" 506#include "wrkspc.inc" 507C CHARACTER*1 XYZ(3) 508C DATA XYZ/'X','Y','Z'/ 509 CHARACTER*8 LABEL, LABEL2, LABELX 510*. Output 511 DIMENSION H1(*) 512* 513* Structure of file 514* 1 : Number of syms 515* 2 : NMO's per sym 516* 3 : NAO's per SYM 517* 4 : Number of elements in CMOAO 518* 4 : CMOAO-expansion matrix (in symmetry packed form) 519* 5 : Number of property AO lists 520* Loop over number of properties 521* Label, offset and length of each proprty list 522* 523* Property integrals for prop1,prop2 ... 524* 525* Note : CMOAO and property integrals written in form 526* given by ONEEL_MAT_DISC 527* 528* Jeppe Olsen, Feb. 98 529* 530 IDUM = 0 531 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'GETH1A') 532* 533*. DIPOLE => DIPLEN 534 IF(LABEL(1:6).EQ.'DIPOLE') THEN 535 LABELX = 'DIPLEN ' 536 ELSE 537 LABELX = LABEL 538 END IF 539* 540 CALL REWINO(LUH) 541*. Skip Number of orbital symmetries 542 READ (LUH,*) NSMOB 543*. Skip Number of MO's per symmetry 544 READ (LUH,*) (NMOS_ENV(ISM),ISM=1,NSMOB) 545*. Skip Number of AO's per symmetry 546 READ (LUH,*) (NAOS_ENV(ISM),ISM=1,NSMOB) 547*. Length of CMO-AO expansion 548 READ(LUH,*) LENGTH 549*. And skip 550 DO IJ = 1, LENGTH 551 READ(LUH,'(E22.15)') 552 END DO 553*. Total number of properties ( 3 for each rank1, 6 for each rank 2) 554 READ(LUH,*) NPROP_COMP 555 IFOUND = 0 556 WRITE(lupri,*) ' NPROP_COMP = ', NPROP_COMP 557 DO IPROP_COM = 1, NPROP_COMP 558 READ(LUH,'(A,I6,I6)') LABEL2,IOFF,LENGTH 559 IF(LABEL2.EQ.LABELX) THEN 560 IOFFA = IOFF 561 LENGTHA = LENGTH 562 IFOUND = 1 563 END IF 564 END DO 565 IF(IFOUND.EQ.0) THEN 566 WRITE(lupri,*) ' Label not found on file 91' 567 WRITE(lupri,'(2A)' ) ' Label = ', LABELX 568 Call Abend2( ' Label not found on file 91' ) 569 END IF 570*. Skip to start of integrals 571 WRITE(lupri,*) ' IOFFA, LENGTHA ', IOFFA,LENGTHA 572 DO IJ = 1, IOFFA - 1 573 READ(LUH,*) 574 END DO 575*. and read 576 CALL SYM_FOR_OP(LABEL,IXYZSYM,IOPSM) 577 CALL ONEEL_MAT_DISC(H1,IOPSM,NSMOB, 578 & NAOS_ENV,NAOS_ENV,LUH,1) 579C ONEEL_MAT_DISC(H,IHSM,NSM,NRPSM,NCPSM,LUH,IFT) 580* 581 NTEST = 000 582 IF(NTEST.GE.100) THEN 583 WRITE(lupri,*) ' Property integrals read in ' 584 CALL PRHONE(H1,NAOS_ENV,IOPSM,NSMOB,0) 585C PRHONE(H,NFUNC,IHSM,NSM,IPACK) 586 END IF 587* 588 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'GETH1A') 589 RETURN 590 END 591*********************************************************************** 592 SUBROUTINE GET_ORB_DIM_ENV(ECORE_ENV) 593* 594* Obtain number of orbitals and basis functions from the 595* programming environment. 596* results stored in NAOS_ENV, NMOS_ENV 597* 598* Obtain environments CORE energy, ECORE_ENV 599* 600* Jeppe Olsen, December 97 601* 602#ifdef VAR_MPI 603 use dalton_mpi 604#endif 605#include "implicit.h" 606#include "priunit.h" 607#include "maxorb.h" 608#include "infpar.h" 609#ifdef VAR_MPI 610#include "mpif.h" 611 integer(kind=MPI_INTEGER_KIND) my_STATUS(MPI_STATUS_SIZE) 612#endif 613#include "mxpdim.inc" 614#include "orbinp.inc" 615#include "lucinp.inc" 616#include "parluci.h" 617! potnuc on infinp.h could also deliver ecore_env 618! 619 NTEST = 000 620! NTEST = 9999 ! debug 621 622 call izero(naos_env,MXPOBS) 623 call izero(nmos_env,MXPOBS) 624 if(luci_myproc .eq. luci_master)then 625 CALL GETOBS_DALTON(ECORE_ENV,NAOS_ENV,NMOS_ENV,NSMOB) 626 end if 627#ifdef VAR_MPI 628 if(luci_nmproc .gt. 1)then 629 call dalton_mpi_bcast(nsmob, luci_master, mpi_comm_world) 630 call dalton_mpi_bcast(ecore_env, luci_master, mpi_comm_world) 631 call dalton_mpi_bcast(naos_env, luci_master, mpi_comm_world) 632 call dalton_mpi_bcast(nmos_env, luci_master, mpi_comm_world) 633 end if 634#endif 635! 636 IF(NTEST.GE.100) THEN 637 WRITE(lupri,*) ' From GET_ORB_FROM_ENV : ' 638 WRITE(lupri,*) ' =======================' 639 WRITE(lupri,*) ' NSMOB ',NSMOB 640 WRITE(lupri,*) ' NAOS_ENV' 641 CALL IWRTMA(NAOS_ENV,1,NSMOB,1,NSMOB) 642 WRITE(lupri,*) ' NMOS_ENV' 643 CALL IWRTMA(NMOS_ENV,1,NSMOB,1,NSMOB) 644 WRITE(lupri,*) ' ECORE_ENV=', ECORE_ENV 645 END IF 646! 647 END 648*********************************************************************** 649* 650 SUBROUTINE GET_PROPINT(H,IHSM,LABEL,SCR,NMO,NBAS,NSM,ILOW) 651* 652*. Obtain Property integrals in MO basis for operator with 653* label LABEL. 654* 655* If ILOW = 1, only the elements below the diagonal are 656* obtained. 657* 658* Jeppe Olsen, June 1997 659* September 97 : ILOW added 660 IMPLICIT REAL*8(A-H,O-Z) 661* 662#include "priunit.h" 663* 664*. Input 665 DIMENSION NMO(*),NBAS(*) 666#include "multd2h.inc" 667*. Output 668 DIMENSION H(*) 669*. Scratch 670 DIMENSION SCR(*) 671*. Scratch should atleaest be of length ** 672* 673 NTEST = 000 674*. Integrals in AO basis, neglect symmetry 675 NBAST = 0 676 NMOT = 0 677 DO ISM = 1, NSM 678 NBAST = NBAST + NBAS(ISM) 679 NMOT = NMOT + NMO(ISM) 680 END DO 681C? WRITE(lupri,*) ' Total number of basis functions ',NBAST 682 LINTMX = NBAST*NBAST 683* 684 KLH1AO = 1 685 KLFREE = KLH1AO + LINTMX 686* 687 KLC = KLFREE 688 KLFREE = KLC + LINTMX 689* 690 KLSCR = KLFREE 691*. Currently only DALTON route is working 692 IDALTON = 1 693 IF(IDALTON.EQ.1) THEN 694C? WRITE(lupri,*) ' Dalton route in action' 695*. Obtain AO property integrals 696C GET_H1AO_DALTON(LABEL,H1AO,IHSM,SCR,NBAS) 697C GET_H1AO(LABEL,H1AO,IHSM,NBAS) 698 CALL GET_H1AO(LABEL,SCR(KLH1AO),IHSM,NBAS) 699C CALL GET_H1AO_DALTON(LABEL,SCR(KLH1AO),IHSM, 700C & SCR(KLSCR),NBAS,NSM) 701*. Obtain MO-AO transformation matrix 702 CALL GET_CMOAO(SCR(KLC)) 703*. Transform from AO to MO basis 704C TRAH1(NBAS,NORB,NSYM,HAO,C,HMO,IHSM,SCR) 705 CALL TRAH1(NBAS,NMO,NSM,SCR(KLH1AO),SCR(KLC),H,IHSM, 706 & SCR(KLSCR)) 707 END IF 708* 709 IF(NTEST .GE. 100 ) THEN 710 WRITE(lupri,*) 'electron integrals in MO basis, full format ' 711 CALL PRHONE(H,NMO,IHSM,NSM,0) 712 END IF 713 IF(ILOW.EQ.1) THEN 714*. Complete to lower half form 715 IOFF_IN = 1 716 IOFF_OUT = 1 717 DO ISM = 1, NSM 718 JSM = MULTD2H(ISM,IHSM) 719 IF(ISM.EQ.JSM) THEN 720*. Copy lower half 721 LDIM = NMO(ISM) 722 NELMNT_IN = LDIM * LDIM 723 NELMNT_OUT = LDIM * (LDIM + 1)/2 724 CALL COPVEC(H(IOFF_IN),SCR(KLSCR),NELMNT_IN) 725 SIGN = 1.0D0 726 CALL TRIPAK_LUCI(SCR(KLSCR),H(IOFF_OUT),1,LDIM,LDIM,SIGN) 727 IOFF_IN = IOFF_IN + NELMNT_IN 728 IOFF_OUT = IOFF_OUT + NELMNT_OUT 729 ELSE IF(ISM.LT.JSM) THEN 730*. Just skip block in input matrix 731 LIDIM = NMO(ISM) 732 LJDIM = NMO(JSM) 733 IOFF_IN = IOFF_IN + LIDIM*LJDIM 734 ELSE IF(ISM.GT.JSM) THEN 735*. Copy block to block 736 LIDIM = NMO(ISM) 737 LJDIM = NMO(JSM) 738 NELMNT = LIDIM*LJDIM 739C CALL TRPMAT(H(IOFF_IN),LIDIM,LJDIM,H(IOFF_OUT)) 740 CALL COPVEC(H(IOFF_IN),H(IOFF_OUT),NELMNT) 741 IOFF_IN = IOFF_IN + NELMNT 742 IOFF_OUT = IOFF_OUT + NELMNT 743 END IF 744 END DO 745 END IF 746*. The one-electron integrals reside in a NMOT X NMOT matrix. 747*. Zero trivial integrals 748 IF(ILOW.EQ.1) THEN 749 NELMNT = IOFF_OUT-1 750 ELSE 751 LENGTH = 0 752 DO ISM = 1, NSM 753 JSM = MULTD2H(ISM,IHSM) 754 NELMNT = NELMNT + NMO(ISM)*NMO(JSM) 755 END DO 756 IFREE = NELMNT + 1 757 END IF 758C? WRITE(lupri,*) ' GET_PROP : NELMNT= ', NELMNT 759 ZERO = 0.0D0 760 NZERO = NMOT*NMOT - NELMNT 761 IFREE = NELMNT + 1 762 CALL SETVEC(H(IFREE),ZERO,NZERO) 763 764 IF(NTEST .GE. 50 ) THEN 765 WRITE(lupri,*) 'electron integrals in MO basis ' 766 CALL PRHONE(H,NMO,IHSM,NSM,ILOW) 767 END IF 768* 769 RETURN 770 END 771*********************************************************************** 772 SUBROUTINE GETH1(H,ISM,ITP,JSM,JTP) 773* 774* One-electron integrals over orbitals belonging to 775* given OS class 776* 777* 778* The orbital symmetries are used to obtain the total 779* symmetry of the one-electron integrals. 780* It is therefore assumed that ISM, JSM represents 781* a correct symmetry block 782* of the integrals 783* 784* Jeppe Olsen, Version of fall 97 785* Summer of 98 : CC options added 786* 787 IMPLICIT REAL*8(A-H,O-Z) 788* 789#include "priunit.h" 790* 791#include "mxpdim.inc" 792#include "wrkspc.inc" 793*.Global pointers 794#include "glbbas.inc" 795#include "lucinp.inc" 796#include "orbinp.inc" 797#include "cc_exc.inc" 798*.Output 799 DIMENSION H(*) 800* 801 NI = NOBPTS(ITP,ISM) 802 NJ = NOBPTS(JTP,JSM) 803* 804 IF(ICC_EXC.EQ.0) THEN 805* 806* Normal one-electron integrals 807* 808 IJ = 0 809 DO J = 1, NJ 810 DO I = 1, NI 811 IJ = IJ+1 812 H(IJ) = GETH1E(I,ITP,ISM,J,JTP,JSM) 813 END DO 814 END DO 815 ELSE 816* 817* Single excitation coefficients dressed up as integrals 818* taken from KCC 819C GET_SX_BLK(HBLK,H,IGAS,ISM,JGAS,JSM) 820*. Note : WORK(KCC1) not perfect choice 821! CALL GET_SX_BLK(H,WORK(KCC1),ITP,ISM,JTP,JSM) 822 END IF 823* 824 NTEST = 0 825 IF(NTEST.NE.0) THEN 826 WRITE(lupri,*) ' H1 for itp ism jtp jsm ',ITP,ISM,JTP,JSM 827 CALL WRTMT_LU(H,NI,NJ,NI,NJ) 828 END IF 829* 830 RETURN 831 END 832*********************************************************************** 833 FUNCTION GETH1E(IORB,ITP,ISM,JORB,JTP,JSM) 834* 835* One-electron integral for active 836* orbitals (IORB,ITP,ISM),(JORB,JTP,JSM) 837* 838* The orbital symmetries are used to obtain the 839* total symmetry of the operator 840 IMPLICIT REAL*8(A-H,O-Z) 841#include "mxpdim.inc" 842C COMMON/BIGGY/WORK(MXPWRD) 843#include "wrkspc.inc" 844* 845#include "glbbas.inc" 846#include "lucinp.inc" 847#include "orbinp.inc" 848#include "multd2h.inc" 849#include "intform.inc" 850* 851 IJSM = MULTD2H(ISM,JSM) 852 IF(IH1FORM.EQ.1) THEN 853*. Normal integrals, lower triangular packed 854 IF(IJSM.EQ.1) THEN 855 GETH1E = 856 & GTH1ES(IREOTS,WORK(KPINT1),WORK(KINT1),IBSO,MXPNGAS, 857 & IOBPTS,NACOBS,IORB,ITP,ISM,JORB,JTP,JSM,1) 858 ELSE 859 GETH1E = 860 & GTH1ES(IREOTS,WORK(KPGINT1(IJSM)),WORK(KINT1),IBSO,MXPNGAS, 861 & IOBPTS,NACOBS,IORB,ITP,ISM,JORB,JTP,JSM,1) 862 END IF 863 ELSE 864*. Integrals are in full blocked form 865 GETH1E = 866 & GTH1ES(IREOTS,WORK(KPGINT1A(IJSM)),WORK(KINT1),IBSO,MXPNGAS, 867 & IOBPTS,NACOBS,IORB,ITP,ISM,JORB,JTP,JSM,0) 868 END IF 869* 870 RETURN 871 END 872*********************************************************************** 873 FUNCTION GETH1I(IORB,JORB) 874* 875* Obtain one -electron integral H(IORB,JOB) 876* 877* Interface from EXPHAM to LUCIA 878 IMPLICIT REAL*8 (A-H,O-Z) 879* 880#include "priunit.h" 881* 882#include "mxpdim.inc" 883#include "orbinp.inc" 884* 885 ISM = ISMFTO(IORB) 886 ITP = ITPFSO(IREOTS(IORB)) 887 IREL = IORB - IOBPTS(ITP,ISM) + 1 888* 889 JSM = ISMFTO(JORB) 890 JTP = ITPFSO(IREOTS(JORB)) 891 JREL = JORB - IOBPTS(JTP,JSM) + 1 892* 893 GETH1I = GETH1E(IREL,ITP,ISM,JREL,JTP,JSM) 894* 895 NTEST = 0 896 IF( NTEST .NE. 0 ) THEN 897 WRITE(lupri,*) ' GETH1I : IORB JORB ', IORB, JORB 898 WRITE(lupri,*) ' ISM ITP IREL ', ISM,ITP,IREL 899 WRITE(lupri,*) ' JSM JTP JREL ', JSM,JTP,JREL 900 WRITE(lupri,*) ' GETH1I = ', GETH1I 901 END IF 902* 903 RETURN 904 END 905*********************************************************************** 906 SUBROUTINE GETINCN2(XINT,ITP,ISM,JTP,JSM,KTP,KSM,LTP,LSM, 907 & IXCHNG,IKSM,JLSM,INTLST,IJKLOF,NSMOB,I2INDX, 908 & ICOUL) 909* 910* Obtain integrals 911* 912* ICOUL = 0 : 913* XINT(IK,JL) = (IJ!KL) for IXCHNG = 0 914* = (IJ!KL)-(IL!KJ) for IXCHNG = 1 915* 916* ICOUL = 1 : 917* XINT(IJ,KL) = (IJ!KL) for IXCHNG = 0 918* = (IJ!KL)-(IL!KJ) for IXCHNG = 1 919* 920* ICOUL = 2 : XINT(IL,JK) = (IJ!KL) for IXCHNG = 0 921* = (IJ!KL)-(IL!KJ) for IXCHNG = 1 922* 923* Storing for ICOUL = 1 not working if IKSM or JLSM .ne. 0 924* 925* 926* Version for integrals stored in INTLST 927* 928* If type equals zero, all integrals of given type are fetched 929* ( added aug8, 98) 930* 931 IMPLICIT REAL*8(A-H,O-Z) 932* 933#include "mxpdim.inc" 934#include "orbinp.inc" 935*. Integral list 936 Real * 8 Intlst(*) 937 Dimension IJKLof(NsmOB,NsmOb,NsmOB) 938*. Pair of orbital indices ( symmetry ordered ) => address in symmetry packed 939*. matrix 940 Dimension I2INDX(*) 941*.Output 942 DIMENSION XINT(*) 943*. Local scratch 944 DIMENSION IJARR(MXPORB) 945* 946 IF(ITP.GE.1) THEN 947 iOrb=NOBPTS(ITP,ISM) 948 ELSE 949 IORB = NTOOBS(ISM) 950 END IF 951* 952 IF(JTP.GE.1) THEN 953 jOrb=NOBPTS(JTP,JSM) 954 ELSE 955 JORB = NTOOBS(JSM) 956 END IF 957* 958 IF(KTP.GE.1) THEN 959 kOrb=NOBPTS(KTP,KSM) 960 ELSE 961 KORB = NTOOBS(KSM) 962 END IF 963* 964 IF(LTP.GE.1) THEN 965 lOrb=NOBPTS(LTP,LSM) 966 ELSE 967 LORB = NTOOBS(LSM) 968 END IF 969* 970*. Offsets relative to start of all orbitals, symmetry ordered 971 IOFF = IBSO(ISM) 972 DO IITP = 1, ITP -1 973 IOFF = IOFF + NOBPTS(IITP,ISM) 974 END DO 975* 976 JOFF = IBSO(JSM) 977 DO JJTP = 1, JTP -1 978 JOFF = JOFF + NOBPTS(JJTP,JSM) 979 END DO 980* 981 KOFF = IBSO(KSM) 982 DO KKTP = 1, KTP -1 983 KOFF = KOFF + NOBPTS(KKTP,KSM) 984 END DO 985* 986 LOFF = IBSO(LSM) 987 DO LLTP = 1, LTP -1 988 LOFF = LOFF + NOBPTS(LLTP,LSM) 989 END DO 990 991* 992* Collect Coulomb terms 993* 994 ijblk = max(ism,jsm)*(max(ism,jsm)-1)/2 + min(ism,jsm) 995 klblk = max(ksm,lsm)*(max(ksm,lsm)-1)/2 + min(ksm,lsm) 996* 997 IF(IJBLK.GT.KLBLK) THEN 998 IJRELKL = 1 999 IBLOFF=IJKLOF(MAX(ISM,JSM),MIN(ISM,JSM),MAX(KSM,LSM)) 1000 ELSE IF (IJBLK.EQ.KLBLK) THEN 1001 IJRELKL = 0 1002 IBLOFF=IJKLOF(MAX(ISM,JSM),MIN(ISM,JSM),MAX(KSM,LSM)) 1003 ELSE IF (IJBLK.LT.KLBLK) THEN 1004 IJRELKL = -1 1005 IBLOFF = IJKLOF(MAX(KSM,LSM),MIN(KSM,LSM),MAX(ISM,JSM)) 1006 END IF 1007* 1008 itOrb=NTOOBS(iSm) 1009 jtOrb=NTOOBS(jSm) 1010 ktOrb=NTOOBS(kSm) 1011 ltOrb=NTOOBS(lSm) 1012* 1013 If(ISM.EQ.JSM) THEN 1014 IJPAIRS = ITORB*(ITORB+1)/2 1015 ELSE 1016 IJPAIRS = ITORB*JTORB 1017 END IF 1018* 1019 IF(KSM.EQ.LSM) THEN 1020 KLPAIRS = KTORB*(KTORB+1)/2 1021 ELSE 1022 KLPAIRS = KTORB*LTORB 1023 END IF 1024* 1025 iInt=0 1026 Do lJeppe=lOff,lOff+lOrb-1 1027 jMin=jOff 1028 If ( JLSM.ne.0 ) jMin=lJeppe 1029 Do jJeppe=jMin,jOff+jOrb-1 1030* 1031* 1032*. Set up array IJ*(IJ-1)/2 1033 IF(IJRELKL.EQ.0) THEN 1034 DO II = IOFF,IOFF+IORB-1 1035 IJ = I2INDX((JJEPPE-1)*NTOOB+II) 1036 IJARR(II) = IJ*(IJ-1)/2 1037 END DO 1038 END IF 1039* 1040 Do kJeppe=kOff,kOff+kOrb-1 1041 iMin = iOff 1042 kl = I2INDX(KJEPPE+(LJEPPE-1)*NTOOB) 1043 If(IKSM.ne.0) iMin = kJeppe 1044 IF(ICOUL.EQ.1) THEN 1045*. Address before integral (1,j!k,l) 1046 IINT = (LJEPPE-LOFF)*Jorb*Korb*Iorb 1047 & + (KJEPPE-KOFF)*Jorb*Iorb 1048 & + (JJEPPE-JOFF)*Iorb 1049 ELSE IF (ICOUL.EQ.2) THEN 1050* Address before (1L,JK) 1051 IINT = (KJEPPE-KOFF)*JORB*LORB*IORB 1052 & + (JJEPPE-JOFF) *LORB*IORB 1053 & + (LJEPPE-LOFF) *IORB 1054 END IF 1055* 1056 IF(IJRELKL.EQ.1) THEN 1057*. Block (ISM JSM ! KSM LSM ) with (Ism,jsm) > (ksm,lsm) 1058 IJKL0 = IBLOFF-1+(kl-1)*ijPairs 1059 IJ0 = (JJEPPE-1)*NTOOB 1060 Do iJeppe=iMin,iOff+iOrb-1 1061 ijkl = ijkl0 + I2INDX(IJEPPE+IJ0) 1062 iInt=iInt+1 1063 Xint(iInt) = Intlst(ijkl) 1064 End Do 1065 END IF 1066* 1067*. block (ISM JSM !ISM JSM) 1068 IF(IJRELKL.EQ.0) THEN 1069 IJ0 = (JJEPPE-1)*NTOOB 1070 KLOFF = KL*(KL-1)/2 1071 IJKL0 = (KL-1)*IJPAIRS-KLOFF 1072 Do iJeppe=iMin,iOff+iOrb-1 1073 ij = I2INDX(IJEPPE+IJ0 ) 1074 If ( ij.ge.kl ) Then 1075C ijkl=ij+(kl-1)*ijPairs-klOff 1076 IJKL = IJKL0 + IJ 1077 Else 1078 IJOFF = IJARR(IJEPPE) 1079 ijkl=kl+(ij-1)*klPairs-ijOff 1080 End If 1081 iInt=iInt+1 1082 Xint(iInt) = Intlst(iblOff-1+ijkl) 1083 End Do 1084 END IF 1085* 1086*. Block (ISM JSM ! KSM LSM ) with (Ism,jsm) < (ksm,lsm) 1087 IF(IJRELKL.EQ.-1) THEN 1088 ijkl0 = IBLOFF-1+KL - KLPAIRS 1089 IJ0 = (JJEPPE-1)*NTOOB 1090 Do iJeppe=iMin,iOff+iOrb-1 1091 IJKL = IJKL0 + I2INDX(IJEPPE + IJ0)*KLPAIRS 1092 iInt=iInt+1 1093 Xint(iInt) = Intlst(ijkl) 1094 End Do 1095 END IF 1096* 1097 End Do 1098 End Do 1099 End Do 1100* 1101* Collect Exchange terms 1102* 1103 If ( IXCHNG.ne.0 ) Then 1104* 1105 IF(ISM.EQ.LSM) THEN 1106 ILPAIRS = ITORB*(ITORB+1)/2 1107 ELSE 1108 ILPAIRS = ITORB*LTORB 1109 END IF 1110* 1111 IF(KSM.EQ.JSM) THEN 1112 KJPAIRS = KTORB*(KTORB+1)/2 1113 ELSE 1114 KJPAIRS = KTORB*JTORB 1115 END IF 1116* 1117 ilblk = max(ism,lsm)*(max(ism,lsm)-1)/2 + min(ism,lsm) 1118 kjblk = max(ksm,jsm)*(max(ksm,jsm)-1)/2 + min(ksm,jsm) 1119 IF(ILBLK.GT.KJBLK) THEN 1120 ILRELKJ = 1 1121 IBLOFF = IJKLOF(MAX(ISM,LSM),MIN(ISM,LSM),MAX(KSM,JSM)) 1122 ELSE IF(ILBLK.EQ.KJBLK) THEN 1123 ILRELKJ = 0 1124 IBLOFF = IJKLOF(MAX(ISM,LSM),MIN(ISM,LSM),MAX(KSM,JSM)) 1125 ELSE IF(ILBLK.LT.KJBLK) THEN 1126 ILRELKJ = -1 1127 IBLOFF = IJKLOF(MAX(KSM,JSM),MIN(KSM,JSM),MAX(ISM,LSM)) 1128 END IF 1129* 1130 iInt=0 1131 Do lJeppe=lOff,lOff+lOrb-1 1132 jMin=jOff 1133 If ( JLSM.ne.0 ) jMin=lJeppe 1134* 1135 IF(ILRELKJ.EQ.0) THEN 1136 DO II = IOFF,IOFF+IORB-1 1137 IL = I2INDX(II+(LJEPPE-1)*NTOOB) 1138 IJARR(II) = IL*(IL-1)/2 1139 END DO 1140 END IF 1141* 1142 Do jJeppe=jMin,jOff+jOrb-1 1143 Do kJeppe=kOff,kOff+kOrb-1 1144 KJ = I2INDX(KJEPPE+(JJEPPE-1)*NTOOB) 1145 KJOFF = KJ*(KJ-1)/2 1146 iMin = iOff 1147* 1148 IF(ICOUL.EQ.1) THEN 1149*. Address before integral (1,j!k,l) 1150 IINT = (LJEPPE-LOFF)*Jorb*Korb*Iorb 1151 & + (KJEPPE-KOFF)*Jorb*Iorb 1152 & + (JJEPPE-JOFF)*Iorb 1153 ELSE IF (ICOUL.EQ.2) THEN 1154* Address before (1L,JK) 1155 IINT = (KJEPPE-KOFF)*JORB*LORB*IORB 1156 & + (JJEPPE-JOFF) *LORB*IORB 1157 & + (LJEPPE-LOFF) *IORB 1158 END IF 1159* 1160 If(IKSM.ne.0) iMin = kJeppe 1161* 1162 IF(ILRELKJ.EQ.1) THEN 1163 ILKJ0 = IBLOFF-1+( kj-1)*ilpairs 1164 IL0 = (LJEPPE-1)*NTOOB 1165 Do iJeppe=iMin,iOff+iOrb-1 1166 ILKJ = ILKJ0 + I2INDX(IJEPPE + IL0) 1167 iInt=iInt+1 1168 XInt(iInt)=XInt(iInt)-Intlst(ilkj) 1169 End Do 1170 END IF 1171* 1172 IF(ILRELKJ.EQ.0) THEN 1173 IL0 = (LJEPPE-1)*NTOOB 1174 ILKJ0 = (kj-1)*ilPairs-kjOff 1175 Do iJeppe=iMin,iOff+iOrb-1 1176 IL = I2INDX(IJEPPE + IL0 ) 1177 If ( il.ge.kj ) Then 1178C ilkj=il+(kj-1)*ilPairs-kjOff 1179 ILKJ = IL + ILKJ0 1180 Else 1181 ILOFF = IJARR(IJEPPE) 1182 ilkj=kj+(il-1)*kjPairs-ilOff 1183 End If 1184 iInt=iInt+1 1185 XInt(iInt)=XInt(iInt)-Intlst(iBLoff-1+ilkj) 1186 End Do 1187 END IF 1188* 1189 IF(ILRELKJ.EQ.-1) THEN 1190 ILKJ0 = IBLOFF-1+KJ-KJPAIRS 1191 IL0 = (LJEPPE-1)*NTOOB 1192 Do iJeppe=iMin,iOff+iOrb-1 1193 ILKJ = ILKJ0 + I2INDX(IJEPPE+ IL0)*KJPAIRS 1194 iInt=iInt+1 1195 XInt(iInt)=XInt(iInt)-Intlst(ilkj) 1196 End Do 1197 END IF 1198* 1199 End Do 1200 End Do 1201 End Do 1202 End If 1203* 1204 RETURN 1205 END 1206*********************************************************************** 1207 SUBROUTINE LGETINT(XINT,ITP,ISM,JTP,JSM,KTP,KSM,LTP,LSM, 1208 & IXCHNG,IKSM,JLSM,ICOUL) 1209 1210* 1211* Outer routine for accessing integral block 1212* 1213 IMPLICIT REAL*8(A-H,O-Z) 1214 REAL*8 XINT(*) 1215* 1216#include "priunit.h" 1217* 1218* 1219#include "mxpdim.inc" 1220#include "lucinp.inc" 1221#include "orbinp.inc" 1222#include "csm.inc" 1223#include "cc_exc.inc" 1224#include "crun.inc" 1225#include "wrkspc.inc" 1226#include "glbbas.inc" 1227* 1228 CALL QENTER('LGETINT') 1229 NTEST = 00 1230* 1231 IF(NTEST.GE.5) 1232 &WRITE(lupri,*) ' LGETINT : ICC_EXC and ICOUL = ', ICC_EXC, ICOUL 1233 1234* ======================= 1235* Usual/Normal integrals 1236* ======================= 1237* 1238*. Integrals in core in internal LUCIA format 1239 CALL GETINCN2(XINT,ITP,ISM,JTP,JSM,KTP,KSM,LTP,LSM, 1240 & IXCHNG,IKSM,JLSM,WORK(KINT2), 1241 & WORK(KPINT2),NSMOB,WORK(KINH1),ICOUL) 1242 1243 IF(NTEST.NE.0) THEN 1244 IF(ITP.EQ.0) THEN 1245 NI = NTOOBS(ISM) 1246 ELSE 1247 NI = NOBPTS(ITP,ISM) 1248 END IF 1249* 1250 IF(KTP.EQ.0) THEN 1251 NK = NTOOBS(KSM) 1252 ELSE 1253 NK = NOBPTS(KTP,KSM) 1254 END IF 1255* 1256 IF(IKSM.EQ.0) THEN 1257 NIK = NI * NK 1258 ELSE 1259 NIK = NI*(NI+1)/2 1260 END IF 1261* 1262 IF(JTP.EQ.0) THEN 1263 NJ = NTOOBS(JSM) 1264 ELSE 1265 NJ = NOBPTS(JTP,JSM) 1266 END IF 1267* 1268 IF(LTP.EQ.0) THEN 1269 NL = NTOOBS(LSM) 1270 ELSE 1271 NL = NOBPTS(LTP,LSM) 1272 END IF 1273* 1274 IF(JLSM.EQ.0) THEN 1275 NJL = NJ * NL 1276 ELSE 1277 NJL = NJ*(NJ+1)/2 1278 END IF 1279 WRITE(lupri,*) ' 2 electron integral block for TS blocks ' 1280 WRITE(lupri,*) ' Ixchng :', IXCHNG 1281 WRITE(lupri,*) ' After GETINC ' 1282 WRITE(lupri,'(1H ,4(A,I2,A,I2,A))') 1283 & '(',ITP,',',ISM,')','(',JTP,',',JSM,')', 1284 & '(',KTP,',',KSM,')','(',LTP,',',LSM,')' 1285 CALL WRTMT_LU(XINT,NIK,NJL,NIK,NJL) 1286 END IF 1287* 1288 CALL QEXIT('LGETINT') 1289C Call Abend2( ' Jeppe forced me to stop in LGETINT ' ) 1290 RETURN 1291 END 1292*********************************************************************** 1293 1294 SUBROUTINE GETOBS_DALTON(ECORE_ENV,NAOS_ENV,NMOS_ENV,NSYM_ENV) 1295 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1296#include "priunit.h" 1297*. Scratch 1298 DIMENSION NBAS(8), NOCC(8), NORB(8) 1299*. Output 1300 DIMENSION NAOS_ENV(*), NMOS_ENV(*) 1301 1302C 1303C Read AO and MO information on file SIRIFC written from SIRIUS. 1304C 1305 ITAP30 = 16 1306 OPEN(ITAP30,STATUS='OLD',FORM='UNFORMATTED',FILE='SIRIFC') 1307 REWIND ITAP30 1308 CALL MOLLAB('TRCCINT ',ITAP30,lupri) 1309 READ (ITAP30) NSYM,NORBT,NBAST,NCMOT,(NOCC(I),I=1,NSYM), 1310 * (NORB(I),I=1,NSYM),(NBAS(I),I=1,NSYM), 1311 * POTNUC,EMCSCF 1312 1313! transfer to output variables 1314 NSYM_ENV = NSYM 1315 ECORE_ENV = POTNUC 1316 CALL ICOPY(NSYM,NORB,1,NMOS_ENV,1) 1317 CALL ICOPY(NSYM,NBAS,1,NAOS_ENV,1) 1318 1319#ifdef LUCI_DEBUG 1320 1321 write(lupri,*) 'GETOBS_DALTON' !hjaaj DEBUG 1322 WRITE(lupri,*) ' Number of basis functions per sym ' 1323 CALL IWRTMA(NBAS,1,8,1,8) 1324* 1325 WRITE(lupri,*) ' Norb as delivered from environment ' 1326 CALL IWRTMA(NORB,1,8,1,8) 1327* 1328 WRITE(lupri,*) ' NOCC as delivered from DALTON (discarded)' 1329 CALL IWRTMA(NOCC,1,8,1,8) 1330*. 1331 WRITE(lupri,*) ' from DALTON: NORBT, NBAST, NCMOT :', 1332 & NORBT,NBAST,NCMOT 1333 1334#endif 1335 1336 END 1337*********************************************************************** 1338 SUBROUTINE GETOBS_LUCIA(NAOS_ENV,NMOS_ENV) 1339* 1340* Obtain info on orbital dimensions from LU91 - LUCIA format 1341* 1342* Jeppe Olsen, Feb. 98 1343* 1344 IMPLICIT REAL*8(A-H,O-Z) 1345*. Output 1346 INTEGER NMOS_ENV(*),NAOS_ENV(*) 1347* 1348 LUH = 91 1349 CALL REWINO(LUH) 1350*. 1351 READ(LUH,*) NSMOB 1352*. 1353 READ(LUH,*) (NMOS_ENV(ISM),ISM=1, NSMOB) 1354* 1355 READ(LUH,*) (NAOS_ENV(ISM),ISM=1, NSMOB) 1356* 1357 RETURN 1358 END 1359*********************************************************************** 1360 FUNCTION GTIJKL(I,J,K,L) 1361* 1362* Obtain integral (I J ! K L ) 1363* where I,J,K and l refers to active orbitals in 1364* Type ordering 1365* 1366 IMPLICIT REAL*8(A-H,O-Z) 1367#include "priunit.h" 1368#include "mxpdim.inc" 1369C COMMON/BIGGY/WORK(MXPWRD) 1370#include "wrkspc.inc" 1371*.GLobal pointers 1372C COMMON/GLBBAS/KINT1,KINT2,KPINT1,KPINT2,KLSM1,KLSM2,KRHO1 1373#include "glbbas.inc" 1374#include "lucinp.inc" 1375#include "orbinp.inc" 1376#include "crun.inc" 1377 IF(INTIMP .EQ. 2 ) THEN 1378*. LUCAS ordering 1379 I12S = 0 1380 I34S = 0 1381 I1234S = 1 1382 GTIJKL = GIJKLL(IREOTS(1),WORK(KPINT2),WORK(KLSM2), 1383 & WORK(KINT2),ISMFTO,IBSO,NACOB,NSMOB, 1384 & NOCOBS,I,J,K,L) 1385 ELSE IF (INTIMP.EQ.1.OR.INTIMP.EQ.5.or.INTIMP.eq.6) THEN 1386! sirius or dirac integral format 1387 GTIJKL = GMIJKL(I,J,K,L,WORK(KINT2),WORK(KPINT2)) 1388 END IF 1389 1390 END 1391 1392*********************************************************************** 1393 FUNCTION GMIJKL(IORB,JORB,KORB,LORB,INTLST,IJKLOF) 1394* 1395* Obtain integral (IORB JORB ! KORB LORB) MOLCAS version 1396* Integrals assumed in core 1397* 1398* Version for integrals stored in INTLST 1399* 1400 IMPLICIT REAL*8(A-H,O-Z) 1401#include "priunit.h" 1402#include "mxpdim.inc" 1403#include "orbinp.inc" 1404#include "lucinp.inc" 1405*. Integral list 1406 Real * 8 Intlst(*) 1407 Dimension IJKLOF(NsmOB,NsmOb,NsmOB) 1408 Logical iSymj,kSyml,ISYMK,JSYML,ijSymkl,IKSYMJL 1409 Logical ijklPerm 1410*. 1411 NTEST = 000 1412* 1413*. The orbital list corresponds to type ordered indices, reform to 1414*. symmetry ordering 1415* 1416 IABS = IREOTS(IORB) 1417 ISM = ISMFTO(IORB) 1418 IOFF = IBSO(ISM) 1419* 1420 JABS = IREOTS(JORB) 1421 JSM = ISMFTO(JORB) 1422 JOFF = IBSO(JSM) 1423* 1424 KABS = IREOTS(KORB) 1425 KSM = ISMFTO(KORB) 1426 KOFF = IBSO(KSM) 1427* 1428 LABS = IREOTS(LORB) 1429 LSM = ISMFTO(LORB) 1430 LOFF = IBSO(LSM) 1431* 1432 If( Ntest.ge. 100) THEN 1433 write(lupri,*) ' GMIJKL at your service ' 1434 WRITE(lupri,*) ' IORB IABS ISM IOFF ',IORB,IABS,ISM,IOFF 1435 WRITE(lupri,*) ' JORB JABS JSM JOFF ',JORB,JABS,JSM,JOFF 1436 WRITE(lupri,*) ' KORB KABS KSM KOFF ',KORB,KABS,KSM,KOFF 1437 WRITE(lupri,*) ' LORB LABS LSM LOFF ',LORB,LABS,LSM,LOFF 1438 call flshfo(lupri) 1439 END IF 1440* 1441 If ( jSm.gt.iSm .or. ( iSm.eq.jSm .and. JABS.gt.IABS)) Then 1442 iSym=jSm 1443 jSym=iSm 1444 I = JABS - JOFF + 1 1445 J = IABS - IOFF + 1 1446 Else 1447 iSym=iSm 1448 jSym=jSm 1449 I = IABS - IOFF + 1 1450 J = JABS - JOFF + 1 1451 End If 1452 ijBlk=jSym+iSym*(iSym-1)/2 1453 If ( lSm.gt.kSm .or. ( kSm.eq.lSm .and. LABS.gt.KABS)) Then 1454 kSym=lSm 1455 lSym=kSm 1456 K = LABS -LOFF + 1 1457 L = KABS - KOFF + 1 1458 Else 1459 kSym=kSm 1460 lSym=lSm 1461 K = KABS - KOFF + 1 1462 L = LABS -LOFF + 1 1463 End If 1464 klBlk=lSym+kSym*(kSym-1)/2 1465* 1466 ijklPerm=.false. 1467 If ( klBlk.gt.ijBlk ) Then 1468 iTemp=iSym 1469 iSym=kSym 1470 kSym=iTemp 1471 iTemp=jSym 1472 jSym=lSym 1473 lSym=iTemp 1474 iTemp=ijBlk 1475 ijBlk=klBlk 1476 klBlk=iTemp 1477 ijklPerm=.true. 1478* 1479 iTemp = i 1480 i = k 1481 k = itemp 1482 iTemp = j 1483 j = l 1484 l = iTemp 1485 End If 1486 If(Ntest .ge. 100 ) then 1487 write(lupri,*) ' i j k l ',i,j,k,l 1488 write(lupri,*) ' Isym,Jsym,Ksym,Lsym',Isym,Jsym,Ksym,Lsym 1489 call flshfo(lupri) 1490 End if 1491* 1492* Define offset for given symmetry block 1493 IBLoff = IJKLof(Isym,Jsym,Ksym) 1494 If(ntest .ge. 100 ) 1495 &WRITE(lupri,*) ' IBLoff Isym Jsym Ksym ', IBLoff,ISym,Jsym,Ksym 1496 iSymj=iSym.eq.jSym 1497 kSyml=kSym.eq.lSym 1498 iSymk=iSym.eq.kSym 1499 jSyml=jSym.eq.lSym 1500 ikSymjl=iSymk.and.jSyml 1501 ijSymkl=iSymj.and.kSyml 1502* 1503 itOrb=NTOOBS(iSym) 1504 jtOrb=NTOOBS(jSym) 1505 ktOrb=NTOOBS(kSym) 1506 ltOrb=NTOOBS(lSym) 1507C? print *,' itOrb,jtOrb,ktOrb,ltOrb',itOrb,jtOrb,ktOrb,ltOrb 1508 If ( iSymj ) Then 1509 ijPairs=itOrb*(itOrb+1)/2 1510 ij=j+i*(i-1)/2 1511 Else 1512 ijPairs=itOrb*jtOrb 1513 ij=j + (i-1)*jtOrb 1514 End if 1515* 1516 IF(KSYML ) THEN 1517 klPairs=ktOrb*(ktOrb+1)/2 1518 kl=l+k*(k-1)/2 1519 ELSE 1520 klPairs=ktOrb*ltOrb 1521 kl=l+(k-1)*ltOrb 1522 End If 1523C? print *,' ijPairs,klPairs',ijPairs,klPairs 1524* 1525 If ( ikSymjl ) Then 1526 If ( ij.gt.kl ) Then 1527 klOff=kl+(kl-1)*(kl-2)/2-1 1528 ijkl=ij+(kl-1)*ijPairs-klOff 1529 Else 1530 ijOff=ij+(ij-1)*(ij-2)/2-1 1531 ijkl=kl+(ij-1)*klPairs-ijOff 1532 End If 1533 Else 1534 ijkl=ij+(kl-1)*ijPairs 1535 End If 1536 If( ntest .ge. 100 )then 1537 write(lupri,*) ' ijkl ', ijkl 1538 write(lupri,*) ' iblOff ', iblOff 1539 write(lupri,*) ' iblOff-1+ijkl', iblOff-1+ijkl 1540 call flshfo(lupri) 1541 end if 1542* 1543 GMIJKL = Intlst(iblOff-1+ijkl) 1544 If( ntest .ge. 100 ) 1545 & write(lupri,*) ' GMIJKL ', GMIJKL 1546* 1547 END 1548*********************************************************************** 1549 SUBROUTINE GTJK(RJ,RK,NTOOB) 1550* 1551* Interface routine for obtaining Coulomb (RJ) and 1552* Exchange integrals (RK) 1553* 1554* Ordering of intgrals is in the internal order 1555 IMPLICIT REAL*8(A-H,O-Z) 1556* 1557*.CRUN 1558C COMMON/CRUN/MAXIT,IRESTR,INTIMP,NP1,NP2,NQ,INCORE,MXCIV,ICISTR, 1559C & NOCSF,IDIAG 1560#include "mxpdim.inc" 1561#include "parluci.h" 1562#include "crun.inc" 1563*.Output 1564 DIMENSION RJ(NTOOB,NTOOB),RK(NTOOB,NTOOB) 1565 1566 IF(INTIMP.EQ.1.OR.INTIMP.EQ.5.or.INTIMP.eq.6) THEN 1567*. Interface to SIRIUS 1568 CALL GTJKS(RJ,RK,NTOOB) 1569 ELSE 1570*. Interface to LUCAS integrals 1571 CALL GTJKL(RJ,RK,NTOOB) 1572 END IF 1573* 1574 NTEST = 0 1575 IF(NTEST.NE.0) THEN 1576 WRITE(luwrt,*) ' RJ and RK from GTJK ' 1577 CALL WRTMT_LU(RJ,NTOOB,NTOOB,NTOOB,NTOOB) 1578 CALL WRTMT_LU(RK,NTOOB,NTOOB,NTOOB,NTOOB) 1579 END IF 1580* 1581 RETURN 1582 END 1583*********************************************************************** 1584 SUBROUTINE GTJKL(RJ,RK,NTOOB) 1585* 1586* Obtain Coulomb integrals (II!JJ) 1587* exchange integrals (IJ!JI) 1588* 1589 IMPLICIT REAL*8(A-H,O-Z) 1590#include "priunit.h" 1591 DIMENSION RJ(NTOOB,NTOOB),RK(NTOOB,NTOOB) 1592* 1593 DO 100 IORB = 1, NTOOB 1594 DO 50 JORB = 1, NTOOB 1595 RJ(IORB,JORB) = GTIJKL(IORB,IORB,JORB,JORB) 1596 RK(IORB,JORB) = GTIJKL(IORB,JORB,JORB,IORB) 1597 50 CONTINUE 1598 100 CONTINUE 1599* 1600 NTEST = 0 1601 IF(NTEST.NE.0) THEN 1602 WRITE(lupri,*) ' RJ and RK from GTJK ' 1603 CALL WRTMT_LU(RJ,NTOOB,NTOOB,NTOOB,NTOOB) 1604 CALL WRTMT_LU(RK,NTOOB,NTOOB,NTOOB,NTOOB) 1605 END IF 1606* 1607 RETURN 1608 END 1609*********************************************************************** 1610 1611* Working on EXPHAM 1612* some known problems : 1613* 1 : if CSF are used diagonal is not delivered to H0mat 1614 SUBROUTINE GTJKS(J,K,NORB) 1615* 1616* Obtain Coulomb and Exchange integrals 1617* from complete integral list stored in core 1618* 1619 IMPLICIT REAL*8 (A-H,O-Z) 1620#include "priunit.h" 1621 REAL*8 J(NORB,NORB),K(NORB,NORB) 1622* 1623 DO 200 IORB = 1, NORB 1624 DO 100 JORB = 1, NORB 1625! write(lupri,*) 'iorb,jorb ==> j',iorb,jorb 1626 J(IORB,JORB) = GTIJKL(IORB,IORB,JORB,JORB) 1627! write(lupri,*) 'iorb,jorb ==> k',iorb,jorb 1628 K(IORB,JORB) = GTIJKL(IORB,JORB,JORB,IORB) 1629 100 CONTINUE 1630 200 CONTINUE 1631* 1632 END 1633*********************************************************************** 1634 subroutine hello_dalton_lucita 1635************************************************************************ 1636* * 1637* Print the program banner, date and time of execution * 1638* * 1639*----------------------------------------------------------------------* 1640* * 1641* written by: * 1642* M.P. Fuelscher * 1643* University of Lund, Sweden, 1993 * 1644* Modified, Timo Fleig, Dec 2001, for DIRAC * 1645* Aug 2004 * 1646* Aug 2006 * 1647* HJAaJ, May 2008, for DALTON * 1648* * 1649************************************************************************ 1650#include "priunit.h" 1651 Character*8 Fmt 1652 Character*70 Line,StLine 1653 integer :: lpaper = 120 1654*----------------------------------------------------------------------* 1655* Start and define the paper width * 1656* Initialize blank and header lines * 1657*----------------------------------------------------------------------* 1658 lLine=Len(Line) 1659 Do i=1,lLine 1660 StLine(i:i)='*' 1661 End Do 1662 left=(lPaper-lLine)/2 1663 Write(Fmt,'(A,I3.3,A)') '(',left,'X,A)' 1664*----------------------------------------------------------------------* 1665* Print the program header * 1666*----------------------------------------------------------------------* 1667 write(lupri,'(/1x,a)') StLine 1668 write(lupri,'( 1x,a)') StLine,' ', 1669 & 'D A L T O N - L U C I T A', 1670 & 'An interface section for LUCITA under DALTON', ' ', 1671 & 'Authors: J. Olsen, Univ. Aarhus', 1672 & ' H. J. Aa. Jensen, Univ. Southern Denmark ', 1673 & ' S. Knecht, Univ. Southern Denmark',' ', 1674! & 'Based on LUCITA-DIRAC interface', 1675! & ' Author: Timo Fleig, Univ. Toulouse', ' ', 1676 & 'Using LUCIA version 1999', ' ', 1677 & ' Author: J. Olsen, Lund/Aarhus',' ', 1678 & 'Parallelization of LUCITA, Duesseldorf/Odense: ', 1679 & ' S. Knecht, Univ. Southern Denmark', ' ', 1680 & 'Citations:', 1681 & ' J. Olsen, P. Joergensen, J. Simons, ', 1682 & ' Chem. Phys. Lett. 169 (1990) 463 ', 1683 & ' S. Knecht, H. J. Aa. Jensen and T. Fleig, ', 1684 & ' J. Chem. Phys., 128 (2008) 014108 ',' ', 1685 & StLine, Stline, Stline 1686 1687 end 1688*********************************************************************** 1689 Function I2EAD(IORB,JORB,KORB,LORB) 1690* 1691* Find adress of integral in LUCIA order 1692* 1693 IMPLICIT REAL*8 (A-H,O-Z) 1694* 1695#include "mxpdim.inc" 1696 1697#include "glbbas.inc" 1698* 1699#include "wrkspc.inc" 1700* 1701 I2EAD = I2EADS(IORB,JORB,KORB,LORB,WORK(KPINT2)) 1702* 1703 RETURN 1704 END 1705*********************************************************************** 1706 FUNCTION I2EADS(IORB,JORB,KORB,LORB,IJKLOF) 1707* 1708* Obtain address of integral (IORB JORB ! KORB LORB) in MOLCAS order 1709* IORB JORB KORB LORB corresponds to SYMMETRY ordered indices !! 1710* Integrals assumed in core 1711* 1712* 1713 IMPLICIT REAL*8(A-H,O-Z) 1714* 1715#include "priunit.h" 1716* 1717* 1718#include "mxpdim.inc" 1719#include "orbinp.inc" 1720#include "lucinp.inc" 1721* 1722 Dimension IJKLOF(NsmOB,NsmOb,NsmOB) 1723 Logical iSymj,kSyml,ISYMK,JSYML,ijSymkl,IKSYMJL 1724 Logical ijklPerm 1725*. 1726 NTEST = 000 1727* 1728 IABS = IORB 1729 ISM = ISMFTO(IREOST(IORB)) 1730 IOFF = IBSO(ISM) 1731* 1732 JABS = JORB 1733 JSM = ISMFTO(IREOST(JORB)) 1734 JOFF = IBSO(JSM) 1735* 1736 KABS = KORB 1737 KSM = ISMFTO(IREOST(KORB)) 1738 KOFF = IBSO(KSM) 1739* 1740 LABS = LORB 1741 LSM = ISMFTO(IREOST(LORB)) 1742 LOFF = IBSO(LSM) 1743* 1744 If( Ntest.ge. 100) THEN 1745 write(lupri,*) ' I2EADS at your service ' 1746 WRITE(lupri,*) ' IORB IABS ISM IOFF ',IORB,IABS,ISM,IOFF 1747 WRITE(lupri,*) ' JORB JABS JSM JOFF ',JORB,JABS,JSM,JOFF 1748 WRITE(lupri,*) ' KORB KABS KSM KOFF ',KORB,KABS,KSM,KOFF 1749 WRITE(lupri,*) ' LORB LABS LSM LOFF ',LORB,LABS,LSM,LOFF 1750 END IF 1751* 1752 If ( jSm.gt.iSm .or. ( iSm.eq.jSm .and. JABS.gt.IABS)) Then 1753 iSym=jSm 1754 jSym=iSm 1755 I = JABS - JOFF + 1 1756 J = IABS - IOFF + 1 1757 Else 1758 iSym=iSm 1759 jSym=jSm 1760 I = IABS - IOFF + 1 1761 J = JABS - JOFF + 1 1762 End If 1763 ijBlk=jSym+iSym*(iSym-1)/2 1764 If ( lSm.gt.kSm .or. ( kSm.eq.lSm .and. LABS.gt.KABS)) Then 1765 kSym=lSm 1766 lSym=kSm 1767 K = LABS -LOFF + 1 1768 L = KABS - KOFF + 1 1769 Else 1770 kSym=kSm 1771 lSym=lSm 1772 K = KABS - KOFF + 1 1773 L = LABS -LOFF + 1 1774 End If 1775 klBlk=lSym+kSym*(kSym-1)/2 1776* 1777 ijklPerm=.false. 1778 If ( klBlk.gt.ijBlk ) Then 1779 iTemp=iSym 1780 iSym=kSym 1781 kSym=iTemp 1782 iTemp=jSym 1783 jSym=lSym 1784 lSym=iTemp 1785 iTemp=ijBlk 1786 ijBlk=klBlk 1787 klBlk=iTemp 1788 ijklPerm=.true. 1789* 1790 iTemp = i 1791 i = k 1792 k = itemp 1793 iTemp = j 1794 j = l 1795 l = iTemp 1796 End If 1797 If(Ntest .ge. 100 ) then 1798 write(lupri,*) ' i j k l ',i,j,k,l 1799 write(lupri,*) ' Isym,Jsym,Ksym,Lsym',Isym,Jsym,Ksym,Lsym 1800 End if 1801* 1802* Define offset for given symmetry block 1803 IBLoff = IJKLof(Isym,Jsym,Ksym) 1804 If(ntest .ge. 100 ) 1805 &WRITE(lupri,*) ' IBLoff Isym Jsym Ksym ', IBLoff,ISym,Jsym,Ksym 1806 iSymj=iSym.eq.jSym 1807 kSyml=kSym.eq.lSym 1808 iSymk=iSym.eq.kSym 1809 jSyml=jSym.eq.lSym 1810 ikSymjl=iSymk.and.jSyml 1811 ijSymkl=iSymj.and.kSyml 1812* 1813 itOrb=NTOOBS(iSym) 1814 jtOrb=NTOOBS(jSym) 1815 ktOrb=NTOOBS(kSym) 1816 ltOrb=NTOOBS(lSym) 1817C? print *,' itOrb,jtOrb,ktOrb,ltOrb',itOrb,jtOrb,ktOrb,ltOrb 1818 If ( iSymj ) Then 1819 ijPairs=itOrb*(itOrb+1)/2 1820 ij=j+i*(i-1)/2 1821 Else 1822 ijPairs=itOrb*jtOrb 1823 ij=j + (i-1)*jtOrb 1824 End if 1825* 1826 IF(KSYML ) THEN 1827 klPairs=ktOrb*(ktOrb+1)/2 1828 kl=l+k*(k-1)/2 1829 ELSE 1830 klPairs=ktOrb*ltOrb 1831 kl=l+(k-1)*ltOrb 1832 End If 1833C? print *,' ijPairs,klPairs',ijPairs,klPairs 1834* 1835 If ( ikSymjl ) Then 1836 If ( ij.gt.kl ) Then 1837 klOff=kl+(kl-1)*(kl-2)/2-1 1838 ijkl=ij+(kl-1)*ijPairs-klOff 1839 Else 1840 ijOff=ij+(ij-1)*(ij-2)/2-1 1841 ijkl=kl+(ij-1)*klPairs-ijOff 1842 End If 1843 Else 1844 ijkl=ij+(kl-1)*ijPairs 1845 End If 1846 If( ntest .ge. 100 ) 1847 & write(lupri,*) ' ijkl ', ijkl 1848* 1849 I2EADS = iblOff-1+ijkl 1850 If( ntest .ge. 100 ) then 1851 write(lupri,*) 'i j k l ', i,j,k,l 1852 write(lupri,*) ' ibloff ijkl ',ibloff,ijkl 1853 write(lupri,*) ' I2EADS = ', I2EADS 1854 END IF 1855* 1856 END 1857 1858*********************************************************************** 1859 subroutine intim(citask_id,int1_mcscf,int2_mcscf, 1860 & update_ijkl_from_env, 1861 & orbital_trial_vector, 1862 & ci_trial_vector, 1863 & mcscf_provides_integrals) 1864* 1865* Interface to external integrals 1866* 1867* If NOINT .ne. 0, only pointers are constructed 1868* Jeppe Olsen, Winter of 1991 1869* 1870* Version : Fall 97 1871* 1872* 1873#ifdef VAR_MPI 1874 use sync_coworkers 1875#endif 1876 use lucita_cfg 1877 use lucita_energy_types 1878 use lucita_mcscf_ci_cfg 1879#include "implicit.h" 1880 ! for addressing of WORK 1881 INTEGER*8 current_free_mem, kintdal_interface 1882 character (len=12), intent(in) :: citask_id 1883 real(8), intent(inout) :: int1_mcscf(*) 1884 real(8), intent(inout) :: int2_mcscf(*) 1885 logical, intent(inout) :: update_ijkl_from_env 1886 logical, intent(in) :: orbital_trial_vector 1887 logical, intent(in) :: ci_trial_vector 1888 logical, intent(in) :: mcscf_provides_integrals 1889#include "priunit.h" 1890#include "mxpdim.inc" 1891#include "wrkspc.inc" 1892#include "crun.inc" 1893#include "glbbas.inc" 1894#include "clunit.inc" 1895#include "lucinp.inc" 1896#include "csm.inc" 1897#include "orbinp.inc" 1898#include "parluci.h" 1899#include "oper.inc" 1900*./CINTFO/ 1901 COMMON/CINTFO/I12S,I34S,I1234S,NINT1,NINT2,NBINT1,NBINT2 1902 character (len=25) :: file_name 1903 logical :: update_ijkl, integral_available 1904 integer :: ngsh_lucita_tmp(max_number_of_gas_spaces, 1905 & max_number_of_ptg_irreps) 1906 integral_available = .false. 1907 1908! pointers for symmetry blocks of integrals 1909 CALL INTPNT(WORK(KPINT1),WORK(KLSM1), 1910 & WORK(KPINT2),WORK(KLSM2)) 1911 1912! pointer for orbital indices for symmetry blocked matrices 1913 CALL ORBINH1(WORK(KINH1),NTOOBS,NTOOB,NSMOB) 1914 1915! 1916! load one-electron integrals and two-electron integrals 1917 if(citask_id.eq.'return CIdia' .or. 1918 & citask_id.eq.'perform Davi' .or. 1919 & citask_id.eq.'ijkl resort ' .or. 1920 & citask_id.eq.'return sigma' )then 1921 1922 IF(INCORE.EQ.1)THEN 1923 write(file_name,'(a17,i3,a1,i4)') 1924 & 'IJKL_REOD_LUCITA-',1,'.',luci_myproc 1925 do i = 18, 20 1926 if(file_name(i:i).eq. ' ') file_name(i:i) = '0' 1927 end do 1928 do i = 22, 25 1929 if(file_name(i:i).eq. ' ') file_name(i:i) = '0' 1930 end do 1931 LUINT_LUCITA = -1 1932! 1933! check for existing file with reordered integrals - 1934! if available read from it... 1935 if(citask_id.eq.'perform Davi' .or. 1936 & citask_id.eq.'ijkl resort ')then 1937 if(luci_myproc .eq. luci_master)then 1938 inquire(file=file_name,exist=integral_available) 1939 end if 1940 end if 1941! 1942 CALL GPOPEN(LUINT_LUCITA,file_name,'UNKNOWN',' ', 1943 & 'UNFORMATTED',IDUMMY,.FALSE.) 1944 is_2e_int_active = i12 1945 1946 update_ijkl = 1947 & update_ijkl_from_env.or.ci_trial_vector.or. 1948 & orbital_trial_vector 1949 1950! default for reordering task 1951 if(citask_id.eq.'ijkl resort ') update_ijkl = .true. 1952 1953 if(update_ijkl)then 1954 1955 if(luci_myproc .eq. luci_master)then 1956 if(.not.integral_available)then 1957 CALL MEMMAN(current_free_mem,0,'SFREEM',2,'SHOWMM') 1958 LSCR = current_free_mem - 1000 1959 idum = 0 1960! set local marker + allocate space for lucita-dalton integral interface 1961 call memman(idum,idum,'MARK ',idum,'dalint') 1962 call memman(kintdal_interface,lscr,'ADDL ',2,'dalint') 1963 1964 CALL LUCITA_GETINT_DALTON(einact_mc2lu,WORK(KINT1), 1965 & WORK(KINT2), 1966 & WORK(kintdal_interface),LSCR, 1967 & int1_mcscf,int2_mcscf, 1968 & mcscf_provides_integrals,1) 1969 ecore = ecore_orig + einact_mc2lu 1970! release memory 1971 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'dalint') 1972 else 1973! read 1e- and 2e-integrals from file 1974! ----------------------------------- 1975 CALL MOLLAB('LUCIIJKL',LUINT_LUCITA,LUPRI) 1976 READ (LUINT_LUCITA) nint1, nint2, ecore, ecore_orig , 1977 & einact_mc2lu 1978 read(LUINT_LUCITA) 1979 & ((ngsh_lucita_tmp(i,j),i=1,max_number_of_gas_spaces), 1980 & j=1,max_number_of_ptg_irreps) 1981 CALL READT(LUINT_LUCITA, nint1, work(kint1)) 1982 if(is_2e_int_active .gt. 1) 1983 & CALL READT(LUINT_LUCITA, nint2, work(kint2)) 1984! check consistency in active spaces for integral run 1985! and actual run 1986 do i = 1, max_number_of_gas_spaces 1987 do j = 1, max_number_of_ptg_irreps 1988 if(ngsh_lucita_tmp(i,j) /= ngsh_lucita(i,j))then 1989 write(lupri,*) '*** error in intim: active spaces'// 1990 & ' in integral run and present CI run do not match.' 1991 write(lupri,*) 'integral run ',ngsh_lucita_tmp(:,:) 1992 write(lupri,*) 'actual CI ',ngsh_lucita(:,:) 1993 call quit('inconsistency in active spaces for the'// 1994 & ' integral CI-run and the present CI-run.') 1995 end if 1996 end do 1997 end do 1998 end if 1999!#define LUCI_DEBUG 2000#ifdef LUCI_DEBUG 2001 write(lupri,*) 2002 & '*** master ',luci_myproc,'reports 1-ints:', 2003 & nint1 2004 call wrtmatmn(work(kint1),1,nint1,1,nint1,lupri) 2005 write(lupri,*) 2006 & '*** master ',luci_myproc,'reports 2-ints:', 2007 & nint2 2008 call wrtmatmn(work(kint2),1,nint2,1,nint2,lupri) 2009 write(lupri,*) 2010 & '*** master ',luci_myproc,'reports ecore ', 2011 & ecore, ecore_orig, einact_mc2lu 2012#undef LUCI_DEBUG 2013#endif 2014 end if ! luci_master only 2015#ifdef VAR_MPI 2016 if(luci_nmproc .gt. 1)then 2017! set sync_ctrl option 2018 sync_ctrl_array(3) = update_ijkl 2019! synchronize the co-workers with the 1-/2-electron integrals 2020 call sync_coworkers_cfg(work(kint1),work(kint2)) 2021 2022#ifdef LUCI_DEBUG 2023 if(luci_myproc .eq. luci_master+1)then 2024 print *, '*** co-worker',luci_myproc,'reports 1-ints:', 2025 & nint1 2026 call wrtmatmn(work(kint1),1,nint1,1,nint1,lupri) 2027 print *, '*** co-worker',luci_myproc,'reports 2-ints:', 2028 & nint2 2029 call wrtmatmn(work(kint2),1,nint2,1,nint2,lupri) 2030 end if ! debug print for co-worker-id luci_master+1 2031!#undef LUCI_DEBUG 2032#endif 2033 end if 2034#endif 2035 2036 if(.not.orbital_trial_vector.and..not.ci_trial_vector)then 2037! put 1e- and 2e-integrals to file 2038! -------------------------------- 2039 WRITE (LUINT_LUCITA) ('********',I=1,3), 'LUCIIJKL' 2040 WRITE(LUINT_LUCITA) nint1, nint2, ecore, ecore_orig , 2041 & einact_mc2lu 2042 write(LUINT_LUCITA) 2043 & ((ngsh_lucita(i,j),i=1,max_number_of_gas_spaces), 2044 & j=1,max_number_of_ptg_irreps) 2045 CALL WRITT(LUINT_LUCITA, nint1, work(kint1)) 2046 if(is_2e_int_active .gt. 1) 2047 & CALL WRITT(LUINT_LUCITA, nint2, work(kint2)) 2048 WRITE (LUINT_LUCITA) ('********',I=1,3), 'EODATA ' 2049 update_ijkl_from_env = .false. 2050 end if 2051! 1e-integrals from MCSCF environment, 2e-integrals from file 2052 if(ci_trial_vector)then 2053 if(is_2e_int_active .gt. 1)then 2054 CALL MOLLAB('LUCIIJKL',LUINT_LUCITA,LUPRI) 2055 READ (LUINT_LUCITA) nint1, nint2, ecore, ecore_orig , 2056 & einact_mc2lu 2057 read(LUINT_LUCITA) 2058 & ((ngsh_lucita_tmp(i,j),i=1,max_number_of_gas_spaces), 2059 & j=1,max_number_of_ptg_irreps) 2060 CALL READT(LUINT_LUCITA, nint1, work(kint2)) 2061 CALL READT(LUINT_LUCITA, nint2, work(kint2)) 2062 end if 2063 end if 2064 else 2065! read 1e- and 2e-integrals from file 2066! ----------------------------------- 2067 CALL MOLLAB('LUCIIJKL',LUINT_LUCITA,LUPRI) 2068 READ (LUINT_LUCITA) nint1, nint2, ecore, ecore_orig , 2069 & einact_mc2lu 2070 read(LUINT_LUCITA) 2071 & ((ngsh_lucita_tmp(i,j),i=1,max_number_of_gas_spaces), 2072 & j=1,max_number_of_ptg_irreps) 2073 CALL READT(LUINT_LUCITA, nint1, work(kint1)) 2074 if(is_2e_int_active .gt. 1) 2075 & CALL READT(LUINT_LUCITA, nint2, work(kint2)) 2076 end if 2077 CALL GPCLOSE(LUINT_LUCITA, 'KEEP') 2078 ELSE 2079 CALL QUIT('LUCITA out-of-core not implemented yet, sorry!') 2080 END IF 2081 2082 call dcopy(nint1,work(kint1),1,work(kint1o),1) 2083 2084 IF(IUSE_PH.EQ.1) CALL FI(WORK(KINT1),ECORE_HEX,1) 2085 2086 ECORE_ORIG = ECORE 2087 ECORE = ECORE + ECORE_HEX 2088 2089#ifdef LUCI_DEBUG 2090 if(luci_myproc.eq.luci_master)then 2091 write(lupri,'(/2x,a,e15.8)') 'Updated core energy: ',ECORE 2092 end if 2093#endif 2094 else 2095 2096 end if ! citask_id switch 2097 END 2098*********************************************************************** 2099 2100 SUBROUTINE PUTINT(XINT,ITP,ISM,JTP,JSM,KTP,KSM,LTP,LSM) 2101* 2102* Put integrals in permanent integral list 2103* 2104* Jeppe Olsen, Jan. 1999 2105* 2106 IMPLICIT REAL*8(A-H,O-Z) 2107* 2108#include "mxpdim.inc" 2109#include "lucinp.inc" 2110#include "orbinp.inc" 2111#include "wrkspc.inc" 2112#include "glbbas.inc" 2113*. Specific input 2114 DIMENSION XINT(*) 2115* 2116 CALL QENTER('PUTIN') 2117*. Offset and number of integrals 2118* 2119 IF(ITP.EQ.0) THEN 2120 NI = NTOOBS(ISM) 2121 ELSE 2122 NI = NOBPTS(ITP,ISM) 2123 END IF 2124* 2125 IOFF = IBSO(ISM) 2126 DO IITP = 1, ITP -1 2127 IOFF = IOFF + NOBPTS(IITP,ISM) 2128 END DO 2129* 2130 IF(JTP.EQ.0) THEN 2131 NJ = NTOOBS(JSM) 2132 ELSE 2133 NJ = NOBPTS(JTP,JSM) 2134 END IF 2135* 2136 JOFF = IBSO(JSM) 2137 DO JJTP = 1, JTP -1 2138 JOFF = JOFF + NOBPTS(JJTP,JSM) 2139 END DO 2140* 2141 IF(KTP.EQ.0) THEN 2142 NK = NTOOBS(KSM) 2143 ELSE 2144 NK = NOBPTS(KTP,KSM) 2145 END IF 2146* 2147 KOFF = IBSO(KSM) 2148 DO KKTP = 1, KTP -1 2149 KOFF = KOFF + NOBPTS(KKTP,KSM) 2150 END DO 2151* 2152 IF(LTP.EQ.0) THEN 2153 NL = NTOOBS(LSM) 2154 ELSE 2155 NL = NOBPTS(LTP,LSM) 2156 END IF 2157* 2158 LOFF = IBSO(LSM) 2159 DO LLTP = 1, LTP -1 2160 LOFF = LOFF + NOBPTS(LLTP,LSM) 2161 END DO 2162* 2163 INT_IN = 0 2164 DO LOB = LOFF,LOFF+NL-1 2165 DO KOB = KOFF,KOFF+NK-1 2166 DO JOB = JOFF,JOFF+NJ-1 2167 DO IOB = IOFF,IOFF+NI-1 2168C? WRITE(6,*) ' IOB, JOB, KOB, LOB', IOB,JOB,KOB,LOB 2169 INT_OUT = I2EAD(IOB,JOB,KOB,LOB) 2170 INT_IN = INT_IN + 1 2171C? WRITE(6,*) ' INT_OUT, INT_IN ', INT_OUT,INT_IN 2172C? WRITE(6,*) ' KINT2-1+INT_OUT = ',KINT2-1+INT_OUT 2173 WORK(KINT2-1+INT_OUT) = XINT(INT_IN) 2174 END DO 2175 END DO 2176 END DO 2177 END DO 2178* 2179 CALL QEXIT('PUTIN') 2180 END 2181 SUBROUTINE MAKE_LUCITA_INTEGRALS(CMO,WORK,LWORK) 2182C 2183C Dec 2010 : get EMY and get FCAC and H2AC matrices in Dalton 2184C 2185 use lucita_cfg 2186#include "implicit.h" 2187 COMMON/CECORE/ECORE,ECORE_ORIG,ECORE_H,ECORE_HEX 2188#include "maxorb.h" 2189#include "infinp.h" 2190 2191 DIMENSION CMO(*), WORK(LWORK) 2192 2193C Used from include files: 2194C priunit: LUPRI 2195C inforb: NCMOT, nsym 2196C inftap: LUINTM, FNINTM 2197C mxpdim: MXPOBS 2198C orbinp: ntoobs 2199#include "priunit.h" 2200#include "inforb.h" 2201#include "inftap.h" 2202 2203 2204 LOGICAL CLOSE_MOTWOINT, DISKH2 2205C 2206C ------------------------------------------------------- 2207C 2208 KFRSAV = 1 2209 KFREE = KFRSAV 2210 LFREE = LWORK 2211C 2212 IF (LUINTM .LE. 0) THEN 2213 CALL GPOPEN(LUINTM,FNINTM,'OLD',' ', 2214 & 'UNFORMATTED',IDUMMY,.FALSE.) 2215 CLOSE_MOTWOINT = .TRUE. 2216 ELSE 2217 CLOSE_MOTWOINT = .FALSE. 2218 END IF 2219 2220! allocate scratch memory and read-in 1/2-electron integrals from dalton-sirius files 2221 call memget('REAL',kfcac,nnashx ,WORK,kfree,lfree) 2222 call memget('REAL',kh2ac,nnashx**2,WORK,kfree,lfree) 2223 2224 DISKH2 = .FALSE. 2225 CALL CIINTS(DISKH2,CMO,EMY,WORK(kfcac), 2226 & WORK(kh2ac),WORK(kfree),lfree) 2227 2228 IF (CLOSE_MOTWOINT) THEN 2229 CALL GPCLOSE(LUINTM,'KEEP') 2230 END IF 2231 2232 2233#ifdef LUCI_DEBUG 2234 WRITE(LUPRI,'(//A)') 'MAKE_LUCITA_INTEGRALS: FCAC matrix' 2235 CALL OUTPAK(work(kfcac),NASHT,-1,LUPRI) 2236 WRITE(LUPRI,'(//A)') 'MAKE_LUCITA_INTEGRALS: H2AC matrix' 2237 CALL OUTPUT(work(kh2ac),1,NNASHX,1,NNASHX,NNASHX, 2238 & NNASHX,-1,LUPRI) 2239#endif 2240 2241! Write one- and two-electron integrals to file FCIDUMP 2242 2243 if (lucita_cfg_fci_dump) 2244 & call write_fcidump_file(work(kfcac),NASHT,work(kh2ac),EMY) 2245 2246 LUINT_LUCITA = -1 2247 CALL GPOPEN(LUINT_LUCITA,'INTEGRALS_LUCITA','NEW',' ', 2248 & 'UNFORMATTED',IDUMMY,.FALSE.) 2249 WRITE (LUINT_LUCITA) ('********',I=1,3), 'LUCIINTS' 2250 WRITE (LUINT_LUCITA) EMY, NNASHX, NSYM, DUM, DUM 2251 CALL WRITT(LUINT_LUCITA, NNASHX, WORK(KFCAC)) 2252 CALL WRITT(LUINT_LUCITA, NNASHX*NNASHX, WORK(KH2AC)) 2253 WRITE (LUINT_LUCITA) ('********',I=1,3), 'EODATA ' 2254 CALL GPCLOSE(LUINT_LUCITA, 'KEEP') 2255! release scratch space 2256 call memrel('lucita-dalton-1/2int-write.done',WORK, 2257 & kfrsav,kfrsav,kfree,lfree) 2258 2259 END ! SUBROUTINE MAKE_LUCITA_INTEGRALS 2260*********************************************************************** 2261 2262 SUBROUTINE LUCITA_GETINT_DALTON(EMY,X1INT,X2INT, 2263 & tmp_dal_work,LWORK, 2264 & int1_mcscf, 2265 & int2_mcscf, 2266 & mcscf_provides_integrals,iintden) 2267C 2268C Dec 2010 : get EMY and get FCAC and H2AC matrices in Dalton 2269C 2270 use lucita_mc_response_cfg 2271#include "implicit.h" 2272 2273 DIMENSION X1INT(*), X2INT(*), tmp_dal_work(lwork) 2274 real(8), intent(inout) :: int1_mcscf(*) 2275 real(8), intent(inout) :: int2_mcscf(*) 2276 integer, intent(in) :: iintden 2277 logical, intent(in) :: mcscf_provides_integrals 2278 2279C Used from include files: 2280C priunit: LUPRI 2281C inforb: NCMOT, nsym 2282C inftap: LUINTM, FNINTM 2283C mxpdim: MXPOBS 2284C orbinp: ntoobs 2285C cgas: ngas 2286#include "priunit.h" 2287#include "inforb.h" 2288#include "inftap.h" 2289#include "mxpdim.inc" 2290#include "orbinp.inc" 2291#include "cgas.inc" 2292#include "oper.inc" 2293 COMMON/CINTFO/I12S,I34S,I1234S,NINT1,NINT2,NBINT1,NBINT2 2294 INTEGER ADASX,ASXAD,ADSXA,SXSXDX,SXDXSX 2295 COMMON/CSMPRD/ADASX(MXPOBS,MXPOBS),ASXAD(MXPOBS,2*MXPOBS), 2296 & ADSXA(MXPOBS,2*MXPOBS), 2297 & SXSXDX(2*MXPOBS,2*MXPOBS),SXDXSX(2*MXPOBS,4*MXPOBS) 2298 2299 LOGICAL CLOSE_MOTWOINT, DISKH2 2300! scratch space 2301 integer, allocatable :: mult_ptg(:,:) 2302 integer :: is2e_active 2303C 2304C ------------------------------------------------------- 2305C 2306! allocate scratch matrix 2307 allocate(mult_ptg(nsym,nsym)) 2308 2309! set multiplication table 2310 mult_ptg = 0 2311 call set_ptg_multiplication_table(mult_ptg,nsym) 2312 2313 is2e_active = i12 2314 2315 KFRSAV = 1 2316 KFREE = KFRSAV 2317 LFREE = LWORK 2318 2319! write(lupri,*) 'mcscf_provides_integrals = ', 2320! & mcscf_provides_integrals 2321 2322 if(.not.mcscf_provides_integrals)then 2323! allocate scratch memory and read-in 1/2-electron integrals from dalton-sirius files 2324 call memget('REAL',kfcac,nnashx ,tmp_dal_work,kfree,lfree) 2325 call memget('REAL',kh2ac,nnashx**2,tmp_dal_work,kfree,lfree) 2326 2327 LUINT_LUCITA = -1 2328 CALL GPOPEN(LUINT_LUCITA,'INTEGRALS_LUCITA','OLD',' ', 2329 & 'UNFORMATTED',IDUMMY,.FALSE.) 2330 CALL MOLLAB('LUCIINTS',LUINT_LUCITA,LUPRI) 2331 READ (LUINT_LUCITA) EMY, NNASHX_x, NSYM_x 2332 IF (NNASHX .NE. NNASHX_x) THEN 2333 call quit('NNASHX in common .ne. NNASHX on INTEGRALS_LUCITA') 2334 END IF 2335 IF (NSYM .NE. NSYM_x) THEN 2336 call quit('NSYM in common .ne. NSYM on INTEGRALS_LUCITA') 2337 END IF 2338 CALL READT(LUINT_LUCITA, NNASHX, tmp_dal_WORK(KFCAC)) 2339 CALL READT(LUINT_LUCITA, NNASHX*NNASHX, tmp_dal_WORK(KH2AC)) 2340 CALL GPCLOSE(LUINT_LUCITA, 'KEEP') 2341 2342!#define LUCI_DEBUG 2343#ifdef LUCI_DEBUG 2344 WRITE(LUPRI,'(//A)') 'LUCITA_GETINT_DALTON: FCAC matrix' 2345 WRITE(LUPRI,'(A,2i8)') 2346 & 'dimension: NASHT,NASHT * NASHT',NASHT,NASHT**2 2347 CALL OUTPAK(tmp_dal_work(kfcac),NASHT,-1,LUPRI) 2348 WRITE(LUPRI,'(//A)') 'LUCITA_GETINT_DALTON: H2AC matrix' 2349 WRITE(LUPRI,'(A,2i8)') 2350 & 'dimension: NNASHX, NNASHX**2',NNASHX,NNASHX**2 2351 CALL OUTPUT(tmp_dal_work(kh2ac),1,NNASHX,1,NNASHX,NNASHX, 2352 & NNASHX,-1,LUPRI) 2353#endif 2354!#undef LUCI_DEBUG 2355 2356!-------------------------------------------------------------- 2357! transfer 1e- and 2e-integrals from Dalton to LUCITA order 2358!-------------------------------------------------------------- 2359 call sirint2luint(tmp_dal_work(kfcac),tmp_dal_work(kh2ac), 2360 & dummy,x1int,x2int,mult_ptg,adsxa,nasht,nnashx, 2361 & ntoobs,nsym,mxpobs,nint1,nint2,ireost,ireots, 2362 & ngas,is2e_active,iintden,-1) 2363! release scratch space 2364 call memrel('lucita-dalton-1/2int-interface.done',tmp_dal_work, 2365 & kfrsav,kfrsav,kfree,lfree) 2366 2367 else ! mcscf_provides_integrals 2368! write(lupri,*) 'distribute ints from mc environment...' 2369 call sirint2luint(int1_mcscf,int2_mcscf, 2370 & dummy,x1int,x2int,mult_ptg,adsxa,nasht,nnashx, 2371 & ntoobs,nsym,mxpobs,nint1,nint2,ireost,ireots, 2372 & ngas,is2e_active,iintden, 2373 & lucita_mc_response_prp_int) 2374 end if 2375 2376 deallocate(mult_ptg) 2377 2378 END 2379*********************************************************************** 2380 2381 subroutine lucita_srdft_h1_adapt(X1INT,x1int_scr,int1_mcscf) 2382C 2383#include "implicit.h" 2384 2385 real(8), intent(inout) :: x1int(*) 2386 real(8), intent(inout) :: x1int_scr(*) 2387 real(8), intent(inout) :: int1_mcscf(*) 2388 2389C Used from include files: 2390C priunit: LUPRI 2391C inforb: NCMOT, nsym 2392C inftap: LUINTM, FNINTM 2393C mxpdim: MXPOBS 2394C orbinp: ntoobs 2395C cgas: ngas 2396#include "priunit.h" 2397#include "inforb.h" 2398#include "inftap.h" 2399#include "mxpdim.inc" 2400#include "orbinp.inc" 2401#include "cgas.inc" 2402 COMMON/CINTFO/I12S,I34S,I1234S,NINT1,NINT2,NBINT1,NBINT2 2403 INTEGER ADASX,ASXAD,ADSXA,SXSXDX,SXDXSX 2404 COMMON/CSMPRD/ADASX(MXPOBS,MXPOBS),ASXAD(MXPOBS,2*MXPOBS), 2405 & ADSXA(MXPOBS,2*MXPOBS), 2406 & SXSXDX(2*MXPOBS,2*MXPOBS),SXDXSX(2*MXPOBS,4*MXPOBS) 2407 2408! scratch space 2409 integer, allocatable :: mult_ptg(:,:) 2410 integer :: is2e_active 2411 integer :: iintden = 1 2412C 2413C ------------------------------------------------------- 2414C 2415! allocate scratch matrix 2416 allocate(mult_ptg(nsym,nsym)) 2417 2418! set multiplication table 2419 mult_ptg = 0 2420 call set_ptg_multiplication_table(mult_ptg,nsym) 2421 2422 is2e_active = 1 ! no 2-electron integrals 2423 2424!----------------------------------------------------------------- 2425! transfer 1e- srDFT contributions from Dalton to LUCITA order 2426!----------------------------------------------------------------- 2427 call sirint2luint(x1int,dummy, 2428 & dummy,x1int_scr,dummy,mult_ptg, 2429 & adsxa,nasht,nnashx, 2430 & ntoobs,nsym,mxpobs,nint1,nxxxx,ireost,ireots, 2431 & ngas,is2e_active,iintden,-1) 2432 2433! scale 1-electron integrals with the sr-dft contributions 2434 call daxpy(nint1,1.0d0,x1int_scr,1,int1_mcscf,1) 2435 2436 deallocate(mult_ptg) 2437 2438 END 2439*********************************************************************** 2440 2441 subroutine lucita_putdens_generic(rho1lu,rho2lu,rho1gen,rho2gen, 2442 & pv_full_scratch, 2443 & isrho2_active,iintden,rhotype, 2444 & state_id) 2445! 2446! Sep 2011 : driver for sorting the 1-/2-particle density matrices from 2447! LUCITA to a generic format (here: SIRIUS) 2448! 2449#ifdef VAR_MPI 2450#ifdef USE_MPI_MOD_F90 2451 use mpi 2452#include "implicit.h" 2453#else 2454#include "implicit.h" 2455#include "mpif.h" 2456#endif 2457#endif 2458#include "parluci.h" 2459 2460 real(8), intent(in) :: rho1lu(*) 2461 real(8), intent(in) :: rho2lu(*) 2462 real(8), intent(out) :: rho1gen(*) 2463 real(8), intent(out) :: rho2gen(*) 2464 real(8), intent(inout) :: pv_full_scratch(*) 2465 integer, intent(in) :: iintden 2466 integer, intent(in) :: isrho2_active 2467 integer, intent(in) :: rhotype 2468 integer, intent(in) :: state_id 2469 2470! Used from include files: 2471! priunit: LUPRI 2472! inforb: NCMOT, nsym 2473! inftap: LUINTM, FNINTM 2474! mxpdim: MXPOBS 2475! orbinp: ntoobs 2476 2477#include "priunit.h" 2478#include "inforb.h" 2479#include "inftap.h" 2480#include "mxpdim.inc" 2481#include "orbinp.inc" 2482#include "cgas.inc" 2483 COMMON/CINTFO/I12S,I34S,I1234S,NINT1,NINT2,NBINT1,NBINT2 2484 INTEGER ADASX,ASXAD,ADSXA,SXSXDX,SXDXSX 2485 COMMON/CSMPRD/ADASX(MXPOBS,MXPOBS),ASXAD(MXPOBS,2*MXPOBS), 2486 & ADSXA(MXPOBS,2*MXPOBS), 2487 & SXSXDX(2*MXPOBS,2*MXPOBS),SXDXSX(2*MXPOBS,4*MXPOBS) 2488 2489! scratch space 2490 integer, allocatable :: mult_ptg(:,:) 2491 integer :: noint_tmp 2492 character (len=25) :: file_name 2493! 2494! ------------------------------------------------------- 2495! 2496! allocate scratch matrix 2497 allocate(mult_ptg(nsym,nsym)) 2498 2499!--------------------------------------------------------------------- 2500! transfer 1e- and 2e-density matrices from LUCITA to SIRIUS order 2501!--------------------------------------------------------------------- 2502 2503! set multiplication table 2504 mult_ptg = 0 2505 call set_ptg_multiplication_table(mult_ptg,nsym) 2506 2507 call dzero(rho1gen,nnashx**1) 2508 if(isrho2_active .gt. 1) call dzero(rho2gen,nnashx**2) 2509 2510 call sirint2luint(rho1gen,rho2gen,pv_full_scratch, 2511 & rho1lu,rho2lu,mult_ptg,adsxa,nasht,nnashx, 2512 & ntoobs,nsym,mxpobs,nint1,nint2,ireost, 2513 & ireots,ngas,isrho2_active,iintden, 2514 & rhotype) 2515 2516 deallocate(mult_ptg) 2517 2518! put 1e- and 2e-density matrices to file 2519! --------------------------------------- 2520 2521 write(file_name,'(a17,i3,a1,i4)') 2522 & 'DENSITIES_LUCITA-',state_id,'.',luci_myproc 2523 do i = 18, 20 2524 if(file_name(i:i).eq. ' ') file_name(i:i) = '0' 2525 end do 2526 do i = 22, 25 2527 if(file_name(i:i).eq. ' ') file_name(i:i) = '0' 2528 end do 2529 2530 LUINT_LUCITA = -1 2531 CALL GPOPEN(LUINT_LUCITA,file_name,'UNKNOWN',' ', 2532 & 'UNFORMATTED',IDUMMY,.FALSE.) 2533 WRITE (LUINT_LUCITA) ('********',I=1,3), 'LUCIDENS' 2534 CALL WRITT(LUINT_LUCITA, NNASHX, rho1gen) 2535 if(isrho2_active .gt. 1) 2536 &CALL WRITT(LUINT_LUCITA, NNASHX*NNASHX, rho2gen) 2537 WRITE (LUINT_LUCITA) ('********',I=1,3), 'EODATA ' 2538 CALL GPCLOSE(LUINT_LUCITA, 'KEEP') 2539 2540!#define LUCI_DEBUG 2541#ifdef LUCI_DEBUG 2542 WRITE(LUPRI,'(//A)') 2543 & 'lucita_putdens_generic: 1-particle density matrix' 2544 WRITE(LUPRI,'(A,2i8)') 2545 & 'dimension: NASHT,NASHT * NASHT',NASHT,NASHT**2 2546 CALL OUTPAK(rho1gen,NASHT,-1,LUPRI) 2547 WRITE(LUPRI,'(//A)') 2548 & 'lucita_putdens_generic: 2-particle density matrix' 2549 WRITE(LUPRI,'(A,2i8)') 2550 & 'dimension: NNASHX, NNASHX**2',NNASHX,NNASHX**2 2551 CALL OUTPUT(rho2gen,1,NNASHX,1,NNASHX,NNASHX, 2552 & NNASHX,-1,LUPRI) 2553#endif 2554!#undef LUCI_DEBUG 2555 2556 END 2557*********************************************************************** 2558 2559 subroutine lucita_spinputdens_1p(rho1slu,rho2lu,rho1gen,rho2gen, 2560 & pv_full_scratch, 2561 & isrho2_active,iintden,rhotype, 2562 & state_id,iab) 2563! 2564! Aug 2012 : driver for sorting the 1-particle spin-density matrices from 2565! LUCITA to a generic format (here: SIRIUS) 2566! 2567#ifdef VAR_MPI 2568#ifdef USE_MPI_MOD_F90 2569 use mpi 2570#include "implicit.h" 2571#else 2572#include "implicit.h" 2573#include "mpif.h" 2574#endif 2575#endif 2576#include "parluci.h" 2577 2578 real(8), intent(in) :: rho1slu(*) 2579 real(8), intent(in) :: rho2lu(*) 2580 real(8), intent(out) :: rho1gen(*) 2581 real(8), intent(out) :: rho2gen(*) 2582 real(8), intent(inout) :: pv_full_scratch(*) 2583 integer, intent(in) :: iintden 2584 integer, intent(in) :: isrho2_active 2585 integer, intent(in) :: rhotype 2586 integer, intent(in) :: state_id 2587 integer, intent(in) :: iab 2588 2589! Used from include files: 2590! priunit: LUPRI 2591! inforb: NCMOT, nsym 2592! inftap: LUINTM, FNINTM 2593! mxpdim: MXPOBS 2594! orbinp: ntoobs 2595 2596#include "priunit.h" 2597#include "inforb.h" 2598#include "inftap.h" 2599#include "mxpdim.inc" 2600#include "orbinp.inc" 2601#include "cgas.inc" 2602 COMMON/CINTFO/I12S,I34S,I1234S,NINT1,NINT2,NBINT1,NBINT2 2603 INTEGER ADASX,ASXAD,ADSXA,SXSXDX,SXDXSX 2604 COMMON/CSMPRD/ADASX(MXPOBS,MXPOBS),ASXAD(MXPOBS,2*MXPOBS), 2605 & ADSXA(MXPOBS,2*MXPOBS), 2606 & SXSXDX(2*MXPOBS,2*MXPOBS),SXDXSX(2*MXPOBS,4*MXPOBS) 2607 2608! scratch space 2609 integer, allocatable :: mult_ptg(:,:) 2610 integer :: noint_tmp 2611 character (len=25) :: file_name 2612! 2613! ------------------------------------------------------- 2614! 2615! allocate scratch matrix 2616 allocate(mult_ptg(nsym,nsym)) 2617 2618!--------------------------------------------------------------------- 2619! transfer 1e-spin density matrix from LUCITA to Sirius order 2620!--------------------------------------------------------------------- 2621 2622! set multiplication table 2623 mult_ptg = 0 2624 call set_ptg_multiplication_table(mult_ptg,nsym) 2625 2626 call dzero(rho1gen,nnashx**1) 2627 if(isrho2_active .gt. 1) call dzero(rho2gen,nnashx**2) 2628 2629 call sirint2luint(rho1gen,rho2gen,pv_full_scratch, 2630 & rho1slu,rho2lu,mult_ptg,adsxa,nasht,nnashx, 2631 & ntoobs,nsym,mxpobs,nint1,nint2,ireost, 2632 & ireots,ngas,isrho2_active,iintden, 2633 & rhotype) 2634 2635 deallocate(mult_ptg) 2636 2637! put 1e- and 2e-density matrices to file 2638! --------------------------------------- 2639 2640 select case(iab) 2641 case(1) ! alpha spin 2642 write(file_name,'(a17,i3,a1,i4)') 2643 & 'SPINDENSITY-A-1p-',state_id,'.',luci_myproc 2644 case(2) ! beta spin 2645 write(file_name,'(a17,i3,a1,i4)') 2646 & 'SPINDENSITY-B-1p-',state_id,'.',luci_myproc 2647 case default ! ??? ask einstein or dirac 2648 stop ' wrong program/hamiltonian - spin is a good quantum num' 2649 end select 2650 2651 do i = 18, 20 2652 if(file_name(i:i).eq. ' ') file_name(i:i) = '0' 2653 end do 2654 do i = 22, 25 2655 if(file_name(i:i).eq. ' ') file_name(i:i) = '0' 2656 end do 2657 2658 LUINT_LUCITA = -1 2659 CALL GPOPEN(LUINT_LUCITA,file_name,'UNKNOWN',' ', 2660 & 'UNFORMATTED',IDUMMY,.FALSE.) 2661 WRITE (LUINT_LUCITA) ('********',I=1,3), 'LUCIDENS' 2662 CALL WRITT(LUINT_LUCITA, NNASHX, rho1gen) 2663 WRITE (LUINT_LUCITA) ('********',I=1,3), 'EODATA ' 2664 CALL GPCLOSE(LUINT_LUCITA, 'KEEP') 2665 2666!#define LUCI_DEBUG 2667#ifdef LUCI_DEBUG 2668 WRITE(LUPRI,'(//A,i2)') 2669 & 'lucita_putdens_generic: 1-particle spin density matrix :',iab 2670 WRITE(LUPRI,'(A,2i8)') 2671 & 'dimension: NASHT, triangular packed',NASHT 2672 CALL OUTPAK(rho1gen,NASHT,-1,LUPRI) 2673#undef LUCI_DEBUG 2674#endif 2675 2676 END 2677*********************************************************************** 2678 2679 subroutine lucita_putdens_ensemble(rho1_ensemble) 2680! 2681#include "implicit.h" 2682 2683 real(8), intent(in) :: rho1_ensemble(*) 2684 2685! Used from include files: 2686! priunit: LUPRI 2687! inforb: NCMOT, nsym 2688! inftap: LUINTM, FNINTM 2689! mxpdim: MXPOBS 2690! orbinp: ntoobs 2691 2692#include "priunit.h" 2693#include "inforb.h" 2694#include "inftap.h" 2695#include "mxpdim.inc" 2696#include "orbinp.inc" 2697#include "cgas.inc" 2698 COMMON/CINTFO/I12S,I34S,I1234S,NINT1,NINT2,NBINT1,NBINT2 2699 INTEGER ADASX,ASXAD,ADSXA,SXSXDX,SXDXSX 2700 COMMON/CSMPRD/ADASX(MXPOBS,MXPOBS),ASXAD(MXPOBS,2*MXPOBS), 2701 & ADSXA(MXPOBS,2*MXPOBS), 2702 & SXSXDX(2*MXPOBS,2*MXPOBS),SXDXSX(2*MXPOBS,4*MXPOBS) 2703 2704! scratch space 2705 character (len=25) :: file_name 2706! 2707! ------------------------------------------------------- 2708! 2709! put 1e-density matrix to file 2710! ----------------------------- 2711 2712 LUINT_LUCITA = 99 2713 write(file_name,'(a16)') 2714 & 'ENSEMBLE-DENSITY' 2715 2716 open(file=file_name(1:16),unit=luint_lucita,status='replace', 2717 & form='unformatted',action='write',position='rewind') 2718 2719 WRITE (LUINT_LUCITA) ('********',I=1,3), 'LUCIDENS' 2720 CALL WRITT(LUINT_LUCITA, NNASHX, rho1_ensemble) 2721 WRITE (LUINT_LUCITA) ('********',I=1,3), 'EODATA ' 2722 2723 close(luint_lucita,status='keep') 2724 2725 END 2726*********************************************************************** 2727 2728 subroutine sirint2luint(h1_in,h2_in,pv_full_scratch,h1_out,h2_out, 2729 & mult_ptg,return_sym_of_jorb,nr_actorb_tot, 2730 & nnashx,nr_orb_per_sym,nr_sym, 2731 & max_orb_per_sym,nr_1int,nr_2int, 2732 & ireost,ireots,nr_gas,is2e_active,iintden, 2733 & rhotype) 2734! 2735! purpose: iintden == 1: reorder integrals from dalton-sirius format to internal 2736! lucita format. 2737! 1-electron integrals: on output symmetry-reduced list 2738! 2-electron integrals: on output symmetry-reduced list 2739! iintden == 2: reorder particle-density matrices from to internal lucita format to 2740! dalton-sirius format. 2741! rhotype controls the format: 2742! ==> 1: rho2(ij,kl) density matrix 2743! ==> 2: rho2(i,j,k,l) used for response 2744! ==> 3: rho2(ij,kl) symmetrized transition density matrix 2745! 2746! ----------------------------------------------------------------- 2747 implicit none 2748#include "priunit.h" 2749 integer :: nr_actorb_tot 2750 integer :: nnashx 2751 integer :: nr_sym 2752 integer :: max_orb_per_sym 2753 integer :: nr_1int 2754 integer :: nr_2int 2755 integer :: return_sym_of_jorb(max_orb_per_sym,max_orb_per_sym*2) 2756 integer :: nr_orb_per_sym(max_orb_per_sym) 2757 integer :: mult_ptg(nr_sym,nr_sym) 2758 integer :: iintden 2759 integer :: nr_gas 2760 integer :: is2e_active 2761 integer :: ireost(*) 2762 integer :: ireots(*) 2763 integer :: rhotype 2764! real(8) :: h1_in(nnashx) 2765! real(8) :: h2_in(nnashx,nnashx) 2766 real(8) :: h1_in(*) 2767 real(8) :: h2_in(*) 2768 real(8) :: pv_full_scratch(nr_actorb_tot,nr_actorb_tot, 2769 & nr_actorb_tot,nr_actorb_tot) 2770 real(8) :: h1_out(*) 2771 real(8) :: h2_out(*) 2772 integer :: mytest 2773! ------------------------------------------------------------------ 2774 2775!#define LUCI_DEBUG 2776#ifdef LUCI_DEBUG 2777 if(iintden .eq. 2)then 2778 write(lupri,*) ' 1-el density matrix' 2779 call wrtmatmn(h1_out,1,nr_actorb_tot**2,1, 2780 & nr_actorb_tot**2,lupri) 2781 if(is2e_active .gt. 1)then 2782 write(lupri,*) ' 2-el density matrix' 2783 call wrtmatmn(h2_out,1,nr_actorb_tot**2*(nr_actorb_tot**2+1)/2, 2784 & 1,nr_actorb_tot**2*(nr_actorb_tot**2+1)/2,lupri) 2785 end if 2786 end if 2787#endif 2788 2789! 1-electron integrals/density matrix 2790! ----------------------------------- 2791 call sir1int2lu1int(h1_in,h1_out,return_sym_of_jorb,nnashx, 2792 & nr_orb_per_sym,nr_sym,max_orb_per_sym, 2793 & nr_1int,nr_actorb_tot,ireost,iintden,rhotype) 2794 2795! 2-electron integrals/density matrix 2796! ----------------------------------- 2797 if(is2e_active .gt. 1)then 2798 if(iintden .eq. 1)then 2799 call sir2int2lu2int(h2_in,h2_out,mult_ptg,nr_actorb_tot, 2800 & nnashx,nr_orb_per_sym,nr_sym, 2801 & max_orb_per_sym,nr_2int) 2802 else 2803 call lu2dens2sir2dens(h2_in,h2_out,pv_full_scratch, 2804 & nr_actorb_tot,ireost,ireots, 2805 & nr_sym,nr_gas,rhotype) 2806 end if 2807 end if 2808 2809#ifdef LUCI_DEBUG 2810 if(iintden .eq. 1)then 2811 write(lupri,*) ' symmetry reduced 1-el' 2812 call wrtmatmn(h1_out,1,nr_1int,1,nr_1int,lupri) 2813 if(is2e_active .gt. 1)then 2814 write(lupri,*) ' symmetry reduced 2-el' 2815 call wrtmatmn(h2_out,1,nr_2int,1,nr_2int,lupri) 2816 end if 2817 end if 2818#undef LUCI_DEBUG 2819#endif 2820!#define LUCI_DEBUG 2821#ifdef LUCI_DEBUG 2822 mytest = 0 2823 if(iintden > 1)then 2824 write(lupri,*) ' 1-el density matrix' 2825 call outpak(h1_in,nr_actorb_tot,-1,lupri) 2826 if(is2e_active .gt. 1 .and. mytest == 1)then 2827 write(lupri,*) ' 2-el density matrix' 2828 call output(h2_in,1,nnashx,1,nnashx,nnashx, 2829 & nnashx,-1,lupri) 2830 end if 2831 end if 2832#undef LUCI_DEBUG 2833#endif 2834 2835 end 2836 2837*********************************************************************** 2838 subroutine sir1int2lu1int(h1_in,h1_out,return_sym_of_jorb,nnashx, 2839 & nr_orb_per_sym,nr_sym,max_orb_per_sym, 2840 & nr_1int,nr_actorb_tot,ireost,iintden, 2841 & rhotype) 2842! 2843! purpose: reorder integrals from dalton-sirius format to internal 2844! lucita format. 2845! 2846! 1-electron integrals: on output symmetry-reduced list 2847! 2848! note from stefan: i am not sure whether my code actually 2849! works for non-totally symmetric 1-electron operators as 2850! well, i.e. if we do property calculations. 2851! this should be checked with due care. 2852! 2853! ----------------------------------------------------------------- 2854 implicit none 2855#include "priunit.h" 2856 integer :: nnashx 2857 integer :: nr_sym 2858 integer :: max_orb_per_sym 2859 integer :: nr_1int 2860 integer :: return_sym_of_jorb(max_orb_per_sym,max_orb_per_sym*2) 2861 integer :: nr_orb_per_sym(max_orb_per_sym) 2862 integer :: nr_actorb_tot 2863 integer :: ireost(*) 2864 integer :: iintden 2865 integer :: rhotype 2866 real(8) :: h1_in(*) 2867 real(8) :: h1_out(*) 2868! ----------------------------------------------------------------- 2869 integer :: isym, jsym, ijsym 2870 integer :: offset_ij, offset_internal, isorb, jsorb 2871 integer :: nr_act_isym, nr_act_jsym, orb_tmp 2872 integer :: i_index, j_index 2873 integer :: loop_counter_i 2874 character (len=90) :: error_message 2875! ----------------------------------------------------------------- 2876 2877!#define LUCI_DEBUG 2878#ifdef LUCI_DEBUG 2879 if(rhotype > 0)then 2880 write(lupri,*) ' distributing 1p-dens mat elements from h1_out:' 2881 call wrtmatmn(h1_out,nr_actorb_tot,nr_actorb_tot, 2882 & nr_actorb_tot,nr_actorb_tot,lupri) 2883 else 2884 write(lupri,*) ' distributing 1e-ints from h1_in:' 2885 do isym = 1, nnashx 2886 write(lupri,*) ' h1_in(',isym,') = ',h1_in(isym) 2887 end do 2888 end if 2889!#undef LUCI_DEBUG 2890#endif 2891 2892! 1-electron integrals/density matrix elements 2893! -------------------------------------------- 2894 orb_tmp = 0 2895 offset_ij = 0 2896 offset_internal = 1 2897 ijsym = 1 2898 2899! insert check for symmetry of 1-electron operator - quit if not 2900! totally symmetric for the time being 2901 if(ijsym .gt. 1)then 2902 error_message = '*** error in integral resorting for dalton'// 2903 & ' 1-e ints. operator is not totally symmetric.' 2904 call quit(error_message) 2905 end if 2906 2907 do isym = 1, nr_sym 2908 2909 2910 nr_act_isym = nr_orb_per_sym(isym) 2911 jsym = return_sym_of_jorb(isym,ijsym) 2912 2913#ifdef LUCI_DEBUG 2914 write(lupri,*) ' isym, jsym',isym, jsym 2915 write(lupri,*) ' nr_act_isym', nr_act_isym 2916#endif 2917 2918 if(jsym >= isym)then 2919 2920 nr_act_jsym = nr_orb_per_sym(jsym) 2921 2922#ifdef LUCI_DEBUG 2923 write(lupri,*) ' nr_act_jsym', nr_act_jsym 2924#endif 2925 2926 do jsorb = 1, nr_act_jsym 2927 loop_counter_i = jsorb 2928 if(rhotype > 0) loop_counter_i = nr_act_isym 2929 do isorb = 1, loop_counter_i 2930! calculate offset on dalton 1e-array (lower triangular matrix) 2931 i_index = (orb_tmp+isorb) 2932 j_index = ((orb_tmp+jsorb)*(orb_tmp+jsorb-1)/2) 2933! no triangular packing if rhotype > 1 2934 if(rhotype > 0) j_index = ((orb_tmp+jsorb-1)) 2935! & * ( orb_tmp+jsorb )) 2936 & * nr_actorb_tot 2937 offset_ij = i_index + j_index 2938 2939 if(iintden .eq. 1)then ! 1e integral 2940 h1_out(offset_internal) = h1_in(offset_ij) 2941 else ! 1p density matrix element 2942! 2943! LUCITA: type-symmetry ordering <==> SIRIUS (or other programs) symmetry-type 2944! array ireost takes care ot that... 2945! 2946 offset_internal = 2947 & (ireost(isorb+orb_tmp)-1) * nr_actorb_tot + 2948 & ireost(jsorb+orb_tmp) 2949 if(rhotype > 1) offset_internal = 2950 & (ireost(jsorb+orb_tmp)-1) * nr_actorb_tot + 2951 & ireost(isorb+orb_tmp) 2952 2953 h1_in(offset_ij) = h1_out(offset_internal) 2954 end if 2955 2956#ifdef LUCI_DEBUG 2957 write(lupri,*)'offset_ij, offset_int,jsorb,isorb,h1,switch', 2958 & offset_ij, offset_internal,jsorb,isorb, 2959 & h1_out(offset_internal), h1_in(offset_ij), iintden 2960#endif 2961 2962! count the number of non-vanishing 1-electron integrals 2963 offset_internal = offset_internal + 1 2964 end do 2965 end do 2966 end if ! jsym >= isym 2967 2968! total # of active orbitals for each isym 2969 orb_tmp = orb_tmp + nr_act_isym 2970 end do 2971 2972! final # of 1-el ints 2973 nr_1int = offset_internal - 1 2974 2975!#define LUCI_DEBUG 2976#ifdef LUCI_DEBUG 2977 if(rhotype <= 0)then 2978 write(lupri,*) ' final 1e- list: # of ints =', nr_1int 2979 call wrtmatmn(h1_out,1,nr_1int,1,nr_1int,lupri) 2980 end if 2981#undef LUCI_DEBUG 2982#endif 2983 2984! symmetrize 1p-dens matrix 2985 if(rhotype > 0)then 2986 2987! write(lupri,*) ' temp 1p-tdm: ' 2988! call wrtmatmn(h1_in,nr_actorb_tot,nr_actorb_tot, 2989! & nr_actorb_tot,nr_actorb_tot,lupri) 2990 2991 call trpad(h1_in,1.0d0,nr_actorb_tot) 2992 2993! write(lupri,*) ' temp 1p-tdm: part 2' 2994! call wrtmatmn(h1_in,nr_actorb_tot,nr_actorb_tot, 2995! & nr_actorb_tot,nr_actorb_tot,lupri) 2996 2997 if(rhotype == 1) call dscal(nr_actorb_tot**2,0.5d0,h1_in,1) 2998 2999! write(lupri,*) ' temp 1p-tdm: part 3' 3000! call wrtmatmn(h1_in,nr_actorb_tot,nr_actorb_tot, 3001! & nr_actorb_tot,nr_actorb_tot,lupri) 3002 3003! pack symmetrized 1p-dens matrix (h1_in) as triangular matrix h1_out 3004 call dsitsp(nr_actorb_tot,h1_in,h1_out) 3005! save triangular matrix in h1_in 3006 call dzero(h1_in,nr_actorb_tot**2) 3007 call dcopy((nr_actorb_tot*(nr_actorb_tot+1))/2,h1_out,1,h1_in,1) 3008 end if 3009 3010 end 3011*********************************************************************** 3012 3013 subroutine sir2int2lu2int(h2_in,h2_out,mult_ptg,nr_actorb_tot, 3014 & nnashx,nr_orb_per_sym,nr_sym, 3015 & max_orb_per_sym,nr_2int) 3016! 3017! purpose: reorder integrals from dalton-sirius format to internal 3018! lucita format. 3019! 3020! 2-electron integrals: on output symmetry-reduced list 3021! 3022! ----------------------------------------------------------------- 3023 implicit none 3024#include "priunit.h" 3025 integer :: nr_actorb_tot 3026 integer :: nnashx 3027 integer :: nr_sym 3028 integer :: max_orb_per_sym 3029 integer :: nr_2int 3030 integer :: nr_orb_per_sym(max_orb_per_sym) 3031 integer :: mult_ptg(nr_sym,nr_sym) 3032 real(8) :: h2_in(nnashx,nnashx) 3033 real(8) :: h2_out(*) 3034! ----------------------------------------------------------------- 3035 integer :: isym, jsym, ksym, lsym, llsym, ijsym, klsym, ijklsym 3036 integer :: offset_internal, offset_ijkl 3037 integer :: nr_act_isym, nr_act_jsym, nr_act_ksym, nr_act_lsym 3038 integer :: ioff, joff, koff, loff 3039 integer :: ij_res, kl_res, ijkl_res 3040 integer :: i_first, j_first, j_last, l_last 3041 integer :: i, j, k, l, ij, kl 3042! ----------------------------------------------------------------- 3043 3044!#define LUCI_DEBUG 3045#ifdef LUCI_DEBUG 3046 print *,' h2_in:' 3047 do isym = 1, nnashx 3048 do jsym = 1, nnashx 3049 print *,' h2_in(',isym,jsym,') = ', 3050 & h2_in(isym,jsym) 3051 end do 3052 end do 3053 print *,' mult_ptg:' 3054 do isym = 1, nr_sym 3055 do jsym = 1, nr_sym 3056 print *,' mult_ptg(',isym,jsym,') = ', 3057 & mult_ptg(isym,jsym) 3058 end do 3059 end do 3060#endif 3061 3062! 2-electron integrals/2-particle density matrix elements 3063! ------------------------------------------------------- 3064 offset_internal = 1 3065 ioff = 1 3066! 4-index loop 3067 do isym = 1, nr_sym 3068 nr_act_isym = nr_orb_per_sym(isym) 3069 joff = 1 3070 do jsym = 1, isym 3071 nr_act_jsym = nr_orb_per_sym(jsym) 3072 koff = 1 3073 do ksym = 1, isym 3074 nr_act_ksym = nr_orb_per_sym(ksym) 3075 llsym = ksym 3076 3077 if(ksym .eq. isym) llsym = jsym 3078 3079 loff = 1 3080 do lsym = 1, llsym 3081 nr_act_lsym = nr_orb_per_sym(lsym) 3082 3083 ijsym = mult_ptg(isym,jsym) 3084! print *, ' isym, jsym',isym,jsym 3085 klsym = mult_ptg(ksym,lsym) 3086! print *, ' ksym, lsym',ksym,lsym 3087 ijklsym = mult_ptg(ijsym,klsym) 3088 3089! restrict loop to non-vanishing symmetry combinations of (ij|kl) 3090 3091 if(ijklsym .eq. 1)then 3092 3093#ifdef LUCI_DEBUG 3094 print *, ' check current symmetries' 3095 print('(/a,4i6)'), 'isym,jsym,ksym,lsym', 3096 & isym,jsym,ksym,lsym 3097#endif 3098 3099! index restrictions 3100 if(isym .eq. jsym)then 3101 ij_res = 1 3102 else 3103 ij_res = 0 3104 end if 3105 if(ksym .eq. lsym)then 3106 kl_res = 1 3107 else 3108 kl_res = 0 3109 end if 3110 3111! particle symmetry 3112 if(isym.eq.ksym.and.jsym.eq.lsym)then 3113 ijkl_res = 1 3114 else 3115 ijkl_res = 0 3116 end if 3117! 3118! K 3119! L 3120! I 3121! J 3122! loop over non-redundant output indices including symmetry offsets 3123 do k = koff, koff + nr_act_ksym - 1 3124 3125 if(kl_res .eq. 1)then 3126 l_last = k 3127 else 3128 l_last = loff + nr_act_lsym - 1 3129 end if 3130 3131 do l = loff, l_last 3132 3133 if(ijkl_res .eq. 1)then 3134 i_first = k 3135 else 3136 i_first = ioff 3137 end if 3138 3139 do i = i_first, ioff + nr_act_isym - 1 3140 3141 if(ij_res .eq. 1)then 3142 j_last = i 3143 else 3144 j_last = joff + nr_act_jsym - 1 3145 end if 3146 3147 if(ijkl_res .eq. 1 .and. i .eq. k)then 3148 j_first = l 3149 else 3150 j_first = joff 3151 end if 3152 3153 do j = j_first, j_last 3154 3155! indices on SIRIUS array 3156 ij = (i*(i-1))/2 + j 3157 kl = (k*(k-1))/2 + l 3158 3159! pick element 3160 h2_out(offset_internal) = h2_in(ij,kl) 3161 3162#ifdef LUCI_DEBUG 3163 print ('(/2x,a)'), 3164 & '-------------------------------------------------------' 3165 print ('(2x,a,2i4,a,2i4,a)'),'next integral:'// 3166 & ' (ij|kl) = (',i,j,'|',k,l,')' 3167 print *, 'ij, kl =',ij, kl 3168 print *, 'offset_internal =',offset_internal 3169 print *, 'h2_out(offset_internal) =', 3170 & h2_out(offset_internal) 3171 print('(2x,a/)'), 3172 & '-------------------------------------------------------' 3173#endif 3174 offset_internal = offset_internal + 1 3175 end do ! j 3176 end do ! i 3177 end do ! l 3178 end do ! k 3179 end if ! non-vanishing non-redundant (wrt symmetry) integrals 3180 loff = loff + nr_act_lsym 3181 end do ! lsym 3182 koff = koff + nr_act_ksym 3183 end do ! ksym 3184 joff = joff + nr_act_jsym 3185 end do ! jsym 3186 ioff = ioff + nr_act_isym 3187 end do ! isym 3188 3189! final # of 2-el ints 3190 nr_2int = offset_internal - 1 3191 3192!#define LUCI_DEBUG 3193#ifdef LUCI_DEBUG 3194 write(lupri,*) ' final 2e- list: # of ints =', nr_2int 3195 call wrtmatmn(h2_out,1,nr_2int,1,nr_2int,lupri) 3196#undef LUCI_DEBUG 3197#endif 3198 3199 end 3200*********************************************************************** 3201 3202 subroutine lu2dens2sir2dens(pv_sirius,pv_lu,pv_full_scratch,nasht, 3203 & ireost,ireots,nr_sym,nr_gas,rhotype) 3204! 3205! purpose: reorder 2e-RDM from dalton-sirius format to internal 3206! lucita format. 3207! 3208! 2e-RDM: on output symmetry-reduced list 3209! 3210! ----------------------------------------------------------------- 3211 implicit none 3212#include "priunit.h" 3213 integer, intent(in) :: nasht 3214 integer, intent(in) :: nr_sym 3215 integer, intent(in) :: nr_gas 3216 integer, intent(in) :: rhotype 3217 integer, intent(in) :: ireost(*) 3218 integer, intent(in) :: ireots(*) 3219 real(8), intent(out) :: pv_sirius((nasht*(nasht+1))/2, 3220 & (nasht*(nasht+1))/2) 3221 real(8), intent(in) :: pv_lu(*) 3222 real(8), intent(inout) :: pv_full_scratch(nasht,nasht,nasht,nasht) 3223! ----------------------------------------------------------------- 3224 integer :: nnashx 3225 integer :: n2ashx 3226 integer :: kl_ST_max, kl_ST_min, kl_ST_index 3227 integer :: l, k, k_ST, l_ST 3228 integer :: i, j, ijkl,klij, ij,kl 3229 real(8), allocatable :: temp_mat(:) 3230 logical :: do_reorder 3231 real(8) :: value 3232! ----------------------------------------------------------------- 3233 3234 nnashx = (nasht*(nasht+1))/2 3235 n2ashx = nasht**2 3236 3237! expand the lower-triangle 2-particle density matrix (lu format) to the full one 3238 call dsptsi(n2ashx,pv_lu,pv_full_scratch) 3239 3240!#define LUCI_DEBUG 3241#ifdef LUCI_DEBUG 3242 write(lupri,*) ' final rho2 (full)' 3243 call wrtmatmn(pv_full_scratch,n2ashx,n2ashx,n2ashx,n2ashx,lupri) 3244#endif 3245 3246 if(rhotype > 1)then 3247 call symmetrize_transition_density_matrix(pv_full_scratch,nasht,1) 3248 end if 3249 3250#ifdef LUCI_DEBUG 3251 write(lupri,*) ' final rho2 (full) - symm attempt' 3252 call wrtmatmn(pv_full_scratch,n2ashx,n2ashx,n2ashx,n2ashx,lupri) 3253 do i = 1, nasht 3254 do j = 1, nasht 3255 do k = 1, nasht 3256 do l = 1, nasht 3257 write(lupri,'(a,i2,a1,i2,a1,i2,a1,i2,a,f16.10)') 3258 & ' pv(',i,',',j,',',k,',',l,') = ', 3259 & pv_full_scratch(i,j,k,l) 3260 end do 3261 end do 3262 end do 3263 end do 3264#endif 3265 3266 do_reorder = (nr_gas .gt. 1 .and. nr_sym .gt. 1) 3267 if (do_reorder) allocate(temp_mat(n2ashx)) 3268 3269 do l = 1, nasht 3270 l_ST = ireots(l) 3271 do k = 1, l 3272 k_ST = ireots(k) 3273 kl_ST_max = max(k_ST,l_ST) 3274 kl_ST_min = min(k_ST,l_ST) 3275 kl_ST_index = (kl_ST_max*(kl_ST_max-1))/2 + kl_ST_min 3276 if(do_reorder)then 3277! reorder from type-symmetry ordering to a symmetry-type one 3278 call reormt(pv_full_scratch(1,1,k,l),temp_mat, 3279 & nasht,nasht,ireost,ireost) 3280! distribute elements on lower triangle 3281 call dgetsp(nasht,temp_mat,pv_sirius(1,kl_ST_index)) 3282 else 3283! distribute elements on lower triangle 3284 call dgetsp(nasht,pv_full_scratch(1,1,k,l), 3285 & pv_sirius(1,kl_ST_index)) 3286 end if 3287 end do 3288 end do 3289 3290 3291 if(do_reorder) deallocate(temp_mat) 3292 3293!#define LUCI_DEBUG 3294#ifdef LUCI_DEBUG 3295 write(lupri,*) ' final PV in SIRIUS format' 3296 call output(pv_sirius,1,nnashx,1,nnashx,nnashx,nnashx,-1,lupri) 3297#undef LUCI_DEBUG 3298#endif 3299 3300 end 3301*********************************************************************** 3302 3303 subroutine set_ptg_multiplication_table(mult_ptg,nr_ptg_irreps) 3304! 3305! purpose: set point-group multiplication table for 3306! dalton -> lucita integral transfer. 3307! 3308! ----------------------------------------------------------------- 3309 implicit none 3310#include "multd2h.inc" 3311 integer :: nr_ptg_irreps, i, j 3312 integer :: mult_ptg(nr_ptg_irreps,nr_ptg_irreps) 3313! ----------------------------------------------------------------- 3314 3315 do j = 1, nr_ptg_irreps 3316 do i = 1, nr_ptg_irreps 3317 mult_ptg(i,j) = multd2h(i,j) 3318 end do 3319 end do 3320 3321 end 3322*********************************************************************** 3323 3324 subroutine write_fcidump_file(AMATRIX,nrow,BMATRIX,corenergy) 3325! ----------------------------------------------------------------- 3326! Written by Katharina Boguslawski and Pawel Tecmer 3327! ETH Zurich, March 2013 3328! 3329! purpose: write one- and two-electron integrals for a given active 3330! space to disk 3331! 3332! filename: FCIDUMP 3333! ----------------------------------------------------------------- 3334 3335 use lucita_cfg 3336 3337#include "implicit.h" 3338#include "maxorb.h" 3339#include "infinp.h" 3340#include "priunit.h" 3341#include "inforb.h" 3342#include "inftap.h" 3343#include "clunit.inc" 3344!#include "orbinp.inc" 3345 3346 character(len=30) :: form1 3347 character(len=30) :: form2 3348 character(len=30) :: form3 3349 DIMENSION AMATRIX(*) 3350 DIMENSION BMATRIX(*) 3351 INTEGER NDUMMY, start 3352 INTEGER NTRIANGLE, norbtot, noccend, noccendi, noccendj 3353 integer, allocatable :: mult_ptg(:,:) 3354 integer nactiveorb(NSYM) 3355 integer :: offseti, offsetj, offsetk, offsetl 3356 real*8 :: corenergy 3357! ----------------------------------------------------------------- 3358 integer :: ISYM, KSYM, syml, JSYM, ijsym, klsym, ijklsym 3359 integer :: i, j, k, l 3360! ----------------------------------------------------------------- 3361 3362 COMMON/CECORE/ECORE,ECORE_ORIG,ECORE_H,ECORE_HEX 3363 3364! allocate memory for multiplication table to write only 3365! symmetry-allowed integrals 3366 allocate(mult_ptg(NSYM,NSYM)) 3367 mult_ptg = 0 3368 call set_ptg_multiplication_table(mult_ptg,NSYM) 3369 3370 LUFCID = -1 3371 3372! define printing format 3373 form1="(F17.12, 4X, I6, I6, I6, I6)" 3374 form2="(A11,I3,A7,I2,A5,I2,A1)" 3375 form3="(F17.8, 4X, I6, I6, I6, I6)" 3376 3377 3378! FCIDUMP info 3379 write(lupri,*) 3380 write(lupri,*)" FCIDUMP file details " 3381 write(lupri,*)"======================" 3382 write(lupri,*) 3383 3384 3385! calculate number of active electrons 3386 norbtot = 0 3387 do KSYM=1,NSYM,1 3388 norbtot = norbtot + NOCC(KSYM) - NISH(KSYM) 3389 nactiveorb(KSYM) = NASH(KSYM) 3390 end do 3391 3392! Print header of FCIDUMP file 3393! using MOLPRO format 3394 if (LUFCID .le. 0) then 3395 CALL GPOPEN(LUFCID,'FCIDUMP','new',' ', 3396 & 'FORMATTED',IDUMMY,.false.) 3397 end if 3398 write(LUFCID,form2) ' &FCI NORB=', norbtot , 3399 & ',NELEC=', lucita_cfg_nr_active_e, ',MS2=', 3400 & lucita_cfg_is_spin_multiplett-1, ',' 3401 write(LUFCID,"(A)",advance='no') ' ORBSYM=' 3402 do ISYM=1,NSYM,1 3403 if (nactiveorb(ISYM).ne.0) then 3404 do i=1,nactiveorb(ISYM),1 3405 write(LUFCID,"(I1,A1)",advance='no') ISYM,',' 3406 end do 3407 end if 3408 end do 3409 3410 3411 write(lupri,form2) ' &FCI NORB=', norbtot , 3412 & ',NELEC=', lucita_cfg_nr_active_e, ',MS2=', 3413 & lucita_cfg_is_spin_multiplett-1, ',' 3414 3415 3416 write(LUFCID,*) 3417 write(LUFCID,"(A7,I1)") ' ISYM=',lucita_cfg_ptg_symmetry 3418 write(LUFCID,"(A5)") ' &END' 3419 3420! NOW: 3421! Print two-electron integrals: 3422! integrals are sorted as KL|IJ 3423! order is different as in MOLPRO, where integrals are sorted as 3424! IJ|KL 3425! Integrals are stored in a set of triangular matrices, define index 3426! which accesses matrix elements 3427 IELE = 0 3428! Determine number of elements in one triangular matrix: 3429 NTRIANGLE = norbtot*(norbtot+1)/2 3430 3431! KL IJ 3432! define orbital offset for K 3433 offsetk = 0 3434! sum over all irreducible representations 3435 do KSYM=1,NSYM,1 3436 if (nactiveorb(KSYM).ne.0) then 3437 do k=1,nactiveorb(KSYM),1 ! iterate through all elements 3438 offsetl = 0 !define orbital offset for L 3439! restrict summation to prevent double counting 3440 do syml=1,KSYM,1 3441 if (nactiveorb(syml).ne.0) then 3442! set upper summation bound for orbital index (prevent double 3443! counting): 3444! if not the same irrep l goes from 1 to number of orbitals 3445 if (KSYM.eq.syml) then 3446 noccend = k 3447 else 3448 noccend = nactiveorb(syml) 3449 end if 3450 do l=1,noccend,1 3451! orbital offset for I 3452 offseti = 0 3453! restrict summation to prevent double counting for both irrep ISYM and 3454! orbital indices i 3455 do ISYM=1,KSYM,1 3456 if (nactiveorb(ISYM).ne.0) then 3457 ijksyml = 0 3458 if (ISYM.eq.KSYM) then 3459 noccendi = k 3460 else 3461 noccendi = nactiveorb(ISYM) 3462 end if 3463 do i=1,noccendi,1 3464! orbital offset J 3465 offsetj = 0 3466! double counting problem: irrep of J must be smaller or equal to 3467! irrep of I 3468 do JSYM=1,ISYM,1 3469! Collect only integrals which are allowed by symmetry, all others 3470! are neglected. 3471! Two cases have to be distinguished: IJ|KL and IK|JL 3472 if (nactiveorb(JSYM).ne.0) then 3473 if(ISYM.eq.JSYM.and.KSYM.eq.syml) then ! first case 3474 ijsym = mult_ptg(ISYM,JSYM) 3475 ksyml = mult_ptg(KSYM,syml) 3476 ijksyml = mult_ptg(ijsym,ksyml) 3477 else ! second case 3478 ijsym = mult_ptg(ISYM,KSYM) 3479 ksyml = mult_ptg(JSYM,syml) 3480 ijksyml = mult_ptg(ijsym,ksyml) 3481 3482 end if 3483 3484 if(ijksyml .eq. 1) then 3485 3486! Set matric index: first determine first index of block KL 3487 IELE = (offsetk+k)*(offsetk+k-1)/2*NTRIANGLE 3488 & +(offsetl+l-1)*NTRIANGLE+1 3489! Prevent double printing of symmetry redundant indices: 3490! set upper summation index 3491! if IJKL in same irrep, restrict j to at most i 3492 if (JSYM.eq.ISYM.and.KSYM.eq.syml) then !.and.ISYM.eq.KSYM) then 3493 noccendj = i 3494!redundant else if (JSYM.eq.ISYM.and.KSYM.eq.syml.and.ISYM.ne.KSYM) then 3495! noccendj = i 3496! if LJ in irrep1 and IK in irrep2 3497 else if (syml.eq.JSYM.and.KSYM.eq.ISYM 3498 & .and.syml.ne.ISYM) then 3499! second restriction to prevent double printing, J<=L in KL IJ 3500 if (k.eq.i) then 3501 noccendj = l 3502 else ! otherwise all J are needed 3503 noccendj = nactiveorb(JSYM) 3504 end if 3505 else 3506 noccendj = nactiveorb(JSYM) 3507 end if 3508 IELE = IELE+(offseti+i)*(offseti+i-1)/2+offsetj 3509 do j=1,noccendj,1 3510! check for redundant integrals 3511 if (JSYM.eq.ISYM.and.KSYM.eq.syml 3512 & .and.ISYM.eq.KSYM) then 3513 if (k.eq.i.and.l.eq.i.and.j.lt.i) then 3514 3515 else if (k.eq.i.and.j.lt.l) then 3516 3517 else 3518 write(LUFCID,form1) BMATRIX(IELE), 3519 & i+offseti, j+offsetj, k+offsetk, l+offsetl 3520 end if 3521 else 3522 write(LUFCID,form1) BMATRIX(IELE), 3523 & i+offseti, j+offsetj, k+offsetk, l+offsetl 3524 end if 3525 IELE = IELE + 1 3526 end do ! do j 3527 offsetj = offsetj + nactiveorb(JSYM) 3528 3529 else ! if not symmetry-allowed go to next irrep and update 3530! orbital offset J 3531 offsetj = offsetj + nactiveorb(JSYM) 3532 end if ! if ijksyml 3533 end if ! if nocc(jsym) 3534 end do ! do jsym 3535 3536 end do ! do i 3537 offseti = offseti + nactiveorb(ISYM) !update orbital 3538! offset I 3539 end if ! if nocc(isym) 3540 end do ! do isym 3541 end do ! do l 3542 offsetl = offsetl + nactiveorb(syml) 3543 end if ! if nocc(syml) 3544 end do ! do syml 3545 end do ! do k 3546 end if ! if nocc(ksym) 3547 offsetk = offsetk + nactiveorb(KSYM) 3548 end do ! do ksym 3549 3550 3551! Print one-electron integrals: 3552! Reset matrix index 3553 IELE = 1 3554! number of dummy indices to be ignored in one-electron 3555! integrals because of symmetry ISYM < actual ISYM 3556 NDUMMY = 0 3557 3558 do ISYM=1,NSYM,1 3559 if (nactiveorb(ISYM).ne.0) then 3560 IELE = IELE + NDUMMY 3561 do i=1,nactiveorb(ISYM),1 ! only loop through same irrep 3562 do j=1,i,1 3563 write(LUFCID,form1) AMATRIX(IELE), i+NDUMMY, j+NDUMMY, 3564 & 0, 0 3565 IELE = IELE + 1 3566 end do 3567! add orbital offset if irrep changes 3568 if (i.lt.nactiveorb(ISYM)) then 3569 IELE = IELE + NDUMMY 3570 end if 3571 end do 3572! update orbital offset 3573 NDUMMY = NDUMMY + nactiveorb(ISYM) 3574 end if 3575 end do 3576 3577 write(LUFCID,form3) POTNUC + corenergy, 0,0 ,0,0 3578 3579 call GPCLOSE(LUFCID,'KEEP') 3580 3581 deallocate(mult_ptg) 3582 3583 return 3584 3585 end 3586