TABLE OF CONTENTS
- ABINIT/m_gwls_GenerateEpsilon
- m_gwls_GenerateEpsilon/driver_generate_dielectric_matrix
- m_gwls_GenerateEpsilon/Driver_GeneratePrintDielectricEigenvalues
- m_gwls_GenerateEpsilon/GeneratePrintDielectricEigenvalues
ABINIT/m_gwls_GenerateEpsilon [ Modules ]
NAME
m_gwls_GenerateEpsilon
FUNCTION
.
COPYRIGHT
Copyright (C) 2009-2022 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 23 module m_gwls_GenerateEpsilon 24 !---------------------------------------------------------------------------------------------------- 25 ! This module contains routines to compute and store the dielectric matrix, which plays a 26 ! central role in the self energy computations. In particular, global arrays are used to store 27 ! the static dielectric matrix. 28 !---------------------------------------------------------------------------------------------------- 29 ! local modules 30 use m_gwls_utility 31 use m_gwls_wf 32 use m_gwls_hamiltonian 33 use m_gwls_lineqsolver 34 use m_gwls_TimingLog 35 use m_gwls_polarisability 36 use m_gwls_model_polarisability 37 use m_gwls_GWlanczos 38 ! Abinit modules 39 use m_abicore 40 use defs_basis 41 use m_dtset 42 43 use m_io_tools, only : get_unit 44 45 implicit none 46 save 47 private
m_gwls_GenerateEpsilon/driver_generate_dielectric_matrix [ Functions ]
[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]
NAME
driver_generate_dielectric_matrix
FUNCTION
.
INPUTS
OUTPUT
SOURCE
80 subroutine driver_generate_dielectric_matrix(epsilon_matrix_function,nseeds,kmax,& 81 epsilon_eigenvalues,Lbasis,debug) 82 !---------------------------------------------------------------------- 83 ! This routine computes the Lanczos approximate representation of the 84 ! implicit dielectic operator and then diagonalizes the banded 85 ! Lanczos matrix. 86 !---------------------------------------------------------------------- 87 interface 88 subroutine epsilon_matrix_function(v_out,v_in,l) 89 use defs_basis 90 91 integer, intent(in) :: l 92 complex(dp), intent(out) :: v_out(l) 93 complex(dp), intent(in) :: v_in(l) 94 95 end subroutine epsilon_matrix_function 96 end interface 97 98 integer, intent(in) :: nseeds, kmax 99 logical, intent(in) :: debug 100 101 real (dp), intent(out) :: epsilon_eigenvalues(nseeds*kmax) 102 complex(dpc), intent(out) :: Lbasis(npw_k,nseeds*kmax) ! array containing the Lanczos basis 103 104 105 ! local variables 106 107 complex(dpc), allocatable :: seeds(:,:) 108 109 complex(dpc),allocatable :: alpha(:,:,:) 110 complex(dpc),allocatable :: beta (:,:,:) 111 112 integer :: mpi_communicator 113 114 ! ************************************************************************* 115 116 117 ! The epsilon operator will act in LA mode. 118 mpi_communicator = mpi_enreg%comm_bandfft 119 120 121 !Create seeds 122 ABI_MALLOC(seeds,(npw_k,nseeds)) 123 call get_seeds(first_seed, nseeds, seeds) 124 125 ! compute the Lanczos basis 126 ABI_MALLOC(alpha,(nseeds,nseeds,kmax)) 127 ABI_MALLOC(beta ,(nseeds,nseeds,kmax)) 128 129 call block_lanczos_algorithm(mpi_communicator,epsilon_matrix_function,kmax,nseeds,npw_k, & 130 seeds,alpha,beta,Lbasis) 131 132 ! Diagonalize the epsilon matrix, which is banded 133 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug) 134 135 if (debug) then 136 call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues) 137 end if 138 139 ABI_FREE(seeds) 140 ABI_FREE(alpha) 141 ABI_FREE(beta) 142 143 end subroutine driver_generate_dielectric_matrix
m_gwls_GenerateEpsilon/Driver_GeneratePrintDielectricEigenvalues [ Functions ]
[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]
NAME
Driver_GeneratePrintDielectricEigenvalues
FUNCTION
.
INPUTS
OUTPUT
SOURCE
371 subroutine Driver_GeneratePrintDielectricEigenvalues(dtset) 372 !---------------------------------------------------------------------- 373 ! Compute the eigenvalues of the various dielectric operators 374 !---------------------------------------------------------------------- 375 type(dataset_type),intent(in) :: dtset 376 377 integer :: kmax_exact, kmax_model, kmax 378 real(dp) :: second_model_parameter 379 380 381 integer :: lm, k, lmax, l1, l2 382 integer :: io_unit 383 integer :: io_unit2 384 385 real(dp) :: time1, time2, time 386 ! local variables 387 character(128) :: output_filename 388 character(256) :: timing_string 389 390 391 complex(dpc), allocatable :: Lbasis_exact(:,:) 392 complex(dpc), allocatable :: Lbasis_model(:,:) 393 394 complex(dpc), allocatable :: sub_Lbasis_exact(:,:) 395 complex(dpc), allocatable :: sub_Lbasis_model(:,:) 396 397 complex(dpc), allocatable :: dummy(:,:) 398 complex(dpc), allocatable :: dummy2(:,:) 399 complex(dpc), allocatable :: dummy3(:,:) 400 401 complex(dpc), allocatable :: alpha_exact(:,:,:) 402 complex(dpc), allocatable :: beta_exact (:,:,:) 403 404 complex(dpc), allocatable :: alpha_model(:,:,:) 405 complex(dpc), allocatable :: beta_model (:,:,:) 406 407 real(dp), allocatable :: eig_exact(:) 408 real(dp), allocatable :: eig_model(:) 409 410 complex(dpc), allocatable :: model_epsilon_matrix(:,:) 411 complex(dpc), allocatable :: vector(:) 412 413 real(dpc) :: tr_eps_1, tr_eps_2, tr_eps_3 414 415 416 integer :: lwork, lrwork, liwork, info 417 complex(dpc), allocatable :: work(:) 418 real(dp) , allocatable :: rwork(:) 419 integer , allocatable :: iwork(:) 420 421 integer :: debug_unit 422 character(50) :: debug_filename 423 424 ! ************************************************************************* 425 426 427 428 429 kmax_exact = dtset%gwls_stern_kmax 430 kmax_model = dtset%gwls_kmax_complement 431 432 !second_model_parameter = dtset%gwls_second_model_parameter 433 second_model_parameter = zero 434 435 436 ! global stuff 437 nseeds = dtset%gwls_nseeds 438 first_seed = dtset%gwls_first_seed 439 e = dtset%gwls_band_index 440 441 442 443 ABI_MALLOC(Lbasis_exact,(npw_k,kmax_exact*nseeds)) 444 ABI_MALLOC(Lbasis_model,(npw_k,kmax_model*nseeds)) 445 446 ABI_MALLOC(alpha_exact, (nseeds,nseeds,kmax_exact)) 447 ABI_MALLOC(beta_exact , (nseeds,nseeds,kmax_exact)) 448 ABI_MALLOC(alpha_model, (nseeds,nseeds,kmax_model)) 449 ABI_MALLOC(beta_model , (nseeds,nseeds,kmax_model)) 450 451 452 ! set omega=0 for exact dielectric operator 453 call set_dielectric_function_frequency([zero,zero]) 454 455 456 457 call cpu_time(time1) 458 output_filename = 'EIGENVALUES_EXACT.dat' 459 call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_k, nseeds, kmax_exact, & 460 output_filename, Lbasis_exact, alpha_exact, beta_exact) 461 462 463 464 call cpu_time(time2) 465 time = time2-time1 466 write(timing_string,'(A)') "Time to compute the EXACT Static Dielectric Matrix : " 467 call write_timing_log(timing_string,time) 468 469 470 471 call cpu_time(time1) 472 call setup_Pk_model(zero,second_model_parameter) 473 output_filename = 'EIGENVALUES_MODEL.dat' 474 call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_model_operator, nseeds, kmax_model, & 475 &output_filename, Lbasis_model, alpha_model, beta_model) 476 call cpu_time(time2) 477 time = time2-time1 478 write(timing_string,'(A)') "Time to compute the MODEL Static Dielectric Matrix : " 479 call write_timing_log(timing_string,time) 480 481 482 call cpu_time(time1) 483 if (kmax_exact <= kmax_model) then 484 kmax = kmax_exact 485 else 486 kmax = kmax_model 487 end if 488 lmax = nseeds*kmax 489 490 ! Build model operator matrix elements in the exact basis 491 ABI_MALLOC(model_epsilon_matrix, (lmax,lmax)) 492 ABI_MALLOC(vector, (npw_k)) 493 494 model_epsilon_matrix = cmplx_0 495 496 do l2 =1 , lmax 497 call matrix_function_epsilon_model_operator(vector ,Lbasis_exact(:,l2),npw_k) 498 499 500 do l1 =1, lmax 501 502 model_epsilon_matrix(l1, l2) = complex_vector_product(Lbasis_exact(:,l1),vector,npw_k) 503 504 end do 505 506 507 end do 508 509 ABI_FREE(vector) 510 511 call cpu_time(time2) 512 time = time2-time1 513 write(timing_string,'(A)') "Compute MODEL matrix elements in EXACT basis : " 514 call write_timing_log(timing_string,time) 515 516 517 518 519 ! Compare the traces 520 521 522 io_unit = get_unit() 523 open(file='DIELECTRIC_TRACE.dat',status=files_status_new,unit=io_unit) 524 write(io_unit,10) '#----------------------------------------------------------------------------' 525 write(io_unit,10) '# ' 526 write(io_unit,10) '# Partial traces ' 527 write(io_unit,10) '# ========================================== ' 528 write(io_unit,10) '# ' 529 write(io_unit,10) '# Tabulate the trace of various operators, as function of the number ' 530 write(io_unit,10) '# of lanczos steps performed. ' 531 write(io_unit,10) '# ' 532 write(io_unit,10) '# ' 533 write(io_unit,10) '# NOTES: ' 534 write(io_unit,10) '# Tr[1-eps^{-1}] is evaluated in the Lanczos basis of eps ' 535 write(io_unit,10) '# Tr[1-eps_m^{-1}] is evaluated in the Lanczos basis of eps_m ' 536 write(io_unit,10) '# Tr[eps_m^{-1}-eps^{-1}] is evaluated in the Lanczos basis of eps ' 537 write(io_unit,10) '# ' 538 write(io_unit,10) '# ' 539 write(io_unit,10) '# k Tr[1-eps^{-1}] Tr[1-eps_m^{-1}] Tr[eps_m^{-1}-eps^{-1}]' 540 write(io_unit,10) '#----------------------------------------------------------------------------' 541 flush(io_unit) 542 543 544 io_unit2 = get_unit() 545 open(file='RPA_ENERGY.dat',status=files_status_new,unit=io_unit2) 546 write(io_unit2,10) '#----------------------------------------------------------------------------' 547 write(io_unit2,10) '# ' 548 write(io_unit2,10) '# RPA TOTAL ENERGY ' 549 write(io_unit2,10) '# ========================================== ' 550 write(io_unit2,10) '# ' 551 write(io_unit2,10) '# It can be shown that the correlation energy, within the RPA, is given ' 552 write(io_unit2,10) '# by: ' 553 write(io_unit2,10) '# E_c = int_0^{infty} dw/(2pi) Tr[ ln(eps{iw)}+1-eps(iw)] ' 554 write(io_unit2,10) '# ' 555 write(io_unit2,10) '# As a gauge of what can be expected as far as convergence is concerned, ' 556 write(io_unit2,10) '# the following will be printed. ' 557 write(io_unit2,10) '# ' 558 write(io_unit2,10) '# I_1 = Tr[ ln(eps) + 1 - eps ] ' 559 write(io_unit2,10) '# I_2 = Tr[ ln(eps_m) + 1 - eps_m ] ' 560 write(io_unit2,10) '# I_3 = Tr[ ln(eps^{-1}_m . eps ) + eps_m - eps ] ' 561 write(io_unit2,10) '# ' 562 write(io_unit2,10) '# ' 563 write(io_unit2,10) '# ' 564 write(io_unit2,10) '# k I_1 I_2 I_3 ' 565 write(io_unit2,10) '#----------------------------------------------------------------------------' 566 flush(io_unit2) 567 568 569 ! Iterate every 10 values of k max, or else the linear algebra gets too expensive... 570 do k = 4, kmax, 4 571 572 ABI_MALLOC(sub_Lbasis_exact,(npw_k,k*nseeds)) 573 ABI_MALLOC(sub_Lbasis_model,(npw_k,k*nseeds)) 574 575 576 ABI_MALLOC(eig_exact,(k*nseeds)) 577 ABI_MALLOC(eig_model,(k*nseeds)) 578 ABI_MALLOC(dummy,(k*nseeds,k*nseeds)) 579 ABI_MALLOC(dummy2,(k*nseeds,k*nseeds)) 580 581 sub_Lbasis_exact(:,:) = Lbasis_exact(:,1:k*nseeds) 582 sub_Lbasis_model(:,:) = Lbasis_model(:,1:k*nseeds) 583 584 ! Diagonalize the epsilon matrix, which is banded 585 call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_exact(:,:,1:k),beta_exact(:,:,1:k),sub_Lbasis_exact,eig_exact,.false.) 586 call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_model(:,:,1:k),beta_model(:,:,1:k),sub_Lbasis_model,eig_model,.false.) 587 588 tr_eps_1 = sum(one-one/eig_exact(:)) 589 tr_eps_2 = sum(one-one/eig_model(:)) 590 591 tr_eps_3 = -sum(one/eig_exact(:)) 592 593 dummy(:,:) = model_epsilon_matrix(1:k*nseeds, 1:k*nseeds) 594 595 call driver_invert_positive_definite_hermitian_matrix(dummy,k*nseeds) 596 597 do lm = 1, k*nseeds 598 599 tr_eps_3 = tr_eps_3 + dble(dummy(lm,lm)) 600 end do 601 602 603 write(io_unit,20) k, tr_eps_1, tr_eps_2, tr_eps_3 604 flush(io_unit) 605 606 607 608 tr_eps_1 = sum(log(eig_exact(:))+one-eig_exact(:)) 609 tr_eps_2 = sum(log(eig_model(:))+one-eig_model(:)) 610 611 dummy2(:,:) = zero 612 tr_eps_3 = zero 613 do lm = 1, k*nseeds 614 dummy2(lm,lm) = eig_exact(lm) 615 tr_eps_3 = tr_eps_3 + dble(model_epsilon_matrix(lm,lm)) - dble(eig_exact(lm)) 616 end do 617 618 ABI_MALLOC(dummy3,(k*nseeds,k*nseeds)) 619 call ZGEMM( 'N', & ! Hermitian conjugate the first array 620 'N', & ! Leave second array as is 621 k*nseeds, & ! the number of rows of the matrix op( A ) 622 k*nseeds, & ! the number of columns of the matrix op( B ) 623 k*nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B ) 624 cmplx_1, & ! alpha constant 625 dummy2, & ! matrix A 626 k*nseeds, & ! LDA 627 dummy, & ! matrix B 628 k*nseeds, & ! LDB 629 cmplx_0, & ! beta constant 630 dummy3, & ! matrix C 631 k*nseeds) ! LDC 632 633 dummy2(:,:) = dummy3(:,:) 634 ABI_FREE(dummy3) 635 636 ! find eigenvalues 637 !call heevd(dummy2, eig_exact) 638 639 lwork = k*nseeds+1 640 lrwork = k*nseeds 641 liwork = 1 642 643 ABI_MALLOC(work,(lwork)) 644 ABI_MALLOC(rwork,(lrwork)) 645 ABI_MALLOC(iwork,(liwork)) 646 647 call zheevd('N', 'U',k*nseeds, dummy2, k*nseeds, eig_exact, work, lwork, rwork, lrwork, iwork, liwork, info) 648 if ( info /= 0) then 649 debug_unit = get_unit() 650 write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log' 651 652 open(debug_unit,file=trim(debug_filename),status='unknown') 653 654 write(debug_unit,'(A)') '*************************************************************************************' 655 write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZHEEVD(1), gwls_GenerateEpsilon' 656 write(debug_unit,'(A)') '*************************************************************************************' 657 658 close(debug_unit) 659 660 end if 661 662 663 664 665 666 tr_eps_3 = tr_eps_3 + sum(log(eig_exact)) 667 668 write(io_unit2,20) k, tr_eps_1, tr_eps_2, tr_eps_3 669 flush(io_unit2) 670 671 672 673 ABI_FREE(work) 674 ABI_FREE(rwork) 675 ABI_FREE(iwork) 676 677 678 679 ABI_FREE(sub_Lbasis_exact) 680 ABI_FREE(sub_Lbasis_model) 681 682 ABI_FREE(eig_exact) 683 ABI_FREE(eig_model) 684 ABI_FREE(dummy) 685 ABI_FREE(dummy2) 686 end do 687 688 close(io_unit) 689 close(io_unit2) 690 691 call cpu_time(time2) 692 time = time2-time1 693 write(timing_string,'(A)') "Time to compute the TRACES of the Dielectric Matrices: " 694 call write_timing_log(timing_string,time) 695 696 ABI_FREE(Lbasis_exact) 697 ABI_FREE(Lbasis_model) 698 ABI_FREE(alpha_exact) 699 ABI_FREE(beta_exact ) 700 ABI_FREE(alpha_model) 701 ABI_FREE(beta_model ) 702 ABI_FREE(model_epsilon_matrix) 703 704 705 706 10 format(A) 707 20 format(I5,3ES24.16) 708 709 end subroutine Driver_GeneratePrintDielectricEigenvalues
m_gwls_GenerateEpsilon/GeneratePrintDielectricEigenvalues [ Functions ]
[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]
NAME
GeneratePrintDielectricEigenvalues
FUNCTION
.
INPUTS
OUTPUT
SOURCE
159 subroutine GeneratePrintDielectricEigenvalues(epsilon_matrix_function,nseeds,kmax,output_filename,Lbasis,alpha,beta) 160 !---------------------------------------------------------------------- 161 ! This routine computes the Lanczos approximate representation of the 162 ! implicit dielectic operator and then diagonalizes the banded 163 ! Lanczos matrix. 164 !---------------------------------------------------------------------- 165 interface 166 subroutine epsilon_matrix_function(v_out,v_in,l) 167 use defs_basis 168 169 integer, intent(in) :: l 170 complex(dp), intent(out) :: v_out(l) 171 complex(dp), intent(in) :: v_in(l) 172 173 end subroutine epsilon_matrix_function 174 end interface 175 176 integer, intent(in) :: nseeds, kmax 177 178 character(*), intent(in) :: output_filename 179 180 181 complex(dpc), intent(out) :: Lbasis(npw_k,nseeds*kmax) 182 complex(dpc), intent(out) :: alpha(nseeds,nseeds,kmax) 183 complex(dpc), intent(out) :: beta (nseeds,nseeds,kmax) 184 185 186 ! local variables 187 188 189 complex(dpc),allocatable :: seeds(:,:) 190 191 complex(dpc),allocatable :: Lbasis_diag(:,:) 192 193 194 real(dp), allocatable :: psik(:,:) 195 real(dp), allocatable :: psir(:,:,:,:) 196 197 real(dp), allocatable :: epsilon_eigenvalues(:) 198 199 200 integer :: mpi_communicator 201 integer :: io_unit 202 integer :: lmax 203 integer :: l 204 integer :: ir1, ir2, ir3 205 integer :: n1, n2, n3 206 207 real(dp) :: R, G 208 real(dp) :: sigma_R, sigma_G 209 real(dp) :: x, y, z 210 211 real(dp),allocatable :: G_array(:) 212 real(dp),allocatable :: R_array(:,:,:) 213 214 logical :: debug 215 216 ! ************************************************************************* 217 218 219 debug = .false. 220 lmax = kmax*nseeds 221 mpi_communicator = mpi_enreg%comm_bandfft 222 !Create seeds 223 ABI_MALLOC(seeds,(npw_k,nseeds)) 224 call get_seeds(first_seed, nseeds, seeds) 225 226 ! compute the Lanczos basis 227 ABI_MALLOC(Lbasis_diag,(npw_k,lmax)) 228 ABI_MALLOC(epsilon_eigenvalues,(lmax)) 229 230 ABI_MALLOC(psik,(2,npw_k)) 231 ABI_MALLOC(psir,(2,n4,n5,n6)) 232 ABI_MALLOC(G_array,(npw_k)) 233 ABI_MALLOC(R_array,(n4,n5,n6)) 234 235 psir = zero 236 R_array = zero 237 238 n1 = n4-1 239 n2 = n5-1 240 n3 = n6 241 242 ! Generate the Lanczos basis and banded eigenvalue representation 243 call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k, seeds,alpha,beta,Lbasis) 244 245 Lbasis_diag = Lbasis 246 247 ! Diagonalize the epsilon matrix, which is banded 248 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis_diag,epsilon_eigenvalues,debug) 249 250 call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,lmax,npw_k,Lbasis_diag,epsilon_eigenvalues) 251 252 io_unit = get_unit() 253 open(file=output_filename,status=files_status_new,unit=io_unit) 254 write(io_unit,10) '#----------------------------------------------------------------------------' 255 write(io_unit,10) '# ' 256 write(io_unit,10) '# Partial eigenvalues ' 257 write(io_unit,10) '# ========================================== ' 258 write(io_unit,10) '# ' 259 write(io_unit,10) '# Tabulate the eigenvalues of the dielectic matrix, as well as some ' 260 write(io_unit,10) '# information regarding the eigenstates. ' 261 write(io_unit,10) '# ' 262 write(io_unit,10) '# ' 263 write(io_unit,10) '# definitions: ' 264 write(io_unit,10) '# l index of the eigenvalue ' 265 write(io_unit,10) '# eig eigenvalue ' 266 write(io_unit,10) '# ' 267 write(io_unit,10) '# ' 268 write(io_unit,10) '# (the following vectors are expressed in crystal units) ' 269 write(io_unit,10) '# ' 270 write(io_unit,10) '# R = < l | |r| | l > ' 271 write(io_unit,10) '# sigma_R = sqrt{< l | (|r|-R)^2 | l > } ' 272 write(io_unit,10) '# ' 273 write(io_unit,10) '# G = < l | |G| | l > ' 274 write(io_unit,10) '# sigma_G = sqrt{< l | (|G|-G)^2 | l > } ' 275 write(io_unit,10) '# ' 276 write(io_unit,10) '# ' 277 write(io_unit,10) '# l eig r sigma_r G sigma_G ' 278 write(io_unit,10) '#----------------------------------------------------------------------------' 279 flush(io_unit) 280 281 282 G_array(:) = kg_k(1,:)**2+ kg_k(2,:)**2+ kg_k(3,:)**2 283 284 G_array(:) = sqrt(G_array(:)) 285 286 R_array = zero 287 288 do ir3=1,n3 289 290 if (ir3 <= n3/2 ) then 291 z = (one*ir3)/(one*n3) 292 else 293 z = (one*ir3)/(one*n3)-one 294 end if 295 296 do ir2=1,n2 297 298 if (ir2 <= n2/2 ) then 299 y = (one*ir2)/(one*n2) 300 else 301 y = (one*ir2)/(one*n2)-one 302 end if 303 304 do ir1=1,n1 305 306 if (ir1 <= n1/2 ) then 307 x = (one*ir1)/(one*n1) 308 else 309 x = (one*ir1)/(one*n1)-one 310 end if 311 312 313 R_array(ir1,ir2,ir3) = sqrt(x**2+y**2+z**2) 314 end do 315 end do 316 end do 317 318 319 320 do l=1, lmax 321 322 psik(1,:) = dble (Lbasis_diag(:,l)) 323 psik(2,:) = dimag(Lbasis_diag(:,l)) 324 325 call g_to_r(psir ,psik) 326 327 328 G = sum(G_array(:)*(psik(1,:)**2+psik(2,:)**2)) 329 R = sum(R_array(:,:,:)*(psir(1,:,:,:)**2+psir(2,:,:,:)**2) )*ucvol/nfft 330 331 sigma_G = sqrt(sum((G_array(:) -G)**2*(psik(1,:)**2 +psik(2,:)**2))) 332 sigma_R = sqrt(sum((R_array(:,:,:)-R)**2*(psir(1,:,:,:)**2+psir(2,:,:,:)**2))*ucvol/nfft) 333 334 335 336 write(io_unit,20) l, epsilon_eigenvalues(l), R,sigma_R, G,sigma_G 337 338 end do 339 340 close(io_unit) 341 342 ABI_FREE(seeds) 343 ABI_FREE(Lbasis_diag) 344 ABI_FREE(psik) 345 ABI_FREE(psir) 346 ABI_FREE(G_array) 347 ABI_FREE(R_array) 348 349 ABI_FREE(epsilon_eigenvalues) 350 351 352 10 format(A) 353 20 format(I5,ES24.16,4F12.8) 354 355 end subroutine GeneratePrintDielectricEigenvalues