TABLE OF CONTENTS
ABINIT/recursion [ Functions ]
NAME
recursion
FUNCTION
This routine computes the recursion coefficients and the corresponding continued fraction to get the density at a point from a fixed potential.
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group (SLeroux,MMancini). 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
exppot=exponential of -1/tsmear*vtrial (computed only once in vtorhorec) coordx, coordy, coordz=coordonnees of the computed point nrec=order of recursion fermie=fermi energy (Hartree) tsmear=temperature (Hartree) dim_trott=dimension of the partial fraction decomposition rtrotter=trotter parameter (real) ZT_p=fourier transform of the Green krenel (computed only once in vtorhorec) typat(:)=type of psp associated to any atom tol=tolerance criteria for stopping recursion debug=debugging variable mpi_enreg=information about MPI paralelisation nfft=number of points in FFT grid ngfft=information about FFT metrec<type(metricrec_type)>=information concerning the infinitesimal metrics inf_ucvol=infinitesimal unit cell volume tim_fourdp=time counter for fourdp natom=number of atoms projec(ngfftrec(1),ngfftrec(2),ngfftrec(3),lmnmax,natom) is the vector, on the ngfftrec grid containing the non-lacal projector $Y_{lm}(r-R_A)f_{lk}(r-R_A) tim= 0 if the time spent in the routine is not taken into account,1 otherwise. For example when measuring time for loading balancing, we don't want to add the time spent in this to the total time calculation
OUTPUT
rho_out=result of the continued fraction multiplied by a multiplicator an, bn2 : coefficient given by recursion.
SIDE EFFECTS
PARENTS
first_rec,vtorhorec
CHILDREN
fourdp,timab,trottersum,vn_nl_rec
NOTES
at this time :
- exppot should be replaced by ?
- coord should be replaced by ?
- need a rectangular box (rmet diagonal matrix)
SOURCE
62 #if defined HAVE_CONFIG_H 63 #include "config.h" 64 #endif 65 66 #include "abi_common.h" 67 68 69 subroutine recursion(exppot,coordx,coordy,coordz,an,bn2,rho_out, & 70 & nrec,fermie,tsmear,rtrotter,dim_trott, & 71 & ZT_p, tol,typat, & 72 & nlrec,mpi_enreg,& 73 & nfft,ngfft,metrec,& 74 & tim_fourdp,natom,projec,tim) 75 76 use m_profiling_abi 77 78 use defs_basis 79 use defs_abitypes 80 use defs_rectypes 81 use m_rec_tools,only : trottersum 82 use m_linalg_interfaces 83 84 !This section has been created automatically by the script Abilint (TD). 85 !Do not modify the following lines by hand. 86 #undef ABI_FUNC 87 #define ABI_FUNC 'recursion' 88 use interfaces_18_timing 89 use interfaces_53_ffts 90 use interfaces_68_recursion, except_this_one => recursion 91 !End of the abilint section 92 93 implicit none 94 95 !Arguments ------------------------------- 96 !scalars 97 integer,intent(in) :: coordx,coordy,coordz,nfft,nrec,tim 98 integer,intent(in) :: tim_fourdp,natom,dim_trott 99 real(dp),intent(in) :: fermie,tol,tsmear,rtrotter 100 real(dp), intent(out) :: rho_out 101 type(MPI_type),intent(in) :: mpi_enreg 102 type(nlpsprec_type),intent(in) :: nlrec 103 type(metricrec_type),intent(in) :: metrec 104 !arrays 105 integer, intent(in) :: ngfft(18) 106 integer, intent(in) :: typat(natom) 107 real(dp), intent(in) :: ZT_p(1:2, 0:nfft-1) 108 real(dp), intent(in) :: exppot(0:nfft-1) 109 real(dp), intent(in) :: projec(0:,0:,0:,1:,1:) 110 real(dp), intent(out) :: an(0:nrec),bn2(0:nrec) 111 !Local variables------------------------------- 112 !not used, debugging purpose only 113 !for debugging purpose, detailled printing only once for density and ekin 114 !scalars 115 integer, parameter :: level = 7, minrec = 3 116 integer :: irec,isign,timab_id,ii 117 real(dp) :: switchimu,switchu 118 real(dp) :: bb,beta,mult,prod_b2,error,errold 119 real(dp) :: inf_ucvol,pi_on_rtrotter,twortrotter,exp1,exp2 120 complex(dpc) :: cinv2rtrotter,coeef_mu,facrec0 121 ! character(len=500) :: msg 122 !arrays 123 real(dp) :: tsec(2) 124 real(dp) :: inf_tr(3) 125 real(dp) :: Zvtempo(1:2, 0:nfft-1) 126 real(dp) :: unold(0:nfft-1),vn(0:nfft-1),un(0:nfft-1) 127 complex(dpc) :: acc_rho(0:nrec) 128 complex(dpc) :: D(0:dim_trott),Dold(0:dim_trott) 129 complex(dpc) :: N(0:dim_trott),Nold(0:dim_trott) 130 131 ! ************************************************************************* 132 133 !--If count time or not 134 timab_id = 616; if(tim/=0) timab_id = 606; 135 136 call timab(timab_id,1,tsec) 137 138 !############################################################## 139 !--Initialisation of metrics 140 inf_ucvol = metrec%ucvol 141 inf_tr = metrec%tr 142 mult = two/inf_ucvol !non-spined system 143 144 beta = one/tsmear 145 !--Variables for optimisation 146 pi_on_rtrotter = pi/rtrotter 147 twortrotter = two*rtrotter 148 exp1 = exp((beta*fermie)/(rtrotter)) 149 exp2 = exp(beta*fermie/(twortrotter)) 150 cinv2rtrotter = cmplx(one/twortrotter,zero,dp) 151 coeef_mu = cmplx(one/exp2,zero,dp) 152 153 !--Initialisation of an,bn,un.... 154 N = czero; D = cone 155 facrec0 = cone 156 Nold = czero; Dold = czero 157 158 an = zero; bn2 = zero; bn2(0) = one 159 bb = zero; vn = zero; unold = zero 160 !--u0 is a Dirac function 161 un = zero 162 un(coordx+ngfft(1)*(coordy+ngfft(2)*coordz)) = one/sqrt(inf_ucvol) 163 164 !--Initialisation of accumulated density 165 acc_rho = czero 166 !--Initialisation of estimated error 167 prod_b2 = twortrotter/exp1 168 errold = zero 169 170 !############################################################## 171 !--Main loop 172 maindo : do irec = 0, nrec 173 174 ! --Get an and bn2 coef by the lanczos method 175 176 ! --Computation of exp(-beta*V/8*p)*un or exp(-beta*V/4*p)*un 177 ! depending on if nl part has to be calculated or not. 178 vn = exppot * un 179 180 ! --First Non-local psp contribution: (Id+sum_atom int dr1(E(r,r1))vn(r1)) 181 ! --Computation of exp(-beta*V_NL/4*p)*vn 182 if(nlrec%nlpsp) then 183 call timab(timab_id,2,tsec) 184 call vn_nl_rec(vn,natom,typat,ngfft(:3),inf_ucvol,nlrec,projec) 185 call timab(timab_id,1,tsec) 186 187 ! --Computation of exp(-beta*V/8*p)*vn in nonlocal case 188 vn = exppot * vn 189 end if !--End if on nlrec%nlpsp 190 191 ! --Convolution with the Green kernel 192 ! --FFT of vn 193 isign = -1 194 call fourdp(1,Zvtempo,vn,isign,mpi_enreg,nfft,ngfft,1,tim_fourdp) 195 196 ! --F(T)F(vn) 197 do ii = 0,nfft-1 198 switchu = Zvtempo(1,ii) 199 switchimu = Zvtempo(2,ii) 200 Zvtempo(1,ii) = switchu*ZT_p(1,ii) - switchimu*ZT_p(2,ii) 201 Zvtempo(2,ii) = switchu*ZT_p(2,ii) + switchimu*ZT_p(1,ii) 202 end do 203 204 ! --F^-1(F(T)F(vn)) 205 isign = 1 206 call fourdp(1,Zvtempo,vn,isign,mpi_enreg,nfft,ngfft,1,tim_fourdp) 207 208 ! --Computation of exp(-beta*V/8*p)*un or exp(-beta*V/4*p)*un 209 ! depending on if nl part has to be calculated or not. 210 211 vn = inf_ucvol * exppot * vn 212 213 ! --Second Non-local psp contribution: (Id+sum_atom E(r,r1))vn 214 if(nlrec%nlpsp) then 215 call timab(timab_id,2,tsec) 216 call vn_nl_rec(vn,natom,typat,ngfft(:3),inf_ucvol,nlrec,projec) 217 call timab(timab_id,1,tsec) 218 219 ! --Computation of exp(-beta*V/8*p)*vn in nonlocal case 220 vn = exppot * vn 221 end if !--End if on nlrec%nlpsp 222 223 ! --Multiplication of a and b2 coef by exp(beta*fermie/(two*rtrotter)) must be done in the continued fraction computation 224 ! --Computation of a and b2 225 an(irec) = inf_ucvol*ddot(nfft,vn,1,un,1) 226 227 ! --an must be positive real 228 ! --We must compute bn2 and prepare for the next iteration 229 if(irec<nrec)then 230 do ii = 0,nfft-1 231 switchu = un(ii) 232 un(ii) = vn(ii)-an(irec)*un(ii)-bb*unold(ii) 233 unold(ii) = switchu 234 bn2(irec+1) = bn2(irec+1)+inf_ucvol*un(ii)*un(ii) 235 end do 236 bb = sqrt(bn2(irec+1)) 237 un = (one/bb)*un 238 end if 239 240 ! ###################################################### 241 ! --Density computation 242 ! density computation is done inside the main looping, juste after the calculus of a and b2, in order to make 243 ! it possible to stop the recursion at the needed accuracy, without doing more recursion loop than needed - 244 ! further developpement 245 246 ! !--using the property that: sum_i(bi*c)^2|(z-ai*c)=1/c*sum_i(bi)^2|(z/c-ai) 247 ! !and for c =exp(-beta*fermie/(two*rtrotter) 248 249 250 call trottersum(dim_trott,error,prod_b2,pi_on_rtrotter,& 251 & facrec0,coeef_mu,exp1,& 252 & an(irec),bn2(irec),& 253 & N,D,Nold,Dold) 254 255 256 if(irec/=nrec .and. irec>=minrec)then 257 if((bn2(irec+1)<tol14).or.(mult*error<tol.and.errold<tol)) exit 258 end if 259 errold = mult*error 260 end do maindo 261 !--Accumulated density 262 rho_out = mult*real(cone-sum(N/D)*cinv2rtrotter,dp) 263 264 265 call timab(timab_id,2,tsec) 266 267 end subroutine recursion