!#define __DEBUG_V_H_vs_SPHEROPOLE ! ! Copyright (C) 2004 PWSCF group ! This file is distributed under the terms of the ! GNU General Public License. See the file `License' ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! MODULE atomic_paw ! !============================================================================ ! ! Module for Projector Augmented Wave (PAW) calculations assuming ! spherical symmetry. Kresse's notations are adopted. ! Contains the type describing a PAW dataset, the routine for ! generating it from the ld1 code, and the routines to compute ! the hamiltonian and energy. ! GGA and LSD are implemented, relativistic calculations are not. ! ! References: ! P.E.Blochl, Phys. Rev. B 50, 17953 (1994) ! G.Kresse, D.Joubert, Phys. Rev. B 59, 1758 (1999) ! ! WARNINGS: ! Still work in progress on many aspects. ! The PAW dataset is written in a temporary format which is not supposed ! to be used any longer. ! Consistency with the input of the ld1 code yet to be checked ! ! Written by Guido Fratesi, february 2005 ! Modified by Riccardo Mazzarello, july 2006 ! Other people involved: Lorenzo Paulatto and Stefano de Gironcoli ! !============================================================================ ! USE kinds, ONLY: dp USE ld1_parameters, ONLY: nwfsx USE parameters, ONLY: lmaxx USE constants, ONLY: pi, fpi, e2, eps8 USE radial_grids, ONLY: ndmx, radial_grid_type USE paw_type, ONLY: paw_t, nullify_pseudo_paw, allocate_pseudo_paw ! IMPLICIT NONE PRIVATE SAVE ! REAL(dp), PARAMETER :: ZERO=0._dp, ONE=1._dp, TWO=2._dp, HALF=0.5_dp ! ! TEMP, to be substituted by module constants ! REAL(dp), PARAMETER :: & ! PI=3.14159265358979323846_dp, FPI=4._dp*PI, E2=2._dp, EPS8=1.0e-8_dp ! !============================================================================ ! PUBLIC :: paw_t PUBLIC :: us2paw PUBLIC :: paw2us PUBLIC :: check_multipole PUBLIC :: new_paw_hamiltonian PUBLIC :: find_bes_qi PUBLIC :: compute_nonlocal_coeff_ion ! CONTAINS ! !============================================================================ ! PUBLIC ROUTINES !!! !============================================================================ ! ! Compute the values of the local pseudopotential and the NL coefficients ! Compute also the total energy ! SUBROUTINE new_paw_hamiltonian (veffps_, ddd_, etot_, & pawset_, nwfc_, l_, j_, nspin_, spin_, oc_, pswfc_, eig_, paw_energy,dddion_) IMPLICIT NONE REAL(dp), INTENT(OUT) :: veffps_(ndmx,2) REAL(dp), INTENT(OUT) :: ddd_(nwfsx,nwfsx,2) REAL(dp), INTENT(OUT) :: etot_ TYPE(paw_t), INTENT(IN) :: pawset_ INTEGER, INTENT(IN) :: nwfc_ INTEGER, INTENT(IN) :: l_(nwfsx) INTEGER, INTENT(IN) :: nspin_ INTEGER, INTENT(IN) :: spin_(nwfsx) REAL(dp), INTENT(IN) :: j_(nwfsx) REAL(dp), INTENT(IN) :: oc_(nwfsx) REAL(dp), INTENT(IN) :: pswfc_(ndmx,nwfsx) REAL(dp), INTENT(IN) :: eig_(nwfsx) REAL(dp), OPTIONAL :: dddion_(nwfsx,nwfsx) REAL(dp), INTENT(OUT), OPTIONAL :: paw_energy(5,3) ! REAL(dp) :: & ! one center: eps, e1, e1ps, & ! energies; veff1(ndmx,2), veff1ps(ndmx,2), & ! eff potentials; chargeps(ndmx,2), charge1(ndmx,2), charge1ps(ndmx,2), & ! charges. projsum(nwfsx,nwfsx,2), eigsum ! sum of projections, sum of eigenval. ! INTEGER :: ns, is, n REAL(dp) :: energy(5,3) ! ! Compute the valence charges CALL compute_charges(projsum, chargeps, charge1, charge1ps, & pawset_, nwfc_, l_, j_, nspin_, spin_, oc_, pswfc_, 1 ) ! ! Check for negative charge ! do is=1,nspin_ do n=2,pawset_%grid%mesh !if (chargeps(n,is)<-1.d-12) & ! call errore('new_paw_hamiltonian','negative rho',1) enddo enddo ! write(766,"(4f12.6)") (pawset_%grid%r(ns), chargeps(ns,1), charge1(ns,1), charge1ps(ns,1), ns=1,pawset_%grid%mesh) ! write(767,"(4f12.6)") (pawset_%grid%r(ns), chargeps(ns,2), charge1(ns,2), charge1ps(ns,2), ns=1,pawset_%grid%mesh) ! ! Compute the one-center energy and effective potential: ! E = Eh[n_v] + Exc[n_v+n_c] - Int[veff*n_v], ! veff = vh[n_v] + v_xc[n_v+n_c], ! where n_v can be n~+n^ or n1 or n1~+n^ and n_c can be nc or nc~ ! n~+n^ , nc~ CALL compute_onecenter_energy ( eps, veffps_, & pawset_, chargeps, pawset_%nlcc, pawset_%psccharge, nspin_,& pawset_%grid%mesh, pawset_%psloc, energy(1,1)) ! n1 , nc CALL compute_onecenter_energy ( e1, veff1, & pawset_, charge1, .TRUE., pawset_%aeccharge, nspin_,& pawset_%irc, pawset_%aeloc, energy(1,2)) ! n1~+n^ , nc~ CALL compute_onecenter_energy ( e1ps, veff1ps, & pawset_, charge1ps, pawset_%nlcc, pawset_%psccharge, nspin_,& pawset_%irc, pawset_%psloc, energy(1,3)) ! Add the local part DO is=1,nspin_ veffps_(1:pawset_%grid%mesh,is) = veffps_(1:pawset_%grid%mesh,is) + & pawset_%psloc(1:pawset_%grid%mesh) veff1 (1:pawset_%grid%mesh,is) = veff1 (1:pawset_%grid%mesh,is) + & pawset_%aeloc(1:pawset_%grid%mesh) veff1ps(1:pawset_%grid%mesh,is) = veff1ps(1:pawset_%grid%mesh,is) + & pawset_%psloc(1:pawset_%grid%mesh) END DO ! ! Compute the nonlocal D coefficients CALL compute_nonlocal_coeff (ddd_,pawset_,nspin_,veffps_,veff1,veff1ps) IF (PRESENT(dddion_)) THEN CALL compute_nonlocal_coeff_ion (dddion_,pawset_) END IF ! ! Compute total energy eigsum=ZERO DO ns=1,nwfc_ IF (oc_(ns)>ZERO) eigsum = eigsum + oc_(ns)*eig_(ns) END DO etot_ = eigsum + eps + e1 - e1ps if (PRESENT(paw_energy)) paw_energy=energy ! END SUBROUTINE new_paw_hamiltonian ! !============================================================================ ! ! Convert the USPP to a PAW dataset ! SUBROUTINE us2paw (pawset_, & zval, grid, rmatch_augfun, ikk, & nbeta, lls, jjs, ocs, enls, els, rcutus, psipaw, phis, betas, & qvan, kindiff, & nlcc, aerhoc, psrhoc, aevtot, psvtot, which_paw_augfun ) USE funct, ONLY : dft_name, get_iexch, get_icorr, get_igcx, get_igcc USE ld1inc, ONLY : zed, file_screen USE paw_type, ONLY : nullify_pseudo_paw, allocate_pseudo_paw USE io_global, ONLY : stdout, ionode, ionode_id USE radial_grids, ONLY : allocate_radial_grid USE mp, only : mp_bcast IMPLICIT NONE TYPE(paw_t), INTENT(OUT) :: pawset_ REAL(dp), INTENT(IN) :: zval type(radial_grid_type), INTENT(IN) :: grid REAL(dp), INTENT(IN) :: rmatch_augfun INTEGER, INTENT(IN) :: ikk(nwfsx) ! INTEGER, INTENT(IN) :: nbeta INTEGER, INTENT(IN) :: lls(nwfsx) REAL(dp), INTENT(IN) :: ocs(nwfsx) CHARACTER(LEN=2), INTENT(IN) :: els(nwfsx) REAL(dp), INTENT(IN) :: jjs(nwfsx) REAL(dp), INTENT(IN) :: rcutus(nwfsx) REAL(dp), INTENT(IN) :: enls(nwfsx) REAL(dp), INTENT(IN) :: psipaw(ndmx,nwfsx) REAL(dp), INTENT(IN) :: phis(ndmx,nwfsx) REAL(dp), INTENT(IN) :: betas(ndmx,nwfsx) ! REAL(dp), INTENT(IN) :: qvan(ndmx,nwfsx,nwfsx) REAL(dp), INTENT(IN) :: kindiff(nwfsx,nwfsx) LOGICAL, INTENT(IN) :: nlcc REAL(dp), INTENT(IN) :: aerhoc(ndmx) REAL(dp), INTENT(IN) :: psrhoc(ndmx) REAL(dp), INTENT(IN) :: aevtot(ndmx) REAL(dp), INTENT(IN) :: psvtot(ndmx) CHARACTER(20), INTENT(IN) :: which_paw_augfun ! REAL(DP), EXTERNAL :: int_0_inf_dr CHARACTER, EXTERNAL :: atom_name*2 REAL(dp) :: vps(ndmx,2), projsum(nwfsx,nwfsx,2), ddd(nwfsx,nwfsx,2), dddion(nwfsx,nwfsx) INTEGER :: irc, ns, ns1, n, leading_power, mesh, ios REAL(dp) :: aux(ndmx), aux2(ndmx,2), raux REAL(dp) :: aecharge(ndmx,2), pscharge(ndmx,2) REAL(dp) :: etot INTEGER :: nspin=1, spin(nwfsx)=1 ! PAW generat. from spin-less calculation CHARACTER(LEN=6) :: shortname ! ! variables for aug. functions generation ! REAL(DP) :: energy(5,3), max_aug_cutoff = -1._dp INTEGER :: nc, iok ! index Bessel function, ... INTEGER :: l1,l2,l3, lll, ircm, ir, ir0 REAL(dp) :: twosigma2, rm ! needed for "gaussian" augfun REAL(dp) :: zeta, resi,usigma=4._dp ! needed for "Bloechl" augfun REAL(DP),ALLOCATABLE :: gaussian(:) ! needed for gaussian and Bloechl REAL(DP),ALLOCATABLE :: aug_real(:,:,:) ! needed for PSQ augfun REAL(DP) :: qc(2), xc(2), b1(2), b2(2) ! needed for BESSEL augfun REAL(DP), ALLOCATABLE :: j1(:,:) ! needed for BESSEL augfun ! mesh = grid%mesh irc = maxval(ikk(1:nbeta))+8 CALL nullify_pseudo_paw(pawset_) CALL allocate_pseudo_paw(pawset_,ndmx,nwfsx,lmaxx) CALL allocate_radial_grid(pawset_%grid, mesh) pawset_%symbol = atom_name(nint(zed)) pawset_%zval = zval ! ! Copy the mesh pawset_%grid%mesh = grid%mesh pawset_%grid%r(:) = grid%r(:) pawset_%grid%r2(:) = grid%r2(:) pawset_%grid%rab(:)= grid%rab(:) pawset_%grid%sqr(:)= grid%sqr(:) pawset_%grid%xmin = grid%xmin pawset_%grid%rmax = grid%rmax pawset_%grid%zmesh = grid%zmesh pawset_%grid%dx = grid%dx ! pawset_%rmatch_augfun = rmatch_augfun if (rmatch_augfun <= 0.0_dp) pawset_%rmatch_augfun = grid%r(irc) pawset_%irc = irc pawset_%ikk(1:nbeta)=ikk(1:nbeta) ! ! Copy the wavefunctions pawset_%nwfc = nbeta pawset_%l(1:nbeta) = lls(1:nbeta) pawset_%lmax = MAXVAL(pawset_%l(1:nbeta)) pawset_%oc(1:nbeta) = ocs(1:nbeta) pawset_%jj(1:nbeta) = jjs(1:nbeta) pawset_%rcutus(1:nbeta) = rcutus(1:nbeta) pawset_%els(1:nbeta)= els(1:nbeta) pawset_%enl(1:nbeta)= enls(1:nbeta) pawset_%aewfc(1:mesh,1:nbeta) = psipaw(1:mesh,1:nbeta) pawset_%pswfc(1:mesh,1:nbeta) = phis (1:mesh,1:nbeta) pawset_%proj (1:mesh,1:nbeta) = betas (1:mesh,1:nbeta) ! pawset_%augfun = 0._dp ! ! ! Augmentation functions: ! The specific radial part is not important, as long as the ! multipole moments of the exact augmentation functions are ! reproduced. So we write on the file the exact augmentation ! functions and their multipole moments, keeping in mind that ! the PW program should not use the radial part as is but ! substitute it with some smoothened analogue. ! ! Sdg>> not quite sure this is correct because the shape of augfun, ! arbitrary as it may be, determines the shape of vloc since this ! is obtained by unscreening vscf with the GIVEN augfun... ! ! moreover in order to give the right electrostatics in the solid the ! augmentation functions should not overlap. ! ! Compute the exact augmentation functions and their moments ! write(stdout,'(//,4x,a)') 'multipoles (all-electron charge) - (pseudo charge)' write(stdout,'(7x,2a,":",2a,2x,6a)') ' ns',' l1','ns1',' l2',& ' l=0 ',' l=1 ',' l=2 ',' l=3 ',' l=4 ', ' l=5 ' pawset_%augfun(:,:,:,:) = 0.0_dp pawset_%augmom(:,:,:) = 0.0_dp DO ns=1,nbeta l1=pawset_%l(ns) DO ns1=1,ns l2=pawset_%l(ns1) do l3 = max(l1-l2,l2-l1), l1+l2, 2 pawset_%augfun(1:mesh,ns,ns1,l3) = & pawset_%aewfc(1:mesh,ns) * pawset_%aewfc(1:mesh,ns1) - & pawset_%pswfc(1:mesh,ns) * pawset_%pswfc(1:mesh,ns1) pawset_%augfun(1:mesh,ns1,ns,l3) = pawset_%augfun(1:mesh,ns,ns1,l3) aux(1:irc) = pawset_%augfun(1:irc,ns,ns1,l3) * pawset_%grid%r(1:irc)**l3 lll = l1 + l2 + 2 + l3 pawset_%augmom(ns,ns1,l3)=int_0_inf_dr(aux(1:irc),pawset_%grid,irc,lll) pawset_%augmom(ns1,ns,l3)=pawset_%augmom(ns,ns1,l3) end do write (stdout,'(7x,2i3,a,2i3,10f8.4)') ns,l1,":",ns1,l2,& (pawset_%augmom(ns,ns1,l3), l3=0,l1+l2) END DO END DO ! ! IF (which_paw_augfun/='AE') THEN ! The following lines enables to test the idependence on the ! specific shape of the radial part of the augmentation functions ! in a spherically averager system (as is the case in atoms) only ! the zero-th moment contribute to the scf charge write(stdout, '(5x,2a)') 'Required augmentation: ',TRIM(which_paw_augfun) ! 101 pawset_%augfun(:,:,:,:) = 0.0_dp DO ns=1,nbeta l1 = pawset_%l(ns) DO ns1=1,ns l2 = pawset_%l(ns1) ! SELECT CASE (which_paw_augfun) CASE ('PSQ') continue ! the work is done at the end CASE ('QVAN') IF(ns==1 .and. ns1==1) & CALL infomsg('us2paw', 'WARNING: QVAN augmentation function are for testing ONLY: '//& 'they will not work in pw!') ! Choose the shape of the augmentation functions: NC Q ... pawset_%augfun(1:mesh,ns,ns1,0) = qvan(1:mesh,ns,ns1) ! Renormalize the aug. functions so that their integral is correct leading_power = l1 + l2 + 2 raux=int_0_inf_dr(pawset_%augfun(1:irc,ns,ns1,0),pawset_%grid,irc,leading_power) IF (ABS(raux) < eps8) CALL errore & ('ld1_to_paw','norm of augm.func. is too small',ns*100+ns1) raux=pawset_%augmom(ns,ns1,0)/raux pawset_%augfun(1:mesh,ns,ns1,0)=raux*pawset_%augfun(1:mesh,ns,ns1,0) ! CASE ('BG') IF(ns==1 .and. ns1==1) & CALL infomsg('us2paw', 'WARNING: using Bloechl style augmentation functions '//& 'is not a good idea, as analytical overlap are not '//& 'implemented in pwscf: use BESSEL or GAUSS instead.') IF(.not. allocated(gaussian)) ALLOCATE(gaussian(mesh)) ! use Bloechl style augmentation functions, as linear combinations of ! gaussians functions (this is quite pointless if the the plane-wave ! code doesn't use double-augmentation with analytical gaussian overlaps) DO ir0 = 1,mesh IF(grid%r(ir0) > pawset_%rmatch_augfun) & exit ENDDO ! look for a sigma large enough to keep (almost) all the gaussian in the aug sphere zeta = (usigma/pawset_%rmatch_augfun)**2 DO ir = 1,mesh gaussian(ir) = exp( - zeta*(grid%r(ir))**2 ) * grid%r2(ir) ENDDO DO l3 = max (l1-l2,l2-l1), l1+l2 ! Functions has to be normalized later, so I can use a constant factor ! = rc**l3 to prevent very large numbers when integrating: aux(1:grid%mesh) = gaussian(1:grid%mesh) * grid%r(1:grid%mesh)**l3 ! Normalization to have unitary multipole l3 ! and check norm of augmentation functions. raux = int_0_inf_dr(aux,pawset_%grid,mesh,l3+2) IF (abs(raux) < eps8) CALL errore & ('ld1_to_paw','norm of augm. func. is too small',ns*100+ns1) ! Check if gaussians are localized enough into the augmentation sphere, ! otherwise try again with a larger sigma. resi = int_0_inf_dr(aux,pawset_%grid,ir0, l3+2) resi = (raux-resi)/raux IF (abs(resi) > 1.e-10_dp) THEN usigma = usigma + .0625_dp GOTO 101 ENDIF raux=pawset_%augmom(ns,ns1,l3)/raux pawset_%augfun(1:mesh,ns,ns1,l3) = raux*gaussian(1:mesh) pawset_%augfun(1:mesh,ns1,ns,l3) = raux*gaussian(1:mesh) ENDDO DEALLOCATE(gaussian) ! CASE ('GAUSS') ! use linear combinations of gaussians functions, not the Bloechl style ! but it looks a bit alike... (used for testing, now obsolete) IF(ns==1 .and. ns1==1) & CALL infomsg('us2paw', 'GAUSS augmentation functions are ususally '//& 'harder than BESSEL; use BESSEL instead unless'//& ' you have discontinuity in local potential') ALLOCATE(gaussian(mesh)) ! rm = pawset_%rmatch_augfun twosigma2 = TWO*(rm/3.0_dp)**2 do ir=1,mesh if (grid%r(ir) <= rm) then gaussian(ir) = ( EXP(-grid%r(ir)**2/twosigma2) + & EXP(-(grid%r(ir)-TWO*rm)**2/twosigma2 ) - & TWO*EXP(-rm**2/twosigma2) ) * grid%r2(ir) else gaussian(ir) = 0.0_dp endif end do DO l3 = max (l1-l2,l2-l1), l1+l2 ! aux(1:irc) = gaussian(1:irc) * pawset_%grid%r(1:irc)**l3 ! calculate the corresponding multipole raux=int_0_inf_dr(aux,pawset_%grid,irc,l3+2) IF (ABS(raux) < eps8) CALL errore & ('ld1_to_paw','norm of augm. func. is too small',ns*100+ns1) raux=pawset_%augmom(ns,ns1,l3)/raux pawset_%augfun(1:mesh,ns,ns1,l3) = raux*gaussian(1:mesh) pawset_%augfun(1:mesh,ns1,ns,l3) = raux*gaussian(1:mesh) ! END DO DEALLOCATE(gaussian) ! CASE ('BESSEL') ! Use Kresse-Joubert style augmentation functions ! Defined as linear combination of Bessel functions. ALLOCATE (j1 (pawset_%grid%mesh,2)) do ir=1,irc if (grid%r(ir) 0 qq(ns,ns1)=ZERO END IF END DO END DO ! vpsloc(1:mesh)=pawset_%psloc(1:mesh) bmat(1:nbeta,1:nbeta)=pawset_%dion(1:nbeta,1:nbeta) ! rhos(1:mesh)=pawset_%pscharge(1:mesh) ! betas(1:mesh,1:nbeta)=pawset_%proj(1:mesh,1:nbeta) pseudotype=3 ! CALL set_dft_from_name (pawset_%dft) END SUBROUTINE paw2us ! !============================================================================ ! ! ... ! SUBROUTINE check_multipole (pawset_) USE radial_grids, ONLY: hartree USE io_global, ONLY : stdout IMPLICIT NONE TYPE(paw_t), INTENT(IN) :: pawset_ INTEGER:: mesh REAL(dp) :: r(ndmx), r2(ndmx), sqr(ndmx), dx INTEGER :: nbeta INTEGER :: lls(nwfsx) INTEGER :: ir, ns1, ns2, l1, l2, l3, irc, ir0 REAL(dp) :: auxpot(ndmx,0:2*lmaxx+2), auxrho(ndmx) ! ! set a few internal variables write (stdout,*) "check_multipole : lmaxx =",lmaxx mesh=pawset_%grid%mesh r(1:mesh)=pawset_%grid%r(1:mesh) r2(1:mesh)=pawset_%grid%r2(1:mesh) sqr(1:mesh)=pawset_%grid%sqr(1:mesh) dx=pawset_%grid%dx irc = pawset_%irc ! nbeta=pawset_%nwfc lls(1:nbeta)=pawset_%l(1:nbeta) ! do ns1=1,nbeta l1 = lls(ns1) do ns2=1,nbeta l2 = lls(ns2) auxpot(:,:) = 0.0_dp do l3 = max(l1-l2,l2-l1), l1+l2 auxrho(1:mesh) = & pawset_%aewfc(1:mesh,ns1) * pawset_%aewfc(1:mesh,ns2) - & pawset_%pswfc(1:mesh,ns1) * pawset_%pswfc(1:mesh,ns2) - & pawset_%augfun(1:mesh,ns1,ns2,l3) call hartree(l3,l1+l2+2,mesh,pawset_%grid,auxrho,auxpot(1,l3)) end do write (stdout,'(a,2i3,a,2i3)') " MULTIPOLO DI ",ns1,l1,":",ns2, l2 do ir=1,irc if (r(ir) < 1.0_dp) ir0 = ir end do do ir=ir0,irc+30, 3 write (stdout,'(10f8.4)') r(ir),(auxpot(ir,l3), l3=0,l1+l2) end do end do end do return END SUBROUTINE check_multipole ! !============================================================================ ! PRIVATE ROUTINES !!! !============================================================================ ! SUBROUTINE compute_charges (projsum_, chargeps_, charge1_, charge1ps_, & pawset_, nwfc_, l_, j_, nspin_, spin_, oc_, pswfc_ , iflag, unit_) USE io_global, ONLY : ionode IMPLICIT NONE REAL(dp), INTENT(OUT) :: projsum_(nwfsx,nwfsx,2) REAL(dp), INTENT(OUT) :: chargeps_(ndmx,2) REAL(dp), INTENT(OUT) :: charge1_(ndmx,2) REAL(dp), INTENT(OUT) :: charge1ps_(ndmx,2) TYPE(paw_t), INTENT(IN) :: pawset_ INTEGER, INTENT(IN) :: nwfc_ INTEGER, INTENT(IN) :: l_(nwfsx) INTEGER, INTENT(IN) :: nspin_ INTEGER, INTENT(IN) :: spin_(nwfsx) REAL(dp), INTENT(IN) :: j_(nwfsx) REAL(dp), INTENT(IN) :: oc_(nwfsx) REAL(dp), INTENT(IN) :: pswfc_(ndmx,nwfsx) INTEGER, OPTIONAL :: unit_, iflag REAL(dp) :: augcharge(ndmx,2), chargetot INTEGER :: i, n, iflag0 iflag0=0 if (present(iflag)) iflag0=iflag CALL compute_sumwfc2(chargeps_,pawset_,nwfc_,pswfc_,oc_,spin_) CALL compute_projsum(projsum_,pawset_,nwfc_,l_,j_,spin_,pswfc_,oc_) ! WRITE (6200,'(20e20.10)') ((projsum_(ns,ns1,1),ns=1,ns1),ns1=1,pawset_%nwfc) CALL compute_onecenter_charge(charge1ps_,pawset_,projsum_,nspin_,"PS") CALL compute_onecenter_charge(charge1_ ,pawset_,projsum_,nspin_,"AE") ! add augmentation charges CALL compute_augcharge(augcharge,pawset_,projsum_,nspin_) if (present(unit_).and.ionode) then write(unit_,*) write(unit_,*) "#" do i=1,pawset_%grid%mesh write (unit_,'(4f12.8)') pawset_%grid%r(i), augcharge(i,1), chargeps_(i,1), charge1ps_(i,1) end do end if chargeps_ (1:pawset_%grid%mesh,1:nspin_) = chargeps_ (1:pawset_%grid%mesh,1:nspin_) & + augcharge(1:pawset_%grid%mesh,1:nspin_) charge1ps_(1:pawset_%grid%mesh,1:nspin_) = charge1ps_(1:pawset_%grid%mesh,1:nspin_) & + augcharge(1:pawset_%grid%mesh,1:nspin_) ! ! If there are unbounded scattering states in the pseudopotential generation, ! n1 and n~1 separately diverge. This makes the hartree and exchange and ! correlation energies separately diverging. The cancellation becomes ! very difficult numerically. Outside the sphere however the two charges ! should be equal and opposite and we set them to zero. ! Is this the right solution? ! do n=pawset_%irc+1,pawset_%grid%mesh chargetot=chargeps_(n,1) if (nspin_==2) chargetot=chargetot+chargeps_(n,2) if (chargetot<1.d-11.or.iflag0==1) then charge1_(n,1:nspin_)=0.0_DP charge1ps_(n,1:nspin_)=0.0_DP endif enddo ! do n=1,pawset_%grid%mesh ! write(6,'(4f20.15)') pawset_%grid%r(n), chargeps_(n,1), charge1_(n,1), ! charge1ps_(n,1) END SUBROUTINE compute_charges ! !============================================================================ ! ! Compute the one-center energy and effective potential: ! E = Eh[n_v] + Exc[n_v+n_c] - Int[veff*n_v], ! veff = vh[n_v] + v_xc[n_v+n_c], ! where n_v can be n~+n^ or n1 or n1~+n^ and n_c can be nc or n~c ! SUBROUTINE compute_onecenter_energy ( totenergy_, veff_, & pawset_, vcharge_, nlcc_, ccharge_, nspin_, iint, vloc, energies_ , unit_) USE funct, ONLY: dft_is_gradient, exc_t, vxc_t !igcx, igcc USE radial_grids, ONLY: hartree USE io_global, ONLY : stdout, ionode IMPLICIT NONE REAL(dp), INTENT(OUT) :: totenergy_ ! H+XC+DC REAL(dp), INTENT(OUT) :: veff_(ndmx,2) ! effective potential TYPE(paw_t), INTENT(IN) :: pawset_ ! the PAW dataset REAL(dp), INTENT(IN) :: vcharge_(ndmx,2) ! valence charge LOGICAL, INTENT(IN) :: nlcc_ ! non-linear core correction REAL(dp), INTENT(IN) :: ccharge_(ndmx) ! core charge INTEGER, INTENT(IN) :: nspin_ ! 1 for LDA, 2 for LSDA INTEGER, INTENT(IN) :: iint ! integrals up to iint REAL(dp), INTENT(IN), OPTIONAL :: vloc(ndmx) ! REAL(dp), INTENT(OUT), OPTIONAL :: energies_(5)! Etot,H,XC,DC terms INTEGER, OPTIONAL :: unit_ ! REAL(dp), PARAMETER :: rho_eq_0(ndmx) = ZERO ! ccharge=0 when nlcc=.f. ! REAL(dp) :: & eh, exc, edc, & ! hartree, xc and double counting energies eloc, & ! local energy rhovtot(ndmx), & ! total valence charge aux(ndmx), & ! auxiliary to compute integrals vh(ndmx), & ! hartree potential vxc(ndmx,2), & ! exchange-correlation potential (LDA+GGA) vgc(ndmx,2), & ! exchange-correlation potential (GGA only) egc(ndmx), & ! exchange correlation energy density (GGA only) rh(2), & ! valence charge at a given point without 4 pi r^2 rhc, & ! core charge at a given point without 4 pi r^2 vxcr(2) ! exchange-correlation potential at a given point ! INTEGER :: i, is INTEGER :: lsd REAL(DP), EXTERNAL :: int_0_inf_dr #if defined __DEBUG_V_H_vs_SPHEROPOLE REAL(DP) :: dummy_charge,aux1(ndmx),aux2(ndmx),res1,res2 #endif ! ! Set up total valence charge rhovtot(1:pawset_%grid%mesh) = vcharge_(1:pawset_%grid%mesh,1) IF (nspin_==2) rhovtot(1:pawset_%grid%mesh) = rhovtot(1:pawset_%grid%mesh) + & vcharge_(1:pawset_%grid%mesh,2) ! ! Hartree CALL hartree(0,2,pawset_%grid%mesh,pawset_%grid,rhovtot,vh) if (PRESENT(unit_).and.ionode) then write (unit_,*) " " write (unit_,*) "#" do i=1,pawset_%grid%mesh write (unit_,'(3f12.7)') pawset_%grid%r(i),rhovtot(i),vh(i) end do end if #if defined __DEBUG_V_H_vs_SPHEROPOLE dummy_charge=int_0_inf_dr(rhovtot,pawset_%grid,pawset_%grid%mesh,2) aux1(1:pawset_%grid%mesh) = FPI*pawset_%grid%r2(1:pawset_%grid%mesh)*vh(1:pawset_%grid%mesh) - & FPI*dummy_charge * pawset_%grid%r(1:pawset_%grid%mesh) aux2(1:pawset_%grid%mesh) = rhovtot(1:pawset_%grid%mesh)*pawset_%grid%r2(1:pawset_%grid%mesh) res1 = int_0_inf_dr(aux1,pawset_%grid,pawset_%grid%mesh,1) res2 = int_0_inf_dr(aux2,pawset_%grid,pawset_%grid%mesh,4) WRITE (stdout,'(4(A,1e15.7))') ' INT rho', dummy_charge,' INT V_H', & res1, ' INT r^2*rho', res2, ' ERR:', (1.d0- res1/ (-res2 * (2.d0*PI/3.d0))) #endif vh(1:pawset_%grid%mesh) = e2 * vh(1:pawset_%grid%mesh) aux(1:pawset_%grid%mesh) = vh(1:pawset_%grid%mesh) * rhovtot(1:pawset_%grid%mesh) eh = HALF * int_0_inf_dr(aux,pawset_%grid,iint,2) ! ! Exhange-Correlation rh=(/ZERO,ZERO/) rhc=ZERO lsd=nspin_-1 DO i=1,pawset_%grid%mesh DO is=1,nspin_ rh(is) = vcharge_(i,is)/pawset_%grid%r2(i)/FPI ENDDO IF (nlcc_) rhc = ccharge_(i)/pawset_%grid%r2(i)/FPI CALL vxc_t(rh,rhc,lsd,vxcr) vxc(i,1:nspin_)=vxcr(1:nspin_) IF (nlcc_) THEN aux(i)=exc_t(rh,rhc,lsd) * (rhovtot(i)+ccharge_(i)) ELSE aux(i)=exc_t(rh,rhc,lsd) * rhovtot(i) END IF END DO IF (dft_is_gradient()) THEN IF (nlcc_) THEN CALL vxcgc(ndmx,pawset_%grid%mesh,nspin_,pawset_%grid%r,& pawset_%grid%r2,vcharge_,ccharge_,vgc,egc,1) ELSE CALL vxcgc(ndmx,pawset_%grid%mesh,nspin_,pawset_%grid%r,& pawset_%grid%r2,vcharge_,rho_eq_0,vgc,egc,1) END IF vxc(1:pawset_%grid%mesh,1:nspin_) = vxc(1:pawset_%grid%mesh,1:nspin_) + & vgc(1:pawset_%grid%mesh,1:nspin_) aux(1:pawset_%grid%mesh) = aux(1:pawset_%grid%mesh) + & egc(1:pawset_%grid%mesh) * pawset_%grid%r2(1:pawset_%grid%mesh) * FPI END IF exc = int_0_inf_dr(aux,pawset_%grid,iint,2) ! ! Double counting edc=ZERO DO is=1,nspin_ veff_(1:pawset_%grid%mesh,is)=vxc(1:pawset_%grid%mesh,is)+vh(1:pawset_%grid%mesh) aux(1:pawset_%grid%mesh)=veff_(1:pawset_%grid%mesh,is)*vcharge_(1:pawset_%grid%mesh,is) edc=edc+int_0_inf_dr(aux,pawset_%grid,iint,2) END DO ! eloc=ZERO ! IF (present(vloc)) THEN DO is=1,nspin_ aux(1:pawset_%grid%mesh)=vloc(1:pawset_%grid%mesh) & *vcharge_(1:pawset_%grid%mesh,is) eloc=eloc+int_0_inf_dr(aux,pawset_%grid,iint,2) ENDDO ENDIF ! ! Total totenergy_ = eh + exc - edc ! IF (PRESENT(energies_)) THEN energies_(1)=totenergy_ energies_(2)=eh energies_(3)=exc energies_(4)=edc energies_(5)=eloc END IF ! END SUBROUTINE compute_onecenter_energy ! !============================================================================ ! ! Compute NL 'D' coefficients = D^ + D1 - D1~ ! SUBROUTINE compute_nonlocal_coeff(ddd_, pawset_, nspin_, veffps_, veff1_, veff1ps_) IMPLICIT NONE REAL(dp), INTENT(OUT) :: ddd_(nwfsx,nwfsx,2) TYPE(paw_t), INTENT(IN) :: pawset_ INTEGER, INTENT(IN) :: nspin_ REAL(dp), INTENT(IN) :: veffps_(ndmx,2) REAL(dp), INTENT(IN) :: veff1_(ndmx,2) REAL(dp), INTENT(IN) :: veff1ps_(ndmx,2) INTEGER :: is, ns, ns1 REAL(dp) :: aux(ndmx), dd REAL(DP), EXTERNAL :: int_0_inf_dr ! REAL(dp):: dddd(nwfsx,nwfsx,3) = 0.d0 ! ! D^ = Int Q*v~ ! D1-D1~ = kindiff + Int[ae*v1*ae - ps*v1~*ps - Q*v1~] ddd_(:,:,:)=ZERO DO is=1,nspin_ DO ns=1,pawset_%nwfc DO ns1=1,ns IF (pawset_%l(ns)==pawset_%l(ns1).and.& ABS(pawset_%jj(ns)-pawset_%jj(ns1))<1.d-8) THEN ! Int[Q*v~] aux(1:pawset_%grid%mesh) = & pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) * & veffps_(1:pawset_%grid%mesh,is) dd = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) ! dddd(ns,ns1,1) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) ! Int[ae*v1*ae] aux(1:pawset_%grid%mesh) = & pawset_%aewfc(1:pawset_%grid%mesh,ns ) * & pawset_%aewfc(1:pawset_%grid%mesh,ns1) * & veff1_(1:pawset_%grid%mesh,is) dd = dd + & int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) ! dddd(ns,ns1,2) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) ! Int[ps*v1~*ps + aufun*v1~] aux(1:pawset_%grid%mesh) = & ( pawset_%pswfc(1:pawset_%grid%mesh,ns ) * & pawset_%pswfc(1:pawset_%grid%mesh,ns1) + & pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) ) * & veff1ps_(1:pawset_%grid%mesh,is) dd = dd - & int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) ! dddd(ns,ns1,3) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) ! collect ddd_(ns,ns1,is) = pawset_%kdiff(ns,ns1) + dd ddd_(ns1,ns,is) = ddd_(ns,ns1,is) END IF END DO END DO END DO ! write(*,*) 'deeq', pawset_%irc ! write(*,'(4f13.8)') dddd(1:4,1:4,1) ! write(*,*) 'ddd ae' ! write(*,'(4f13.8)') dddd(1:4,1:4,2) ! write(*,*) 'ddd ps' ! write(*,'(4f13.8)') dddd(1:4,1:4,3) END SUBROUTINE compute_nonlocal_coeff ! !============================================================================ ! ! 'D_ion' coefficients = D1 - D1~ ! SUBROUTINE compute_nonlocal_coeff_ion(ddd_, pawset_) IMPLICIT NONE REAL(dp), INTENT(OUT) :: ddd_(nwfsx,nwfsx) TYPE(paw_t), INTENT(IN) :: pawset_ INTEGER :: ns, ns1 REAL(dp) :: aux(ndmx), dd REAL(DP), EXTERNAL :: int_0_inf_dr ! ! D^ = Int Q*v~ ! D1-D1~ = kindiff + Int[ae*v1*ae - ps*v1~*ps - Q*v1~] ! write(666,"(4f12.6)") (pawset_%grid%r(ns), veffps_(ns,1), veff1_(ns,1), veff1ps_(ns,1), ns=1,pawset_%grid%mesh) ! write(667,"(4f12.6)") (pawset_%grid%r(ns), veffps_(ns,2), veff1_(ns,2), veff1ps_(ns,2), ns=1,pawset_%grid%mesh) ddd_(:,:)=ZERO DO ns=1,pawset_%nwfc DO ns1=1,ns IF (pawset_%l(ns)==pawset_%l(ns1).and. & ABS(pawset_%jj(ns)-pawset_%jj(ns1))<1.d-8 ) THEN ! Int[ae*v1*ae] aux(1:pawset_%grid%mesh) = & pawset_%aewfc(1:pawset_%grid%mesh,ns ) * & pawset_%aewfc(1:pawset_%grid%mesh,ns1) * & pawset_%aeloc(1:pawset_%grid%mesh) dd = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) ! Int[ps*v1~*ps + Q*v1~] aux(1:pawset_%grid%mesh) = & ( pawset_%pswfc(1:pawset_%grid%mesh,ns ) * & pawset_%pswfc(1:pawset_%grid%mesh,ns1) + & pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) ) * & pawset_%psloc(1:pawset_%grid%mesh) dd = dd - & int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) ! ddd_(ns,ns1) = pawset_%kdiff(ns,ns1) + dd ddd_(ns1,ns)=ddd_(ns,ns1) END IF END DO END DO END SUBROUTINE compute_nonlocal_coeff_ion ! !============================================================================ ! ! Write PAW dataset wfc and potentials on files ! SUBROUTINE human_write_paw(pawset_) USE io_global, ONLY : ionode IMPLICIT NONE TYPE(paw_t), INTENT(In) :: pawset_ INTEGER :: n,ns IF (.not.ionode) return DO ns=1,pawset_%nwfc WRITE (5000+ns,'(A)') "# r AEwfc PSwfc projector" DO n=1,pawset_%grid%mesh WRITE (5000+ns,'(4f20.12)') pawset_%grid%r(n), pawset_%aewfc(n,ns), pawset_%pswfc(n,ns),pawset_%proj(n,ns) END DO END DO ! WRITE (6000,'(A)') "# r AEpot PSpot" DO n=1,pawset_%grid%mesh WRITE (6000,'(3f20.8)') pawset_%grid%r(n), pawset_%aeloc(n), pawset_%psloc(n) END DO END SUBROUTINE human_write_paw ! !============================================================================ ! LOWER-LEVEL ROUTINES !!! !============================================================================ ! ! Compute Sum_i oc_i * wfc_i^2 ! SUBROUTINE compute_sumwfc2(charge_, pawset_, nwfc_, wfc_, oc_, spin_) IMPLICIT NONE REAL(dp), INTENT(OUT) :: charge_(ndmx,2) TYPE(paw_t), INTENT(IN) :: pawset_ INTEGER, INTENT(IN) :: nwfc_ REAL(dp), INTENT(IN) :: wfc_(ndmx,nwfsx) REAL(dp), INTENT(IN) :: oc_(nwfsx) INTEGER, INTENT(IN) :: spin_(nwfsx) INTEGER :: nf charge_(1:pawset_%grid%mesh,:)=ZERO DO nf=1,nwfc_ IF (oc_(nf)>ZERO) charge_(1:pawset_%grid%mesh,spin_(nf)) = & charge_(1:pawset_%grid%mesh,spin_(nf)) + oc_(nf) * & wfc_(1:pawset_%grid%mesh,nf) * wfc_(1:pawset_%grid%mesh,nf) END DO END SUBROUTINE compute_sumwfc2 ! !============================================================================ ! ! Compute Sum_n oc_n ! SUBROUTINE compute_projsum (projsum_, pawset_, nwfc_, l_, j_, spin_, pswfc_, oc_) REAL(dp), INTENT(OUT) :: projsum_(nwfsx,nwfsx,2) TYPE(paw_t), INTENT(IN) :: pawset_ INTEGER, INTENT(IN) :: nwfc_ INTEGER, INTENT(IN) :: l_(nwfsx) INTEGER, INTENT(IN) :: spin_(nwfsx) REAL(dp), INTENT(IN) :: j_(nwfsx) REAL(dp), INTENT(IN) :: pswfc_(ndmx,nwfsx) REAL(dp), INTENT(IN) :: oc_(nwfsx) REAL(dp) :: proj_dot_wfc(nwfsx,nwfsx), aux(ndmx) INTEGER :: ns, ns1, nf, nr, is, nst REAL(DP), EXTERNAL :: int_0_inf_dr ! Compute DO ns=1,pawset_%nwfc DO nf=1,nwfc_ IF (pawset_%l(ns)==l_(nf).AND.pawset_%jj(ns)==j_(nf)) THEN DO nr=1,pawset_%grid%mesh aux(nr)=pawset_%proj(nr,ns)*pswfc_(nr,nf) END DO nst=(l_(nf)+1)*2 proj_dot_wfc(ns,nf)=int_0_inf_dr(aux,pawset_%grid,pawset_%ikk(ns),nst) ELSE proj_dot_wfc(ns,nf)=ZERO END IF END DO END DO ! Do the sum over the wavefunctions projsum_(:,:,:)=ZERO DO ns=1,pawset_%nwfc DO ns1=1,ns DO nf=1,nwfc_ is=spin_(nf) IF (oc_(nf)>ZERO) THEN projsum_(ns,ns1,is) = projsum_(ns,ns1,is) + oc_(nf) * & proj_dot_wfc(ns,nf) * proj_dot_wfc(ns1,nf) END IF END DO projsum_(ns1,ns,:)=projsum_(ns,ns1,:) END DO END DO END SUBROUTINE compute_projsum ! !============================================================================ ! ! ! Compute augmentation charge as Sum_ij W_ij * Q_ij ! SUBROUTINE compute_augcharge(augcharge_, pawset_, projsum_, nspin_) IMPLICIT NONE REAL(dp), INTENT(OUT) :: augcharge_(ndmx,2) TYPE(paw_t), INTENT(IN) :: pawset_ REAL(dp), INTENT(IN) :: projsum_(nwfsx,nwfsx,2) INTEGER, INTENT(IN) :: nspin_ INTEGER :: ns, ns1, is REAL(dp) :: factor augcharge_=ZERO DO is=1,nspin_ DO ns=1,pawset_%nwfc DO ns1=1,ns ! multiply by TWO off-diagonal terms factor=TWO IF (ns1==ns) factor=ONE ! ! NB: in a spherically averaged system only the l=0 component ! of the augmentation functions is present ! augcharge_(1:pawset_%grid%mesh,is) = augcharge_(1:pawset_%grid%mesh,is) + & factor * projsum_(ns,ns1,is) * & pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) END DO END DO END DO END SUBROUTINE compute_augcharge ! !============================================================================ ! ! Compute n1 and n1~, as Sum_ij W_ij wfc_i(r)*wfc_j(r) ! SUBROUTINE compute_onecenter_charge(charge1_, pawset_, projsum_, nspin_, which_wfc) IMPLICIT NONE REAL(dp), INTENT(OUT) :: charge1_(ndmx,2) TYPE(paw_t), INTENT(IN) :: pawset_ REAL(dp), INTENT(IN) :: projsum_(nwfsx,nwfsx,2) INTEGER, INTENT(IN) :: nspin_ CHARACTER(LEN=2),INTENT(IN) :: which_wfc INTEGER :: ns, ns1, is REAL(dp) :: factor charge1_=ZERO DO is=1,nspin_ DO ns=1,pawset_%nwfc DO ns1=1,ns ! multiply times TWO off-diagonal terms IF (ns1==ns) THEN factor=ONE ELSE factor=TWO END IF SELECT CASE (which_wfc) CASE ("AE") charge1_(1:pawset_%grid%mesh,is) = charge1_(1:pawset_%grid%mesh,is) + factor * & projsum_(ns,ns1,is) * pawset_%aewfc(1:pawset_%grid%mesh,ns) * & pawset_%aewfc(1:pawset_%grid%mesh,ns1) CASE ("PS") charge1_(1:pawset_%grid%mesh,is) = charge1_(1:pawset_%grid%mesh,is) + factor * & projsum_(ns,ns1,is) * pawset_%pswfc(1:pawset_%grid%mesh,ns) * & pawset_%pswfc(1:pawset_%grid%mesh,ns1) CASE DEFAULT call errore ('compute_onecenter_charge','specify AE or PS wavefunctions',1) END SELECT END DO END DO END DO END SUBROUTINE compute_onecenter_charge ! !-------------------------------------------------------------------------- SUBROUTINE find_bes_qi(qc,rmatch,lam,ncn,iok) !-------------------------------------------------------------------------- ! ! This routine finds two values of q such that the ! functions f_l have a derivative equal to 0 at rmatch ! IMPLICIT NONE INTEGER, INTENT(IN) :: & lam, & ! input: the angular momentum ncn ! input: the number of qi to compute INTEGER, INTENT(INOUT) :: & iok ! output: if 0 the calculation in this routine is ok REAL (dp), INTENT(OUT) :: & qc(ncn) ! output: the values of qi REAL (dp), INTENT(IN) :: rmatch REAL (dp) :: & zeroderjl (2,7) ! first two zeros of the first derivative of ! spherical Bessel function j_l for l = 0,...,6 INTEGER :: & nc ! counter on the q found data zeroderjl / 0.0_dp, 4.4934094579614_dp, & 2.0815759780862_dp, 5.9403699890844_dp, & 3.3420936578747_dp, 7.2899322987026_dp, & 4.5140996477983_dp, 8.5837549433127_dp, & 5.6467036213923_dp, 9.8404460168549_dp, & 6.7564563311363_dp, 11.0702068269176_dp, & 7.8510776799611_dp, 12.2793339053177_dp / iok=0 IF (ncn.gt.2) & CALL errore('find_aug_qi','ncn is too large',1) IF (lam.gt.6) & CALL errore('find_aug_qi','l not programmed',1) ! ! fix deltaq and the maximum step number ! DO nc = 1, ncn ! qc(nc) = zeroderjl (nc, lam + 1) / rmatch ! ENDDO RETURN END SUBROUTINE find_bes_qi ! !============================================================================ ! END MODULE atomic_paw