1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck ccsd_sortao */ 20 SUBROUTINE CCSD_SORTAO(WORK,LWORK) 21C 22C Written by Henrik Koch 25-Sep-1993 23C 24C Purpose: Read destribution of AO integrals. 25C 26 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 27#include "priunit.h" 28#include "iratdef.h" 29#include "maxorb.h" 30C 31 LOGICAL FIRST, ENABLED 32 DIMENSION WORK(LWORK) 33C 34 CHARACTER*8 NAME(8) 35 CHARACTER*8 LBLSAV 36C 37C SAVE FIRST 38C 39 DATA FIRST /.TRUE./ 40 DATA NAME /'CCAOIN_1','CCAOIN_2','CCAOIN_3','CCAOIN_4', 41 * 'CCAOIN_5','CCAOIN_6','CCAOIN_7','CCAOIN_8'/ 42 COMMON/SORTIO/LUAOIN(8) 43#include "inftap.h" 44#include "ccorb.h" 45#include "ccsdsym.h" 46#include "ccfop.h" 47#include "ccsdinp.h" 48#include "ccpack.h" 49#include "r12int.h" 50C 51C-------------------- 52C Only sort once. 53C-------------------- 54C 55C IF (FIRST) THEN 56C FIRST = .FALSE. 57C ELSE 58C RETURN 59C ENDIF 60C 61C------------------------- 62C Open integral files. 63C------------------------- 64C 65 IF (LUINTA.LE.0) THEN 66 LUINTA = -1 67 CALL MAKE_AOTWOINT(WORK,LWORK) 68 CALL GPOPEN(LUINTA,'AOTWOINT','UNKNOWN',' ','UNFORMATTED', 69 & IDUMMY,.FALSE.) 70 END IF 71C 72 DO 50 ISYM = 1,NSYM 73C 74 NFILE = 0 75 CALL WOPEN2(NFILE,NAME(ISYM),64,0) 76 LUAOIN(ISYM) = NFILE 77C 78 50 CONTINUE 79C 80C------------------------------ 81C Skip sorting if required. 82C------------------------------ 83C 84 IF (NOSORT) THEN 85C 86 DO ISYMD = 1,NSYM 87 ISYDIS = MULD2H(ISYMD,ISYMOP) 88 LENGTH = NDISAO(ISYDIS) 89 IOFFID = 1 90 DO ID = 1,NBAS(ISYMD) 91 IDEL = ID + IBAS(ISYMD) 92 IOFFINT(IDEL) = IOFFID 93 IOFFID = IOFFID + LENGTH 94 END DO 95 END DO 96C 97 RETURN 98C 99 END IF 100C 101C------------------------ 102C Buffer information. 103C------------------------ 104C 105 IF (CCR12) THEN 106 NALLBAS = 0 107 DO I = 1, NSYM 108 NALLBAS = NALLBAS + MBAS1(I) + MBAS2(I) 109 END DO 110 ELSE 111 NALLBAS = NBAST 112 ENDIF 113C 114 LBUF = 600 115 IF (NALLBAS .LE. 255) THEN 116 NIBUF = 1 117 NBITS = 8 118 IBIT1 = 2**8 - 1 119 IBIT2 = 2**16 - 1 120 ELSE 121 NIBUF = 2 122 NBITS = 16 123 IBIT1 = 2**16 - 1 124 IBIT2 = 0 ! not used when NIBUF .eq. 2 125 END IF 126C 127C----------------------- 128C Buffer allocation. 129C----------------------- 130C 131 KRBUF = 1 132 KIBUF = KRBUF + LBUF 133 KAOAB = KIBUF + (NIBUF*LBUF + 1)/2 + 1 ! IBUF always integer*4 134 KAOG = KAOAB + (N2BASX + 1)/IRAT + 1 135 KEND1 = KAOG + (NBAST*NSYM + 1)/IRAT + 1 136 LWRK1 = LWORK - KEND1 137C 138 IF (LWRK1 .LT. 0) CALL QUIT('Insufficient work space in CCRDAO') 139C 140C------------------------------------------------------ 141C Calculate in the index arrays needed in the sort. 142C------------------------------------------------------ 143C 144 CALL CCSD_INIT2(WORK(KAOAB),WORK(KAOG)) 145C 146C------------------------------------------- 147C set up table for packing of integrals: 148C------------------------------------------- 149C 150 IF (LPACKINT) THEN 151 DTIME = SECOND() 152 CALL INITPCKR8(THRPCKINT,IPCKTABINT,ENABLED) 153 154 IF (.NOT.ENABLED) THEN 155 WRITE(LUPRI,'(A)') 156 & 'packing routines not enabled for this architecture...', 157 & '...the integral packing is switched off...' 158 LPACKINT = .FALSE. 159 END IF 160 161 NTOTPCK = 0 162 NTOTINT = 0 163 PCKRATIO = 1.0D0 164 PCKTIME = SECOND() - DTIME 165 END IF 166C 167C------------------------------------ 168C Loop over batches of integrals. 169C------------------------------------ 170C 171C 172 DO 100 ISYMD = 1,NSYM 173C 174 IOFF2 = 1 175C 176 ISYDIS = MULD2H(ISYMD,ISYMOP) 177 NTOTD = NBAS(ISYMD) 178 IF (NTOTD .EQ. 0) GOTO 100 179 LENGTH = NDISAO(ISYDIS) 180C 181 NUMBAT = MIN(NTOTD,LWRK1/LENGTH) 182C 183 IF (NUMBAT .EQ. 0) THEN 184 WRITE(LUPRI,*) 'In CCSD_SORTAO NUMBAT is zero' 185 CALL QUIT('Insufficient work space in CCRDAO') 186 ENDIF 187C 188 ITOTBA = (NTOTD-1)/NUMBAT + 1 189C 190 ID1 = IBAS(ISYMD) + 1 191 ID2 = IBAS(ISYMD) 192 IOFF1 = IBAS(ISYMD) 193C 194 DO 200 I = 1,ITOTBA 195C 196 INUMBA = NUMBAT 197 IF (NUMBAT*I .GT. NTOTD) THEN 198 INUMBA = NTOTD - NUMBAT*(I-1) 199 ENDIF 200C 201 ID2 = ID2 + INUMBA 202C 203 CALL DZERO(WORK(KEND1),LENGTH*INUMBA) 204C 205 CALL CCSD_SORT1(LUINTA,WORK(KEND1),WORK(KIBUF),WORK(KRBUF), 206 * WORK(KAOAB),WORK(KAOG),ISYDIS,LENGTH, 207 * IOFF1,ID1,ID2,NIBUF,LBUF,NBITS,IBIT1) 208C 209 CALL CCSD_SORT2(WORK(KEND1),IOFF2,INUMBA,LENGTH, 210 * ID1,ID2,ISYMD) 211C 212 IF (IPRINT .GT. 50) THEN 213 CALL AROUND('Integral distribution') 214 IPRC = KEND1 215 DO 210 IPRD = ID1,ID2 216 WRITE(LUPRI,*) 'D distribution',IPRD 217 DO 220 IPSYMG = 1,NSYM 218 WRITE(LUPRI,*) 'Gamma symmetry',IPSYMG 219 ISYMAB = MULD2H(IPSYMG,ISYDIS) 220 CALL OUTPUT(WORK(IPRC),1,NNBST(ISYMAB),1, 221 * NBAS(IPSYMG),NNBST(ISYMAB), 222 * NBAS(IPSYMG),1,LUPRI) 223 IPRC = IPRC + NNBST(ISYMAB)*NBAS(IPSYMG) 224 220 CONTINUE 225 210 CONTINUE 226 END IF 227C 228 ID1 = ID1 + INUMBA 229 IOFF1 = IOFF1 + INUMBA 230C 231 200 CONTINUE 232C 233 100 CONTINUE 234C 235C------------------------------------- 236C Print packing statistics: 237C------------------------------------- 238C 239 IF (IPRINT.GT.0 .AND. LPACKINT) THEN 240 WRITE (LUPRI,'(//10X,A,F9.2,A)') 241 & 'Time needed to pack integrals: ', 242 & PCKTIME, ' seconds' 243 WRITE (LUPRI,'(10X,A,G9.2)') 244 & 'Threshold used for packing: ', 245 & THRPCKINT 246 NTOTINT = MAX(NTOTINT,1) 247 PCKRATIO = DBLE(NTOTPCK)/DBLE(NTOTINT) 248 WRITE (LUPRI,'(10X,A,F9.2,A)') 249 & 'Reduction obtained by packing: ', 250 & 100.0D0*(1.0D0 - PCKRATIO),' %' 251 END IF 252C 253C------------------------------------- 254C Close integral files and delete. 255C------------------------------------- 256C 257 IF (KEEPAOTWO .EQ. 0) THEN 258 CALL GPCLOSE(LUINTA,'DELETE') 259 ELSE 260 CALL GPCLOSE(LUINTA,'KEEP') 261 ENDIF 262C 263 RETURN 264 END 265C /* Deck ccsd_sort1 */ 266 SUBROUTINE CCSD_SORT1(LUINTA,XINT,IBUF4,RBUF,KAOAB,KAOG,ISYDIS, 267 * LENGTH,IOFF,ID1,ID2,NIBUF,LBUF,NBITS,IBIT1) 268C 269C Written by Henrik Koch 25-Sep-1993 270C 271#include "implicit.h" 272#include "priunit.h" 273#include "ibtpar.h" 274#include "ccorb.h" 275 REAL*8 XINT(LENGTH,*),RBUF(LBUF) 276 INTEGER*4 IBUF4(NIBUF*LBUF), LENGTH4 277 INTEGER KAOAB(NBAST,NBAST),KAOG(NBAST,NSYM) 278#include "ccsdsym.h" 279#include "r12int.h" 280 281C 282C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 283C 284 REWIND (LUINTA) 285 CALL MOLLAB('BASTWOEL',LUINTA,LUPRI) 286C 287 IF (NIBUF .EQ. 1) THEN 288C 289 10 READ(LUINTA,ERR=2000) RBUF,IBUF4,LENGTH4 290C 291 IF (LENGTH4 .EQ. -1) GOTO 100 292C 293 DO I = 1,LENGTH4 294C 295 LABLE = IBUF4(I) 296 VALUE = RBUF(I) 297C 298 IP = IAND( LABLE ,IBIT1) 299 IQ = IAND(ISHFT(LABLE, -NBITS),IBIT1) 300 IR = IAND(ISHFT(LABLE,-2*NBITS),IBIT1) 301 IS = IAND(ISHFT(LABLE,-3*NBITS),IBIT1) 302 IF (NOAUXB) CALL IJKAUX(IP,IQ,IR,IS) 303C 304 IF ((IS .GE. ID1) .AND. (IS .LE. ID2)) THEN 305 IADR = KAOG(IR,ISYDIS) + KAOAB(IP,IQ) 306 XINT(IADR,IS-IOFF) = VALUE 307 ENDIF 308 IF ((IR .GE. ID1) .AND. (IR .LE. ID2)) THEN 309 IADR = KAOG(IS,ISYDIS) + KAOAB(IP,IQ) 310 XINT(IADR,IR-IOFF) = VALUE 311 ENDIF 312 IF ((IP .GE. ID1) .AND. (IP .LE. ID2)) THEN 313 IADR = KAOG(IQ,ISYDIS) + KAOAB(IR,IS) 314 XINT(IADR,IP-IOFF) = VALUE 315 ENDIF 316 IF ((IQ .GE. ID1) .AND. (IQ .LE. ID2)) THEN 317 IADR = KAOG(IP,ISYDIS) + KAOAB(IR,IS) 318 XINT(IADR,IQ-IOFF) = VALUE 319 ENDIF 320C 321 END DO 322C 323 GOTO 10 324 100 CONTINUE 325C 326 ELSE 327C 328 30 READ(LUINTA,ERR=2000) RBUF,IBUF4,LENGTH4 329C 330 IF (LENGTH4 .EQ. -1) GOTO 200 331C 332 DO 40 I = 1,LENGTH4 333C 334 LABLE1 = IBUF4(2*I-1) 335 LABLE2 = IBUF4(2*I ) 336 VALUE = RBUF(I) 337C 338 IP = IAND( LABLE1 ,IBIT1) 339 IQ = IAND(ISHFT(LABLE1,-NBITS),IBIT1) 340 IR = IAND( LABLE2, IBIT1) 341 IS = IAND(ISHFT(LABLE2,-NBITS),IBIT1) 342 IF (NOAUXB) CALL IJKAUX(IP,IQ,IR,IS) 343C 344 IF ((IS .GE. ID1) .AND. (IS .LE. ID2)) THEN 345 IADR = KAOG(IR,ISYDIS) + KAOAB(IP,IQ) 346 XINT(IADR,IS-IOFF) = VALUE 347 ENDIF 348 IF ((IR .GE. ID1) .AND. (IR .LE. ID2)) THEN 349 IADR = KAOG(IS,ISYDIS) + KAOAB(IP,IQ) 350 XINT(IADR,IR-IOFF) = VALUE 351 ENDIF 352 IF ((IP .GE. ID1) .AND. (IP .LE. ID2)) THEN 353 IADR = KAOG(IQ,ISYDIS) + KAOAB(IR,IS) 354 XINT(IADR,IP-IOFF) = VALUE 355 ENDIF 356 IF ((IQ .GE. ID1) .AND. (IQ .LE. ID2)) THEN 357 IADR = KAOG(IP,ISYDIS) + KAOAB(IR,IS) 358 XINT(IADR,IQ-IOFF) = VALUE 359 ENDIF 360C 361 40 CONTINUE 362C 363 GOTO 30 364 200 CONTINUE 365C 366 ENDIF 367C 368 RETURN 369 2000 CALL QUIT('Error reading AOTWOINT in CCSD_SORT1') 370 END 371C /* Deck ccsd_sort2 */ 372 SUBROUTINE CCSD_SORT2(XINT,IOFF,INUMBA,LENGTH,ID1,ID2,ISYM) 373C 374C Written by Henrik Koch 25-Sep-1993 375C 376#include "implicit.h" 377#include "priunit.h" 378 DIMENSION XINT(*) 379C 380 CHARACTER*8 NAME(8) 381C 382 DATA NAME /'CCAOIN_1','CCAOIN_2','CCAOIN_3','CCAOIN_4', 383 * 'CCAOIN_5','CCAOIN_6','CCAOIN_7','CCAOIN_8'/ 384 385 COMMON/SORTIO/LUAOIN(8) 386#include "ccorb.h" 387#include "maxorb.h" 388#include "ccpack.h" 389C 390 NFILE = LUAOIN(ISYM) 391 KOFF = 1 392C 393 DO IDEL = ID1, ID2 394C 395 NBYTE = LENGTH*8 396 NDWORDS = LENGTH 397C 398 IF (LPACKINT) THEN 399 400 DTIME = SECOND() 401 402 CALL PCKR8(XINT(KOFF),LENGTH,XINT(KOFF),NBYTE, 403 & IPCKTABINT,LPACKINT) 404 405 NDWORDS = (NBYTE+7)/8 406 NTOTINT = NTOTINT + LENGTH / NBAST 407 NTOTPCK = NTOTPCK + NDWORDS / NBAST 408 PCKTIME = PCKTIME + SECOND() - DTIME 409 410 END IF 411C 412 CALL PUTWA2(NFILE,NAME(ISYM),XINT(KOFF),IOFF,NDWORDS) 413C 414 IOFFINT(IDEL) = IOFF 415 NPCKINT(IDEL) = NBYTE 416C 417 IOFF = IOFF + NDWORDS 418 KOFF = KOFF + LENGTH 419C 420 END DO 421C 422 RETURN 423 END 424C /* Deck ccsd_init2 */ 425 SUBROUTINE CCSD_INIT2(KAOAB,KAOG) 426C 427C Henrik Koch and Alfredo Sanchez. 29-Jun-1994 428C 429C Set up indexing arrays 430C 431#include "implicit.h" 432#include "priunit.h" 433#include "ccorb.h" 434 DIMENSION KAOAB(NBAST,NBAST),KAOG(NBAST,NSYM) 435#include "ccsdsym.h" 436C 437C 438 DO 100 ISYMAB = 1,NSYM 439 ICOUNT = 0 440 DO 110 ISYMB = 1,NSYM 441 ISYMA = MULD2H(ISYMB,ISYMAB) 442 IF (ISYMB .GT. ISYMA) THEN 443 DO 120 B = 1,NBAS(ISYMB) 444 IB = IBAS(ISYMB) + B 445 DO 130 A = 1,NBAS(ISYMA) 446 IA = IBAS(ISYMA) + A 447 ICOUNT = ICOUNT + 1 448 KAOAB(IA,IB) = ICOUNT 449 KAOAB(IB,IA) = ICOUNT 450 130 CONTINUE 451 120 CONTINUE 452 ELSE IF (ISYMA .EQ. ISYMB) THEN 453 DO 140 B = 1,NBAS(ISYMB) 454 IB = IBAS(ISYMB) + B 455 DO 150 A = 1,B 456 IA = IBAS(ISYMA) + A 457 ICOUNT = ICOUNT + 1 458 KAOAB(IA,IB) = ICOUNT 459 KAOAB(IB,IA) = ICOUNT 460 150 CONTINUE 461 140 CONTINUE 462 END IF 463 110 CONTINUE 464 100 CONTINUE 465C 466 DO 200 ISYMD = 1,NSYM 467 ISYDIS = MULD2H(ISYMD,ISYMOP) 468 ICOUNT = 0 469 DO 210 ISYMG = 1,NSYM 470 ISYMAB = MULD2H(ISYMG,ISYDIS) 471 DO 220 G = 1,NBAS(ISYMG) 472 IG = IBAS(ISYMG) + G 473 KAOG(IG,ISYMD) = ICOUNT 474 ICOUNT = ICOUNT + NNBST(ISYMAB) 475 220 CONTINUE 476 210 CONTINUE 477 200 CONTINUE 478C 479 RETURN 480 END 481