TABLE OF CONTENTS


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

  Copyright (C) 1998-2018 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 .

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
  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.
  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