TABLE OF CONTENTS
- ABINIT/m_paw_dfptnl
- ABINIT/paw_dfptnl_accrhoij
- m_paw_dfptnl/paw_dfptnl_energy
- m_paw_dfptnl/paw_dfptnl_xc
ABINIT/m_paw_dfptnl [ Modules ]
NAME
m_paw_dfptnl
FUNCTION
This module contains several routines used to compute PAW contributions to a 3rd-order energy or 2nd-order PAW occupancies.
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (LB) 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
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_paw_dfptnl 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_xmpi 29 30 use m_pawang, only : pawang_type 31 use m_pawrad, only : pawrad_type,simp_gen 32 use m_pawtab, only : pawtab_type 33 use m_paw_an, only : paw_an_type 34 use m_pawrhoij, only : pawrhoij_type 35 use m_pawcprj, only : pawcprj_type 36 use m_paw_denpot, only : pawdensities 37 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 38 39 implicit none 40 41 private 42 43 !public procedures. 44 public :: paw_dfptnl_energy ! Compute the XC PAW on-site contributions to a 3rd-order energy 45 public :: paw_dfptnl_xc ! Compute a contribution of the 3rd-derivative of XC energy of ONE PAW sphere 46 public :: paw_dfptnl_accrhoij ! Accumulate the 2nd order PAW quantities rhoij^(2) 47 48 CONTAINS !========================================================================================
ABINIT/paw_dfptnl_accrhoij [ Functions ]
NAME
paw_dfptnl_accrhoij
FUNCTION
Accumulate the 2nd order PAW quantities rhoij^(2) (augmentation occupancies) This routine is similar to pawaccrhoij.F90 but is implemented independently in order to not overload the original routine.
INPUTS
atindx(natom)=index table for atoms (sorted-->random), inverse of atindx.
cplex: if 1, WFs (or 1st-order WFs) are REAL, if 2, COMPLEX
cwaveprj0_pert1(natom,nspinor) = wave function at given n,k projected with non-local projectors:
cwaveprj0%cp =<p_i|
cwaveprj0%dcp(1)=<p_i^(pert1)|Cnk>
cwaveprj0_pert2(natom,nspinor) = wave function at given n,k projected with non-local projectors:
cwaveprj0%cp =<p_i|Cnk>
cwaveprj0%dcp(1)=<p_i^(pert1)|Cnk>
cwaveprj1_pert12(natom,nspinor)= 1st order wave function at given n,k projected with non-local projectors:
cwaveprj1%cp =<p_i|Cnk^(pert1)>
cwaveprj1%dcp(1)=<p_i^(pert2)|Cnk^(pert1)>
cwaveprj1_pert21(natom,nspinor)= 1st order wave function at given n,k projected with non-local projectors:
cwaveprj1%cp =<p_i|Cnk^(pert2)>
cwaveprj1%dcp(1)=<p_i^(pert1)|Cnk^(pert2)>
ipert1=index of the first perturbation
ipert2=index of the second perturbation
isppol=index of current spin component
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
nspinor=number of spinorial components (on current proc)
occ_k=occupation number for current band n,k
wtk_k=weight assigned to current k-point
SIDE EFFECTS
pawrhoij(natom) <type(pawrhoij_type)>= 2-nd order paw rhoij occupancies and related data On output, has been updated with the contribution of current n,k pawrhoij(:)%rhoij_(lmn2_size,nspden) (non symetrized)
SOURCE
505 subroutine paw_dfptnl_accrhoij(atindx,cplex,cwaveprj0_pert1,cwaveprj0_pert2,& 506 & cwaveprj1_pert12,cwaveprj1_pert21,ipert1,ipert2,isppol,my_natom,natom,& 507 & nspinor,occ_k,pawrhoij,wtk_k,& 508 & comm_atom,mpi_atmtab ) ! optional (parallelism) 509 510 !Arguments --------------------------------------------- 511 !scalars 512 integer,intent(in) :: cplex,ipert1,ipert2,isppol,my_natom,natom,nspinor 513 integer,optional,intent(in) :: comm_atom 514 real(dp),intent(in) :: occ_k,wtk_k 515 !arrays 516 integer,intent(in) :: atindx(natom) 517 integer,optional,target,intent(in) :: mpi_atmtab(:) 518 type(pawcprj_type),intent(in) :: cwaveprj0_pert1(natom,nspinor),cwaveprj0_pert2(natom,nspinor) 519 type(pawcprj_type),intent(in) :: cwaveprj1_pert12(natom,nspinor),cwaveprj1_pert21(natom,nspinor) 520 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 521 522 !Local variables --------------------------------------- 523 !scalars 524 integer :: cplex_rhoij,iatm,iatom,iatom1,ilmn,iplex,j0lmn,jlmn,klmn,klmn_im,klmn_re 525 integer :: my_comm_atom,ncpgr 526 logical :: compute_impart,compute_impart_cplex 527 logical :: my_atmtab_allocated,paral_atom 528 real(dp) :: ro11_im,ro11_re,weight 529 character(len=500) :: message 530 !arrays 531 integer,pointer :: my_atmtab(:) 532 real(dp) :: cpi0(2,nspinor),d1cpi0(2,nspinor),d2cpi0(2,nspinor) 533 real(dp) :: cpj0(2,nspinor),d1cpj0(2,nspinor),d2cpj0(2,nspinor) 534 real(dp) :: cpi1(2,nspinor),d2cpi1(2,nspinor) 535 real(dp) :: cpj1(2,nspinor),d2cpj1(2,nspinor) 536 real(dp) :: cpi2(2,nspinor),d1cpi2(2,nspinor) 537 real(dp) :: cpj2(2,nspinor),d1cpj2(2,nspinor) 538 539 ! *********************************************************************** 540 541 DBG_ENTER("COLL") 542 543 if (my_natom==0) return 544 545 ncpgr=1 546 if (ipert1<0.or.ipert1>natom+2.or.ipert2<0.or.ipert2>natom+2) then 547 message = 'paw_dfptnl_accrhoij: Necessary conditions on ipert1 or ipert2: 0<=ipert<=natom+2' 548 ABI_BUG(message) 549 end if 550 if (pawrhoij(1)%qphase/=1) then 551 message="paw_dfptnl_accrhoij not supposed to be called with q/=0!" 552 ABI_BUG(message) 553 end if 554 555 !Set up parallelism over atoms 556 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 557 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 558 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 559 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 560 & my_natom_ref=my_natom) 561 562 weight=wtk_k*occ_k 563 if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1) weight=half*weight 564 565 !! ================================================================== 566 !! === Accumulate (n,k) contribution to partial 2nd-order rhoij === 567 !! ================================================================== 568 569 compute_impart=(pawrhoij(1)%cplex_rhoij==2) 570 compute_impart_cplex=((pawrhoij(1)%cplex_rhoij==2).and.(cplex==2)) 571 572 ! NOT USED FOR PAWRHO21! => only for PAWRHO2 (full second derivative) 573 !!Accumulate : < Psi^(pert1) | p_i^(0) > < p_j^(0) | Psi^(pert2) > 574 !! + < Psi^(pert2) | p_i^(0) > < p_j^(0) | Psi^(pert1) > 575 ! if (nspinor==1) then 576 ! do iatom=1,my_natom 577 ! iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 578 ! iatm=atindx(iatom1) 579 ! cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 580 ! do jlmn=1,pawrhoij(iatom)%lmn_size 581 ! j0lmn=jlmn*(jlmn-1)/2 582 ! cpj1(1:2,1)=cwaveprj1_pert12(iatm,1)%cp(1:2,jlmn) ! < p_j^(0) | Psi^(pert1) > 583 ! cpj2(1:2,1)=cwaveprj1_pert21(iatm,1)%cp(1:2,jlmn) ! < p_j^(0) | Psi^(pert2) > 584 ! do ilmn=1,jlmn 585 ! klmn=j0lmn+ilmn 586 ! klmn_re=cplex_rhoij*(klmn-1)+1 587 ! cpi1(1:2,1)=cwaveprj1_pert12(iatm,1)%cp(1:2,ilmn) ! < p_i^(0) | Psi^(pert1) > 588 ! cpi2(1:2,1)=cwaveprj1_pert21(iatm,1)%cp(1:2,ilmn) ! < p_i^(0) | Psi^(pert2) > 589 ! ro11_re=zero 590 ! do iplex=1,cplex 591 ! ro11_re=ro11_re+cpi1(iplex,1)*cpj2(iplex,1)+cpj1(iplex,1)*cpi2(iplex,1) 592 ! end do 593 ! pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 594 ! if (compute_impart_cplex) then 595 ! klmn_im=klmn_re+1 596 ! ro11_im= cpi1(1,1)*cpj2(2,1)-cpi1(2,1)*cpj2(1,1) 597 ! ro11_im=ro11_im+cpj1(1,1)*cpi2(2,1)-cpj1(2,1)*cpi2(1,1) 598 ! pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 599 ! end if 600 ! end do 601 ! end do 602 ! end do 603 ! else ! nspinor=2 604 ! ABI_BUG("paw_dfptnl_accrhoij is not implemented for nspinor=2") 605 ! end if 606 607 !Accumulate : < Psi^(pert1) | p_i^(pert2) > < p_j^(0) | Psi^(0) > 608 ! + < Psi^(pert1) | p_i^(0) > < p_j^(pert2) | Psi^(0) > 609 ! + < Psi^(0) | p_i^(pert2) > < p_j^(0) | Psi^(pert1) > 610 ! + < Psi^(0) | p_i^(0) > < p_j^(pert2) | Psi^(pert1) > 611 if (ipert2>0.and.ipert2<=natom) then 612 if (nspinor==1) then 613 do iatom=1,my_natom 614 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 615 iatm=atindx(iatom1) 616 if (iatom/=ipert2) cycle ! To move atom "ipert2" does not change projectors of other atoms 617 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 618 do jlmn=1,pawrhoij(iatom)%lmn_size 619 j0lmn=jlmn*(jlmn-1)/2 620 cpj0(1:2,1) =cwaveprj0_pert2 (iatm,1)% cp(1:2, jlmn) ! < p_j^(0) | Psi^(0) > 621 d2cpj0(1:2,1)=cwaveprj0_pert2 (iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert2) | Psi^(0) > 622 cpj1(1:2,1) =cwaveprj1_pert12(iatm,1)% cp(1:2, jlmn) ! < p_j^(0) | Psi^(pert1) > 623 d2cpj1(1:2,1)=cwaveprj1_pert12(iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert2) | Psi^(pert1) > 624 do ilmn=1,jlmn 625 klmn=j0lmn+ilmn 626 klmn_re=cplex_rhoij*(klmn-1)+1 627 cpi0(1:2,1) =cwaveprj0_pert2 (iatm,1)% cp(1:2, ilmn) ! < p_i^(0) | Psi^(0) > 628 d2cpi0(1:2,1)=cwaveprj0_pert2 (iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(0) > 629 cpi1(1:2,1) =cwaveprj1_pert12(iatm,1)% cp(1:2, ilmn) ! < p_i^(0) | Psi^(pert1) > 630 d2cpi1(1:2,1)=cwaveprj1_pert12(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(pert1) > 631 ro11_re=zero 632 do iplex=1,cplex 633 ro11_re=ro11_re+d2cpi1(iplex,1)* cpj0(iplex,1) 634 ro11_re=ro11_re+ cpi1(iplex,1)*d2cpj0(iplex,1) 635 ro11_re=ro11_re+d2cpi0(iplex,1)* cpj1(iplex,1) 636 ro11_re=ro11_re+ cpi0(iplex,1)*d2cpj1(iplex,1) 637 end do 638 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 639 if (compute_impart_cplex) then 640 klmn_im=klmn_re+1 641 ro11_im= d2cpi1(1,1)* cpj0(2,1)-d2cpi1(2,1)* cpj0(1,1) 642 ro11_im=ro11_im+ cpi1(1,1)*d2cpj0(2,1)- cpi1(2,1)*d2cpj0(1,1) 643 ro11_im=ro11_im+d2cpi0(1,1)* cpj1(2,1)-d2cpi0(2,1)* cpj1(1,1) 644 ro11_im=ro11_im+ cpi0(1,1)*d2cpj1(2,1)- cpi0(2,1)*d2cpj1(1,1) 645 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 646 end if 647 end do 648 end do 649 end do 650 else ! nspinor=2 651 ABI_BUG("paw_dfptnl_accrhoij is not implemented for nspinor=2") 652 end if 653 end if 654 655 !Accumulate : < Psi^(pert2) | p_i^(pert1) > < p_j^(0) | Psi^(0) > 656 ! + < Psi^(pert2) | p_i^(0) > < p_j^(pert1) | Psi^(0) > 657 ! + < Psi^(0) | p_i^(pert1) > < p_j^(0) | Psi^(pert2) > 658 ! + < Psi^(0) | p_i^(0) > < p_j^(pert1) | Psi^(pert2) > 659 if (ipert1>0.and.ipert1<=natom) then 660 if (nspinor==1) then 661 do iatom=1,my_natom 662 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 663 iatm=atindx(iatom1) 664 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 665 if (iatom/=ipert1) cycle ! To move atom "ipert1" does not change projectors of other atoms 666 do jlmn=1,pawrhoij(iatom)%lmn_size 667 j0lmn=jlmn*(jlmn-1)/2 668 cpj0(1:2,1) =cwaveprj0_pert1 (iatm,1)% cp(1:2, jlmn) ! < p_j^(0) | Psi^(0) > 669 d1cpj0(1:2,1)=cwaveprj0_pert1 (iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert1) | Psi^(0) > 670 cpj2(1:2,1) =cwaveprj1_pert21(iatm,1)% cp(1:2, jlmn) ! < p_j^(0) | Psi^(pert2) > 671 d1cpj2(1:2,1)=cwaveprj1_pert21(iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert1) | Psi^(pert2) > 672 do ilmn=1,jlmn 673 klmn=j0lmn+ilmn 674 klmn_re=cplex_rhoij*(klmn-1)+1 675 cpi0(1:2,1) =cwaveprj0_pert1 (iatm,1)% cp(1:2, ilmn) ! < p_i^(0) | Psi^(0) > 676 d1cpi0(1:2,1)=cwaveprj0_pert1 (iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(0) > 677 cpi2(1:2,1) =cwaveprj1_pert21(iatm,1)% cp(1:2, ilmn) ! < p_i^(0) | Psi^(pert2) > 678 d1cpi2(1:2,1)=cwaveprj1_pert21(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(pert2) > 679 ro11_re=zero 680 do iplex=1,cplex 681 ro11_re=ro11_re+d1cpi2(iplex,1)* cpj0(iplex,1) 682 ro11_re=ro11_re+ cpi2(iplex,1)*d1cpj0(iplex,1) 683 ro11_re=ro11_re+d1cpi0(iplex,1)* cpj2(iplex,1) 684 ro11_re=ro11_re+ cpi0(iplex,1)*d1cpj2(iplex,1) 685 end do 686 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 687 if (compute_impart_cplex) then 688 klmn_im=klmn_re+1 689 ro11_im= d1cpi2(1,1)* cpj0(2,1)-d1cpi2(2,1)* cpj0(1,1) 690 ro11_im=ro11_im+ cpi2(1,1)*d1cpj0(2,1)- cpi2(2,1)*d1cpj0(1,1) 691 ro11_im=ro11_im+d1cpi0(1,1)* cpj2(2,1)-d1cpi0(2,1)* cpj2(1,1) 692 ro11_im=ro11_im+ cpi0(1,1)*d1cpj2(2,1)- cpi0(2,1)*d1cpj2(1,1) 693 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 694 end if 695 end do 696 end do 697 end do 698 else ! nspinor=2 699 message="paw_dfptnl_accrhoij is not implemented for nspinor=2!" 700 ABI_BUG(message) 701 end if 702 end if 703 ! End 704 705 !Accumulate : < Psi^(0) | p_i^(pert1) > < p_j^(pert2) | Psi^(0) > 706 ! + < Psi^(0) | p_i^(pert2) > < p_j^(pert1) | Psi^(0) > 707 if (ipert1>0.and.ipert1<=natom.and.ipert2>0.and.ipert2<=natom) then 708 if (nspinor==1) then 709 do iatom=1,my_natom 710 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 711 iatm=atindx(iatom1) 712 if (iatom/=ipert1.or.iatom/=ipert2) cycle ! To move atom "ipert" does not change projectors of other atoms 713 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 714 do jlmn=1,pawrhoij(iatom)%lmn_size 715 j0lmn=jlmn*(jlmn-1)/2 716 d1cpj0(1:2,1)=cwaveprj0_pert1(iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert1) | Psi^(0) > 717 d2cpj0(1:2,1)=cwaveprj0_pert2(iatm,1)%dcp(1:2,1,jlmn) ! < p_j^(pert2) | Psi^(0) > 718 do ilmn=1,jlmn 719 klmn=j0lmn+ilmn 720 klmn_re=cplex_rhoij*(klmn-1)+1 721 d1cpi0(1:2,1)=cwaveprj0_pert1(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(0) > 722 d2cpi0(1:2,1)=cwaveprj0_pert2(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(0) > 723 ro11_re=zero 724 do iplex=1,cplex 725 ro11_re=ro11_re+d1cpi0(iplex,1)*d2cpj0(iplex,1)+d2cpi0(iplex,1)*d1cpj0(iplex,1) 726 end do 727 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 728 if (compute_impart_cplex) then 729 klmn_im=klmn_re+1 730 ro11_im= d1cpi0(1,1)*d2cpj0(2,1)-d1cpi0(2,1)*d2cpj0(1,1) 731 ro11_im=ro11_im+d2cpi0(1,1)*d1cpj0(2,1)-d2cpi0(2,1)*d1cpj0(1,1) 732 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 733 end if 734 end do 735 end do 736 end do 737 else ! nspinor=2 738 message="paw_dfptnl_accrhoij is not implemented for nspinor=2!" 739 ABI_BUG(message) 740 end if 741 end if 742 743 !Destroy atom table used for parallelism 744 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 745 746 DBG_EXIT("COLL") 747 748 end subroutine paw_dfptnl_accrhoij
m_paw_dfptnl/paw_dfptnl_energy [ Functions ]
[ Top ] [ m_paw_dfptnl ] [ Functions ]
NAME
paw_dfptnl_energy
FUNCTION
Compute the XC PAW on-site contributions to a 3rd-order energy. It is equal to: E_onsite= \sum_at [ E_at(kxc,rho1,rho2,rho3) - E_at(tkxc,trho1,trho2,trho3) ] where E_at(...) is computed in paw_dfptnl_xc.F90. The atomic densities are computed from pawrhoij1,pawrhoij2 and pawrhoij3. This routine is similar to pawdfptenergy.F90 but is implemented independently in order to not overload the original routine. LDA ONLY - USE THE DENSITY OVER A WHOLE SPHERICAL GRID (r,theta,phi)
INPUTS
ixc= choice of exchange-correlation scheme 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=total number of atoms in cell ntypat=number of types of atoms in unit cell. paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh paw_an1(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh This corresponds to (j1) perturbation pawang <type(pawang_type)>=paw angular mesh and related data pawprtvol=control print volume and debugging output for PAW pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawrhoij_1-2-3(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
OUTPUT
d3exc= real and imaginary parts of the contribution to the third derivative of the total energy
SIDE EFFECTS
SOURCE
91 subroutine paw_dfptnl_energy(d3exc,ixc,my_natom,natom,ntypat,& 92 & paw_an0,pawang,pawprtvol,pawrad,& 93 & pawrhoij_1,pawrhoij_2,pawrhoij_3,& 94 & pawtab,pawxcdev,& 95 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 96 97 !Arguments --------------------------------------------- 98 !scalars 99 integer,intent(in) :: ixc,my_natom,natom,ntypat 100 integer,intent(in) :: pawprtvol,pawxcdev 101 integer,optional,intent(in) :: comm_atom 102 type(pawang_type),intent(in) :: pawang 103 !arrays 104 integer,optional,target,intent(in) :: mpi_atmtab(:) 105 real(dp),intent(out) :: d3exc(2) 106 type(paw_an_type),intent(in) :: paw_an0(my_natom) 107 type(pawrad_type),intent(in) :: pawrad(ntypat) 108 type(pawrhoij_type),intent(in) :: pawrhoij_1(my_natom) 109 type(pawrhoij_type),intent(in) :: pawrhoij_2(my_natom) 110 type(pawrhoij_type),intent(in) :: pawrhoij_3(my_natom) 111 type(pawtab_type),intent(in) :: pawtab(ntypat) 112 113 !Local variables --------------------------------------- 114 !scalars 115 integer :: cplex_1,cplex_2,cplex_3,iatom,iatom_tot,itypat 116 integer :: lm_size_all,mesh_size,my_comm_atom,npts,nspden,nzlmopt 117 integer :: opt_compch,usecore,usetcore,usexcnhat 118 logical :: my_atmtab_allocated,paral_atom 119 real(dp) :: compch,d3exc1_iat(2) 120 character(len=500) :: msg 121 !arrays 122 integer,pointer :: my_atmtab(:) 123 logical,allocatable :: lmselect_1(:),lmselect_2(:),lmselect_3(:),lmselect_tmp(:) 124 real(dp),allocatable :: nhat1_1(:,:,:),rho1_1(:,:,:),trho1_1(:,:,:) 125 real(dp),allocatable :: nhat1_2(:,:,:),rho1_2(:,:,:),trho1_2(:,:,:) 126 real(dp),allocatable :: nhat1_3(:,:,:),rho1_3(:,:,:),trho1_3(:,:,:) 127 128 ! ************************************************************************* 129 130 DBG_ENTER("COLL") 131 132 nzlmopt = 0 ! compute all LM-moments of the density and use all LM-moments 133 134 if (pawxcdev/=0) then 135 msg="paw_dfptnl_energy is not implemented for pawxcdev/=0" 136 ABI_BUG(msg) 137 end if 138 if (my_natom>0) then 139 if (pawrhoij_1(1)%qphase/=1.or.pawrhoij_2(1)%qphase/=1.or.pawrhoij_3(1)%qphase/=1) then 140 msg="paw_dfptnl_energy not supposed to be called with q/=0!" 141 ABI_BUG(msg) 142 end if 143 end if 144 145 !Set up parallelism over atoms 146 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 147 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 148 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 149 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 150 151 !!Various inits 152 opt_compch=0; !optvxc=1;optexc=3 153 usecore=0;usetcore=0 ! This is true for phonons and Efield pert. 154 usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat) 155 156 npts=pawang%angl_size 157 158 d3exc = zero 159 160 !================ Loop on atomic sites ======================= 161 do iatom=1,my_natom 162 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 163 164 itypat=pawrhoij_1(iatom)%itypat 165 mesh_size=pawtab(itypat)%mesh_size 166 nspden=pawrhoij_1(iatom)%nspden 167 cplex_1=pawrhoij_1(iatom)%cplex_rhoij 168 cplex_2=pawrhoij_2(iatom)%cplex_rhoij 169 cplex_3=pawrhoij_3(iatom)%cplex_rhoij 170 lm_size_all=paw_an0(iatom)%lm_size 171 172 ABI_MALLOC(lmselect_tmp,(lm_size_all)) 173 lmselect_tmp(:)=.true. 174 175 ! Compute on-site 1st-order densities (pert1) 176 ABI_MALLOC(lmselect_1,(lm_size_all)) 177 lmselect_1(:)=paw_an0(iatom)%lmselect(:) 178 ABI_MALLOC(rho1_1,(cplex_1*mesh_size,lm_size_all,nspden)) 179 ABI_MALLOC(trho1_1,(cplex_1*mesh_size,lm_size_all,nspden)) 180 ABI_MALLOC(nhat1_1,(cplex_1*mesh_size,lm_size_all,nspden*usexcnhat)) 181 call pawdensities(compch,cplex_1,iatom_tot,lmselect_tmp,lmselect_1,& 182 & lm_size_all,nhat1_1,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,& 183 & pawrad(itypat),pawrhoij_1(iatom),pawtab(itypat),rho1_1,trho1_1) 184 ! Compute on-site 1st-order densities (pert2) 185 ABI_MALLOC(lmselect_2,(lm_size_all)) 186 lmselect_2(:)=paw_an0(iatom)%lmselect(:) 187 ABI_MALLOC(rho1_2,(cplex_2*mesh_size,lm_size_all,nspden)) 188 ABI_MALLOC(trho1_2,(cplex_2*mesh_size,lm_size_all,nspden)) 189 ABI_MALLOC(nhat1_2,(cplex_2*mesh_size,lm_size_all,nspden*usexcnhat)) 190 call pawdensities(compch,cplex_2,iatom_tot,lmselect_tmp,lmselect_2,& 191 & lm_size_all,nhat1_2,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,& 192 & pawrad(itypat),pawrhoij_2(iatom),pawtab(itypat),rho1_2,trho1_2) 193 ! Compute on-site 1st-order densities (pert3) 194 ABI_MALLOC(lmselect_3,(lm_size_all)) 195 lmselect_3(:)=paw_an0(iatom)%lmselect(:) 196 ABI_MALLOC(rho1_3,(cplex_3*mesh_size,lm_size_all,nspden)) 197 ABI_MALLOC(trho1_3,(cplex_3*mesh_size,lm_size_all,nspden)) 198 ABI_MALLOC(nhat1_3,(cplex_3*mesh_size,lm_size_all,nspden*usexcnhat)) 199 call pawdensities(compch,cplex_3,iatom_tot,lmselect_tmp,lmselect_3,& 200 & lm_size_all,nhat1_3,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,& 201 & pawrad(itypat),pawrhoij_3(iatom),pawtab(itypat),rho1_3,trho1_3) 202 ABI_FREE(lmselect_tmp) 203 204 call paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,paw_an0(iatom)%k3xc1,lm_size_all,& 205 & lmselect_1,lmselect_2,lmselect_3,nhat1_1,nhat1_2,nhat1_3,& 206 & paw_an0(iatom)%nk3xc1,mesh_size,nspden,pawang,pawrad(itypat),& 207 & rho1_1,rho1_2,rho1_3,0) 208 d3exc = d3exc + d3exc1_iat 209 210 call paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,paw_an0(iatom)%k3xct1,lm_size_all,& 211 & lmselect_1,lmselect_2,lmselect_3,nhat1_1,nhat1_2,nhat1_3,& 212 & paw_an0(iatom)%nk3xc1,mesh_size,nspden,pawang,pawrad(itypat),& 213 & trho1_1,trho1_2,trho1_3,usexcnhat) 214 d3exc = d3exc - d3exc1_iat 215 216 ABI_FREE(lmselect_1) 217 ABI_FREE(lmselect_2) 218 ABI_FREE(lmselect_3) 219 ABI_FREE(nhat1_1) 220 ABI_FREE(nhat1_2) 221 ABI_FREE(nhat1_3) 222 ABI_FREE(rho1_1) 223 ABI_FREE(rho1_2) 224 ABI_FREE(rho1_3) 225 ABI_FREE(trho1_1) 226 ABI_FREE(trho1_2) 227 ABI_FREE(trho1_3) 228 229 ! ================ End loop oon atomic sites ======================= 230 end do 231 232 !!Reduction in case of parallelism 233 ! if (paral_atom) then 234 ! call xmpi_sum(delta_energy,my_comm_atom,ierr) 235 ! end if 236 237 !!Destroy atom table used for parallelism 238 ! call free_my_atmtab(my_atmtab,my_atmtab_allocated) 239 240 ! call timab(567,2,tsec) 241 242 DBG_EXIT("COLL") 243 244 end subroutine paw_dfptnl_energy
m_paw_dfptnl/paw_dfptnl_xc [ Functions ]
[ Top ] [ m_paw_dfptnl ] [ Functions ]
NAME
paw_dfptnl_xc
FUNCTION
Compute a contribution of the third derivative of XC energy of ONE PAW sphere Om_a. It is equal to: E_at(kxc,rho1,rho2,rho3) = Int_{Om_a} dr kxc(r) * rho1(r) * rho2(r) * rho3(r) where kxc,rho1,rho2 and rho3 are inputs. This routine is similar to m_pawxc.F90:pawxc_dfpt(...) but is implemented independently in order to not overload the original routine. LDA ONLY - USE THE DENSITY OVER A WHOLE SPHERICAL GRID (r,theta,phi)
INPUTS
cplex_1-2-3= if 1, 1st-order densities are REAL, if 2, COMPLEX d3exc1_iat=third-order derivative to compute ixc= choice of exchange-correlation scheme kxc(nrad,pawang%angl_size,nkxc)=GS xc kernel lm_size=size of density array rhor (see below) lmselect1-2-3(lm_size)=select the non-zero LM-moments of input density rhor1-2-3 nhat1-2-3(cplex_den*nrad,lm_size,nspden)=first-order change of compensation density (total in 1st half and spin-up in 2nd half if nspden=2) nkxc=second dimension of the kxc array nrad=size of radial mesh for densities/potentials (might be different from pawrad%mesh_size) nspden=number of spin-density components option=0 compute both 2nd-order XC energy and 1st-order potential 1 compute only 1st-order XC potential 2 compute only 2nd-order XC energy, XC potential is temporary computed here 3 compute only 2nd-order XC energy, XC potential is input in vxc1(:) pawang <type(pawang_type)>=paw angular mesh and related data pawrad <type(pawrad_type)>=paw radial mesh and related data rhor1-2-3(cplex_den*nrad,lm_size,nspden)=first-order change of density usexcnhat= 0 if compensation density does not have to be used 1 if compensation density has to be used in d2Exc only
OUTPUT
d3exc1_iat = E_at(kxc,rho1,rho2,rho3) (see FUNCTION above)
SIDE EFFECTS
SOURCE
291 subroutine paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,kxc,lm_size,lmselect1,lmselect2,lmselect3,& 292 & nhat1,nhat2,nhat3,nkxc,nrad,nspden,pawang,pawrad,rhor1,rhor2,rhor3,usexcnhat) 293 294 !Arguments ------------------------------------ 295 !scalars 296 integer,intent(in) :: cplex_1,cplex_2,cplex_3,ixc,lm_size,nkxc,nrad,nspden,usexcnhat 297 type(pawang_type),intent(in) :: pawang 298 type(pawrad_type),intent(in) :: pawrad 299 !arrays 300 logical,intent(in) :: lmselect1(lm_size),lmselect2(lm_size),lmselect3(lm_size) 301 real(dp),intent(out) :: d3exc1_iat(2) 302 real(dp),intent(in) :: kxc(nrad,pawang%angl_size,nkxc) 303 real(dp),intent(in) :: nhat1(cplex_1*nrad,lm_size,nspden*((usexcnhat+1)/2)) 304 real(dp),intent(in) :: nhat2(cplex_2*nrad,lm_size,nspden*((usexcnhat+1)/2)) 305 real(dp),intent(in) :: nhat3(cplex_3*nrad,lm_size,nspden*((usexcnhat+1)/2)) 306 real(dp),intent(in) :: rhor1(cplex_1*nrad,lm_size,nspden) 307 real(dp),intent(in) :: rhor2(cplex_2*nrad,lm_size,nspden) 308 real(dp),intent(in) :: rhor3(cplex_3*nrad,lm_size,nspden) 309 310 !Local variables------------------------------- 311 !scalars 312 integer :: ii,ilm,ipts,ispden,lm_size_eff,npts 313 real(dp) :: d3exc1_int,rho1u,rho1d,rho2u,rho2d,rho3u,rho3d 314 character(len=500) :: msg 315 !arrays 316 ! real(dp) :: tsec(2) 317 real(dp),allocatable :: ff(:),rho1arr(:,:),rho2arr(:,:),rho3arr(:,:) 318 319 ! ************************************************************************* 320 321 !---------------------------------------------------------------------- 322 !----- Initializations 323 !---------------------------------------------------------------------- 324 325 npts=pawang%angl_size 326 lm_size_eff=min(lm_size,pawang%ylm_size) 327 328 d3exc1_iat(:) = zero 329 330 !Special case: no XC applied 331 if (ixc==0) then 332 msg='Note that no xc is applied (ixc=0). Returning' 333 ABI_WARNING(msg) 334 return 335 end if 336 337 ABI_MALLOC(rho1arr,(cplex_1*nrad,nspden)) 338 ABI_MALLOC(rho2arr,(cplex_2*nrad,nspden)) 339 ABI_MALLOC(rho3arr,(cplex_3*nrad,nspden)) 340 341 !Restriction : all cplex must be 1 342 if (cplex_1/=1.or.cplex_2/=1.or.cplex_3/=1) then 343 msg='All cplex must be one (for the moment...)' 344 ABI_BUG(msg) 345 end if 346 !Restriction : nspden must be 1 347 ! if (nkxc>1) then 348 ! msg='nkxc must be one (<=> nspden=1) (for the moment...)' 349 ! ABI_BUG(msg) 350 ! end if 351 352 ABI_MALLOC(ff,(nrad)) 353 354 !!---------------------------------------------------------------------- 355 !!----- Loop on the angular part and inits 356 !!---------------------------------------------------------------------- 357 358 !Do loop on the angular part (theta,phi) 359 do ipts=1,npts 360 361 ! Copy the input 1st-order density for this (theta,phi) - PERT1 362 rho1arr(:,:)=zero 363 if (usexcnhat==0) then 364 do ispden=1,nspden 365 do ilm=1,lm_size_eff 366 if (lmselect1(ilm)) rho1arr(:,ispden)=rho1arr(:,ispden) & 367 & +rhor1(:,ilm,ispden)*pawang%ylmr(ilm,ipts) 368 end do 369 end do 370 else 371 do ispden=1,nspden 372 do ilm=1,lm_size_eff 373 if (lmselect1(ilm)) rho1arr(:,ispden)=rho1arr(:,ispden) & 374 & +(rhor1(:,ilm,ispden)+nhat1(:,ilm,ispden))*pawang%ylmr(ilm,ipts) 375 end do 376 end do 377 end if 378 379 ! Copy the input 1st-order density for this (theta,phi) - PERT2 380 rho2arr(:,:)=zero 381 if (usexcnhat==0) then 382 do ispden=1,nspden 383 do ilm=1,lm_size_eff 384 if (lmselect2(ilm)) rho2arr(:,ispden)=rho2arr(:,ispden) & 385 & +rhor2(:,ilm,ispden)*pawang%ylmr(ilm,ipts) 386 end do 387 end do 388 else 389 do ispden=1,nspden 390 do ilm=1,lm_size_eff 391 if (lmselect2(ilm)) rho2arr(:,ispden)=rho2arr(:,ispden) & 392 & +(rhor2(:,ilm,ispden)+nhat2(:,ilm,ispden))*pawang%ylmr(ilm,ipts) 393 end do 394 end do 395 end if 396 397 ! Copy the input 1st-order density for this (theta,phi) - PERT3 398 rho3arr(:,:)=zero 399 if (usexcnhat==0) then 400 do ispden=1,nspden 401 do ilm=1,lm_size_eff 402 if (lmselect3(ilm)) rho3arr(:,ispden)=rho3arr(:,ispden) & 403 & +rhor3(:,ilm,ispden)*pawang%ylmr(ilm,ipts) 404 end do 405 end do 406 else 407 do ispden=1,nspden 408 do ilm=1,lm_size_eff 409 if (lmselect3(ilm)) rho3arr(:,ispden)=rho3arr(:,ispden) & 410 & +(rhor3(:,ilm,ispden)+nhat3(:,ilm,ispden))*pawang%ylmr(ilm,ipts) 411 end do 412 end do 413 end if 414 415 ! ---------------------------------------------------------------------- 416 ! ----- Accumulate and store 3nd-order change of XC energy 417 ! ---------------------------------------------------------------------- 418 419 if (cplex_1==1.and.cplex_2==1.and.cplex_3==1) then ! all cplex are 1 : 420 if (nspden==1) then 421 ff(:)=kxc(:,ipts,1)*rho1arr(:,1)*rho2arr(:,1)*rho3arr(:,1) 422 else if (nspden==2) then 423 do ii=1,nrad 424 rho1u=rho1arr(ii,2) 425 rho1d=rho1arr(ii,1)-rho1arr(ii,2) 426 rho2u=rho2arr(ii,2) 427 rho2d=rho2arr(ii,1)-rho2arr(ii,2) 428 rho3u=rho3arr(ii,2) 429 rho3d=rho3arr(ii,1)-rho3arr(ii,2) 430 ff(ii)=& 431 ! uuu uud 432 & kxc(ii,ipts,1)*rho1u*rho2u*rho3u + kxc(ii,ipts,2)*rho1u*rho2u*rho3d + & 433 ! udu udd 434 & kxc(ii,ipts,2)*rho1u*rho2d*rho3u + kxc(ii,ipts,3)*rho1u*rho2d*rho3d + & 435 ! duu dud 436 & kxc(ii,ipts,2)*rho1d*rho2u*rho3u + kxc(ii,ipts,3)*rho1d*rho2u*rho3d + & 437 ! ddu ddd 438 & kxc(ii,ipts,3)*rho1d*rho2d*rho3u + kxc(ii,ipts,4)*rho1d*rho2d*rho3d 439 end do 440 end if 441 end if 442 443 ff(1:nrad)=ff(1:nrad)*pawrad%rad(1:nrad)**2 444 call simp_gen(d3exc1_int,ff,pawrad) 445 d3exc1_iat(1)=d3exc1_iat(1)+d3exc1_int*pawang%angwgth(ipts) 446 447 ! ----- End of the loop on npts (angular part) 448 end do 449 450 d3exc1_iat = d3exc1_iat*four_pi 451 452 ABI_FREE(ff) 453 ABI_FREE(rho1arr) 454 ABI_FREE(rho2arr) 455 ABI_FREE(rho3arr) 456 457 end subroutine paw_dfptnl_xc