TABLE OF CONTENTS
ABINIT/fermisolverec [ Functions ]
NAME
fermisolverec
FUNCTION
This routine computes the fermi energy in order to have a given number of valence electrons in the recursion method, using a Ridder s Method
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group ( ). 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
debug_rec=debugging variable nb_rec=order of recursion nb_point=number of discretization point in one dimension (=n1=n2=n3) temperature=temperature (Hartree) trotter=trotter parameter nelect=number of valence electrons (dtset%nelect) acc=accuracy for the fermi energy max_it=maximum number of iteration for the Ridder's Method long_tranche=number of point computed by thi proc mpi_enreg=information about MPI parallelization inf_ucvol=infinitesimal unit cell volume gputopo=true if topology gpu-cpu= 2 or 3
OUTPUT
SIDE EFFECTS
fermie=fermi energy rho=density, recomputed for the new fermi energy a, b2 : coefficient given by recursion recomputed for the new fermi energy
PARENTS
vtorhorec
CHILDREN
alloc_dens_cuda,dealloc_dens_cuda,density_cuda,density_rec,timab,wrtout xmpi_barrier,xmpi_sum
NOTES
at this time :
SOURCE
50 #if defined HAVE_CONFIG_H 51 #include "config.h" 52 #endif 53 54 #include "abi_common.h" 55 56 57 subroutine fermisolverec(fermie,rho,a,b2,debug_rec,nb_rec, & 58 & temperature,trotter,nelect, & 59 & acc, max_it, & 60 & long_tranche,mpi_enreg,& 61 & inf_ucvol,gputopo) 62 63 use defs_basis 64 use defs_abitypes 65 use defs_rectypes 66 use m_xmpi 67 use m_errors 68 use m_profiling_abi 69 70 #ifdef HAVE_GPU_CUDA 71 use m_initcuda,only : cudap 72 #endif 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'fermisolverec' 78 use interfaces_14_hidewrite 79 use interfaces_18_timing 80 use interfaces_68_recursion, except_this_one => fermisolverec 81 !End of the abilint section 82 83 implicit none 84 85 !Arguments ------------------------------- 86 !scalars 87 integer,intent(in) :: long_tranche,max_it,nb_rec,trotter 88 logical,intent(in) :: debug_rec,gputopo 89 real(dp),intent(in) :: acc,inf_ucvol,nelect,temperature 90 real(dp), intent(inout) :: fermie 91 type(MPI_type),intent(in) :: mpi_enreg 92 !arrays 93 real(dp), intent(inout) :: a(0:nb_rec,long_tranche), b2(0:nb_rec,long_tranche) 94 real(dp), intent(inout) :: rho(long_tranche) 95 96 !Local variables------------------------------- 97 !scalars 98 integer :: ierr,ii,ipointlocal,nn,dim_trott 99 real(dp) :: beta,fermieh,fermiel,fermiem,fermienew,nelecth,nelectl,nelectm 100 real(dp) :: nelectnew,res_nelecth,res_nelectl,res_nelectm,res_nelectnew 101 real(dp) :: rtrotter,ss,fermitol 102 character(len=500) :: msg 103 !arrays 104 real(dp) :: tsec(2) 105 real(dp) :: rhotry(long_tranche) 106 !no_abirules 107 #ifdef HAVE_GPU_CUDA 108 integer :: swt_tm,npitch 109 real(cudap) :: rhocu(long_tranche) 110 real(dp) :: tsec2(2) 111 #endif 112 113 ! ************************************************************************* 114 115 #ifdef HAVE_GPU_CUDA 116 swt_tm = 0 117 #endif 118 119 call timab(609,1,tsec) 120 121 beta = one/temperature 122 rtrotter = max(half,real(trotter,dp)) 123 dim_trott = max(0,2*trotter-1) 124 125 write(msg,'(a)')' -- fermisolverec ---------------------------------' 126 call wrtout(std_out,msg,'COLL') 127 if(debug_rec) then 128 write (msg,'(a,d10.3)')' nelect= ',nelect 129 call wrtout(std_out,msg,'COLL') 130 end if 131 !initialisation of fermiel 132 fermiel = fermie 133 call timab(609,2,tsec) 134 135 !initialisation fermitol 136 fermitol = acc 137 #ifdef HAVE_GPU_CUDA_SP 138 if(gputopo) fermitol = 1.d-3 139 #endif 140 141 if(gputopo) then 142 #ifdef HAVE_GPU_CUDA 143 swt_tm = 1 144 ! allocate array an and bn2 on gpu for computation of trotter formula 145 call alloc_dens_cuda(long_tranche,nb_rec,dim_trott,npitch,& 146 & real(a,cudap),real(b2,cudap)) 147 148 call timab(617,1,tsec) 149 call density_cuda(npitch,long_tranche,nb_rec,dim_trott,& 150 & real(fermiel,cudap),real(temperature,cudap),& 151 & real(rtrotter,cudap),real(inf_ucvol,cudap),& 152 & real(tol14,cudap),& 153 & rhocu) 154 rhotry = real(rhocu,dp) 155 call timab(617,2,tsec) 156 #endif 157 else 158 do ipointlocal = 1,long_tranche 159 call density_rec(a(:,ipointlocal),& 160 & b2(:,ipointlocal),& 161 & rhotry(ipointlocal),& 162 & nb_rec,fermiel,temperature,rtrotter,dim_trott, & 163 & tol14,inf_ucvol) 164 end do 165 end if 166 167 call timab(609,1,tsec) 168 nelectl = sum(rhotry) 169 call xmpi_sum( nelectl,mpi_enreg%comm_bandfft,ierr) 170 res_nelectl = inf_ucvol*nelectl - nelect 171 172 if (res_nelectl /= zero) then 173 ! initialisation of fermih 174 ! excess of electrons -> smaller fermi 175 res_nelecth = zero 176 ii = 1 177 fermieh = fermie - ten*sign(one,res_nelectl)*temperature 178 do while(ii<6 .and. res_nelecth*res_nelectl>=0) 179 fermieh = fermieh - ten*sign(one,res_nelectl)*temperature 180 call timab(609,2,tsec) 181 182 if(gputopo) then 183 #ifdef HAVE_GPU_CUDA 184 call timab(617,1,tsec) 185 call density_cuda(npitch,long_tranche,nb_rec,dim_trott,& 186 & real(fermieh,cudap),real(temperature,cudap),& 187 & real(rtrotter,cudap),real(inf_ucvol,cudap),& 188 & real(tol14,cudap),& 189 & rhocu) 190 rhotry = real(rhocu,dp) 191 call timab(617,2,tsec) 192 #endif 193 else 194 do ipointlocal = 1,long_tranche 195 call density_rec(a(:,ipointlocal), & 196 & b2(:,ipointlocal), & 197 & rhotry(ipointlocal), & 198 & nb_rec,fermieh,temperature,rtrotter,dim_trott, & 199 & tol14,inf_ucvol) 200 end do 201 end if 202 call timab(609,1,tsec) 203 nelecth = sum(rhotry) 204 call xmpi_sum( nelecth,mpi_enreg%comm_bandfft ,ierr); 205 res_nelecth = inf_ucvol*nelecth - nelect 206 207 if(debug_rec) then 208 write (msg,'(a,es11.4e2,a,es11.4e2)') ' Fermi energy interval',fermieh,' ',fermiel 209 call wrtout(std_out,msg,'COLL') 210 end if 211 ii = ii +1 212 end do 213 214 if (res_nelecth*res_nelectl>0) then 215 write (msg,'(4a)')' fermisolverec : ERROR- ',ch10,& 216 & ' initial guess for fermi energy doesnt permit to find solutions in solver',ch10 217 MSG_ERROR(msg) 218 end if 219 220 ! MAIN LOOP ------------------------------------------------------ 221 main : do nn=1,max_it 222 ! fermiem computation 223 fermiem = 0.5d0*(fermiel+fermieh) 224 225 ! nelectm = zero 226 call timab(609,2,tsec) 227 228 if(gputopo) then 229 #ifdef HAVE_GPU_CUDA 230 call timab(617,1,tsec) 231 call density_cuda(npitch,long_tranche,nb_rec,dim_trott,& 232 & real(fermiem,cudap),real(temperature,cudap),& 233 & real(rtrotter,cudap),real(inf_ucvol,cudap),& 234 & real(tol14,cudap),& 235 & rhocu) 236 rhotry = real(rhocu,dp) 237 call timab(617,2,tsec) 238 #endif 239 else 240 do ipointlocal = 1,long_tranche 241 call density_rec(a(:,ipointlocal), & 242 & b2(:,ipointlocal), & 243 & rhotry(ipointlocal), & 244 & nb_rec,fermiem,temperature,rtrotter,dim_trott, & 245 & tol14,inf_ucvol) 246 end do 247 end if 248 249 call timab(609,1,tsec) 250 nelectm = sum(rhotry) 251 call xmpi_sum( nelectm,mpi_enreg%comm_bandfft,ierr) 252 res_nelectm = inf_ucvol*nelectm - nelect 253 254 ! new guess 255 ss = sqrt(res_nelectm**two-res_nelectl*res_nelecth) 256 fermienew = fermiem + (fermiem-fermiel)*sign(one, res_nelectl-res_nelecth)*res_nelectm/ss 257 258 call timab(609,2,tsec) 259 if(gputopo) then 260 #ifdef HAVE_GPU_CUDA 261 call timab(617,1,tsec) 262 call density_cuda(npitch,long_tranche,nb_rec,dim_trott,& 263 & real(fermienew,cudap),real(temperature,cudap),& 264 & real(rtrotter,cudap),real(inf_ucvol,cudap),& 265 & real(tol14,cudap),& 266 & rhocu) 267 rhotry = real(rhocu,dp) 268 call timab(617,2,tsec) 269 #endif 270 else 271 do ipointlocal = 1,long_tranche 272 call density_rec(a(:,ipointlocal), & 273 & b2(:,ipointlocal), & 274 & rhotry(ipointlocal), & 275 & nb_rec,fermienew,temperature,rtrotter,dim_trott, & 276 & tol14,inf_ucvol) 277 end do 278 end if 279 280 call timab(609,1,tsec) 281 nelectnew = sum(rhotry) 282 call xmpi_sum( nelectnew,mpi_enreg%comm_bandfft ,ierr); 283 res_nelectnew = inf_ucvol*nelectnew - nelect 284 285 ! fermiel et fermieh for new iteration 286 if (sign(res_nelectm,res_nelectnew) /= res_nelectm) then 287 fermiel = fermiem 288 res_nelectl = res_nelectm 289 fermieh = fermienew 290 res_nelecth = res_nelectnew 291 else if (sign(res_nelectl,res_nelectnew) /= res_nelectl) then 292 fermieh = fermienew 293 res_nelecth = res_nelectnew 294 else if (sign(res_nelecth,res_nelectnew) /= res_nelecth) then 295 fermiel = fermienew 296 res_nelectl = res_nelectnew 297 end if 298 299 ! are we within the tolerance ? 300 if ((abs(res_nelectnew) < fermitol).or.(nn == max_it)) then 301 fermie = fermienew 302 rho = rhotry 303 if(debug_rec) then 304 write (msg,'(a,es11.4e2,a,i4)')' err, num_iter ', res_nelectnew, ' ',nn 305 call wrtout(std_out,msg,'COLL') 306 write(msg,'(a,50a)')' ',('-',ii=1,50) 307 call wrtout(std_out,msg,'COLL') 308 end if 309 exit main 310 end if 311 312 end do main 313 314 end if 315 316 #ifdef HAVE_GPU_CUDA 317 !deallocate array on GPU 318 if(gputopo) then 319 call dealloc_dens_cuda() 320 end if 321 call timab(613+swt_tm,1,tsec2) !!--start time-counter: sync gpu-cpu 322 call xmpi_barrier(mpi_enreg%comm_bandfft) 323 call timab(613+swt_tm,2,tsec2) !!--stop time-counter: sync gpu-cpu 324 #endif 325 326 327 call timab(609,2,tsec) 328 end subroutine fermisolverec