TABLE OF CONTENTS
ABINIT/linemin [ Functions ]
NAME
linemin
FUNCTION
Performs the "line minimization" w.r.t. the angle theta on a unit circle to update the wavefunction associated with the current k-point and band label. This routine is used only when the electric field is on (otherwise it could in principle also be used, but there is a simpler procedure, as originally coded in abinit).
COPYRIGHT
Copyright (C) 2000-2017 ABINIT group (MVeithen,ISouza,JIniguez) 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
chc = <C|H_0|C> where |C> is the wavefunction of the current band
detovc = determinant of the overlap matrix S
detovd = determinant of the overlap matrix where for the band
that is being updated <C| is replaced by <D| (search direction)
dhc = Re[<D|H_0|C>]
dhd = <D|H_0|D>
efield_dot = reciprocal lattice coordinates of the electric field
iline = index of the current line minimization
nkpt = number of k-points
nstr(idir) = number of strings along the idir-th direction
sdeg = spin degeneracy
OUTPUT
bcut(ifor,idir) = branch cut of the ellipse associated with (ifor,idir) costh = cos(thetam) hel(ifor,idir) = helicity of the ellipse associated with (ifor,idir) phase_end = total change in Zak phase, must be equal to dphase_aux1 + n*two_pi sinth = sin(thetam) thetam = optimal angle theta in line_minimization
SIDE EFFECTS
Input/Output dphase_aux1 = change in Zak phase accumulated during the loop over iline (can be used for debugging in cgwf.f) phase_init = initial Zak phase (before doing the first line minimization)
TODO
NOTES
We are making the "frozen Hamiltonian approximation", i.e., the Hamiltonian does not change with theta (we are neglecting the dependence of the Hartree and exchange-correlation terms on theta; the original abinit routine does the same)
PARENTS
cgwf
CHILDREN
etheta,rhophi,wrtout
SOURCE
67 #if defined HAVE_CONFIG_H 68 #include "config.h" 69 #endif 70 71 #include "abi_common.h" 72 73 74 subroutine linemin(bcut,chc,costh,detovc,detovd,dhc,dhd,dphase_aux1,& 75 & efield_dot,iline,nkpt,nstr,hel,phase_end,phase_init,sdeg,sinth,thetam) 76 77 use defs_basis 78 use m_errors 79 use m_profiling_abi 80 81 use m_numeric_tools, only : rhophi 82 83 !This section has been created automatically by the script Abilint (TD). 84 !Do not modify the following lines by hand. 85 #undef ABI_FUNC 86 #define ABI_FUNC 'linemin' 87 use interfaces_14_hidewrite 88 use interfaces_67_common, except_this_one => linemin 89 !End of the abilint section 90 91 implicit none 92 93 !Arguments ------------------------------------ 94 !scalars 95 integer,intent(in) :: iline,nkpt 96 real(dp),intent(in) :: chc,dhc,dhd,sdeg 97 real(dp),intent(out) :: costh,sinth,thetam 98 !arrays 99 integer,intent(in) :: nstr(3) 100 integer,intent(out) :: hel(2,3) 101 real(dp),intent(in) :: detovc(2,2,3),detovd(2,2,3),efield_dot(3) 102 real(dp),intent(inout) :: dphase_aux1(3),phase_init(3) 103 real(dp),intent(out) :: bcut(2,3),phase_end(3) 104 105 !Local variables ------------------------- 106 !scalars 107 integer :: idir,ifor,igrid,iter,maxiter,ngrid 108 real(dp) :: aa,angle,bb,big_axis,cc,cphi_0,delta_theta,e0,e1 109 real(dp) :: excentr,iab,phase0,phase_min,phi_0,rdum,sgn,small_axis,sphi_0 110 real(dp) :: theta,theta_0,val 111 logical :: flag_neg 112 character(len=500) :: message 113 !arrays 114 real(dp) :: g_theta(2),theta_min(2),theta_try(2) 115 real(dp) :: esave(251),e1save(251) !!REC 116 117 ! *********************************************************************** 118 119 !Compute the helicity and the branch cut of the ellipse in the complex 120 !plane associated with the overlap between a k-point and one of its neighbours 121 122 do idir = 1, 3 123 124 if (abs(efield_dot(idir)) < tol12) cycle 125 126 do ifor = 1, 2 127 128 aa = half*(detovc(1,ifor,idir)*detovc(1,ifor,idir) + & 129 & detovc(2,ifor,idir)*detovc(2,ifor,idir) + & 130 & detovd(1,ifor,idir)*detovd(1,ifor,idir) + & 131 & detovd(2,ifor,idir)*detovd(2,ifor,idir)) 132 133 bb = half*(detovc(1,ifor,idir)*detovc(1,ifor,idir) + & 134 & detovc(2,ifor,idir)*detovc(2,ifor,idir) - & 135 & detovd(1,ifor,idir)*detovd(1,ifor,idir) - & 136 & detovd(2,ifor,idir)*detovd(2,ifor,idir)) 137 138 cc = detovc(1,ifor,idir)*detovd(1,ifor,idir) + & 139 & detovc(2,ifor,idir)*detovd(2,ifor,idir) 140 141 iab = detovc(1,ifor,idir)*detovd(2,ifor,idir) - & 142 & detovc(2,ifor,idir)*detovd(1,ifor,idir) 143 144 if (iab >= zero) then 145 hel(ifor,idir) = 1 146 else 147 hel(ifor,idir) = -1 148 end if 149 150 if (abs(bb) > tol8) then 151 theta_0 = half*atan(cc/bb) 152 else 153 theta_0 = quarter*pi 154 end if 155 156 if (bb < zero) theta_0 = theta_0 + pi*half 157 158 g_theta(:) = cos(theta_0)*detovc(:,ifor,idir) + & 159 & sin(theta_0)*detovd(:,ifor,idir) 160 ! DEBUG 161 ! write(std_out,*)'before rhophi, g_theta =',g_theta 162 ! ENDDEBUG 163 call rhophi(g_theta,phi_0,rdum) 164 ! DEBUG 165 ! write(std_out,*)'after rhophi, phi_0 = ',phi_0 166 ! ENDDEBUG 167 168 cphi_0 = cos(phi_0) 169 sphi_0 = sin(phi_0) 170 171 rdum = aa - sqrt(bb*bb + cc*cc) 172 if (rdum < zero) rdum = zero 173 small_axis = sqrt(rdum) 174 big_axis = sqrt(aa + sqrt(bb*bb + cc*cc)) 175 excentr = hel(ifor,idir)*small_axis/big_axis 176 177 ! Find angle for which phi = pi 178 if (abs(excentr) > tol8) then 179 angle = atan(tan(pi-phi_0)/excentr) 180 else 181 if (tan(pi-phi_0)*hel(ifor,idir) > zero) then 182 angle = half*pi 183 else 184 angle = -0.5_dp*pi 185 end if 186 end if 187 bcut(ifor,idir) = angle + theta_0 188 189 190 ! Compute the branch-cut angle 191 if (hel(ifor,idir) == 1) then 192 if ((sphi_0 > 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) + pi 193 if ((sphi_0 < 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) - pi 194 else 195 if ((sphi_0 > 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) - pi 196 if ((sphi_0 < 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) + pi 197 end if 198 199 if (bcut(ifor,idir) > pi) bcut(ifor,idir) = bcut(ifor,idir) - two_pi 200 if (bcut(ifor,idir) < -1_dp*pi) bcut(ifor,idir) = bcut(ifor,idir) + two_pi 201 202 ! DEBUG 203 ! write(std_out,'(a,2x,i3,2x,i3,5x,f16.9,5x,i2)')'linemin: ifor,idir,bcut,hel',& 204 ! & ifor,idir,bcut(ifor,idir),hel(ifor,idir) 205 ! write(std_out,'(a,5x,f16.9,5x,f16.9)')'linemin: big_axis,small_axis ',& 206 ! & big_axis,small_axis 207 ! ENDDEBUG 208 209 end do ! ifor 210 end do ! idir 211 212 !--------------------------------------------------------------------------- 213 214 !Perform the "line minimization" w.r.t. the angle theta on a unit circle 215 !to update the wavefunction associated with the current k-point and band label. 216 217 ngrid = 250 ! initial number of subdivisions in [-pi/2,pi/2] 218 !for finding extrema 219 maxiter = 100 220 delta_theta = pi/ngrid 221 222 !DEBUG 223 !write(std_out,*)'linemin: theta, e0, e1, e1fdiff' 224 !ENDDEBUG 225 226 227 !Get the interval where the absolute minimum of E(theta) is located 228 229 val = huge(one) ! large number 230 flag_neg=.false. 231 theta_min(:) = ten 232 do igrid = 1, ngrid+1 233 234 theta = (igrid - 1)*delta_theta - pi*half 235 call etheta(bcut,chc,detovc,detovd,dhc,dhd,efield_dot,e0,e1,& 236 & hel,nkpt,nstr,sdeg,theta) 237 238 esave(igrid)=e0 !!REC 239 e1save(igrid)=e1 !!REC 240 241 ! It is important to detect when the slope changes from negative to positive 242 ! Moreover, a slope being extremely close to zero must be ignored 243 244 ! DEBUG 245 ! write(std_out,*)' igrid,e0,e1,val,theta_min(:)=',igrid,theta,e0,e1,val,theta_min(:) 246 ! ENDDEBUG 247 248 ! Store e1 and theta if negative ... 249 if(e1 < -tol10)then 250 theta_try(1)=theta 251 flag_neg=.true. 252 end if 253 ! A change of sign is just happening 254 if(e1 > tol10 .and. flag_neg)then 255 theta_try(2)=theta 256 flag_neg=.false. 257 ! Still, must be better than the previous minimum in order to succeed 258 if (e0 < val-tol10) then 259 val=e0 260 theta_min(:)=theta_try(:) 261 end if 262 end if 263 end do 264 265 !In case the minimum was not found 266 267 if (abs(theta_min(1) - ten) < tol10) then 268 ! REC start 269 write(message,'(a,a)')ch10,& 270 & ' linemin: ERROR- cannot find theta_min.' 271 call wrtout(std_out,message,'COLL') 272 write(message,'(a,a)')ch10,& 273 & ' igrid theta esave(igrid) e1save(igrid) ' 274 call wrtout(std_out,message,'COLL') 275 do igrid = 1, ngrid+1 276 theta = (igrid - 1)*delta_theta - pi*half 277 write(std_out,'(i6,3f16.9)')igrid,theta,esave(igrid),e1save(igrid) 278 !write(101,'(i6,3f16.9)')igrid,theta,esave(igrid),e1save(igrid) 279 end do 280 write(message,'(6a)')ch10,& 281 & ' linemin : ERROR - ',ch10,& 282 & ' Cannot find theta_min. No minimum exists : the field is too strong ! ',ch10,& 283 & ' Try decreasing difference between D and 4 Pi P by changing structure or D (only for fixed D calculation)' 284 call wrtout(std_out,message,'COLL') 285 message = ' linemin cannot find theta_min' 286 MSG_ERROR(message) 287 end if 288 289 !Compute the mimum of E(theta) 290 291 292 iter = 0 293 do while ((delta_theta > tol8).and.(iter < maxiter)) 294 delta_theta = half*(theta_min(2) - theta_min(1)) 295 theta = theta_min(1) + delta_theta 296 call etheta(bcut,chc,detovc,detovd,dhc,dhd,efield_dot,e0,e1,& 297 & hel,nkpt,nstr,sdeg,theta) 298 if (e1 > zero) then 299 theta_min(2) = theta 300 else 301 theta_min(1) = theta 302 end if 303 iter = iter + 1 304 305 ! DEBUG 306 ! write(std_out,'(a,2x,i3,2(2x,f16.9))')'iter,e0,e1 = ',iter,e0,e1 307 ! ENDDEBUG 308 309 end do 310 311 costh = cos(theta) 312 sinth = sin(theta) 313 314 thetam = theta 315 316 !DEBUG 317 !write(std_out,*)'linemin : thetam = ',thetam 318 !ENDDEBUG 319 320 321 !--------------------------------------------------------------------------- 322 323 !Compute and store the change in electronic polarization 324 325 sgn = one 326 do idir = 1, 3 327 328 if (abs(efield_dot(idir)) < tol12) cycle 329 330 phase_end(idir) = zero 331 do ifor = 1, 2 332 333 g_theta(:) = detovc(:,ifor,idir) 334 ! DEBUG 335 ! write(std_out,*)'before rhophi (2nd call), g_theta =',g_theta 336 ! ENDDEBUG 337 call rhophi(g_theta,phase0,rdum) 338 ! DEBUG 339 ! write(std_out,*)'after rhophi, phase0 = ',phase0 340 ! ENDDEBUG 341 342 if(iline == 1) phase_init(idir) = phase_init(idir) + sgn*phase0 343 344 g_theta(:) = costh*detovc(:,ifor,idir) + sinth*detovd(:,ifor,idir) 345 call rhophi(g_theta,phase_min,rdum) 346 347 phase_end(idir) = phase_end(idir) + sgn*phase_min 348 349 ! Correct for branch cuts (remove jumps) 350 if (bcut(ifor,idir) <= zero) phase0 = phase0 + hel(ifor,idir)*two_pi 351 if(thetam >= bcut(ifor,idir)) phase_min = phase_min + hel(ifor,idir)*two_pi 352 353 dphase_aux1(idir) = dphase_aux1(idir) + sgn*(phase_min - phase0) 354 355 sgn = -1_dp*sgn 356 357 end do ! idir 358 end do ! ifor 359 360 !DEBUG 361 !write(std_out,'(a,3(2x,f16.9))')'dphase_aux1 = ',(dphase_aux1(idir),idir = 1, 3) 362 !write(std_out,*)' linemin: debug, exit.' 363 !ENDDEBUG 364 365 end subroutine linemin