TABLE OF CONTENTS
- ABINIT/fm12a1
- ABINIT/fm12a1t
- ABINIT/fp12a1
- ABINIT/fp32a1
- ABINIT/ifermi12
- ABINIT/ifermi32
- ABINIT/ifermi52
- ABINIT/ifermim12
- ABINIT/m_vtorhotf
- ABINIT/vtorhotf
- ABINIT/xp12a1
- ABINIT/zfermi1
- ABINIT/zfermi12
- ABINIT/zfermi2
- ABINIT/zfermi3
- ABINIT/zfermi32
- ABINIT/zfermi52
- ABINIT/zfermim12
- vtorhotf/tf
- vtorhotf/tfek
ABINIT/fm12a1 [ 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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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