TABLE OF CONTENTS

ABINIT/fm12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

 fm12a1

FUNCTION

INPUTS

OUTPUT

SOURCE

1269  function fm12a1 (x)
1270 
1271 !Arguments -------------------------------
1272  real(dp),intent(in) :: x
1273  real(dp) :: fm12a1
1274 
1275  real(dp) :: y
1276 
1277 !
1278 !**********************************************************************
1279 !*                                                                    *
1280 !*               Integrale de Fermi d'ordre -1/2                      *
1281 !*    Fm12(x) = somme de 0 a l'infini de (dt*t**-1/2)/(1+exp(t-x))    *
1282 !*                                                                    *
1283 !**********************************************************************
1284 !
1285 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1286 !Erreur relative maximum annoncee 4.75 e-5
1287 !
1288  if (x.lt.2._dp) then
1289    y=exp(x)
1290    fm12a1=y*(23.1456_dp+y*(13.7820_dp+y))&
1291 &   /(13.0586_dp+y*(17.0048_dp+y*(5.07527_dp+y*0.236620_dp)))
1292  else
1293    y=1./(x*x)
1294    fm12a1=sqrt(x)*(0.0153602_dp+y*(0.146815_dp+y))&
1295 &   /(0.00768015_dp+y*(0.0763700_dp+y*0.570485_dp))
1296  end if
1297 !
1298 !**********************************************************************
1299  end function fm12a1

ABINIT/fm12a1t [ Functions ]

[ Top ] [ Functions ]

NAME

 fm12a1t

FUNCTION

INPUTS

OUTPUT

SOURCE

1314  subroutine fm12a1t (cktf,rtnewt,tphysel,vtrial,rhor_middx,rhor_mid,nfft)
1315 
1316  integer,intent(in) :: nfft
1317  real(dp),intent(in) :: tphysel,rtnewt,cktf
1318  real(dp),intent(in) :: vtrial(nfft)
1319  real(dp),intent(out) :: rhor_middx(nfft),rhor_mid(nfft)
1320 
1321  !intrinsic exp,sqrt
1322  integer :: ifft
1323  real(dp) :: x,y,sqrtx
1324 
1325 !
1326 !**********************************************************************
1327 !*                                                                    *
1328 !*               Integrale de Fermi d'ordre -1/2                      *
1329 !*    Fm12(x) = somme de 0 a l'infini de (dt*t**-1/2)/(1+exp(t-x))    *
1330 !*                      ....                                              *
1331 !**********************************************************************
1332 !
1333 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1334 !Erreur relative maximum annoncee 4.75 e-5
1335 !
1336  do ifft=1,nfft
1337    x=(rtnewt-vtrial(ifft))/tphysel
1338    if (x.lt.2._dp) then
1339      y=exp(x)
1340      rhor_middx(ifft)=cktf*y*(23.1456e0_dp+y*(13.7820e0_dp+y))&
1341 &     /(13.0586e0_dp+y*(17.0048e0_dp+y*(5.07527e0_dp+y*0.236620e0_dp)))
1342      rhor_mid(ifft)=cktf*y*(21.8168_dp+y*(13.1693_dp+y))&
1343 &     /(24.6180+y*(23.5546_dp+y*(4.76290_dp+y*0.134481_dp)))
1344    else
1345      y=1._dp/(x*x)
1346      sqrtx=sqrt(x)
1347      rhor_middx(ifft)=cktf*sqrtx*(0.0153602e0_dp+y*(0.146815e0_dp+y))&
1348 &     /(0.00768015e0_dp+y*(0.0763700e0_dp+y*0.570485e0_dp))
1349      rhor_mid(ifft)=cktf*x*sqrtx*(0.0473011_dp+y*(0.548433_dp+y))&
1350 &     /(0.0709478_dp+y*(0.737041_dp+y*0.382065_dp))
1351    end if
1352  end do
1353 !
1354 !**********************************************************************
1355  end subroutine fm12a1t

ABINIT/fp12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

 fp12a1

FUNCTION

INPUTS

OUTPUT

SOURCE

1131  function fp12a1 (x)
1132 
1133 ! Arguments -------------------------------
1134  real(dp),intent(in) :: x
1135  real(dp) :: fp12a1
1136 
1137  real(dp) :: y
1138 
1139 !**********************************************************************
1140 !*                                                                    *
1141 !*               Integrale de Fermi d'ordre 1/2                       *
1142 !*    Fp12(x) = somme de 0 a l'infini de (dt*t**1/2)/(1+exp(t-x))     *
1143 !*                                                                    *
1144 !**********************************************************************
1145 !
1146 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1147 !Erreur relative maximum annoncee 5.54 e-5
1148 !Erreur relative maximum constatee : -5.53e-5 pour eta = 2
1149 !
1150  if (x.lt.2._dp) then
1151    y=exp(x)
1152    fp12a1=y*(21.8168_dp+y*(13.1693_dp+y))&
1153 &   /(24.6180_dp+y*(23.5546_dp+y*(4.76290_dp+y*0.134481_dp)))
1154  else
1155    y=one/(x*x)
1156    fp12a1=x*sqrt(x)*(0.0473011_dp+y*(0.548433_dp+y))&
1157 &   /(0.0709478_dp+y*(0.737041_dp+y*0.382065_dp))
1158  end if
1159 !
1160 !**********************************************************************
1161  end function fp12a1

ABINIT/fp32a1 [ Functions ]

[ Top ] [ Functions ]

NAME

 fp32a1

FUNCTION

INPUTS

OUTPUT

SOURCE

1176  function fp32a1 (x)
1177 
1178 !Arguments -------------------------------
1179  real(dp),intent(in) :: x
1180  real(dp) :: fp32a1
1181 
1182  real(dp) :: y,x2
1183 
1184 !
1185 !**********************************************************************
1186 !*                                                                    *
1187 !*               Integrale de Fermi d'ordre 3/2                       *
1188 !*    Fp32(x) = somme de 0 a l'infini de (dt*t**3/2)/(1+exp(t-x))     *
1189 !*                                                                    *
1190 !**********************************************************************
1191 !
1192 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1193 !Erreur relative maximum annoncee 6.54 e-5
1194 !Erreur relative maximum constatee : -5.84e-5 pour eta = -5
1195 !
1196  if (x.lt.two) then
1197    y=exp(x)
1198    fp32a1=y*(135.863_dp+y*(49.2764_dp+y))/(102.210_dp+y*(55.0312_dp+y*4.23365_dp))
1199  else
1200    x2=x*x
1201    y=1._dp/x2
1202    fp32a1=x2*sqrt(x)*(0.154699_dp+y*(1.20037_dp+y))&
1203 &   /(0.386765_dp+y*(0.608119_dp-y*0.165665_dp))
1204  end if
1205 !
1206 !**********************************************************************
1207  end function fp32a1

ABINIT/ifermi12 [ Functions ]

[ Top ] [ Functions ]

NAME

 ifermi12

FUNCTION

   this routine applies a rational function expansion to get the inverse
   fermi-dirac integral of order 1/2 when it is equal to f.
   maximum error is 4.19d-9.   reference: antia apjs 84,101 1993

INPUTS

OUTPUT

SOURCE

929  function ifermi12(ff)
930 
931 !Arguments -------------------------------
932  real(dp), intent(in) :: ff
933  real(dp) :: ifermi12
934 !Local variables-------------------------------
935  integer :: ii,m1,k1,m2,k2
936  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
937 
938 !..load the coefficients of the expansion
939  data  an,m1,k1,m2,k2 /0.5e0_dp, 4, 3, 6, 5/
940  data  (a1(ii),ii=1,5)/ 1.999266880833e4_dp,   5.702479099336e3_dp,&
941 & 6.610132843877e2_dp,   3.818838129486e1_dp,&
942 & 1.0e0_dp/
943  data  (b1(ii),ii=1,4)/ 1.771804140488e4_dp,  -2.014785161019e3_dp,&
944 & 9.130355392717e1_dp,  -1.670718177489e0_dp/
945  data  (a2(ii),ii=1,7)/-1.277060388085e-2_dp,  7.187946804945e-2_dp,&
946 & -4.262314235106e-1_dp,  4.997559426872e-1_dp,&
947 & -1.285579118012e0_dp,  -3.930805454272e-1_dp,&
948 & 1.0e0_dp/
949  data  (b2(ii),ii=1,6)/-9.745794806288e-3_dp,  5.485432756838e-2_dp,&
950 & -3.299466243260e-1_dp,  4.077841975923e-1_dp,&
951 & -1.145531476975e0_dp,  -6.067091689181e-2_dp/
952 
953 ! *************************************************************************
954 
955  if (ff .lt. 4.0e0_dp) then
956    rn = ff + a1(m1)
957    do ii=m1-1,1,-1
958      rn = rn*ff + a1(ii)
959    end do
960    den = b1(k1+1)
961    do ii=k1,1,-1
962      den = den*ff + b1(ii)
963    end do
964    ifermi12 = log(ff * rn/den)
965 
966  else
967    ff1 = one/ff**(one/(one + an))
968    rn = ff1 + a2(m2)
969    do ii=m2-1,1,-1
970      rn = rn*ff1 + a2(ii)
971    end do
972    den = b2(k2+1)
973    do ii=k2,1,-1
974      den = den*ff1 + b2(ii)
975    end do
976    ifermi12 = rn/(den*ff1)
977  end if
978 
979 end function ifermi12

ABINIT/ifermi32 [ Functions ]

[ Top ] [ Functions ]

NAME

 ifermi32

FUNCTION

   this routine applies a rational function expansion to get the inverse
   fermi-dirac integral of order 3/2 when it is equal to f.
   maximum error is 2.26d-9.   reference: antia apjs 84,101 1993

INPUTS

OUTPUT

SOURCE

 997  function ifermi32(ff)
 998 
 999 !Arguments -------------------------------
1000  real(dp), intent(in) :: ff
1001  real(dp) :: ifermi32
1002 !Local variables-------------------------------
1003  integer :: ii,m1,k1,m2,k2
1004  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
1005 
1006 !..load the coefficients of the expansion
1007  data  an,m1,k1,m2,k2 /1.5e0_dp, 3, 4, 6, 5/
1008  data  (a1(ii),ii=1,4)/ 1.715627994191e2_dp,   1.125926232897e2_dp,&
1009 & 2.056296753055e1_dp,   1.0e0_dp/
1010  data  (b1(ii),ii=1,5)/ 2.280653583157e2_dp,   1.193456203021e2_dp,&
1011 & 1.167743113540e1_dp,  -3.226808804038e-1_dp,&
1012 & 3.519268762788e-3_dp/
1013  data  (a2(ii),ii=1,7)/-6.321828169799e-3_dp, -2.183147266896e-2_dp,&
1014 & -1.057562799320e-1_dp, -4.657944387545e-1_dp,&
1015 & -5.951932864088e-1_dp,  3.684471177100e-1_dp,&
1016 & 1.0e0_dp/
1017  data  (b2(ii),ii=1,6)/-4.381942605018e-3_dp, -1.513236504100e-2_dp,&
1018 & -7.850001283886e-2_dp, -3.407561772612e-1_dp,&
1019 & -5.074812565486e-1_dp, -1.387107009074e-1_dp/
1020 
1021 ! *************************************************************************
1022 
1023  if (ff .lt. 4.0e0_dp) then
1024    rn = ff + a1(m1)
1025    do ii=m1-1,1,-1
1026      rn = rn*ff + a1(ii)
1027    end do
1028    den = b1(k1+1)
1029    do ii=k1,1,-1
1030      den = den*ff + b1(ii)
1031    end do
1032    ifermi32 = log(ff * rn/den)
1033 
1034  else
1035    ff1 = one/ff**(one/(one + an))
1036    rn = ff1 + a2(m2)
1037    do ii=m2-1,1,-1
1038      rn = rn*ff1 + a2(ii)
1039    end do
1040    den = b2(k2+1)
1041    do ii=k2,1,-1
1042      den = den*ff1 + b2(ii)
1043    end do
1044    ifermi32 = rn/(den*ff1)
1045  end if
1046 
1047 end function ifermi32

ABINIT/ifermi52 [ Functions ]

[ Top ] [ Functions ]

NAME

 ifermi52

FUNCTION

   this routine applies a rational function expansion to get the inverse
   fermi-dirac integral of order 5/2 when it is equal to f.
   maximum error is 6.17d-9.   reference: antia apjs 84,101 1993

INPUTS

OUTPUT

SOURCE

1065  function ifermi52(ff)
1066 
1067 !Arguments -------------------------------
1068  real(dp), intent(in) :: ff
1069  real(dp) :: ifermi52
1070 
1071 !Local variables-------------------------------
1072  integer :: ii,m1,k1,m2,k2
1073  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
1074 
1075 !..load the coefficients of the expansion
1076  data  an,m1,k1,m2,k2 /2.5e0_dp, 2, 3, 6, 6/
1077  data  (a1(ii),ii=1,3)/ 2.138969250409e2_dp,   3.539903493971e1_dp,&
1078 & 1.0e0_dp/
1079  data  (b1(ii),ii=1,4)/ 7.108545512710e2_dp,   9.873746988121e1_dp,&
1080 & 1.067755522895e0_dp,  -1.182798726503e-2_dp/
1081  data  (a2(ii),ii=1,7)/-3.312041011227e-2_dp,  1.315763372315e-1_dp,&
1082 & -4.820942898296e-1_dp,  5.099038074944e-1_dp,&
1083 & 5.495613498630e-1_dp, -1.498867562255e0_dp,&
1084 & 1.0e0_dp/
1085  data  (b2(ii),ii=1,7)/-2.315515517515e-2_dp,  9.198776585252e-2_dp,&
1086 & -3.835879295548e-1_dp,  5.415026856351e-1_dp,&
1087 & -3.847241692193e-1_dp,  3.739781456585e-2_dp,&
1088 & -3.008504449098e-2_dp/
1089 
1090 ! *************************************************************************
1091 
1092  if (ff .lt. 4.0e0_dp) then
1093    rn = ff + a1(m1)
1094    do ii=m1-1,1,-1
1095      rn = rn*ff + a1(ii)
1096    end do
1097    den = b1(k1+1)
1098    do ii=k1,1,-1
1099      den = den*ff + b1(ii)
1100    end do
1101    ifermi52 = log(ff * rn/den)
1102 
1103  else
1104    ff1 = one/ff**(one/(one + an))
1105    rn = ff1 + a2(m2)
1106    do ii=m2-1,1,-1
1107      rn = rn*ff1 + a2(ii)
1108    end do
1109    den = b2(k2+1)
1110    do ii=k2,1,-1
1111      den = den*ff1 + b2(ii)
1112    end do
1113    ifermi52 = rn/(den*ff1)
1114  end if
1115 
1116 end function ifermi52

ABINIT/ifermim12 [ Functions ]

[ Top ] [ Functions ]

NAME

 ifermim12

FUNCTION

  this routine applies a rational function expansion to get the inverse
  fermi-dirac integral of order -1/2 when it is equal to f.
  maximum error is 3.03d-9.   reference: antia apjs 84,101 1993

INPUTS

OUTPUT

SOURCE

858  function ifermim12(ff)
859 
860 !Arguments -------------------------------
861  real(dp), intent(in) :: ff
862  real(dp) :: ifermim12
863 !Local variables-------------------------------
864  integer :: ii,m1,k1,m2,k2
865  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
866 
867 !..load the coefficients of the expansion
868  data  an,m1,k1,m2,k2 /-0.5e0_dp, 5, 6, 6, 6/
869  data  (a1(ii),ii=1,6)/-1.570044577033e4_dp,   1.001958278442e4_dp,&
870 & -2.805343454951e3_dp,   4.121170498099e2_dp,&
871 & -3.174780572961e1_dp,   1.0e0_dp/
872  data  (b1(ii),ii=1,7)/-2.782831558471e4_dp,   2.886114034012e4_dp,&
873 & -1.274243093149e4_dp,   3.063252215963e3_dp,&
874 & -4.225615045074e2_dp,   3.168918168284e1_dp,&
875 & -1.008561571363e0_dp/
876  data  (a2(ii),ii=1,7)/ 2.206779160034e-8_dp,  -1.437701234283e-6_dp,&
877 & 6.103116850636e-5_dp,  -1.169411057416e-3_dp,&
878 & 1.814141021608e-2_dp,  -9.588603457639e-2_dp,&
879 & 1.0e0_dp/
880  data  (b2(ii),ii=1,7)/ 8.827116613576e-8_dp,  -5.750804196059e-6_dp,&
881 & 2.429627688357e-4_dp,  -4.601959491394e-3_dp,&
882 & 6.932122275919e-2_dp,  -3.217372489776e-1_dp,&
883 & 3.124344749296e0_dp/
884 
885 ! *************************************************************************
886 
887  if (ff .lt. 4.0e0_dp) then
888    rn = ff + a1(m1)
889    do ii=m1-1,1,-1
890      rn = rn*ff + a1(ii)
891    end do
892    den = b1(k1+1)
893    do ii=k1,1,-1
894      den = den*ff + b1(ii)
895    end do
896    ifermim12 = log(ff * rn/den)
897 
898  else
899    ff1 = one/ff**(one/(one + an))
900    rn = ff1 + a2(m2)
901    do ii=m2-1,1,-1
902      rn = rn*ff1 + a2(ii)
903    end do
904    den = b2(k2+1)
905    do ii=k2,1,-1
906      den = den*ff1 + b2(ii)
907    end do
908    ifermim12 = rn/(den*ff1)
909  end if
910 
911 end function ifermim12

ABINIT/m_vtorhotf [ Modules ]

[ Top ] [ Modules ]

NAME

  m_vtorhotf

FUNCTION

 Computes the new density from a fixed potential (vtrial) using the Thomas-Fermi functional

COPYRIGHT

  Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, MF, AR, MM)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_vtorhotf
23 
24  use defs_basis
25  use m_abicore
26  use m_errors
27  use m_xmpi
28  use m_dtset
29 
30  use defs_abitypes, only : MPI_type
31  use m_time,     only : timab
32  use m_spacepar,  only : symrhg
33 
34  implicit none
35 
36  private

ABINIT/vtorhotf [ Functions ]

[ Top ] [ Functions ]

NAME

 vtorhotf

FUNCTION

 This routine computes the new density from a fixed potential (vtrial)
 using the Thomas-Fermi functional

INPUTS

  dtset <type(dataset_type)>=all input variables for this dataset
  gprimd(3,3)=dimensional reciprocal space primitive translations
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell.
  nfft=number of fft grid points
  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym=number of symmetry elements in space group
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  gprimd(3,3)=dimensional real space primitive translations
  ucvol=unit cell volume in bohr**3.
  vtrial(nfft,nspden)=INPUT Vtrial(r).

OUTPUT

  ek=kinetic energy part of total energy.
  enlx=nonlocal psp + potential Fock ACE part of total energy.
  entropy=entropy due to the occupation number smearing (if metal)
  fermie=fermi energy (Hartree)
  grnl(3*natom)=stores grads of nonlocal energy wrt length scales
   (3x3 tensor) and grads wrt atomic coordinates (3*natom)

SIDE EFFECTS

  rhog(2,nfft)=array for Fourier transform of electron density
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3.

SOURCE

 82 subroutine vtorhotf(dtset,ek,enlx,entropy,fermie,gprimd,grnl,&
 83 &  irrzon,mpi_enreg,natom,nfft,nspden,nsppol,nsym,phnons,rhog,rhor,rprimd,ucvol,vtrial)
 84 
 85 !Arguments ------------------------------------
 86 !scalars
 87  integer,intent(in) :: natom,nfft,nspden,nsppol,nsym
 88  real(dp),intent(in) :: ucvol
 89  real(dp),intent(out) :: ek,enlx,entropy,fermie
 90  type(MPI_type),intent(in) :: mpi_enreg
 91  type(dataset_type),intent(in) :: dtset
 92 !arrays
 93  integer,intent(in) :: irrzon((dtset%ngfft(1)*dtset%ngfft(1)*dtset%ngfft(1))**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))
 94  real(dp),intent(in) :: gprimd(3,3)
 95  real(dp),intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(1)*dtset%ngfft(1))**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))
 96  real(dp),intent(in) :: rprimd(3,3),vtrial(nfft,nspden)
 97  real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,nspden)
 98  real(dp),intent(out) :: grnl(3*natom)
 99 
100 !Local variables-------------------------------
101 !scalars
102  integer,parameter :: jdichomax=20,level=111
103  integer :: i1,i2,i3,ierr,ifft,ii,ir,iscf,jdicho
104  integer :: me_fft,n1,n2,n3,nfftot,nproc_fft,prtvol
105  real(dp),save :: cktf,fermie_tol,nelect_mid
106  real(dp) :: dnelect_mid_dx,dxrtnewt,eektemp,eektf,feektemp,feektf
107  real(dp) :: rtnewt,sum_rhor_mid,sum_rhor_middx
108  logical,save :: lfirst_time_tf=.true.
109  logical :: lnewtonraphson
110  character(len=500) :: message
111 !arrays
112  real(dp) :: tsec(2)
113  real(dp),allocatable :: betamumoinsV(:),rhor_mid(:),rhor_middx(:)
114 
115 ! *************************************************************************
116 
117 !Keep track of total time spent in vtorho
118  call timab(21,1,tsec)
119 
120 !Structured debugging if prtvol==-level
121  prtvol=dtset%prtvol
122  if(prtvol==-level)then
123    write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' vtorho : enter '
124    call wrtout(std_out,message,'COLL')
125  end if
126 
127  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
128  me_fft=dtset%ngfft(11) ; nproc_fft=dtset%ngfft(10)
129  iscf=dtset%iscf
130 !Debugging : print vtrial and rhor
131  if(prtvol==-level)then
132    write(message,'(a)') '   ir              vtrial(ir)     rhor(ir) '
133    call wrtout(std_out,message,'COLL')
134    do ir=1,nfft
135 !    if(ir<=11 .or. mod(ir,301)==0 )then
136      i3=(ir-1)/n1/(n2/nproc_fft)
137      i2=(ir-1-i3*n1*n2/nproc_fft)/n1
138      i1=ir-1-i3*n1*n2/nproc_fft-(i2-me_fft)*n1
139      write(message,'(i5,3i3,a,2d13.6)')ir,i1,i2,i3,' ',vtrial(ir,1),rhor(ir,1)
140      call wrtout(std_out,message,'COLL')
141      if(nspden>=2)then
142        write(message,'(a,2d13.6)')'               ',vtrial(ir,2),rhor(ir,2)
143        call wrtout(std_out,message,'COLL')
144      end if
145 !    end if
146    end do
147  end if
148 
149  ek=zero
150  enlx=zero
151  grnl(:)=zero
152 
153 !Initialize rhor if needed
154  if(iscf>0) rhor(:,:)=zero
155 
156 !call Thomas-Fermi for the density
157  call tf
158 !Compute energy terms
159  call tfek
160 
161  call timab(21,2,tsec)
162 !End thomas fermi
163 
164  contains

ABINIT/xp12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

 xp12a1

FUNCTION

INPUTS

OUTPUT

SOURCE

1222  function xp12a1 (y)
1223 
1224 !Arguments -------------------------------
1225  real(dp) :: xp12a1
1226  real(dp),intent(in) :: y
1227 
1228  real(dp),parameter :: deux=2._dp,deuxs3=deux/3._dp
1229  real(dp) :: top,den,z
1230 
1231 !
1232 !**********************************************************************
1233 !*                                                                    *
1234 !*              Calcul de eta tel que fp12 (eta) = y                  *
1235 !*          ou fp12 est l'integrale de Fermi d'ordre +1/2             *
1236 !*                                                                    *
1237 !**********************************************************************
1238 !
1239 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1240 !Erreur relative maximum annoncee sur exp(eta) : 3.02 e-5
1241 !
1242  if (y.lt.4._dp) then
1243    top=44.593646_dp+y*(11.288764_dp+y)
1244    den=39.519346_dp+y*(-5.7517464_dp+y*0.26594291_dp)
1245    xp12a1=log(y*top/den)
1246  else
1247    z=y**(-deuxs3)
1248    top=34.873722_dp+z*(-26.922515_dp+z)
1249    den=26.612832_dp+z*(-20.452930_dp+z*11.808945_dp)
1250    xp12a1=top/(z*den)
1251  end if
1252 !
1253 !**********************************************************************
1254  end function xp12a1

ABINIT/zfermi1 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermi1

FUNCTION

INPUTS

OUTPUT

SOURCE

480  function zfermi1(xx)
481 !..
482 !..this routine applies a rational function expansion to get the fermi-dirac
483 !..integral of order 1 evaluated at x. maximum error is 1.0e-8.
484 !..reference: antia  priv comm. 11sep94
485 !..
486 !..declare
487 
488 !Arguments -------------------------------
489  real(dp), intent(in) :: xx
490  real(dp):: zfermi1
491 !Local variables-------------------------------
492  integer ::  ii,m1,k1,m2,k2
493  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
494 
495 !..load the coefficients of the expansion
496  data  an,m1,k1,m2,k2 /1.0_dp, 7, 4, 9, 5/
497  data  (a1(ii),ii=1,8)/-7.606458638543e7_dp,  -1.143519707857e8_dp,&
498 & -5.167289383236e7_dp,  -7.304766495775e6_dp,&
499 & -1.630563622280e5_dp,   3.145920924780e3_dp,&
500 & -7.156354090495e1_dp,   1.0_dp/
501  data  (b1(ii),ii=1,5)/-7.606458639561e7_dp,  -1.333681162517e8_dp,&
502 & -7.656332234147e7_dp,  -1.638081306504e7_dp,&
503 & -1.044683266663e6_dp/
504  data (a2(ii),ii=1,10)/-3.493105157219e-7_dp, -5.628286279892e-5_dp,&
505 & -5.188757767899e-3_dp, -2.097205947730e-1_dp,&
506 & -3.353243201574_dp,    -1.682094530855e1_dp,&
507 & -2.042542575231e1_dp,   3.551366939795_dp,&
508 & -2.400826804233_dp,     1.0_dp/
509  data  (b2(ii),ii=1,6)/-6.986210315105e-7_dp, -1.102673536040e-4_dp,&
510 & -1.001475250797e-2_dp, -3.864923270059e-1_dp,&
511 & -5.435619477378_dp,    -1.563274262745e1_dp/
512 
513 ! *************************************************************************
514 
515  if (xx .lt. 2.0_dp) then
516    xx1 = exp(xx)
517    rn = xx1 + a1(m1)
518    do ii=m1-1,1,-1
519      rn = rn*xx1 + a1(ii)
520    end do
521    den = b1(k1+1)
522    do ii=k1,1,-1
523      den = den*xx1 + b1(ii)
524    end do
525    zfermi1 = xx1 * rn/den
526 
527  else
528    xx1 = 1.0_dp/(xx*xx)
529    rn = xx1 + a2(m2)
530    do ii=m2-1,1,-1
531      rn = rn*xx1 + a2(ii)
532    end do
533    den = b2(k2+1)
534    do ii=k2,1,-1
535      den = den*xx1 + b2(ii)
536    end do
537    zfermi1 = xx*xx*rn/den
538  end if
539 
540 end function zfermi1

ABINIT/zfermi12 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermi12

FUNCTION

INPUTS

OUTPUT

SOURCE

399  function zfermi12(xx)
400 !..
401 !..this routine applies a rational function expansion to get the fermi-dirac
402 !..integral of order 1/2 evaluated at x. maximum error is 5.47d-13.
403 !..reference: antia apjs 84,101 1993
404 !..
405 !..declare
406 
407 !Arguments -------------------------------
408  real(dp), intent(in) :: xx
409  real(dp):: zfermi12
410 
411 !Local variables-------------------------------
412  integer ::         ii,m1,k1,m2,k2
413  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
414 
415 !..load the coefficients of the expansion
416  data  an,m1,k1,m2,k2 /0.5e0_dp, 7, 7, 10, 11/
417  data  (a1(ii),ii=1,8)/5.75834152995465e6_dp,   1.30964880355883e7_dp,&
418 & 1.07608632249013e7_dp,   3.93536421893014e6_dp,&
419 & 6.42493233715640e5_dp,   4.16031909245777e4_dp,&
420 & 7.77238678539648e2_dp,   1.0e0_dp/
421  data  (b1(ii),ii=1,8)/6.49759261942269e6_dp,   1.70750501625775e7_dp,&
422 & 1.69288134856160e7_dp,   7.95192647756086e6_dp,&
423 & 1.83167424554505e6_dp,   1.95155948326832e5_dp,&
424 & 8.17922106644547e3_dp,   9.02129136642157e1_dp/
425  data (a2(ii),ii=1,11)/4.85378381173415e-14_dp, 1.64429113030738e-11_dp,&
426 & 3.76794942277806e-9_dp,  4.69233883900644e-7_dp,&
427 & 3.40679845803144e-5_dp,  1.32212995937796e-3_dp,&
428 & 2.60768398973913e-2_dp,  2.48653216266227e-1_dp,&
429 & 1.08037861921488e0_dp,   1.91247528779676e0_dp,&
430 & 1.0e0_dp/
431  data (b2(ii),ii=1,12)/7.28067571760518e-14_dp, 2.45745452167585e-11_dp,&
432 & 5.62152894375277e-9_dp,  6.96888634549649e-7_dp,&
433 & 5.02360015186394e-5_dp,  1.92040136756592e-3_dp,&
434 & 3.66887808002874e-2_dp,  3.24095226486468e-1_dp,&
435 & 1.16434871200131e0_dp,   1.34981244060549e0_dp,&
436 & 2.01311836975930e-1_dp, -2.14562434782759e-2_dp/
437 
438 ! *************************************************************************
439 
440  if (xx .lt. two) then
441    xx1 = exp(xx)
442    rn = xx1 + a1(m1)
443    do ii=m1-1,1,-1
444      rn = rn*xx1 + a1(ii)
445    end do
446    den = b1(k1+1)
447    do ii=k1,1,-1
448      den = den*xx1 + b1(ii)
449    end do
450    zfermi12 = xx1 * rn/den
451 
452  else
453    xx1 = one/(xx*xx)
454    rn = xx1 + a2(m2)
455    do ii=m2-1,1,-1
456      rn = rn*xx1 + a2(ii)
457    end do
458    den = b2(k2+1)
459    do ii=k2,1,-1
460      den = den*xx1 + b2(ii)
461    end do
462    zfermi12 = xx*sqrt(xx)*rn/den
463  end if
464 
465 end function zfermi12

ABINIT/zfermi2 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermi2

FUNCTION

  this routine applies a rational function expansion to get the fermi-dirac
  integral of order 2 evaluated at x. maximum error is 1.0e-8.
  reference: antia  priv comm. 11sep94

INPUTS

OUTPUT

SOURCE

636  function zfermi2(xx)
637 
638 !Arguments -------------------------------
639  real(dp), intent(in) :: xx
640  real(dp) :: zfermi2
641 !Local variables-------------------------------
642  integer ::  ii,m1,k1,m2,k2
643  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
644 
645 !..load the coefficients of the expansion
646  data  an,m1,k1,m2,k2 /2.0_dp, 7, 4, 5, 9/
647  data  (a1(ii),ii=1,8)/-1.434885992395e8_dp,  -2.001711155617e8_dp,&
648 & -8.507067153428e7_dp,  -1.175118281976e7_dp,&
649 & -3.145120854293e5_dp,   4.275771034579e3_dp,&
650 & -8.069902926891e1_dp,   1.0e0_dp/
651  data  (b1(ii),ii=1,5)/-7.174429962316e7_dp,  -1.090535948744e8_dp,&
652 & -5.350984486022e7_dp,  -9.646265123816e6_dp,&
653 & -5.113415562845e5_dp/
654  data  (a2(ii),ii=1,6)/ 6.919705180051e-8_dp,  1.134026972699e-5_dp,&
655 & 7.967092675369e-4_dp,  2.432500578301e-2_dp,&
656 & 2.784751844942e-1_dp,  1.0e0_dp/
657  data (b2(ii),ii=1,10)/ 2.075911553728e-7_dp,  3.197196691324e-5_dp,&
658 & 2.074576609543e-3_dp,  5.250009686722e-2_dp,&
659 & 3.171705130118e-1_dp, -1.147237720706e-1_dp,&
660 & 6.638430718056e-2_dp, -1.356814647640e-2_dp,&
661 & -3.648576227388e-2_dp,  3.621098757460e-2_dp/
662 
663 ! *************************************************************************
664 
665  if (xx .lt. 2.0e0_dp) then
666    xx1 = exp(xx)
667    rn = xx1 + a1(m1)
668    do ii=m1-1,1,-1
669      rn = rn*xx1 + a1(ii)
670    end do
671    den = b1(k1+1)
672    do ii=k1,1,-1
673      den = den*xx1 + b1(ii)
674    end do
675    zfermi2 = xx1 * rn/den
676 
677  else
678    xx1 = one/(xx*xx)
679    rn = xx1 + a2(m2)
680    do ii=m2-1,1,-1
681      rn = rn*xx1 + a2(ii)
682    end do
683    den = b2(k2+1)
684    do ii=k2,1,-1
685      den = den*xx1 + b2(ii)
686    end do
687    zfermi2 = xx*xx*xx*rn/den
688  end if
689 
690 end function zfermi2

ABINIT/zfermi3 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermi3

FUNCTION

  this routine applies a rational function expansion to get the fermi-dirac
  integral of order 3 evaluated at x. maximum error is 1.0e-8.
  reference: antia  priv comm. 11sep94

INPUTS

OUTPUT

SOURCE

785  function zfermi3(xx)
786 
787 !Arguments -------------------------------
788  real(dp), intent(in) :: xx
789  real(dp):: zfermi3
790 
791 !Local variables-------------------------------
792  integer :: ii,m1,k1,m2,k2
793  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
794 
795 !..load the coefficients of the expansion
796  data  an,m1,k1,m2,k2 /3.0, 4, 6, 7, 7/
797  data  (a1(ii),ii=1,5)/ 6.317036716422e2_dp,    7.514163924637e2_dp,&
798 & 2.711961035750e2_dp,    3.274540902317e1_dp,&
799 & 1.0_dp/
800  data  (b1(ii),ii=1,7)/ 1.052839452797e2_dp,    1.318163114785e2_dp,&
801 & 5.213807524405e1_dp,    7.500064111991_dp,&
802 & 3.383020205492e-1_dp,   2.342176749453e-3_dp,&
803 & -8.445226098359e-6_dp/
804  data  (a2(ii),ii=1,8)/ 1.360999428425e-8_dp,   1.651419468084e-6_dp,&
805 & 1.021455604288e-4_dp,   3.041270709839e-3_dp,&
806 & 4.584298418374e-2_dp,   3.440523212512e-1_dp,&
807 & 1.077505444383_dp,    1.0_dp/
808  data  (b2(ii),ii=1,8)/ 5.443997714076e-8_dp,   5.531075760054e-6_dp,&
809 & 2.969285281294e-4_dp,   6.052488134435e-3_dp,&
810 & 5.041144894964e-2_dp,   1.048282487684e-1_dp,&
811 & 1.280969214096e-2_dp,  -2.851555446444e-3_dp/
812 
813 ! *************************************************************************
814 
815  if (xx .lt. two) then
816    xx1 = exp(xx)
817    rn = xx1 + a1(m1)
818    do ii=m1-1,1,-1
819      rn = rn*xx1 + a1(ii)
820    end do
821    den = b1(k1+1)
822    do ii=k1,1,-1
823      den = den*xx1 + b1(ii)
824    end do
825    zfermi3 = xx1 * rn/den
826 
827  else
828    xx1 = one/(xx*xx)
829    rn = xx1 + a2(m2)
830    do ii=m2-1,1,-1
831      rn = rn*xx1 + a2(ii)
832    end do
833    den = b2(k2+1)
834    do ii=k2,1,-1
835      den = den*xx1 + b2(ii)
836    end do
837    zfermi3 = xx*xx*xx*xx*rn/den
838  end if
839 
840 end function zfermi3

ABINIT/zfermi32 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermi32

FUNCTION

  this routine applies a rational function expansion to get the fermi-dirac
  integral of order 3/2 evaluated at x. maximum error is 5.07d-13.
  reference: antia apjs 84,101 1993

INPUTS

OUTPUT

SOURCE

558  function zfermi32(xx)
559 
560 !Arguments -------------------------------
561  real(dp), intent(in) :: xx
562  real(dp) :: zfermi32
563 
564 !Local variables-------------------------------
565  integer :: ii,m1,k1,m2,k2
566  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
567 
568 !..load the coefficients of the expansion
569  data  an,m1,k1,m2,k2 /1.5e0_dp, 6, 7, 9, 10/
570  data  (a1(ii),ii=1,7)/4.32326386604283e4_dp,   8.55472308218786e4_dp,&
571 & 5.95275291210962e4_dp,   1.77294861572005e4_dp,&
572 & 2.21876607796460e3_dp,   9.90562948053193e1_dp,&
573 & 1.0e0_dp/
574  data  (b1(ii),ii=1,8)/3.25218725353467e4_dp,   7.01022511904373e4_dp,&
575 & 5.50859144223638e4_dp,   1.95942074576400e4_dp,&
576 & 3.20803912586318e3_dp,   2.20853967067789e2_dp,&
577 & 5.05580641737527e0_dp,   1.99507945223266e-2_dp/
578  data (a2(ii),ii=1,10)/2.80452693148553e-13_dp, 8.60096863656367e-11_dp,&
579 & 1.62974620742993e-8_dp,  1.63598843752050e-6_dp,&
580 & 9.12915407846722e-5_dp,  2.62988766922117e-3_dp,&
581 & 3.85682997219346e-2_dp,  2.78383256609605e-1_dp,&
582 & 9.02250179334496e-1_dp,  1.0e0_dp/
583  data (b2(ii),ii=1,11)/7.01131732871184e-13_dp, 2.10699282897576e-10_dp,&
584 & 3.94452010378723e-8_dp,  3.84703231868724e-6_dp,&
585 & 2.04569943213216e-4_dp,  5.31999109566385e-3_dp,&
586 & 6.39899717779153e-2_dp,  3.14236143831882e-1_dp,&
587 & 4.70252591891375e-1_dp, -2.15540156936373e-2_dp,&
588 & 2.34829436438087e-3_dp/
589 
590 ! *************************************************************************
591 
592  if (xx .lt. 2.0e0_dp) then
593    xx1 = exp(xx)
594    rn = xx1 + a1(m1)
595    do ii=m1-1,1,-1
596      rn = rn*xx1 + a1(ii)
597    end do
598    den = b1(k1+1)
599    do ii=k1,1,-1
600      den = den*xx1 + b1(ii)
601    end do
602    zfermi32 = xx1 * rn/den
603 
604  else
605    xx1 = one/(xx*xx)
606    rn = xx1 + a2(m2)
607    do ii=m2-1,1,-1
608      rn = rn*xx1 + a2(ii)
609    end do
610    den = b2(k2+1)
611    do ii=k2,1,-1
612      den = den*xx1 + b2(ii)
613    end do
614    zfermi32 = xx*xx*sqrt(xx)*rn/den
615  end if
616 
617 end function zfermi32

ABINIT/zfermi52 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermi52

FUNCTION

  this routine applies a rational function expansion to get the fermi-dirac
  integral of order 5/2 evaluated at x. maximum error is 2.47d-13.
  reference: antia apjs 84,101 1993

INPUTS

OUTPUT

SOURCE

708  function zfermi52(xx)
709 
710 !Arguments -------------------------------
711  real(dp), intent(in) :: xx
712  real(dp) :: zfermi52
713 
714 !Local variables-------------------------------
715  integer :: ii,m1,k1,m2,k2
716  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
717 
718 !..load the coefficients of the expansion
719  data  an,m1,k1,m2,k2 /2.5e0_dp, 6, 7, 10, 9/
720  data  (a1(ii),ii=1,7)/6.61606300631656e4_dp,   1.20132462801652e5_dp,&
721 & 7.67255995316812e4_dp,   2.10427138842443e4_dp,&
722 & 2.44325236813275e3_dp,   1.02589947781696e2_dp,&
723 & 1.0e0_dp/
724  data  (b1(ii),ii=1,8)/1.99078071053871e4_dp,   3.79076097261066e4_dp,&
725 & 2.60117136841197e4_dp,   7.97584657659364e3_dp,&
726 & 1.10886130159658e3_dp,   6.35483623268093e1_dp,&
727 & 1.16951072617142e0_dp,   3.31482978240026e-3_dp/
728  data (a2(ii),ii=1,11)/8.42667076131315e-12_dp, 2.31618876821567e-9_dp,&
729 & 3.54323824923987e-7_dp,  2.77981736000034e-5_dp,&
730 & 1.14008027400645e-3_dp,  2.32779790773633e-2_dp,&
731 & 2.39564845938301e-1_dp,  1.24415366126179e0_dp,&
732 & 3.18831203950106e0_dp,   3.42040216997894e0_dp,&
733 & 1.0e0_dp/
734  data (b2(ii),ii=1,10)/2.94933476646033e-11_dp, 7.68215783076936e-9_dp,&
735 & 1.12919616415947e-6_dp,  8.09451165406274e-5_dp,&
736 & 2.81111224925648e-3_dp,  3.99937801931919e-2_dp,&
737 & 2.27132567866839e-1_dp,  5.31886045222680e-1_dp,&
738 & 3.70866321410385e-1_dp,  2.27326643192516e-2_dp/
739 
740 ! *************************************************************************
741 
742  if (xx .lt. two) then
743    xx1 = exp(xx)
744    rn = xx1 + a1(m1)
745    do ii=m1-1,1,-1
746      rn = rn*xx1 + a1(ii)
747    end do
748    den = b1(k1+1)
749    do ii=k1,1,-1
750      den = den*xx1 + b1(ii)
751    end do
752    zfermi52 = xx1 * rn/den
753 
754  else
755    xx1 = one/(xx*xx)
756    rn = xx1 + a2(m2)
757    do ii=m2-1,1,-1
758      rn = rn*xx1 + a2(ii)
759    end do
760    den = b2(k2+1)
761    do ii=k2,1,-1
762      den = den*xx1 + b2(ii)
763    end do
764    zfermi52 = xx*xx*xx*sqrt(xx)*rn/den
765  end if
766 
767 end function zfermi52

ABINIT/zfermim12 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermim12

FUNCTION



INPUTS


OUTPUT


SOURCE


325  function zfermim12(xx)
326 
327 !Arguments -------------------------------
328  real(dp), intent(in) :: xx
329  real(dp) :: zfermim12
330 
331 !Local variables-------------------------------
332  integer ::  ii,m1,k1,m2,k2
333  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
334 !..load the coefficients of the expansion
335  data  an,m1,k1,m2,k2 /-0.5e0_dp, 7, 7, 11, 11/
336  data  (a1(ii),ii=1,8)/ 1.71446374704454e7_dp,    3.88148302324068e7_dp,&
337 & 3.16743385304962e7_dp,    1.14587609192151e7_dp,&
338 & 1.83696370756153E6_dp,    1.14980998186874e5_dp,&
339 & 1.98276889924768e3_dp,    1.0e0_dp/
340  data  (b1(ii),ii=1,8)/ 9.67282587452899e6_dp,    2.87386436731785e7_dp,&
341 & 3.26070130734158e7_dp,    1.77657027846367e7_dp,&
342 & 4.81648022267831e6_dp,    6.13709569333207e5_dp,&
343 & 3.13595854332114e4_dp,    4.35061725080755e2_dp/
344  data (a2(ii),ii=1,12)/-4.46620341924942e-15_dp, -1.58654991146236e-12_dp,&
345 & -4.44467627042232e-10_dp, -6.84738791621745e-8_dp,&
346 & -6.64932238528105e-6_dp,  -3.69976170193942e-4_dp,&
347 & -1.12295393687006e-2_dp,  -1.60926102124442e-1_dp,&
348 & -8.52408612877447e-1_dp,  -7.45519953763928e-1_dp,&
349 & 2.98435207466372e0_dp,    1.0e0_dp/
350  data (b2(ii),ii=1,12)/-2.23310170962369e-15_dp, -7.94193282071464e-13_dp,&
351 & -2.22564376956228e-10_dp, -3.43299431079845e-8_dp,&
352 & -3.33919612678907e-6_dp,  -1.86432212187088e-4_dp,&
353 & -5.69764436880529e-3_dp,  -8.34904593067194e-2_dp,&
354 & -4.78770844009440e-1_dp,  -4.99759250374148e-1_dp,&
355 & 1.86795964993052e0_dp,    4.16485970495288e-1_dp/
356 
357 ! *************************************************************************
358 
359  if (xx .lt. 2.0e0_dp) then
360    xx1 = exp(xx)
361    rn = xx1 + a1(m1)
362    do ii=m1-1,1,-1
363      rn = rn*xx1 + a1(ii)
364    end do
365    den = b1(k1+1)
366    do ii=k1,1,-1
367      den = den*xx1 + b1(ii)
368    end do
369    zfermim12 = xx1 * rn/den
370 !  ..
371  else
372    xx1 = one/(xx*xx)
373    rn = xx1 + a2(m2)
374    do ii=m2-1,1,-1
375      rn = rn*xx1 + a2(ii)
376    end do
377    den = b2(k2+1)
378    do ii=k2,1,-1
379      den = den*xx1 + b2(ii)
380    end do
381    zfermim12 = sqrt(xx)*rn/den
382  end if
383 
384 end function zfermim12







vtorhotf/tf [ Functions ]



[ Top ] [ vtorhotf ] [ Functions ]


NAME


 tf



FUNCTION


INPUTS


OUTPUT


SOURCE


178   subroutine tf()
179 
180 ! *************************************************************************
181 
182    ABI_MALLOC(rhor_mid,(nfft))
183    ABI_MALLOC(rhor_middx,(nfft))
184    fermie_tol=1.e-10_dp
185    cktf=one/two/pi**2*(two*dtset%tphysel)**1.5_dp
186 
187 !  Should be made an input variable, if TF really needed for production
188 !  rtnewt=dtset%userra
189    rtnewt=zero
190 
191 !  Newton Raphson
192    if (lfirst_time_tf) then
193      lfirst_time_tf=.false.
194    end if
195    jdicho=0
196    lnewtonraphson=.false.
197    do while (.not.lnewtonraphson)
198      jdicho=jdicho+1
199 !    do ifft=1,nfft
200 !    rhor_mid(ifft)=cktf*zfermi12((rtnewt-vtrial(ifft,1))/dtset%tphysel)
201 !    rhor_middx(ifft)=cktf*zfermim12((rtnewt-vtrial(ifft,1))/dtset%tphysel)
202 !    end do
203      call fm12a1t(cktf,rtnewt,dtset%tphysel,vtrial(:,1),rhor_middx,rhor_mid,&
204 &     nfft)
205      sum_rhor_mid=sum(rhor_mid(:))
206      sum_rhor_middx=sum(rhor_middx(:))
207      call xmpi_sum(sum_rhor_mid,mpi_enreg%comm_fft ,ierr)
208      call xmpi_sum(sum_rhor_middx,mpi_enreg%comm_fft ,ierr)
209      nelect_mid=sum_rhor_mid*ucvol/(nfft*nproc_fft)-dtset%nelect
210      dnelect_mid_dx=sum_rhor_middx*ucvol/(nfft*nproc_fft)/dtset%tphysel/two
211      dxrtnewt=nelect_mid/dnelect_mid_dx
212      rtnewt=rtnewt-dxrtnewt
213      if (abs(nelect_mid) < fermie_tol/2._dp) then
214        lnewtonraphson=.true.
215      end if
216      if (jdicho > jdichomax) then
217        ABI_ERROR('NEWTON RAPHSON NOT CONVERGED')
218      end if
219    end do
220    fermie=rtnewt
221    rhor(:,1)=rhor_mid(:)
222    ABI_FREE(rhor_mid)
223    ABI_FREE(rhor_middx)
224 
225 !  DEBUG
226 !  write(std_out,*)'fmid,nmid,jdicho',fermie,nelect_mid,jdicho
227 !  ENDDEBUG
228 
229 !  Compute rhog
230    call timab(70,1,tsec)
231 
232    nfftot=dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)
233    call symrhg(1,gprimd,irrzon,mpi_enreg,nfft,nfftot,dtset%ngfft,nspden,nsppol,nsym,phnons,&
234 &   rhog,rhor,rprimd,dtset%symafm,dtset%symrel,dtset%tnons)
235 
236 !  We now have both rho(r) and rho(G), symmetrized, and if nsppol=2
237 !  we also have the spin-up density, symmetrized, in rhor(:,2).
238 
239    call timab(70,2,tsec)
240 
241 end subroutine tf







vtorhotf/tfek [ Functions ]



[ Top ] [ vtorhotf ] [ Functions ]


NAME


 tfek



FUNCTION


 This is the calculation of the kinetic energy for Thomas Fermi
 Energy and free energy must be distinguished



INPUTS


OUTPUT


SOURCE


258   subroutine tfek()
259 
260 ! *************************************************************************
261 
262    ABI_MALLOC(betamumoinsV,(nfft))
263    cktf=one/two/pi**2*(two*dtset%tphysel)**1.5_dp
264    eektf=zero
265    feektf=zero
266    do ifft=1,nfft
267 
268 !    betamumoinsv(ifft)=ifermi12(rhor(ifft,1)/cktf)
269      betamumoinsv(ifft)=(rtnewt-vtrial(ifft,1))/dtset%tphysel
270 !    eektemp=zfermi32(betamumoinsV(ifft))/zfermi12(betamumoinsV(ifft))
271      eektemp=fp32a1(betamumoinsV(ifft))/rhor(ifft,1)*cktf
272      feektemp=betamumoinsV(ifft)-two/three*eektemp
273      feektf=feektf+feektemp*rhor(ifft,1)
274      eektf=eektf+eektemp*rhor(ifft,1)
275    end do
276 !  Init mpi_comm
277    call timab(48,1,tsec)
278    call xmpi_sum(eektf,mpi_enreg%comm_fft ,ierr)
279    call xmpi_sum(feektf,mpi_enreg%comm_fft ,ierr)
280    call timab(48,2,tsec)
281    eektf=eektf*dtset%tphysel
282    eektf=eektf*ucvol/dble(nfft*nproc_fft)
283    feektf=feektf*dtset%tphysel
284    feektf=feektf*ucvol/dble(nfft*nproc_fft)
285 !  DEBUG
286 !  write(std_out,*)'eektf',eektf
287 !  stop ('vtorhotf')
288 !  ENDDEBUG
289    ek=eektf
290    entropy=(eektf-feektf)/dtset%tphysel
291    ABI_FREE(betamumoinsV)
292  end subroutine tfek