TABLE OF CONTENTS
ABINIT/recursion_nl [ Functions ]
NAME
recursion_nl
FUNCTION
Given a $|un>$ vector on the real-space grid this routine calculates the density in $|un>$ by recursion 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
exppot=exponential of -1/tsmear*vtrial (computed only once in vtorhorec) trotter=trotter parameter dim_trott=dimension of the partial fraction decomposition tsmear=temperature (Hartree) tol=tolerance criteria for stopping recursion_nl ngfft=information about FFT(dtset%ngfft a priori different from ngfftrec) rset<recursion_type> contains all parameter of recursion typat(natom)=type of pseudo potential associated to any atom 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)
OUTPUT
rho_out=result of the continued fraction multiplied by a multiplicator
SIDE EFFECTS
un(:,:,:)=initial vector on the grid. it is changed in output
PARENTS
nlenergyrec
CHILDREN
fourdp,timab,trottersum,vn_nl_rec,wrtout
NOTES
at this time :
- need a rectangular box (rmet diagonal matrix)
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 55 subroutine recursion_nl(exppot,un,rho_out,rset,ngfft, & 56 & tsmear,trotter,dim_trott,tol,typat,& 57 & natom,projec) 58 59 use m_profiling_abi 60 61 use defs_basis 62 use defs_abitypes 63 use defs_rectypes 64 use m_rec_tools,only : trottersum 65 use m_linalg_interfaces 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 'recursion_nl' 71 use interfaces_14_hidewrite 72 use interfaces_18_timing 73 use interfaces_53_ffts 74 use interfaces_68_recursion, except_this_one => recursion_nl 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------- 80 !scalars 81 integer,intent(in) :: trotter,natom,dim_trott 82 real(dp),intent(in) :: tol,tsmear 83 type(recursion_type),intent(in) :: rset 84 real(dp), intent(out) :: rho_out 85 !arrays 86 integer,intent(in) :: typat(natom),ngfft(18) 87 real(dp),intent(in) :: exppot(0:ngfft(1)*ngfft(2)*ngfft(3)-1) 88 real(dp),intent(inout) :: un(0:rset%nfftrec-1) 89 real(dp),pointer :: projec(:,:,:,:,:) 90 !Local variables------------------------------- 91 !scalars 92 integer, parameter :: minrec = 3 93 integer :: irec,isign,ii 94 real(dp) :: bb,beta,mult,prod_b2,rtrotter 95 real(dp) :: inf_ucvol,pi_on_rtrotter,twortrotter,exp1 96 real(dp) :: exp2,error,errold 97 real(dp) :: switchu,switchimu 98 complex(dpc) :: facrec0,cinv2rtrotter,coeef_mu 99 character(len=500) :: msg 100 type(mpi_type),pointer:: mpi_loc 101 !arrays 102 real(dp):: tsec(2) 103 real(dp):: inf_tr(3) 104 real(dp):: an(0:rset%min_nrec),bn2(0:rset%min_nrec) 105 real(dp):: vn(0:rset%nfftrec-1) 106 real(dp):: unold(0:rset%nfftrec-1) 107 real(dp):: Zvtempo(1:2,0:rset%nfftrec-1) 108 complex(dpc) :: acc_rho(0:rset%min_nrec) 109 complex(dpc) :: D(0:dim_trott),Dold(0:dim_trott) 110 complex(dpc) :: N(0:dim_trott),Nold(0:dim_trott) 111 112 ! ************************************************************************* 113 114 call timab(608,1,tsec) !--start time-counter: recursion_nl 115 if(rset%debug)then 116 msg=' ' 117 call wrtout(std_out,msg,'COLL') 118 end if 119 120 !############################################################## 121 beta = one/tsmear 122 123 !--Rewriting the trotter parameter 124 rtrotter = max(half,real(trotter,dp)) 125 126 !--Initialisation of mpi 127 mpi_loc => rset%mpi 128 129 !--Initialisation of metrics 130 inf_ucvol = rset%inf%ucvol 131 inf_tr = rset%inf%tr 132 mult = one !--In the case of the calculus of the NL-energy 133 134 !--Initialisation of an,bn,un.... 135 N = czero; D = cone 136 facrec0 = cone 137 Nold = czero; Dold = czero 138 139 an = zero; bn2 = zero; bn2(0) = one 140 bb = zero; vn = zero; unold = zero 141 142 !--Variables for optimisation 143 pi_on_rtrotter = pi/rtrotter 144 twortrotter = two*rtrotter 145 exp1 = exp((beta*rset%efermi)/(rtrotter)) 146 exp2 = exp(beta*rset%efermi/(twortrotter)) 147 cinv2rtrotter = cmplx(one/twortrotter,zero,dp) 148 coeef_mu = cmplx(one/exp2,zero,dp) 149 150 !--Initialisation of accumulated density 151 acc_rho = czero 152 !--Initialisation of estimated error 153 prod_b2 = twortrotter/exp1 154 errold = zero 155 156 !############################################################## 157 !--Main loop 158 maindo : do irec = 0, rset%min_nrec 159 ! --Get an and bn2 coef by the lanczos method 160 161 ! --Computation of exp(-beta*V/8*p)*un 162 vn = exppot * un 163 164 ! --First Non-local psp contribution: (Id+sum_atom E(r,r1))vn 165 call timab(608,2,tsec) 166 call vn_nl_rec(vn,natom,typat,rset%ngfftrec(:3),inf_ucvol,rset%nl,projec) 167 call timab(608,1,tsec) 168 169 ! --Computation of exp(-beta*V/8*p)*un 170 vn = exppot * vn 171 172 ! --Convolution with the Green kernel 173 ! --FFT of vn 174 isign = -1 175 call fourdp(1,Zvtempo,vn,isign,mpi_loc,rset%nfftrec,rset%ngfftrec,1,6) 176 177 ! --F(T)F(vn) 178 do ii = 0,rset%nfftrec-1 179 switchu = Zvtempo(1,ii) 180 switchimu = Zvtempo(2,ii) 181 Zvtempo(1,ii) = switchu*rset%ZT_p(1,ii) - switchimu*rset%ZT_p(2,ii) 182 Zvtempo(2,ii) = switchu*rset%ZT_p(2,ii) + switchimu*rset%ZT_p(1,ii) 183 end do 184 185 ! --F^-1(F(T)F(vn)) 186 isign = 1 187 call fourdp(1,Zvtempo,vn,isign,mpi_loc,rset%nfftrec,rset%ngfftrec,1,6) 188 189 ! --Computation of exp(-beta*V/2*p)*vn 190 vn = inf_ucvol * exppot * vn 191 192 ! --Second Non-local psp contribution: (Id+sum_atom E(r,r1))vn 193 call timab(608,2,tsec) 194 call vn_nl_rec(vn,natom,typat,rset%ngfftrec(:3),inf_ucvol,rset%nl,projec) 195 call timab(608,1,tsec) 196 197 ! --Computation of exp(-beta*V/8*p)*vn 198 vn = exppot * vn 199 200 201 ! --Multiplication of a and b2 coef by exp(beta*fermie/(2.d0*rtrotter)) must be done in the continued fraction computation 202 ! --Computation of a and b2 203 an(irec) = inf_ucvol*ddot(rset%nfftrec,vn,1,un,1) !--an must be positive real 204 205 ! --We must compute bn2 and prepare for the next iteration 206 if(irec<rset%min_nrec)then 207 do ii = 0,rset%nfftrec-1 208 switchu = un(ii) 209 un(ii) = vn(ii)-an(irec)*un(ii)-bb*unold(ii) 210 unold(ii) = switchu 211 bn2(irec+1) = bn2(irec+1)+inf_ucvol*un(ii)*un(ii) 212 end do 213 bb = sqrt(bn2(irec+1)) 214 un = (one/bb)*un 215 end if 216 217 ! ###################################################### 218 ! --Density computation 219 ! in order to make it possible to stop the recursion_nl at the 220 ! needed accuracy, without doing more recursion_nl loop than needed further developpement 221 222 call trottersum(dim_trott,error,& 223 & prod_b2,pi_on_rtrotter,& 224 & facrec0,coeef_mu,exp1,& 225 & an(irec),bn2(irec),& 226 & N,D,Nold,Dold) 227 228 229 if(irec/=rset%min_nrec .and. irec>=minrec)then 230 if((bn2(irec+1)<tol14).or.(mult*error<tol.and.errold<tol)) exit 231 end if 232 errold = mult*error 233 end do maindo 234 !--Accumulated density 235 rho_out = mult*real(cone-sum(N/D)*cinv2rtrotter,dp) 236 237 call timab(608,2,tsec) !--stop time-counter: recursion_nl 238 239 end subroutine recursion_nl