TABLE OF CONTENTS
- ABINIT/gwls_GenerateEpsilon
- gwls_GenerateEpsilon/driver_generate_dielectric_matrix
- gwls_GenerateEpsilon/Driver_GeneratePrintDielectricEigenvalues
- gwls_GenerateEpsilon/GeneratePrintDielectricEigenvalues
ABINIT/gwls_GenerateEpsilon [ Modules ]
NAME
gwls_GenerateEpsilon
FUNCTION
.
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (JLJ, BR, MC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
PARENTS
CHILDREN
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 28 29 module gwls_GenerateEpsilon 30 !---------------------------------------------------------------------------------------------------- 31 ! This module contains routines to compute and store the dielectric matrix, which plays a 32 ! central role in the self energy computations. In particular, global arrays are used to store 33 ! the static dielectric matrix. 34 !---------------------------------------------------------------------------------------------------- 35 ! local modules 36 use m_gwls_utility 37 use gwls_wf 38 use gwls_hamiltonian 39 use gwls_lineqsolver 40 use gwls_TimingLog 41 use gwls_polarisability 42 use gwls_model_polarisability 43 use gwls_GWlanczos 44 ! Abinit modules 45 use m_profiling_abi 46 use defs_basis 47 use defs_datatypes 48 use defs_abitypes 49 50 use m_io_tools, only : get_unit 51 52 implicit none 53 save 54 private
gwls_GenerateEpsilon/driver_generate_dielectric_matrix [ Functions ]
[ Top ] [ gwls_GenerateEpsilon ] [ Functions ]
NAME
driver_generate_dielectric_matrix
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_ComputeCorrelationEnergy
CHILDREN
cpu_time,diagonalize_lanczos_banded driver_invert_positive_definite_hermitian_matrix generateprintdielectriceigenvalues matrix_function_epsilon_model_operator set_dielectric_function_frequency,setup_pk_model,write_timing_log,zgemm zheevd
SOURCE
98 subroutine driver_generate_dielectric_matrix(epsilon_matrix_function,nseeds,kmax,& 99 epsilon_eigenvalues,Lbasis,debug) 100 !---------------------------------------------------------------------- 101 ! This routine computes the Lanczos approximate representation of the 102 ! implicit dielectic operator and then diagonalizes the banded 103 ! Lanczos matrix. 104 !---------------------------------------------------------------------- 105 106 !This section has been created automatically by the script Abilint (TD). 107 !Do not modify the following lines by hand. 108 #undef ABI_FUNC 109 #define ABI_FUNC 'driver_generate_dielectric_matrix' 110 !End of the abilint section 111 112 implicit none 113 interface 114 subroutine epsilon_matrix_function(v_out,v_in,l) 115 use defs_basis 116 117 integer, intent(in) :: l 118 complex(dp), intent(out) :: v_out(l) 119 complex(dp), intent(in) :: v_in(l) 120 121 end subroutine epsilon_matrix_function 122 end interface 123 124 integer, intent(in) :: nseeds, kmax 125 logical, intent(in) :: debug 126 127 real (dp), intent(out) :: epsilon_eigenvalues(nseeds*kmax) 128 complex(dpc), intent(out) :: Lbasis(npw_k,nseeds*kmax) ! array containing the Lanczos basis 129 130 131 ! local variables 132 133 complex(dpc), allocatable :: seeds(:,:) 134 135 complex(dpc),allocatable :: alpha(:,:,:) 136 complex(dpc),allocatable :: beta (:,:,:) 137 138 integer :: mpi_communicator 139 140 ! ************************************************************************* 141 142 143 ! The epsilon operator will act in LA mode. 144 mpi_communicator = mpi_enreg%comm_bandfft 145 146 147 !Create seeds 148 ABI_ALLOCATE(seeds,(npw_k,nseeds)) 149 call get_seeds(first_seed, nseeds, seeds) 150 151 ! compute the Lanczos basis 152 ABI_ALLOCATE(alpha,(nseeds,nseeds,kmax)) 153 ABI_ALLOCATE(beta ,(nseeds,nseeds,kmax)) 154 155 call block_lanczos_algorithm(mpi_communicator,epsilon_matrix_function,kmax,nseeds,npw_k, & 156 seeds,alpha,beta,Lbasis) 157 158 ! Diagonalize the epsilon matrix, which is banded 159 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug) 160 161 if (debug) then 162 call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues) 163 end if 164 165 ABI_DEALLOCATE(seeds) 166 ABI_DEALLOCATE(alpha) 167 ABI_DEALLOCATE(beta) 168 169 end subroutine driver_generate_dielectric_matrix
gwls_GenerateEpsilon/Driver_GeneratePrintDielectricEigenvalues [ Functions ]
[ Top ] [ gwls_GenerateEpsilon ] [ Functions ]
NAME
Driver_GeneratePrintDielectricEigenvalues
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_sternheimer
CHILDREN
cpu_time,diagonalize_lanczos_banded driver_invert_positive_definite_hermitian_matrix generateprintdielectriceigenvalues matrix_function_epsilon_model_operator set_dielectric_function_frequency,setup_pk_model,write_timing_log,zgemm zheevd
SOURCE
427 subroutine Driver_GeneratePrintDielectricEigenvalues(dtset) 428 !---------------------------------------------------------------------- 429 ! Compute the eigenvalues of the various dielectric operators 430 !---------------------------------------------------------------------- 431 432 !This section has been created automatically by the script Abilint (TD). 433 !Do not modify the following lines by hand. 434 #undef ABI_FUNC 435 #define ABI_FUNC 'Driver_GeneratePrintDielectricEigenvalues' 436 !End of the abilint section 437 438 type(dataset_type),intent(in) :: dtset 439 440 integer :: kmax_exact, kmax_model, kmax 441 real(dp) :: second_model_parameter 442 443 444 integer :: lm, k, lmax, l1, l2 445 integer :: io_unit 446 integer :: io_unit2 447 448 real(dp) :: time1, time2, time 449 ! local variables 450 character(128) :: output_filename 451 character(256) :: timing_string 452 453 454 complex(dpc), allocatable :: Lbasis_exact(:,:) 455 complex(dpc), allocatable :: Lbasis_model(:,:) 456 457 complex(dpc), allocatable :: sub_Lbasis_exact(:,:) 458 complex(dpc), allocatable :: sub_Lbasis_model(:,:) 459 460 complex(dpc), allocatable :: dummy(:,:) 461 complex(dpc), allocatable :: dummy2(:,:) 462 complex(dpc), allocatable :: dummy3(:,:) 463 464 complex(dpc), allocatable :: alpha_exact(:,:,:) 465 complex(dpc), allocatable :: beta_exact (:,:,:) 466 467 complex(dpc), allocatable :: alpha_model(:,:,:) 468 complex(dpc), allocatable :: beta_model (:,:,:) 469 470 real(dp), allocatable :: eig_exact(:) 471 real(dp), allocatable :: eig_model(:) 472 473 complex(dpc), allocatable :: model_epsilon_matrix(:,:) 474 complex(dpc), allocatable :: vector(:) 475 476 real(dpc) :: tr_eps_1, tr_eps_2, tr_eps_3 477 478 479 integer :: lwork, lrwork, liwork, info 480 complex(dpc), allocatable :: work(:) 481 real(dp) , allocatable :: rwork(:) 482 integer , allocatable :: iwork(:) 483 484 integer :: debug_unit 485 character(50) :: debug_filename 486 487 ! ************************************************************************* 488 489 490 491 492 kmax_exact = dtset%gwls_stern_kmax 493 kmax_model = dtset%gwls_kmax_complement 494 495 !second_model_parameter = dtset%gwls_second_model_parameter 496 second_model_parameter = zero 497 498 499 ! global stuff 500 nseeds = dtset%gwls_nseeds 501 first_seed = dtset%gwls_first_seed 502 e = dtset%gwls_band_index 503 504 505 506 ABI_ALLOCATE(Lbasis_exact,(npw_k,kmax_exact*nseeds)) 507 ABI_ALLOCATE(Lbasis_model,(npw_k,kmax_model*nseeds)) 508 509 ABI_ALLOCATE(alpha_exact, (nseeds,nseeds,kmax_exact)) 510 ABI_ALLOCATE(beta_exact , (nseeds,nseeds,kmax_exact)) 511 ABI_ALLOCATE(alpha_model, (nseeds,nseeds,kmax_model)) 512 ABI_ALLOCATE(beta_model , (nseeds,nseeds,kmax_model)) 513 514 515 ! set omega=0 for exact dielectric operator 516 call set_dielectric_function_frequency([zero,zero]) 517 518 519 520 call cpu_time(time1) 521 output_filename = 'EIGENVALUES_EXACT.dat' 522 call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_k, kmax_exact, & 523 output_filename, Lbasis_exact, alpha_exact, beta_exact) 524 525 526 527 call cpu_time(time2) 528 time = time2-time1 529 write(timing_string,'(A)') "Time to compute the EXACT Static Dielectric Matrix : " 530 call write_timing_log(timing_string,time) 531 532 533 534 call cpu_time(time1) 535 call setup_Pk_model(zero,second_model_parameter) 536 output_filename = 'EIGENVALUES_MODEL.dat' 537 call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_model_operator, kmax_model, output_filename,& 538 Lbasis_model, alpha_model, beta_model) 539 call cpu_time(time2) 540 time = time2-time1 541 write(timing_string,'(A)') "Time to compute the MODEL Static Dielectric Matrix : " 542 call write_timing_log(timing_string,time) 543 544 545 call cpu_time(time1) 546 if (kmax_exact <= kmax_model) then 547 kmax = kmax_exact 548 else 549 kmax = kmax_model 550 end if 551 lmax = nseeds*kmax 552 553 ! Build model operator matrix elements in the exact basis 554 ABI_ALLOCATE(model_epsilon_matrix, (lmax,lmax)) 555 ABI_ALLOCATE(vector, (npw_k)) 556 557 model_epsilon_matrix = cmplx_0 558 559 do l2 =1 , lmax 560 call matrix_function_epsilon_model_operator(vector ,Lbasis_exact(:,l2),npw_k) 561 562 563 do l1 =1, lmax 564 565 model_epsilon_matrix(l1, l2) = complex_vector_product(Lbasis_exact(:,l1),vector,npw_k) 566 567 end do 568 569 570 end do 571 572 ABI_DEALLOCATE(vector) 573 574 call cpu_time(time2) 575 time = time2-time1 576 write(timing_string,'(A)') "Compute MODEL matrix elements in EXACT basis : " 577 call write_timing_log(timing_string,time) 578 579 580 581 582 ! Compare the traces 583 584 585 io_unit = get_unit() 586 open(file='DIELECTRIC_TRACE.dat',status=files_status_new,unit=io_unit) 587 write(io_unit,10) '#----------------------------------------------------------------------------' 588 write(io_unit,10) '# ' 589 write(io_unit,10) '# Partial traces ' 590 write(io_unit,10) '# ========================================== ' 591 write(io_unit,10) '# ' 592 write(io_unit,10) '# Tabulate the trace of various operators, as function of the number ' 593 write(io_unit,10) '# of lanczos steps performed. ' 594 write(io_unit,10) '# ' 595 write(io_unit,10) '# ' 596 write(io_unit,10) '# NOTES: ' 597 write(io_unit,10) '# Tr[1-eps^{-1}] is evaluated in the Lanczos basis of eps ' 598 write(io_unit,10) '# Tr[1-eps_m^{-1}] is evaluated in the Lanczos basis of eps_m ' 599 write(io_unit,10) '# Tr[eps_m^{-1}-eps^{-1}] is evaluated in the Lanczos basis of eps ' 600 write(io_unit,10) '# ' 601 write(io_unit,10) '# ' 602 write(io_unit,10) '# k Tr[1-eps^{-1}] Tr[1-eps_m^{-1}] Tr[eps_m^{-1}-eps^{-1}]' 603 write(io_unit,10) '#----------------------------------------------------------------------------' 604 flush(io_unit) 605 606 607 io_unit2 = get_unit() 608 open(file='RPA_ENERGY.dat',status=files_status_new,unit=io_unit2) 609 write(io_unit2,10) '#----------------------------------------------------------------------------' 610 write(io_unit2,10) '# ' 611 write(io_unit2,10) '# RPA TOTAL ENERGY ' 612 write(io_unit2,10) '# ========================================== ' 613 write(io_unit2,10) '# ' 614 write(io_unit2,10) '# It can be shown that the correlation energy, within the RPA, is given ' 615 write(io_unit2,10) '# by: ' 616 write(io_unit2,10) '# E_c = int_0^{infty} dw/(2pi) Tr[ ln(eps{iw)}+1-eps(iw)] ' 617 write(io_unit2,10) '# ' 618 write(io_unit2,10) '# As a gauge of what can be expected as far as convergence is concerned, ' 619 write(io_unit2,10) '# the following will be printed. ' 620 write(io_unit2,10) '# ' 621 write(io_unit2,10) '# I_1 = Tr[ ln(eps) + 1 - eps ] ' 622 write(io_unit2,10) '# I_2 = Tr[ ln(eps_m) + 1 - eps_m ] ' 623 write(io_unit2,10) '# I_3 = Tr[ ln(eps^{-1}_m . eps ) + eps_m - eps ] ' 624 write(io_unit2,10) '# ' 625 write(io_unit2,10) '# ' 626 write(io_unit2,10) '# ' 627 write(io_unit2,10) '# k I_1 I_2 I_3 ' 628 write(io_unit2,10) '#----------------------------------------------------------------------------' 629 flush(io_unit2) 630 631 632 ! Iterate every 10 values of k max, or else the linear algebra gets too expensive... 633 do k = 4, kmax, 4 634 635 ABI_ALLOCATE(sub_Lbasis_exact,(npw_k,k*nseeds)) 636 ABI_ALLOCATE(sub_Lbasis_model,(npw_k,k*nseeds)) 637 638 639 ABI_ALLOCATE(eig_exact,(k*nseeds)) 640 ABI_ALLOCATE(eig_model,(k*nseeds)) 641 ABI_ALLOCATE(dummy,(k*nseeds,k*nseeds)) 642 ABI_ALLOCATE(dummy2,(k*nseeds,k*nseeds)) 643 644 sub_Lbasis_exact(:,:) = Lbasis_exact(:,1:k*nseeds) 645 sub_Lbasis_model(:,:) = Lbasis_model(:,1:k*nseeds) 646 647 ! Diagonalize the epsilon matrix, which is banded 648 call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_exact(:,:,1:k),beta_exact(:,:,1:k),sub_Lbasis_exact,eig_exact,.false.) 649 call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_model(:,:,1:k),beta_model(:,:,1:k),sub_Lbasis_model,eig_model,.false.) 650 651 tr_eps_1 = sum(one-one/eig_exact(:)) 652 tr_eps_2 = sum(one-one/eig_model(:)) 653 654 tr_eps_3 = -sum(one/eig_exact(:)) 655 656 dummy(:,:) = model_epsilon_matrix(1:k*nseeds, 1:k*nseeds) 657 658 call driver_invert_positive_definite_hermitian_matrix(dummy,k*nseeds) 659 660 do lm = 1, k*nseeds 661 662 tr_eps_3 = tr_eps_3 + dble(dummy(lm,lm)) 663 end do 664 665 666 write(io_unit,20) k, tr_eps_1, tr_eps_2, tr_eps_3 667 flush(io_unit) 668 669 670 671 tr_eps_1 = sum(log(eig_exact(:))+one-eig_exact(:)) 672 tr_eps_2 = sum(log(eig_model(:))+one-eig_model(:)) 673 674 dummy2(:,:) = zero 675 tr_eps_3 = zero 676 do lm = 1, k*nseeds 677 dummy2(lm,lm) = eig_exact(lm) 678 tr_eps_3 = tr_eps_3 + dble(model_epsilon_matrix(lm,lm)) - dble(eig_exact(lm)) 679 end do 680 681 ABI_ALLOCATE(dummy3,(k*nseeds,k*nseeds)) 682 call ZGEMM( 'N', & ! Hermitian conjugate the first array 683 'N', & ! Leave second array as is 684 k*nseeds, & ! the number of rows of the matrix op( A ) 685 k*nseeds, & ! the number of columns of the matrix op( B ) 686 k*nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B ) 687 cmplx_1, & ! alpha constant 688 dummy2, & ! matrix A 689 k*nseeds, & ! LDA 690 dummy, & ! matrix B 691 k*nseeds, & ! LDB 692 cmplx_0, & ! beta constant 693 dummy3, & ! matrix C 694 k*nseeds) ! LDC 695 696 dummy2(:,:) = dummy3(:,:) 697 ABI_DEALLOCATE(dummy3) 698 699 ! find eigenvalues 700 !call heevd(dummy2, eig_exact) 701 702 lwork = k*nseeds+1 703 lrwork = k*nseeds 704 liwork = 1 705 706 ABI_ALLOCATE(work,(lwork)) 707 ABI_ALLOCATE(rwork,(lrwork)) 708 ABI_ALLOCATE(iwork,(liwork)) 709 710 call zheevd('N', 'U',k*nseeds, dummy2, k*nseeds, eig_exact, work, lwork, rwork, lrwork, iwork, liwork, info) 711 if ( info /= 0) then 712 debug_unit = get_unit() 713 write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log' 714 715 open(debug_unit,file=trim(debug_filename),status='unknown') 716 717 write(debug_unit,'(A)') '*************************************************************************************' 718 write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZHEEVD(1), gwls_GenerateEpsilon' 719 write(debug_unit,'(A)') '*************************************************************************************' 720 721 close(debug_unit) 722 723 end if 724 725 726 727 728 729 tr_eps_3 = tr_eps_3 + sum(log(eig_exact)) 730 731 write(io_unit2,20) k, tr_eps_1, tr_eps_2, tr_eps_3 732 flush(io_unit2) 733 734 735 736 ABI_DEALLOCATE(work) 737 ABI_DEALLOCATE(rwork) 738 ABI_DEALLOCATE(iwork) 739 740 741 742 ABI_DEALLOCATE(sub_Lbasis_exact) 743 ABI_DEALLOCATE(sub_Lbasis_model) 744 745 ABI_DEALLOCATE(eig_exact) 746 ABI_DEALLOCATE(eig_model) 747 ABI_DEALLOCATE(dummy) 748 ABI_DEALLOCATE(dummy2) 749 end do 750 751 close(io_unit) 752 close(io_unit2) 753 754 call cpu_time(time2) 755 time = time2-time1 756 write(timing_string,'(A)') "Time to compute the TRACES of the Dielectric Matrices: " 757 call write_timing_log(timing_string,time) 758 759 ABI_DEALLOCATE(Lbasis_exact) 760 ABI_DEALLOCATE(Lbasis_model) 761 ABI_DEALLOCATE(alpha_exact) 762 ABI_DEALLOCATE(beta_exact ) 763 ABI_DEALLOCATE(alpha_model) 764 ABI_DEALLOCATE(beta_model ) 765 ABI_DEALLOCATE(model_epsilon_matrix) 766 767 768 769 10 format(A) 770 20 format(I5,3ES24.16) 771 772 end subroutine Driver_GeneratePrintDielectricEigenvalues
gwls_GenerateEpsilon/GeneratePrintDielectricEigenvalues [ Functions ]
[ Top ] [ gwls_GenerateEpsilon ] [ Functions ]
NAME
GeneratePrintDielectricEigenvalues
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_GenerateEpsilon
CHILDREN
cpu_time,diagonalize_lanczos_banded driver_invert_positive_definite_hermitian_matrix generateprintdielectriceigenvalues matrix_function_epsilon_model_operator set_dielectric_function_frequency,setup_pk_model,write_timing_log,zgemm zheevd
SOURCE
196 subroutine GeneratePrintDielectricEigenvalues(epsilon_matrix_function,kmax,output_filename,Lbasis,alpha,beta) 197 !---------------------------------------------------------------------- 198 ! This routine computes the Lanczos approximate representation of the 199 ! implicit dielectic operator and then diagonalizes the banded 200 ! Lanczos matrix. 201 !---------------------------------------------------------------------- 202 203 !This section has been created automatically by the script Abilint (TD). 204 !Do not modify the following lines by hand. 205 #undef ABI_FUNC 206 #define ABI_FUNC 'GeneratePrintDielectricEigenvalues' 207 !End of the abilint section 208 209 implicit none 210 interface 211 subroutine epsilon_matrix_function(v_out,v_in,l) 212 use defs_basis 213 214 integer, intent(in) :: l 215 complex(dp), intent(out) :: v_out(l) 216 complex(dp), intent(in) :: v_in(l) 217 218 end subroutine epsilon_matrix_function 219 end interface 220 221 integer, intent(in) :: kmax 222 223 character(*), intent(in) :: output_filename 224 225 226 complex(dpc), intent(out) :: Lbasis(npw_k,nseeds*kmax) 227 complex(dpc), intent(out) :: alpha(nseeds,nseeds,kmax) 228 complex(dpc), intent(out) :: beta (nseeds,nseeds,kmax) 229 230 231 ! local variables 232 233 234 complex(dpc),allocatable :: seeds(:,:) 235 236 complex(dpc),allocatable :: Lbasis_diag(:,:) 237 238 239 real(dp), allocatable :: psik(:,:) 240 real(dp), allocatable :: psir(:,:,:,:) 241 242 real(dp), allocatable :: epsilon_eigenvalues(:) 243 244 245 integer :: mpi_communicator 246 integer :: io_unit 247 integer :: lmax 248 integer :: l 249 integer :: ir1, ir2, ir3 250 integer :: n1, n2, n3 251 252 real(dp) :: R, G 253 real(dp) :: sigma_R, sigma_G 254 real(dp) :: x, y, z 255 256 real(dp),allocatable :: G_array(:) 257 real(dp),allocatable :: R_array(:,:,:) 258 259 logical :: debug 260 261 ! ************************************************************************* 262 263 264 debug = .false. 265 lmax = kmax*nseeds 266 mpi_communicator = mpi_enreg%comm_bandfft 267 !Create seeds 268 ABI_ALLOCATE(seeds,(npw_k,nseeds)) 269 call get_seeds(first_seed, nseeds, seeds) 270 271 ! compute the Lanczos basis 272 ABI_ALLOCATE(Lbasis_diag,(npw_k,lmax)) 273 ABI_ALLOCATE(epsilon_eigenvalues,(lmax)) 274 275 ABI_ALLOCATE(psik,(2,npw_k)) 276 ABI_ALLOCATE(psir,(2,n4,n5,n6)) 277 ABI_ALLOCATE(G_array,(npw_k)) 278 ABI_ALLOCATE(R_array,(n4,n5,n6)) 279 280 psir = zero 281 R_array = zero 282 283 n1 = n4-1 284 n2 = n5-1 285 n3 = n6 286 287 ! Generate the Lanczos basis and banded eigenvalue representation 288 call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k, seeds,alpha,beta,Lbasis) 289 290 Lbasis_diag = Lbasis 291 292 ! Diagonalize the epsilon matrix, which is banded 293 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis_diag,epsilon_eigenvalues,debug) 294 295 call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,lmax,npw_k,Lbasis_diag,epsilon_eigenvalues) 296 297 io_unit = get_unit() 298 open(file=output_filename,status=files_status_new,unit=io_unit) 299 write(io_unit,10) '#----------------------------------------------------------------------------' 300 write(io_unit,10) '# ' 301 write(io_unit,10) '# Partial eigenvalues ' 302 write(io_unit,10) '# ========================================== ' 303 write(io_unit,10) '# ' 304 write(io_unit,10) '# Tabulate the eigenvalues of the dielectic matrix, as well as some ' 305 write(io_unit,10) '# information regarding the eigenstates. ' 306 write(io_unit,10) '# ' 307 write(io_unit,10) '# ' 308 write(io_unit,10) '# definitions: ' 309 write(io_unit,10) '# l index of the eigenvalue ' 310 write(io_unit,10) '# eig eigenvalue ' 311 write(io_unit,10) '# ' 312 write(io_unit,10) '# ' 313 write(io_unit,10) '# (the following vectors are expressed in crystal units) ' 314 write(io_unit,10) '# ' 315 write(io_unit,10) '# R = < l | |r| | l > ' 316 write(io_unit,10) '# sigma_R = sqrt{< l | (|r|-R)^2 | l > } ' 317 write(io_unit,10) '# ' 318 write(io_unit,10) '# G = < l | |G| | l > ' 319 write(io_unit,10) '# sigma_G = sqrt{< l | (|G|-G)^2 | l > } ' 320 write(io_unit,10) '# ' 321 write(io_unit,10) '# ' 322 write(io_unit,10) '# l eig r sigma_r G sigma_G ' 323 write(io_unit,10) '#----------------------------------------------------------------------------' 324 flush(io_unit) 325 326 327 G_array(:) = kg_k(1,:)**2+ kg_k(2,:)**2+ kg_k(3,:)**2 328 329 G_array(:) = sqrt(G_array(:)) 330 331 R_array = zero 332 333 do ir3=1,n3 334 335 if (ir3 <= n3/2 ) then 336 z = (one*ir3)/(one*n3) 337 else 338 z = (one*ir3)/(one*n3)-one 339 end if 340 341 do ir2=1,n2 342 343 if (ir2 <= n2/2 ) then 344 y = (one*ir2)/(one*n2) 345 else 346 y = (one*ir2)/(one*n2)-one 347 end if 348 349 do ir1=1,n1 350 351 if (ir1 <= n1/2 ) then 352 x = (one*ir1)/(one*n1) 353 else 354 x = (one*ir1)/(one*n1)-one 355 end if 356 357 358 R_array(ir1,ir2,ir3) = sqrt(x**2+y**2+z**2) 359 end do 360 end do 361 end do 362 363 364 365 do l=1, lmax 366 367 psik(1,:) = dble (Lbasis_diag(:,l)) 368 psik(2,:) = dimag(Lbasis_diag(:,l)) 369 370 call g_to_r(psir ,psik) 371 372 373 G = sum(G_array(:)*(psik(1,:)**2+psik(2,:)**2)) 374 R = sum(R_array(:,:,:)*(psir(1,:,:,:)**2+psir(2,:,:,:)**2) )*ucvol/nfft 375 376 sigma_G = sqrt(sum((G_array(:) -G)**2*(psik(1,:)**2 +psik(2,:)**2))) 377 sigma_R = sqrt(sum((R_array(:,:,:)-R)**2*(psir(1,:,:,:)**2+psir(2,:,:,:)**2))*ucvol/nfft) 378 379 380 381 write(io_unit,20) l, epsilon_eigenvalues(l), R,sigma_R, G,sigma_G 382 383 end do 384 385 close(io_unit) 386 387 ABI_DEALLOCATE(seeds) 388 ABI_DEALLOCATE(Lbasis_diag) 389 ABI_DEALLOCATE(psik) 390 ABI_DEALLOCATE(psir) 391 ABI_DEALLOCATE(G_array) 392 ABI_DEALLOCATE(R_array) 393 394 ABI_DEALLOCATE(epsilon_eigenvalues) 395 396 397 10 format(A) 398 20 format(I5,ES24.16,4F12.8) 399 400 end subroutine GeneratePrintDielectricEigenvalues