TABLE OF CONTENTS
ABINIT/m_paw_mkrho [ Modules ]
NAME
m_paw_mkrho
FUNCTION
This module contains routines used to compute PAW density on the real space fine grid.
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (MT, JWZ) 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 MODULE m_paw_mkrho 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 29 use defs_abitypes, only : MPI_type 30 use m_time, only : timab 31 use m_pawang, only : pawang_type 32 use m_pawrad, only : pawrad_type,pawrad_deducer0 33 use m_pawtab, only : pawtab_type,pawtab_get_lsize 34 use m_paw_sphharm, only : initylmr 35 use m_pawfgrtab, only : pawfgrtab_type,pawfgrtab_init,pawfgrtab_free 36 use m_pawrhoij, only : pawrhoij_type,pawrhoij_copy,pawrhoij_free_unpacked, & 37 & pawrhoij_nullify,pawrhoij_free,pawrhoij_symrhoij 38 use m_pawfgr, only : pawfgr_type 39 use m_paw_nhat, only : pawmknhat,nhatgrid 40 use m_paral_atom, only : get_my_atmtab,free_my_atmtab 41 use m_fourier_interpol, only : transgrid 42 43 use m_sort, only : sort_dp 44 use m_splines, only : spline,splint 45 use m_io_tools, only : open_file 46 use m_geometry, only : xred2xcart 47 use m_pptools, only : printxsf 48 use m_fft, only : fourdp 49 50 implicit none 51 52 private 53 54 !public procedures. 55 public :: pawmkrho ! Build PAW electronic density on fine grid, including compensation charge density 56 public :: denfgr ! Build complete PAW electronic density on fine grid, including on-site contributions 57 58 CONTAINS !========================================================================================
m_paw_mkrho/denfgr [ Functions ]
[ Top ] [ m_paw_mkrho ] [ Functions ]
NAME
denfgr
FUNCTION
Construct complete electron density on fine grid, by removing nhat and adding PAW corrections
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) gmet(3,3)=reciprocal space metric tensor in bohr**-2. mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms my_natom=number of atoms treated by current processor natom= number of atoms in cell nattyp(ntypat)= # atoms of each type. ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv) nhat(pawfgr%nfft,nspden)= compensation charge density used in PAW nspinor=Number of spinor components nsppol=Number of independent spin components. nspden= number of spin densities ntypat= number of types of atoms in the cell pawfgr <type(pawfgr_type)>= data about the fine grid pawrad(ntypat) <type(pawrad_type)>= radial mesh data for each type of atom pawrhoij(natom) <type(pawrhoij_type)>= rho_ij data for each atom pawtab(ntypat) <type(pawtab_type)>= PAW functions around each type of atom rhor(pawfgr%nfft,nspden)= input density ($\tilde{n}+\hat{n}$ in PAW case) rprimd(3,3)=dimensional primitive translations for real space (bohr) typat(natom)= list of atom types ucvol=unit cell volume (bohr**3) xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
rhor_paw(pawfgr%nfft,nspden)= full electron density on the fine grid
NOTES
In PAW calculations, the valence density present in rhor includes the
compensation charge density $\hat{n}$, and also doesn't include the on-site
PAW contributions. For post-processing and proper visualization it is necessary
to use the full electronic density, which is what this subroutine constructs.
Specifically, it removes $\hat{n}$ from rhor, and also computes the on-site PAW
terms. This is nothing other than the proper PAW treatment of the density
operator $|\mathbf{r}\rangle\langle\mathbf{r}|$, and yields the formula
$$\tilde{n}+\sum_{ij}\rho_ij\left[\varphi_i(\mathbf{r})\varphi_j(\mathbf{r})-
\tilde{\varphi}_i(\mathbf{r})\tilde{\varphi}_j(\mathbf{r})\right]$$
Notice that this formula is expressed on the fine grid, and requires
interpolating the PAW radial functions onto this grid, as well as calling
initylmr in order to get the angular functions on the grid points.
SOURCE
332 subroutine denfgr(atindx1,gmet,spaceComm_in,my_natom,natom,nattyp,ngfft,nhat,nspinor,nsppol,nspden,ntypat, & 333 & pawfgr,pawrad,pawrhoij,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,rhor_nt_one,rprimd,typat,ucvol,xred,& 334 & abs_n_tilde_nt_diff,znucl,mpi_atmtab,comm_atom) ! Optional arguments 335 336 !Arguments ------------------------------------ 337 !scalars 338 integer,intent(in) :: my_natom,natom,nspden,ntypat,prtvol,nsppol,nspinor 339 integer,optional,intent(in) :: comm_atom 340 real(dp),intent(in) :: ucvol 341 type(pawfgr_type),intent(in) :: pawfgr 342 !arrays 343 integer,intent(in) :: spaceComm_in 344 integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),typat(natom) 345 integer,optional,target,intent(in) :: mpi_atmtab(:) 346 real(dp),intent(in) :: gmet(3,3),nhat(pawfgr%nfft,nspden) 347 real(dp),intent(in) :: rhor(pawfgr%nfft,nspden),rprimd(3,3) 348 real(dp),intent(inout) :: xred(3,natom) 349 real(dp),intent(out) :: rhor_paw(pawfgr%nfft,nspden) 350 real(dp),intent(out) :: rhor_n_one(pawfgr%nfft,nspden) 351 real(dp),intent(out) :: rhor_nt_one(pawfgr%nfft,nspden) 352 real(dp),optional,intent(out) :: abs_n_tilde_nt_diff(nspden) 353 real(dp),optional,intent(in) :: znucl(ntypat) 354 type(pawrad_type),intent(in) :: pawrad(ntypat) 355 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 356 type(pawtab_type),target,intent(in) :: pawtab(ntypat) 357 358 !Local variables------------------------------- 359 !scalars 360 integer,parameter :: master=0 361 integer :: delta,iatom,ierr,ifgd,ifftsph,inl,inrm,ipsang,irhoij 362 integer :: ispden,itypat,il,im,ilm,iln,ilmn 363 integer :: jl,jlm,jln,jm,j0lmn,jlmn 364 integer :: klmn,my_comm_atom,my_start_indx,my_end_indx 365 integer :: nfgd,nnl,normchoice,nprocs,optcut,optgr0,optgr1,optgr2 366 integer :: optrad,option,my_rank,remainder,tmp_unt 367 real(dp) :: phj,phi,rR,tphj,tphi,ybcbeg,ybcend 368 logical :: my_atmtab_allocated,paral_atom 369 character(len=500) :: message 370 !arrays 371 integer,allocatable :: l_size_atm(:),nrm_ifftsph(:) 372 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 373 integer,pointer :: my_atmtab(:) 374 real(dp) :: ylmgr(3,3,0) 375 real(dp) :: yvals(4),xcart(3,natom) 376 real(dp),allocatable :: diag(:),nrm(:),phigrd(:,:),tphigrd(:,:),ylm(:,:),ypp(:) 377 real(dp),allocatable :: phi_at_zero(:),tphi_at_zero(:) 378 real(dp),allocatable :: rhor_tmp(:,:),tot_rhor(:) 379 character(len=fnlen) :: xsf_fname 380 type(pawfgrtab_type) :: local_pawfgrtab(my_natom) 381 382 ! ************************************************************************ 383 384 DBG_ENTER("COLL") 385 386 if (my_natom>0) then 387 ABI_CHECK(pawrhoij(1)%qphase==1,'denfgr not supposed to be called with qphase/=1!') 388 end if 389 390 !Set up parallelism over atoms (compatible only with band-FFT parallelism) 391 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 392 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 393 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 394 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 395 & my_natom_ref=my_natom) 396 397 !MG: FIXME It won't work if atom-parallelism is used 398 !but even the loop over atoms below should be rewritten in this case. 399 400 !use a local copy of pawfgrtab to make sure we use the correction in the paw spheres 401 !the usual pawfgrtab uses r_shape which may not be the same as r_paw 402 if (my_natom>0) then 403 if (paral_atom) then 404 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,typat,mpi_atmtab=my_atmtab) 405 call pawfgrtab_init(local_pawfgrtab,pawrhoij(1)%qphase,l_size_atm,nspden,typat,& 406 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 407 else 408 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,typat) 409 call pawfgrtab_init(local_pawfgrtab,pawrhoij(1)%qphase,l_size_atm,nspden,typat) 410 end if 411 ABI_FREE(l_size_atm) 412 end if 413 414 !Note: call to nhatgrid: comm_fft not used because FFT parallelism 415 !is done manually below 416 optcut = 1 ! use rpaw to construct local_pawfgrtab 417 optgr0 = 0; optgr1 = 0; optgr2 = 0 ! dont need gY terms locally 418 optrad = 1 ! do store r-R 419 if (paral_atom) then 420 call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,& 421 & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,rprimd,typat,ucvol,xred,& 422 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 423 else 424 call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,& 425 & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,rprimd,typat,ucvol,xred) 426 end if 427 !now local_pawfgrtab is ready to use 428 429 !Initialise output arrays. 430 rhor_paw=zero; rhor_n_one=zero; rhor_nt_one=zero 431 432 !Initialise and check parallell execution 433 my_rank = xmpi_comm_rank(spaceComm_in) 434 nprocs = xmpi_comm_size(spaceComm_in) 435 436 437 !loop over atoms in cell 438 do iatom = 1, my_natom 439 itypat = pawrhoij(iatom)%itypat 440 indlmn => pawtab(itypat)%indlmn 441 nfgd = local_pawfgrtab(iatom)%nfgd ! number of points in the fine grid for this PAW sphere 442 nnl = pawtab(itypat)%basis_size ! number of nl elements in PAW basis 443 444 ! Division of fine grid points among processors 445 if (nprocs==1) then ! Make sure everything runs with one proc 446 write(message,'(a)') ' In denfgr - number of processors: 1' 447 call wrtout(std_out,message,'COLL') 448 write(message,'(a)') ' Calculation of PAW density done in serial' 449 call wrtout(std_out,message,'COLL') 450 write(message,'(a,I6)') ' Number of fine grid points:',nfgd 451 call wrtout(std_out,message,'COLL') 452 my_start_indx = 1 453 my_end_indx = nfgd 454 else ! Divide up the fine grid points among the processors 455 write(message,'(a,I4)') ' In denfgr - number of processors: ',nprocs 456 call wrtout(std_out,message,'COLL') 457 write(message,'(a)') ' Calculation of PAW density done in parallel' 458 call wrtout(std_out,message,'COLL') 459 write(message,'(a,I6)') ' Number of fine grid points:',nfgd 460 call wrtout(std_out,message,'COLL') 461 ! Divide the fine grid points among the processors 462 delta = int(floor(real(nfgd)/real(nprocs))) 463 remainder = nfgd-nprocs*delta 464 my_start_indx = 1+my_rank*delta 465 my_end_indx = (my_rank+1)*delta 466 ! Divide the remainder points among the processors 467 ! by shuffling indices 468 if ((my_rank+1)>remainder) then 469 my_start_indx = my_start_indx + remainder 470 my_end_indx = my_end_indx + remainder 471 else 472 my_start_indx = my_start_indx + my_rank 473 my_end_indx = my_end_indx + my_rank + 1 474 end if 475 if (prtvol>9) then 476 write(message,'(a,I6)') ' My index Starts at: ',my_start_indx 477 call wrtout(std_out,message,'PERS') 478 write(message,'(a,I6)') ' Ends at: ',my_end_indx 479 call wrtout(std_out,message,'PERS') 480 write(message,'(a,I6)') ' # pts: ',my_end_indx+1-my_start_indx 481 call wrtout(std_out,message,'PERS') 482 end if 483 end if 484 485 write(message,'(a,I3,a,I3)') ' Entered loop for atom: ',iatom,' of:',natom 486 call wrtout(std_out,message,'PERS') 487 488 ! obtain |r-R| values on fine grid 489 ABI_MALLOC(nrm,(nfgd)) 490 do ifgd=1, nfgd 491 nrm(ifgd) = sqrt(dot_product(local_pawfgrtab(iatom)%rfgd(:,ifgd),local_pawfgrtab(iatom)%rfgd(:,ifgd))) 492 end do ! these are the |r-R| values 493 494 ! compute Ylm for each r-R vector. 495 ! ---- 496 ipsang = 1 + (pawtab(itypat)%l_size - 1)/2 ! recall l_size=2*l_max+1 497 ABI_MALLOC(ylm,(ipsang*ipsang,nfgd)) 498 option = 1 ! compute Ylm(r-R) for vectors 499 normchoice = 1 ! use computed norms of input vectors 500 call initylmr(ipsang,normchoice,nfgd,nrm,option,local_pawfgrtab(iatom)%rfgd,ylm,ylmgr) 501 502 ! in order to do spline fits, the |r-R| data must be sorted 503 ! ---- 504 ABI_MALLOC(nrm_ifftsph,(nfgd)) 505 nrm_ifftsph(:) = local_pawfgrtab(iatom)%ifftsph(:) ! copy of indices of points, to be rearranged by sort_dp 506 call sort_dp(nfgd,nrm,nrm_ifftsph,tol8) ! sort the nrm points, keeping track of which goes where 507 508 ! now make spline fits of phi and tphi onto the fine grid around the atom 509 ! ---- 510 ABI_MALLOC(phigrd,(nfgd,nnl)) 511 ABI_MALLOC(tphigrd,(nfgd,nnl)) 512 ABI_MALLOC(phi_at_zero,(nnl)) 513 ABI_MALLOC(tphi_at_zero,(nnl)) 514 ABI_MALLOC(ypp,(pawtab(itypat)%mesh_size)) 515 ABI_MALLOC(diag,(pawtab(itypat)%mesh_size)) 516 517 do inl = 1, nnl 518 519 ! spline phi onto points 520 ypp(:) = zero; diag(:) = zero; ybcbeg = zero; ybcend = zero; 521 call spline(pawrad(itypat)%rad,pawtab(itypat)%phi(:,inl),pawtab(itypat)%mesh_size,ybcbeg,ybcend,ypp) 522 call splint(pawtab(itypat)%mesh_size,pawrad(itypat)%rad,pawtab(itypat)%phi(:,inl),ypp,nfgd,nrm,phigrd(:,inl)) 523 524 ! next splint tphi onto points 525 ypp(:) = zero; diag(:) = zero; ybcbeg = zero; ybcend = zero; 526 call spline(pawrad(itypat)%rad,pawtab(itypat)%tphi(:,inl),pawtab(itypat)%mesh_size,ybcbeg,ybcend,ypp) 527 call splint(pawtab(itypat)%mesh_size,pawrad(itypat)%rad,pawtab(itypat)%tphi(:,inl),ypp,nfgd,nrm,tphigrd(:,inl)) 528 529 ! Find out the value of the basis function at zero using extrapolation 530 yvals = zero 531 ! Extrapolate only if this is an s-state (l=0) 532 if (indlmn(1,inl)==0) then 533 yvals(2:4) = pawtab(itypat)%phi(2:4,inl)/pawrad(itypat)%rad(2:4) 534 call pawrad_deducer0(yvals,4,pawrad(itypat)) 535 write(std_out,*) 'phi_at_zero: ',yvals(1),' from:',yvals(2:4) 536 end if 537 phi_at_zero(inl) = yvals(1) 538 539 yvals = zero 540 ! Extrapolate only if this is an s-state (l=0) 541 if (indlmn(1,inl)==0) then 542 yvals(2:4) = pawtab(itypat)%tphi(2:4,inl)/pawrad(itypat)%rad(2:4) 543 call pawrad_deducer0(yvals,4,pawrad(itypat)) 544 write(std_out,*) 'tphi_at_zero: ',yvals(1),' from:',yvals(2:4) 545 end if 546 tphi_at_zero(inl) = yvals(1) 547 548 end do ! end loop over nnl basis functions 549 ABI_FREE(ypp) 550 ABI_FREE(diag) 551 552 ! loop over basis elements for this atom 553 ! because we have to store things like <phi|r'><r'|phi>-<tphi|r'><r'|tphi> at each point of the 554 ! fine grid, there is no integration, and hence no simplifications of the Y_lm's. That's why 555 ! we have to loop through the basis elements in exhaustive detail, rather than just a loop over 556 ! lmn2_size or something comparable. 557 ! ---- 558 if (prtvol>9) then 559 write(message,'(a,I3)') ' Entering j-loop over basis elements for atom:',iatom 560 call wrtout(std_out,message,'PERS') 561 end if 562 563 do jlmn=1,pawtab(itypat)%lmn_size 564 565 if (prtvol>9) then 566 write(message,'(2(a,I3))') ' Element:',jlmn,' of:',pawtab(itypat)%lmn_size 567 call wrtout(std_out,message,'PERS') 568 end if 569 570 jl=indlmn(1,jlmn) 571 jm=indlmn(2,jlmn) 572 jlm=indlmn(4,jlmn) 573 jln=indlmn(5,jlmn) 574 j0lmn=jlmn*(jlmn-1)/2 575 576 if (prtvol>9) then 577 write(message,'(a,I3)') ' Entering i-loop for j:',jlmn 578 call wrtout(std_out,message,'PERS') 579 end if 580 581 do ilmn=1,jlmn 582 583 if (prtvol>9) then 584 write(message,'(2(a,I3))') ' Element:',ilmn,' of:',jlmn 585 call wrtout(std_out,message,'PERS') 586 end if 587 588 il=indlmn(1,ilmn) 589 im=indlmn(2,ilmn) 590 iln=indlmn(5,ilmn) 591 ilm=indlmn(4,ilmn) 592 klmn=j0lmn+ilmn 593 594 if (prtvol>9) then 595 write(message,'(a)') ' Entering loop over nonzero elems of rhoij' 596 call wrtout(std_out,message,'PERS') 597 end if 598 599 ! Loop over non-zero elements of rhoij 600 do irhoij=1,pawrhoij(iatom)%nrhoijsel 601 if (klmn==pawrhoij(iatom)%rhoijselect(irhoij)) then ! rho_ij /= 0 for this klmn 602 603 do ifgd=my_start_indx, my_end_indx ! loop over fine grid points in current PAW sphere 604 ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid 605 606 ! have to retrieve the spline point to use since these were sorted 607 do inrm=1, nfgd 608 if(nrm_ifftsph(inrm) == ifftsph) exit ! have found nrm point corresponding to nfgd point 609 end do ! now inrm is the index of the sorted nrm vector to use 610 611 ! avoid division by zero 612 if(nrm(inrm) > zero) then 613 rR = nrm(inrm) ! value of |r-R| in the following 614 ! recall that <r|phi>=u(r)*Slm(r^)/r 615 phj = phigrd(inrm,jln)*ylm(jlm,ifgd)/rR 616 phi = phigrd(inrm,iln)*ylm(ilm,ifgd)/rR 617 tphj = tphigrd(inrm,jln)*ylm(jlm,ifgd)/rR 618 tphi = tphigrd(inrm,iln)*ylm(ilm,ifgd)/rR 619 else 620 ! use precalculated <r|phi>=u(r)*Slm(r^)/r at r=0 621 phj = phi_at_zero(jln)*ylm(jlm,ifgd) 622 phi = phi_at_zero(iln)*ylm(ilm,ifgd) 623 tphj = tphi_at_zero(jln)*ylm(jlm,ifgd) 624 tphi = tphi_at_zero(iln)*ylm(ilm,ifgd) 625 end if ! check if |r-R| = 0 626 627 do ispden=1,nspden 628 if (pawrhoij(iatom)%cplex_rhoij == 1) then 629 rhor_paw(ifftsph,ispden) = rhor_paw(ifftsph,ispden) + & 630 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*(phj*phi - tphj*tphi) 631 632 rhor_n_one(ifftsph,ispden) = rhor_n_one(ifftsph,ispden) + & 633 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*phj*phi 634 635 rhor_nt_one(ifftsph,ispden) = rhor_nt_one(ifftsph,ispden) + & 636 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*tphj*tphi 637 else 638 rhor_paw(ifftsph,ispden) = rhor_paw(ifftsph,ispden) + & 639 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*(phj*phi - tphj*tphi) 640 641 rhor_n_one(ifftsph,ispden) = rhor_n_one(ifftsph,ispden) + & 642 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*phj*phi 643 644 rhor_nt_one(ifftsph,ispden) = rhor_nt_one(ifftsph,ispden) + & 645 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*tphj*tphi 646 end if ! end check on cplex rhoij 647 648 end do ! end loop over nsdpen 649 end do ! end loop over nfgd 650 end if ! end selection on rhoij /= 0 651 end do ! end loop over non-zero rhoij 652 end do ! end loop over ilmn atomic basis states 653 end do ! end loop over jlmn atomic basis states 654 655 ABI_FREE(nrm) 656 ABI_FREE(nrm_ifftsph) 657 ABI_FREE(phigrd) 658 ABI_FREE(tphigrd) 659 ABI_FREE(ylm) 660 ABI_FREE(phi_at_zero) 661 ABI_FREE(tphi_at_zero) 662 end do ! Loop on atoms 663 664 !MPI sum on each node the different contributions to the PAW densities. 665 call xmpi_sum(rhor_paw,spaceComm_in,ierr) 666 call xmpi_sum(rhor_n_one,spaceComm_in,ierr) 667 call xmpi_sum(rhor_nt_one,spaceComm_in,ierr) 668 if (paral_atom) then 669 call xmpi_sum(rhor_paw,my_comm_atom,ierr) 670 call xmpi_sum(rhor_n_one,my_comm_atom,ierr) 671 call xmpi_sum(rhor_nt_one,my_comm_atom,ierr) 672 end if 673 674 call wrtout(std_out,' *** Partial contributions to PAW rhor summed ***','PERS') 675 call xmpi_barrier(spaceComm_in) 676 677 !Add the plane-wave contribution \tilde{n} and remove \hat{n} 678 !BE careful here since the storage mode of rhoij and rhor is different. 679 select case (nspinor) 680 case (1) 681 if (nsppol==1) then 682 rhor_paw = rhor_paw + rhor - nhat 683 else ! Spin-polarised case: rhor_paw contains rhor_paw(spin_up,spin_down) but we need rhor_paw(total,spin_up) 684 ABI_MALLOC(tot_rhor,(pawfgr%nfft)) 685 ! 686 ! AE rhor 687 tot_rhor(:) = SUM(rhor_paw,DIM=2) 688 rhor_paw(:,2) = rhor_paw(:,1) 689 rhor_paw(:,1) = tot_rhor 690 rhor_paw = rhor_paw + rhor - nhat 691 ! 692 ! onsite AE rhor 693 tot_rhor(:) = SUM(rhor_n_one,DIM=2) 694 rhor_n_one(:,2) = rhor_n_one(:,1) 695 rhor_n_one(:,1) = tot_rhor 696 ! 697 ! onsite PS rhor 698 tot_rhor(:) = SUM(rhor_nt_one,DIM=2) 699 rhor_nt_one(:,2) = rhor_nt_one(:,1) 700 rhor_nt_one(:,1) = tot_rhor 701 702 ABI_FREE(tot_rhor) 703 end if 704 705 case (2) 706 ! * if nspden==4, rhor contains (n^11, n^22, Re[n^12], Im[n^12]. 707 ! Storage mode for rhoij is different, See pawaccrhoij. 708 ABI_ERROR("nspinor 2 not coded") 709 case default 710 write(message,'(a,i0)')" Wrong value for nspinor=",nspinor 711 ABI_ERROR(message) 712 end select 713 714 !if (prtvol>9) then ! Check normalisation 715 !write(message,'(a,F8.4)') ' PAWDEN - NORM OF DENSITY: ',SUM(rhor_paw(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3)) 716 !call wrtout(std_out,message,'COLL') 717 !end if 718 719 if (present(abs_n_tilde_nt_diff).AND.present(znucl)) then 720 ABI_MALLOC(rhor_tmp,(pawfgr%nfft,nspden)) 721 do ispden=1,nspden 722 rhor_tmp(:,ispden) = zero 723 do iatom=1,my_natom 724 do ifgd=1,local_pawfgrtab(iatom)%nfgd ! loop over fine grid points in current PAW sphere 725 ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid 726 rhor_tmp(ifftsph,ispden) = rhor(ifftsph,ispden) - nhat(ifftsph,ispden) & 727 & - rhor_nt_one(ifftsph,ispden) 728 end do !ifgd 729 end do ! iatom 730 end do ! ispden 731 if (paral_atom) then 732 call xmpi_sum(rhor_tmp,my_comm_atom,ierr) 733 end if 734 735 if (my_rank==master) then 736 do ispden=1,nspden 737 ! Write to xsf file 738 call xred2xcart(natom,rprimd,xcart,xred) 739 write(xsf_fname,'(a,I0,a)') 'N_tilde_onsite_diff_sp',ispden,'.xsf' 740 if (open_file(xsf_fname,message, unit=tmp_unt,status='unknown',form='formatted') /= 0) then 741 ABI_ERROR(message) 742 end if 743 call printxsf(ngfft(1),ngfft(2),ngfft(3),rhor_tmp(:,ispden),rprimd,& 744 & (/zero,zero,zero/),natom,ntypat,typat,xcart,znucl,tmp_unt,0) 745 close(tmp_unt) 746 abs_n_tilde_nt_diff(ispden) = SUM(ABS(rhor_tmp(:,ispden)))/pawfgr%nfft 747 write(message,'(4(a),F16.9,2(a,I0),a)') ch10,' Wrote xsf file with \tilde{n}-\tilde{n}^1.',ch10,& 748 & ' Value of norm |\tilde{n}-\tilde{n}^1|:',& 749 & abs_n_tilde_nt_diff(ispden),' spin: ',ispden,' of ',nspden,ch10 750 call wrtout(std_out,message,'COLL') 751 end do 752 end if 753 ABI_FREE(rhor_tmp) 754 755 else if ((present(abs_n_tilde_nt_diff).AND.(.NOT.present(znucl))) & 756 & .OR.(.NOT.present(abs_n_tilde_nt_diff).AND.(present(znucl)))) then 757 write(message,'(a)') ' Both abs_n_tilde_nt_diff *and* znucl must be passed',ch10,& 758 & 'to denfgr for |\tilde{n}-\tilde{n}^1| norm evaluation.' 759 ABI_ERROR(message) 760 end if 761 762 call pawfgrtab_free(local_pawfgrtab) 763 764 !Destroy atom table used for parallelism 765 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 766 767 call xmpi_barrier(spaceComm_in) 768 769 DBG_EXIT("COLL") 770 771 end subroutine denfgr
m_paw_mkrho/pawmkrho [ Functions ]
[ Top ] [ m_paw_mkrho ] [ Functions ]
NAME
pawmkrho
FUNCTION
PAW only: Build total pseudo (compensated) density (\tild_rho + \hat_rho) Build compensation charge density (\hat_rho) Build occupation matrix (packed storage)
INPUTS
compute_rhor_rhog: if 1: set the computation of rhor and rhog in addition to the compensating charge.
if 0: compute only the compensating charge
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
1 for GS calculations
gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
indsym(4,nsym,natom)=indirect indexing array for atom labels
ipert=index of perturbation if pawrhoij is a pertubed rhoij
no meaning for ground-state calculations (should be 0)
idir=direction of atomic displacement (in case of atomic displ. perturb.)
mpi_enreg=information about MPI parallelization
my_natom=number of atoms treated by current processor
natom=number of atoms in cell
nspden=number of spin-density components
nsym=number of symmetry elements in space group
ntypat=number of types of atoms in unit cell.
paral_kgb=option for (kpt,g vectors,bands) parallelism
pawang <type(pawang_type)>=angular mesh discretization and related data
pawang_sym <type(pawang_type)>=angular data used for symmetrization
optional parameter only needed for RF calculations
pawfgr <type(paw_fgr_type)>=fine rectangular grid parameters
pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
pawprtvol=control print volume and debugging output for PAW
pawrhoij0(natom) <type(pawrhoij_type)>= GS paw rhoij occupancies and related data (used only if ipert>0)
optional parameter only needed for RF calculations
pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
qphon(3)=wavevector of the phonon (RF only)
rhopsg(2,pawfgr%nfftc)= pseudo density given on the coarse grid in reciprocal space
rhopsr(pawfgr%nfftc,nspden)= pseudo density given on the coarse grid in real space
rprimd(3,3)=real space primitive translations.
symafm(nsym)=(anti)ferromagnetic part of symmetry operations
symrec(3,3,nsym)=symmetries of group in terms of operations on
reciprocal space primitive translations
typat(natom)=type for each atom
ucvol=volume of the unit cell
xred(3,natom)= reduced atomic coordinates
OUTPUT
compch_fft=compensation charge inside spheres integrated over fine fft grid pawnhat(pawfgr%nfft,nspden)=compensation charge density on fine rectangular grid (optional argument) rhog(2,pawfgr%nfft)= compensated pseudo density given on the fine grid in reciprocal space This output is optional rhor(pawfgr%nfft,nspden)= compensated pseudo density given on the fine grid in real space
SIDE EFFECTS
pawrhoij(my_natom)= PAW occupancies
At input : values at previous step in packed storage (pawrhoij()%rhoijp)
At output: values (symmetrized) in packed storage (pawrhoij()%rhoijp)
pawrhoij_unsym(:)= unsymmetrized PAW occupancies
At input : values (unsymmetrized) in unpacked storage (pawrhoij()%rhoij_)
At output: values in unpacked storage (pawrhoij()%rhoij_) are destroyed
NOTES
pawrhoij and pawrhoij_unsym can be identical (refer to the same pawrhoij datastructure). They should be different only if pawrhoij is distributed over atomic sites (in that case pawrhoij_unsym should not be distributed over atomic sites).
SOURCE
132 subroutine pawmkrho(compute_rhor_rhog,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,& 133 & my_natom,natom,nspden,nsym,ntypat,paral_kgb,pawang,pawfgr,pawfgrtab,pawprtvol,& 134 & pawrhoij,pawrhoij_unsym,& 135 & pawtab,qphon,rhopsg,rhopsr,rhor,rprimd,symafm,symrec,typat,ucvol,usewvl,xred,& 136 & pawang_sym,pawnhat,pawnhatgr,pawrhoij0,rhog) ! optional arguments 137 138 !Arguments ------------------------------------ 139 !scalars 140 integer,intent(in) :: compute_rhor_rhog,cplex,idir,ipert,my_natom,natom,nspden,nsym,ntypat,paral_kgb,pawprtvol 141 integer,intent(in) :: usewvl 142 real(dp),intent(in) :: ucvol 143 real(dp),intent(out) :: compch_fft 144 type(MPI_type),intent(in) :: mpi_enreg 145 type(pawang_type),intent(in) :: pawang 146 type(pawang_type),intent(in),optional :: pawang_sym 147 type(pawfgr_type),intent(in) :: pawfgr 148 !arrays 149 integer,intent(in) :: indsym(4,nsym,natom) 150 integer,intent(in) :: symafm(nsym),symrec(3,3,nsym),typat(natom) 151 real(dp),intent(in) :: gprimd(3,3),qphon(3),rprimd(3,3),xred(3,natom) 152 real(dp),intent(inout),target,optional :: pawnhat(cplex*pawfgr%nfft,nspden) !vz_i 153 real(dp),intent(inout),target,optional :: pawnhatgr(:,:,:) !vz_i 154 real(dp),intent(inout) :: rhor(cplex*pawfgr%nfft,nspden*compute_rhor_rhog) 155 real(dp),intent(out),optional :: rhog(2,pawfgr%nfft*compute_rhor_rhog) 156 real(dp),intent(inout) :: rhopsg(2,pawfgr%nfftc*compute_rhor_rhog),rhopsr(cplex*pawfgr%nfftc,nspden*compute_rhor_rhog) 157 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 158 type(pawrhoij_type),intent(inout),target :: pawrhoij(:) 159 type(pawrhoij_type),intent(inout) :: pawrhoij_unsym(:) 160 type(pawrhoij_type),intent(in),target,optional :: pawrhoij0(my_natom) 161 type(pawtab_type),intent(in) :: pawtab(ntypat) 162 163 !Local variables------------------------------- 164 !scalars 165 integer :: choice,ider,izero,option 166 character(len=500) :: msg 167 !arrays 168 real(dp) :: tsec(2) 169 real(dp),target :: rhodum(0,0,0) 170 real(dp),pointer :: pawnhat_ptr(:,:) 171 real(dp),pointer :: pawnhatgr_ptr(:,:,:) 172 type(pawrhoij_type),pointer :: pawrhoij_ptr(:),pawrhoij0_ptr(:) 173 174 ! *********************************************************************** 175 176 DBG_ENTER("COLL") 177 178 call timab(556,1,tsec) 179 180 !Compatibility tests 181 if (size(pawrhoij_unsym)>0) then 182 if (pawrhoij_unsym(1)%use_rhoij_==0) then 183 msg=' rhoij_ field must be allocated in pawrhoij_unsym !' 184 ABI_BUG(msg) 185 end if 186 end if 187 if (ipert>0.and.(.not.present(pawrhoij0))) then 188 msg=' pawrhoij0 must be present when ipert>0 !' 189 ABI_BUG(msg) 190 end if 191 192 !Symetrize PAW occupation matrix and store it in packed storage 193 call timab(557,1,tsec) 194 option=1;choice=1 195 if (present(pawang_sym)) then 196 call pawrhoij_symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,& 197 & natom,nsym,ntypat,option,pawang_sym,pawprtvol,pawtab,rprimd,symafm,& 198 & symrec,typat,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 199 & qphon=qphon) 200 else 201 call pawrhoij_symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,& 202 & natom,nsym,ntypat,option,pawang,pawprtvol,pawtab,rprimd,symafm,& 203 & symrec,typat,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 204 & qphon=qphon) 205 end if 206 call pawrhoij_free_unpacked(pawrhoij_unsym) 207 call timab(557,2,tsec) 208 209 !In somes cases (parallelism), has to distribute the PAW occupation matrix 210 if (size(pawrhoij)==natom.and.(my_natom/=natom)) then 211 ABI_MALLOC(pawrhoij_ptr,(my_natom)) 212 call pawrhoij_nullify(pawrhoij_ptr) 213 call pawrhoij_copy(pawrhoij,pawrhoij_ptr,& 214 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom, & 215 & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.) 216 else 217 pawrhoij_ptr=>pawrhoij 218 end if 219 220 !Compute compensation charge density 221 ider=0;izero=0 222 if (present(pawnhat)) then 223 pawnhat_ptr => pawnhat 224 else 225 ABI_MALLOC(pawnhat_ptr,(pawfgr%nfft,nspden)) 226 end if 227 if (present(pawnhatgr)) then 228 pawnhatgr_ptr => pawnhatgr 229 else 230 pawnhatgr_ptr => rhodum 231 end if 232 if (present(pawrhoij0)) then 233 pawrhoij0_ptr => pawrhoij0 234 else 235 pawrhoij0_ptr => pawrhoij_ptr 236 end if 237 238 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,& 239 & pawfgr%nfft,pawfgr%ngfft,ider,nspden,ntypat,pawang,pawfgrtab,& 240 & pawnhatgr_ptr,pawnhat_ptr,pawrhoij_ptr,pawrhoij0_ptr,pawtab,qphon,rprimd,ucvol,usewvl,xred,& 241 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 242 & comm_fft=mpi_enreg%comm_fft,paral_kgb=paral_kgb,me_g0=mpi_enreg%me_g0,& 243 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 244 245 if (compute_rhor_rhog/=0) then 246 ! Transfer pseudo density from coarse grid to fine grid 247 if(usewvl==0) then 248 call transgrid(cplex,mpi_enreg,nspden,+1,1,0,paral_kgb,pawfgr,rhopsg,rhodum,rhopsr,rhor) 249 end if 250 251 ! Add pseudo density and compensation charge density (on fine grid) 252 rhor(:,:)=rhor(:,:)+pawnhat_ptr(:,:) 253 254 ! Compute compensated pseudo density in reciprocal space 255 if (present(rhog)) then 256 call fourdp(cplex,rhog,rhor(:,1),-1,mpi_enreg,pawfgr%nfft,1,pawfgr%ngfft,0) 257 end if 258 end if 259 260 !Free temporary memory spaces 261 if (.not.present(pawnhat)) then 262 ABI_FREE(pawnhat_ptr) 263 end if 264 if (size(pawrhoij)==natom.and.(my_natom/=natom)) then 265 call pawrhoij_free(pawrhoij_ptr) 266 ABI_FREE(pawrhoij_ptr) 267 end if 268 nullify(pawnhat_ptr) 269 nullify(pawnhatgr_ptr) 270 nullify(pawrhoij_ptr) 271 272 call timab(556,2,tsec) 273 274 DBG_EXIT("COLL") 275 276 end subroutine pawmkrho