./64_psp/vhtnzc.F90

ABINIT/vhtnzc

ABINIT/vhtnzc [ Functions ]

NAME

 vhtnzc

FUNCTION

 Compute VHartree(tild[n_Z+n_core]) from input ncore

 Copyright (C) 1998-2017 ABINIT group (GJ, FJ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  mesh=dimension of radial mesh
  nc= core density (to be pseudized)
  rad(mesh)=radial mesh
  rc=cut-off radius
  znucl=nuclear number of atom as specified in psp file

OUTPUT

  vhtnzc(mesh) = hartree potential induced by density tild[n_Z+n_core] (pseudo core density + nucleus)

PARENTS

      psp6cc

CHILDREN

      ctrap

SOURCE

 34 #if defined HAVE_CONFIG_H
 35 #include "config.h"
 36 #endif
 37 
 38 #include "abi_common.h"
 39 
 40 subroutine vhtnzc(nc,rc,vh_tnzc,mesh,rad,znucl)
 41 
 42  use defs_basis
 43  use m_profiling_abi
 44 
 45 !This section has been created automatically by the script Abilint (TD).
 46 !Do not modify the following lines by hand.
 47 #undef ABI_FUNC
 48 #define ABI_FUNC 'vhtnzc'
 49  use interfaces_32_util
 50 !End of the abilint section
 51 
 52  implicit none
 53 
 54 !Arguments ------------------------------------
 55 !scalars
 56  integer,intent(in) :: mesh
 57  real(dp),intent(in) :: znucl
 58  real(dp),intent(in) :: rc
 59 !arrays
 60  real(dp),intent(in) :: nc(mesh),rad(mesh)
 61  real(dp),intent(out) :: vh_tnzc(mesh)
 62 
 63 !Local variables-------------------------------
 64 !scalars
 65  integer :: ir,nc1
 66  real(dp) :: gnorm,rc1,step,yp1,yp2,yp3
 67 !arrays
 68  real(dp),allocatable :: den1(:),den2(:),den3(:),den4(:),nzc(:),rvhn(:),shapefunc(:)
 69 
 70 ! *************************************************************************
 71 
 72  rc1=rc/four
 73 
 74  step=log(rad(2)/rad(1))
 75  nc1=int(log(rc1/rad(1))/step)+1
 76  rc1=rad(nc1)
 77 
 78  ABI_ALLOCATE(shapefunc,(mesh))
 79  shapefunc(1)=one
 80  shapefunc(2:nc1)=(sin(pi*rad(2:nc1)/rc1)/(pi*rad(2:nc1)/rc1))**2
 81  if (nc1<mesh) shapefunc(nc1+1:mesh)=zero
 82 
 83  ABI_ALLOCATE(den1,(mesh))
 84  den1(1:mesh)=four_pi*shapefunc(1:mesh)*rad(1:mesh)**3
 85  call ctrap(mesh,den1,step,gnorm)
 86  gnorm =one/gnorm
 87  ABI_DEALLOCATE(den1)
 88 
 89  ABI_ALLOCATE(nzc,(mesh))
 90  nzc(1:mesh)=four*pi*nc(1:mesh)*rad(1:mesh)**2-four_pi*shapefunc(1:mesh)*rad(1:mesh)**2*znucl*gnorm
 91  ABI_DEALLOCATE(shapefunc)
 92 
 93  ABI_ALLOCATE(rvhn,(mesh))
 94  rvhn(1)=zero
 95 
 96  ABI_ALLOCATE(den1,(mesh))
 97  ABI_ALLOCATE(den2,(mesh))
 98  ABI_ALLOCATE(den3,(mesh))
 99  ABI_ALLOCATE(den4,(mesh))
100 
101  den1(1)=zero;den2(1)=zero
102  do ir=2,mesh
103    den1(ir)= rad(ir)*nzc(ir)
104    den2(ir)= den1(ir)/rad(ir)
105  end do
106 
107 !For first few points do stupid integral
108  den3(1)=zero;den4(1)=zero
109  do ir=2,mesh
110    call ctrap(ir,den1(1:ir),step,den3(ir))
111    call ctrap(ir,den2(1:ir),step,den4(ir))
112  end do
113 
114  do ir=1,mesh
115    rvhn(ir)=den3(ir)+rad(ir)*(den4(mesh)-den4(ir))
116  end do
117 
118  ABI_DEALLOCATE(den1)
119  ABI_DEALLOCATE(den2)
120  ABI_DEALLOCATE(den3)
121  ABI_DEALLOCATE(den4)
122 
123  vh_tnzc(2:mesh)=rvhn(2:mesh)/rad(2:mesh)
124  yp2=(vh_tnzc(3)-vh_tnzc(2))/(rad(3)-rad(2))
125  yp3=(vh_tnzc(4)-vh_tnzc(3))/(rad(4)-rad(3))
126  yp1=yp2+(yp2-yp3)*rad(2)/(rad(3)-rad(2))
127  vh_tnzc(1)=vh_tnzc(2)-(yp1+yp2)*rad(2)
128 
129  ABI_DEALLOCATE(nzc)
130  ABI_DEALLOCATE(rvhn)
131 
132 end subroutine vhtnzc

TABLE OF CONTENTS

ABINIT/vhtnzc [ Functions ]