1! 2! Copyright (C) 2007-2011 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! This module contains the variables and routines necessary to the implementation 9! of the two-dimensional Coulomb cutoff. Details of the implementation can be found in: 10! 11! Sohier, T., Calandra, M., & Mauri, F. (2017), 12! "Density functional perturbation theory for gated two-dimensional heterostructures: 13! Theoretical developments and application to flexural phonons in graphene." 14! Physical Review B, 96(7), 75448. https://doi.org/10.1103/PhysRevB.96.075448 15! 16!---------------------------------------------------------------------------- 17MODULE Coul_cut_2D 18 !---------------------------------------------------------------------------- 19 !! This module contains the variables and subroutines needed for the 20 !! 2D Coulomb cutoff. 21 ! 22 USE kinds, ONLY : DP 23 USE constants, ONLY : tpi, pi 24 ! 25 SAVE 26 ! 27 LOGICAL :: do_cutoff_2D = .FALSE. 28 !! flag for 2D cutoff. If true, the cutoff is active. 29 REAL(DP) :: lz 30 !! The distance in the out-plne direction after which potential are cut off. 31 ! 32 REAL(DP), ALLOCATABLE :: cutoff_2D(:) 33 !! The factor appended to the Coulomb Kernel to cut off potentials 34 REAL(DP), ALLOCATABLE :: lr_Vloc(:,:) 35 !! The long-range part of the local part of the ionic potential 36 ! 37CONTAINS 38! 39!---------------------------------------------------------------------- 40SUBROUTINE cutoff_fact() 41 !---------------------------------------------------------------------- 42 !! This routine calculates the cutoff factor in G-space and stores it in 43 !! a vector called \(\text{cutoff_2D}(:)\), to be re-used in various routines. 44 !! See Eq.(24) of PRB 96, 075448 45 ! 46 USE kinds 47 USE io_global, ONLY : stdout 48 USE gvect, ONLY : g, ngm, ngmx 49 USE cell_base, ONLY : alat, celldm, at 50 ! 51 IMPLICIT NONE 52 ! 53 ! ... local variables 54 ! 55 INTEGER :: ng, i 56 ! counter over G vectors, cartesian coord. 57 REAL(DP) :: Gzlz, Gplz 58 ! 59 ALLOCATE( cutoff_2D(ngmx) ) 60 ! 61 ! Message to indicate that the cutoff is active. 62 WRITE(stdout, *) "----2D----2D----2D----2D----2D----2D----2D----2D----2D----2D----2D----2D" 63 WRITE(stdout, *) " The code is running with the 2D cutoff" 64 WRITE(stdout, *) " Please refer to:" 65 WRITE(stdout, *) " Sohier, T., Calandra, M., & Mauri, F. (2017), " 66 WRITE(stdout, *) " Density functional perturbation theory for gated two-dimensional heterostructures:" 67 WRITE(stdout, *) " Theoretical developments and application to flexural phonons in graphene." 68 WRITE(stdout, *) " Physical Review B, 96(7), 75448. https://doi.org/10.1103/PhysRevB.96.075448" 69 WRITE(stdout, *) "----2D----2D----2D----2D----2D----2D----2D----2D----2D----2D----2D----2D" 70 ! at(:,i) are the lattice vectors of the simulation cell, a_i, 71 ! in alat units: a_i(:) = at(:,i)/alat 72 ! Check that material is in the x-y plane 73 DO i = 1, 2 74 IF (ABS(at(3,i))>1d-8) WRITE(stdout, *) "2D CODE WILL NOT WORK, 2D MATERIAL NOT IN X-Y PLANE!!" 75 ENDDO 76 ! define cutoff distnce and compute cutoff factor 77 lz = 0.5d0*at(3,3)*alat 78 DO ng = 1, ngm 79 Gplz = SQRT( g(1,ng)**2 + g(2,ng)**2 )*tpi*lz/alat 80 Gzlz = g(3,ng)*tpi*lz/alat 81 cutoff_2D(ng) = 1.0d0 - EXP(-Gplz)*COS(Gzlz) 82 ENDDO 83 ! 84 RETURN 85 ! 86END SUBROUTINE cutoff_fact 87! 88!---------------------------------------------------------------------- 89SUBROUTINE cutoff_lr_Vloc( ) 90 !---------------------------------------------------------------------- 91 !! This routine calculates the long-range part of \(\text{vloc}(g)\) for 92 !! 2D calculations. 93 !! See Eq. (32) of PRB 96, 075448. 94 ! 95 USE kinds 96 USE constants, ONLY : fpi, e2, eps8 97 USE fft_base, ONLY : dfftp 98 USE gvect, ONLY : ngm, gg, g, ngmx 99 ! gg is G^2 in increasing order (in units of tpiba2=(2pi/a)^2) 100 USE ions_base, ONLY : zv, nsp 101 USE uspp_param, ONLY : upf 102 USE cell_base, ONLY : omega, tpiba2 103 ! 104 IMPLICIT NONE 105 ! 106 ! ... local variables 107 ! 108 INTEGER :: ng, nt, ng0 109 REAL(DP) ::fac 110 ! 111 IF (.NOT. ALLOCATED(lr_Vloc)) ALLOCATE( lr_Vloc(ngmx,nsp) ) 112 ! 113 lr_Vloc(:,:) = 0.0d0 114 ! set G=0 value to zero 115 IF (gg(1)<eps8) THEN 116 lr_Vloc(1,:) = 0.0d0 117 ng0 = 2 118 ELSE 119 ng0 = 1 120 ENDIF 121 ! Set g.neq.0 values 122 DO nt = 1, nsp 123 fac = upf(nt)%zp * e2 / tpiba2 124 DO ng = ng0, ngm 125 lr_Vloc(ng,nt) = - fpi / omega* fac * cutoff_2D(ng)* & 126 & EXP( -gg(ng) * tpiba2 * 0.25d0) / gg(ng) 127 ENDDO 128 ENDDO 129 ! 130 RETURN 131 ! 132END SUBROUTINE cutoff_lr_Vloc 133! 134!---------------------------------------------------------------------- 135SUBROUTINE cutoff_local( aux ) 136 !---------------------------------------------------------------------- 137 !! This subroutine is called to re-add the long-range part of the local 138 !! part of the ionic potential, using \(\text{lr_Vloc}\) computed in 139 !! routine \(\texttt{cutoff_lr_Vloc}\). 140 !! See Eq. (33) of PRB 96, 075448 141 ! 142 USE kinds 143 USE fft_base, ONLY : dfftp 144 USE gvect, ONLY : ngm 145 USE vlocal, ONLY : strf 146 USE ions_base, ONLY : nsp 147 ! 148 IMPLICIT NONE 149 ! 150 COMPLEX(DP), INTENT(INOUT):: aux(dfftp%nnr) 151 !! input: local part of ionic potential 152 ! 153 ! ... local variables 154 ! 155 INTEGER :: ng, nt 156 ! 157 DO nt = 1, nsp 158 DO ng = 1, ngm 159 aux(dfftp%nl(ng)) = aux(dfftp%nl(ng)) + lr_Vloc(ng,nt) * strf(ng,nt) 160 ENDDO 161 ENDDO 162 ! 163 RETURN 164 ! 165END SUBROUTINE cutoff_local 166! 167! 168!---------------------------------------------------------------------- 169SUBROUTINE cutoff_hartree( rhog, aux1, ehart ) 170 !---------------------------------------------------------------------- 171 !! This subroutine cuts off the Hartree potential and defines Hartree 172 !! energy accordingly in G-space. 173 !! See Eq. (34) and (41) of PRB 96, 075448 174 ! 175 USE kinds 176 USE gvect, ONLY : ngm, gg , gstart 177 USE io_global, ONLY : stdout 178 ! 179 IMPLICIT NONE 180 ! 181 COMPLEX(DP), INTENT(IN) :: rhog(ngm) 182 !! local potential 183 REAL(DP), INTENT(INOUT) :: aux1(2,ngm) 184 !! Hartree potential 185 REAL(DP), INTENT(INOUT) :: ehart 186 !! Hartree energy 187 ! 188 ! ... local variables 189 ! 190 INTEGER :: ig 191 REAL(DP) :: fac 192 REAL(DP) :: rgtot_re, rgtot_im 193 ! 194 DO ig = gstart, ngm 195 ! 196 fac = 1.D0 / gg(ig) * cutoff_2D(ig) 197 ! 198 rgtot_re = REAL( rhog(ig) ) 199 rgtot_im = AIMAG( rhog(ig) ) 200 ! 201 ehart = ehart + ( rgtot_re**2 + rgtot_im**2 ) * fac 202 ! 203 aux1(1,ig) = rgtot_re * fac 204 aux1(2,ig) = rgtot_im * fac 205 ! 206 ENDDO 207 ! 208 RETURN 209 ! 210END SUBROUTINE cutoff_hartree 211! 212! 213!---------------------------------------------------------------------- 214SUBROUTINE cutoff_ewald( alpha, ewaldg, omega ) 215 !---------------------------------------------------------------------- 216 !! This subroutine defines computes the cutoff version of the 217 !! Ewald sum in G space. 218 !! See Eq. (46) of PRB 96, 075448 219 ! 220 USE kinds 221 USE gvect, ONLY : ngm, gg, gstart 222 USE ions_base, ONLY : zv, nsp, nat, ityp 223 USE cell_base, ONLY : tpiba2, alat 224 USE vlocal, ONLY : strf 225 USE io_global, ONLY : stdout 226 ! 227 IMPLICIT NONE 228 ! 229 REAL(DP), INTENT(IN) :: alpha 230 !! tuning parameter for Ewald LR/SR separation 231 REAL(DP), INTENT(INOUT) :: ewaldg 232 !! Ewald sum 233 REAL(DP), INTENT(IN) :: omega 234 !! unit-cell volume 235 ! 236 ! ... local variables 237 ! 238 INTEGER :: ng, nt, na, nr, ir, iz, nz, rmax 239 COMPLEX(DP) :: rhon 240 REAL(DP) :: rp, z 241 ! 242 ! The G=0 component of the long-ranged local part of the 243 ! pseudopotential minus the Hartree potential is set to 0. 244 ! This is equivalent to substracting the finite non-singular 245 ! part of the ionic potential at G=0. See Appendix D.2 of PRB 96, 075448. 246 ! In practice, with respect to the 3D ewald sum, we must subtract the 247 ! G=0 energy term - 4 pi/omega * e**2/2 * charge**2 / alpha / 4.0d0. 248 ! That is, we simply set setting ewaldg(G=0)=0. 249 ! 250 ewaldg = 0.0d0 251 ! now the G.neq.0 terms 252 DO ng = gstart, ngm 253 rhon = (0.d0, 0.d0) 254 DO nt = 1, nsp 255 rhon = rhon + zv(nt) * CONJG(strf(ng,nt) ) 256 ENDDO 257 ewaldg = ewaldg + ABS(rhon)**2 * EXP( - gg(ng) * tpiba2 /& 258 alpha / 4.d0) / gg(ng)*cutoff_2D(ng) / tpiba2 259 ENDDO 260 ewaldg = 2.d0 * tpi / omega * ewaldg 261 ! 262 ! ... here add the other constant term (Phi_self) 263 ! 264 IF (gstart == 2) THEN 265 DO na = 1, nat 266 ewaldg = ewaldg - zv(ityp(na))**2 * SQRT(8.d0 / tpi * & 267 alpha) 268 ENDDO 269 ENDIF 270 ! 271 RETURN 272 ! 273END SUBROUTINE cutoff_ewald 274! 275! 276!---------------------------------------------------------------------- 277SUBROUTINE cutoff_force_ew( aux, alpha ) 278 !---------------------------------------------------------------------- 279 !! This subroutine cuts off the Ewald contribution to the forces. More 280 !! precisely, it cuts off the LR ion-ion potential that is then used to 281 !! compute the Ewald forces. 282 !! See Eq. (55) of PRB 96, 075448 (note that Eq. (56), derived from Eq. (55), 283 !! looks somewhat different from what is implemented in the code, but it is 284 !! equivalent). 285 ! 286 USE kinds 287 USE gvect, ONLY : ngm, gg , gstart 288 USE cell_base, ONLY : tpiba2, alat 289 ! 290 IMPLICIT NONE 291 ! 292 COMPLEX(DP), INTENT(INOUT) :: aux(ngm) 293 !! long-range part of the ionic potential 294 REAL(DP), INTENT(IN) :: alpha 295 !! tuning parameter for the LR/SR separation 296 ! 297 ! ... local variables 298 ! 299 INTEGER :: ig 300 ! 301 DO ig = gstart, ngm 302 aux(ig) = aux(ig) * EXP( - gg(ig) * tpiba2 / alpha / 4.d0) & 303 / (gg(ig) * tpiba2) * cutoff_2D(ig) 304 ENDDO 305 ! 306 RETURN 307 ! 308END SUBROUTINE cutoff_force_ew 309! 310! 311!---------------------------------------------------------------------- 312SUBROUTINE cutoff_force_lc( aux, forcelc ) 313 !---------------------------------------------------------------------- 314 !! This subroutine re-adds the cutoff contribution from the long-range 315 !! local part of the ionic potential to the forces. In the 2D code, this 316 !! contribution is missing from the \(\text{Vloc}\). 317 !! See Eq. (54) of PRB 96, 075448. 318 ! 319 USE kinds 320 USE gvect, ONLY : ngm, gg, g , gstart 321 USE constants, ONLY : fpi, e2, eps8, tpi 322 USE uspp_param, ONLY : upf 323 USE cell_base, ONLY : tpiba2, alat, omega 324 USE ions_base, ONLY : nat, zv, tau, ityp 325 USE io_global, ONLY : stdout 326 USE fft_base, ONLY : dfftp 327 ! 328 IMPLICIT NONE 329 ! 330 COMPLEX(DP), INTENT(IN) :: aux(dfftp%nnr) 331 !! local ionic potential 332 REAL(DP), INTENT(INOUT) :: forcelc(3,nat) 333 !! corresponding force contribution 334 ! 335 ! ... local variables 336 ! 337 REAL(DP) :: arg 338 INTEGER :: ig, na, ipol 339 ! 340 DO na = 1, nat 341 DO ig = gstart, ngm 342 arg = (g(1,ig) * tau(1,na) + g(2,ig) * tau(2,na) + & 343 g(3,ig) * tau(3,na) ) * tpi 344 DO ipol = 1, 3 345 forcelc(ipol,na) = forcelc (ipol,na) + tpi / alat * & 346 g(ipol,ig) * lr_Vloc(ig, ityp(na)) * omega * & 347 ( SIN(arg)*DBLE(aux(dfftp%nl(ig))) + COS(arg)*AIMAG(aux(dfftp%nl(ig))) ) 348 ENDDO 349 ENDDO 350 ENDDO 351 ! 352 RETURN 353 ! 354END SUBROUTINE cutoff_force_lc 355! 356! 357!---------------------------------------------------------------------- 358SUBROUTINE cutoff_stres_evloc( psic_G, strf, evloc ) 359 !---------------------------------------------------------------------- 360 !! This subroutine adds the contribution from the cutoff long-range part 361 !! of the local part of the ionic potential to \(\text{evloc}\). 362 !! evloc corresponds to the delta term in Eq. (63) of PRB 96, 075448. 363 !! It is the energy of the electrons in the local ionic potential. 364 !! Note that it is not calculated as such (by itself) in the standard code. 365 !! Indeed, it is "hidden" in the sum of KS eigenvalues. That is why we need 366 !! to re-compute it here for the stress. 367 ! 368 USE kinds 369 USE ions_base, ONLY : ntyp => nsp 370 USE gvect, ONLY : ngm , gstart 371 USE io_global, ONLY : stdout 372 USE fft_base, ONLY : dfftp 373 ! 374 IMPLICIT NONE 375 ! 376 COMPLEX(DP), INTENT(IN) :: psic_G(dfftp%nnr) 377 !! charge density in G space 378 COMPLEX(DP), INTENT(IN) :: strf(ngm,ntyp) 379 !! the structure factor 380 REAL(DP), INTENT(INOUT) :: evloc 381 !! the energy of the electrons in the local ionic potential 382 ! 383 ! ... local variables 384 ! 385 INTEGER :: ng, nt 386 ! 387 ! If gstart=2, it means g(1) is G=0, but we have nothing to add for G=0 388 ! So we start at gstart. 389 DO nt = 1, ntyp 390 DO ng = gstart, ngm 391 evloc = evloc + DBLE( CONJG(psic_G(dfftp%nl(ng))) * strf(ng,nt) ) & 392 * lr_Vloc(ng,nt) 393 ENDDO 394 ENDDO 395 ! 396 RETURN 397 ! 398END SUBROUTINE cutoff_stres_evloc 399! 400! 401!---------------------------------------------------------------------- 402SUBROUTINE cutoff_stres_sigmaloc( psic_G, strf, sigmaloc ) 403 !---------------------------------------------------------------------- 404 !! This subroutine adds the contribution from the cutoff long-range part 405 !! of the local part of the ionic potential to the rest of the 406 !! \(\text{sigmaloc}\). That is, the rest of Eq. (63) of PRB 96, 075448. 407 ! 408 USE kinds 409 USE ions_base, ONLY : ntyp => nsp 410 USE constants, ONLY : eps8 411 USE gvect, ONLY : ngm, g, gg, gstart 412 USE cell_base, ONLY : tpiba, tpiba2, alat, omega 413 USE io_global, ONLY : stdout 414 USE fft_base, ONLY : dfftp 415 ! 416 IMPLICIT NONE 417 ! 418 COMPLEX(DP), INTENT(IN) :: psic_G(dfftp%nnr) 419 !! charge density in G space 420 COMPLEX(DP), INTENT(IN) :: strf(ngm,ntyp) 421 !! the structure factor 422 REAL(DP), INTENT(INOUT) :: sigmaloc(3,3) 423 !! stress contribution for the local ionic potential 424 ! 425 ! ... local variables 426 ! 427 INTEGER :: ng, nt, l, m 428 REAL(DP) :: Gp, G2lzo2Gp, beta, dlr_Vloc 429 ! 430 ! no G=0 contribution 431 DO nt = 1, ntyp 432 DO ng = gstart, ngm 433 ! 434 Gp = SQRT( g(1,ng)**2 + g(2,ng)**2 )*tpiba 435 ! below is a somewhat cumbersome way to define beta of Eq. (61) of PRB 96, 075448 436 IF (Gp < eps8) THEN 437 ! G^2*lz/2|Gp| 438 G2lzo2Gp = 0.0d0 439 beta = 0.0d0 440 ELSE 441 G2lzo2Gp = gg(ng)*tpiba2*lz/2.0d0/Gp 442 beta = G2lzo2Gp*(1.0d0-cutoff_2D(ng))/cutoff_2D(ng) 443 ENDIF 444 ! dlr_vloc corresponds to the derivative of the long-range local ionic potential 445 ! with respect to G 446 DO l = 1, 3 447 IF (l == 3) THEN 448 dlr_Vloc = - 1.0d0/ (gg(ng)*tpiba2) * lr_Vloc(ng,nt) & 449 * (1.0d0+ gg(ng)*tpiba2/4.0d0) 450 ELSE 451 dlr_Vloc = - 1.0d0/ (gg(ng)*tpiba2) * lr_Vloc(ng,nt) & 452 * (1.0d0- beta + gg(ng)*tpiba2/4.0d0) 453 ENDIF 454 ! 455 DO m = 1, l 456 sigmaloc(l,m) = sigmaloc(l,m) + DBLE( CONJG( psic_G(dfftp%nl(ng) ) ) & 457 * strf(ng,nt) ) * 2.0d0 * dlr_Vloc & 458 * tpiba2 * g(l,ng) * g(m,ng) 459 ENDDO 460 ENDDO 461 ! 462 ENDDO 463 ENDDO 464 ! 465 RETURN 466 ! 467END SUBROUTINE cutoff_stres_sigmaloc 468! 469! 470!---------------------------------------------------------------------- 471SUBROUTINE cutoff_stres_sigmahar( psic_G, sigmahar ) 472 !---------------------------------------------------------------------- 473 !! This subroutine cuts off the Hartree part of the stress. 474 !! See Eq. (62) of PRB 96, 075448. 475 ! 476 USE kinds 477 USE gvect, ONLY : ngm, g, gg, gstart 478 USE constants, ONLY : eps8 479 USE cell_base, ONLY : tpiba2, alat, tpiba 480 USE io_global, ONLY : stdout 481 USE fft_base, ONLY : dfftp 482 ! 483 IMPLICIT NONE 484 ! 485 COMPLEX(DP), INTENT(IN) :: psic_G(dfftp%nnr) 486 !! charge density in G-space 487 REAL(DP), INTENT(INOUT) :: sigmahar(3,3) 488 !! hartree contribution to stress 489 ! 490 ! ... local variables 491 ! 492 INTEGER :: ng, nt, l, m 493 REAL(DP) :: Gp, G2lzo2Gp, beta, shart, g2, fact 494 ! 495 DO ng = gstart, ngm 496 Gp = SQRT( g(1,ng)**2 + g(2,ng)**2 )*tpiba 497 IF (Gp < eps8) THEN 498 G2lzo2Gp = 0.0d0 499 beta = 0.0d0 500 ELSE 501 G2lzo2Gp = gg(ng)*tpiba2*lz/2.0d0/Gp 502 beta = G2lzo2Gp*(1.0d0-cutoff_2D(ng))/cutoff_2D(ng) 503 ENDIF 504 g2 = gg (ng) * tpiba2 505 shart = psic_G(dfftp%nl(ng)) * CONJG(psic_G(dfftp%nl(ng))) / g2 * cutoff_2D(ng) 506 DO l = 1, 3 507 IF (l == 3) THEN 508 fact = 1.0d0 509 ELSE 510 fact = 1.0d0 - beta 511 ENDIF 512 DO m = 1, l 513 sigmahar(l,m) = sigmahar(l,m) + shart * tpiba2 * 2 * & 514 g(l,ng) * g(m,ng) / g2 * fact 515 ENDDO 516 ENDDO 517 ENDDO 518 !sigma is multiplied by 0.5*fpi*e2 after 519 ! 520 RETURN 521 ! 522END SUBROUTINE cutoff_stres_sigmahar 523! 524! 525!---------------------------------------------------------------------- 526SUBROUTINE cutoff_stres_sigmaewa( alpha, sdewald, sigmaewa ) 527 !---------------------------------------------------------------------- 528 !! This subroutine cuts off the Ewald part of the stress. 529 !! See Eq. (64) in PRB 96 075448 530 ! 531 USE kinds 532 USE ions_base, ONLY : nat, zv, tau, ityp 533 USE constants, ONLY : e2, eps8 534 USE gvect, ONLY : ngm, g, gg, gstart 535 USE cell_base, ONLY : tpiba2, alat, omega, tpiba 536 USE io_global, ONLY : stdout 537 ! 538 IMPLICIT NONE 539 ! 540 REAL(DP), INTENT(IN) :: alpha 541 !! tuning param for LR/SR separation 542 REAL(DP), INTENT(INOUT) :: sigmaewa(3,3) 543 !! ewald contribution to stress 544 REAL(DP), INTENT(INOUT) :: sdewald 545 !! constant and diagonal terms 546 ! 547 ! ... local variables 548 ! 549 INTEGER :: ng, na, l, m 550 REAL(DP) :: Gp, G2lzo2Gp, beta, sewald, g2, g2a, arg, fact 551 COMPLEX(DP) :: rhostar 552 ! 553 ! g(1) is a problem if it's G=0, because we divide by G^2. 554 ! So start at gstart. 555 ! fact=1.0d0, gamma_only not implemented 556 ! G=0 componenent of the long-range part of the local part of the 557 ! pseudopotminus the Hartree potential is set to 0. 558 ! in other words, sdewald=0. 559 ! sdewald is the last term in equation B1 of PRB 32 3792. 560 ! See also similar comment for ewaldg in cutoff_ewald routine 561 ! 562 sdewald = 0.0D0 563 DO ng = gstart, ngm 564 Gp = SQRT( g(1,ng)**2 + g(2,ng)**2 )*tpiba 565 IF (Gp < eps8) THEN 566 G2lzo2Gp = 0.0d0 567 beta = 0.0d0 568 ELSE 569 G2lzo2Gp = gg(ng)*tpiba2*lz/2.0d0/Gp 570 beta = G2lzo2Gp*(1.0d0-cutoff_2D(ng))/cutoff_2D(ng) 571 ENDIF 572 g2 = gg(ng) * tpiba2 573 g2a = g2 / 4.d0 / alpha 574 rhostar = (0.d0, 0.d0) 575 DO na = 1, nat 576 arg = (g(1,ng) * tau(1,na) + g(2,ng) * tau(2,na) + & 577 g(3,ng) * tau(3,na) ) * tpi 578 rhostar = rhostar + zv (ityp(na) ) * CMPLX(COS(arg), SIN(arg), KIND=DP) 579 ENDDO 580 rhostar = rhostar / omega 581 sewald = tpi * e2 * EXP(-g2a) / g2* cutoff_2D(ng) * ABS(rhostar)**2 582 ! ... sewald is an other diagonal term that is similar to the diagonal terms 583 ! in the other stress contributions. It basically gives a term prop to 584 ! the ewald energy 585 sdewald = sdewald-sewald 586 DO l = 1, 3 587 IF (l == 3) THEN 588 fact = (g2a + 1.0d0) 589 ELSE 590 fact = (1.0d0+g2a-beta) 591 ENDIF 592 ! 593 DO m = 1, l 594 sigmaewa(l,m) = sigmaewa(l,m) + sewald * tpiba2 * 2.d0 * & 595 g(l,ng) * g(m,ng) / g2 * fact 596 ENDDO 597 ENDDO 598 ! 599 ENDDO 600 ! 601 RETURN 602 ! 603END SUBROUTINE cutoff_stres_sigmaewa 604! 605! 606END MODULE Coul_cut_2D 607