## ABINIT/fm12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

``` fm12a1
```

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

```1509  function fm12a1 (x)
1510
1511
1512 !This section has been created automatically by the script Abilint (TD).
1513 !Do not modify the following lines by hand.
1514 #undef ABI_FUNC
1515 #define ABI_FUNC 'fm12a1'
1516 !End of the abilint section
1517
1518  implicit none
1519
1520 !Arguments -------------------------------
1521  real(dp),intent(in) :: x
1522  real(dp) :: fm12a1
1523
1524  real(dp) :: y
1525
1526 !
1527 !**********************************************************************
1528 !*                                                                    *
1529 !*               Integrale de Fermi d'ordre -1/2                      *
1530 !*    Fm12(x) = somme de 0 a l'infini de (dt*t**-1/2)/(1+exp(t-x))    *
1531 !*                                                                    *
1532 !**********************************************************************
1533 !
1534 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1535 !Erreur relative maximum annoncee 4.75 e-5
1536 !
1537  if (x.lt.2._dp) then
1538    y=exp(x)
1539    fm12a1=y*(23.1456_dp+y*(13.7820_dp+y))&
1540 &   /(13.0586_dp+y*(17.0048_dp+y*(5.07527_dp+y*0.236620_dp)))
1541  else
1542    y=1./(x*x)
1543    fm12a1=sqrt(x)*(0.0153602_dp+y*(0.146815_dp+y))&
1544 &   /(0.00768015_dp+y*(0.0763700_dp+y*0.570485_dp))
1545  end if
1546 !
1547 !**********************************************************************
1548  end function fm12a1
```

## ABINIT/fm12a1t [ Functions ]

[ Top ] [ Functions ]

NAME

``` fm12a1t
```

FUNCTION

INPUTS

OUTPUT

PARENTS

```      vtorhotf
```

CHILDREN

SOURCE

```1568  subroutine fm12a1t (cktf,rtnewt,tphysel,vtrial,rhor_middx,rhor_mid,nfft)
1569
1570
1571 !This section has been created automatically by the script Abilint (TD).
1572 !Do not modify the following lines by hand.
1573 #undef ABI_FUNC
1574 #define ABI_FUNC 'fm12a1t'
1575 !End of the abilint section
1576
1577  implicit none
1578
1579  integer,intent(in) :: nfft
1580  real(dp),intent(in) :: tphysel,rtnewt,cktf
1581  real(dp),intent(in) :: vtrial(nfft)
1582  real(dp),intent(out) :: rhor_middx(nfft),rhor_mid(nfft)
1583
1584  !intrinsic exp,sqrt
1585  integer :: ifft
1586  real(dp) :: x,y,sqrtx
1587
1588 !
1589 !**********************************************************************
1590 !*                                                                    *
1591 !*               Integrale de Fermi d'ordre -1/2                      *
1592 !*    Fm12(x) = somme de 0 a l'infini de (dt*t**-1/2)/(1+exp(t-x))    *
1593 !*                      ....                                              *
1594 !**********************************************************************
1595 !
1596 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1597 !Erreur relative maximum annoncee 4.75 e-5
1598 !
1599  do ifft=1,nfft
1600    x=(rtnewt-vtrial(ifft))/tphysel
1601    if (x.lt.2._dp) then
1602      y=exp(x)
1603      rhor_middx(ifft)=cktf*y*(23.1456e0_dp+y*(13.7820e0_dp+y))&
1604 &     /(13.0586e0_dp+y*(17.0048e0_dp+y*(5.07527e0_dp+y*0.236620e0_dp)))
1605      rhor_mid(ifft)=cktf*y*(21.8168_dp+y*(13.1693_dp+y))&
1606 &     /(24.6180+y*(23.5546_dp+y*(4.76290_dp+y*0.134481_dp)))
1607    else
1608      y=1._dp/(x*x)
1609      sqrtx=sqrt(x)
1610      rhor_middx(ifft)=cktf*sqrtx*(0.0153602e0_dp+y*(0.146815e0_dp+y))&
1611 &     /(0.00768015e0_dp+y*(0.0763700e0_dp+y*0.570485e0_dp))
1612      rhor_mid(ifft)=cktf*x*sqrtx*(0.0473011_dp+y*(0.548433_dp+y))&
1613 &     /(0.0709478_dp+y*(0.737041_dp+y*0.382065_dp))
1614    end if
1615  end do
1616 !
1617 !**********************************************************************
1618  end subroutine fm12a1t
```

## ABINIT/fp12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

``` fp12a1
```

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

```1332  function fp12a1 (x)
1333
1334
1335 !This section has been created automatically by the script Abilint (TD).
1336 !Do not modify the following lines by hand.
1337 #undef ABI_FUNC
1338 #define ABI_FUNC 'fp12a1'
1339 !End of the abilint section
1340
1341  implicit none
1342
1343 ! Arguments -------------------------------
1344  real(dp),intent(in) :: x
1345  real(dp) :: fp12a1
1346
1347  real(dp) :: y
1348
1349 !**********************************************************************
1350 !*                                                                    *
1351 !*               Integrale de Fermi d'ordre 1/2                       *
1352 !*    Fp12(x) = somme de 0 a l'infini de (dt*t**1/2)/(1+exp(t-x))     *
1353 !*                                                                    *
1354 !**********************************************************************
1355 !
1356 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1357 !Erreur relative maximum annoncee 5.54 e-5
1358 !Erreur relative maximum constatee : -5.53e-5 pour eta = 2
1359 !
1360  if (x.lt.2._dp) then
1361    y=exp(x)
1362    fp12a1=y*(21.8168_dp+y*(13.1693_dp+y))&
1363 &   /(24.6180_dp+y*(23.5546_dp+y*(4.76290_dp+y*0.134481_dp)))
1364  else
1365    y=one/(x*x)
1366    fp12a1=x*sqrt(x)*(0.0473011_dp+y*(0.548433_dp+y))&
1367 &   /(0.0709478_dp+y*(0.737041_dp+y*0.382065_dp))
1368  end if
1369 !
1370 !**********************************************************************
1371  end function fp12a1
```

## ABINIT/fp32a1 [ Functions ]

[ Top ] [ Functions ]

NAME

``` fp32a1
```

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

```1390  function fp32a1 (x)
1391
1392
1393 !This section has been created automatically by the script Abilint (TD).
1394 !Do not modify the following lines by hand.
1395 #undef ABI_FUNC
1396 #define ABI_FUNC 'fp32a1'
1397 !End of the abilint section
1398
1399  implicit none
1400
1401 !Arguments -------------------------------
1402  real(dp),intent(in) :: x
1403  real(dp) :: fp32a1
1404
1405  real(dp) :: y,x2
1406
1407 !
1408 !**********************************************************************
1409 !*                                                                    *
1410 !*               Integrale de Fermi d'ordre 3/2                       *
1411 !*    Fp32(x) = somme de 0 a l'infini de (dt*t**3/2)/(1+exp(t-x))     *
1412 !*                                                                    *
1413 !**********************************************************************
1414 !
1415 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1416 !Erreur relative maximum annoncee 6.54 e-5
1417 !Erreur relative maximum constatee : -5.84e-5 pour eta = -5
1418 !
1419  if (x.lt.two) then
1420    y=exp(x)
1421    fp32a1=y*(135.863_dp+y*(49.2764_dp+y))/(102.210_dp+y*(55.0312_dp+y*4.23365_dp))
1422  else
1423    x2=x*x
1424    y=1._dp/x2
1425    fp32a1=x2*sqrt(x)*(0.154699_dp+y*(1.20037_dp+y))&
1426 &   /(0.386765_dp+y*(0.608119_dp-y*0.165665_dp))
1427  end if
1428 !
1429 !**********************************************************************
1430  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

PARENTS

CHILDREN

SOURCE

```1091  function ifermi12(ff)
1092
1093
1094 !This section has been created automatically by the script Abilint (TD).
1095 !Do not modify the following lines by hand.
1096 #undef ABI_FUNC
1097 #define ABI_FUNC 'ifermi12'
1098 !End of the abilint section
1099
1100  implicit none
1101
1102 !Arguments -------------------------------
1103  real(dp), intent(in) :: ff
1104  real(dp) :: ifermi12
1105 !Local variables-------------------------------
1106  integer :: ii,m1,k1,m2,k2
1107  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
1108
1109 !..load the coefficients of the expansion
1110  data  an,m1,k1,m2,k2 /0.5e0_dp, 4, 3, 6, 5/
1111  data  (a1(ii),ii=1,5)/ 1.999266880833e4_dp,   5.702479099336e3_dp,&
1112 & 6.610132843877e2_dp,   3.818838129486e1_dp,&
1113 & 1.0e0_dp/
1114  data  (b1(ii),ii=1,4)/ 1.771804140488e4_dp,  -2.014785161019e3_dp,&
1115 & 9.130355392717e1_dp,  -1.670718177489e0_dp/
1116  data  (a2(ii),ii=1,7)/-1.277060388085e-2_dp,  7.187946804945e-2_dp,&
1117 & -4.262314235106e-1_dp,  4.997559426872e-1_dp,&
1118 & -1.285579118012e0_dp,  -3.930805454272e-1_dp,&
1119 & 1.0e0_dp/
1120  data  (b2(ii),ii=1,6)/-9.745794806288e-3_dp,  5.485432756838e-2_dp,&
1121 & -3.299466243260e-1_dp,  4.077841975923e-1_dp,&
1122 & -1.145531476975e0_dp,  -6.067091689181e-2_dp/
1123
1124 ! *************************************************************************
1125
1126  if (ff .lt. 4.0e0_dp) then
1127    rn = ff + a1(m1)
1128    do ii=m1-1,1,-1
1129      rn = rn*ff + a1(ii)
1130    end do
1131    den = b1(k1+1)
1132    do ii=k1,1,-1
1133      den = den*ff + b1(ii)
1134    end do
1135    ifermi12 = log(ff * rn/den)
1136
1137  else
1138    ff1 = one/ff**(one/(one + an))
1139    rn = ff1 + a2(m2)
1140    do ii=m2-1,1,-1
1141      rn = rn*ff1 + a2(ii)
1142    end do
1143    den = b2(k2+1)
1144    do ii=k2,1,-1
1145      den = den*ff1 + b2(ii)
1146    end do
1147    ifermi12 = rn/(den*ff1)
1148  end if
1149
1150 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

PARENTS

CHILDREN

SOURCE

```1172  function ifermi32(ff)
1173
1174
1175 !This section has been created automatically by the script Abilint (TD).
1176 !Do not modify the following lines by hand.
1177 #undef ABI_FUNC
1178 #define ABI_FUNC 'ifermi32'
1179 !End of the abilint section
1180
1181  implicit none
1182
1183 !Arguments -------------------------------
1184  real(dp), intent(in) :: ff
1185  real(dp) :: ifermi32
1186 !Local variables-------------------------------
1187  integer :: ii,m1,k1,m2,k2
1188  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
1189
1190 !..load the coefficients of the expansion
1191  data  an,m1,k1,m2,k2 /1.5e0_dp, 3, 4, 6, 5/
1192  data  (a1(ii),ii=1,4)/ 1.715627994191e2_dp,   1.125926232897e2_dp,&
1193 & 2.056296753055e1_dp,   1.0e0_dp/
1194  data  (b1(ii),ii=1,5)/ 2.280653583157e2_dp,   1.193456203021e2_dp,&
1195 & 1.167743113540e1_dp,  -3.226808804038e-1_dp,&
1196 & 3.519268762788e-3_dp/
1197  data  (a2(ii),ii=1,7)/-6.321828169799e-3_dp, -2.183147266896e-2_dp,&
1198 & -1.057562799320e-1_dp, -4.657944387545e-1_dp,&
1199 & -5.951932864088e-1_dp,  3.684471177100e-1_dp,&
1200 & 1.0e0_dp/
1201  data  (b2(ii),ii=1,6)/-4.381942605018e-3_dp, -1.513236504100e-2_dp,&
1202 & -7.850001283886e-2_dp, -3.407561772612e-1_dp,&
1203 & -5.074812565486e-1_dp, -1.387107009074e-1_dp/
1204
1205 ! *************************************************************************
1206
1207  if (ff .lt. 4.0e0_dp) then
1208    rn = ff + a1(m1)
1209    do ii=m1-1,1,-1
1210      rn = rn*ff + a1(ii)
1211    end do
1212    den = b1(k1+1)
1213    do ii=k1,1,-1
1214      den = den*ff + b1(ii)
1215    end do
1216    ifermi32 = log(ff * rn/den)
1217
1218  else
1219    ff1 = one/ff**(one/(one + an))
1220    rn = ff1 + a2(m2)
1221    do ii=m2-1,1,-1
1222      rn = rn*ff1 + a2(ii)
1223    end do
1224    den = b2(k2+1)
1225    do ii=k2,1,-1
1226      den = den*ff1 + b2(ii)
1227    end do
1228    ifermi32 = rn/(den*ff1)
1229  end if
1230
1231 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

PARENTS

CHILDREN

SOURCE

```1253  function ifermi52(ff)
1254
1255
1256 !This section has been created automatically by the script Abilint (TD).
1257 !Do not modify the following lines by hand.
1258 #undef ABI_FUNC
1259 #define ABI_FUNC 'ifermi52'
1260 !End of the abilint section
1261
1262  implicit none
1263
1264 !Arguments -------------------------------
1265  real(dp), intent(in) :: ff
1266  real(dp) :: ifermi52
1267
1268 !Local variables-------------------------------
1269  integer :: ii,m1,k1,m2,k2
1270  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
1271
1272 !..load the coefficients of the expansion
1273  data  an,m1,k1,m2,k2 /2.5e0_dp, 2, 3, 6, 6/
1274  data  (a1(ii),ii=1,3)/ 2.138969250409e2_dp,   3.539903493971e1_dp,&
1275 & 1.0e0_dp/
1276  data  (b1(ii),ii=1,4)/ 7.108545512710e2_dp,   9.873746988121e1_dp,&
1277 & 1.067755522895e0_dp,  -1.182798726503e-2_dp/
1278  data  (a2(ii),ii=1,7)/-3.312041011227e-2_dp,  1.315763372315e-1_dp,&
1279 & -4.820942898296e-1_dp,  5.099038074944e-1_dp,&
1280 & 5.495613498630e-1_dp, -1.498867562255e0_dp,&
1281 & 1.0e0_dp/
1282  data  (b2(ii),ii=1,7)/-2.315515517515e-2_dp,  9.198776585252e-2_dp,&
1283 & -3.835879295548e-1_dp,  5.415026856351e-1_dp,&
1284 & -3.847241692193e-1_dp,  3.739781456585e-2_dp,&
1285 & -3.008504449098e-2_dp/
1286
1287 ! *************************************************************************
1288
1289  if (ff .lt. 4.0e0_dp) then
1290    rn = ff + a1(m1)
1291    do ii=m1-1,1,-1
1292      rn = rn*ff + a1(ii)
1293    end do
1294    den = b1(k1+1)
1295    do ii=k1,1,-1
1296      den = den*ff + b1(ii)
1297    end do
1298    ifermi52 = log(ff * rn/den)
1299
1300  else
1301    ff1 = one/ff**(one/(one + an))
1302    rn = ff1 + a2(m2)
1303    do ii=m2-1,1,-1
1304      rn = rn*ff1 + a2(ii)
1305    end do
1306    den = b2(k2+1)
1307    do ii=k2,1,-1
1308      den = den*ff1 + b2(ii)
1309    end do
1310    ifermi52 = rn/(den*ff1)
1311  end if
1312
1313 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

PARENTS

CHILDREN

SOURCE

```1007  function ifermim12(ff)
1008
1009
1010 !This section has been created automatically by the script Abilint (TD).
1011 !Do not modify the following lines by hand.
1012 #undef ABI_FUNC
1013 #define ABI_FUNC 'ifermim12'
1014 !End of the abilint section
1015
1016  implicit none
1017
1018 !Arguments -------------------------------
1019  real(dp), intent(in) :: ff
1020  real(dp) :: ifermim12
1021 !Local variables-------------------------------
1022  integer :: ii,m1,k1,m2,k2
1023  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
1024
1025 !..load the coefficients of the expansion
1026  data  an,m1,k1,m2,k2 /-0.5e0_dp, 5, 6, 6, 6/
1027  data  (a1(ii),ii=1,6)/-1.570044577033e4_dp,   1.001958278442e4_dp,&
1028 & -2.805343454951e3_dp,   4.121170498099e2_dp,&
1029 & -3.174780572961e1_dp,   1.0e0_dp/
1030  data  (b1(ii),ii=1,7)/-2.782831558471e4_dp,   2.886114034012e4_dp,&
1031 & -1.274243093149e4_dp,   3.063252215963e3_dp,&
1032 & -4.225615045074e2_dp,   3.168918168284e1_dp,&
1033 & -1.008561571363e0_dp/
1034  data  (a2(ii),ii=1,7)/ 2.206779160034e-8_dp,  -1.437701234283e-6_dp,&
1035 & 6.103116850636e-5_dp,  -1.169411057416e-3_dp,&
1036 & 1.814141021608e-2_dp,  -9.588603457639e-2_dp,&
1037 & 1.0e0_dp/
1038  data  (b2(ii),ii=1,7)/ 8.827116613576e-8_dp,  -5.750804196059e-6_dp,&
1039 & 2.429627688357e-4_dp,  -4.601959491394e-3_dp,&
1040 & 6.932122275919e-2_dp,  -3.217372489776e-1_dp,&
1041 & 3.124344749296e0_dp/
1042
1043 ! *************************************************************************
1044
1045  if (ff .lt. 4.0e0_dp) then
1046    rn = ff + a1(m1)
1047    do ii=m1-1,1,-1
1048      rn = rn*ff + a1(ii)
1049    end do
1050    den = b1(k1+1)
1051    do ii=k1,1,-1
1052      den = den*ff + b1(ii)
1053    end do
1054    ifermim12 = log(ff * rn/den)
1055
1056  else
1057    ff1 = one/ff**(one/(one + an))
1058    rn = ff1 + a2(m2)
1059    do ii=m2-1,1,-1
1060      rn = rn*ff1 + a2(ii)
1061    end do
1062    den = b2(k2+1)
1063    do ii=k2,1,-1
1064      den = den*ff1 + b2(ii)
1065    end do
1066    ifermim12 = rn/(den*ff1)
1067  end if
1068
1069 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 (C) 1998-2018 ABINIT group (DCA, XG, GMR, MF, AR, MM)
GNU General Public License, see ~abinit/COPYING
or http://www.gnu.org/copyleft/gpl.txt .
```

PARENTS

CHILDREN

SOURCE

```21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24
25 #include "abi_common.h"
26
27 module m_vtorhotf
28
29  use defs_basis
30  use defs_abitypes
31  use m_abicore
32  use m_errors
33  use m_xmpi
34
35  use m_time,     only : timab
36  use m_dtfil,    only : status
37  use m_spacepar,  only : symrhg
38
39  implicit none
40
41  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

```  dtfil <type(datafiles_type)>=variables related to files
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
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.
enl=nonlocal pseudopotential 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.
```

PARENTS

```      scfcv
```

CHILDREN

SOURCE

``` 93 subroutine vtorhotf(dtfil,dtset,ek,enl,entropy,fermie,gprimd,grnl,&
94 &  irrzon,mpi_enreg,natom,nfft,nspden,nsppol,nsym,phnons,rhog,rhor,rprimd,ucvol,vtrial)
95
96
97 !This section has been created automatically by the script Abilint (TD).
98 !Do not modify the following lines by hand.
99 #undef ABI_FUNC
100 #define ABI_FUNC 'vtorhotf'
101 !End of the abilint section
102
103  implicit none
104
105 !Arguments ------------------------------------
106 !scalars
107  integer,intent(in) :: natom,nfft,nspden,nsppol,nsym
108  real(dp),intent(in) :: ucvol
109  real(dp),intent(out) :: ek,enl,entropy,fermie
110  type(MPI_type),intent(in) :: mpi_enreg
111  type(datafiles_type),intent(in) :: dtfil
112  type(dataset_type),intent(in) :: dtset
113 !arrays
114  integer,intent(in) :: irrzon((dtset%ngfft(1)*dtset%ngfft(1)*dtset%ngfft(1))**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))
115  real(dp),intent(in) :: gprimd(3,3)
116  real(dp),intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(1)*dtset%ngfft(1))**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))
117  real(dp),intent(in) :: rprimd(3,3),vtrial(nfft,nspden)
118  real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,nspden)
119  real(dp),intent(out) :: grnl(3*natom)
120
121 !Local variables-------------------------------
122 !scalars
123  integer,parameter :: jdichomax=20,level=111
124  integer :: i1,i2,i3,ierr,iexit,ifft,ii,ir,iscf,jdicho
125  integer :: me_fft,n1,n2,n3,nfftot,nproc_fft,prtvol
126  real(dp),save :: cktf,fermie_tol,nelect_mid
127  real(dp) :: dnelect_mid_dx,dxrtnewt,eektemp,eektf,feektemp,feektf
128  real(dp) :: rtnewt,sum_rhor_mid,sum_rhor_middx
129  logical,save :: lfirst_time_tf=.true.
130  logical :: lnewtonraphson
131  character(len=500) :: message
132 !arrays
133  real(dp) :: tsec(2)
134  real(dp),allocatable :: betamumoinsV(:),rhor_mid(:),rhor_middx(:)
135 !no_abirules
136
137 ! *************************************************************************
138
139 !Keep track of total time spent in vtorho
140  call timab(21,1,tsec)
141
142
143  call status(0,dtfil%filstat,iexit,level,'enter         ')
144
145 !Structured debugging if prtvol==-level
146  prtvol=dtset%prtvol
147  if(prtvol==-level)then
148    write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' vtorho : enter '
149    call wrtout(std_out,message,'COLL')
150  end if
151
152  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
153  me_fft=dtset%ngfft(11) ; nproc_fft=dtset%ngfft(10)
154  iscf=dtset%iscf
155 !Debugging : print vtrial and rhor
156  if(prtvol==-level)then
157    write(message,'(a)') '   ir              vtrial(ir)     rhor(ir) '
158    call wrtout(std_out,message,'COLL')
159    do ir=1,nfft
160 !    if(ir<=11 .or. mod(ir,301)==0 )then
161      i3=(ir-1)/n1/(n2/nproc_fft)
162      i2=(ir-1-i3*n1*n2/nproc_fft)/n1
163      i1=ir-1-i3*n1*n2/nproc_fft-(i2-me_fft)*n1
164      write(message,'(i5,3i3,a,2d13.6)')ir,i1,i2,i3,' ',vtrial(ir,1),rhor(ir,1)
165      call wrtout(std_out,message,'COLL')
166      if(nspden>=2)then
167        write(message,'(a,2d13.6)')'               ',vtrial(ir,2),rhor(ir,2)
168        call wrtout(std_out,message,'COLL')
169      end if
170 !    end if
171    end do
172  end if
173
174  ek=zero
175  enl=zero
176  grnl(:)=zero
177
178 !Initialize rhor if needed
179  if(iscf>0) rhor(:,:)=zero
180
181 !call Thomas-Fermi for the density
182  call tf
183 !Compute energy terms
184  call tfek
185
186  call status(0,dtfil%filstat,iexit,level,'exit          ')
187
188  call timab(21,2,tsec)
189 !End thomas fermi
190
191  contains
```

## ABINIT/xp12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

``` xp12a1
```

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

```1449  function xp12a1 (y)
1450
1451
1452 !This section has been created automatically by the script Abilint (TD).
1453 !Do not modify the following lines by hand.
1454 #undef ABI_FUNC
1455 #define ABI_FUNC 'xp12a1'
1456 !End of the abilint section
1457
1458  implicit none
1459
1460 !Arguments -------------------------------
1461  real(dp) :: xp12a1
1462  real(dp),intent(in) :: y
1463
1464  real(dp),parameter :: deux=2._dp,deuxs3=deux/3._dp
1465  real(dp) :: top,den,z
1466
1467 !
1468 !**********************************************************************
1469 !*                                                                    *
1470 !*              Calcul de eta tel que fp12 (eta) = y                  *
1471 !*          ou fp12 est l'integrale de Fermi d'ordre +1/2             *
1472 !*                                                                    *
1473 !**********************************************************************
1474 !
1475 !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1476 !Erreur relative maximum annoncee sur exp(eta) : 3.02 e-5
1477 !
1478  if (y.lt.4._dp) then
1479    top=44.593646_dp+y*(11.288764_dp+y)
1480    den=39.519346_dp+y*(-5.7517464_dp+y*0.26594291_dp)
1481    xp12a1=log(y*top/den)
1482  else
1483    z=y**(-deuxs3)
1484    top=34.873722_dp+z*(-26.922515_dp+z)
1485    den=26.612832_dp+z*(-20.452930_dp+z*11.808945_dp)
1486    xp12a1=top/(z*den)
1487  end if
1488 !
1489 !**********************************************************************
1490  end function xp12a1
```

## ABINIT/zfermi1 [ Functions ]

[ Top ] [ Functions ]

NAME

``` zfermi1
```

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

```565  function zfermi1(xx)
566 !..
567 !..this routine applies a rational function expansion to get the fermi-dirac
568 !..integral of order 1 evaluated at x. maximum error is 1.0e-8.
569 !..reference: antia  priv comm. 11sep94
570 !..
571 !..declare
572
573 !This section has been created automatically by the script Abilint (TD).
574 !Do not modify the following lines by hand.
575 #undef ABI_FUNC
576 #define ABI_FUNC 'zfermi1'
577 !End of the abilint section
578
579  implicit none
580
581 !Arguments -------------------------------
582  real(dp), intent(in) :: xx
583  real(dp):: zfermi1
584 !Local variables-------------------------------
585  integer ::  ii,m1,k1,m2,k2
586  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
587
588 !..load the coefficients of the expansion
589  data  an,m1,k1,m2,k2 /1.0_dp, 7, 4, 9, 5/
590  data  (a1(ii),ii=1,8)/-7.606458638543e7_dp,  -1.143519707857e8_dp,&
591 & -5.167289383236e7_dp,  -7.304766495775e6_dp,&
592 & -1.630563622280e5_dp,   3.145920924780e3_dp,&
593 & -7.156354090495e1_dp,   1.0_dp/
594  data  (b1(ii),ii=1,5)/-7.606458639561e7_dp,  -1.333681162517e8_dp,&
595 & -7.656332234147e7_dp,  -1.638081306504e7_dp,&
596 & -1.044683266663e6_dp/
597  data (a2(ii),ii=1,10)/-3.493105157219e-7_dp, -5.628286279892e-5_dp,&
598 & -5.188757767899e-3_dp, -2.097205947730e-1_dp,&
599 & -3.353243201574_dp,    -1.682094530855e1_dp,&
600 & -2.042542575231e1_dp,   3.551366939795_dp,&
601 & -2.400826804233_dp,     1.0_dp/
602  data  (b2(ii),ii=1,6)/-6.986210315105e-7_dp, -1.102673536040e-4_dp,&
603 & -1.001475250797e-2_dp, -3.864923270059e-1_dp,&
604 & -5.435619477378_dp,    -1.563274262745e1_dp/
605
606 ! *************************************************************************
607
608  if (xx .lt. 2.0_dp) then
609    xx1 = exp(xx)
610    rn = xx1 + a1(m1)
611    do ii=m1-1,1,-1
612      rn = rn*xx1 + a1(ii)
613    end do
614    den = b1(k1+1)
615    do ii=k1,1,-1
616      den = den*xx1 + b1(ii)
617    end do
618    zfermi1 = xx1 * rn/den
619
620  else
621    xx1 = 1.0_dp/(xx*xx)
622    rn = xx1 + a2(m2)
623    do ii=m2-1,1,-1
624      rn = rn*xx1 + a2(ii)
625    end do
626    den = b2(k2+1)
627    do ii=k2,1,-1
628      den = den*xx1 + b2(ii)
629    end do
630    zfermi1 = xx*xx*rn/den
631  end if
632
633 end function zfermi1
```

## ABINIT/zfermi12 [ Functions ]

[ Top ] [ Functions ]

NAME

``` zfermi12
```

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

```472  function zfermi12(xx)
473 !..
474 !..this routine applies a rational function expansion to get the fermi-dirac
475 !..integral of order 1/2 evaluated at x. maximum error is 5.47d-13.
476 !..reference: antia apjs 84,101 1993
477 !..
478 !..declare
479
480 !This section has been created automatically by the script Abilint (TD).
481 !Do not modify the following lines by hand.
482 #undef ABI_FUNC
483 #define ABI_FUNC 'zfermi12'
484 !End of the abilint section
485
486  implicit none
487
488 !Arguments -------------------------------
489  real(dp), intent(in) :: xx
490  real(dp):: zfermi12
491
492 !Local variables-------------------------------
493  integer ::         ii,m1,k1,m2,k2
494  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
495
496 !..load the coefficients of the expansion
497  data  an,m1,k1,m2,k2 /0.5e0_dp, 7, 7, 10, 11/
498  data  (a1(ii),ii=1,8)/5.75834152995465e6_dp,   1.30964880355883e7_dp,&
499 & 1.07608632249013e7_dp,   3.93536421893014e6_dp,&
500 & 6.42493233715640e5_dp,   4.16031909245777e4_dp,&
501 & 7.77238678539648e2_dp,   1.0e0_dp/
502  data  (b1(ii),ii=1,8)/6.49759261942269e6_dp,   1.70750501625775e7_dp,&
503 & 1.69288134856160e7_dp,   7.95192647756086e6_dp,&
504 & 1.83167424554505e6_dp,   1.95155948326832e5_dp,&
505 & 8.17922106644547e3_dp,   9.02129136642157e1_dp/
506  data (a2(ii),ii=1,11)/4.85378381173415e-14_dp, 1.64429113030738e-11_dp,&
507 & 3.76794942277806e-9_dp,  4.69233883900644e-7_dp,&
508 & 3.40679845803144e-5_dp,  1.32212995937796e-3_dp,&
509 & 2.60768398973913e-2_dp,  2.48653216266227e-1_dp,&
510 & 1.08037861921488e0_dp,   1.91247528779676e0_dp,&
511 & 1.0e0_dp/
512  data (b2(ii),ii=1,12)/7.28067571760518e-14_dp, 2.45745452167585e-11_dp,&
513 & 5.62152894375277e-9_dp,  6.96888634549649e-7_dp,&
514 & 5.02360015186394e-5_dp,  1.92040136756592e-3_dp,&
515 & 3.66887808002874e-2_dp,  3.24095226486468e-1_dp,&
516 & 1.16434871200131e0_dp,   1.34981244060549e0_dp,&
517 & 2.01311836975930e-1_dp, -2.14562434782759e-2_dp/
518
519 ! *************************************************************************
520
521  if (xx .lt. two) then
522    xx1 = exp(xx)
523    rn = xx1 + a1(m1)
524    do ii=m1-1,1,-1
525      rn = rn*xx1 + a1(ii)
526    end do
527    den = b1(k1+1)
528    do ii=k1,1,-1
529      den = den*xx1 + b1(ii)
530    end do
531    zfermi12 = xx1 * rn/den
532
533  else
534    xx1 = one/(xx*xx)
535    rn = xx1 + a2(m2)
536    do ii=m2-1,1,-1
537      rn = rn*xx1 + a2(ii)
538    end do
539    den = b2(k2+1)
540    do ii=k2,1,-1
541      den = den*xx1 + b2(ii)
542    end do
543    zfermi12 = xx*sqrt(xx)*rn/den
544  end if
545
546 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

PARENTS

CHILDREN

SOURCE

```746  function zfermi2(xx)
747
748
749 !This section has been created automatically by the script Abilint (TD).
750 !Do not modify the following lines by hand.
751 #undef ABI_FUNC
752 #define ABI_FUNC 'zfermi2'
753 !End of the abilint section
754
755  implicit none
756
757 !Arguments -------------------------------
758  real(dp), intent(in) :: xx
759  real(dp) :: zfermi2
760 !Local variables-------------------------------
761  integer ::  ii,m1,k1,m2,k2
762  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
763
764 !..load the coefficients of the expansion
765  data  an,m1,k1,m2,k2 /2.0_dp, 7, 4, 5, 9/
766  data  (a1(ii),ii=1,8)/-1.434885992395e8_dp,  -2.001711155617e8_dp,&
767 & -8.507067153428e7_dp,  -1.175118281976e7_dp,&
768 & -3.145120854293e5_dp,   4.275771034579e3_dp,&
769 & -8.069902926891e1_dp,   1.0e0_dp/
770  data  (b1(ii),ii=1,5)/-7.174429962316e7_dp,  -1.090535948744e8_dp,&
771 & -5.350984486022e7_dp,  -9.646265123816e6_dp,&
772 & -5.113415562845e5_dp/
773  data  (a2(ii),ii=1,6)/ 6.919705180051e-8_dp,  1.134026972699e-5_dp,&
774 & 7.967092675369e-4_dp,  2.432500578301e-2_dp,&
775 & 2.784751844942e-1_dp,  1.0e0_dp/
776  data (b2(ii),ii=1,10)/ 2.075911553728e-7_dp,  3.197196691324e-5_dp,&
777 & 2.074576609543e-3_dp,  5.250009686722e-2_dp,&
778 & 3.171705130118e-1_dp, -1.147237720706e-1_dp,&
779 & 6.638430718056e-2_dp, -1.356814647640e-2_dp,&
780 & -3.648576227388e-2_dp,  3.621098757460e-2_dp/
781
782 ! *************************************************************************
783
784  if (xx .lt. 2.0e0_dp) then
785    xx1 = exp(xx)
786    rn = xx1 + a1(m1)
787    do ii=m1-1,1,-1
788      rn = rn*xx1 + a1(ii)
789    end do
790    den = b1(k1+1)
791    do ii=k1,1,-1
792      den = den*xx1 + b1(ii)
793    end do
794    zfermi2 = xx1 * rn/den
795
796  else
797    xx1 = one/(xx*xx)
798    rn = xx1 + a2(m2)
799    do ii=m2-1,1,-1
800      rn = rn*xx1 + a2(ii)
801    end do
802    den = b2(k2+1)
803    do ii=k2,1,-1
804      den = den*xx1 + b2(ii)
805    end do
806    zfermi2 = xx*xx*xx*rn/den
807  end if
808
809 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

PARENTS

CHILDREN

SOURCE

```921  function zfermi3(xx)
922
923
924 !This section has been created automatically by the script Abilint (TD).
925 !Do not modify the following lines by hand.
926 #undef ABI_FUNC
927 #define ABI_FUNC 'zfermi3'
928 !End of the abilint section
929
930  implicit none
931
932 !Arguments -------------------------------
933  real(dp), intent(in) :: xx
934  real(dp):: zfermi3
935
936 !Local variables-------------------------------
937  integer :: ii,m1,k1,m2,k2
938  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
939
940 !..load the coefficients of the expansion
941  data  an,m1,k1,m2,k2 /3.0, 4, 6, 7, 7/
942  data  (a1(ii),ii=1,5)/ 6.317036716422e2_dp,    7.514163924637e2_dp,&
943 & 2.711961035750e2_dp,    3.274540902317e1_dp,&
944 & 1.0_dp/
945  data  (b1(ii),ii=1,7)/ 1.052839452797e2_dp,    1.318163114785e2_dp,&
946 & 5.213807524405e1_dp,    7.500064111991_dp,&
947 & 3.383020205492e-1_dp,   2.342176749453e-3_dp,&
948 & -8.445226098359e-6_dp/
949  data  (a2(ii),ii=1,8)/ 1.360999428425e-8_dp,   1.651419468084e-6_dp,&
950 & 1.021455604288e-4_dp,   3.041270709839e-3_dp,&
951 & 4.584298418374e-2_dp,   3.440523212512e-1_dp,&
952 & 1.077505444383_dp,    1.0_dp/
953  data  (b2(ii),ii=1,8)/ 5.443997714076e-8_dp,   5.531075760054e-6_dp,&
954 & 2.969285281294e-4_dp,   6.052488134435e-3_dp,&
955 & 5.041144894964e-2_dp,   1.048282487684e-1_dp,&
956 & 1.280969214096e-2_dp,  -2.851555446444e-3_dp/
957
958 ! *************************************************************************
959
960  if (xx .lt. two) then
961    xx1 = exp(xx)
962    rn = xx1 + a1(m1)
963    do ii=m1-1,1,-1
964      rn = rn*xx1 + a1(ii)
965    end do
966    den = b1(k1+1)
967    do ii=k1,1,-1
968      den = den*xx1 + b1(ii)
969    end do
970    zfermi3 = xx1 * rn/den
971
972  else
973    xx1 = one/(xx*xx)
974    rn = xx1 + a2(m2)
975    do ii=m2-1,1,-1
976      rn = rn*xx1 + a2(ii)
977    end do
978    den = b2(k2+1)
979    do ii=k2,1,-1
980      den = den*xx1 + b2(ii)
981    end do
982    zfermi3 = xx*xx*xx*xx*rn/den
983  end if
984
985 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

PARENTS

CHILDREN

SOURCE

```655  function zfermi32(xx)
656
657
658 !This section has been created automatically by the script Abilint (TD).
659 !Do not modify the following lines by hand.
660 #undef ABI_FUNC
661 #define ABI_FUNC 'zfermi32'
662 !End of the abilint section
663
664  implicit none
665
666 !Arguments -------------------------------
667  real(dp), intent(in) :: xx
668  real(dp) :: zfermi32
669
670 !Local variables-------------------------------
671  integer :: ii,m1,k1,m2,k2
672  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
673
674 !..load the coefficients of the expansion
675  data  an,m1,k1,m2,k2 /1.5e0_dp, 6, 7, 9, 10/
676  data  (a1(ii),ii=1,7)/4.32326386604283e4_dp,   8.55472308218786e4_dp,&
677 & 5.95275291210962e4_dp,   1.77294861572005e4_dp,&
678 & 2.21876607796460e3_dp,   9.90562948053193e1_dp,&
679 & 1.0e0_dp/
680  data  (b1(ii),ii=1,8)/3.25218725353467e4_dp,   7.01022511904373e4_dp,&
681 & 5.50859144223638e4_dp,   1.95942074576400e4_dp,&
682 & 3.20803912586318e3_dp,   2.20853967067789e2_dp,&
683 & 5.05580641737527e0_dp,   1.99507945223266e-2_dp/
684  data (a2(ii),ii=1,10)/2.80452693148553e-13_dp, 8.60096863656367e-11_dp,&
685 & 1.62974620742993e-8_dp,  1.63598843752050e-6_dp,&
686 & 9.12915407846722e-5_dp,  2.62988766922117e-3_dp,&
687 & 3.85682997219346e-2_dp,  2.78383256609605e-1_dp,&
688 & 9.02250179334496e-1_dp,  1.0e0_dp/
689  data (b2(ii),ii=1,11)/7.01131732871184e-13_dp, 2.10699282897576e-10_dp,&
690 & 3.94452010378723e-8_dp,  3.84703231868724e-6_dp,&
691 & 2.04569943213216e-4_dp,  5.31999109566385e-3_dp,&
692 & 6.39899717779153e-2_dp,  3.14236143831882e-1_dp,&
693 & 4.70252591891375e-1_dp, -2.15540156936373e-2_dp,&
694 & 2.34829436438087e-3_dp/
695
696 ! *************************************************************************
697
698  if (xx .lt. 2.0e0_dp) then
699    xx1 = exp(xx)
700    rn = xx1 + a1(m1)
701    do ii=m1-1,1,-1
702      rn = rn*xx1 + a1(ii)
703    end do
704    den = b1(k1+1)
705    do ii=k1,1,-1
706      den = den*xx1 + b1(ii)
707    end do
708    zfermi32 = xx1 * rn/den
709
710  else
711    xx1 = one/(xx*xx)
712    rn = xx1 + a2(m2)
713    do ii=m2-1,1,-1
714      rn = rn*xx1 + a2(ii)
715    end do
716    den = b2(k2+1)
717    do ii=k2,1,-1
718      den = den*xx1 + b2(ii)
719    end do
720    zfermi32 = xx*xx*sqrt(xx)*rn/den
721  end if
722
723 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

PARENTS

CHILDREN

SOURCE

```831  function zfermi52(xx)
832
833
834 !This section has been created automatically by the script Abilint (TD).
835 !Do not modify the following lines by hand.
836 #undef ABI_FUNC
837 #define ABI_FUNC 'zfermi52'
838 !End of the abilint section
839
840  implicit none
841
842 !Arguments -------------------------------
843  real(dp), intent(in) :: xx
844  real(dp) :: zfermi52
845
846 !Local variables-------------------------------
847  integer :: ii,m1,k1,m2,k2
848  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
849
850 !..load the coefficients of the expansion
851  data  an,m1,k1,m2,k2 /2.5e0_dp, 6, 7, 10, 9/
852  data  (a1(ii),ii=1,7)/6.61606300631656e4_dp,   1.20132462801652e5_dp,&
853 & 7.67255995316812e4_dp,   2.10427138842443e4_dp,&
854 & 2.44325236813275e3_dp,   1.02589947781696e2_dp,&
855 & 1.0e0_dp/
856  data  (b1(ii),ii=1,8)/1.99078071053871e4_dp,   3.79076097261066e4_dp,&
857 & 2.60117136841197e4_dp,   7.97584657659364e3_dp,&
858 & 1.10886130159658e3_dp,   6.35483623268093e1_dp,&
859 & 1.16951072617142e0_dp,   3.31482978240026e-3_dp/
860  data (a2(ii),ii=1,11)/8.42667076131315e-12_dp, 2.31618876821567e-9_dp,&
861 & 3.54323824923987e-7_dp,  2.77981736000034e-5_dp,&
862 & 1.14008027400645e-3_dp,  2.32779790773633e-2_dp,&
863 & 2.39564845938301e-1_dp,  1.24415366126179e0_dp,&
864 & 3.18831203950106e0_dp,   3.42040216997894e0_dp,&
865 & 1.0e0_dp/
866  data (b2(ii),ii=1,10)/2.94933476646033e-11_dp, 7.68215783076936e-9_dp,&
867 & 1.12919616415947e-6_dp,  8.09451165406274e-5_dp,&
868 & 2.81111224925648e-3_dp,  3.99937801931919e-2_dp,&
869 & 2.27132567866839e-1_dp,  5.31886045222680e-1_dp,&
870 & 3.70866321410385e-1_dp,  2.27326643192516e-2_dp/
871
872 ! *************************************************************************
873
874  if (xx .lt. two) then
875    xx1 = exp(xx)
876    rn = xx1 + a1(m1)
877    do ii=m1-1,1,-1
878      rn = rn*xx1 + a1(ii)
879    end do
880    den = b1(k1+1)
881    do ii=k1,1,-1
882      den = den*xx1 + b1(ii)
883    end do
884    zfermi52 = xx1 * rn/den
885
886  else
887    xx1 = one/(xx*xx)
888    rn = xx1 + a2(m2)
889    do ii=m2-1,1,-1
890      rn = rn*xx1 + a2(ii)
891    end do
892    den = b2(k2+1)
893    do ii=k2,1,-1
894      den = den*xx1 + b2(ii)
895    end do
896    zfermi52 = xx*xx*xx*sqrt(xx)*rn/den
897  end if
898
899 end function zfermi52
```

## ABINIT/zfermim12 [ Functions ]

[ Top ] [ Functions ]

NAME

``` zfermim12
```

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

```385  function zfermim12(xx)
386
387
388 !This section has been created automatically by the script Abilint (TD).
389 !Do not modify the following lines by hand.
390 #undef ABI_FUNC
391 #define ABI_FUNC 'zfermim12'
392 !End of the abilint section
393
394  implicit none
395
396 !Arguments -------------------------------
397  real(dp), intent(in) :: xx
398  real(dp) :: zfermim12
399
400 !Local variables-------------------------------
401  integer ::  ii,m1,k1,m2,k2
402  real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
403 !..load the coefficients of the expansion
404  data  an,m1,k1,m2,k2 /-0.5e0_dp, 7, 7, 11, 11/
405  data  (a1(ii),ii=1,8)/ 1.71446374704454e7_dp,    3.88148302324068e7_dp,&
406 & 3.16743385304962e7_dp,    1.14587609192151e7_dp,&
407 & 1.83696370756153E6_dp,    1.14980998186874e5_dp,&
408 & 1.98276889924768e3_dp,    1.0e0_dp/
409  data  (b1(ii),ii=1,8)/ 9.67282587452899e6_dp,    2.87386436731785e7_dp,&
410 & 3.26070130734158e7_dp,    1.77657027846367e7_dp,&
411 & 4.81648022267831e6_dp,    6.13709569333207e5_dp,&
412 & 3.13595854332114e4_dp,    4.35061725080755e2_dp/
413  data (a2(ii),ii=1,12)/-4.46620341924942e-15_dp, -1.58654991146236e-12_dp,&
414 & -4.44467627042232e-10_dp, -6.84738791621745e-8_dp,&
415 & -6.64932238528105e-6_dp,  -3.69976170193942e-4_dp,&
416 & -1.12295393687006e-2_dp,  -1.60926102124442e-1_dp,&
417 & -8.52408612877447e-1_dp,  -7.45519953763928e-1_dp,&
418 & 2.98435207466372e0_dp,    1.0e0_dp/
419  data (b2(ii),ii=1,12)/-2.23310170962369e-15_dp, -7.94193282071464e-13_dp,&
420 & -2.22564376956228e-10_dp, -3.43299431079845e-8_dp,&
421 & -3.33919612678907e-6_dp,  -1.86432212187088e-4_dp,&
422 & -5.69764436880529e-3_dp,  -8.34904593067194e-2_dp,&
423 & -4.78770844009440e-1_dp,  -4.99759250374148e-1_dp,&
424 & 1.86795964993052e0_dp,    4.16485970495288e-1_dp/
425
426 ! *************************************************************************
427
428  if (xx .lt. 2.0e0_dp) then
429    xx1 = exp(xx)
430    rn = xx1 + a1(m1)
431    do ii=m1-1,1,-1
432      rn = rn*xx1 + a1(ii)
433    end do
434    den = b1(k1+1)
435    do ii=k1,1,-1
436      den = den*xx1 + b1(ii)
437    end do
438    zfermim12 = xx1 * rn/den
439 !  ..
440  else
441    xx1 = one/(xx*xx)
442    rn = xx1 + a2(m2)
443    do ii=m2-1,1,-1
444      rn = rn*xx1 + a2(ii)
445    end do
446    den = b2(k2+1)
447    do ii=k2,1,-1
448      den = den*xx1 + b2(ii)
449    end do
450    zfermim12 = sqrt(xx)*rn/den
451  end if
452
453 end function zfermim12
```

## vtorhotf/tf [ Functions ]

[ Top ] [ vtorhotf ] [ Functions ]

NAME

``` tf
```

FUNCTION

INPUTS

OUTPUT

PARENTS

```      vtorhotf
```

CHILDREN

SOURCE

```210   subroutine tf()
211
212
213 !This section has been created automatically by the script Abilint (TD).
214 !Do not modify the following lines by hand.
215 #undef ABI_FUNC
216 #define ABI_FUNC 'tf'
217 !End of the abilint section
218
219   implicit none
220
221 ! *************************************************************************
222
223    ABI_ALLOCATE(rhor_mid,(nfft))
224    ABI_ALLOCATE(rhor_middx,(nfft))
225    fermie_tol=1.e-10_dp
226    cktf=one/two/pi**2*(two*dtset%tphysel)**1.5_dp
227
228 !  Should be made an input variable, if TF really needed for production
229 !  rtnewt=dtset%userra
230    rtnewt=zero
231
232 !  Newton Raphson
233    if (lfirst_time_tf) then
234      lfirst_time_tf=.false.
235    end if
236    jdicho=0
237    lnewtonraphson=.false.
238    do while (.not.lnewtonraphson)
239      jdicho=jdicho+1
240 !    do ifft=1,nfft
241 !    rhor_mid(ifft)=cktf*zfermi12((rtnewt-vtrial(ifft,1))/dtset%tphysel)
242 !    rhor_middx(ifft)=cktf*zfermim12((rtnewt-vtrial(ifft,1))/dtset%tphysel)
243 !    end do
244      call fm12a1t(cktf,rtnewt,dtset%tphysel,vtrial(:,1),rhor_middx,rhor_mid,&
245 &     nfft)
246      sum_rhor_mid=sum(rhor_mid(:))
247      sum_rhor_middx=sum(rhor_middx(:))
248      call xmpi_sum(sum_rhor_mid,mpi_enreg%comm_fft ,ierr)
249      call xmpi_sum(sum_rhor_middx,mpi_enreg%comm_fft ,ierr)
250      nelect_mid=sum_rhor_mid*ucvol/(nfft*nproc_fft)-dtset%nelect
251      dnelect_mid_dx=sum_rhor_middx*ucvol/(nfft*nproc_fft)/dtset%tphysel/two
252      dxrtnewt=nelect_mid/dnelect_mid_dx
253      rtnewt=rtnewt-dxrtnewt
254      if (abs(nelect_mid) < fermie_tol/2._dp) then
255        lnewtonraphson=.true.
256      end if
257      if (jdicho > jdichomax) then
258        MSG_ERROR('NEWTON RAPHSON NOT CONVERGED')
259      end if
260    end do
261    fermie=rtnewt
262    rhor(:,1)=rhor_mid(:)
263    ABI_DEALLOCATE(rhor_mid)
264    ABI_DEALLOCATE(rhor_middx)
265
266 !  DEBUG
267 !  write(std_out,*)'fmid,nmid,jdicho',fermie,nelect_mid,jdicho
268 !  ENDDEBUG
269
270 !  Compute rhog
271    call timab(70,1,tsec)
272
273    call status(0,dtfil%filstat,iexit,level,'compute rhog  ')
274    nfftot=dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)
275    call symrhg(1,gprimd,irrzon,mpi_enreg,nfft,nfftot,dtset%ngfft,nspden,nsppol,nsym,dtset%paral_kgb,phnons,&
276 &   rhog,rhor,rprimd,dtset%symafm,dtset%symrel)
277
278 !  We now have both rho(r) and rho(G), symmetrized, and if nsppol=2
279 !  we also have the spin-up density, symmetrized, in rhor(:,2).
280
281    call timab(70,2,tsec)
282
283 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

PARENTS

```      vtorhotf
```

CHILDREN

SOURCE

```305   subroutine tfek()
306
307
308 !This section has been created automatically by the script Abilint (TD).
309 !Do not modify the following lines by hand.
310 #undef ABI_FUNC
311 #define ABI_FUNC 'tfek'
312 !End of the abilint section
313
314   implicit none
315
316 ! *************************************************************************
317
318    ABI_ALLOCATE(betamumoinsV,(nfft))
319    cktf=one/two/pi**2*(two*dtset%tphysel)**1.5_dp
320    eektf=zero
321    feektf=zero
322    do ifft=1,nfft
323
324 !    betamumoinsv(ifft)=ifermi12(rhor(ifft,1)/cktf)
325      betamumoinsv(ifft)=(rtnewt-vtrial(ifft,1))/dtset%tphysel
326 !    eektemp=zfermi32(betamumoinsV(ifft))/zfermi12(betamumoinsV(ifft))
327      eektemp=fp32a1(betamumoinsV(ifft))/rhor(ifft,1)*cktf
328      feektemp=betamumoinsV(ifft)-two/three*eektemp
329      feektf=feektf+feektemp*rhor(ifft,1)
330      eektf=eektf+eektemp*rhor(ifft,1)
331    end do
332 !  Init mpi_comm
333    call timab(48,1,tsec)
334    call xmpi_sum(eektf,mpi_enreg%comm_fft ,ierr)
335    call xmpi_sum(feektf,mpi_enreg%comm_fft ,ierr)
336    call timab(48,2,tsec)
337    eektf=eektf*dtset%tphysel
338    eektf=eektf*ucvol/dble(nfft*nproc_fft)
339    feektf=feektf*dtset%tphysel
340    feektf=feektf*ucvol/dble(nfft*nproc_fft)
341 !  DEBUG
342 !  write(std_out,*)'eektf',eektf
343 !  stop ('vtorhotf')
344 !  ENDDEBUG
345    ek=eektf
346    entropy=(eektf-feektf)/dtset%tphysel
347    ABI_DEALLOCATE(betamumoinsV)
348  end subroutine tfek
```