./64_psp/psden.F90

ABINIT/psden

ABINIT/psden [ Functions ]

NAME

 psden

FUNCTION

 Calculate a pseudo-density from an original density on a radial grid (regular or logarithmic)

 Copyright (C) 1998-2017 ABINIT group (GJ,FJ,MT)
 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

  ilog=1 if grid is logarithmic, else 0
  mesh= dimension of nc
  nc(mesh)= density to be pseudized
  rc= cut-off radius
  rad(mesh) = radial mesh

OUTPUT

  ff(mesh)= pseudized density

  Optional:
    ff1(mesh)= 1st derivative of pseudo density (only r<rc modified)
    ff2(mesh)= 2nd derivative of pseudo density (only r<rc modified)

NOTES

    ff=exp(-(a+b.r^2+c.r^4))

PARENTS

      psp6cc

CHILDREN

      ctrap

SOURCE

 41 #if defined HAVE_CONFIG_H
 42 #include "config.h"
 43 #endif
 44 
 45 #include "abi_common.h"
 46 
 47 subroutine psden(ilog,ff,mesh,nc,rc,rad,ff1,ff2)
 48 
 49  use defs_basis
 50  use m_profiling_abi
 51  use m_errors
 52 
 53 !This section has been created automatically by the script Abilint (TD).
 54 !Do not modify the following lines by hand.
 55 #undef ABI_FUNC
 56 #define ABI_FUNC 'psden'
 57  use interfaces_32_util
 58 !End of the abilint section
 59 
 60  implicit none
 61 
 62 !Arguments ------------------------------------
 63 !scalars
 64  integer,intent(in) :: ilog,mesh
 65  real(dp),intent(in) :: rc
 66 !arrays
 67  real(dp),intent(in) :: nc(mesh),rad(mesh)
 68  real(dp),intent(out) :: ff(mesh)
 69  real(dp),intent(inout),optional :: ff1(mesh),ff2(mesh)
 70 
 71 !Local variables-------------------------------
 72 !scalars
 73  integer :: ii,nc1
 74  real(dp) :: aa,aa1,aa2,bb,cc,c1,c3,f0,f0p,norm1,norm2,rc1,step
 75 !arrays
 76  real(dp),allocatable :: fpir(:),gg(:)
 77 
 78 ! *************************************************************************
 79 
 80  rc1=rc/four
 81 
 82  ABI_ALLOCATE(fpir,(mesh))
 83  fpir(1:mesh)=four_pi*rad(1:mesh)**2
 84  if (ilog==1) fpir(1:mesh)=fpir(1:mesh)*rad(1:mesh)
 85 
 86  if (ilog==0) then
 87    step=rad(2)-rad(1)
 88    nc1=int(rc1/step)+1
 89    rc1=(nc1-1)*step
 90  else if (ilog==1) then
 91    step=log(rad(2)/rad(1))
 92    nc1=int(log(rc1/rad(1))/step)+1
 93    rc1=rad(nc1)
 94  end if
 95  ff(1:nc1)=nc(1:nc1)*fpir(1:nc1)
 96  call ctrap(nc1,ff(1:nc1),step,c3)
 97  if (ilog==1) c3=c3+half*ff(1)
 98  f0=nc(nc1);c1=-log(f0)
 99  f0p=half*(nc(nc1+1)-nc(nc1-1))/step
100 
101  ii=0;aa1=zero;norm1=c3+one
102  do while (norm1>c3.and.ii<100)
103    ii=ii+1;aa1=aa1+one
104    aa=c1-aa1*rc1**4+rc1*(f0p/f0+four*aa1*rc1**3)*half
105    bb=-half*(f0p/f0+four*aa1*rc1**3)/rc1
106    ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa1*rad(1:nc1)**4)
107    call ctrap(nc1,ff(1:nc1),step,norm1)
108    if (ilog==1) norm1=norm1+half*ff(1)
109  end do
110  if (ii==100) then 
111    MSG_ERROR('Big pb 1 in psden !')
112  end if
113 
114  ii=0;aa2=zero;norm2=c3-one
115  do while (norm2<c3.and.ii<100)
116    ii=ii+1;aa2=aa2-one
117    aa=c1-aa2*rc1**4+rc1*(f0p/f0+four*aa2*rc1**3)*half
118    bb=-half*(f0p/f0+four*aa2*rc1**3)/rc1
119    ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa2*rad(1:nc1)**4)
120    call ctrap(nc1,ff(1:nc1),step,norm2)
121    if (ilog==1) norm2=norm2+half*ff(1)
122  end do
123  if (ii==100) then 
124    MSG_ERROR('Big pb 2 in psden !')
125  end if
126 
127  do while (abs(norm2-c3)>tol10)
128 
129    cc=(aa1+aa2)*half
130    aa=c1-cc*rc1**4+rc1*(f0p/f0+four*cc*rc1**3)*half
131    bb=-half*(f0p/f0+four*cc*rc1**3)/rc1
132    ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-cc*rad(1:nc1)**4)
133    call ctrap (nc1,ff(1:nc1),step,norm2)
134    if (ilog==1) norm2=norm2+half*ff(1)
135    if ((norm1-c3)*(norm2-c3)>zero) then
136      aa1=cc
137      norm1=norm2
138    else
139      aa2=cc
140    end if
141 
142  end do ! while
143 
144  ff(1)=exp(-aa);if (ilog==1) ff(1)=ff(1)*exp(-bb*rad(1)**2-cc*rad(1)**4)
145  ff(2:nc1)=ff(2:nc1)/fpir(2:nc1)
146  if (nc1<mesh) ff(nc1+1:mesh)=nc(nc1+1:mesh)
147  if (present(ff1)) ff1(1:nc1)=-(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)*ff(1:nc1)
148  if (present(ff2)) ff2(1:nc1)=-(two*bb+12.0_dp*cc*rad(1:nc1)**2)*ff(1:nc1) &
149 & +(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)**2*ff(1:nc1)
150 
151  ABI_ALLOCATE(gg,(mesh))
152  gg(1:mesh)=fpir(1:mesh)*ff(1:mesh)
153  call ctrap(mesh,gg(1:mesh),step,norm1)
154  if (ilog==1) norm1=norm1+half*gg(1)
155  write(std_out,*) 'psden: tild_nc integral= ',norm1
156  ABI_DEALLOCATE(gg)
157 
158  ABI_DEALLOCATE(fpir)
159 
160 end subroutine psden

TABLE OF CONTENTS

ABINIT/psden [ Functions ]