1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 subroutine bsc_cellxc( irel, ider, cell, 9 $ NMesh, NSpan, maxp, mtype, 10 & XMesh, nspin, Dens, EX, EC, DX, DC, VXC, 11 & DVXCDN, stressl ) 12 13C ******************************************************************* 14C Finds total exchange-correlation energy and potential in a 15C periodic cell. 16C This version implements the Local (spin) Density Approximation and 17C the Generalized-Gradient-Aproximation with the 'explicit mesh 18C functional' approach of White & Bird, PRB 50, 4954 (1994). 19C Gradients are 'defined' by numerical derivatives, using 2*NN+1 mesh 20C points, where NN is a parameter defined below 21C Ref: L.C.Balbas et al, PRB 64, 165110 (2001) 22C Wrtten by J.M.Soler using algorithms developed by 23C L.C.Balbas, J.L.Martins and J.M.Soler, Dec.1996 - Aug.1997 24C Parallel version written by J.Gale. June 1999. 25C Argument DVXCDN added by J.Junquera. November 2000. 26C Adapted for multiple functionals in the same run by J.Gale 2005 27C ************************* INPUT *********************************** 28C integer irel : Relativistic exchange? (0=>no, 1=>yes) 29C integer ider : Return dVxc/drho in DVXCDN? 30C 0=>no, 1=>yes (available only for LDA) 31C real*8 cell(3,3) : Unit cell vectors cell(ixyz,ivector) 32C integer NMesh(3) : Number of mesh divisions of each vector 33C integer NSpan(3) : Physical dimensions of arrays XMesh, Dens and 34C VXC (or memory span between array elements) 35C See usage section for more information 36C integer maxp : Physical dimension of XMesh, Dens, and VXC 37C integer mtype : Mesh type: 38C 0 => Uniform mesh 39C 1 => Adaptive mesh, given in cartesian coord 40C 2 => Adaptive mesh, given in cell-vector coord 41C real XMesh(3,maxp): Mesh point coordinates (not used if mtype=0) 42C When mtype=2, cartesian coordinates are 43C Xcart(ix,im) = Sum_iv(cell(ix,iv)*XMesh(iv,ip)) 44C Notice single precision in this version 45C integer nspin : nspin=1 => unpolarized; nspin=2 => polarized; 46C nspin=4 => non-collinear polarization 47C real Dens(maxp,nspin) : Total (nspin=1) or spin (nspin=2) electron 48C density at mesh points, ordered as 49C IP = I1+NSpan(1)*((I2-1)+NSpan(2)*(I3-1)), 50C with I1=1,...,NMesh(1), etc 51C For non-collinear polarization, the density 52C matrix is given by: Dens(1)=D11, Dens(2)=D22, 53C Dens(3)=Real(D12), Dens(4)=Im(D12) 54C Notice single precision in this version 55C ************************* OUTPUT ********************************** 56C real*8 EX : Total exchange energy 57C real*8 EC : Total correlation energy 58C real*8 DX : IntegralOf( rho * (eps_x - v_x) ) 59C real*8 DC : IntegralOf( rho * (eps_c - v_c) ) 60C real VXC(maxp,nspin) : (Spin) exch-corr potential 61C Notice single precision in this version 62C real DVXCDN(maxp,nspin,nspin) : Derivatives of exchange-correlation 63C potential respect the charge density 64C Not used unless ider=1. Available only for LDA 65C real*8 stressl(3,3) : xc contribution to the stress tensor, 66C assuming constant density (not charge), 67C i.e. r->r' => rho'(r') = rho(r) 68C For plane-wave and grid (finite diff) basis 69C sets, density rescaling gives an extra term 70C (not included) (DX+DC-EX-EC)/cell_volume for 71C the diagonal elements of stress. For other 72C basis sets, the extra term is, in general: 73C IntegralOf(v_xc * d_rho/d_strain) / cell_vol 74C ************************ UNITS ************************************ 75C Distances in atomic units (Bohr). 76C Densities in atomic units (electrons/Bohr**3) 77C Energy unit depending of parameter EUnit below 78C Stress in EUnit/Bohr**3 79C ************************ USAGE ************************************ 80C Typical calls for different array dimensions: 81C parameter ( maxp = 1000000 ) 82C DIMENSION NMesh(3), Dens(maxp,2), VXC(maxp,2) 83C Find cell vectors 84C Find density at N1*N2*N3 mesh points (less than maxp) and place 85C them consecutively in array Dens 86C NMesh(1) = N1 87C NMesh(2) = N2 88C NMesh(3) = N3 89C CALL cellXC( 0, 0, cell, NMesh, NMesh, maxp, 0, XMesh, 90C . 2, Dens, EX, EC, DX, DC, VXC, DVXCDN, stress ) 91C Or alternatively: 92C parameter ( M1=100, M2=100, M3=100 ) 93C DIMENSION NMesh(3), NSpan(3), Dens(M1,M2,M3,2), VXC(M1,M2,M3,2) 94C DATA NSpan / M1, M2, M3 / 95C Find cell vectors 96C Find Dens at N1*N2*N3 mesh points 97C NMesh(1) = N1 98C NMesh(2) = N2 99C NMesh(3) = N3 100C CALL cellXC( 0, 0, cell, NMesh, NSpan, M1*M2*M3, 0, XMesh, 101C . 2, Dens, EX, EC, DX, DC, VXC, DVXCDN, stress ) 102C ********* BEHAVIOUR *********************************************** 103C - Notice that XMesh is not used if mtype=0, and DVXCDN is not 104C used if ider=0. This means that you do not even need to dimension 105C them in the calling program. See usage section for calling examples. 106C - If FNCTNL='LDA', Dens and VXC may be the same physical array. 107C - Stops and prints a warning if arguments functl or XCauth are not 108C one of the allowed possibilities. 109C - Since the exchange and correlation part is usually a small fraction 110C of a typical electronic structure calculation, this routine has 111C been coded with emphasis on simplicity and functionality, not in 112C eficiency. The parallel version was written by J.Gale. 113C ********* ROUTINES CALLED ***************************************** 114C GGAXC, LDAXC, RECLAT, TIMER, VOLCEL 115C ******************************************************************* 116 117 use precision, only : dp, grid_p 118 use sys, only : die 119 use alloc, only : re_alloc, de_alloc 120 use mesh, only : nsm, meshLim 121 use siestaxc, only : ggaxc, ldaxc 122 use siestaxc, only : getXC 123#ifdef MPI 124C Modules 125 use mesh, only : nmeshg 126 use parallel, only : Node, Nodes, ProcessorY 127 use parallelsubs, only : HowManyMeshPerNode 128 use mpi_siesta 129 use moreMeshSubs, only : distExtMeshData, gathExtMeshData 130#endif 131 132 implicit none 133 134C Argument types and dimensions 135 integer, intent(in) :: ider 136 integer, intent(in) :: irel 137 integer, intent(in) :: maxp 138 integer, intent(in) :: mtype 139 integer, intent(in) :: NMesh(3) 140 integer, intent(in) :: NSpan(3) 141 integer, intent(in) :: nspin 142 real(dp), intent(in) :: cell(3,3) 143 real(dp), intent(out) :: DC 144 real(dp), intent(out) :: DX 145 real(dp), intent(out) :: EC 146 real(dp), intent(out) :: EX 147 real(dp), intent(inout) :: stressl(3,3) 148 149C If you change next line, please change also the argument 150C explanations in the header comments 151!!! Pending (AG) 152* real(dp) 153 real(grid_p), intent(in) :: Dens(maxp,nspin) 154 real(grid_p), intent(out) :: DVXCDN(maxp,nspin,nspin) 155 real(grid_p), intent(out) :: VXC(maxp,nspin) 156 real(grid_p), intent(in) :: XMesh(3,*) 157 158C Fix the order of the numerical derivatives 159C NN is the number of points used in each coordinate and direction, 160C i.e. a total of 6*NN neighbour points is used to find the gradients 161 integer, parameter :: NN = 3 162 163C Fix energy unit: EUnit=1.0 => Hartrees, 164C EUnit=0.5 => Rydbergs, 165C EUnit=0.03674903 => eV 166 real(dp), parameter :: EUnit = 0.5_dp 167 168C Fix switch to skip points with zero density 169 logical, parameter :: skip0 = .true. 170 real(dp), parameter :: denmin = 1.0e-15_dp 171 172C Parameter mspin must be equal or larger than nspin 173 integer, parameter :: mspin = 4 174 175#ifdef MPI 176C MPI related variables 177 integer :: IOut, INN, JPNN(3,-NN:NN), MPIerror 178#endif 179 180C Local variables and arrays 181 logical GGA, GGAfunctl 182 integer I1, I2, I3, IC, IN, IP, IS, IX, 183 & J1, J2, J3, JN, JP(3,-NN:NN), JS, JX, 184 & KS, NPG, nf 185 real(dp) D(mspin), DECDD(mspin), DECDGD(3,mspin), 186 & dentot, DEXDD(mspin), DEXDGD(3,mspin), 187 & DGDM(-NN:NN), DGIDFJ(3,3,-NN:NN), 188 & DMDX(3,3), DVol, DXDM(3,3), 189 & DVCDN(mspin*mspin), DVXDN(mspin*mspin), 190 & EPSC, EPSX, F1, F2, GD(3,mspin), 191 & VOLCEL, volume, stress(3,3) 192 external RECLAT, VOLCEL 193 194 integer BS(3), iDistr 195 real(grid_p), pointer :: bdensX(:,:,:), bdensY(:,:,:), 196 & bdensZ(:,:,:), bvxcX(:,:,:), 197 & bvxcY(:,:,:), bvxcZ(:,:,:) 198 199 integer :: nXCfunc 200 character(len=20) :: XCauth(10), XCfunc(10) 201 real(dp) :: XCweightX(10), XCweightC(10) 202 203#ifdef DEBUG 204 call write_debug( ' PRE cellxc' ) 205#endif 206 207C Start time counter (intended only for debugging and development) 208#ifdef _TRACE_ 209 call MPI_Barrier( MPI_Comm_World, MPIerror ) 210 call MPItrace_event( 1000, 3 ) 211#endif 212 call timer( 'cellXC', 1 ) 213 214 call getXC( nXCfunc, XCfunc, XCauth, XCweightX, XCweightC) 215 216 GGA = .false. 217 do nf = 1,nXCfunc 218 if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga' ) then 219 GGA = .true. 220 endif 221 if ( XCfunc(nf).eq.'VDW' .or. XCfunc(nf).eq.'vdw') then 222 call die("BSC's cellxc cannot handle VDW functionals") 223 endif 224 enddo 225 226 if (ider.ne.0 .and. GGA) then 227 call die('BSC cellXC: ider=1 available only for LDA') 228 endif 229 230C Check value of mspin 231 if (mspin.lt.nspin) then 232 write(6,*) 'cellXC: parameter mspin must be at least ', nspin 233 call die() 234 endif 235 236 BS(1) = (meshLim(2,1)-meshLim(1,1)+1)*NSM 237 BS(2) = (meshLim(2,2)-meshLim(1,2)+1)*NSM 238 BS(3) = (meshLim(2,3)-meshLim(1,3)+1)*NSM 239#ifdef MPI 240C If GGA and the number of processors is greater than 1 then we need 241C to exchange border densities so that the density derivatives can 242C be calculated. 243 if (GGA.and.Nodes.gt.1) then 244 if (NN.gt.BS(1).or.NN.gt.BS(2).or.NN.gt.BS(3)) then 245 if (Node.eq.0) then 246 write(6,'('' ERROR - number of fine mesh points per '', 247 & ''Node must be greater than finite difference order '')') 248 endif 249 call die() 250 endif 251 252 iDistr = 3 253 if (BS(1).ne.NMeshG(1)) then 254 nullify(bdensX) 255 call re_alloc( bdensX, 1, BS(2)*BS(3), 1, 2*NN, 256 & 1, nspin, 'bdensX', 'cellxc' ) 257 nullify(bvxcX) 258 call re_alloc( bvxcX, 1, BS(2)*BS(3), 1, 2*NN, 259 & 1, nspin, 'bvxcX', 'cellxc' ) 260 call distExtMeshData( iDistr, 1, BS(2)*BS(3), NSM, NN, NSPIN, 261 & maxp, NMeshG, dens, bdensX ) 262 endif 263 264 if (BS(2).ne.NMeshG(2)) then 265 nullify(bdensY) 266 call re_alloc( bdensY, 1, BS(1)*BS(3), 1, 2*NN, 267 & 1, nspin, 'bdensY', 'cellxc' ) 268 nullify(bvxcY) 269 call re_alloc( bvxcY, 1, BS(1)*BS(3), 1, 2*NN, 270 & 1, nspin, 'bvxcY', 'cellxc' ) 271 call distExtMeshData( iDistr, 2, BS(1)*BS(3), NSM, NN, NSPIN, 272 & maxp, NMeshG, dens, bdensY ) 273 endif 274 275 if (BS(3).ne.NMeshG(3)) then 276 nullify(bdensZ) 277 call re_alloc( bdensZ, 1, BS(1)*BS(2), 1, 2*NN, 278 & 1, nspin, 'bdensZ', 'cellxc' ) 279 nullify(bvxcZ) 280 call re_alloc( bvxcZ, 1, BS(1)*BS(2), 1, 2*NN, 281 & 1, nspin, 'bvxcZ', 'cellxc' ) 282 call distExtMeshData( iDistr, 3, BS(1)*BS(2), NSM, NN, NSPIN, 283 & maxp, NMeshG, dens, bdensZ ) 284 endif 285 endif 286#endif 287 288C Find weights of numerical derivation from Lagrange interp. formula 289 do IN = -NN,NN 290 F1 = 1.0_dp 291 F2 = 1.0_dp 292 do JN = -NN,NN 293 if (JN.ne.IN .and. JN.ne.0) F1 = F1 * (0 - JN) 294 if (JN.ne.IN) F2 = F2 * (IN - JN) 295 enddo 296 DGDM(IN) = F1 / F2 297 enddo 298 DGDM(0) = 0.0_dp 299 300C Find total number of mesh points 301#ifdef MPI 302 NPG = NMeshG(1) * NMeshG(2) * NMeshG(3) 303#else 304 NPG = NMesh(1) * NMesh(2) * NMesh(3) 305#endif 306 307C Find Jacobian matrix dx/dmesh for uniform mesh 308 if (mtype.eq.0) then 309 310C Find mesh cell volume 311 DVol = VOLCEL( cell ) / DBLE(NPG) 312 313 if (GGA) then 314 315C Find mesh unit vectors and reciprocal mesh vectors 316 do IC = 1,3 317 do IX = 1,3 318#ifdef MPI 319 DXDM(IX,IC) = cell(IX,IC) / NMeshG(IC) 320#else 321 DXDM(IX,IC) = cell(IX,IC) / NMesh(IC) 322#endif 323 enddo 324 enddo 325 call RECLAT( DXDM, DMDX, 0 ) 326 327C Find the weights for the derivative d(gradF(i))/d(F(j)) of 328C the gradient at point i with respect to the value at point j 329 do IN = -NN,NN 330 do IC = 1,3 331 do IX = 1,3 332 DGIDFJ(IX,IC,IN) = DMDX(IX,IC) * DGDM(IN) 333 enddo 334 enddo 335 enddo 336 337 endif 338 339 endif 340 341C Initialize output 342 EX = 0.0_dp 343 EC = 0.0_dp 344 DX = 0.0_dp 345 DC = 0.0_dp 346 VXC(1:maxp,1:nspin) = 0.0_grid_p 347 if (ider.eq.1) then 348 DVXCDN(1:maxp,1:nspin,1:nspin) = 0.0_grid_p 349 endif 350 351#ifdef MPI 352C Initialise buffer regions of Vxc 353 if (GGA.and.Nodes.gt.1) then 354 if (BS(1).ne.NMeshG(1)) bvxcX = 0.0_grid_p 355 if (BS(2).ne.NMeshG(2)) bvxcY = 0.0_grid_p 356 if (BS(3).ne.NMeshG(3)) bvxcZ = 0.0_grid_p 357 endif 358#endif 359 stress(1:3,1:3) = 0.0_dp 360 361C Loop on mesh points 362 do I3 = 0,NMesh(3)-1 363 do I2 = 0,NMesh(2)-1 364 do I1 = 0,NMesh(1)-1 365 366C Find mesh index of this point 367 IP = 1 + I1 + NSpan(1) * I2 + NSpan(1) * NSpan(2) * I3 368 369C Skip point if density=0 370 if (skip0) then 371 dentot = 0.0_dp 372 do IS = 1,MIN(nspin,2) 373 dentot = dentot + MAX(0.0_grid_p,Dens(IP,IS)) 374 enddo 375 if (dentot .lt. denmin) then 376 do IS = 1,nspin 377 VXC(IP,IS) = 0.0_grid_p 378 enddo 379 goto 210 380 endif 381 endif 382 383C Find mesh indexes of neighbour points 384C Note : a negative index indicates a point from the buffer region 385 if (GGA .or. mtype.ne.0) then 386 387C X-direction 388#ifdef MPI 389 if (NMesh(1).eq.NMeshG(1)) then 390#endif 391 do IN = -NN,NN 392 J1 = MOD( I1+IN+100*BS(1), BS(1) ) 393 JP(1,IN) = 1+J1+BS(1)*I2+BS(1)*BS(2)*I3 394 enddo 395#ifdef MPI 396 else 397 do IN = -NN,NN 398 J1 = I1+IN 399 if (J1.lt.0) then 400C Out of block to the left - negative index 401 IOut = -J1 402 JP(1,IN) = -(1+I2+BS(2)*I3) 403 JPNN(1,IN) = J1 404 elseif (J1.gt.(BS(1)-1)) then 405C Out of block to the right - negative index 406 IOut = J1 - BS(1) + 1 407 JP(1,IN) = -(1+I2+BS(2)*I3) 408 JPNN(1,IN) = IOut 409 else 410C In block - positive index 411 JP(1,IN) = 1+J1+BS(1)*I2+BS(1)*BS(2)*I3 412 endif 413 enddo 414 endif 415#endif 416 417 418C Y-direction 419#ifdef MPI 420 if (NMesh(2).eq.NMeshG(2)) then 421#endif 422 do IN = -NN,NN 423 J2 = MOD( I2+IN+100*BS(2), BS(2) ) 424 JP(2,IN) = 1+I1+BS(1)*J2+BS(1)*BS(2)*I3 425 enddo 426#ifdef MPI 427 else 428 do IN = -NN,NN 429 J2 = I2+IN 430 if (J2.lt.0) then 431C Out of block to the left - negative index 432 IOut = -J2 433 JP(2,IN) = -(1+I1+BS(1)*I3) 434 JPNN(2,IN) = J2 435 elseif (J2.gt.(BS(2)-1)) then 436C Out of block to the right - negative index 437 IOut = J2 - BS(2) + 1 438 JP(2,IN) = -(1+I1+BS(1)*I3) 439 JPNN(2,IN) = IOut 440 else 441C In block - positive index 442 JP(2,IN) = 1+I1+BS(1)*J2+BS(1)*BS(2)*I3 443 endif 444 enddo 445 endif 446#endif 447 448C Z-direction 449#ifdef MPI 450 if (NMesh(3).eq.NMeshG(3)) then 451#endif 452 do IN = -NN,NN 453 J3 = MOD( I3+IN+100*BS(3), BS(3) ) 454 JP(3,IN) = 1+I1+BS(1)*I2+BS(1)*BS(2)*J3 455 enddo 456#ifdef MPI 457 else 458 do IN = -NN,NN 459 J3 = I3+IN 460 if (J3.lt.0) then 461C Out of block to the left - negative index 462 IOut = -J3 463 JP(3,IN) = -(1+I1+BS(1)*I2) 464 JPNN(3,IN) = J3 465 elseif (J3.gt.(BS(3)-1)) then 466C Out of block to the right - negative index 467 IOut = J3 - BS(3) + 1 468 JP(3,IN) = -(1+I1+BS(1)*I2) 469 JPNN(3,IN) = IOut 470 else 471C In block - positive index 472 JP(3,IN) = 1+I1+BS(1)*I2+BS(1)*BS(2)*J3 473 endif 474 enddo 475 endif 476#endif 477 endif 478 479C Find Jacobian matrix dx/dmesh for adaptative mesh 480 if (mtype .ne. 0) then 481 482C Find dx/dmesh 483 do IC = 1,3 484 do IX = 1,3 485 DXDM(IX,IC) = 0.0_dp 486 do IN = -NN,NN 487 if (mtype .eq. 1) then 488 DXDM(IX,IC) = DXDM(IX,IC) + 489 & XMesh(IX,JP(IC,IN)) * DGDM(IN) 490 else 491 DXDM(IX,IC) = DXDM(IX,IC) + 492 & ( cell(IX,1) * XMesh(1,JP(IC,IN)) + 493 & cell(IX,2) * XMesh(2,JP(IC,IN)) + 494 & cell(IX,3) * XMesh(3,JP(IC,IN)) ) * DGDM(IN) 495 endif 496 enddo 497 enddo 498 enddo 499 500C Find inverse of matrix dx/dmesh 501 call reclat( DXDM, DMDX, 0 ) 502 503C Find differential of volume = determinant of Jacobian matrix 504 DVol = VOLCEL( DXDM ) 505 506C Find the weights for the derivative d(gradF(i))/d(F(j)), of 507C the gradient at point i with respect to the value at point j 508 if (GGA) then 509 do IN = -NN,NN 510 do IC = 1,3 511 do IX = 1,3 512 DGIDFJ(IX,IC,IN) = DMDX(IX,IC) * DGDM(IN) 513 enddo 514 enddo 515 enddo 516 endif 517 518 endif 519 520C Find density and gradient of density at this point 521 do IS = 1,nspin 522 D(IS) = Dens(IP,IS) 523 enddo 524C Test to ensure that densities are always > 0 added to 525C avoid floating point errors in ggaxc. JDG & JMS 526 do IS = 1,min(nspin,2) 527 D(IS) = max(0.0_dp,D(IS)) 528* D(IS) = max(denmin,D(IS)) 529 enddo 530 if (GGA) then 531#ifdef MPI 532 if (Nodes.eq.1) then 533#endif 534 do IS = 1,nspin 535 do IX = 1,3 536 GD(IX,IS) = 0.0_dp 537 do IN = -NN,NN 538 GD(IX,IS) = GD(IX,IS) + 539 & DGIDFJ(IX,1,IN) * Dens(JP(1,IN),IS) + 540 & DGIDFJ(IX,2,IN) * Dens(JP(2,IN),IS) + 541 & DGIDFJ(IX,3,IN) * Dens(JP(3,IN),IS) 542 enddo 543 enddo 544 enddo 545#ifdef MPI 546 else 547 do IS = 1,nspin 548 do IX = 1,3 549 GD(IX,IS) = 0.0_dp 550 do IN = -NN,NN 551 552 if (JP(1,IN).gt.0) then 553 GD(IX,IS) = GD(IX,IS) + 554 & DGIDFJ(IX,1,IN) * Dens(JP(1,IN),IS) 555 else 556 INN = JPNN(1,IN) 557 if (INN.lt.0) then 558 GD(IX,IS) = GD(IX,IS) + 559 & DGIDFJ(IX,1,IN) * bdensX(-JP(1,IN),-INN,IS) 560 else 561 GD(IX,IS) = GD(IX,IS) + 562 & DGIDFJ(IX,1,IN) * bdensX(-JP(1,IN),NN+INN,IS) 563 endif 564 endif 565 566 enddo 567 do IN = -NN,NN 568 if (JP(2,IN).gt.0) then 569 GD(IX,IS) = GD(IX,IS) + 570 & DGIDFJ(IX,2,IN) * Dens(JP(2,IN),IS) 571 else 572 INN = JPNN(2,IN) 573 if (INN.lt.0) then 574 GD(IX,IS) = GD(IX,IS) + 575 & DGIDFJ(IX,2,IN) * bdensY(-JP(2,IN),-INN,IS) 576 else 577 GD(IX,IS) = GD(IX,IS) + 578 & DGIDFJ(IX,2,IN) * bdensY(-JP(2,IN),NN+INN,IS) 579 endif 580 endif 581 enddo 582 do IN = -NN,NN 583 if (JP(3,IN).gt.0) then 584 GD(IX,IS) = GD(IX,IS) + 585 & DGIDFJ(IX,3,IN) * Dens(JP(3,IN),IS) 586 else 587 INN = JPNN(3,IN) 588 if (INN.lt.0) then 589 GD(IX,IS) = GD(IX,IS) + 590 & DGIDFJ(IX,3,IN) * bdensZ(-JP(3,IN),-INN,IS) 591 else 592 GD(IX,IS) = GD(IX,IS) + 593 & DGIDFJ(IX,3,IN) * bdensZ(-JP(3,IN),NN+INN,IS) 594 endif 595 endif 596 enddo 597 enddo 598 enddo 599 endif 600#endif 601 endif 602 603 604C Loop over all functionals 605 do nf = 1,nXCfunc 606 607C Is this a GGA? 608 if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then 609 GGAfunctl = .true. 610 else 611 GGAfunctl = .false. 612 endif 613 614C Find exchange and correlation energy densities and their 615C derivatives with respect to density and density gradient 616 if (GGAfunctl) then 617 call ggaxc( XCauth(nf), irel, nspin, D, GD, 618 & EPSX, EPSC, DEXDD, DECDD, DEXDGD, DECDGD ) 619 else 620 call ldaxc( XCauth(nf), irel, nspin, D, EPSX, EPSC, DEXDD, 621 & DECDD, DVXDN, DVCDN ) 622 endif 623 624C Scale return values by weight for this functional 625 EPSX = XCweightX(nf)*EPSX 626 EPSC = XCweightC(nf)*EPSC 627 do IS = 1,nspin 628 DEXDD(IS) = XCweightX(nf)*DEXDD(IS) 629 DECDD(IS) = XCweightC(nf)*DECDD(IS) 630 enddo 631 if (GGAfunctl) then 632 do IS = 1,nspin 633 DEXDGD(1:3,IS) = XCweightX(nf)*DEXDGD(1:3,IS) 634 DECDGD(1:3,IS) = XCweightC(nf)*DECDGD(1:3,IS) 635 enddo 636 else 637 DVXDN(1:nspin*nspin) = XCweightX(nf)*DVXDN(1:nspin*nspin) 638 DVCDN(1:nspin*nspin) = XCweightC(nf)*DVCDN(1:nspin*nspin) 639 endif 640 641C Add contributions to exchange-correlation energy and its 642C derivatives with respect to density at all points 643 do IS = 1,min(nspin,2) 644 EX = EX + DVol * D(IS) * EPSX 645 EC = EC + DVol * D(IS) * EPSC 646 DX = DX + DVol * D(IS) * EPSX 647 DC = DC + DVol * D(IS) * EPSC 648 enddo 649 do IS = 1,nspin 650 DX = DX - DVol * D(IS) * DEXDD(IS) 651 DC = DC - DVol * D(IS) * DECDD(IS) 652 if (GGAfunctl) then 653 VXC(IP,IS) = VXC(IP,IS) + DEXDD(IS) + DECDD(IS) 654#ifdef MPI 655 if (Nodes.eq.1) then 656#endif 657 do IN = -NN,NN 658 do IC = 1,3 659 do IX = 1,3 660 DX = DX - DVol * Dens(JP(IC,IN),IS) * 661 & DEXDGD(IX,IS) * DGIDFJ(IX,IC,IN) 662 DC = DC - DVol * Dens(JP(IC,IN),IS) * 663 & DECDGD(IX,IS) * DGIDFJ(IX,IC,IN) 664 VXC(JP(IC,IN),IS) = VXC(JP(IC,IN),IS) + 665 & (DEXDGD(IX,IS)+DECDGD(IX,IS))*DGIDFJ(IX,IC,IN) 666 enddo 667 enddo 668 enddo 669#ifdef MPI 670 else 671 do IN = -NN,NN 672 673C X-direction 674 if (JP(1,IN).gt.0) then 675 do IX = 1,3 676 DX = DX - DVol * Dens(JP(1,IN),IS) * 677 & DEXDGD(IX,IS) * DGIDFJ(IX,1,IN) 678 DC = DC - DVol * Dens(JP(1,IN),IS) * 679 & DECDGD(IX,IS) * DGIDFJ(IX,1,IN) 680 VXC(JP(1,IN),IS) = VXC(JP(1,IN),IS) + 681 & (DEXDGD(IX,IS)+DECDGD(IX,IS))*DGIDFJ(IX,1,IN) 682 enddo 683 else 684 INN = JPNN(1,IN) 685 if (INN.lt.0) then 686 do IX = 1,3 687 DX = DX - DVol * bdensX(-JP(1,IN),-INN,IS) * 688 & DEXDGD(IX,IS) * DGIDFJ(IX,1,IN) 689 DC = DC - DVol * bdensX(-JP(1,IN),-INN,IS) * 690 & DECDGD(IX,IS) * DGIDFJ(IX,1,IN) 691 bvxcX(-JP(1,IN),-INN,IS) = (DEXDGD(IX,IS) + 692 & DECDGD(IX,IS))*DGIDFJ(IX,1,IN) + 693 & bvxcX(-JP(1,IN),-INN,IS) 694 enddo 695 else 696 do IX = 1,3 697 DX = DX - DVol * bdensX(-JP(1,IN),NN+INN,IS) * 698 & DEXDGD(IX,IS) * DGIDFJ(IX,1,IN) 699 DC = DC - DVol * bdensX(-JP(1,IN),NN+INN,IS) * 700 & DECDGD(IX,IS) * DGIDFJ(IX,1,IN) 701 bvxcX(-JP(1,IN),NN+INN,IS) = (DEXDGD(IX,IS) + 702 & DECDGD(IX,IS))*DGIDFJ(IX,1,IN) + 703 & bvxcX(-JP(1,IN),NN+INN,IS) 704 enddo 705 endif 706 endif 707C Y-direction 708 if (JP(2,IN).gt.0) then 709 do IX = 1,3 710 DX = DX - DVol * Dens(JP(2,IN),IS) * 711 & DEXDGD(IX,IS) * DGIDFJ(IX,2,IN) 712 DC = DC - DVol * Dens(JP(2,IN),IS) * 713 & DECDGD(IX,IS) * DGIDFJ(IX,2,IN) 714 VXC(JP(2,IN),IS) = VXC(JP(2,IN),IS) + 715 & (DEXDGD(IX,IS)+DECDGD(IX,IS))*DGIDFJ(IX,2,IN) 716 enddo 717 else 718 INN = JPNN(2,IN) 719 if (INN.lt.0) then 720 do IX = 1,3 721 DX = DX - DVol * bdensY(-JP(2,IN),-INN,IS) * 722 & DEXDGD(IX,IS) * DGIDFJ(IX,2,IN) 723 DC = DC - DVol * bdensY(-JP(2,IN),-INN,IS) * 724 & DECDGD(IX,IS) * DGIDFJ(IX,2,IN) 725 bvxcY(-JP(2,IN),-INN,IS) = (DEXDGD(IX,IS) + 726 & DECDGD(IX,IS))*DGIDFJ(IX,2,IN) + 727 & bvxcY(-JP(2,IN),-INN,IS) 728 enddo 729 else 730 do IX = 1,3 731 DX = DX - DVol * bdensY(-JP(2,IN),NN+INN,IS) * 732 & DEXDGD(IX,IS) * DGIDFJ(IX,2,IN) 733 DC = DC - DVol * bdensY(-JP(2,IN),NN+INN,IS) * 734 & DECDGD(IX,IS) * DGIDFJ(IX,2,IN) 735 bvxcY(-JP(2,IN),NN+INN,IS) = (DEXDGD(IX,IS) + 736 & DECDGD(IX,IS))*DGIDFJ(IX,2,IN) + 737 & bvxcY(-JP(2,IN),NN+INN,IS) 738 enddo 739 endif 740 endif 741 742C Z-direction 743 if (JP(3,IN).gt.0) then 744 do IX = 1,3 745 DX = DX - DVol * Dens(JP(3,IN),IS) * 746 & DEXDGD(IX,IS) * DGIDFJ(IX,3,IN) 747 DC = DC - DVol * Dens(JP(3,IN),IS) * 748 & DECDGD(IX,IS) * DGIDFJ(IX,3,IN) 749 VXC(JP(3,IN),IS) = VXC(JP(3,IN),IS) + 750 & (DEXDGD(IX,IS)+DECDGD(IX,IS))*DGIDFJ(IX,3,IN) 751 enddo 752 else 753 INN = JPNN(3,IN) 754 if (INN.lt.0) then 755 do IX = 1,3 756 DX = DX - DVol * bdensZ(-JP(3,IN),-INN,IS) * 757 & DEXDGD(IX,IS) * DGIDFJ(IX,3,IN) 758 DC = DC - DVol * bdensZ(-JP(3,IN),-INN,IS) * 759 & DECDGD(IX,IS) * DGIDFJ(IX,3,IN) 760 bvxcZ(-JP(3,IN),-INN,IS) = (DEXDGD(IX,IS) + 761 & DECDGD(IX,IS)) * DGIDFJ(IX,3,IN) + 762 & bvxcZ(-JP(3,IN),-INN,IS) 763 enddo 764 else 765 do IX = 1,3 766 DX = DX - DVol * bdensZ(-JP(3,IN),NN+INN,IS) * 767 & DEXDGD(IX,IS) * DGIDFJ(IX,3,IN) 768 DC = DC - DVol * bdensZ(-JP(3,IN),NN+INN,IS) * 769 & DECDGD(IX,IS) * DGIDFJ(IX,3,IN) 770 bvxcZ(-JP(3,IN),NN+INN,IS) = (DEXDGD(IX,IS) + 771 & DECDGD(IX,IS))*DGIDFJ(IX,3,IN) + 772 & bvxcZ(-JP(3,IN),NN+INN,IS) 773 enddo 774 endif 775 endif 776 777 enddo 778 endif 779 780#endif 781 else 782 VXC(IP,IS) = VXC(IP,IS) + DEXDD(IS) + DECDD(IS) 783 if (ider .eq. 1) then 784 do JS = 1, nspin 785 KS = JS + (IS-1)*nspin 786 DVXCDN(IP,JS,IS) = DVXCDN(IP,JS,IS) + 787 & DVXDN(KS) + DVCDN(KS) 788 enddo 789 endif 790 endif 791 enddo 792 793C Add contribution to stress due to change in gradient of density 794C originated by the deformation of the mesh with strain 795 if (GGAfunctl) then 796 do JX = 1,3 797 do IX = 1,3 798 do IS = 1,nspin 799 stress(IX,JX) = stress(IX,JX) - DVol * GD(IX,IS) * 800 & ( DEXDGD(JX,IS) + DECDGD(JX,IS) ) 801 enddo 802 enddo 803 enddo 804 endif 805 806C End of loop over functionals 807 enddo 808 809 210 continue 810 811 enddo 812 enddo 813 enddo 814 815#ifdef MPI 816C Return buffer region contributions to VXC to their correct nodes 817 if (GGA.and.Nodes.gt.1) then 818 if (BS(1).ne.NMeshG(1)) then 819 call gathExtMeshData( iDistr, 1, BS(2)*BS(3), NSM, NN, NSPIN, 820 & maxp, NMeshG, bvxcX, VXC ) 821 endif 822 if (BS(2).ne.NMeshG(2)) then 823 call gathExtMeshData( iDistr, 2, BS(1)*BS(3), NSM, NN, NSPIN, 824 & maxp, NMeshG, bvxcY, VXC ) 825 endif 826 if (BS(3).ne.NMeshG(3)) then 827 call gathExtMeshData( iDistr, 3, BS(1)*BS(2), NSM, NN, NSPIN, 828 & maxp, NMeshG, bvxcZ, VXC ) 829 endif 830 endif 831#endif 832 833C Add contribution to stress from the change of volume with strain and 834C divide by volume to get correct stress definition (dE/dStrain)/Vol 835 volume = VOLCEL( cell ) 836 do JX = 1,3 837 stress(JX,JX) = stress(JX,JX) + EX + EC 838 do IX = 1,3 839 stress(IX,JX) = stress(IX,JX) / volume 840 enddo 841 enddo 842C Divide by energy unit 843 EX = EX / EUnit 844 EC = EC / EUnit 845 DX = DX / EUnit 846 DC = DC / EUnit 847 do IS = 1,nspin 848 do I3 = 0,NMesh(3)-1 849 do I2 = 0,NMesh(2)-1 850 do I1 = 0,NMesh(1)-1 851 IP = 1 + I1 + BS(1) * I2 + BS(1) * BS(2) * I3 852 VXC(IP,IS) = VXC(IP,IS) / EUnit 853 enddo 854 enddo 855 enddo 856 enddo 857 do JX = 1,3 858 do IX = 1,3 859 stressl(IX,JX) = stressl(IX,JX) + (stress(IX,JX) / EUnit) 860 enddo 861 enddo 862 863 if (ider.eq.1 .and. .not.GGA) then 864 do IS = 1,nspin 865 do JS = 1,nspin 866 do I3 = 0,NMesh(3)-1 867 do I2 = 0,NMesh(2)-1 868 do I1 = 0,NMesh(1)-1 869 IP = 1 + I1 + BS(1) * I2 + BS(1) * BS(2) * I3 870 DVXCDN(IP,JS,IS) = DVXCDN(IP,JS,IS) / EUnit 871 enddo 872 enddo 873 enddo 874 enddo 875 enddo 876 endif 877 878#ifdef MPI 879C Free memory 880 if (GGA.and.Nodes.gt.1) then 881 if (BS(1).ne.NMeshG(1)) then 882 call de_alloc( bdensX, 'bdensX', 'cellxc' ) 883 call de_alloc( bvxcX, 'bvxcX', 'cellxc' ) 884 endif 885 if (BS(2).ne.NMeshG(2)) then 886 call de_alloc( bdensY, 'bdensY', 'cellxc' ) 887 call de_alloc( bvxcY, 'bvxcY', 'cellxc' ) 888 endif 889 if (BS(3).ne.NMeshG(3)) then 890 call de_alloc( bvxcZ, 'bvxcZ', 'cellxc' ) 891 call de_alloc( bdensZ, 'bdensZ', 'cellxc' ) 892 endif 893 endif 894#endif 895 896C Stop time counter 897#ifdef _TRACE_ 898 call MPI_Barrier( MPI_Comm_World, MPIerror ) 899 call MPItrace_event( 1000, 0 ) 900#endif 901 call timer( 'cellXC', 2 ) 902 903 end subroutine bsc_cellxc 904