1!#define __DEBUG_V_H_vs_SPHEROPOLE 2! 3! Copyright (C) 2004 PWSCF group 4! This file is distributed under the terms of the 5! GNU General Public License. See the file `License' 6! in the root directory of the present distribution, 7! or http://www.gnu.org/copyleft/gpl.txt . 8! 9MODULE atomic_paw 10 ! 11 !============================================================================ 12 ! 13 ! Module for Projector Augmented Wave (PAW) calculations assuming 14 ! spherical symmetry. Kresse's notations are adopted. 15 ! Contains the type describing a PAW dataset, the routine for 16 ! generating it from the ld1 code, and the routines to compute 17 ! the hamiltonian and energy. 18 ! GGA and LSD are implemented, relativistic calculations are not. 19 ! 20 ! References: 21 ! P.E.Blochl, Phys. Rev. B 50, 17953 (1994) 22 ! G.Kresse, D.Joubert, Phys. Rev. B 59, 1758 (1999) 23 ! 24 ! WARNINGS: 25 ! Still work in progress on many aspects. 26 ! The PAW dataset is written in a temporary format which is not supposed 27 ! to be used any longer. 28 ! Consistency with the input of the ld1 code yet to be checked 29 ! 30 ! Written by Guido Fratesi, february 2005 31 ! Modified by Riccardo Mazzarello, july 2006 32 ! Fully Relativistic generalization by Andrea Dal Corso, november 2009 33 ! Other people involved: Lorenzo Paulatto and Stefano de Gironcoli 34 ! 35 !============================================================================ 36 ! 37 USE kinds, ONLY: dp 38 USE ld1_parameters, ONLY: nwfsx 39 USE upf_params, ONLY: lmaxx 40 USE constants, ONLY: pi, fpi, e2, eps8 41 USE radial_grids, ONLY: ndmx, radial_grid_type 42 USE paw_type, ONLY: paw_t, nullify_pseudo_paw, allocate_pseudo_paw 43 ! 44 IMPLICIT NONE 45 PRIVATE 46 SAVE 47 ! 48 REAL(dp), PARAMETER :: ZERO=0._dp, ONE=1._dp, TWO=2._dp, HALF=0.5_dp 49 ! 50 ! TEMP, to be substituted by module constants 51! REAL(dp), PARAMETER :: & 52! PI=3.14159265358979323846_dp, FPI=4._dp*PI, E2=2._dp, EPS8=1.0e-8_dp 53 ! 54 !============================================================================ 55 ! 56 PUBLIC :: paw_t 57 PUBLIC :: us2paw 58 PUBLIC :: paw2us 59 PUBLIC :: check_multipole 60 PUBLIC :: new_paw_hamiltonian 61 PUBLIC :: find_bes_qi 62 PUBLIC :: compute_nonlocal_coeff_ion 63 ! 64CONTAINS 65 ! 66 !============================================================================ 67 ! PUBLIC ROUTINES !!! 68 !============================================================================ 69 ! 70 ! Compute the values of the local pseudopotential and the NL coefficients 71 ! Compute also the total energy 72 ! 73 SUBROUTINE new_paw_hamiltonian (veffps_, ddd_, etot_, & 74 pawset_, nwfc_, l_, j_, nspin_, spin_, oc_, pswfc_, eig_, paw_energy,dddion_) 75 IMPLICIT NONE 76 REAL(dp), INTENT(OUT) :: veffps_(ndmx,2) 77 REAL(dp), INTENT(OUT) :: ddd_(nwfsx,nwfsx,2) 78 REAL(dp), INTENT(OUT) :: etot_ 79 TYPE(paw_t), INTENT(IN) :: pawset_ 80 INTEGER, INTENT(IN) :: nwfc_ 81 INTEGER, INTENT(IN) :: l_(nwfsx) 82 INTEGER, INTENT(IN) :: nspin_ 83 INTEGER, INTENT(IN) :: spin_(nwfsx) 84 REAL(dp), INTENT(IN) :: j_(nwfsx) 85 REAL(dp), INTENT(IN) :: oc_(nwfsx) 86 REAL(dp), INTENT(IN) :: pswfc_(ndmx,nwfsx) 87 REAL(dp), INTENT(IN) :: eig_(nwfsx) 88 REAL(dp), OPTIONAL :: dddion_(nwfsx,nwfsx) 89 REAL(dp), INTENT(OUT), OPTIONAL :: paw_energy(5,3) 90 ! 91 REAL(dp) :: & ! one center: 92 eps, e1, e1ps, & ! energies; 93 veff1(ndmx,2), veff1ps(ndmx,2), & ! eff potentials; 94 chargeps(ndmx,2), charge1(ndmx,2), charge1ps(ndmx,2), & ! charges. 95 projsum(nwfsx,nwfsx,2), eigsum ! sum of projections, sum of eigenval. 96 ! 97 INTEGER :: ns, is, n 98 REAL(dp) :: energy(5,3) 99 ! 100 ! Compute the valence charges 101 CALL compute_charges(projsum, chargeps, charge1, charge1ps, & 102 pawset_, nwfc_, l_, j_, nspin_, spin_, oc_, pswfc_, 1 ) 103 ! 104 ! Check for negative charge 105 ! 106 do is=1,nspin_ 107 do n=2,pawset_%grid%mesh 108! write(6,*) n, pawset_%grid%r(n), chargeps(n,is) 109 if (chargeps(n,is)<-1.d-12) & 110 call errore('new_paw_hamiltonian','negative rho',1) 111 enddo 112 enddo 113 114! write(766,"(4f12.6)") (pawset_%grid%r(ns), chargeps(ns,1), charge1(ns,1), charge1ps(ns,1), ns=1,pawset_%grid%mesh) 115! write(767,"(4f12.6)") (pawset_%grid%r(ns), chargeps(ns,2), charge1(ns,2), charge1ps(ns,2), ns=1,pawset_%grid%mesh) 116 ! 117 ! Compute the one-center energy and effective potential: 118 ! E = Eh[n_v] + Exc[n_v+n_c] - Int[veff*n_v], 119 ! veff = vh[n_v] + v_xc[n_v+n_c], 120 ! where n_v can be n~+n^ or n1 or n1~+n^ and n_c can be nc or nc~ 121 ! n~+n^ , nc~ 122 CALL compute_onecenter_energy ( eps, veffps_, & 123 pawset_, chargeps, pawset_%nlcc, pawset_%psccharge, nspin_,& 124 pawset_%grid%mesh, pawset_%psloc, energy(1,1)) 125 ! n1 , nc 126 CALL compute_onecenter_energy ( e1, veff1, & 127 pawset_, charge1, .TRUE., pawset_%aeccharge, nspin_,& 128 pawset_%irc, pawset_%aeloc, energy(1,2)) 129 ! n1~+n^ , nc~ 130 CALL compute_onecenter_energy ( e1ps, veff1ps, & 131 pawset_, charge1ps, pawset_%nlcc, pawset_%psccharge, nspin_,& 132 pawset_%irc, pawset_%psloc, energy(1,3)) 133 ! Add the local part 134 DO is=1,nspin_ 135 veffps_(1:pawset_%grid%mesh,is) = veffps_(1:pawset_%grid%mesh,is) + & 136 pawset_%psloc(1:pawset_%grid%mesh) 137 veff1 (1:pawset_%grid%mesh,is) = veff1 (1:pawset_%grid%mesh,is) + & 138 pawset_%aeloc(1:pawset_%grid%mesh) 139 veff1ps(1:pawset_%grid%mesh,is) = veff1ps(1:pawset_%grid%mesh,is) + & 140 pawset_%psloc(1:pawset_%grid%mesh) 141 END DO 142 ! 143 ! Compute the nonlocal D coefficients 144 CALL compute_nonlocal_coeff (ddd_,pawset_,nspin_,veffps_,veff1,veff1ps) 145 IF (PRESENT(dddion_)) THEN 146 CALL compute_nonlocal_coeff_ion (dddion_,pawset_) 147 END IF 148 ! 149 ! Compute total energy 150 eigsum=ZERO 151 DO ns=1,nwfc_ 152 IF (oc_(ns)>ZERO) eigsum = eigsum + oc_(ns)*eig_(ns) 153 END DO 154 etot_ = eigsum + eps + e1 - e1ps 155 156 if (PRESENT(paw_energy)) paw_energy=energy 157 ! 158 END SUBROUTINE new_paw_hamiltonian 159 ! 160 !============================================================================ 161 ! 162 ! Convert the USPP to a PAW dataset 163 ! 164 SUBROUTINE us2paw (pawset_, & 165 zval, grid, rmatch_augfun, ikk, & 166 nbeta, lls, jjs, ocs, enls, els, rcutus, psipaw, psipaw_rel, & 167 phis, betas, qvan, kindiff, & 168 nlcc, aerhoc, psrhoc, aevtot, psvtot, which_paw_augfun,rel ) 169 170 USE funct, ONLY : get_dft_long 171 USE ld1inc, ONLY : zed, file_screen 172 USE paw_type, ONLY : nullify_pseudo_paw, allocate_pseudo_paw 173 USE io_global, ONLY : stdout, ionode, ionode_id 174 USE radial_grids, ONLY : allocate_radial_grid 175 USE mp, only : mp_bcast 176 USE mp_world, only : world_comm 177 IMPLICIT NONE 178 TYPE(paw_t), INTENT(OUT) :: pawset_ 179 REAL(dp), INTENT(IN) :: zval 180 type(radial_grid_type), INTENT(IN) :: grid 181 REAL(dp), INTENT(IN) :: rmatch_augfun 182 INTEGER, INTENT(IN) :: ikk(nwfsx) 183 ! 184 INTEGER, INTENT(IN) :: nbeta, rel 185 INTEGER, INTENT(IN) :: lls(nwfsx) 186 REAL(dp), INTENT(IN) :: ocs(nwfsx) 187 CHARACTER(LEN=2), INTENT(IN) :: els(nwfsx) 188 REAL(dp), INTENT(IN) :: jjs(nwfsx) 189 REAL(dp), INTENT(IN) :: rcutus(nwfsx) 190 REAL(dp), INTENT(IN) :: enls(nwfsx) 191 REAL(dp), INTENT(IN) :: psipaw(ndmx,nwfsx) 192 REAL(dp), INTENT(IN) :: psipaw_rel(ndmx,nwfsx) 193 REAL(dp), INTENT(IN) :: phis(ndmx,nwfsx) 194 REAL(dp), INTENT(IN) :: betas(ndmx,nwfsx) 195 ! 196 REAL(dp), INTENT(IN) :: qvan(ndmx,nwfsx,nwfsx) 197 REAL(dp), INTENT(IN) :: kindiff(nwfsx,nwfsx) 198 LOGICAL, INTENT(IN) :: nlcc 199 REAL(dp), INTENT(IN) :: aerhoc(ndmx) 200 REAL(dp), INTENT(IN) :: psrhoc(ndmx) 201 REAL(dp), INTENT(IN) :: aevtot(ndmx) 202 REAL(dp), INTENT(IN) :: psvtot(ndmx) 203 CHARACTER(20), INTENT(IN) :: which_paw_augfun 204 ! 205 REAL(DP), EXTERNAL :: int_0_inf_dr 206 CHARACTER(len=2), EXTERNAL :: atom_name 207 REAL(dp) :: vps(ndmx,2), projsum(nwfsx,nwfsx,2), ddd(nwfsx,nwfsx,2), dddion(nwfsx,nwfsx) 208 INTEGER :: irc, ns, ns1, n, leading_power, mesh, ios 209 REAL(dp) :: aux(ndmx), aux2(ndmx,2), raux 210 REAL(dp) :: aecharge(ndmx,2), pscharge(ndmx,2) 211 REAL(dp) :: etot 212 INTEGER :: nspin=1, spin(nwfsx)=1 ! PAW generat. from spin-less calculation 213 ! 214 ! variables for aug. functions generation 215 ! 216 REAL(DP) :: energy(5,3), max_aug_cutoff = -1._dp 217 INTEGER :: nc, iok ! index Bessel function, ... 218 INTEGER :: l1,l2,l3, lll, ircm, ir, ir0 219 REAL(dp) :: twosigma2, rm ! needed for "gaussian" augfun 220 REAL(dp) :: zeta, resi,usigma=4._dp ! needed for "Bloechl" augfun 221 REAL(DP),ALLOCATABLE :: gaussian(:) ! needed for gaussian and Bloechl 222 REAL(DP),ALLOCATABLE :: aug_real(:,:,:) ! needed for PSQ augfun 223 REAL(DP) :: qc(2), xc(2), b1(2), b2(2) ! needed for BESSEL augfun 224 REAL(DP), ALLOCATABLE :: j1(:,:) ! needed for BESSEL augfun 225 ! 226 mesh = grid%mesh 227 irc = maxval(ikk(1:nbeta))+8 228 CALL nullify_pseudo_paw(pawset_) 229 CALL allocate_pseudo_paw(pawset_,ndmx,nwfsx,lmaxx) 230 CALL allocate_radial_grid(pawset_%grid, mesh) 231 pawset_%symbol = atom_name(nint(zed)) 232 pawset_%zval = zval 233 ! 234 ! Copy the mesh 235 pawset_%grid%mesh = grid%mesh 236 pawset_%grid%r(:) = grid%r(:) 237 pawset_%grid%r2(:) = grid%r2(:) 238 pawset_%grid%rab(:)= grid%rab(:) 239 pawset_%grid%sqr(:)= grid%sqr(:) 240 pawset_%grid%xmin = grid%xmin 241 pawset_%grid%rmax = grid%rmax 242 pawset_%grid%zmesh = grid%zmesh 243 pawset_%grid%dx = grid%dx 244 ! 245 pawset_%rmatch_augfun = rmatch_augfun 246 !if (rmatch_augfun <= 0.0_dp) pawset_%rmatch_augfun = grid%r(irc) 247 if (rmatch_augfun <= 0.0_dp) pawset_%rmatch_augfun = MAXVAL(rcutus(1:nbeta)) 248 pawset_%rel = rel 249 pawset_%irc = irc 250 pawset_%ikk(1:nbeta)=ikk(1:nbeta) 251 ! 252 ! Copy the wavefunctions 253 pawset_%nwfc = nbeta 254 pawset_%l(1:nbeta) = lls(1:nbeta) 255 pawset_%lmax = MAXVAL(pawset_%l(1:nbeta)) 256 pawset_%oc(1:nbeta) = ocs(1:nbeta) 257 pawset_%jj(1:nbeta) = jjs(1:nbeta) 258 pawset_%rcutus(1:nbeta) = rcutus(1:nbeta) 259 pawset_%els(1:nbeta)= els(1:nbeta) 260 pawset_%enl(1:nbeta)= enls(1:nbeta) 261 pawset_%aewfc(1:mesh,1:nbeta) = psipaw(1:mesh,1:nbeta) 262 pawset_%aewfc_rel(1:mesh,1:nbeta) = psipaw_rel(1:mesh,1:nbeta) 263 pawset_%pswfc(1:mesh,1:nbeta) = phis (1:mesh,1:nbeta) 264 pawset_%proj (1:mesh,1:nbeta) = betas (1:mesh,1:nbeta) 265 ! 266 pawset_%augfun = 0._dp 267 ! 268 ! 269 ! Augmentation functions: 270 ! The specific radial part is not important, as long as the 271 ! multipole moments of the exact augmentation functions are 272 ! reproduced. So we write on the file the exact augmentation 273 ! functions and their multipole moments, keeping in mind that 274 ! the PW program should not use the radial part as is but 275 ! substitute it with some smoothened analogue. 276 ! 277 ! Sdg>> not quite sure this is correct because the shape of augfun, 278 ! arbitrary as it may be, determines the shape of vloc since this 279 ! is obtained by unscreening vscf with the GIVEN augfun... 280 ! 281 ! moreover in order to give the right electrostatics in the solid the 282 ! augmentation functions should not overlap. 283 ! 284 ! Compute the exact augmentation functions and their moments 285 ! 286 write(stdout,'(//,4x,a)') 'multipoles (all-electron charge) - (pseudo charge)' 287 write(stdout,'(7x,2a,":",2a,2x,6a)') ' ns',' l1','ns1',' l2',& 288 ' l=0 ',' l=1 ',' l=2 ',' l=3 ',' l=4 ', ' l=5 ' 289 pawset_%augfun(:,:,:,:) = 0.0_dp 290 pawset_%augmom(:,:,:) = 0.0_dp 291 DO ns=1,nbeta 292 l1=pawset_%l(ns) 293 DO ns1=1,ns 294 l2=pawset_%l(ns1) 295 do l3 = max(l1-l2,l2-l1), l1+l2, 2 296 pawset_%augfun(1:mesh,ns,ns1,l3) = & 297 pawset_%aewfc(1:mesh,ns) * pawset_%aewfc(1:mesh,ns1) - & 298 pawset_%pswfc(1:mesh,ns) * pawset_%pswfc(1:mesh,ns1) 299 IF (pawset_%rel==2) & 300 pawset_%augfun(1:irc,ns,ns1,l3) =pawset_%augfun(1:irc,ns,ns1,l3)& 301 +pawset_%aewfc_rel(1:irc,ns) * pawset_%aewfc_rel(1:irc,ns1) 302 pawset_%augfun(1:mesh,ns1,ns,l3) = pawset_%augfun(1:mesh,ns,ns1,l3) 303 aux(1:irc) = pawset_%augfun(1:irc,ns,ns1,l3) * pawset_%grid%r(1:irc)**l3 304 lll = l1 + l2 + 2 + l3 305 pawset_%augmom(ns,ns1,l3)=int_0_inf_dr(aux(1:irc),pawset_%grid,irc,lll) 306 pawset_%augmom(ns1,ns,l3)=pawset_%augmom(ns,ns1,l3) 307 end do 308 write (stdout,'(7x,2i3,a,2i3,10f8.4)') ns,l1,":",ns1,l2,& 309 (pawset_%augmom(ns,ns1,l3), l3=0,l1+l2) 310 END DO 311 END DO 312 ! 313 ! 314 IF (which_paw_augfun/='AE') THEN 315 ! The following lines enables to test the idependence on the 316 ! specific shape of the radial part of the augmentation functions 317 ! in a spherically averager system (as is the case in atoms) only 318 ! the zero-th moment contribute to the scf charge 319 write(stdout, '(5x,3a,f12.6)') 'Required augmentation: ',TRIM(which_paw_augfun), " radius:", pawset_%rmatch_augfun 320 ! 321101 pawset_%augfun(:,:,:,:) = 0.0_dp 322 DO ns=1,nbeta 323 l1 = pawset_%l(ns) 324 DO ns1=1,ns 325 l2 = pawset_%l(ns1) 326 ! 327 SELECT CASE (which_paw_augfun) 328 CASE ('PSQ') 329 continue ! the work is done at the end 330 CASE ('QVAN') 331 IF(ns==1 .and. ns1==1) & 332 CALL infomsg('us2paw', 'WARNING: QVAN augmentation function are for testing ONLY: '//& 333 'they will not work in pw!') 334 ! Choose the shape of the augmentation functions: NC Q ... 335 pawset_%augfun(1:mesh,ns,ns1,0) = qvan(1:mesh,ns,ns1) 336 ! Renormalize the aug. functions so that their integral is correct 337 leading_power = l1 + l2 + 2 338 raux=int_0_inf_dr(pawset_%augfun(1:irc,ns,ns1,0),pawset_%grid,irc,leading_power) 339 IF (ABS(raux) < eps8) CALL errore & 340 ('ld1_to_paw','norm of augm.func. is too small',ns*100+ns1) 341 raux=pawset_%augmom(ns,ns1,0)/raux 342 pawset_%augfun(1:mesh,ns,ns1,0)=raux*pawset_%augfun(1:mesh,ns,ns1,0) 343 ! 344 CASE ('BG') 345 IF(ns==1 .and. ns1==1) & 346 CALL infomsg('us2paw', 'WARNING: using Bloechl style augmentation functions '//& 347 'is not a good idea, as analytical overlap are not '//& 348 'implemented in pwscf: use BESSEL or GAUSS instead.') 349 IF(.not. allocated(gaussian)) ALLOCATE(gaussian(mesh)) 350 ! use Bloechl style augmentation functions, as linear combinations of 351 ! gaussians functions (this is quite pointless if the the plane-wave 352 ! code doesn't use double-augmentation with analytical gaussian overlaps) 353 DO ir0 = 1,mesh 354 IF(grid%r(ir0) > pawset_%rmatch_augfun) & 355 exit 356 ENDDO 357 ! look for a sigma large enough to keep (almost) all the gaussian in the aug sphere 358 zeta = (usigma/pawset_%rmatch_augfun)**2 359 DO ir = 1,mesh 360 gaussian(ir) = exp( - zeta*(grid%r(ir))**2 ) * grid%r2(ir) 361 ENDDO 362 DO l3 = max (l1-l2,l2-l1), l1+l2 363 ! Functions has to be normalized later, so I can use a constant factor 364 ! = rc**l3 to prevent very large numbers when integrating: 365 aux(1:grid%mesh) = gaussian(1:grid%mesh) * grid%r(1:grid%mesh)**l3 366 ! Normalization to have unitary multipole l3 367 ! and check norm of augmentation functions. 368 raux = int_0_inf_dr(aux,pawset_%grid,mesh,l3+2) 369 IF (abs(raux) < eps8) CALL errore & 370 ('ld1_to_paw','norm of augm. func. is too small',ns*100+ns1) 371 ! Check if gaussians are localized enough into the augmentation sphere, 372 ! otherwise try again with a larger sigma. 373 resi = int_0_inf_dr(aux,pawset_%grid,ir0, l3+2) 374 resi = (raux-resi)/raux 375 IF (abs(resi) > 1.e-10_dp) THEN 376 usigma = usigma + .0625_dp 377 GOTO 101 378 ENDIF 379 raux=pawset_%augmom(ns,ns1,l3)/raux 380 381 pawset_%augfun(1:mesh,ns,ns1,l3) = raux*gaussian(1:mesh) 382 pawset_%augfun(1:mesh,ns1,ns,l3) = raux*gaussian(1:mesh) 383 ENDDO 384 DEALLOCATE(gaussian) 385 ! 386 CASE ('GAUSS') 387 ! use linear combinations of gaussians functions, not the Bloechl style 388 ! but it looks a bit alike... (used for testing, now obsolete) 389 IF(ns==1 .and. ns1==1) & 390 CALL infomsg('us2paw', 'GAUSS augmentation functions are ususally '//& 391 'harder than BESSEL; use BESSEL instead unless'//& 392 ' you have discontinuity in local potential') 393 ALLOCATE(gaussian(mesh)) 394 ! 395 rm = pawset_%rmatch_augfun 396 twosigma2 = TWO*(rm/3.0_dp)**2 397 do ir=1,mesh 398 if (grid%r(ir) <= rm) then 399 gaussian(ir) = ( EXP(-grid%r(ir)**2/twosigma2) + & 400 EXP(-(grid%r(ir)-TWO*rm)**2/twosigma2 ) - & 401 TWO*EXP(-rm**2/twosigma2) ) * grid%r2(ir) 402 else 403 gaussian(ir) = 0.0_dp 404 endif 405 end do 406 DO l3 = max (l1-l2,l2-l1), l1+l2 407 ! 408 aux(1:irc) = gaussian(1:irc) * pawset_%grid%r(1:irc)**l3 409 ! calculate the corresponding multipole 410 raux=int_0_inf_dr(aux,pawset_%grid,irc,l3+2) 411 IF (ABS(raux) < eps8) CALL errore & 412 ('ld1_to_paw','norm of augm. func. is too small',ns*100+ns1) 413 raux=pawset_%augmom(ns,ns1,l3)/raux 414 pawset_%augfun(1:mesh,ns,ns1,l3) = raux*gaussian(1:mesh) 415 pawset_%augfun(1:mesh,ns1,ns,l3) = raux*gaussian(1:mesh) 416 ! 417 END DO 418 DEALLOCATE(gaussian) 419 ! 420 CASE ('BESSEL') 421 ! Use Kresse-Joubert style augmentation functions 422 ! Defined as linear combination of Bessel functions. 423 ALLOCATE (j1 (pawset_%grid%mesh,2)) 424 do ir=1,irc 425 !if (grid%r(ir)<pawset_%rmatch_augfun) ircm=ir 426 if (grid%r(ir)<pawset_%rmatch_augfun) ircm=ir 427 end do 428 do l3 = max(l1-l2,l2-l1), l1+l2 429 ! 430 CALL find_bes_qi(qc,pawset_%grid%r(ircm),l3,2,iok) 431 IF (iok.ne.0) CALL errore('compute_augfun', & 432 'problems with find_bes_qi',1) 433 DO nc = 1, 2 434 ! 435 CALL sph_bes(irc,grid%r(1),qc(nc),l3,j1(1,nc)) 436 aux(1:irc) = j1(1:irc,nc) * grid%r(1:irc)**(l3+2) 437 b1(nc) = j1(ircm,nc) 438 b2(nc) = int_0_inf_dr(aux,pawset_%grid,ircm,l3+2) 439 ! 440 ENDDO 441 xc(1) = b1(2) / (b1(2) * b2(1) - b1(1) * b2(2)) 442 xc(2) = - b1(1) * xc(1) / b1(2) 443 pawset_%augfun(1:ircm,ns,ns1,l3) = & 444 pawset_%augmom(ns,ns1,l3) * grid%r2(1:ircm) * & 445 (xc(1) * j1(1:ircm,1) + xc(2) * j1(1:ircm,2)) 446 ! Symmetrize augmentation functions: 447 pawset_%augfun(1:mesh,ns1,ns,l3)=pawset_%augfun(1:mesh,ns,ns1,l3) 448 ! 449 ! Save higher bessel coefficient to compute suggested cutoff 450 max_aug_cutoff=MAX( max_aug_cutoff, MAXVAL(qc(1:2))**2) 451 ! 452 end do 453 DEALLOCATE (j1) 454 ! 455 CASE DEFAULT 456 ! 457 CALL errore ('ld1_to_paw','Specified augmentation functions ('// & 458 TRIM(which_paw_augfun)//') not allowed or coded',1) 459 ! 460 END SELECT 461 END DO 462 END DO 463 IF ( which_paw_augfun == 'BG') & 464 write(stdout,"(5x,a,f12.6)") "Gaussians generated with zeta: ", zeta 465 IF ( which_paw_augfun == 'BESSEL') & 466 write(stdout,'(5x, "Suggested rho cutoff for augmentation:",f7.2," Ry")') max_aug_cutoff 467 ! 468 IF ( which_paw_augfun == 'PSQ') THEN 469 ALLOCATE(aug_real(ndmx,nwfsx,nwfsx)) 470 DO ns=1,nbeta 471 DO ns1=ns,nbeta 472 aug_real(1:mesh,ns,ns1) = & 473 pawset_%aewfc(1:mesh,ns) * pawset_%aewfc(1:mesh,ns1) - & 474 pawset_%pswfc(1:mesh,ns) * pawset_%pswfc(1:mesh,ns1) 475 IF (pawset_%rel==2) & 476 aug_real(1:irc,ns,ns1) = aug_real(1:irc,ns,ns1) + & 477 pawset_%aewfc_rel(1:irc,ns)*pawset_%aewfc_rel(1:irc,ns1) 478 aug_real(1:mesh,ns1,ns) = aug_real(1:mesh,ns,ns1) 479 ENDDO 480 ENDDO 481 ! 482 CALL pseudo_q(aug_real,pawset_%augfun) 483 ! 484 DEALLOCATE(aug_real) 485 ENDIF 486 487 END IF 488! 489! Outside the PAW spheres augfun should be exactly 0. On some machine 490! it is equal to zero to machine precision and sometimes it is negative, 491! so as to confuse the check for negative charge. So we set it to zero 492! explicitly. 493! 494 DO ns = 1, nbeta 495 l1 = pawset_%l(ns) 496 DO ns1 = ns, nbeta 497 l2 = pawset_%l(ns1) 498 DO l3 = max (l1-l2,l2-l1), l1+l2 499 DO n = pawset_%irc+1, pawset_%grid%mesh 500 pawset_%augfun(n,ns,ns1,l3) = 0.0_DP 501 pawset_%augfun(n,ns1,ns,l3) = 0.0_DP 502 END DO 503 END DO 504 END DO 505 END DO 506 ! 507 ! 508 ! Copy kinetic energy differences 509 pawset_%kdiff(1:nbeta,1:nbeta)=kindiff(1:nbeta,1:nbeta) 510 ! 511 ! Copy the core charge (if not NLCC copy only the AE one) 512 pawset_%nlcc=nlcc 513 IF (pawset_%nlcc) pawset_%psccharge(1:mesh)=psrhoc(1:mesh) 514 pawset_%aeccharge(1:mesh)=aerhoc(1:mesh) 515 ! 516 ! Copy the local potentials 517 pawset_%psloc(1:mesh)=psvtot(1:mesh) 518 pawset_%aeloc(1:mesh)=aevtot(1:mesh) 519 ! and descreen them: 520 CALL compute_charges(projsum, pscharge, aecharge, aux2, & 521 pawset_, nbeta, lls, jjs, nspin, spin, ocs, phis ) 522 pawset_%pscharge(1:mesh)=pscharge(1:mesh,1) 523 ! 524 CALL compute_onecenter_energy ( raux, aux2, & 525 pawset_, pscharge, pawset_%nlcc, pawset_%psccharge, nspin, & 526 pawset_%grid%mesh ) 527 pawset_%psloc(1:mesh)=psvtot(1:mesh)-aux2(1:mesh,1) 528 ! 529 CALL compute_onecenter_energy ( raux, aux2, & 530 pawset_, aecharge, .TRUE., pawset_%aeccharge, nspin, & 531 pawset_%grid%mesh) 532 pawset_%aeloc(1:mesh)=aevtot(1:mesh)-aux2(1:mesh,1) 533 ! 534 if (file_screen .ne.' ') then 535 if (ionode) & 536 open(unit=20,file=file_screen, status='unknown', iostat=ios, err=105 ) 537105 call mp_bcast(ios, ionode_id, world_comm) 538 call errore('descreening','opening file'//file_screen,abs(ios)) 539 if (ionode) then 540 write(20,'(a)') "# n, r(n), aeloc(n), psloc(n), pscharge(n)" 541 do n=1,grid%mesh 542 write(20,'(i5,7e12.4)') n,grid%r(n), pawset_%aeloc(n), & 543 pawset_%psloc(n), pawset_%pscharge(n) 544 enddo 545 close(20) 546 endif 547 endif 548 ! 549 write(pawset_%dft,'(80x)') !fill it with spaces 550 pawset_%dft = get_dft_long ( ) 551 ! 552 ! Generate the paw hamiltonian for test (should be equal to the US one) 553 CALL new_paw_hamiltonian (vps, ddd, etot, & 554 pawset_, pawset_%nwfc, pawset_%l, pawset_%jj, nspin, spin, pawset_%oc, & 555 pawset_%pswfc, pawset_%enl, energy, dddion) 556 pawset_%dion(1:nbeta,1:nbeta)=dddion(1:nbeta,1:nbeta) 557 WRITE(stdout,'(/5x,A,f12.6,A)') 'Estimated PAW energy =',etot,' Ryd' 558 WRITE(stdout,'(/5x,A)') 'The PAW screened D coefficients' 559 DO ns1=1,pawset_%nwfc 560 WRITE(stdout,'(6f12.5)') (ddd(ns1,ns,1),ns=1,pawset_%nwfc) 561 END DO 562 WRITE(stdout,'(/5x,A)') 'The PAW descreened D coefficients (US)' 563 DO ns1=1,pawset_%nwfc 564 WRITE(stdout,'(6f12.5)') (dddion(ns1,ns),ns=1,pawset_%nwfc) 565 END DO 566 ! 567 ! 568 END SUBROUTINE us2paw 569 ! 570 !============================================================================ 571 ! 572 ! ... 573 ! 574 SUBROUTINE paw2us (pawset_,zval,grid,nbeta,lls,jjs,ikk,betas,qq,qvan,& 575 vpsloc, bmat, rhos, els, rcutus, pseudotype,psipaw_rel) 576 USE funct, ONLY : set_dft_from_name 577 use radial_grids, only: radial_grid_type, allocate_radial_grid 578 IMPLICIT NONE 579 TYPE(radial_grid_type), INTENT(OUT) :: grid 580 TYPE(paw_t), INTENT(IN) :: pawset_ 581 REAL(dp), INTENT(OUT) :: zval 582 INTEGER, INTENT(OUT) :: nbeta 583 INTEGER, INTENT(OUT) :: lls(nwfsx) 584 INTEGER, INTENT(OUT) :: ikk(nwfsx) 585 INTEGER, INTENT(OUT) :: pseudotype 586 CHARACTER(LEN=2), INTENT(OUT) :: els(nwfsx) 587 REAL(dp), INTENT(OUT) :: jjs(nwfsx) 588 REAL(dp), INTENT(OUT) :: rcutus(nwfsx) 589 REAL(dp), INTENT(OUT) :: betas(ndmx,nwfsx) 590 REAL(dp), INTENT(OUT) :: psipaw_rel(ndmx,nwfsx) 591 REAL(dp), INTENT(OUT) :: qq(nwfsx,nwfsx) 592 REAL(dp), INTENT(OUT) :: qvan(ndmx,nwfsx,nwfsx) 593 REAL(dp), INTENT(OUT) :: vpsloc(ndmx) ! the local pseudopotential 594 REAL(dp), INTENT(OUT) :: bmat(nwfsx,nwfsx)! the pseudo coefficients (unscreened D) 595 REAL(dp), INTENT(OUT) :: rhos(ndmx) ! the pseudo density 596 597 INTEGER :: ns, ns1, mesh 598 zval=pawset_%zval 599 ! 600 mesh=pawset_%grid%mesh 601 ! Copy the mesh 602 call allocate_radial_grid(grid, mesh) 603 grid%mesh = pawset_%grid%mesh 604 grid%r(1:mesh) = pawset_%grid%r(1:mesh) 605 grid%r2(1:mesh) = pawset_%grid%r2(1:mesh) 606 grid%rab(1:mesh)= pawset_%grid%rab(1:mesh) 607 grid%sqr(1:mesh)= pawset_%grid%sqr(1:mesh) 608 grid%xmin = pawset_%grid%xmin 609 grid%rmax = pawset_%grid%rmax 610 grid%zmesh = pawset_%grid%zmesh 611 grid%dx = pawset_%grid%dx 612 613 ! 614 nbeta=pawset_%nwfc 615 lls(1:nbeta)=pawset_%l(1:nbeta) 616 jjs(1:nbeta)=pawset_%jj(1:nbeta) 617 els(1:nbeta)=pawset_%els(1:nbeta) 618 rcutus(1:nbeta)=pawset_%rcutus(1:nbeta) 619 ikk(1:nbeta)=pawset_%ikk(1:nbeta) 620 ! 621 DO ns=1,nbeta 622 DO ns1=1,nbeta 623 qvan(1:mesh,ns,ns1)=pawset_%augfun(1:mesh,ns,ns1,0) 624 IF (lls(ns)==lls(ns1)) THEN 625 qq(ns,ns1)=pawset_%augmom(ns,ns1,0) 626 ELSE ! different spherical harmonic => 0 627 qq(ns,ns1)=ZERO 628 END IF 629 END DO 630 END DO 631 ! 632 vpsloc(1:mesh)=pawset_%psloc(1:mesh) 633 bmat(1:nbeta,1:nbeta)=pawset_%dion(1:nbeta,1:nbeta) 634 ! 635 rhos(1:mesh)=pawset_%pscharge(1:mesh) 636 ! 637 psipaw_rel(1:mesh,1:nbeta)=pawset_%aewfc_rel(1:mesh,1:nbeta) 638 betas(1:mesh,1:nbeta)=pawset_%proj(1:mesh,1:nbeta) 639 pseudotype=3 640 ! 641 CALL set_dft_from_name (pawset_%dft) 642 643 END SUBROUTINE paw2us 644 ! 645 !============================================================================ 646 ! 647 ! ... 648 ! 649 SUBROUTINE check_multipole (pawset_) 650 USE radial_grids, ONLY: hartree 651 USE io_global, ONLY : stdout 652 IMPLICIT NONE 653 TYPE(paw_t), INTENT(IN) :: pawset_ 654 INTEGER:: mesh 655 REAL(dp) :: r(ndmx), r2(ndmx), sqr(ndmx), dx 656 INTEGER :: nbeta 657 INTEGER :: lls(nwfsx) 658 INTEGER :: ir, ns1, ns2, l1, l2, l3, irc, ir0 659 REAL(dp) :: auxpot(ndmx,0:2*lmaxx+2), auxrho(ndmx) 660 ! 661 ! set a few internal variables 662 write (stdout,*) "check_multipole : lmaxx =",lmaxx 663 mesh=pawset_%grid%mesh 664 r(1:mesh)=pawset_%grid%r(1:mesh) 665 r2(1:mesh)=pawset_%grid%r2(1:mesh) 666 sqr(1:mesh)=pawset_%grid%sqr(1:mesh) 667 dx=pawset_%grid%dx 668 irc = pawset_%irc 669 ! 670 nbeta=pawset_%nwfc 671 lls(1:nbeta)=pawset_%l(1:nbeta) 672 ! 673 do ns1=1,nbeta 674 l1 = lls(ns1) 675 do ns2=1,nbeta 676 l2 = lls(ns2) 677 auxpot(:,:) = 0.0_dp 678 do l3 = max(l1-l2,l2-l1), l1+l2 679 auxrho(1:mesh) = & 680 pawset_%aewfc(1:mesh,ns1) * pawset_%aewfc(1:mesh,ns2) - & 681 pawset_%pswfc(1:mesh,ns1) * pawset_%pswfc(1:mesh,ns2) - & 682 pawset_%augfun(1:mesh,ns1,ns2,l3) 683 call hartree(l3,l1+l2+2,mesh,pawset_%grid,auxrho,auxpot(1,l3)) 684 end do 685 write (stdout,'(a,2i3,a,2i3)') " MULTIPOLO DI ",ns1,l1,":",ns2, l2 686 do ir=1,irc 687 if (r(ir) < 1.0_dp) ir0 = ir 688 end do 689 do ir=ir0,irc+30, 3 690 write (stdout,'(10f8.4)') r(ir),(auxpot(ir,l3), l3=0,l1+l2) 691 end do 692 end do 693 end do 694 return 695 END SUBROUTINE check_multipole 696 ! 697 !============================================================================ 698 ! PRIVATE ROUTINES !!! 699 !============================================================================ 700 ! 701 SUBROUTINE compute_charges (projsum_, chargeps_, charge1_, charge1ps_, & 702 pawset_, nwfc_, l_, j_, nspin_, spin_, oc_, pswfc_ , iflag, unit_) 703 USE io_global, ONLY : ionode 704 IMPLICIT NONE 705 REAL(dp), INTENT(OUT) :: projsum_(nwfsx,nwfsx,2) 706 REAL(dp), INTENT(OUT) :: chargeps_(ndmx,2) 707 REAL(dp), INTENT(OUT) :: charge1_(ndmx,2) 708 REAL(dp), INTENT(OUT) :: charge1ps_(ndmx,2) 709 TYPE(paw_t), INTENT(IN) :: pawset_ 710 INTEGER, INTENT(IN) :: nwfc_ 711 INTEGER, INTENT(IN) :: l_(nwfsx) 712 INTEGER, INTENT(IN) :: nspin_ 713 INTEGER, INTENT(IN) :: spin_(nwfsx) 714 REAL(dp), INTENT(IN) :: j_(nwfsx) 715 REAL(dp), INTENT(IN) :: oc_(nwfsx) 716 REAL(dp), INTENT(IN) :: pswfc_(ndmx,nwfsx) 717 INTEGER, OPTIONAL :: unit_, iflag 718 REAL(dp) :: augcharge(ndmx,2), chargetot 719 INTEGER :: i, n, iflag0 720 721 iflag0=0 722 if (present(iflag)) iflag0=iflag 723 CALL compute_sumwfc2(chargeps_,pawset_,nwfc_,pswfc_,oc_,spin_) 724 CALL compute_projsum(projsum_,pawset_,nwfc_,l_,j_,spin_,pswfc_,oc_) 725! WRITE (6200,'(20e20.10)') ((projsum_(ns,ns1,1),ns=1,ns1),ns1=1,pawset_%nwfc) 726 CALL compute_onecenter_charge(charge1ps_,pawset_,projsum_,nspin_,"PS") 727 CALL compute_onecenter_charge(charge1_ ,pawset_,projsum_,nspin_,"AE") 728 ! add augmentation charges 729 CALL compute_augcharge(augcharge,pawset_,projsum_,nspin_) 730 731 if (present(unit_).and.ionode) then 732 write(unit_,*) 733 write(unit_,*) "#" 734 do i=1,pawset_%grid%mesh 735 write (unit_,'(4f12.8)') pawset_%grid%r(i), augcharge(i,1), chargeps_(i,1), charge1ps_(i,1) 736 end do 737 end if 738 chargeps_ (1:pawset_%grid%mesh,1:nspin_) = chargeps_ (1:pawset_%grid%mesh,1:nspin_) & 739 + augcharge(1:pawset_%grid%mesh,1:nspin_) 740 charge1ps_(1:pawset_%grid%mesh,1:nspin_) = charge1ps_(1:pawset_%grid%mesh,1:nspin_) & 741 + augcharge(1:pawset_%grid%mesh,1:nspin_) 742! 743! If there are unbounded scattering states in the pseudopotential generation, 744! n1 and n~1 separately diverge. This makes the hartree and exchange and 745! correlation energies separately diverging. The cancellation becomes 746! very difficult numerically. Outside the sphere however the two charges 747! should be equal and opposite and we set them to zero. 748! Is this the right solution? 749! 750 do n=pawset_%irc+1,pawset_%grid%mesh 751 chargetot=chargeps_(n,1) 752 if (nspin_==2) chargetot=chargetot+chargeps_(n,2) 753 if (chargetot<1.d-11.or.iflag0==1) then 754 charge1_(n,1:nspin_)=0.0_DP 755 charge1ps_(n,1:nspin_)=0.0_DP 756 endif 757 enddo 758! do n=1,pawset_%grid%mesh 759! write(6,'(4f20.15)') pawset_%grid%r(n), chargeps_(n,1), charge1_(n,1), 760! charge1ps_(n,1) 761 762 END SUBROUTINE compute_charges 763 ! 764 !============================================================================ 765 ! 766 ! Compute the one-center energy and effective potential: 767 ! E = Eh[n_v] + Exc[n_v+n_c] - Int[veff*n_v], 768 ! veff = vh[n_v] + v_xc[n_v+n_c], 769 ! where n_v can be n~+n^ or n1 or n1~+n^ and n_c can be nc or n~c 770 ! 771 SUBROUTINE compute_onecenter_energy ( totenergy_, veff_, & 772 pawset_, vcharge_, nlcc_, ccharge_, nspin_, iint, vloc, energies_ , unit_) 773 USE funct, ONLY: dft_is_gradient 774 USE radial_grids, ONLY: hartree 775 USE io_global, ONLY : stdout, ionode 776 IMPLICIT NONE 777 REAL(dp), INTENT(OUT) :: totenergy_ ! H+XC+DC 778 REAL(dp), INTENT(OUT) :: veff_(ndmx,2) ! effective potential 779 TYPE(paw_t), INTENT(IN) :: pawset_ ! the PAW dataset 780 REAL(dp), INTENT(IN) :: vcharge_(ndmx,2) ! valence charge 781 LOGICAL, INTENT(IN) :: nlcc_ ! non-linear core correction 782 REAL(dp), INTENT(IN) :: ccharge_(ndmx) ! core charge 783 INTEGER, INTENT(IN) :: nspin_ ! 1 for LDA, 2 for LSDA 784 INTEGER, INTENT(IN) :: iint ! integrals up to iint 785 REAL(dp), INTENT(IN), OPTIONAL :: vloc(ndmx) ! 786 REAL(dp), INTENT(OUT), OPTIONAL :: energies_(5)! Etot,H,XC,DC terms 787 INTEGER, OPTIONAL :: unit_ 788 ! 789 REAL(dp), PARAMETER :: rho_eq_0(ndmx) = ZERO ! ccharge=0 when nlcc=.f. 790 ! 791 REAL(dp) :: & 792 eh, exc, edc, & ! hartree, xc and double counting energies 793 eloc, & ! local energy 794 rhovtot(ndmx), & ! total valence charge 795 aux(ndmx), & ! auxiliary to compute integrals 796 vh(ndmx), & ! hartree potential 797 vxc(ndmx,2), & ! exchange-correlation potential (LDA+GGA) 798 vgc(ndmx,2), & ! exchange-correlation potential (GGA only) 799 egc(ndmx), & ! exchange correlation energy density (GGA only) 800 rh(2), & ! valence charge at a given point without 4 pi r^2 801 rhc, & ! core charge at a given point without 4 pi r^2 802 vxcr(2) ! exchange-correlation potential at a given point 803 REAL(dp) :: & ! compatibility with metaGGA - not yet used 804 tau(ndmx) = ZERO, vtau(ndmx) = ZERO ! 805 ! 806 INTEGER :: i, is 807 INTEGER :: lsd 808 REAL(DP), EXTERNAL :: int_0_inf_dr 809#if defined __DEBUG_V_H_vs_SPHEROPOLE 810 REAL(DP) :: dummy_charge,aux1(ndmx),aux2(ndmx),res1,res2 811#endif 812 ! 813 ! Set up total valence charge 814 rhovtot(1:pawset_%grid%mesh) = vcharge_(1:pawset_%grid%mesh,1) 815 IF (nspin_==2) rhovtot(1:pawset_%grid%mesh) = rhovtot(1:pawset_%grid%mesh) + & 816 vcharge_(1:pawset_%grid%mesh,2) 817 ! 818 ! Hartree 819 CALL hartree(0,2,pawset_%grid%mesh,pawset_%grid,rhovtot,vh) 820 if (PRESENT(unit_).and.ionode) then 821 write (unit_,*) " " 822 write (unit_,*) "#" 823 do i=1,pawset_%grid%mesh 824 write (unit_,'(3f12.7)') pawset_%grid%r(i),rhovtot(i),vh(i) 825 end do 826 end if 827#if defined __DEBUG_V_H_vs_SPHEROPOLE 828 dummy_charge=int_0_inf_dr(rhovtot,pawset_%grid,pawset_%grid%mesh,2) 829 aux1(1:pawset_%grid%mesh) = FPI*pawset_%grid%r2(1:pawset_%grid%mesh)*vh(1:pawset_%grid%mesh) - & 830 FPI*dummy_charge * pawset_%grid%r(1:pawset_%grid%mesh) 831 aux2(1:pawset_%grid%mesh) = rhovtot(1:pawset_%grid%mesh)*pawset_%grid%r2(1:pawset_%grid%mesh) 832 res1 = int_0_inf_dr(aux1,pawset_%grid,pawset_%grid%mesh,1) 833 res2 = int_0_inf_dr(aux2,pawset_%grid,pawset_%grid%mesh,4) 834 WRITE (stdout,'(4(A,1e15.7))') ' INT rho', dummy_charge,' INT V_H', & 835 res1, ' INT r^2*rho', res2, ' ERR:', (1.d0- res1/ (-res2 * (2.d0*PI/3.d0))) 836#endif 837 vh(1:pawset_%grid%mesh) = e2 * vh(1:pawset_%grid%mesh) 838 aux(1:pawset_%grid%mesh) = vh(1:pawset_%grid%mesh) * rhovtot(1:pawset_%grid%mesh) 839 eh = HALF * int_0_inf_dr(aux,pawset_%grid,iint,2) 840 ! 841 ! Exhange-Correlation 842 rh=(/ZERO,ZERO/) 843 rhc=ZERO 844 lsd=nspin_-1 845 DO i=1,pawset_%grid%mesh 846 DO is=1,nspin_ 847 rh(is) = vcharge_(i,is)/pawset_%grid%r2(i)/FPI 848 ENDDO 849 IF (nlcc_) rhc = ccharge_(i)/pawset_%grid%r2(i)/FPI 850 CALL vxc_t(lsd,rh,rhc,exc,vxcr) 851 vxc(i,1:nspin_)=vxcr(1:nspin_) 852 IF (nlcc_) THEN 853 aux(i)=exc * (rhovtot(i)+ccharge_(i)) 854 ELSE 855 aux(i)=exc * rhovtot(i) 856 END IF 857 END DO 858 IF (dft_is_gradient()) THEN 859 IF (nlcc_) THEN 860 CALL vxcgc(ndmx,pawset_%grid%mesh,nspin_,pawset_%grid%r,& 861 pawset_%grid%r2,vcharge_,ccharge_,vgc,egc, & 862 tau, vtau, 1) 863 ELSE 864 CALL vxcgc(ndmx,pawset_%grid%mesh,nspin_,pawset_%grid%r,& 865 pawset_%grid%r2,vcharge_,rho_eq_0,vgc,egc, & 866 tau, vtau, 1) 867 END IF 868 vxc(1:pawset_%grid%mesh,1:nspin_) = vxc(1:pawset_%grid%mesh,1:nspin_) + & 869 vgc(1:pawset_%grid%mesh,1:nspin_) 870 aux(1:pawset_%grid%mesh) = aux(1:pawset_%grid%mesh) + & 871 egc(1:pawset_%grid%mesh) * pawset_%grid%r2(1:pawset_%grid%mesh) * FPI 872 END IF 873 exc = int_0_inf_dr(aux,pawset_%grid,iint,2) 874 ! 875 ! Double counting 876 edc=ZERO 877 DO is=1,nspin_ 878 veff_(1:pawset_%grid%mesh,is)=vxc(1:pawset_%grid%mesh,is)+vh(1:pawset_%grid%mesh) 879 aux(1:pawset_%grid%mesh)=veff_(1:pawset_%grid%mesh,is)*vcharge_(1:pawset_%grid%mesh,is) 880 edc=edc+int_0_inf_dr(aux,pawset_%grid,iint,2) 881 END DO 882 ! 883 eloc=ZERO 884 ! 885 IF (present(vloc)) THEN 886 DO is=1,nspin_ 887 aux(1:pawset_%grid%mesh)=vloc(1:pawset_%grid%mesh) & 888 *vcharge_(1:pawset_%grid%mesh,is) 889 eloc=eloc+int_0_inf_dr(aux,pawset_%grid,iint,2) 890 ENDDO 891 ENDIF 892 ! 893 ! Total 894 totenergy_ = eh + exc - edc 895 ! 896 IF (PRESENT(energies_)) THEN 897 energies_(1)=totenergy_ 898 energies_(2)=eh 899 energies_(3)=exc 900 energies_(4)=edc 901 energies_(5)=eloc 902 END IF 903 ! 904 END SUBROUTINE compute_onecenter_energy 905 ! 906 !============================================================================ 907 ! 908 ! Compute NL 'D' coefficients = D^ + D1 - D1~ 909 ! 910 SUBROUTINE compute_nonlocal_coeff(ddd_, pawset_, nspin_, veffps_, veff1_, veff1ps_) 911 IMPLICIT NONE 912 REAL(dp), INTENT(OUT) :: ddd_(nwfsx,nwfsx,2) 913 TYPE(paw_t), INTENT(IN) :: pawset_ 914 INTEGER, INTENT(IN) :: nspin_ 915 REAL(dp), INTENT(IN) :: veffps_(ndmx,2) 916 REAL(dp), INTENT(IN) :: veff1_(ndmx,2) 917 REAL(dp), INTENT(IN) :: veff1ps_(ndmx,2) 918 INTEGER :: is, ns, ns1 919 REAL(dp) :: aux(ndmx), dd 920 REAL(DP), EXTERNAL :: int_0_inf_dr 921! REAL(dp):: dddd(nwfsx,nwfsx,3) = 0.d0 922 923 ! 924 ! D^ = Int Q*v~ 925 ! D1-D1~ = kindiff + Int[ae*v1*ae - ps*v1~*ps - Q*v1~] 926 ddd_(:,:,:)=ZERO 927 DO is=1,nspin_ 928 DO ns=1,pawset_%nwfc 929 DO ns1=1,ns 930 IF (pawset_%l(ns)==pawset_%l(ns1).and.& 931 ABS(pawset_%jj(ns)-pawset_%jj(ns1))<1.d-8) THEN 932 ! Int[Q*v~] 933 aux(1:pawset_%grid%mesh) = & 934 pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) * & 935 veffps_(1:pawset_%grid%mesh,is) 936 dd = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) 937! dddd(ns,ns1,1) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) 938 ! Int[ae*v1*ae] 939 aux(1:pawset_%grid%mesh) = & 940 pawset_%aewfc(1:pawset_%grid%mesh,ns ) * & 941 pawset_%aewfc(1:pawset_%grid%mesh,ns1) * & 942 veff1_(1:pawset_%grid%mesh,is) 943 IF (pawset_%rel==2) & 944 aux(1:pawset_%irc) = aux(1:pawset_%irc) + & 945 pawset_%aewfc_rel(1:pawset_%irc,ns ) * & 946 pawset_%aewfc_rel(1:pawset_%irc,ns1) * & 947 veff1_(1:pawset_%irc,is) 948 dd = dd + & 949 int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) 950! dddd(ns,ns1,2) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) 951 ! Int[ps*v1~*ps + aufun*v1~] 952 aux(1:pawset_%grid%mesh) = & 953 ( pawset_%pswfc(1:pawset_%grid%mesh,ns ) * & 954 pawset_%pswfc(1:pawset_%grid%mesh,ns1) + & 955 pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) ) * & 956 veff1ps_(1:pawset_%grid%mesh,is) 957 dd = dd - & 958 int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) 959! dddd(ns,ns1,3) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) 960 ! collect 961 ddd_(ns,ns1,is) = pawset_%kdiff(ns,ns1) + dd 962 ddd_(ns1,ns,is) = ddd_(ns,ns1,is) 963 END IF 964 END DO 965 END DO 966 END DO 967! write(*,*) 'deeq', pawset_%irc 968! write(*,'(4f13.8)') dddd(1:4,1:4,1) 969! write(*,*) 'ddd ae' 970! write(*,'(4f13.8)') dddd(1:4,1:4,2) 971! write(*,*) 'ddd ps' 972! write(*,'(4f13.8)') dddd(1:4,1:4,3) 973 974 END SUBROUTINE compute_nonlocal_coeff 975 ! 976 !============================================================================ 977 ! 978 ! 'D_ion' coefficients = D1 - D1~ 979 ! 980 SUBROUTINE compute_nonlocal_coeff_ion(ddd_, pawset_) 981 IMPLICIT NONE 982 REAL(dp), INTENT(OUT) :: ddd_(nwfsx,nwfsx) 983 TYPE(paw_t), INTENT(IN) :: pawset_ 984 INTEGER :: ns, ns1 985 REAL(dp) :: aux(ndmx), dd 986 REAL(DP), EXTERNAL :: int_0_inf_dr 987 ! 988 ! D^ = Int Q*v~ 989 ! D1-D1~ = kindiff + Int[ae*v1*ae - ps*v1~*ps - Q*v1~] 990! write(666,"(4f12.6)") (pawset_%grid%r(ns), veffps_(ns,1), veff1_(ns,1), veff1ps_(ns,1), ns=1,pawset_%grid%mesh) 991! write(667,"(4f12.6)") (pawset_%grid%r(ns), veffps_(ns,2), veff1_(ns,2), veff1ps_(ns,2), ns=1,pawset_%grid%mesh) 992 ddd_(:,:)=ZERO 993 DO ns=1,pawset_%nwfc 994 DO ns1=1,ns 995 IF (pawset_%l(ns)==pawset_%l(ns1).and. & 996 ABS(pawset_%jj(ns)-pawset_%jj(ns1))<1.d-8 ) THEN 997 ! Int[ae*v1*ae] 998 aux(1:pawset_%grid%mesh) = & 999 pawset_%aewfc(1:pawset_%grid%mesh,ns ) * & 1000 pawset_%aewfc(1:pawset_%grid%mesh,ns1) * & 1001 pawset_%aeloc(1:pawset_%grid%mesh) 1002 IF (pawset_%rel==2) & 1003 aux(1:pawset_%irc) = aux(1:pawset_%irc) + & 1004 pawset_%aewfc_rel(1:pawset_%irc,ns ) * & 1005 pawset_%aewfc_rel(1:pawset_%irc,ns1) * & 1006 pawset_%aeloc(1:pawset_%irc) 1007 dd = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) 1008 ! Int[ps*v1~*ps + Q*v1~] 1009 aux(1:pawset_%grid%mesh) = & 1010 ( pawset_%pswfc(1:pawset_%grid%mesh,ns ) * & 1011 pawset_%pswfc(1:pawset_%grid%mesh,ns1) + & 1012 pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) ) * & 1013 pawset_%psloc(1:pawset_%grid%mesh) 1014 dd = dd - & 1015 int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2) 1016 ! 1017 ddd_(ns,ns1) = pawset_%kdiff(ns,ns1) + dd 1018 ddd_(ns1,ns)=ddd_(ns,ns1) 1019 END IF 1020 END DO 1021 END DO 1022 END SUBROUTINE compute_nonlocal_coeff_ion 1023 ! 1024 !============================================================================ 1025 ! 1026 ! Write PAW dataset wfc and potentials on files 1027 ! 1028 SUBROUTINE human_write_paw(pawset_) 1029 USE io_global, ONLY : ionode 1030 IMPLICIT NONE 1031 TYPE(paw_t), INTENT(In) :: pawset_ 1032 INTEGER :: n,ns 1033 IF (.not.ionode) return 1034 DO ns=1,pawset_%nwfc 1035 WRITE (5000+ns,'(A)') "# r AEwfc PSwfc projector" 1036 DO n=1,pawset_%grid%mesh 1037 WRITE (5000+ns,'(4f20.12)') pawset_%grid%r(n), pawset_%aewfc(n,ns), pawset_%pswfc(n,ns),pawset_%proj(n,ns) 1038 END DO 1039 END DO 1040 ! 1041 WRITE (6000,'(A)') "# r AEpot PSpot" 1042 DO n=1,pawset_%grid%mesh 1043 WRITE (6000,'(3f20.8)') pawset_%grid%r(n), pawset_%aeloc(n), pawset_%psloc(n) 1044 END DO 1045 END SUBROUTINE human_write_paw 1046 ! 1047 !============================================================================ 1048 ! LOWER-LEVEL ROUTINES !!! 1049 !============================================================================ 1050 ! 1051 ! Compute Sum_i oc_i * wfc_i^2 1052 ! 1053 SUBROUTINE compute_sumwfc2(charge_, pawset_, nwfc_, wfc_, oc_, spin_) 1054 IMPLICIT NONE 1055 REAL(dp), INTENT(OUT) :: charge_(ndmx,2) 1056 TYPE(paw_t), INTENT(IN) :: pawset_ 1057 INTEGER, INTENT(IN) :: nwfc_ 1058 REAL(dp), INTENT(IN) :: wfc_(ndmx,nwfsx) 1059 REAL(dp), INTENT(IN) :: oc_(nwfsx) 1060 INTEGER, INTENT(IN) :: spin_(nwfsx) 1061 INTEGER :: nf 1062 charge_(1:pawset_%grid%mesh,:)=ZERO 1063 DO nf=1,nwfc_ 1064 IF (oc_(nf)>ZERO) charge_(1:pawset_%grid%mesh,spin_(nf)) = & 1065 charge_(1:pawset_%grid%mesh,spin_(nf)) + oc_(nf) * & 1066 wfc_(1:pawset_%grid%mesh,nf) * wfc_(1:pawset_%grid%mesh,nf) 1067 END DO 1068 END SUBROUTINE compute_sumwfc2 1069 ! 1070 !============================================================================ 1071 ! 1072 ! Compute Sum_n oc_n <pswfc_n|proj_i> <proj_j|pswfc_n> 1073 ! 1074 SUBROUTINE compute_projsum (projsum_, pawset_, nwfc_, l_, j_, spin_, pswfc_, oc_) 1075 REAL(dp), INTENT(OUT) :: projsum_(nwfsx,nwfsx,2) 1076 TYPE(paw_t), INTENT(IN) :: pawset_ 1077 INTEGER, INTENT(IN) :: nwfc_ 1078 INTEGER, INTENT(IN) :: l_(nwfsx) 1079 INTEGER, INTENT(IN) :: spin_(nwfsx) 1080 REAL(dp), INTENT(IN) :: j_(nwfsx) 1081 REAL(dp), INTENT(IN) :: pswfc_(ndmx,nwfsx) 1082 REAL(dp), INTENT(IN) :: oc_(nwfsx) 1083 REAL(dp) :: proj_dot_wfc(nwfsx,nwfsx), aux(ndmx) 1084 INTEGER :: ns, ns1, nf, nr, is, nst 1085 REAL(DP), EXTERNAL :: int_0_inf_dr 1086 ! Compute <projector|wavefunction> 1087 DO ns=1,pawset_%nwfc 1088 DO nf=1,nwfc_ 1089 IF (pawset_%l(ns)==l_(nf).AND.pawset_%jj(ns)==j_(nf)) THEN 1090 DO nr=1,pawset_%grid%mesh 1091 aux(nr)=pawset_%proj(nr,ns)*pswfc_(nr,nf) 1092 END DO 1093 nst=(l_(nf)+1)*2 1094 proj_dot_wfc(ns,nf)=int_0_inf_dr(aux,pawset_%grid,pawset_%ikk(ns),nst) 1095 ELSE 1096 proj_dot_wfc(ns,nf)=ZERO 1097 END IF 1098 END DO 1099 END DO 1100 ! Do the sum over the wavefunctions 1101 projsum_(:,:,:)=ZERO 1102 DO ns=1,pawset_%nwfc 1103 DO ns1=1,ns 1104 DO nf=1,nwfc_ 1105 is=spin_(nf) 1106 IF (oc_(nf)>ZERO) THEN 1107 projsum_(ns,ns1,is) = projsum_(ns,ns1,is) + oc_(nf) * & 1108 proj_dot_wfc(ns,nf) * proj_dot_wfc(ns1,nf) 1109 END IF 1110 END DO 1111 projsum_(ns1,ns,:)=projsum_(ns,ns1,:) 1112 END DO 1113 END DO 1114 END SUBROUTINE compute_projsum 1115 ! 1116 !============================================================================ 1117 ! 1118 ! 1119 ! Compute augmentation charge as Sum_ij W_ij * Q_ij 1120 ! 1121 SUBROUTINE compute_augcharge(augcharge_, pawset_, projsum_, nspin_) 1122 IMPLICIT NONE 1123 REAL(dp), INTENT(OUT) :: augcharge_(ndmx,2) 1124 TYPE(paw_t), INTENT(IN) :: pawset_ 1125 REAL(dp), INTENT(IN) :: projsum_(nwfsx,nwfsx,2) 1126 INTEGER, INTENT(IN) :: nspin_ 1127 INTEGER :: ns, ns1, is 1128 REAL(dp) :: factor 1129 augcharge_=ZERO 1130 DO is=1,nspin_ 1131 DO ns=1,pawset_%nwfc 1132 DO ns1=1,ns 1133 ! multiply by TWO off-diagonal terms 1134 factor=TWO 1135 IF (ns1==ns) factor=ONE 1136 ! 1137 ! NB: in a spherically averaged system only the l=0 component 1138 ! of the augmentation functions is present 1139 ! 1140 augcharge_(1:pawset_%grid%mesh,is) = augcharge_(1:pawset_%grid%mesh,is) + & 1141 factor * projsum_(ns,ns1,is) * & 1142 pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) 1143 END DO 1144 END DO 1145 END DO 1146 END SUBROUTINE compute_augcharge 1147 ! 1148 !============================================================================ 1149 ! 1150 ! Compute n1 and n1~, as Sum_ij W_ij wfc_i(r)*wfc_j(r) 1151 ! 1152 SUBROUTINE compute_onecenter_charge(charge1_, pawset_, projsum_, nspin_, which_wfc) 1153 IMPLICIT NONE 1154 REAL(dp), INTENT(OUT) :: charge1_(ndmx,2) 1155 TYPE(paw_t), INTENT(IN) :: pawset_ 1156 REAL(dp), INTENT(IN) :: projsum_(nwfsx,nwfsx,2) 1157 INTEGER, INTENT(IN) :: nspin_ 1158 CHARACTER(LEN=2),INTENT(IN) :: which_wfc 1159 INTEGER :: ns, ns1, is 1160 REAL(dp) :: factor 1161 charge1_=ZERO 1162 DO is=1,nspin_ 1163 DO ns=1,pawset_%nwfc 1164 DO ns1=1,ns 1165 ! multiply times TWO off-diagonal terms 1166 IF (ns1==ns) THEN 1167 factor=ONE 1168 ELSE 1169 factor=TWO 1170 END IF 1171 SELECT CASE (which_wfc) 1172 CASE ("AE") 1173 charge1_(1:pawset_%grid%mesh,is) = charge1_(1:pawset_%grid%mesh,is) + factor * & 1174 projsum_(ns,ns1,is) * pawset_%aewfc(1:pawset_%grid%mesh,ns) * & 1175 pawset_%aewfc(1:pawset_%grid%mesh,ns1) 1176 IF (pawset_%rel==2) & 1177 charge1_(1:pawset_%irc,is) = charge1_(1:pawset_%irc,is) + factor * & 1178 projsum_(ns,ns1,is) * pawset_%aewfc_rel(1:pawset_%irc,ns) * & 1179 pawset_%aewfc_rel(1:pawset_%irc,ns1) 1180 CASE ("PS") 1181 charge1_(1:pawset_%grid%mesh,is) = charge1_(1:pawset_%grid%mesh,is) + factor * & 1182 projsum_(ns,ns1,is) * pawset_%pswfc(1:pawset_%grid%mesh,ns) * & 1183 pawset_%pswfc(1:pawset_%grid%mesh,ns1) 1184 CASE DEFAULT 1185 call errore ('compute_onecenter_charge','specify AE or PS wavefunctions',1) 1186 END SELECT 1187 END DO 1188 END DO 1189 END DO 1190 END SUBROUTINE compute_onecenter_charge 1191! 1192!-------------------------------------------------------------------------- 1193SUBROUTINE find_bes_qi(qc,rmatch,lam,ncn,iok) 1194 !-------------------------------------------------------------------------- 1195 ! 1196 ! This routine finds two values of q such that the 1197 ! functions f_l have a derivative equal to 0 at rmatch 1198 ! 1199 IMPLICIT NONE 1200 1201 INTEGER, INTENT(IN) :: & 1202 lam, & ! input: the angular momentum 1203 ncn ! input: the number of qi to compute 1204 INTEGER, INTENT(INOUT) :: & 1205 iok ! output: if 0 the calculation in this routine is ok 1206 1207 REAL (dp), INTENT(OUT) :: & 1208 qc(ncn) ! output: the values of qi 1209 REAL (dp), INTENT(IN) :: rmatch 1210 1211 REAL (dp) :: & 1212 zeroderjl (2,7) ! first two zeros of the first derivative of 1213 ! spherical Bessel function j_l for l = 0,...,6 1214 1215 INTEGER :: & 1216 nc ! counter on the q found 1217 1218 data zeroderjl / 0.0_dp, 4.4934094579614_dp, & 1219 2.0815759780862_dp, 5.9403699890844_dp, & 1220 3.3420936578747_dp, 7.2899322987026_dp, & 1221 4.5140996477983_dp, 8.5837549433127_dp, & 1222 5.6467036213923_dp, 9.8404460168549_dp, & 1223 6.7564563311363_dp, 11.0702068269176_dp, & 1224 7.8510776799611_dp, 12.2793339053177_dp / 1225 iok=0 1226 IF (ncn.gt.2) & 1227 CALL errore('find_aug_qi','ncn is too large',1) 1228 1229 IF (lam.gt.6) & 1230 CALL errore('find_aug_qi','l not programmed',1) 1231 ! 1232 ! fix deltaq and the maximum step number 1233 ! 1234 DO nc = 1, ncn 1235 ! 1236 qc(nc) = zeroderjl (nc, lam + 1) / rmatch 1237 ! 1238 ENDDO 1239 RETURN 1240END SUBROUTINE find_bes_qi 1241 ! 1242 !============================================================================ 1243 ! 1244END MODULE atomic_paw 1245 1246