Added the possibility to pseudize the core charge with two bessel function.

(Experimental)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@4197 c92efa57-630b-4861-b058-cf58834340f0

atomic/ld1_readin.f90

@@ -33,7 +33,7 @@ subroutine ld1_readin
                          file_core, file_beta, file_chi, author, &
                          nld, rlderiv, eminld, emaxld, deld, &
                          ecutmin, ecutmax, decut, rytoev_fact, verbosity, &
-                         frozen_core, lsdts
+                         frozen_core, lsdts, new_core_ps
 
   use funct, only : set_dft_from_name
   use radial_grids, only: do_mesh, check_mesh
@@ -102,6 +102,7 @@ subroutine ld1_readin
        zval,  &    ! the pseudo valence
        lloc,  &    ! l component considered as local 
        nlcc,  &    ! if true nlcc is set
+       new_core_ps, & ! if true pseudize the core charge with bessel functions
        rcore, &    ! the core radius for nlcc
        rcloc, &    ! the local cut-off for pseudo
        lpaw,  &    ! if true create a PAW dataset
@@ -278,6 +279,7 @@ subroutine ld1_readin
      lloc=-1
      rcloc=-1_dp
      nlcc=.false.
+     new_core_ps=.false.
      rcore=0.0_dp
      rho0=0.0_dp
      tm  = .false.
@@ -587,6 +589,7 @@ implicit none
   call mp_bcast( zval,  ionode_id )
   call mp_bcast( lloc,  ionode_id )
   call mp_bcast( nlcc,  ionode_id )
+  call mp_bcast( new_core_ps,  ionode_id )
   call mp_bcast( rcore, ionode_id )
   call mp_bcast( rcloc, ionode_id )
   call mp_bcast( lpaw,  ionode_id )

atomic/ld1inc.f90

@@ -226,6 +226,7 @@ module ld1inc
        rhoc(ndmx)      ! the core charge
 
   logical :: &
+       new_core_ps, & ! if true pseudize the core charge with bessel functions
        nlcc    ! if true nlcc pseudopotential
   !
   !  the potential for the scf

atomic/set_rho_core.f90

@@ -18,12 +18,12 @@ subroutine set_rho_core
   use io_global, only : stdout, ionode, ionode_id
   use mp,        only : mp_bcast
   use ld1inc, only : nlcc, lpaw, grid, rhoc, aeccharge, psccharge, rcore, &
-                     nwf, oc, rel, core_state, psi, file_core
+                     nwf, oc, rel, core_state, psi, file_core, new_core_ps
   implicit none
 
   real(DP) :: drho, const, br1, br2, &
-       eps1, eps2, br12, a, b, eps12, totrho
-  real(DP), allocatable:: rhov(:), rhoco(:)
+       eps1, eps2, br12, xc(8), a, b, eps12, totrho
+  real(DP), allocatable:: rhov(:)
   real(DP), external :: int_0_inf_dr
   integer :: i, ik, n, ns, ios
 
@@ -32,7 +32,7 @@ subroutine set_rho_core
   else
      if (lpaw) write(stdout,'(/,5x,'' Computing core charge for PAW: '')')
   end if
-  allocate (rhov(grid%mesh), rhoco(grid%mesh))
+  allocate (rhov(grid%mesh))
   !
   !      calculates core charge density
   !
@@ -66,8 +66,7 @@ subroutine set_rho_core
   endif
 
 !  write(stdout,'("Integrated core charge",f15.10)') totrho
-  rhoco(:) = rhoc(1:grid%mesh)
-  if (lpaw) aeccharge(1:grid%mesh) = rhoc(1:grid%mesh)
+  aeccharge(1:grid%mesh) = rhoc(1:grid%mesh)
   !
   if (rcore > 0.0_dp) then
      !      rcore read on input
@@ -91,39 +90,47 @@ subroutine set_rho_core
   !      rhoc(r) = core charge        for r > rcore
   !      rhoc(r) = r^2 a sin(b r)/r   for r < rcore
   !
-  if (drho > 0.0_dp) then
-     call infomsg('set_rho_core','d rho/ d r > 0')
-     return
-  endif
-  const= grid%r(ik)*drho / ( rhoc(ik)/grid%r2(ik) ) + 1.0_dp
-  if (const > 0.0_dp) then
-     br1 = 0.00001_dp
-     br2 = pi/2.0_dp-0.00001_dp
+  if (new_core_ps) then
+     call compute_phius(1,ik,aeccharge,rhoc,xc,0,'  ')
   else
-     br1 = pi/2.0_dp+0.00001_dp
-     br2 = pi
-  end if
-  do n=1, 15
-     eps1 = br1 - const*tan(br1)
-     eps2 = br2 - const*tan(br2)
-     br12 = (br1+br2)/2.0_dp
-     eps12 = br12 - const*tan(br12)
-     if(eps1*eps12 < 0.0_dp) then
-        br2 = br12
-     else if(eps12*eps2 < 0.0_dp) then
-        br1 = br12
+     if (drho > 0.0_dp) then
+        call infomsg('set_rho_core','d rho/ d r > 0')
+        return
+     endif
+     const= grid%r(ik)*drho / ( rhoc(ik)/grid%r2(ik) ) + 1.0_dp
+     if (const > 0.0_dp) then
+        br1 = 0.00001_dp
+        br2 = pi/2.0_dp-0.00001_dp
      else
-        call errore('set_rho_core','error in bisection',n)
+        br1 = pi/2.0_dp+0.00001_dp
+        br2 = pi
      end if
-  end do
-  b = br12/grid%r(ik)
-  a = ( rhoc(ik)/grid%r2(ik) ) * grid%r(ik)/sin(br12)
-  do n=1,ik
-     rhoc(n) = a*sin(b*grid%r(n))/grid%r(n) * grid%r2(n)
-  end do
+     do n=1, 15
+        eps1 = br1 - const*tan(br1)
+        eps2 = br2 - const*tan(br2)
+        br12 = (br1+br2)/2.0_dp
+        eps12 = br12 - const*tan(br12)
+        if(eps1*eps12 < 0.0_dp) then
+           br2 = br12
+        else if(eps12*eps2 < 0.0_dp) then
+           br1 = br12
+        else
+           call errore('set_rho_core','error in bisection',n)
+        end if
+     end do
+     b = br12/grid%r(ik)
+     a = ( rhoc(ik)/grid%r2(ik) ) * grid%r(ik)/sin(br12)
+     do n=1,ik
+        rhoc(n) = a*sin(b*grid%r(n))/grid%r(n) * grid%r2(n)
+     end do
+  endif
   if (lpaw) psccharge(1:grid%mesh) = rhoc(1:grid%mesh)
   write(stdout,'(/,5x,''  r > '',f4.2,'' : true rho core'')') grid%r(ik)
-  write(stdout,110) grid%r(ik), a, b
+  if (new_core_ps) then
+     write(stdout,'(5x,"Core charge pseudized with two Bessel functions")')
+  else
+     write(stdout,110) grid%r(ik), a, b
+  endif
 110 format (5x, '  r < ',f4.2,' : rho core = a sin(br)/r', &
        '    a=',f7.2,'  b=',f7.2/)
 1100 continue
@@ -136,7 +143,7 @@ subroutine set_rho_core
      if (ionode) then
         if (totrho>1.d-6) then
            do n=1,grid%mesh
-              write(26,'(4f20.10)') grid%r(n),rhoc(n),rhov(n),rhoco(n)
+              write(26,'(4f20.10)') grid%r(n),rhoc(n),rhov(n),aeccharge(n)
            enddo
         else
            do n=1,grid%mesh
@@ -149,6 +156,6 @@ subroutine set_rho_core
   totrho = int_0_inf_dr(rhoc,grid,grid%mesh,2)
   write(stdout,'(6x,''Integrated core pseudo-charge : '',f6.2)')  totrho
   if (.not.nlcc) rhoc(1:grid%mesh) = 0.0_dp
-  deallocate (rhoco, rhov)
+  deallocate (rhov)
   return
 end subroutine set_rho_core

atomic_doc/INPUT_LD1

@@ -147,6 +147,9 @@ Needed for pseudo-potential (PP) test. optional for PP generation:
               corresponding value of rcut is used. 
       nlcc =.true. : produce a PP with the nonlinear core    (default .false.) 
              correction of Froyen, Cohen, and Louie
+      new_core_ps                                            (default .false.)
+             if .true. pseudizes the core charge with bessel functions. 
+             requires nlcc=.true.
       rcore= matching radius (a.u.) for the smoothing of the core charge
              If not specified, the matching radius is determined 
              by the condition  rho_core(rcore) = 2*rho_valence(rcore)