added LDA-1/2 correction to LD1. Courtesy of Leonardo Matheus Marion Jorge

from University of San Paolo.


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

atomic/Makefile

@@ -71,6 +71,7 @@ pseudovloc.o \
 pseudo_q.o \
 read_pseudo_ncpp.o \
 read_pseudo_rrkj3.o \
+run_lda_half.o \
 run_pseudo.o \
 run_test.o \
 scf.o \

atomic/ld1.f90

@@ -60,6 +60,13 @@ program ld1
      call run_test ( )
      call ld1_writeout ( )
      !
+  elseif (iswitch.eq.4) then
+     !
+     ! LDA-1/2 correction to the input pseudopotential 
+     !
+     call run_lda_half ( )
+     call ld1_writeout ( )
+     !
   else
      call errore('ld1','iswitch not implemented',1)
   endif

atomic/ld1_readin.f90

@@ -39,7 +39,8 @@ subroutine ld1_readin
                          rmatch_augfun, which_augfun,         & !paw
                          rhos, bmat, lsmall, &              ! extra for paw2us
                          lgipaw_reconstruction, lsave_wfc, &
-                         relpert, noscf
+                         relpert, noscf, &
+                         rcutv ! LDA-1/2
 
   use funct, only : set_dft_from_name
   use radial_grids, only: do_mesh, check_mesh
@@ -105,7 +106,8 @@ subroutine ld1_readin
        ecutmin,       & ! for test with spherical Bessel functions:
        ecutmax,       & ! min and max energy cutoff for j_l(qr),
        decut,         & ! step: ecut = ecutmin, ecutmin+decut, ... , ecutmax
-       rm               ! radius of the box
+       rm,            & ! radius of the box
+       rcutv            ! CUT for LDA-1/2
 
   namelist /inputp/ &
        pseudotype,&! the pseudopotential type
@@ -237,7 +239,8 @@ subroutine ld1_readin
      if (zdum /= zed) call errore &
           ('ld1_readin','inconsistent Z/atom specification',nint(zdum))
   end if
-  if (iswitch < 1 .or. iswitch > 3) &
+! with LDA-1/2 now iswitch <=4
+  if (iswitch < 1 .or. iswitch > 4) &
        call errore('ld1_readin','wrong iswitch',1)
   if (eminld > emaxld) &
        call errore('ld1_readin','eminld or emaxld wrong',1)
@@ -322,7 +325,8 @@ subroutine ld1_readin
   ! generate the radial grid - note that if iswitch = 2 the radial grid
   ! is not generated but read from the pseudopotential file
   !
-  if (iswitch /= 2) then
+! also for LDA-1/2 radial grid has to be read
+  if (iswitch /= 2.or.iswitch/=4) then
      call do_mesh(rmax,zed,xmin,dx,0,grid)
      rhoc=0.0_dp
   endif
@@ -442,10 +446,21 @@ subroutine ld1_readin
   decut   = 5.0_dp
   rm      =30.0_dp
 
+!
+! default value for LDA-1/2
+!
+  rcutv = -1.0
+!
   if (ionode) read(5,test,err=300,iostat=ios)
 300  call mp_bcast(ios, ionode_id)
 
-  if (iswitch==2) call errore('ld1_readin','reading test',abs(ios))
+! LDA-1/2
+  if(iswitch==4) nconf = 2
+  if(iswitch==4.and.rcutv<0.0) call errore('ld1_readin','inconsistent rcutv',1)
+!
+
+! Added iswitch ==4 for LDA-1/2
+  if (iswitch==2.or.iswitch==4) call errore('ld1_readin','reading test',abs(ios))
   call bcast_test()
   call mp_bcast(configts, ionode_id)
   !
@@ -557,7 +572,7 @@ subroutine ld1_readin
   !
   !    PP testing: reading the pseudopotential
   !
-  if (iswitch ==2) then
+  if (iswitch ==2.or.iswitch==4) then
      lpaw=.false.
      !
      if (file_pseudo == ' ') &

atomic/ld1_writeout.f90

@@ -21,7 +21,8 @@ subroutine ld1_writeout
                      nconf , lpaw, rel, pawsetup, pseudotype, &
                      rhoc, vnl, phits, vpsloc, & 
                      elts, llts, octs, rcut, etots, nwfts, &
-                     lmax, lloc, zval, nlc, nnl, alps, alpc, alc, cc, nlcc
+                     lmax, lloc, zval, nlc, nnl, alps, alpc, alc, cc, nlcc, &
+                     iswitch
   use funct, only : get_dft_name
   use paw_type, only : deallocate_pseudo_paw
 
@@ -37,7 +38,7 @@ subroutine ld1_writeout
   
   if (file_pseudopw == ' ') return
 
-  if (nconf > 1) &
+  if (iswitch/=4.and.nconf > 1) &
        call errore('ld1_writeout','more than one test configuration',1)
 
   if ( (( rel == 2) .or. lpaw) &

atomic/ld1inc.f90

@@ -155,6 +155,10 @@ module ld1inc
   character(len=2) ::  &
        elts(nwfsx)           ! the label of the states
   !
+  ! for LDA-1/2
+  !
+  real(DP) :: rcutv ! CUT for LDA-1/2
+  !
   !    The control of the run
   !
   integer ::      &