TABLE OF CONTENTS

ABINIT/fm12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

 fm12a1

FUNCTION

INPUTS

OUTPUT

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

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

ABINIT/fm12a1t [ Functions ]

[ Top ] [ Functions ]

NAME

 fm12a1t

FUNCTION

INPUTS

OUTPUT

PARENTS

      vtorhotf

CHILDREN

SOURCE

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

ABINIT/fp12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

 fp12a1

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

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

ABINIT/fp32a1 [ Functions ]

[ Top ] [ Functions ]

NAME

 fp32a1

FUNCTION

INPUTS

OUTPUT

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

1389  function fp32a1 (x)
1390 
1391    use defs_basis
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

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

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

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

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

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

ABINIT/vtorhotf [ Functions ]

[ Top ] [ Functions ]

NAME

 vtorhotf

FUNCTION

 This routine 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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

 52 #if defined HAVE_CONFIG_H
 53 #include "config.h"
 54 #endif
 55 
 56 #include "abi_common.h"
 57 
 58 subroutine vtorhotf(dtfil,dtset,ek,enl,entropy,fermie,gprimd,grnl,&
 59 &  irrzon,mpi_enreg,natom,nfft,nspden,nsppol,nsym,phnons,rhog,rhor,rprimd,ucvol,vtrial)
 60 
 61  use defs_basis
 62  use defs_abitypes
 63  use m_profiling_abi
 64  use m_errors
 65  use m_xmpi
 66 
 67 !This section has been created automatically by the script Abilint (TD).
 68 !Do not modify the following lines by hand.
 69 #undef ABI_FUNC
 70 #define ABI_FUNC 'vtorhotf'
 71  use interfaces_14_hidewrite
 72  use interfaces_18_timing
 73  use interfaces_32_util
 74 !End of the abilint section
 75 
 76  implicit none
 77 
 78 !Arguments ------------------------------------
 79 !scalars
 80  integer,intent(in) :: natom,nfft,nspden,nsppol,nsym
 81  real(dp),intent(in) :: ucvol
 82  real(dp),intent(out) :: ek,enl,entropy,fermie
 83  type(MPI_type),intent(in) :: mpi_enreg
 84  type(datafiles_type),intent(in) :: dtfil
 85  type(dataset_type),intent(in) :: dtset
 86 !arrays
 87  integer,intent(in) :: irrzon((dtset%ngfft(1)*dtset%ngfft(1)*dtset%ngfft(1))**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))
 88  real(dp),intent(in) :: gprimd(3,3)
 89  real(dp),intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(1)*dtset%ngfft(1))**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))
 90  real(dp),intent(in) :: rprimd(3,3),vtrial(nfft,nspden)
 91  real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,nspden)
 92  real(dp),intent(out) :: grnl(3*natom)
 93 
 94 !Local variables-------------------------------
 95 !scalars
 96  integer,parameter :: jdichomax=20,level=111
 97  integer :: i1,i2,i3,ierr,iexit,ifft,ii,ir,iscf,jdicho
 98  integer :: me_fft,n1,n2,n3,nfftot,nproc_fft,prtvol
 99  real(dp),save :: cktf,fermie_tol,nelect_mid
100  real(dp) :: dnelect_mid_dx,dxrtnewt,eektemp,eektf,feektemp,feektf
101  real(dp) :: rtnewt,sum_rhor_mid,sum_rhor_middx
102  logical,save :: lfirst_time_tf=.true.
103  logical :: lnewtonraphson
104  character(len=500) :: message
105 !arrays
106  real(dp) :: tsec(2)
107  real(dp),allocatable :: betamumoinsV(:),rhor_mid(:),rhor_middx(:)
108 !no_abirules
109 
110 ! *************************************************************************
111 
112 !Keep track of total time spent in vtorho
113  call timab(21,1,tsec)
114 
115 
116  call status(0,dtfil%filstat,iexit,level,'enter         ')
117 
118 !Structured debugging if prtvol==-level
119  prtvol=dtset%prtvol
120  if(prtvol==-level)then
121    write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' vtorho : enter '
122    call wrtout(std_out,message,'COLL')
123  end if
124 
125  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
126  me_fft=dtset%ngfft(11) ; nproc_fft=dtset%ngfft(10)
127  iscf=dtset%iscf
128 !Debugging : print vtrial and rhor
129  if(prtvol==-level)then
130    write(message,'(a)') '   ir              vtrial(ir)     rhor(ir) '
131    call wrtout(std_out,message,'COLL')
132    do ir=1,nfft
133 !    if(ir<=11 .or. mod(ir,301)==0 )then
134      i3=(ir-1)/n1/(n2/nproc_fft)
135      i2=(ir-1-i3*n1*n2/nproc_fft)/n1
136      i1=ir-1-i3*n1*n2/nproc_fft-(i2-me_fft)*n1
137      write(message,'(i5,3i3,a,2d13.6)')ir,i1,i2,i3,' ',vtrial(ir,1),rhor(ir,1)
138      call wrtout(std_out,message,'COLL')
139      if(nspden>=2)then
140        write(message,'(a,2d13.6)')'               ',vtrial(ir,2),rhor(ir,2)
141        call wrtout(std_out,message,'COLL')
142      end if
143 !    end if
144    end do
145  end if
146 
147  ek=zero
148  enl=zero
149  grnl(:)=zero
150 
151 !Initialize rhor if needed
152  if(iscf>0) rhor(:,:)=zero
153 
154 !call Thomas-Fermi for the density
155  call tf
156 !Compute energy terms
157  call tfek
158 
159  call status(0,dtfil%filstat,iexit,level,'exit          ')
160 
161  call timab(21,2,tsec)
162 !End thomas fermi
163 
164  contains

ABINIT/xp12a1 [ Functions ]

[ Top ] [ Functions ]

NAME

 xp12a1

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

1448  function xp12a1 (y)
1449 
1450    use defs_basis
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

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

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

ABINIT/zfermi12 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermi12

FUNCTION

INPUTS

OUTPUT

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

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

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

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

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

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

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

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

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

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

ABINIT/zfermim12 [ Functions ]

[ Top ] [ Functions ]

NAME

 zfermim12

FUNCTION



INPUTS


OUTPUT


PARENTS


  Will be filled automatically by the parent script



CHILDREN


  Will be filled automatically by the parent script



SOURCE


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







vtorhotf/tf [ Functions ]



[ Top ] [ vtorhotf ] [ Functions ]


NAME


 tf



FUNCTION


INPUTS


OUTPUT


PARENTS


      vtorhotf



CHILDREN


SOURCE


183   subroutine tf()
184 
185 
186 !This section has been created automatically by the script Abilint (TD).
187 !Do not modify the following lines by hand.
188 #undef ABI_FUNC
189 #define ABI_FUNC 'tf'
190  use interfaces_18_timing
191  use interfaces_32_util
192  use interfaces_67_common
193 !End of the abilint section
194 
195   implicit none
196 
197 ! *************************************************************************
198 
199    ABI_ALLOCATE(rhor_mid,(nfft))
200    ABI_ALLOCATE(rhor_middx,(nfft))
201    fermie_tol=1.e-10_dp
202    cktf=one/two/pi**2*(two*dtset%tphysel)**1.5_dp
203 
204 !  Should be made an input variable, if TF really needed for production
205 !  rtnewt=dtset%userra
206    rtnewt=zero
207 
208 !  Newton Raphson
209    if (lfirst_time_tf) then
210      lfirst_time_tf=.false.
211    end if
212    jdicho=0
213    lnewtonraphson=.false.
214    do while (.not.lnewtonraphson)
215      jdicho=jdicho+1
216 !    do ifft=1,nfft
217 !    rhor_mid(ifft)=cktf*zfermi12((rtnewt-vtrial(ifft,1))/dtset%tphysel)
218 !    rhor_middx(ifft)=cktf*zfermim12((rtnewt-vtrial(ifft,1))/dtset%tphysel)
219 !    end do
220      call fm12a1t(cktf,rtnewt,dtset%tphysel,vtrial(:,1),rhor_middx,rhor_mid,&
221 &     nfft)
222      sum_rhor_mid=sum(rhor_mid(:))
223      sum_rhor_middx=sum(rhor_middx(:))
224      call xmpi_sum(sum_rhor_mid,mpi_enreg%comm_fft ,ierr)
225      call xmpi_sum(sum_rhor_middx,mpi_enreg%comm_fft ,ierr)
226      nelect_mid=sum_rhor_mid*ucvol/(nfft*nproc_fft)-dtset%nelect
227      dnelect_mid_dx=sum_rhor_middx*ucvol/(nfft*nproc_fft)/dtset%tphysel/two
228      dxrtnewt=nelect_mid/dnelect_mid_dx
229      rtnewt=rtnewt-dxrtnewt
230      if (abs(nelect_mid) < fermie_tol/2._dp) then
231        lnewtonraphson=.true.
232      end if
233      if (jdicho > jdichomax) then
234        MSG_ERROR('NEWTON RAPHSON NOT CONVERGED')
235      end if
236    end do
237    fermie=rtnewt
238    rhor(:,1)=rhor_mid(:)
239    ABI_DEALLOCATE(rhor_mid)
240    ABI_DEALLOCATE(rhor_middx)
241 
242 !  DEBUG
243 !  write(std_out,*)'fmid,nmid,jdicho',fermie,nelect_mid,jdicho
244 !  ENDDEBUG
245 
246 !  Compute rhog
247    call timab(70,1,tsec)
248 
249    call status(0,dtfil%filstat,iexit,level,'compute rhog  ')
250    nfftot=dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)
251    call symrhg(1,gprimd,irrzon,mpi_enreg,nfft,nfftot,dtset%ngfft,nspden,nsppol,nsym,dtset%paral_kgb,phnons,&
252 &   rhog,rhor,rprimd,dtset%symafm,dtset%symrel)
253 
254 !  We now have both rho(r) and rho(G), symmetrized, and if nsppol=2
255 !  we also have the spin-up density, symmetrized, in rhor(:,2).
256 
257    call timab(70,2,tsec)
258  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


280   subroutine tfek()
281 
282 
283 !This section has been created automatically by the script Abilint (TD).
284 !Do not modify the following lines by hand.
285 #undef ABI_FUNC
286 #define ABI_FUNC 'tfek'
287  use interfaces_18_timing
288  use interfaces_67_common
289 !End of the abilint section
290 
291   implicit none
292 
293 ! *************************************************************************
294 
295    ABI_ALLOCATE(betamumoinsV,(nfft))
296    cktf=one/two/pi**2*(two*dtset%tphysel)**1.5_dp
297    eektf=zero
298    feektf=zero
299    do ifft=1,nfft
300 
301 !    betamumoinsv(ifft)=ifermi12(rhor(ifft,1)/cktf)
302      betamumoinsv(ifft)=(rtnewt-vtrial(ifft,1))/dtset%tphysel
303 !    eektemp=zfermi32(betamumoinsV(ifft))/zfermi12(betamumoinsV(ifft))
304      eektemp=fp32a1(betamumoinsV(ifft))/rhor(ifft,1)*cktf
305      feektemp=betamumoinsV(ifft)-two/three*eektemp
306      feektf=feektf+feektemp*rhor(ifft,1)
307      eektf=eektf+eektemp*rhor(ifft,1)
308    end do
309 !  Init mpi_comm
310    call timab(48,1,tsec)
311    call xmpi_sum(eektf,mpi_enreg%comm_fft ,ierr)
312    call xmpi_sum(feektf,mpi_enreg%comm_fft ,ierr)
313    call timab(48,2,tsec)
314    eektf=eektf*dtset%tphysel
315    eektf=eektf*ucvol/dble(nfft*nproc_fft)
316    feektf=feektf*dtset%tphysel
317    feektf=feektf*ucvol/dble(nfft*nproc_fft)
318 !  DEBUG
319 !  write(std_out,*)'eektf',eektf
320 !  stop ('vtorhotf')
321 !  ENDDEBUG
322    ek=eektf
323    entropy=(eektf-feektf)/dtset%tphysel
324    ABI_DEALLOCATE(betamumoinsV)
325  end subroutine tfek
326 
327 end subroutine vtorhotf