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! 18module mlcc3_intermediates 19! 20! 21! mlcc3 routines to calculate fock matrix and MO integrals 22! Authors Henrik Koch and Rolf H. Myhre 23! December 2014 24! 25 use mlcc_typedef 26 use mlcc_block_import 27 use mlcc_work 28 use mlcc3_active_spaces 29 use mlcc3_data 30 use mlcc3_reordering 31 use mlcc3_various 32! 33! AO integrals 34 real(dp), dimension(:), pointer, private :: ao_integrals_pack => null() 35! 36! Unsorted virtual MO integrals 37 real(dp), dimension(:), pointer, private :: b_D_c_k => null() 38 real(dp), dimension(:), pointer, private :: D_b_k_c => null() 39! 40! Unsorted occupied MO integrals 41 real(dp), dimension(:), pointer, private :: L_j_c_k => null() 42 real(dp), dimension(:), pointer, private :: j_L_k_c => null() 43! 44! Unsorted active MO integrals 45 real(dp), dimension(:), pointer, private :: j_b_k_c => null() 46! 47! Sorted virtual MO integrals 48 real(dp), dimension(:), pointer, private :: integral_sort => null() 49! 50! First transformed array 51 real(dp), dimension(:), pointer, private :: alpha_beta_k_pack_hole => null() 52 real(dp), dimension(:), pointer, private :: alpha_beta_k_pack_part => null() 53! 54! n_basis x n_basis array, used many places 55 real(dp), dimension(:), pointer, private :: alpha_beta => null() 56! 57! n_virtual_active x n_basis 58 real(dp), dimension(:), pointer, private :: b_beta => null() 59! 60! n_virtual_active x n_virtual_general 61 real(dp), dimension(:), pointer, private :: b_D => null() 62 real(dp), dimension(:), pointer, private :: D_b => null() 63! 64! n_occupied_active x n_occupied_general 65 real(dp), dimension(:), pointer, private :: j_L => null() 66 real(dp), dimension(:), pointer, private :: L_j => null() 67! 68! n_occupied_active x n_virtual_active 69 real(dp), dimension(:), pointer, private :: j_b => null() 70! 71contains 72! 73 subroutine mlcc3_intermediates_calc(resp) 74! 75! intermediates calculator 76! Authors Henrik Koch and Rolf H. Myhre 77! December 2014 78! 79! Purpose: loop over AO's, read in AO integrals and calculate 80! Fock matrix and MO integrals 81! 82 implicit none 83! 84 logical, intent(in) :: resp 85! 86 integer :: delta 87 integer :: batch, n_batches, batch_size, memory_use, batch_size2 88 logical :: first 89 integer :: integral_size_vir, integral_size_occ 90 integer :: ioerror 91 integer :: D,b,c,k,j,L, k_off 92 integer :: batch_start, batch_end, n_occ_print, n_occ_int 93 integer :: rec_number, n_act_int 94 integer :: work_remains 95! 96 integer :: bDck_print_unit, Dbkc_print_unit 97 integer :: Ljck_print_unit, jLkc_print_unit 98 integer :: jbkc_print_unit 99! 100 real(dp) :: ddot 101 real(dp) :: bDcknorm, Dbkcnorm,Ljcknorm,jLkcnorm,jbkcnorm 102! 103! File names 104 character(len=12) :: bDck_name 105 character(len=12) :: Dbkc_name 106 character(len=12) :: Ljck_name 107 character(len=12) :: jLkc_name 108 character(len=12) :: jbkc_name 109! 110! 111! Timing variables 112 real :: time_tot, time_start, time_1, time_2 113 real :: time_fock, time_int, time_trans, time_open 114 real :: time_read, time_write, time_close 115 real :: time_zero, time_reord 116 real :: time_bDck, time_Dbkc, time_Ljck, time_jLkc, time_jbkc 117! 118 time_fock = 0.0 119 time_int = 0.0 120 time_write = 0.0 121 time_open = 0.0 122 time_read = 0.0 123 time_trans = 0.0 124 time_close = 0.0 125 time_zero = 0.0 126 time_bDck = 0.0 127 time_Dbkc = 0.0 128 time_Ljck = 0.0 129 time_jLkc = 0.0 130 time_jbkc = 0.0 131 time_reord = 0.0 132! 1339999 format(7x,'Time used in',2x,A22,2x,': ',f14.4,' seconds') 134! 135 call cpu_time(time_start) 136! 137! Energy or response integrals 138! ---------------------------- 139! 140 bDcknorm = zero 141 Dbkcnorm = zero 142 Ljcknorm = zero 143 jLkcnorm = zero 144 jbkcnorm = zero 145! 146 if(resp) then 147! 148 bDck_name = bDck_resp_name 149 Dbkc_name = Dbkc_resp_name 150 Ljck_name = Ljck_resp_name 151 jLkc_name = jLkc_resp_name 152! 153 else 154! 155 bDck_name = bDck_file_name 156 Dbkc_name = Dbkc_file_name 157 Ljck_name = Ljck_file_name 158 jLkc_name = jLkc_file_name 159 jbkc_name = jbkc_file_name 160! 161 end if 162! 163! Some printing variables 164! ----------------------- 165 n_act_int = n_vir_act*n_vir_act*n_occ_act 166 n_occ_int = n_occ_gen*n_occ_act*n_vir_act 167! 168 first = .true. 169! 170! Allocate space for packed AO integrals 171! 172 call work_allocator(ao_integrals_pack,n_ao_ints) 173! 174! Alllocate arrays independent of batch size 175! 176 call work_allocator(alpha_beta,n_basis_2) 177 call work_allocator(b_beta,n_vir_act*n_basis) 178! 179 call work_allocator(b_D,n_vir_act*n_vir_gen) 180 call work_allocator(D_b,n_vir_act*n_vir_gen) 181! 182 call work_allocator(L_j,n_occ_act*n_occ_gen) 183 call work_allocator(j_L,n_occ_act*n_occ_gen) 184! 185 call work_allocator(j_b,n_occ_act*n_vir_act) 186! 187 call work_allocator(integral_sort,n_vir_aag) 188! 189! 190! Figure out how many k's we can batch over 191! 192 memory_use = 2*(n_vir_aag + n_occ_int + n_basis_2_pack) + n_act_int 193! 194 work_remains = work_free() 195! 196 if (memory_use .gt. work_remains) then 197 call quit('Not enough memory for a batch in mlcc3_intermediates') 198 end if 199! 200 batch_size = min(work_remains/memory_use,n_occ_act) 201 n_batches = (n_occ_act - 1)/batch_size + 1 202! 203 batch_size2 = batch_size !store for offsets 204! 205 if(print_mlcc3 .ge. 4) then 206 write(lupri,*) 207 write(lupri,*) 'work_remains:', work_remains 208 write(lupri,*) 'batch_size: ', batch_size 209 write(lupri,*) 'memory_use: ', memory_use 210 write(lupri,*) 'n_batches: ', n_batches 211 write(lupri,*) 212 end if 213! 214! Allocate arrays dependent on batchsize 215! 216! Fully transformed virtual integrals 217! 218 call work_allocator(b_D_c_k,n_vir_aag*batch_size) 219 call work_allocator(D_b_k_c,n_vir_aag*batch_size) 220! 221! Fully transformed occupied integrals 222! 223 call work_allocator(L_j_c_k,n_occ_int*batch_size) 224 call work_allocator(j_L_k_c,n_occ_int*batch_size) 225! 226! 227 if(.not. resp) then 228! 229! Fully transformed active integrals 230! 231 call work_allocator(j_b_k_c,n_act_int*batch_size) 232! 233 end if 234! 235! Partially transformed integrals 236! 237 call work_allocator(alpha_beta_k_pack_hole,n_basis_2_pack*batch_size) 238 call work_allocator(alpha_beta_k_pack_part,n_basis_2_pack*batch_size) 239! 240! 241! ----------------------------------------------------------------------------- 242! Open files for writing integrals outside batch loop, sequential direct access 243! We assume we will not reach the Fortran limit of 2GB record length 244! ----------------------------------------------------------------------------- 245! 246 call cpu_time(time_1) 247! 248 n_occ_print = n_occ_gen*n_vir_act 249! 250! (bD|ck) 251! 252 bDck_print_unit = new_unit() 253 open(unit=bDck_print_unit, file=bDck_name, action='write', status='replace', & 254 & access='direct', form='unformatted', recl=dp*n_vir_aag, iostat=ioerror) 255! 256 if(ioerror .ne. 0) then 257 write(lupri,*) "Error opening bDck file" 258 write(lupri,*) "Error code: ", ioerror 259 call quit('Error with opening file in intermediates') 260 end if 261! 262! (Db|kc) 263! 264 Dbkc_print_unit = new_unit() 265 open(unit=Dbkc_print_unit, file=Dbkc_name, action='write', status='replace', & 266 & access='direct', form='unformatted', recl=dp*n_vir_aag, iostat=ioerror) 267! 268 if(ioerror .ne. 0) then 269 write(lupri,*) "Error opening Dbkc file" 270 write(lupri,*) "Error code: ", ioerror 271 call quit('Error with opening file in intermediates') 272 end if 273! 274! (Lj|ck) 275! 276 Ljck_print_unit = new_unit() 277 open(unit=Ljck_print_unit, file=Ljck_name, action='write', status='replace', & 278 & access='direct', form='unformatted', recl=dp*n_occ_print, iostat=ioerror) 279! 280 if(ioerror .ne. 0) then 281 write(lupri,*) "Error opening Ljck file" 282 write(lupri,*) "Error code: ", ioerror 283 call quit('Error with opening file in intermediates') 284 end if 285! 286! (jL|kc) 287! 288 jLkc_print_unit = new_unit() 289 open(unit=jLkc_print_unit, file=jLkc_name, action='write', status='replace', & 290 & access='direct', form='unformatted', recl=dp*n_occ_print, iostat=ioerror) 291! 292 if(ioerror .ne. 0) then 293 write(lupri,*) "Error opening jLkc file" 294 write(lupri,*) "Error code: ", ioerror 295 call quit('Error with opening file in intermediates') 296 end if 297! 298! (jb|kc) 299! 300 if(.not. resp) then 301 jbkc_print_unit = new_unit() 302 open(unit=jbkc_print_unit, file=jbkc_name, action='write', status='replace', & 303 & access='direct', form='unformatted', recl=dp*n_vir_act**2, iostat=ioerror) 304! 305 if(ioerror .ne. 0) then 306 write(lupri,*) "Error opening jbkc file" 307 write(lupri,*) "Error code: ", ioerror 308 call quit('Error with opening file in intermediates') 309 end if 310! 311 end if 312! 313 call cpu_time(time_open) 314 time_open = time_open - time_1 315! 316! ------------------- 317! |Loop over batches| 318! ------------------- 319! 320 do batch = 1,n_batches 321! 322! 323! Set integral arrays to zero 324! --------------------------- 325! 326 call cpu_time(time_1) 327! 328 call dzero(b_D_c_k,n_vir_aag*batch_size) 329 call dzero(D_b_k_c,n_vir_aag*batch_size) 330 call dzero(L_j_c_k,n_occ_int*batch_size) 331 call dzero(j_L_k_c,n_occ_int*batch_size) 332 if(.not. resp) then 333 call dzero(j_b_k_c,n_act_int*batch_size) 334 end if 335! 336 call cpu_time(time_2) 337 time_zero = time_zero + time_2 - time_1 338! 339! 340! Find batch offsets 341! ------------------ 342! 343 batch_start = batch_size*(batch-1) + 1 344! 345 if (batch .eq. n_batches) then 346 batch_size = n_occ_act - batch_size*(n_batches - 1) 347 end if 348! 349 batch_end = batch_start + batch_size - 1 350! 351 if(print_mlcc3 .ge. 4) then 352 write(lupri,*) 353 write(lupri,*) 'n_batches: ', n_batches 354 write(lupri,*) 'batch: ', batch 355 write(lupri,*) 'batch_size: ', batch_size 356 write(lupri,*) 'batch_start: ', batch_start 357 write(lupri,*) 'batch_end: ', batch_end 358 write(lupri,*) 359 end if 360! 361! Loop over the delta index in (alpha beta | gamma delta) 362! 363 do delta = 1,n_basis 364! 365! Read the packed AO integrals into the array ao_integrals_pack 366! ------------------------------------------------------------- 367! 368 call cpu_time(time_1) 369! 370 call mlcc3_int_reader(delta) 371! 372 call cpu_time(time_2) 373 time_read = time_read + time_2 - time_1 374! 375! Fock matrix calculation. Calculate both t1-transformed and regular versions 376! --------------------------------------------------------------------------- 377! 378 call cpu_time(time_1) 379! 380 if (first) then 381! 382 call mlcc3_fock_calc(delta,resp) 383! 384 end if 385! 386 call cpu_time(time_2) 387 time_fock = time_fock + time_2 - time_1 388! 389! 390! Integral transform for triples code 391! ----------------------------------- 392 call cpu_time(time_1) 393! 394 if(.not. resp) then !Energy integrals 395 call mlcc3_integral_transform(batch,batch_size,batch_size2,delta) 396 else !response integrals 397 call mlcc3_integral_trans_resp(batch,batch_size,batch_size2,delta) 398 end if 399! 400 call cpu_time(time_2) 401 time_int = time_int + time_2 - time_1 402! 403 end do 404! 405! Do not calculate Fock matrix again 406 first = .false. 407! 408! 409! Print norms for debug 410! 411 if(print_mlcc3 .ge. 7) then 412! 413 bDcknorm = bDcknorm + ddot(n_vir_aag*batch_size,b_D_c_k,1,b_D_c_k,1) 414 Dbkcnorm = Dbkcnorm + ddot(n_vir_aag*batch_size,D_b_k_c,1,D_b_k_c,1) 415 Ljcknorm = Ljcknorm + ddot(n_occ_int*batch_size,L_j_c_k,1,L_j_c_k,1) 416 jLkcnorm = jLkcnorm + ddot(n_occ_int*batch_size,j_L_k_c,1,j_L_k_c,1) 417 if(.not. resp) then 418 jbkcnorm = jbkcnorm + ddot(n_act_int*batch_size,j_b_k_c,1,j_b_k_c,1) 419 end if 420! 421 write(lupri,*) 422 write(lupri,*) 'bDcknorm',bDcknorm 423 write(lupri,*) 'Dbkcnorm',Dbkcnorm 424 write(lupri,*) 'Ljcknorm',Ljcknorm 425 write(lupri,*) 'jLkcnorm',jLkcnorm 426 if(.not. resp) then 427 write(lupri,*) 'jbkcnorm',jbkcnorm 428 end if 429 write(lupri,*) 430! 431 end if 432! 433! 434! ----------------------------------------- 435! |Reorder the integrals and print to file| 436! ----------------------------------------- 437! 438 do k = batch_start,batch_end 439! 440 k_off = k - batch_start + 1 441! 442! Reorder bDck integrals to D b c, k order 443! 444 call cpu_time(time_1) 445! 446 call bDck_order(b_D_c_k,integral_sort,n_vir_act,n_vir_gen,k_off,batch_size) 447! 448 call cpu_time(time_2) 449 time_bDck = time_bDck + time_2 - time_1 450! 451! Write the integrals to disk. 452! 453 call cpu_time(time_1) 454! 455 write(bDck_print_unit,rec=k,iostat=ioerror) integral_sort(1:n_vir_aag) 456! 457 call cpu_time(time_2) 458 time_write = time_write + time_2 - time_1 459! 460 if(ioerror .ne. 0) then 461 write(lupri,*) "Error writing to bDck file" 462 write(lupri,*) "Error code: ", ioerror 463 call quit('Error writing integrals to file in intermediates') 464 end if 465! 466! Reorder Dbkc integrals to b c D, k order 467! 468 call cpu_time(time_1) 469! 470 call Dbkc_order(D_b_k_c,integral_sort,n_vir_act,n_vir_gen,k_off,batch_size) 471! 472 call cpu_time(time_2) 473 time_Dbkc = time_Dbkc + time_2 - time_1 474! 475! Write the integrals to disk. 476! 477 call cpu_time(time_1) 478! 479 write(Dbkc_print_unit,rec=k,iostat=ioerror) integral_sort(1:n_vir_aag) 480! 481 call cpu_time(time_2) 482 time_write = time_write + time_2 - time_1 483! 484 if(ioerror .ne. 0) then 485 write(lupri,*) "Error writing to Dbkc file" 486 write(lupri,*) "Error code: ", ioerror 487 call quit('Error writing integrals to file in intermediates') 488 end if 489! 490 do j = 1,n_occ_act 491! 492! 493! Find record number 494! 495 rec_number = n_occ_act*(k-1) + j 496! 497! Reorder Ljck integrals to L c, j k order 498! 499 call cpu_time(time_1) 500! 501 call Ljck_order(L_j_c_k,integral_sort,n_vir_act,n_vir_gen,& 502 & n_occ_act,n_occ_gen,j,k_off,batch_size) 503! 504 call cpu_time(time_2) 505 time_Ljck = time_Ljck + time_2 - time_1 506! 507! Write the integrals to disk. 508! 509 call cpu_time(time_1) 510! 511 write(Ljck_print_unit,rec=rec_number,iostat=ioerror) integral_sort(1:n_occ_print) 512! 513 call cpu_time(time_2) 514 time_write = time_write + time_2 - time_1 515! 516 if(ioerror .ne. 0) then 517 write(lupri,*) "Error writing to Ljck file" 518 write(lupri,*) "Error code: ", ioerror 519 call quit('Error writing integrals to file in intermediates') 520 end if 521! 522! Reorder jLkc integrals to c L, j k order 523! 524 call cpu_time(time_1) 525! 526 call jLkc_order(j_L_k_c,integral_sort,n_vir_act,n_vir_gen,& 527 & n_occ_act,n_occ_gen,j,k_off,batch_size) 528! 529 call cpu_time(time_2) 530 time_jLkc = time_jLkc + time_2 - time_1 531! 532! Write the integrals to disk. 533! 534 call cpu_time(time_1) 535! 536 write(jLkc_print_unit,rec=rec_number,iostat=ioerror) integral_sort(1:n_occ_print) 537! 538 call cpu_time(time_2) 539 time_write = time_write + time_2 - time_1 540! 541 if(ioerror .ne. 0) then 542 write(lupri,*) "Error writing to jLkc file" 543 write(lupri,*) "Error code: ", ioerror 544 call quit('Error writing integrals to file in intermediates') 545 end if 546! 547! 548 if(.not. resp) then 549! 550! Reorder jbkc integrals to b c, j k order 551! 552 call cpu_time(time_1) 553! 554 call jbkc_order(j_b_k_c,integral_sort,n_vir_act,n_vir_gen,& 555 & n_occ_act,n_occ_gen,j,k_off,batch_size) 556! 557 call cpu_time(time_2) 558 time_jbkc = time_jbkc + time_2 - time_1 559! 560! Write the integrals to disk. 561! 562 call cpu_time(time_1) 563! 564 write(jbkc_print_unit,rec=rec_number,iostat=ioerror) integral_sort(1:n_vir_act**2) 565! 566 call cpu_time(time_2) 567 time_write = time_write + time_2 - time_1 568! 569 if(ioerror .ne. 0) then 570 write(lupri,*) "Error writing to jbkc file" 571 write(lupri,*) "Error code: ", ioerror 572 call quit('Error writing integrals to file in intermediates') 573 end if 574! 575 end if 576! 577 end do 578! 579 end do 580! 581 call cpu_time(time_write) 582 time_write = time_write - time_1 583! 584! 585 end do 586! 587! ---------------------- 588! |Close integral files| 589! ---------------------- 590! 591! 592 call cpu_time(time_1) 593! 594 close(bDck_print_unit,status='keep',iostat=ioerror) 595! 596 if(ioerror .ne. 0) then 597 write(lupri,*) "Error closing bDck file" 598 write(lupri,*) "Error code: ", ioerror 599 call quit('Error closing integral file in intermediates') 600 end if 601! 602 close(Dbkc_print_unit,status='keep',iostat=ioerror) 603! 604 if(ioerror .ne. 0) then 605 write(lupri,*) "Error closing Dbkc file" 606 write(lupri,*) "Error code: ", ioerror 607 call quit('Error closing integral file in intermediates') 608 end if 609! 610 close(Ljck_print_unit,status='keep',iostat=ioerror) 611! 612 if(ioerror .ne. 0) then 613 write(lupri,*) "Error closing Ljck file" 614 write(lupri,*) "Error code: ", ioerror 615 call quit('Error closing integral file in intermediates') 616 end if 617! 618 close(jLkc_print_unit,status='keep',iostat=ioerror) 619! 620 if(ioerror .ne. 0) then 621 write(lupri,*) "Error closing jLkc file" 622 write(lupri,*) "Error code: ", ioerror 623 call quit('Error closing integral file in intermediates') 624 end if 625! 626! 627 if(.not. resp) then 628! 629 close(jbkc_print_unit,status='keep',iostat=ioerror) 630! 631 if(ioerror .ne. 0) then 632 write(lupri,*) "Error closing jbkc file" 633 write(lupri,*) "Error code: ", ioerror 634 call quit('Error closing integral file in intermediates') 635 end if 636! 637 end if 638! 639 call cpu_time(time_close) 640 time_close = time_close - time_1 641! 642! 643! ------------------------------------------------------------------------------- 644! |Transform Fock matrices to MO basis. Use the alpha_beta array as a help array| 645! |Note that the Fock matrix is calculated as the transpose. Easier to use | 646! ------------------------------------------------------------------------------- 647! 648! 649 call cpu_time(time_1) 650! 651 call mlcc3_fock_ao_mo(resp) 652! 653 call cpu_time(time_trans) 654 time_trans = time_trans - time_1 655! 656! 657! ------------------ 658! |Free memory used| 659! ------------------ 660! 661! Batch dependent allocations 662 call work_deallocator(alpha_beta_k_pack_part) 663 call work_deallocator(alpha_beta_k_pack_hole) 664! 665 if(.not. resp) then 666 call work_deallocator(j_b_k_c) 667 end if 668! 669 call work_deallocator(j_L_k_c) 670 call work_deallocator(L_j_c_k) 671! 672 call work_deallocator(D_b_k_c) 673 call work_deallocator(b_D_c_k) 674! 675! Sorted integrals 676 call work_deallocator(integral_sort) 677! 678! Various help arrays 679 call work_deallocator(j_b) 680 call work_deallocator(j_L) 681 call work_deallocator(L_j) 682! 683 call work_deallocator(D_b) 684 call work_deallocator(b_D) 685! 686 call work_deallocator(b_beta) 687 call work_deallocator(alpha_beta) 688 call work_deallocator(ao_integrals_pack) 689! 690! 691! 692! Total timing 693 call cpu_time(time_tot) 694! 695 time_tot = time_tot - time_start 696 time_reord = time_bDck + time_Dbkc + time_Ljck + time_jLkc + time_jbkc 697! 698! Print timings 699! ------------- 700! 701 if (print_mlcc3 .ge. 4) then 702! 703 write(lupri,*) 704 write(lupri,*) 705 write(lupri,*) 'Timings from mlcc3_int after loop' 706 write(lupri,9999) 'Total', time_tot 707 write(lupri,9999) 'open files', time_open 708 write(lupri,9999) 'zero arrays', time_zero 709 write(lupri,9999) 'read AO ints', time_read 710 write(lupri,9999) 'AO Fock mat', time_fock 711 write(lupri,9999) 'integral trans', time_int 712 write(lupri,9999) 'bDck order', time_bDck 713 write(lupri,9999) 'Dbkc order', time_Dbkc 714 write(lupri,9999) 'Ljck order', time_Ljck 715 write(lupri,9999) 'jLkc order', time_jLkc 716 write(lupri,9999) 'jbkc order', time_jbkc 717 write(lupri,9999) 'total order', time_reord 718 write(lupri,9999) 'write ints', time_write 719 write(lupri,9999) 'close files', time_close 720 write(lupri,9999) 'Fock transform', time_trans 721 write(lupri,*) 722 write(lupri,*) 723! 724 end if 725! 726! 727 end subroutine mlcc3_intermediates_calc 728! 729 subroutine mlcc3_int_reader(delta) 730! 731! integral reader 732! Authors Henrik Koch and Rolf H. Myhre 733! December 2014 734! 735! Purpose: read in (alpha beta | gamma delta) integrals for fixed delta 736! and alpha >= beta 737! 738 implicit none 739! 740 integer, intent(in) :: delta 741 integer :: dumdel = 1 742 logical :: direct_logic = .false. 743 real(dp), pointer, dimension(:) :: work_point_end => null() 744 integer :: work_remains 745! 746 work_point_end => work_end() 747 work_remains = work_free() 748! 749 call ccrdao(ao_integrals_pack,delta,dumdel,work_point_end,work_remains,dumdel,direct_logic) 750 751 end subroutine mlcc3_int_reader 752! 753! 754 subroutine mlcc3_fock_calc(delta,resp) 755! 756! Fock matrix calculator 757! Authors Henrik Koch and Rolf H. Myhre 758! December 2014 759! 760! Purpose: Calculate the various fock matrices using the alpha_beta array as a help array 761! 762! 763 implicit none 764! 765 integer, intent(in) :: delta 766 logical, intent(in) :: resp 767 integer :: gamm, gamm_off, gamm_off_pack, delta_off, gamma_delta 768 real(dp) :: ddot 769! 770! 771! 772! Calculate delta offset 773 delta_off = n_basis*(delta -1) + 1 774! 775 do gamm = 1,n_basis 776! 777! Calculate gamma offset 778 gamm_off = n_basis*(gamm-1) + 1 779 gamm_off_pack = n_basis_2_pack*(gamm-1) + 1 780 gamma_delta = n_basis*(gamm-1) + delta 781! 782! square up integral(alpha,beta)_gamm,delta 783! 784 call mlcc3_square_packed(ao_integrals_pack(gamm_off_pack:n_basis_2_pack),alpha_beta,n_basis) 785! 786! Add coloumb contribution. F_gamma,delta += 2(alpa beta | gamma delta)*D_alpha,beta 787! 788 mo_fock_mat(gamma_delta) = mo_fock_mat(gamma_delta) + & 789 & 2*ddot(n_basis_2,alpha_beta,1,ao_density,1) 790! 791 mo_fock_mat_t1(gamma_delta) = mo_fock_mat_t1(gamma_delta) + & 792 & 2*ddot(n_basis_2,alpha_beta,1,ao_density_t1,1) 793! 794 if(resp) then 795 mo_fock_mat_c1(gamma_delta) = mo_fock_mat_c1(gamma_delta) + & 796 & 2*ddot(n_basis_2,alpha_beta,1,ao_density_c1,1) 797 end if 798! 799! Subtract exchange contribution. F_gamma,delta -= (gamma beta | alpha delta)*D_alpha,beta 800! 801 call dgemv('N',n_basis,n_basis,-one,alpha_beta,n_basis, & 802 & ao_density(gamm_off),1,one,mo_fock_mat(delta_off),1) 803! 804 call dgemv('N',n_basis,n_basis,-one,alpha_beta,n_basis, & 805 & ao_density_t1(gamm_off),1,one,mo_fock_mat_t1(delta_off),1) 806! 807 if(resp) then 808 call dgemv('N',n_basis,n_basis,-one,alpha_beta,n_basis, & 809 & ao_density_c1(gamm_off),1,one,mo_fock_mat_c1(delta_off),1) 810 end if 811! 812 end do 813! 814! 815 end subroutine mlcc3_fock_calc 816! 817! 818 subroutine mlcc3_fock_ao_mo(resp) 819! 820! transform AO Fock matrices 821! Authors Henrik Koch and Rolf H. Myhre 822! December 2014 823! 824! Purpose: Transform Fock matrices from AO to MO basis 825! 826! 827 implicit none 828! 829 logical, intent(in) :: resp 830! 831! 832 if(resp) then 833! 834 call dgemm('T','N',n_orbitals,n_basis,n_basis,one,lambda_hole,n_basis, & 835 & mo_fock_mat_c1,n_basis,zero,alpha_beta,n_orbitals) 836 call dgemm('N','N',n_orbitals,n_orbitals,n_basis,one,alpha_beta,n_orbitals, & 837 & lambda_part,n_basis,zero,mo_fock_mat_c1,n_orbitals) 838! 839 call dgemm('T','N',n_orbitals,n_basis,n_basis,one,lambda_hole,n_basis, & 840 & mo_fock_mat_t1,n_basis,zero,alpha_beta,n_orbitals) 841 call dgemm('N','N',n_orbitals,n_orbitals,n_basis,one,alpha_beta,n_orbitals, & 842 & lambda_part_resp,n_basis,one,mo_fock_mat_c1,n_orbitals) 843! 844 call dgemm('T','N',n_orbitals,n_basis,n_basis,one,lambda_hole_resp,n_basis, & 845 & mo_fock_mat_t1,n_basis,zero,alpha_beta,n_orbitals) 846 call dgemm('N','N',n_orbitals,n_orbitals,n_basis,one,alpha_beta,n_orbitals, & 847 & lambda_part,n_basis,one,mo_fock_mat_c1,n_orbitals) 848! 849 end if 850! 851! 852 call dgemm('T','N',n_orbitals,n_basis,n_basis,one,lambda_hole,n_basis, & 853 & mo_fock_mat_t1,n_basis,zero,alpha_beta,n_orbitals) 854 call dgemm('N','N',n_orbitals,n_orbitals,n_basis,one,alpha_beta,n_orbitals, & 855 & lambda_part,n_basis,zero,mo_fock_mat_t1,n_orbitals) 856! 857! 858 call dgemm('T','N',n_orbitals,n_basis,n_basis,one,orb_coefficients,n_basis, & 859 & mo_fock_mat,n_basis,zero,alpha_beta,n_orbitals) 860 call dgemm('N','N',n_orbitals,n_orbitals,n_basis,one,alpha_beta,n_orbitals, & 861 & orb_coefficients,n_basis,zero,mo_fock_mat,n_orbitals) 862! 863 end subroutine mlcc3_fock_ao_mo 864! 865! 866! 867 subroutine mlcc3_integral_transform(batch,batch_size,batch_size2,delta) 868! 869! integral transformation routine 870! Authors Henrik Koch and Rolf H. Myhre 871! December 2014 872! 873! Purpose: Transform from the ao integrals to the integrals required 874! for the MLCC3 triples and write to disk 875! (alpha beta | gamma delta) -> (bD|ck) 876! (alpha beta | gamma delta) -> (Db|kc) 877! (alpha beta | gamma delta) -> (Lj|ck) 878! (alpha beta | gamma delta) -> (jL|kc) 879! (alpha beta | gamma delta) -> (jb|kc) 880! 881! 882 implicit none 883! 884 integer, intent(in) :: delta, batch, batch_size, batch_size2 885! 886 integer :: k, c 887 integer :: lambda_vir_off, lambda_occ_off, lambda_batch_off, lambda_gen_off 888 integer :: lambda_c_delta, int_off 889 integer :: alpha_beta_k_off,alpha_beta_k_end 890! 891! 892! 893! Calculate lambda offsets. Inactive occupied MOs come first in ordering 894! 895 lambda_occ_off = n_basis*n_occ_inact + 1 896 lambda_batch_off = n_basis*(n_occ_inact + batch_size2*(batch-1)) + 1 897 lambda_vir_off = n_basis*n_occ + 1 898 lambda_gen_off = n_basis*n_gen_inact + 1 899! 900! 901! Start with transforming gamma to the active occupied k, hole and particle 902! ------------------------------------------------------------------------- 903! 904 call dgemm('N','N',n_basis_2_pack,batch_size,n_basis,one,ao_integrals_pack,n_basis_2_pack, & 905 & lambda_hole(lambda_batch_off),n_basis,zero,alpha_beta_k_pack_hole,n_basis_2_pack) 906! 907 call dgemm('N','N',n_basis_2_pack,batch_size,n_basis,one,ao_integrals_pack,n_basis_2_pack, & 908 & lambda_part(lambda_batch_off),n_basis,zero,alpha_beta_k_pack_part,n_basis_2_pack) 909! 910! 911! Loop over the k indices in the batch 912! ------------------------------------ 913! 914 do k = 1,batch_size 915! 916 alpha_beta_k_off = n_basis_2_pack*(k-1) + 1 917 alpha_beta_k_end = n_basis_2_pack*k 918! 919! --------------------------------- 920! |First set of integrals, (bD|ck)| 921! --------------------------------- 922! 923! Square up hole transformed alpha beta 924! ------------------------------------- 925! 926 call mlcc3_square_packed(alpha_beta_k_pack_hole(alpha_beta_k_off:alpha_beta_k_end), & 927 & alpha_beta,n_basis) 928! 929! Transform alpha to b for bDck 930! ----------------------------- 931! 932 call dgemm('T','N',n_vir_act,n_basis,n_basis,one,lambda_part(lambda_vir_off),n_basis, & 933 & alpha_beta,n_basis,zero,b_beta,n_vir_act) 934! 935! Transform beta to D 936! ------------------- 937! 938 call dgemm('N','N',n_vir_act,n_vir_gen,n_basis,one,b_beta,n_vir_act, & 939 & lambda_hole(lambda_vir_off),n_basis,zero,b_D,n_vir_act) 940! 941! ---------------------------------- 942! |Second set of integrals, (Lj|ck)| 943! ---------------------------------- 944! 945! Transform beta to j for Ljck, reuse alpha_beta 946! ---------------------------------------------- 947! 948 call dgemm('N','N',n_basis,n_occ_act,n_basis,one,alpha_beta,n_basis, & 949 & lambda_hole(lambda_occ_off),n_basis,zero,b_beta,n_basis) 950! 951! Transform alpha to L 952! -------------------- 953! 954 call dgemm('T','N', n_occ_gen,n_occ_act,n_basis,one,lambda_part(lambda_gen_off),n_basis, & 955 & b_beta,n_basis,zero,L_j,n_occ_gen) 956! 957! 958! --------------------------------- 959! |Third set of integrals, (Db|kc)| 960! --------------------------------- 961! 962! Square up particle transformed alpha beta 963! ----------------------------------------- 964! 965 call mlcc3_square_packed(alpha_beta_k_pack_part(alpha_beta_k_off:alpha_beta_k_end), & 966 & alpha_beta,n_basis) 967! 968! Transform beta to b for Dbkc 969! ----------------------------- 970! 971 call dgemm('N','N',n_basis,n_vir_act,n_basis,one,alpha_beta,n_basis, & 972 & lambda_hole(lambda_vir_off),n_basis,zero,b_beta,n_basis) 973! 974! Transform alpha to D 975! -------------------- 976! 977 call dgemm('T','N', n_vir_gen,n_vir_act,n_basis,one,lambda_part(lambda_vir_off),n_basis, & 978 & b_beta,n_basis,zero,D_b,n_vir_gen) 979! 980! 981! ---------------------------------- 982! |Fourth set of integrals, (jL|kc)| 983! ---------------------------------- 984! 985! Transform alpha to j for jLkc, reuse alpha_beta 986! ----------------------------------------------- 987! 988 call dgemm('T','N',n_occ_act,n_basis,n_basis,one,lambda_part(lambda_occ_off),n_basis, & 989 & alpha_beta,n_basis,zero,b_beta,n_occ_act) 990! 991! Transform beta to L 992! ------------------- 993! 994 call dgemm('N','N',n_occ_act,n_occ_gen,n_basis,one,b_beta,n_occ_act, & 995 & lambda_hole(lambda_gen_off),n_basis,zero,j_L,n_occ_act) 996! 997! 998! --------------------------------- 999! |Fifth set of integrals, (jb|kc)| 1000! --------------------------------- 1001! 1002! 1003! Transform alpha to j for jbkc, reuse alpha_beta 1004! ----------------------------------------------- 1005! 1006 call dgemm('T','N',n_occ_act,n_basis,n_basis,one,lambda_part(lambda_occ_off),n_basis, & 1007 & alpha_beta,n_basis,zero,b_beta,n_occ_act) 1008! 1009! Transform beta to b for jbkc 1010! ----------------------------- 1011! 1012 call dgemm('N','N',n_occ_act,n_vir_act,n_basis,one,b_beta,n_occ_act, & 1013 & lambda_hole(lambda_vir_off),n_basis,zero,j_b,n_occ_act) 1014! 1015! 1016! Transform delta to c. Start by looping over c index 1017! --------------------------------------------------- 1018! 1019 do c = 1,n_vir_act 1020! 1021 lambda_c_delta = n_basis*(n_occ + c - 1) + delta 1022! 1023! 1024! Transform bDck 1025! -------------- 1026 int_off = n_vir_act*n_vir_gen*n_vir_act*(k-1) + n_vir_act*n_vir_gen*(c-1) + 1 1027! 1028 call daxpy(n_vir_act*n_vir_gen,lambda_part(lambda_c_delta),b_D,1,b_D_c_k(int_off),1) 1029! 1030! 1031! Transform Ljck 1032! -------------- 1033 int_off = n_occ_gen*n_occ_act*n_vir_act*(k-1) + n_occ_gen*n_occ_act*(c-1) + 1 1034 1035 call daxpy(n_occ_act*n_occ_gen,lambda_part(lambda_c_delta),L_j,1,L_j_c_k(int_off),1) 1036! 1037! 1038! Transform Dbkc 1039! -------------- 1040 int_off = n_vir_gen*n_vir_act*batch_size*(c-1) + n_vir_gen*n_vir_act*(k-1) + 1 1041! 1042 call daxpy(n_vir_act*n_vir_gen,lambda_hole(lambda_c_delta),D_b,1,D_b_k_c(int_off),1) 1043! 1044! 1045! Transform jLkc 1046! -------------- 1047 int_off = n_occ_act*n_occ_gen*batch_size*(c-1) + n_occ_act*n_occ_gen*(k-1) + 1 1048 1049 call daxpy(n_occ_act*n_occ_gen,lambda_hole(lambda_c_delta),j_L,1,j_L_k_c(int_off),1) 1050! 1051! 1052! Transform jbkc 1053! -------------- 1054 int_off = n_occ_act*n_vir_act*batch_size*(c-1) + n_occ_act*n_vir_act*(k-1) + 1 1055! 1056 call daxpy(n_occ_act*n_vir_act,lambda_hole(lambda_c_delta),j_b,1,j_b_k_c(int_off),1) 1057! 1058 end do 1059! 1060 end do 1061! 1062! 1063! 1064! 1065 end subroutine mlcc3_integral_transform 1066! 1067! 1068 subroutine mlcc3_integral_trans_resp(batch,batch_size,batch_size2,delta) 1069! 1070! integral transformation routine 1071! Authors Henrik Koch and Rolf H. Myhre 1072! December 2014 1073! 1074! Purpose: Transform from the ao integrals to the integrals required 1075! for the MLCC3 triples excitation energies and write to disk 1076! (alpha beta | gamma delta) -> (b'D|ck) + (bD|c'k) + (bD|ck') 1077! (alpha beta | gamma delta) -> (Db'|kc) 1078! (alpha beta | gamma delta) -> (Lj'|ck) + (Lj|c'k) + (Lj|ck') 1079! (alpha beta | gamma delta) -> (j'L|kc) 1080! 1081! 1082 implicit none 1083! 1084 integer, intent(in) :: delta, batch, batch_size, batch_size2 1085! 1086 integer :: k, c, i 1087 integer :: lambda_vir_off, lambda_occ_off, lambda_batch_off, lambda_gen_off 1088 integer :: lambda_c_delta, int_off 1089 integer :: alpha_beta_k_off, alpha_beta_k_end 1090! 1091 real(dp), dimension(:), pointer :: lambda_hole_1 => null() 1092! 1093 real(dp), dimension(:), pointer :: lambda_hole_2 => null() 1094 real(dp), dimension(:), pointer :: lambda_part_2 => null() 1095! 1096 real(dp), dimension(:), pointer :: lambda_part_3 => null() 1097! 1098! 1099! 1100! Calculate lambda offsets. Inactive occupied MOs come first in ordering 1101! 1102 lambda_occ_off = n_basis*n_occ_inact + 1 1103 lambda_batch_off = n_basis*(n_occ_inact + batch_size2*(batch-1)) + 1 1104 lambda_vir_off = n_basis*n_occ + 1 1105 lambda_gen_off = n_basis*n_gen_inact + 1 1106! 1107! 1108 do i=1,3 !Loop to get all contributions to |ck) integrals 1109! 1110! 1111! Set up lambda matrices for different contributions 1112! -------------------------------------------------- 1113! 1114 if(i .eq. 1) then !b and j C_1 transformed 1115! 1116 lambda_hole_1 => lambda_hole 1117 lambda_hole_2 => lambda_hole_resp 1118 lambda_part_2 => lambda_part_resp 1119 lambda_part_3 => lambda_part 1120! 1121 else if(i .eq. 2) then !k C_1 transformed 1122! 1123 lambda_hole_1 => lambda_hole 1124 lambda_hole_2 => lambda_hole 1125 lambda_part_2 => lambda_part 1126 lambda_part_3 => lambda_part_resp 1127! 1128 else !c C_1 transformed 1129! 1130 lambda_hole_1 => lambda_hole_resp 1131 lambda_hole_2 => lambda_hole 1132 lambda_part_2 => lambda_part 1133 lambda_part_3 => lambda_part 1134! 1135 end if 1136! 1137! Start with transforming gamma to the active occupied k, hole and particle 1138! ------------------------------------------------------------------------- 1139! 1140 if(i .ne. 2) then ! same in 1 and 2 1141! 1142 call dgemm('N','N',n_basis_2_pack,batch_size,n_basis,one,ao_integrals_pack,n_basis_2_pack, & 1143 & lambda_hole_1(lambda_batch_off),n_basis,zero,alpha_beta_k_pack_hole,n_basis_2_pack) 1144! 1145 end if 1146! 1147 if(i .eq. 1) then 1148! 1149 call dgemm('N','N',n_basis_2_pack,batch_size,n_basis,one,ao_integrals_pack,n_basis_2_pack, & 1150 & lambda_part(lambda_batch_off),n_basis,zero,alpha_beta_k_pack_part,n_basis_2_pack) 1151! 1152 end if 1153! 1154! 1155! Loop over the k indices in the batch 1156! ------------------------------------ 1157! 1158 do k = 1,batch_size 1159! 1160 alpha_beta_k_off = n_basis_2_pack*(k-1) + 1 1161 alpha_beta_k_end = n_basis_2_pack*k 1162! 1163! --------------------------------- 1164! |First set of integrals, (bD|ck)| 1165! --------------------------------- 1166! 1167! Square up hole transformed alpha beta 1168! ------------------------------------- 1169! 1170 call mlcc3_square_packed(alpha_beta_k_pack_hole(alpha_beta_k_off:alpha_beta_k_end), & 1171 & alpha_beta,n_basis) 1172! 1173! Transform alpha to b for bDck 1174! ----------------------------- 1175! 1176 call dgemm('T','N',n_vir_act,n_basis,n_basis,one,lambda_part_2(lambda_vir_off),n_basis, & 1177 & alpha_beta,n_basis,zero,b_beta,n_vir_act) 1178! 1179! Transform beta to D 1180! ------------------- 1181! 1182 call dgemm('N','N',n_vir_act,n_vir_gen,n_basis,one,b_beta,n_vir_act, & 1183 & lambda_hole(lambda_vir_off),n_basis,zero,b_D,n_vir_act) 1184! 1185! ---------------------------------- 1186! |Second set of integrals, (Lj|ck)| 1187! ---------------------------------- 1188! 1189! Transform beta to j for Ljck, reuse alpha_beta 1190! ---------------------------------------------- 1191! 1192 call dgemm('N','N',n_basis,n_occ_act,n_basis,one,alpha_beta,n_basis, & 1193 & lambda_hole_2(lambda_occ_off),n_basis,zero,b_beta,n_basis) 1194! 1195! Transform alpha to L 1196! -------------------- 1197! 1198 call dgemm('T','N', n_occ_gen,n_occ_act,n_basis,one,lambda_part(lambda_gen_off),n_basis, & 1199 & b_beta,n_basis,zero,L_j,n_occ_gen) 1200! 1201! 1202 if(i .eq. 1) then 1203! 1204! --------------------------------- 1205! |Third set of integrals, (Db|kc)| 1206! --------------------------------- 1207! 1208! Square up particle transformed alpha beta 1209! ----------------------------------------- 1210! 1211 call mlcc3_square_packed(alpha_beta_k_pack_part(alpha_beta_k_off:alpha_beta_k_end), & 1212 & alpha_beta,n_basis) 1213! 1214! Transform beta to b for Dbkc 1215! ----------------------------- 1216! 1217 call dgemm('N','N',n_basis,n_vir_act,n_basis,one,alpha_beta,n_basis, & 1218 & lambda_hole(lambda_vir_off),n_basis,zero,b_beta,n_basis) 1219! 1220! Transform alpha to D 1221! -------------------- 1222! 1223 call dgemm('T','N', n_vir_gen,n_vir_act,n_basis,one,lambda_part_resp(lambda_vir_off), & 1224 & n_basis,b_beta,n_basis,zero,D_b,n_vir_gen) 1225! 1226! 1227! ---------------------------------- 1228! |Fourth set of integrals, (jL|kc)| 1229! ---------------------------------- 1230! 1231! Transform alpha to j for jLkc, reuse alpha_beta 1232! ----------------------------------------------- 1233! 1234 call dgemm('T','N',n_occ_act,n_basis,n_basis,one,lambda_part(lambda_occ_off), & 1235 & n_basis,alpha_beta,n_basis,zero,b_beta,n_occ_act) 1236! 1237! Transform beta to L 1238! ------------------- 1239! 1240 call dgemm('N','N',n_occ_act,n_occ_gen,n_basis,one,b_beta,n_occ_act, & 1241 & lambda_hole_resp(lambda_gen_off),n_basis,zero,j_L,n_occ_act) 1242! 1243 end if 1244! 1245! 1246! 1247! Transform delta to c. Start by looping over c index 1248! --------------------------------------------------- 1249! 1250 do c = 1,n_vir_act 1251! 1252 lambda_c_delta = n_basis*(n_occ + c - 1) + delta 1253! 1254! 1255! Transform bDck 1256! -------------- 1257 int_off = n_vir_act*n_vir_gen*n_vir_act*(k-1) + n_vir_act*n_vir_gen*(c-1) + 1 1258! 1259 call daxpy(n_vir_act*n_vir_gen,lambda_part_3(lambda_c_delta),b_D,1,b_D_c_k(int_off),1) 1260! 1261! 1262! Transform Ljck 1263! -------------- 1264 int_off = n_occ_gen*n_occ_act*n_vir_act*(k-1) + n_occ_gen*n_occ_act*(c-1) + 1 1265 1266 call daxpy(n_occ_act*n_occ_gen,lambda_part_3(lambda_c_delta),L_j,1,L_j_c_k(int_off),1) 1267! 1268! 1269 if(i .eq. 1) then 1270! 1271! Transform Dbkc 1272! -------------- 1273 int_off = n_vir_gen*n_vir_act*batch_size*(c-1) + n_vir_gen*n_vir_act*(k-1) + 1 1274! 1275 call daxpy(n_vir_act*n_vir_gen,lambda_hole(lambda_c_delta),D_b,1,D_b_k_c(int_off),1) 1276! 1277! 1278! Transform jLkc 1279! -------------- 1280 int_off = n_occ_act*n_occ_gen*batch_size*(c-1) + n_occ_act*n_occ_gen*(k-1) + 1 1281 1282 call daxpy(n_occ_act*n_occ_gen,lambda_hole(lambda_c_delta),j_L,1,j_L_k_c(int_off),1) 1283! 1284 end if 1285! 1286 end do 1287! 1288 end do 1289! 1290 end do 1291! 1292 end subroutine mlcc3_integral_trans_resp 1293! 1294! 1295end module mlcc3_intermediates 1296