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 9! 10! The old dhscf has been split in two parts: an initialization routine 11! (dhscf_init), which is called after every geometry change but before 12! the main scf loop, and a dhscf proper, which is called at every step 13! of the scf loop 14! 15! NOTE that the mesh initialization part is now done *unconditionally* 16! in dhscf_init, i.e, after *every* geometry change, even if the 17! change does not involve a cell change. The reason is to avoid 18! complexity, since now the mesh parallel distributions will depend on 19! the detailed atomic positions even if the cell does not change. 20! 21! Besides, the relative cost of a "mesh only" initialization is negligible. 22! The only real observable effect would be a printout of "initmesh" data 23! at every geometry iteration. 24! 25 module m_dhscf 26 27! 28! To facilitate the communication among dhscf_init and dhscf, 29! some arrays that hold data which do not change during the SCF loop 30! have been made into module variables 31 32! Some others are scratch, such as nmpl, ntpl, etc 33 34 use precision, only : dp, grid_p 35 use m_dfscf, only : dfscf 36 implicit none 37 38 real(grid_p), pointer :: rhopcc(:), rhoatm(:), Vna(:) 39 real(dp) :: Uharrs ! Harris energy 40 logical :: IsDiag 41 integer :: nml(3), ntml(3), npcc, nmpl, ntpl 42 integer :: nbcell 43 real(dp) :: bcell(3,3), cell(3,3), 44 & dvol, field(3), rmax, scell(3,3) 45 real(dp) :: G2mesh = 0.0_dp 46 logical :: debug_dhscf = .false. 47 logical :: use_bsc_cellxc 48 49 character(len=*), parameter :: debug_fmt = 50 & '('' -- dhscf : Node '',i3,tr1,a,4(tr1,e20.7))' 51 52 public :: dhscf_init, dhscf 53 54 CONTAINS 55 56 subroutine dhscf_init(nspin, norb, iaorb, iphorb, 57 & nuo, nuotot, nua, na, 58 & isa, xa, indxua, ucell, 59 & mscell, G2max, ntm, 60 & maxnd, numd, listdptr, listd, datm, 61 & Fal, stressl) 62 63 use precision, only : dp, grid_p 64 use parallel, only : Node, Nodes 65 use atmfuncs, only : rcut, rcore 66 use fdf 67 use sys, only : die 68 use mesh, only : xdsp, nsm, nsp, meshLim 69 use parsing 70 71 use alloc, only : re_alloc, de_alloc 72 use siesta_options, only : harrisfun 73 use meshsubs, only : PartialCoreOnMesh 74 use meshsubs, only : NeutralAtomOnMesh 75 use meshsubs, only : PhiOnMesh 76 use meshsubs, only : InitMesh 77 use meshsubs, only : InitAtomMesh 78 use meshsubs, only : setupExtMesh 79 use meshsubs, only : distriPhiOnMesh 80 81 use moreMeshSubs, only : setMeshDistr, distMeshData 82 use moreMeshSubs, only : UNIFORM, QUADRATIC, LINEAR 83 use moreMeshSubs, only : TO_SEQUENTIAL, TO_CLUSTER, KEEP 84 85 use meshdscf, only : createLocalDscfPointers 86 87 use iogrid_netcdf, only: set_box_limits 88#ifdef NCDF_4 89 use files, only : slabel 90 use m_ncdf_siesta, only : cdf_init_grid 91 use m_ncdf_io, only : cdf_init_mesh 92#endif 93 use m_efield, only : initialize_efield, acting_efield 94 use m_efield, only : user_specified_field 95 use dipole_m, only : init_dipole_correction 96 97 use m_doping_uniform, only: initialize_doping_uniform 98 use m_doping_uniform, only: compute_doping_structs_uniform, 99 $ doping_active 100 use m_rhog, only: rhog, rhog_in 101 use m_rhog, only: order_rhog 102 use siesta_options, only: mix_charge 103#ifdef NCDF_4 104 use siesta_options, only: write_cdf 105#endif 106#ifdef MPI 107 use mpi_siesta 108#endif 109 110 use m_mesh_node, only: init_mesh_node 111 use m_charge_add, only: init_charge_add 112 use m_hartree_add, only: init_hartree_add 113 114 use m_ts_global_vars,only: TSmode 115 use m_ts_options, only : IsVolt, N_Elec, Elecs 116 use m_ts_electype, only : Elec_ortho2semi_inf 117 use m_ts_voltage, only : ts_init_voltage 118 119 implicit none 120 integer, intent(in) :: nspin, norb, iaorb(norb), iphorb(norb), 121 & nuo, nuotot, nua, na, isa(na), 122 & indxua(na), mscell(3,3), maxnd, 123 & numd(nuo), listdptr(nuo), listd(maxnd) 124 real(dp), intent(in) :: xa(3,na), ucell(3,3), datm(norb) 125 real(dp), intent(inout) :: g2max 126 integer, intent(inout) :: ntm(3) 127 real(dp), intent(inout) :: Fal(3,nua), stressl(3,3) 128 129 real(dp), parameter :: tiny = 1.e-12_dp 130 integer :: io, ia, iphi, is, n, i, j 131 integer :: nsc(3), nsd 132 real(dp) :: DStres(3,3), volume 133 real(dp), external :: volcel, ddot 134 real(grid_p) :: dummy_Drho(1,1), dummy_Vaux(1), 135 & dummy_Vscf(1) 136 logical, save :: frstme = .true. ! Keeps state 137 real(grid_p), pointer :: Vscf(:,:), rhoatm_par(:) 138 integer, pointer :: numphi(:), numphi_par(:) 139 character(len=10) :: shape 140 141 integer :: nm(3) ! For call to initMesh 142 logical :: need_gradients_in_xc, is_vdw 143 144 ! Transport direction (unit-cell aligned) 145 integer :: iE 146 logical :: is_ortho 147!--------------------------------------------------------------------- BEGIN 148#ifdef DEBUG 149 call write_debug( ' PRE dhscf_init' ) 150#endif 151C ---------------------------------------------------------------------- 152C General initialisation 153C ---------------------------------------------------------------------- 154C Start time counter 155 call timer( 'DHSCF_Init', 1 ) 156 157 nsd = min(nspin,2) 158 nullify(Vscf,rhoatm_par) 159 160 if (frstme) then 161 debug_dhscf = fdf_get('Debug.DHSCF', .false.) 162 163 nullify( xdsp, rhopcc, Vna, rhoatm ) 164C nsm lives in module m_dhscf now !! AG** 165 nsm = fdf_integer( 'MeshSubDivisions', 2 ) 166 nsm = max( nsm, 1 ) 167 168C Set mesh sub-division variables & perform one off allocation 169 nsp = nsm*nsm*nsm 170 call re_alloc( xdsp, 1, 3, 1, nsp, 'xdsp', 'dhscf_init' ) 171 172 endif ! First time 173 174 175C ---------------------------------------------------------------------- 176C Orbital initialisation : part 1 177C ---------------------------------------------------------------------- 178 179C Find the maximum orbital radius 180 rmax = 0.0_dp 181 do io = 1, norb 182 ia = iaorb(io) ! Atomic index of each orbital 183 iphi = iphorb(io) ! Orbital index of each orbital in its atom 184 is = isa(ia) ! Species index of each atom 185 rmax = max( rmax, rcut(is,iphi) ) 186 enddo 187 188C Start time counter for mesh initialization 189 call timer( 'DHSCF1', 1 ) 190 191C ---------------------------------------------------------------------- 192C Unit cell handling 193C ---------------------------------------------------------------------- 194C Find diagonal unit cell and supercell 195 call digcel( ucell, mscell, cell, scell, nsc, IsDiag ) 196 if (.not.IsDiag) then 197 if (Node.eq.0) then 198 write(6,'(/,a,3(/,a,3f12.6,a,i6))') 199 & 'DHSCF: WARNING: New shape of unit cell and supercell:', 200 & ('DHSCF:',(cell(i,j),i=1,3),' x',nsc(j),j=1,3) 201 endif 202 endif 203 204C Find the system shape 205 call shaper( cell, nua, isa, xa, shape, nbcell, bcell ) 206 207C Find system volume 208 volume = volcel( cell ) 209 210C ---------------------------------------------------------------------- 211C Mesh initialization 212C ---------------------------------------------------------------------- 213 call InitMesh( na, cell, norb, iaorb, iphorb, isa, rmax, 214 & G2max, G2mesh, nsc, nmpl, nm, 215 & nml, ntm, ntml, ntpl, dvol ) 216 217! Setup box descriptors for each processor, 218! held in module iogrid_netcdf 219 call set_box_limits( ntm, nsm ) 220 221 ! Initialize information on local mesh for each node 222 call init_mesh_node( cell, ntm , meshLim , nsm ) 223 224 ! Setup charge additions in the mesh 225 call init_charge_add(cell, ntm) 226 227 ! Setup Hartree additions in the mesh 228 call init_hartree_add(cell, ntm) 229 230#ifdef NCDF_4 231 ! Initialize the box for each node... 232 call cdf_init_mesh( ntm, nsm ) 233 if ( write_cdf ) then 234 call cdf_init_grid( trim(slabel)//'.nc', ntm) 235 end if 236#endif 237 238C Stop time counter for mesh initialization 239 call timer( 'DHSCF1', 2 ) 240C ---------------------------------------------------------------------- 241C End of mesh initialization 242C ---------------------------------------------------------------------- 243 244C ---------------------------------------------------------------------- 245C Initialize atomic orbitals, density and potential 246C ---------------------------------------------------------------------- 247C Start time counter for atomic initializations 248 call timer( 'DHSCF2', 1 ) 249 250 if (nodes.gt.1) then 251 call setMeshDistr( UNIFORM, nsm, nsp, 252 & nml, nmpl, ntml, ntpl ) 253 endif 254 255C Initialise quantities relating to the atom-mesh positioning 256 call InitAtomMesh( UNIFORM, na, xa ) 257 258! Check if we need stencils in cellxc, and detect incompatibility 259! with Harris forces 260 call xc_setup(need_gradients_in_xc, is_vdw) 261 262 if (need_gradients_in_xc .and. harrisfun) then 263 call die("** Harris forces not implemented for non-LDA XC") 264 endif 265 266 ! Give the opportunity to use BSC's version, 267 ! unless we have a vdw functional 268 use_bsc_cellxc = fdf_get("XC.Use.BSC.CellXC", .false.) 269 if (is_vdw .and. use_bsc_cellxc) then 270 call message("INFO", 271 $ "BSC's cellxc cannot be used with vdW functionals") 272 use_bsc_cellxc = .false. 273 endif 274 275C Compute the number of orbitals on the mesh and recompute the 276C partions for every processor in order to have a similar load 277C in each of them. 278 nullify( numphi ) 279 call re_alloc( numphi, 1, nmpl, 'numphi', 'dhscf_init' ) 280!$OMP parallel do default(shared), private(i) 281 do i= 1, nmpl 282 numphi(i) = 0 283 enddo 284!$OMP end parallel do 285 286 call distriPhiOnMesh( nm, nmpl, norb, iaorb, iphorb, 287 & isa, numphi, need_gradients_in_xc ) 288 289C Find if there are partial-core-corrections for any atom 290 npcc = 0 291 do ia = 1,na 292 if (rcore(isa(ia)) .gt. tiny) npcc = 1 293 enddo 294 295C Find partial-core-correction energy density 296C Vscf and Vaux are not used here 297 call re_alloc( rhopcc, 1, ntpl*npcc+1, 'rhopcc', 'dhscf_init' ) 298 if (npcc .eq. 1) then 299 call PartialCoreOnMesh( na, isa, ntpl, rhopcc, indxua, 300 & nsd, dvol, volume, dummy_Vscf, dummy_Vaux, Fal, stressl, 301 & .false., .false. ) 302 call reord( rhopcc, rhopcc, nml, nsm, TO_SEQUENTIAL ) 303 if ( debug_dhscf ) then 304 write(*,debug_fmt) Node,'rhopcc',sum(rhopcc)*dvol 305 end if 306 endif 307 308C Find neutral-atom potential 309C Drho is not used here 310 call re_alloc( Vna, 1, ntpl, 'Vna', 'dhscf_init' ) 311 call NeutralAtomOnMesh( na, isa, ntpl, Vna, indxua, dvol, 312 & volume, dummy_DRho, Fal, stressl, 313 & .false., .false. ) 314 call reord( Vna, Vna, nml, nsm, TO_SEQUENTIAL ) 315 if ( debug_dhscf ) then 316 write(*,debug_fmt) Node,'Vna',sqrt(sum(Vna**2)) 317 end if 318 319 if (nodes.gt.1) then 320 if (node .eq. 0) then 321 write(6,"(a)") "Setting up quadratic distribution..." 322 endif 323 call setMeshDistr( QUADRATIC, nsm, nsp, 324 & nml, nmpl, ntml, ntpl ) 325 326C Create extended mesh arrays for the second data distribution 327 call setupExtMesh( QUADRATIC, rmax ) 328 329C Compute atom positions for the second data distribution 330 call InitAtomMesh( QUADRATIC, na, xa ) 331 endif 332 333C Calculate orbital values on mesh 334! numphi has already been computed in distriPhiOnMesh 335! in the UNIFORM distribution 336 if (nodes.eq.1) then 337 numphi_par => numphi 338 else 339 nullify(numphi_par) 340 call re_alloc( numphi_par, 1, nmpl, 'numphi_par', 341 & 'dhscf_init' ) 342 call distMeshData( UNIFORM, numphi, QUADRATIC, 343 & numphi_par, KEEP ) 344 endif 345 346 call PhiOnMesh( nmpl, norb, iaorb, iphorb, isa, numphi_par ) 347 348 if (nodes.gt.1) then 349 call de_alloc( numphi_par, 'numphi_par', 'dhscf_init' ) 350 endif 351 call de_alloc( numphi, 'numphi', 'dhscf_init' ) 352 353C ---------------------------------------------------------------------- 354C Create sparse indexing for Dscf as needed for local mesh 355C Note that this is done in the QUADRATIC distribution 356C since 'endpht' (computed finally in PhiOnMesh and stored in 357C meshphi module) is in that distribution. 358C ---------------------------------------------------------------------- 359 if (Nodes.gt.1) then 360 call CreateLocalDscfPointers( nmpl, nuotot, numd, listdptr, 361 & listd ) 362 endif 363 364C ---------------------------------------------------------------------- 365C Calculate terms relating to the neutral atoms on the mesh 366C ---------------------------------------------------------------------- 367C Find Harris (sum of atomic) electron density 368 call re_alloc( rhoatm_par, 1, ntpl, 'rhoatm_par', 'dhscf_init' ) 369 call rhooda( norb, nmpl, datm, rhoatm_par, iaorb, iphorb, isa ) 370! rhoatm_par comes out of here in clustered form in QUADRATIC dist 371 372C Routine Poison should use the uniform data distribution 373 if (nodes.gt.1) then 374 call setMeshDistr( UNIFORM, nsm, nsp, 375 & nml, nmpl, ntml, ntpl ) 376 endif 377 378C Create Rhoatm using UNIFORM distr, in sequential form 379 call re_alloc( rhoatm, 1, ntpl, 'rhoatm', 'dhscf_init' ) 380 call distMeshData( QUADRATIC, rhoatm_par, 381 & UNIFORM, rhoatm, TO_SEQUENTIAL ) 382 383 if ( debug_dhscf ) then 384 write(*,debug_fmt) Node,'rhoatm', sum(rhoatm)*dvol 385 end if 386 387! 388! AG: The initialization of doping structs could be done here now, 389! in the uniform distribution, and with a simple loop over 390! rhoatm. 391 392 if (frstme) call initialize_doping_uniform() 393 if (doping_active) then 394 call compute_doping_structs_uniform(ntpl,rhoatm,nsd) 395 ! Will get the global number of hit points 396 ! Then, the doping density to be added can be simply computed 397 endif 398 399C Allocate Temporal array 400 call re_alloc( Vscf, 1, ntpl, 1, 1, 'Vscf', 'dhscf_init' ) 401 402C Vscf is filled here but not used later 403C Uharrs is computed (and saved) 404C DStres is computed but not used later 405 call poison( cell, ntml(1), ntml(2), ntml(3), ntm, rhoatm, 406 & Uharrs, Vscf, DStres, nsm ) 407 call de_alloc( Vscf, 'Vscf', 'dhscf_init' ) 408 409! Always deallocate rhoatm_par, as it was used even if nodes=1 410 call de_alloc( rhoatm_par, 'rhoatm_par', 'dhscf_init' ) 411 412 if (mix_charge) then 413 call re_alloc( rhog, 1, 2, 1, ntpl, 1, nspin, 414 $ 'Rhog', 'dhscf_init' ) 415 call re_alloc( rhog_in, 1, 2, 1, ntpl, 1, nspin, 416 $ 'Rhog_in', 'dhscf_init' ) 417 call order_rhog(cell, ntml(1), ntml(2), ntml(3), ntm, nsm) 418 endif 419 420C Stop time counter for atomic initializations 421 call timer( 'DHSCF2', 2 ) 422 423C ---------------------------------------------------------------------- 424C At the end of initializations: 425! We are in the UNIFORM distribution 426! Rhoatm, Rhopcc and Vna are in UNIFORM dist, sequential form 427! The index array endpht is in the QUADRATIC distribution 428C ---------------------------------------------------------------------- 429 if (frstme) then 430 call initialize_efield() 431 call init_dipole_correction(nbcell, bcell, acting_efield) 432 end if 433 434 ! Check if we need to add the potential 435 ! corresponding to the voltage-drop. 436 if ( TSmode ) then 437 ! These routines are important if there are cell-changes 438 if ( IsVolt ) then 439 call ts_init_voltage( cell, nua, xa, ntm) 440 end if 441 442 if ( acting_efield ) then 443 ! We do not allow the electric field for 444 ! transiesta runs with V = 0, either. 445 ! It does not make sense, only for fields perpendicular 446 ! to the applied bias. 447 448 ! We need to check that the e-field is perpendicular 449 ! to the transport direction, and that the system is 450 ! either a chain, or a slab. 451 ! However, due to the allowance of a dipole correction 452 ! along the transport direction for buffer calculations 453 ! we have to allow all shapes. (atom is not transiesta 454 ! compatible anyway) 455 456 ! check that we do not violate the periodicity 457 if ( Node .eq. 0 ) then 458 write(*,'(/,2(2a,/))') 'ts-WARNING: ', 459 & 'E-field/dipole-correction! ', 460 & 'ts-WARNING: ', 461 & 'I hope you know what you are doing!' 462 end if 463 464 ! This is either dipole or user, or both 465 is_ortho = .true. 466 if ( N_Elec == 1 ) then 467 ! When using a single electrode 468 ! one may always have an electric field. 469 ! For fields along the semi-infinite direction 470 ! it corresponds to a capacitor setup 471 ! For other fields it is simply a gate-like field. 472 ! However, when the electrode uses real space self-energies 473 ! then we do not allow *any* fields along the semi-inf 474 ! directions since that would mean an infinite field. 475 select case ( Elecs(1)%t_dir ) 476 case ( 4 , 5 , 6, 7 ) 477 is_ortho = Elec_ortho2semi_inf(Elecs(1), 478 & user_specified_field) 479 end select 480 else 481 ! For more than 1 electrode we do not allow any electric 482 ! fields along the semi-infinite directions. 483 do iE = 1 , N_Elec 484 is_ortho = is_ortho .and. 485 & Elec_ortho2semi_inf(Elecs(iE), user_specified_field) 486 end do 487 end if 488 if ( .not. is_ortho ) then 489 call die('User defined E-field must be 490 &perpendicular to semi-infinite directions.') 491 end if 492 end if ! acting_efield 493 494 ! We know that we currently allow people to do more than 495 ! they probably should be allowed. However, there are many 496 ! corner cases that may require dipole corrections, or 497 ! electric fields to "correct" an intrinsic dipole. 498 ! For instance, what should we do with a dipole in a transiesta 499 ! calculation? 500 ! Should we apply a field to counter act it in a device 501 ! calculation? 502 503 end if 504 505 frstme = .false. 506 507 call timer( 'DHSCF_Init', 2 ) 508#ifdef DEBUG 509 call write_debug( ' POS dhscf_init' ) 510#endif 511!------------------------------------------------------------------------- END 512 end subroutine dhscf_init 513 514 subroutine dhscf( nspin, norb, iaorb, iphorb, nuo, 515 & nuotot, nua, na, isa, xa, indxua, 516 & ntm, ifa, istr, iHmat, 517 & filesOut, maxnd, numd, 518 & listdptr, listd, Dscf, datm, maxnh, Hmat, 519 & Enaatm, Enascf, Uatm, Uscf, DUscf, DUext, 520 & Exc, Dxc, dipol, stress, Fal, stressl, 521 & use_rhog_in, charge_density_only ) 522 523C 524C Calculates the self-consistent field contributions to Hamiltonian 525C matrix elements, total energy and atomic forces. 526C Coded by J.M. Soler, August 1996. July 1997. 527C Modified by J.D. Gale, February 2000. 528C 529C ---------------------------------------------------------------------- 530C Input : 531C ---------------------------------------------------------------------- 532C integer nspin : Number of different spin polarisations 533C nspin=1 => Unpolarized, nspin=2 => polarized 534C nspin=4 => Noncollinear spin or spin-orbit 535C integer norb : Total number of basis orbitals in supercell 536C integer iaorb(norb) : Atom to which each orbital belongs 537C integer iphorb(norb) : Orbital index (within atom) of each orbital 538C integer nuo : Number of orbitals in a unit cell in this node 539C integer nuotot : Number of orbitals in a unit cell 540C integer nua : Number of atoms in unit cell 541C integer na : Number of atoms in supercell 542C integer isa(na) : Species index of all atoms in supercell 543C real*8 xa(3,na) : Atomic positions of all atoms in supercell 544C integer indxua : Index of equivalent atom in unit cell 545C integer ifa : Switch which fixes whether the SCF contrib. 546C to atomic forces is calculated and added to fa 547C integer istr : Switch which fixes whether the SCF contrib. 548C to stress is calculated and added to stress 549C integer iHmat : Switch which fixes whether the Hmat matrix 550C elements are calculated or not. 551C type(filesOut_t) filesOut : output file names (If blank => not saved) 552C integer maxnd : First dimension of listd and Dscf 553C integer numd(nuo) : Number of nonzero density-matrix 554C elements for each matrix row 555C integer listdptr(nuo) : Pointer to start of rows of density-matrix 556C integer listd(maxnd) : Nonzero-density-matrix-element column 557C indexes for each matrix row 558C real*8 Dscf(maxnd,h_spin_dim): SCF density-matrix elements 559C real*8 datm(norb) : Harris density-matrix diagonal elements 560C (atomic occupation charges of orbitals) 561C integer maxnh : First dimension of listh and Hmat 562C real*8 Hmat(maxnh,h_spin_dim) : Hamiltonian matrix in sparse form, 563C to which are added the matrix elements 564C <ORB_I | DeltaV | ORB_J>, where 565C DeltaV = Vna + Vxc(SCF) + 566C Vhartree(RhoSCF-RhoHarris) 567C ---------------------------------------------------------------------- 568C Input/output : 569C ---------------------------------------------------------------------- 570C integer ntm(3) : Number of mesh divisions of each cell 571C vector, including subgrid. 572C ---------------------------------------------------------------------- 573C Output : 574C ---------------------------------------------------------------------- 575C real*8 Enaatm : Integral of Vna * rhoatm 576C real*8 Enascf : Integral of Vna * rhoscf 577C real*8 Uatm : Harris hartree electron-interaction energy 578C real*8 Uscf : SCF hartree electron-interaction energy 579C real*8 DUscf : Electrostatic (Hartree) energy of 580C (rhoscf - rhoatm) density 581C real*8 DUext : Interaction energy with external electric field 582C real*8 Exc : SCF exchange-correlation energy 583C real*8 Dxc : SCF double-counting correction to Exc 584C Dxc = integral of ( (epsxc - Vxc) * Rho ) 585C All energies in Rydbergs 586C real*8 dipol(3): Electric dipole (in a.u.) 587C only when the system is a molecule/slab 588C ---------------------------------------------------------------------- 589C Input/output : 590C ---------------------------------------------------------------------- 591C real*8 Fal(3,nua) : Atomic forces, to which the SCF contribution 592C is added by this routine when ifa=1. 593C the SCF contribution is minus the derivative 594C of ( Enascf - Enaatm + DUscf + Exc ) with 595C respect to atomic positions, in Ry/Bohr. 596C contributions local to this node 597C real*8 stressl(3,3): Stress tensor, to which the SCF contribution 598C is added by this routine when ifa=1. 599C the SCF contribution is minus the derivative of 600C ( Enascf - Enaatm + DUscf + Exc ) / volume 601C with respect to the strain tensor, in Ry. 602C contributions local to this node 603C ---------------------------------------------------------------------- 604C Units : 605C ---------------------------------------------------------------------- 606C Energies in Rydbergs 607C Distances in Bohr 608C ---------------------------------------------------------------------- 609C Modules 610 use precision, only : dp, grid_p 611 612 use parallel, only : Node, Nodes 613 use atmfuncs, only : rcut, rcore 614 use units, only : Debye, eV, Ang 615 use fdf 616 use sys, only : die, bye 617 use mesh, only : nsm, nsp, meshLim 618 use parsing 619 use m_iorho, only : write_rho 620 use m_forhar, only : forhar 621 use alloc, only : re_alloc, de_alloc 622 use files, only : slabel 623 use files, only : filesOut_t ! derived type for output file names 624 use siesta_options, only : harrisfun, save_initial_charge_density 625 use siesta_options, only : analyze_charge_density_only 626 use meshsubs, only : LocalChargeOnMesh 627 use meshsubs, only : PartialCoreOnMesh 628 use meshsubs, only : NeutralAtomOnMesh 629 630 use moreMeshSubs, only : setMeshDistr, distMeshData 631 use moreMeshSubs, only : UNIFORM, QUADRATIC, LINEAR 632 use moreMeshSubs, only : TO_SEQUENTIAL, TO_CLUSTER, KEEP 633 use m_partial_charges, only: compute_partial_charges 634 use m_partial_charges, only: want_partial_charges 635 636 use siestaXC, only : cellXC ! Finds xc energy and potential 637 use siestaXC, only : myMeshBox ! Returns my processor mesh box 638 use siestaXC, only : jms_setMeshDistr => setMeshDistr 639 ! Sets a distribution of mesh 640 ! points over parallel processors 641 use m_vmat, only: vmat 642 use m_vmatsp, only: vmatsp 643 use m_rhoofd, only: rhoofd 644 use m_rhoofdsp, only: rhoofdsp 645#ifdef MPI 646 use mpi_siesta 647#endif 648 use iogrid_netcdf, only: write_grid_netcdf 649 use iogrid_netcdf, only: read_grid_netcdf 650 use siesta_options, only: read_charge_cdf 651 use siesta_options, only: savebader 652 use siesta_options, only: read_deformation_charge_cdf 653 use siesta_options, only: mix_charge 654 655 use m_efield, only: add_potential_from_field 656 use m_efield, only: user_specified_field, acting_efield 657 use dipole_m, only: get_dipole_origin 658 use dipole_m, only: get_field_from_dipole 659 use dipole_m, only: dipole_charge, dipole_potential 660 use dipole_m, only: dip_vacuum 661 use dipole_m, only: SLABDIPOLECORR 662 use dipole_m, only: SLABDIPOLECORR_NONE 663 use dipole_m, only: SLABDIPOLECORR_CHARGE, SLABDIPOLECORR_VACUUM 664 665 use m_doping_uniform, only: doping_active, doping_uniform 666 667 use m_charge_add, only: charge_add 668 use m_hartree_add, only: hartree_add 669 670#ifdef NCDF_4 671 use siesta_options, only: write_cdf 672 use m_ncdf_siesta, only: cdf_save_grid 673#endif 674 use m_rhofft, only: rhofft, FORWARD, BACKWARD 675 use m_rhog, only: rhog_in, rhog 676 use m_spin, only: spin 677 use m_spin, only: Spiral, qSpiral 678 679 use m_ts_global_vars,only: TSmode, TSrun 680 use m_ts_options, only: IsVolt, Elecs, N_elec 681 use m_ts_voltage, only: ts_voltage 682 use m_ts_hartree, only: ts_hartree_fix 683 684 implicit none 685 686 integer 687 & maxnd, maxnh, nua, na, norb, nspin, nuo, nuotot, 688 & iaorb(norb), ifa, iHmat, 689 & indxua(na), iphorb(norb), isa(na), 690 & istr, listd(*), listdptr(nuo), ntm(3), numd(nuo) 691 692 real(dp) 693 & datm(norb), dipol(3), Dscf(:,:), 694 & DUext, DUscf, Dxc, Enaatm, Enascf, Exc, 695 & Hmat(:,:), Uatm, Uscf, xa(3,na), 696 & stress(3,3), Fal(3,nua), stressl(3,3) 697 698 type(filesOut_t) filesOut 699 700 logical, intent(in), optional :: use_rhog_in 701 logical, intent(in), optional :: charge_density_only 702 703C ---------------------------------------------------------------------- 704C Routines called internally: 705C ---------------------------------------------------------------------- 706C cellxc(...) : Finds total exch-corr energy and potential 707C cross(a,b,c) : Finds the cross product of two vectors 708C dfscf(...) : Finds SCF contribution to atomic forces 709C dipole_charge : Finds electric dipole moment using the charges 710C doping(...) : Adds a background charge for doped systems 711C write_rho(...) : Saves electron density on a file 712C poison(...) : Solves Poisson equation 713C reord(...) : Reorders electron density and potential arrays 714C rhooda(...) : Finds Harris electron density in the mesh 715C rhoofd(...) : Finds SCF electron density in the mesh 716C rhoofdsp(...) : Finds SCF electron density in the mesh for 717C spiral arrangement of spins 718C timer(...) : Finds CPU times 719C vmat(...) : Finds matrix elements of SCF potential 720C vmatsp(...) : Finds matrix elements of SCF potential for 721C spiral arrangement of spins 722C delk(...) : Finds matrix elements of exp(i \vec{k} \cdot \vec{r}) 723C real*8 volcel( cell ) : Returns volume of unit cell 724C ---------------------------------------------------------------------- 725C Internal variables and arrays: 726C ---------------------------------------------------------------------- 727C integer nbcell : Number of bulk lattice vectors 728C real*8 bcell(3,3) : Bulk lattice vectors 729C real*8 cell(3,3) : Auxiliary lattice vectors (same as ucell) 730C real*8 const : Auxiliary variable (constant within a loop) 731C real*8 DEc : Auxiliary variable to call cellxc 732C real*8 DEx : Auxiliary variable to call cellxc 733C real*8 dvol : Mesh-cell volume 734C real*8 Ec : Correlation energy 735C real*8 Ex : Exchange energy 736C real*8 field(3) : External electric field 737C integer i : General-purpose index 738C integer ia : Atom index 739C integer io : Orbital index 740C integer ip : Point index 741C integer is : Species index 742C logical IsDiag : Is supercell diagonal? 743C integer ispin : Spin index 744C integer j : General-purpose index 745 746C integer myBox(2,3) : My processor's mesh box 747 748C integer nbcell : Number of independent bulk lattice vectors 749C integer npcc : Partial core corrections? (0=no, 1=yes) 750C integer nsd : Number of diagonal spin values (1 or 2) 751C integer ntpl : Number of mesh Total Points in unit cell 752C (including subpoints) locally 753C real*4 rhoatm(ntpl) : Harris electron density 754C real*4 rhopcc(ntpl) : Partial-core-correction density for xc 755C real*4 DRho(ntpl) : Selfconsistent electron density difference 756C real*8 rhotot : Total density at one point 757C real*8 rmax : Maximum orbital radius 758C real*8 scell(3,3) : Supercell vectors 759C real*4 Vaux(ntpl) : Auxiliary potential array 760C real*4 Vna(ntpl) : Sum of neutral-atom potentials 761C real*8 volume : Unit cell volume 762C real*4 Vscf(ntpl) : Hartree potential of selfconsistent density 763C real*8 x0(3) : Center of molecule 764C logical harrisfun : Harris functional or Kohn-Sham? 765C ---------------------------------------------------------------------- 766 767C Local variables 768 integer :: i, ia, ip, ispin, nsd, np_vac 769 real(dp) :: b1Xb2(3), const, DEc, DEx, DStres(3,3), 770 & Ec, Ex, rhotot, x0(3), volume, Vmax_vac, Vmean_vac 771 772 logical :: use_rhog 773 774 real(dp), external :: volcel, ddot 775 776 external 777 & cross, 778 & poison, 779 & reord, rhooda, 780 & timer 781 782! Dummy arrays for bsc_cellxc call 783 real(grid_p) :: aux3(3,1) 784 real(grid_p) :: dummy_DVxcdn(1,1,1) 785 real(grid_p), pointer :: Vscf_gga(:,:), DRho_gga(:,:) 786 external bsc_cellxc 787 788! Interface to JMS's SiestaXC 789 integer :: myBox(2,3) 790 real(dp) :: stressXC(3,3) 791 792C Work arrays 793 real(grid_p), pointer :: Vscf(:,:), Vscf_par(:,:), 794 & DRho(:,:), DRho_par(:,:), 795 & Vaux(:), Vaux_par(:), Chlocal(:), 796 & Totchar(:), fsrc(:), fdst(:), 797 & rhoatm_quad(:) => null(), 798 & DRho_quad(:,:) => null() 799 ! Temporary reciprocal spin quantity 800 real(grid_p) :: rnsd 801 802#ifdef MPI 803 integer :: MPIerror 804 real(dp) :: sbuffer(7), rbuffer(7) 805#endif 806 807#ifdef DEBUG 808 call write_debug( ' PRE DHSCF' ) 809#endif 810 811 if ( spin%H /= size(Dscf,dim=2) ) then 812 call die('Spin components is not equal to options.') 813 end if 814 815 if ( debug_dhscf ) then 816 write(*,debug_fmt) Node,'DM', 817 & (sqrt(sum(Dscf(:,ispin)**2)),ispin=1,spin%H) 818 write(*,debug_fmt) Node,'H', 819 & (sqrt(sum(Hmat(:,ispin)**2)),ispin=1,spin%H) 820 end if 821 822!-------------------------------------------------------------------- BEGIN 823C ---------------------------------------------------------------------- 824C Start of SCF iteration part 825C ---------------------------------------------------------------------- 826 827C ---------------------------------------------------------------------- 828C At the end of DHSCF_INIT, and also at the end of any previous 829! call to dhscf, we were in the UNIFORM distribution 830! Rhoatm, Rhopcc and Vna were in UNIFORM dist, sequential form 831! The index array endpht was in the QUADRATIC distribution 832C ---------------------------------------------------------------------- 833 834#ifdef _TRACE_ 835 call MPI_Barrier( MPI_Comm_World, MPIerror ) 836 call MPItrace_restart( ) 837#endif 838 call timer( 'DHSCF', 1 ) 839 call timer( 'DHSCF3', 1 ) 840 841 nullify( Vscf, Vscf_par, DRho, DRho_par, 842 & Vaux, Vaux_par, Chlocal, Totchar ) 843 nullify( Vscf_gga, DRho_gga) 844 845 volume = volcel(cell) 846 847C------------------------------------------------------------------------- 848 849 if (analyze_charge_density_only) then 850 !! Use the functionality in the first block 851 !! of the routine to get charge files and partial charges 852 call setup_analysis_options() 853 endif 854 855 if (filesOut%vna .ne. ' ') then 856 ! Uniform dist, sequential form 857#ifdef NCDF_4 858 if ( write_cdf ) then 859 call cdf_save_grid(trim(slabel)//'.nc','Vna',1,ntml,Vna) 860 else 861 call write_rho( filesOut%vna, 862 & cell, ntm, nsm, ntpl, 1, Vna) 863 call write_grid_netcdf( cell, ntm, 1, ntpl, Vna, "Vna") 864 end if 865#else 866 call write_rho( filesOut%vna, 867 & cell, ntm, nsm, ntpl, 1, Vna ) 868 call write_grid_netcdf( cell, ntm, 1, ntpl, Vna, "Vna") 869#endif 870 endif 871 872C Allocate memory for DRho using the UNIFORM data distribution 873 call re_alloc( DRho, 1, ntpl, 1, nspin, 'DRho','dhscf' ) 874 875C Find number of diagonal spin values 876 nsd = min( nspin, 2 ) 877 if ( nsd == 1 ) then 878 rnsd = 1._grid_p 879 else 880 rnsd = 1._grid_p / nsd 881 end if 882 883C ---------------------------------------------------------------------- 884C Find SCF electron density at mesh points. Store it in array DRho 885C ---------------------------------------------------------------------- 886! 887! The reading routine works in the uniform distribution, in 888! sequential form 889! 890 if (present(use_rhog_in)) then 891 use_rhog = use_rhog_in 892 else 893 use_rhog = .false. 894 endif 895 if (use_rhog) then 896 ! fourier transform back into drho 897 call rhofft(cell, ntml(1), ntml(2), ntml(3), ntm, nspin, 898 $ DRho,rhog_in,BACKWARD) 899 900 else if (read_charge_cdf) then 901 call read_grid_netcdf(ntm(1:3),nspin,ntpl,DRho,"Rho") 902 read_charge_cdf = .false. 903 else if (read_deformation_charge_cdf) then 904 call read_grid_netcdf(ntm(1:3),nspin,ntpl,DRho,"DeltaRho") 905 ! Add to diagonal components only 906 do ispin = 1,nsd 907 do ip= 1, ntpl 908C rhoatm and Drho are in sequential mode 909 DRho(ip,ispin) = DRho(ip,ispin) + rhoatm(ip) * rnsd 910 enddo 911 enddo 912 read_deformation_charge_cdf = .false. 913 else 914 ! Set the QUADRATIC distribution and allocate memory for DRho_par 915 ! since the construction of the density from the DM and orbital 916 ! data needs that distribution 917 if (nodes.gt.1) then 918 call setMeshDistr( QUADRATIC, nsm, nsp, 919 & nml, nmpl, ntml, ntpl ) 920 endif 921 call re_alloc( DRho_par, 1, ntpl, 1, nspin, 922 & 'DRho_par','dhscf' ) 923 if (Spiral) then 924 call rhoofdsp( norb, nml, nmpl, maxnd, numd, listdptr, listd, 925 & nspin, Dscf, DRho_par, nuo, nuotot, iaorb, 926 & iphorb, isa, qspiral ) 927 else 928 call rhoofd( norb, nmpl, maxnd, numd, listdptr, listd, 929 & nspin, Dscf, DRho_par, 930 & nuo, nuotot, iaorb, iphorb, isa ) 931 endif 932 ! DRHO_par is here in QUADRATIC, clustered form 933 934C Set the UNIFORM distribution again and copy DRho to it 935 if (nodes.gt.1) then 936 call setMeshDistr( UNIFORM, nsm, nsp, nml, nmpl, ntml, ntpl ) 937 endif 938 939 do ispin = 1, nspin 940 fsrc => DRho_par(:,ispin) 941 fdst => DRho(:,ispin) 942 ! Sequential to be able to write it out 943 ! if nodes==1, this call will just reorder 944 call distMeshData( QUADRATIC, fsrc, 945 & UNIFORM, fdst, TO_SEQUENTIAL ) 946 enddo 947 call de_alloc( DRho_par, 'DRho_par','dhscf' ) 948 949 if ( debug_dhscf ) then 950 write(*,debug_fmt) Node,'Rho', 951 & (sum(DRho(:,ispin))*dvol,ispin=1,nspin) 952 end if 953 954 if (save_initial_charge_density) then 955 ! This section is to be deprecated in favor 956 ! of "analyze_charge_density_only" 957 ! (except for the special name for the .nc file) 958#ifdef NCDF_4 959 if ( write_cdf ) then 960 call cdf_save_grid(trim(slabel)//'.nc','RhoInit',nspin, 961 & ntml,DRho) 962 else 963 call write_rho( "RHO_INIT", cell, ntm, nsm, ntpl, 964 $ nspin, DRho) 965 call write_grid_netcdf( cell, ntm, nspin, ntpl, DRho, 966 $ "RhoInit") 967 end if 968#else 969 call write_rho( "RHO_INIT", cell, ntm, nsm, ntpl, 970 $ nspin, DRho) 971 call write_grid_netcdf( cell, ntm, nspin, ntpl, DRho, 972 $ "RhoInit") 973#endif 974 call timer('DHSCF3',2) 975 call timer('DHSCF',2) 976 call bye("STOP after producing RHO_INIT from input DM") 977 endif 978 979 endif 980 981 if (mix_charge) then 982 ! Save fourier transform of charge density 983 call rhofft(cell, ntml(1), ntml(2), ntml(3), ntm, nspin, 984 $ DRho,rhog,FORWARD) 985 endif 986! 987! Proper place to integrate Hirshfeld and Voronoi code, 988! since we have just computed rhoatm and Rho. 989 990 if (want_partial_charges) then 991 ! The endpht array is in the quadratic distribution, so 992 ! we need to use it for this... 993 if (nodes.gt.1) then 994 call setMeshDistr( QUADRATIC, nsm, nsp, 995 & nml, nmpl, ntml, ntpl ) 996 endif 997 call re_alloc( DRho_quad, 1, ntpl, 1, nspin, 998 & 'DRho_quad','dhscf' ) 999 call re_alloc( rhoatm_quad, 1, ntpl, 1000 & 'rhoatm_quad','dhscf' ) 1001 ! Redistribute grid-density 1002 do ispin = 1, nspin 1003 fsrc => DRho(:,ispin) 1004 fdst => DRho_quad(:,ispin) 1005 ! if nodes==1, this call will just reorder 1006 call distMeshData( UNIFORM, fsrc, 1007 & QUADRATIC, fdst, TO_CLUSTER ) 1008 enddo 1009 call distMeshData( UNIFORM, rhoatm, 1010 & QUADRATIC, rhoatm_quad, TO_CLUSTER ) 1011 1012 call compute_partial_charges(DRho_quad,rhoatm_quad, 1013 . nspin, iaorb, iphorb, 1014 . isa, nmpl,dvol) 1015 1016 call de_alloc(rhoatm_quad,'rhoatm_quad','dhscf') 1017 call de_alloc(Drho_quad,'DRho_quad','dhscf') 1018 if (nodes.gt.1) then 1019 call setMeshDistr( UNIFORM, nsm, nsp, 1020 & nml, nmpl, ntml, ntpl ) 1021 endif 1022 endif 1023 1024C ---------------------------------------------------------------------- 1025C Save electron density 1026C ---------------------------------------------------------------------- 1027 if (filesOut%rho .ne. ' ') then 1028 ! DRho is already using a uniform, sequential form 1029#ifdef NCDF_4 1030 if ( write_cdf ) then 1031 call cdf_save_grid(trim(slabel)//'.nc','Rho',nspin,ntml, 1032 & DRho) 1033 else 1034 call write_rho( filesOut%rho, cell, ntm, nsm, ntpl, nspin, 1035 & DRho ) 1036 call write_grid_netcdf( cell, ntm, nspin, ntpl, DRho, "Rho") 1037 end if 1038#else 1039 call write_rho( filesOut%rho, cell, ntm, nsm, ntpl, nspin, 1040 & DRho ) 1041 call write_grid_netcdf( cell, ntm, nspin, ntpl, DRho, "Rho") 1042#endif 1043 endif 1044C ---------------------------------------------------------------------- 1045C Save the diffuse ionic charge and/or the total (ionic+electronic) charge 1046C ---------------------------------------------------------------------- 1047 if (filesOut%psch .ne. ' ' .or. filesOut%toch .ne. ' ') then 1048C Find diffuse ionic charge on mesh 1049 ! Note that the *OnMesh routines, except PhiOnMesh, 1050 ! work with any distribution, thanks to the fact that 1051 ! the ipa, idop, and indexp arrays are distro-specific 1052 call re_alloc( Chlocal, 1, ntpl, 'Chlocal', 'dhscf' ) 1053 call LocalChargeOnMesh( na, isa, ntpl, Chlocal, indxua ) 1054 ! Chlocal comes out in clustered form, so we convert it 1055 call reord( Chlocal, Chlocal, nml, nsm, TO_SEQUENTIAL ) 1056 1057 if ( debug_dhscf ) then 1058 write(*,debug_fmt) Node,'Chlocal',sum(Chlocal)*dvol 1059 end if 1060 1061C Save diffuse ionic charge 1062 if (filesOut%psch .ne. ' ') then 1063#ifdef NCDF_4 1064 if ( write_cdf ) then 1065 call cdf_save_grid(trim(slabel)//'.nc','Chlocal',1,ntml, 1066 & Chlocal) 1067 else 1068 call write_rho( filesOut%psch, cell, ntm, nsm, ntpl, 1, 1069 & Chlocal) 1070 call write_grid_netcdf( cell, ntm, 1, ntpl, Chlocal, 1071 & 'Chlocal' ) 1072 end if 1073#else 1074 call write_rho( filesOut%psch, cell, ntm, nsm, ntpl, 1, 1075 & Chlocal) 1076 call write_grid_netcdf( cell, ntm, 1, ntpl, 1077 & Chlocal, 'Chlocal' ) 1078#endif 1079 endif 1080 1081C Save total (ionic+electronic) charge 1082 if ( filesOut%toch .ne. ' ') then 1083 ! ***************** 1084 ! ** IMPORTANT ** 1085 ! The Chlocal array is re-used to minimize memory 1086 ! usage. In the this small snippet the Chlocal 1087 ! array will contain the total charge, and 1088 ! if the logic should change, (i.e. should Chlocal 1089 ! be retained) is the Totchar needed to be re-instantiated. 1090 ! ***************** 1091 1092!$OMP parallel default(shared), private(ispin,ip) 1093 do ispin = 1, nsd 1094!$OMP do 1095 do ip = 1, ntpl 1096 Chlocal(ip) = Chlocal(ip) + DRho(ip,ispin) 1097 end do 1098!$OMP end do 1099 end do 1100!$OMP end parallel 1101 1102 ! See note above 1103#ifdef NCDF_4 1104 if ( write_cdf ) then 1105 call cdf_save_grid(trim(slabel)//'.nc','RhoTot',1,ntml, 1106 & Chlocal) 1107 else 1108 call write_rho(filesOut%toch,cell,ntm,nsm,ntpl,1,Chlocal) 1109 call write_grid_netcdf( cell, ntm, 1, ntpl, Chlocal, 1110 & "TotalCharge") 1111 end if 1112#else 1113 call write_rho( filesOut%toch, cell, ntm, nsm, ntpl, 1, 1114 & Chlocal ) 1115 call write_grid_netcdf( cell, ntm, 1, ntpl, 1116 & Chlocal, "TotalCharge") 1117#endif 1118 end if 1119 call de_alloc( Chlocal, 'Chlocal', 'dhscf' ) 1120 endif 1121 1122C ---------------------------------------------------------------------- 1123C Save the total charge (model core + valence) for Bader analysis 1124C ---------------------------------------------------------------------- 1125 1126 ! The test for toch guarantees that we are in "analysis mode" 1127 if (filesOut%toch .ne. ' ' .and. savebader) then 1128 call save_bader_charge() 1129 endif 1130 1131C Find difference between selfconsistent and atomic densities 1132 1133 !Both DRho and rhoatm are using a UNIFORM, sequential form 1134!$OMP parallel default(shared), private(ispin,ip) 1135 do ispin = 1,nsd 1136!$OMP do 1137 do ip = 1,ntpl 1138 DRho(ip,ispin) = DRho(ip,ispin) - rhoatm(ip) * rnsd 1139 enddo 1140!$OMP end do nowait 1141 enddo 1142!$OMP end parallel 1143 1144C ---------------------------------------------------------------------- 1145C Save electron density difference 1146C ---------------------------------------------------------------------- 1147 if (filesOut%drho .ne. ' ') then 1148#ifdef NCDF_4 1149 if ( write_cdf ) then 1150 call cdf_save_grid(trim(slabel)//'.nc','RhoDelta',nspin,ntml, 1151 & DRho) 1152 else 1153 call write_rho( filesOut%drho, cell, ntm, nsm, ntpl, nspin, 1154 & DRho ) 1155 call write_grid_netcdf( cell, ntm, nspin, ntpl, DRho, 1156 & "DeltaRho") 1157 end if 1158#else 1159 call write_rho( filesOut%drho, cell, ntm, nsm, ntpl, nspin, 1160 & DRho ) 1161 call write_grid_netcdf( cell, ntm, nspin, ntpl, 1162 & DRho, "DeltaRho") 1163#endif 1164 endif 1165 1166 if (present(charge_density_only)) then 1167 if (charge_density_only) then 1168 call timer('DHSCF3',2) 1169 call timer('DHSCF',2) 1170 call de_alloc( DRho, 'DRho', 'dhscf' ) 1171 RETURN 1172 endif 1173 endif 1174 1175! End of analysis section 1176! Can exit now, if requested 1177 1178 if (analyze_charge_density_only) then 1179 call timer('DHSCF3',2) 1180 call timer('DHSCF',2) 1181 call bye("STOP after analyzing charge from input DM") 1182 endif 1183 1184!------------------------------------------------------------- 1185C Transform spin density into sum and difference 1186 1187 ! TODO Check for NC/SO: 1188 ! Should we diagonalize locally at every point first? 1189 if (nsd .eq. 2) then 1190!$OMP parallel do default(shared), private(rhotot,ip) 1191 do ip = 1,ntpl 1192 rhotot = DRho(ip,1) + DRho(ip,2) 1193 DRho(ip,2) = DRho(ip,2) - DRho(ip,1) 1194 DRho(ip,1) = rhotot 1195 enddo 1196!$OMP end parallel do 1197 endif 1198 1199 if ( debug_dhscf ) then 1200 write(*,debug_fmt) Node,'DRho-clean', 1201 & (sum(DRho(:,ispin))*dvol,ispin=1,nspin) 1202 end if 1203 1204C Add a background charge to neutralize the net charge, to 1205C model doped systems. It only adds the charge at points 1206C where there are atoms (i.e., not in vacuum). 1207C First, call with 'task=0' to add background charge 1208 if (doping_active) call doping_uniform(cell,ntpl,0, 1209 $ DRho(:,1),rhoatm) 1210 1211 ! Add doping in cell (from ChargeGeometries/Geometry.Charge) 1212 ! Note that this routine will return immediately if no dopant is present 1213 call charge_add('+',cell, ntpl, DRho(:,1) ) 1214 1215 if ( debug_dhscf ) then 1216 write(*,debug_fmt) Node,'DRho-doped', 1217 & (sum(DRho(:,ispin))*dvol,ispin=1,nspin) 1218 end if 1219 1220C ---------------------------------------------------------------------- 1221C Calculate the dipole moment 1222C ---------------------------------------------------------------------- 1223 dipol(1:3) = 0.0_dp 1224 if ( nbcell /= 3 ) then 1225 ! Find dipole in system 1226 1227 select case ( SLABDIPOLECORR ) 1228 case ( SLABDIPOLECORR_VACUUM ) 1229 ! we should not calculate the dipole here 1230 ! User requests the vacuum dipole calculation 1231 case default 1232 ! Regardless of whether the user requests a 1233 ! dipole correction or not, then we will calculate the dipole 1234 ! using the charge method. 1235 1236 ! Get dipole origin 1237 call get_dipole_origin(cell, nua, xa, x0) 1238 1239 ! This routine is distribution-blind 1240 ! and will reduce over all processors. 1241 call dipole_charge(cell, dvol, ntm, ntml, nsm, DRho, 1242 & x0, dipol) 1243 1244 ! Orthogonalize dipole to bulk directions 1245 if ( nbcell == 1 ) then ! chain 1246 const = ddot(3,dipol,1,bcell,1) / ddot(3,bcell,1,bcell,1) 1247 dipol(1:3) = dipol(1:3) - const * bcell(1:3,1) 1248 else if ( nbcell == 2 ) then ! slab 1249 call cross( bcell(1,1), bcell(1,2), b1Xb2 ) 1250 const = ddot(3,dipol,1,b1Xb2,1) / ddot(3,b1Xb2,1,b1Xb2,1) 1251 dipol(1:3) = const * b1Xb2(1:3) 1252 end if 1253 1254 end select 1255 1256 endif 1257 1258C ---------------------------------------------------------------------- 1259C Find Hartree potential of DRho = rhoscf-rhoatm. Store it in Vaux 1260C ---------------------------------------------------------------------- 1261! Solve Poisson's equation 1262 call re_alloc( Vaux, 1, ntpl, 'Vaux', 'dhscf' ) 1263 1264 call poison( cell, ntml(1), ntml(2), ntml(3), ntm, DRho, 1265 & DUscf, Vaux, DStres, nsm ) 1266 1267 if ( debug_dhscf ) then 1268 write(*,debug_fmt) Node,'Poisson',sqrt(sum(Vaux(:)**2)) 1269 end if 1270 1271 ! Vscf is in the UNIFORM, sequential form, and only using 1272 ! the first spin index 1273 1274 ! We require that even the SIESTA potential is "fixed" 1275 ! NOTE, this will only do something if 1276 ! TS.Hartree.Fix is set 1277 call ts_hartree_fix( ntm, ntml, Vaux) 1278 1279C Add contribution to stress from electrostatic energy of rhoscf-rhoatm 1280 if (istr .eq. 1) then 1281 stressl(1:3,1:3) = stressl(1:3,1:3) + DStres(1:3,1:3) 1282 endif 1283 1284C ---------------------------------------------------------------------- 1285C Find electrostatic (Hartree) energy of full SCF electron density 1286C using the original data distribution 1287C ---------------------------------------------------------------------- 1288 Uatm = Uharrs 1289 Uscf = 0._dp 1290!$OMP parallel do default(shared), private(ip), 1291!$OMP&reduction(+:Uscf) 1292 do ip = 1, ntpl 1293 Uscf = Uscf + Vaux(ip) * rhoatm(ip) 1294 enddo 1295!$OMP end parallel do 1296 Uscf = Uscf * dVol + Uatm + DUscf 1297 1298C Call doping with 'task=1' to remove background charge added previously 1299C The extra charge thus only affects the Hartree energy and potential, 1300C but not the contribution to Enascf ( = \Int_{Vna*\rho}) 1301 if (doping_active) call doping_uniform(cell,ntpl,1, 1302 $ DRho(:,1),rhoatm) 1303 1304 ! Remove doping in cell (from ChargeGeometries/Geometry.Charge) 1305 ! Note that this routine will return immediately if no dopant is present 1306 call charge_add('-',cell, ntpl, DRho(:,1) ) 1307 1308C ---------------------------------------------------------------------- 1309C Add neutral-atom potential to Vaux 1310C ---------------------------------------------------------------------- 1311 Enaatm = 0.0_dp 1312 Enascf = 0.0_dp 1313!$OMP parallel do default(shared), private(ip), 1314!$OMP&reduction(+:Enaatm,Enascf) 1315 do ip = 1, ntpl 1316 Enaatm = Enaatm + Vna(ip) * rhoatm(ip) 1317 Enascf = Enascf + Vna(ip) * DRho(ip,1) 1318 Vaux(ip) = Vaux(ip) + Vna(ip) 1319 enddo 1320!$OMP end parallel do 1321 Enaatm = Enaatm * dVol 1322 Enascf = Enaatm + Enascf * dVol 1323 1324 1325 if ( SLABDIPOLECORR == SLABDIPOLECORR_VACUUM ) then 1326 ! In case the user requested a vacuum dipole, calculate it here 1327 ! REMARK 1328 ! This method will only ever work for 3D-FFT poisson solved 1329 ! potential 1330 ! Note we do this on the potential which will be written 1331 ! to ElectrostaticPotential 1332 call dipole_potential(cell, nua, isa, xa, ntm, ntml, nsm, 1333 & dip_vacuum, Vaux, dipol) 1334 1335 ! Orthogonalize dipole to bulk directions 1336 if ( nbcell == 1 ) then ! chain 1337 const = ddot(3,dipol,1,bcell,1) / ddot(3,bcell,1,bcell,1) 1338 dipol(1:3) = dipol(1:3) - const * bcell(1:3,1) 1339 else if ( nbcell == 2 ) then ! slab 1340 call cross( bcell(1,1), bcell(1,2), b1Xb2 ) 1341 const = ddot(3,dipol,1,b1Xb2,1) / ddot(3,b1Xb2,1,b1Xb2,1) 1342 dipol(1:3) = const * b1Xb2(1:3) 1343 end if 1344 1345 end if 1346 1347C ---------------------------------------------------------------------- 1348C Add potential from external electric field (if present) 1349C ---------------------------------------------------------------------- 1350 if ( acting_efield ) then 1351 select case ( SLABDIPOLECORR ) 1352 case ( SLABDIPOLECORR_NONE ) 1353 ! No dipole correction 1354 field(:) = 0._dp 1355 DUext = 0._dp 1356 1357 case default 1358 ! some kind of dipole correction 1359 field(:) = get_field_from_dipole(dipol, cell) 1360 1361 if ( Node == 0 ) then 1362 write(6,'(a,3(tr1,f12.4),a)') 1363 $ 'Dipole moment in unit cell =', dipol/Debye, ' D' 1364 write(6,'(a,3(tr1,f12.6),a)') 1365 $ 'Electric field for dipole correction =', 1366 $ field/eV*Ang, ' eV/Ang/e' 1367 end if 1368 1369 ! The dipole correction energy has an extra factor 1370 ! of one half because the field involved is internal. 1371 ! See the paper by Bengtsson DOI:10.1103/PhysRevB.59.12301 1372 ! Hence we compute this part separately 1373 DUext = -0.5_dp * ddot(3,field,1,dipol,1) 1374 end select 1375 1376 ! Add the external electric field 1377 field(:) = field(:) + user_specified_field(:) 1378 1379 ! This routine expects a sequential array, 1380 ! but it is distribution-blind 1381 ! Note that we do not need to specify a center 1382 ! Since we *just* need the saw-tooth to be located 1383 ! in the vacuum region (beyond orbitals) 1384 call add_potential_from_field( field, cell, nua, isa, xa, 1385 & ntm, ntml, nsm, Vaux ) 1386 1387 ! Add energy of external electric field 1388 DUext = DUext - ddot(3,user_specified_field,1,dipol,1) 1389 1390 end if 1391 1392! --------------------------------------------------------------------- 1393! Transiesta: 1394! add the potential corresponding to the (possible) voltage-drop. 1395! note that ts_voltage is not sharing the reord wih efield since 1396! we should not encounter both at the same time. 1397! --------------------------------------------------------------------- 1398 if (TSmode.and.IsVolt.and.TSrun) then 1399 ! This routine expects a sequential array, 1400 ! in whatever distribution 1401 call ts_voltage(cell, ntm, ntml, Vaux) 1402 endif 1403 1404! ---------------------------------------------------------------------- 1405! Add potential from user defined geometries (if present) 1406! ---------------------------------------------------------------------- 1407 call hartree_add( cell, ntpl, Vaux ) 1408 1409C ---------------------------------------------------------------------- 1410C Save electrostatic potential 1411C ---------------------------------------------------------------------- 1412 if (filesOut%vh .ne. ' ') then 1413 ! Note that only the first spin component is used 1414#ifdef NCDF_4 1415 if ( write_cdf ) then 1416 call cdf_save_grid(trim(slabel)//'.nc','Vh',1,ntml, 1417 & Vaux) 1418 else 1419 call write_rho( filesOut%vh, cell, ntm, nsm, ntpl, 1, Vaux ) 1420 call write_grid_netcdf( cell, ntm, 1, ntpl, Vaux, 1421 & "ElectrostaticPotential") 1422 end if 1423#else 1424 call write_rho( filesOut%vh, cell, ntm, nsm, ntpl, 1, Vaux ) 1425 call write_grid_netcdf( cell, ntm, 1, ntpl, 1426 & Vaux, "ElectrostaticPotential") 1427#endif 1428 endif 1429 1430C Get back spin density from sum and difference 1431 ! TODO Check for NC/SO: 1432 ! Should we diagonalize locally at every point first? 1433 if (nsd .eq. 2) then 1434!$OMP parallel do default(shared), private(ip,rhotot) 1435 do ip = 1, ntpl 1436 rhotot = DRho(ip,1) 1437 DRho(ip,1) = 0.5_dp * (rhotot - DRho(ip,2)) 1438 DRho(ip,2) = 0.5_dp * (rhotot + DRho(ip,2)) 1439 enddo 1440!$OMP end parallel do 1441 endif 1442 1443C Exchange-correlation energy 1444C ---------------------------------------------------------------------- 1445 call re_alloc( Vscf, 1, ntpl, 1, nspin, 'Vscf', 'dhscf' ) 1446 1447 if (npcc .eq. 1) then 1448 1449!$OMP parallel default(shared), private(ip,ispin) 1450 do ispin = 1,nsd 1451!$OMP do 1452 do ip= 1, ntpl 1453 DRho(ip,ispin) = DRho(ip,ispin) + 1454 & (rhopcc(ip)+rhoatm(ip)) * rnsd 1455 enddo 1456!$OMP end do nowait 1457 enddo 1458!$OMP end parallel 1459 1460 else 1461 1462!$OMP parallel default(shared), private(ip,ispin) 1463 do ispin = 1,nsd 1464!$OMP do 1465 do ip= 1, ntpl 1466 DRho(ip,ispin) = DRho(ip,ispin) + 1467 & rhoatm(ip) * rnsd 1468 enddo 1469!$OMP end do nowait 1470 enddo 1471!$OMP end parallel 1472 1473 end if 1474 1475 ! Write the electron density used by cellxc 1476 if (filesOut%rhoxc .ne. ' ') then 1477#ifdef NCDF_4 1478 if ( write_cdf ) then 1479 call cdf_save_grid(trim(slabel)//'.nc','RhoXC',nspin,ntml, 1480 & DRho ) 1481 else 1482 call write_rho( filesOut%rhoxc, cell, ntm, nsm, ntpl, nspin, 1483 & DRho ) 1484 call write_grid_netcdf( cell, ntm, nspin, ntpl, DRho, 1485 & "RhoXC") 1486 end if 1487#else 1488 call write_rho( filesOut%rhoxc, cell, ntm, nsm, ntpl, nspin, 1489 & DRho ) 1490 call write_grid_netcdf( cell, ntm, nspin, ntpl, DRho, 1491 & "RhoXC") 1492#endif 1493 endif 1494 1495! Everything now is in UNIFORM, sequential form 1496 1497 call timer("XC",1) 1498 1499 ! Switch to "zero/not-zero rho" distribution (miscalled 'LINEAR') 1500 if (nodes.gt.1) then 1501 call setMeshDistr( LINEAR, nsm, nsp, nml, nmpl, ntml, ntpl ) 1502 endif 1503 1504 call re_alloc( Vscf_gga, 1, ntpl, 1, nspin, 'Vscf_gga','dhscf') 1505 call re_alloc( DRho_gga, 1, ntpl, 1, nspin, 'DRho_gga','dhscf') 1506 1507 ! Redistribute all spin densities 1508 do ispin = 1, nspin 1509 fsrc => DRho(:,ispin) 1510 fdst => DRho_gga(:,ispin) 1511 call distMeshData( UNIFORM, fsrc, LINEAR, fdst, KEEP ) 1512 enddo 1513 1514 if (use_bsc_cellxc) then 1515 1516 call timer("BSC-CellXC",1) 1517 call bsc_cellxc( 0, 0, cell, ntml, ntml, ntpl, 0, aux3, nspin, 1518 & DRho_gga, Ex, Ec, DEx, DEc, Vscf_gga, 1519 & dummy_DVxcdn, stressl ) 1520 1521 call timer("BSC-CellXC",2) 1522 1523 else 1524 1525 call timer("GXC-CellXC",1) 1526 1527 ! Note that RG's meshLim is in blocks of nsm points, and 1-based 1528 ! Example. 1529 ! 1:16 17:32 (base 1) 1530 ! 0:15 16:31 (base 0) 1531 ! Times nsm=2: 1532 ! 0:30 32:62 1533 ! Add 1 to lb, and nsm to ub: 1534 ! 1:32 33:64 (base 1, in units of small-points) 1535 1536 myBox(1,:) = (meshLim(1,:)-1)*nsm + 1 1537 myBox(2,:) = (meshLim(2,:)-1)*nsm + nsm 1538 1539 call cellXC( 0, cell, ntm, myBox(1,1), myBox(2,1), 1540 & myBox(1,2), myBox(2,2), 1541 & myBox(1,3), myBox(2,3), nspin, 1542 & DRho_gga, Ex, Ec, DEx, DEc, stressXC, Vscf_gga, 1543 & keep_input_distribution = .true. ) 1544 ! Vscf is still sequential after the call to JMS's cellxc 1545 stress = stress + stressXC 1546 call timer("GXC-CellXC",2) 1547 1548 endif ! use_bsc_cellxc 1549 1550 ! Go back to uniform distribution 1551 if (nodes.gt.1) then 1552 call setMeshDistr( UNIFORM, nsm, nsp, nml, nmpl, ntml, ntpl) 1553 endif 1554 1555 ! Redistribute to the Vxc array 1556 do ispin = 1,nspin 1557 fsrc => Vscf_gga(:,ispin) 1558 fdst => Vscf(:,ispin) 1559 call distMeshData( LINEAR, fsrc, UNIFORM, fdst, KEEP ) 1560 enddo 1561 call de_alloc( DRho_gga, 'DRho_gga', 'dhscf' ) 1562 call de_alloc( Vscf_gga, 'Vscf_gga', 'dhscf' ) 1563 1564 call timer("XC",2) 1565 1566 if ( debug_dhscf ) then 1567 write(*,debug_fmt) Node,'XC', 1568 & (sqrt(sum(Vscf(:,ispin)**2)),ispin=1,nspin) 1569 end if 1570 1571 Exc = Ex + Ec 1572 Dxc = DEx + DEc 1573 1574! Vscf contains only Vxc, and is UNIFORM and sequential 1575! Now we add up the other contributions to it, at 1576! the same time that we get DRho back to true DeltaRho form 1577!$OMP parallel default(shared), private(ip,ispin) 1578 1579 ! Hartree potential only has diagonal components 1580 do ispin = 1,nsd 1581 if (npcc .eq. 1) then 1582!$OMP do 1583 do ip = 1,ntpl 1584 DRho(ip,ispin) = DRho(ip,ispin) - 1585 & (rhoatm(ip)+rhopcc(ip)) * rnsd 1586 Vscf(ip,ispin) = Vscf(ip,ispin) + Vaux(ip) 1587 enddo 1588!$OMP end do 1589 else 1590!$OMP do 1591 do ip = 1,ntpl 1592 DRho(ip,ispin) = DRho(ip,ispin) - rhoatm(ip) * rnsd 1593 Vscf(ip,ispin) = Vscf(ip,ispin) + Vaux(ip) 1594 enddo 1595!$OMP end do 1596 endif 1597 enddo 1598!$OMP end parallel 1599 1600C ---------------------------------------------------------------------- 1601C Save total potential 1602C ---------------------------------------------------------------------- 1603 if (filesOut%vt .ne. ' ') then 1604#ifdef NCDF_4 1605 if ( write_cdf ) then 1606 call cdf_save_grid(trim(slabel)//'.nc','Vt',nspin,ntml, 1607 & Vscf) 1608 else 1609 call write_rho( filesOut%vt, cell, ntm, nsm, ntpl, nspin, 1610 & Vscf ) 1611 call write_grid_netcdf( cell, ntm, nspin, ntpl, Vscf, 1612 & "TotalPotential") 1613 end if 1614#else 1615 call write_rho( filesOut%vt, cell, ntm, nsm, ntpl, nspin, Vscf) 1616 call write_grid_netcdf( cell, ntm, nspin, ntpl, 1617 & Vscf, "TotalPotential") 1618#endif 1619 endif 1620 1621C ---------------------------------------------------------------------- 1622C Print vacuum level 1623C ---------------------------------------------------------------------- 1624 1625 if (filesOut%vt/=' ' .or. filesOut%vh/=' ') then 1626 forall(ispin=1:nsd) 1627 . DRho(:,ispin) = DRho(:,ispin) + rhoatm(:) * rnsd 1628 call vacuum_level( ntpl, nspin, DRho, Vscf, 1629 . np_vac, Vmax_vac, Vmean_vac ) 1630 forall(ispin=1:nsd) 1631 . DRho(:,ispin) = DRho(:,ispin) - rhoatm(:) * rnsd 1632 if (np_vac>0 .and. Node==0) print'(/,a,2f12.6,a)', 1633 . 'dhscf: Vacuum level (max, mean) =', 1634 . Vmax_vac/eV, Vmean_vac/eV, ' eV' 1635 endif 1636 1637C ---------------------------------------------------------------------- 1638C Find SCF contribution to hamiltonian matrix elements 1639C ---------------------------------------------------------------------- 1640 if (iHmat .eq. 1) then 1641 if (nodes.gt.1) then 1642 call setMeshDistr( QUADRATIC, nsm, nsp, 1643 & nml, nmpl, ntml, ntpl ) 1644 endif 1645 1646! This is a work array, to which we copy Vscf 1647 call re_alloc( Vscf_par, 1, ntpl, 1, nspin, 1648 & 'Vscf_par', 'dhscf' ) 1649 1650 do ispin = 1, nspin 1651 fsrc => Vscf(:,ispin) 1652 fdst => Vscf_par(:,ispin) 1653 call distMeshData( UNIFORM, fsrc, 1654 & QUADRATIC, fdst, TO_CLUSTER ) 1655 enddo 1656 1657 if (Spiral) then 1658 call vmatsp( norb, nml, nmpl, dvol, nspin, Vscf_par, maxnd, 1659 & numd, listdptr, listd, Hmat, nuo, 1660 & nuotot, iaorb, iphorb, isa, qspiral ) 1661 else 1662 call vmat( norb, nmpl, dvol, spin, Vscf_par, maxnd, 1663 & numd, listdptr, listd, Hmat, nuo, 1664 & nuotot, iaorb, iphorb, isa ) 1665 endif 1666 1667 call de_alloc( Vscf_par, 'Vscf_par', 'dhscf' ) 1668 if (nodes.gt.1) then 1669! Everything back to UNIFORM, sequential 1670 call setMeshDistr( UNIFORM, nsm, nsp, 1671 & nml, nmpl, ntml, ntpl ) 1672 endif 1673 endif 1674 1675 1676#ifdef MPI 1677C Global reduction of Uscf/DUscf/Uatm/Enaatm/Enascf 1678 sbuffer(1) = Uscf 1679 sbuffer(2) = DUscf 1680 sbuffer(3) = Uatm 1681 sbuffer(4) = Enaatm 1682 sbuffer(5) = Enascf 1683 if (use_bsc_cellxc) then 1684 sbuffer(6) = Exc 1685 sbuffer(7) = Dxc 1686 endif 1687 1688 call MPI_AllReduce( sbuffer, rbuffer, 7, MPI_double_precision, 1689 & MPI_Sum, MPI_Comm_World, MPIerror ) 1690 Uscf = rbuffer(1) 1691 DUscf = rbuffer(2) 1692 Uatm = rbuffer(3) 1693 Enaatm = rbuffer(4) 1694 Enascf = rbuffer(5) 1695 if (use_bsc_cellxc) then 1696 Exc = rbuffer(6) 1697 Dxc = rbuffer(7) 1698 endif 1699#endif /* MPI */ 1700 1701C Add contribution to stress from the derivative of the Jacobian of --- 1702C r->r' (strained r) in the integral of Vna*(rhoscf-rhoatm) 1703 if (istr .eq. 1) then 1704 do i = 1,3 1705 stress(i,i) = stress(i,i) + ( Enascf - Enaatm ) / volume 1706 enddo 1707 endif 1708 1709C Stop time counter for SCF iteration part 1710 call timer( 'DHSCF3', 2 ) 1711 1712C ---------------------------------------------------------------------- 1713C End of SCF iteration part 1714C ---------------------------------------------------------------------- 1715 1716 if (ifa.eq.1 .or. istr.eq.1) then 1717C ---------------------------------------------------------------------- 1718C Forces and stress : SCF contribution 1719C ---------------------------------------------------------------------- 1720C Start time counter for force calculation part 1721 call timer( 'DHSCF4', 1 ) 1722 1723C Find contribution of partial-core-correction 1724 if (npcc .eq. 1) then 1725 call reord( rhopcc, rhopcc, nml, nsm, TO_CLUSTER ) 1726 call reord( Vaux, Vaux, nml, nsm, TO_CLUSTER ) 1727 ! The partial core calculation only acts on 1728 ! the diagonal spin-components (no need to 1729 ! redistribute un-used elements) 1730 do ispin = 1, nsd 1731 call reord( Vscf(:,ispin), Vscf(:,ispin), 1732 & nml, nsm, TO_CLUSTER ) 1733 enddo 1734 1735 call PartialCoreOnMesh( na, isa, ntpl, rhopcc, indxua, nsd, 1736 & dvol, volume, Vscf, Vaux, Fal, 1737 & stressl, ifa.ne.0, istr.ne.0 ) 1738 1739 call reord( rhopcc, rhopcc, nml, nsm, TO_SEQUENTIAL ) 1740 call reord( Vaux, Vaux, nml, nsm, TO_SEQUENTIAL ) 1741 ! ** see above 1742 do ispin = 1, nsd 1743 call reord( Vscf(:,ispin), Vscf(:,ispin), 1744 & nml, nsm, TO_SEQUENTIAL ) 1745 enddo 1746 1747 if ( debug_dhscf ) then 1748 write(*,debug_fmt) Node,'PartialCore', 1749 & (sqrt(sum(Vscf(:,ispin)**2)),ispin=1,nsd) 1750 end if 1751 endif 1752 1753 if ( harrisfun) then 1754! Forhar deals internally with its own needs 1755! for distribution changes 1756! Upon entry, everything is UNIFORM, sequential form 1757 call forhar( ntpl, nspin, nml, ntml, ntm, npcc, cell, 1758 & rhoatm, rhopcc, Vna, DRho, Vscf, Vaux ) 1759! Upon return, everything is UNIFORM, sequential form 1760 endif 1761 1762C Transform spin density into sum and difference 1763 ! TODO NC/SO 1764 ! Should we perform local diagonalization? 1765 if (nsd .eq. 2) then 1766!$OMP parallel do default(shared), private(rhotot,ip) 1767 do ip = 1,ntpl 1768 rhotot = DRho(ip,1) + DRho(ip,2) 1769 DRho(ip,2) = DRho(ip,2) - DRho(ip,1) 1770 DRho(ip,1) = rhotot 1771 enddo 1772!$OMP end parallel do 1773 endif 1774 1775C Find contribution of neutral-atom potential 1776 call reord( Vna, Vna, nml, nsm, TO_CLUSTER ) 1777 call reord( DRho, DRho, nml, nsm, TO_CLUSTER ) 1778 call NeutralAtomOnMesh( na, isa, ntpl, Vna, indxua, dvol, 1779 & volume, DRho, Fal, stressl, 1780 & ifa.ne.0, istr.ne.0 ) 1781 call reord( DRho, DRho, nml, nsm, TO_SEQUENTIAL ) 1782 call reord( Vna, Vna, nml, nsm, TO_SEQUENTIAL ) 1783 1784 if (nodes.gt.1) then 1785 call setMeshDistr( QUADRATIC, nsm, nsp, 1786 & nml, nmpl, ntml, ntpl ) 1787 endif 1788 1789 call re_alloc( Vscf_par, 1, ntpl, 1, nspin, 1790 & 'Vscf_par', 'dhscf' ) 1791 do ispin = 1, nspin 1792 fsrc => Vscf(:,ispin) 1793 fdst => Vscf_par(:,ispin) 1794 call distMeshData( UNIFORM, fsrc, 1795 & QUADRATIC, fdst, TO_CLUSTER ) 1796 enddo 1797 1798! Remember that Vaux contains everything except Vxc 1799 call re_alloc( Vaux_par, 1, ntpl, 'Vaux_par', 'dhscf' ) 1800 call distMeshData( UNIFORM, Vaux, 1801 & QUADRATIC, Vaux_par, TO_CLUSTER ) 1802 1803 call dfscf( ifa, istr, na, norb, nuo, nuotot, nmpl, nspin, 1804 & indxua, isa, iaorb, iphorb, 1805 & maxnd, numd, listdptr, listd, Dscf, datm, 1806 & Vscf_par, Vaux_par, dvol, volume, Fal, stressl ) 1807 1808 call de_alloc( Vaux_par, 'Vaux_par', 'dhscf' ) 1809 call de_alloc( Vscf_par, 'Vscf_par', 'dhscf' ) 1810 1811 if (nodes.gt.1) then 1812 call setMeshDistr( UNIFORM, nsm, nsp, nml, nmpl, ntml, ntpl ) 1813 endif 1814 1815 1816C Stop time counter for force calculation part 1817 call timer( 'DHSCF4', 2 ) 1818C ---------------------------------------------------------------------- 1819C End of force and stress calculation 1820C ---------------------------------------------------------------------- 1821 endif 1822 1823! We are in the UNIFORM distribution 1824! Rhoatm, Rhopcc and Vna are in UNIFORM dist, sequential form 1825! The index array endpht is in the QUADRATIC distribution 1826 1827C Stop time counter 1828 call timer( 'DHSCF', 2 ) 1829 1830C ---------------------------------------------------------------------- 1831C Free locally allocated memory 1832C ---------------------------------------------------------------------- 1833 call de_alloc( Vaux, 'Vaux', 'dhscf' ) 1834 call de_alloc( Vscf, 'Vscf', 'dhscf' ) 1835 call de_alloc( DRho, 'DRho', 'dhscf' ) 1836 1837#ifdef DEBUG 1838 call write_debug( ' POS DHSCF' ) 1839#endif 1840!------------------------------------------------------------------------ END 1841 CONTAINS 1842 1843 subroutine save_bader_charge() 1844 use meshsubs, only: ModelCoreChargeOnMesh 1845#ifdef NCDF_4 1846 use siesta_options, only: write_cdf 1847 use m_ncdf_siesta, only: cdf_save_grid 1848#endif 1849 ! Auxiliary routine to output the Bader Charge 1850 ! 1851 real(grid_p), pointer :: BaderCharge(:) => null() 1852 1853 call re_alloc( BaderCharge, 1, ntpl, name='BaderCharge', 1854 & routine='dhscf' ) 1855 1856 ! Find a model core charge by re-scaling the local charge 1857 call ModelCoreChargeOnMesh( na, isa, ntpl, BaderCharge, indxua ) 1858 ! It comes out in clustered form, so we convert it 1859 call reord( BaderCharge, BaderCharge, nml, nsm, TO_SEQUENTIAL) 1860 do ispin = 1,nsd 1861 BaderCharge(1:ntpl) = BaderCharge(1:ntpl) + DRho(1:ntpl,ispin) 1862 enddo 1863#ifdef NCDF_4 1864 if ( write_cdf ) then 1865 call cdf_save_grid(trim(slabel)//'.nc','RhoBader',1,ntml, 1866 & BaderCharge) 1867 else 1868 call write_rho( trim(slabel)// ".BADER", cell, 1869 $ ntm, nsm, ntpl, 1, BaderCharge ) 1870 call write_grid_netcdf( cell, ntm, 1, ntpl, 1871 $ BaderCharge, "BaderCharge") 1872 end if 1873#else 1874 call write_rho( trim(slabel)// ".BADER", cell, 1875 $ ntm, nsm, ntpl, 1, BaderCharge ) 1876 call write_grid_netcdf( cell, ntm, 1, ntpl, 1877 $ BaderCharge, "BaderCharge") 1878#endif 1879 1880 call de_alloc( BaderCharge, name='BaderCharge' ) 1881 end subroutine save_bader_charge 1882 1883 subroutine setup_analysis_options() 1884 !! For the analyze-charge-density-only case, 1885 !! avoiding any diagonalization 1886 1887 use siesta_options, only: hirshpop, voropop 1888 use siesta_options, only: saverho, savedrho, saverhoxc 1889 use siesta_options, only: savevh, savevt, savevna 1890 use siesta_options, only: savepsch, savetoch 1891 1892 want_partial_charges = (hirshpop .or. voropop) 1893 1894 if (saverho) filesOut%rho = trim(slabel)//'.RHO' 1895 if (savedrho) filesOut%drho = trim(slabel)//'.DRHO' 1896 if (saverhoxc) filesOut%rhoxc = trim(slabel)//'.RHOXC' 1897 if (savevh) filesOut%vh = trim(slabel)//'.VH' 1898 if (savevt) filesOut%vt = trim(slabel)//'.VT' 1899 if (savevna) filesOut%vna = trim(slabel)//'.VNA' 1900 if (savepsch) filesOut%psch = trim(slabel)//'.IOCH' 1901 if (savetoch) filesOut%toch = trim(slabel)//'.TOCH' 1902 1903 end subroutine setup_analysis_options 1904 1905 end subroutine dhscf 1906 1907 subroutine delk_wrapper(isigneikr, norb, maxnd, 1908 & numd, listdptr, listd, 1909 & nuo, nuotot, iaorb, iphorb, isa ) 1910 1911 use m_delk, only : delk ! The real workhorse, similar to vmat 1912 1913 use moreMeshSubs, only : setMeshDistr 1914 use moreMeshSubs, only : UNIFORM, QUADRATIC 1915 use parallel, only : Nodes 1916 use mesh, only : nsm, nsp 1917 1918! 1919! This is a wrapper to call delk, using some of the module 1920! variables of m_dhscf, but from outside dhscf itself. 1921! 1922 integer :: isigneikr, 1923 & norb, nuo, nuotot, maxnd, 1924 & iaorb(*), iphorb(*), isa(*), 1925 & numd(nuo), 1926 & listdptr(nuo), listd(maxnd) 1927 1928! The dhscf module variables used are: 1929! nmpl 1930! dvol 1931! nml 1932! nmpl 1933! ntml 1934! ntpl 1935! 1936! Some of them might be put somewhere else (mesh?) to allow some 1937! of the kitchen-sink functionality of dhscf to be made more modular. 1938! For example, this wrapper might live independently if enough mesh 1939! information is made available to it. 1940 1941C ---------------------------------------------------------------------- 1942C Calculate matrix elements of exp(i \vec{k} \cdot \vec{r}) 1943C ---------------------------------------------------------------------- 1944 if (isigneikr .eq. 1 .or. isigneikr .eq. -1) then 1945 1946 if (nodes.gt.1) then 1947 call setMeshDistr( QUADRATIC, nsm, nsp, 1948 & nml, nmpl, ntml, ntpl ) 1949 endif 1950 1951 call delk( isigneikr, norb, nmpl, dvol, maxnd, 1952 & numd, listdptr, listd, 1953 & nuo, nuotot, iaorb, iphorb, isa ) 1954 1955 if (nodes.gt.1) then 1956! Everything back to UNIFORM, sequential 1957 call setMeshDistr( UNIFORM, nsm, nsp, 1958 & nml, nmpl, ntml, ntpl ) 1959 endif 1960 1961 endif 1962 end subroutine delk_wrapper 1963 1964 !> Chech whether we need to initialize stencils for 1965 !> BSC's version of cellxc, and whether we have a vdw 1966 !> functional 1967 subroutine xc_setup(need_grads,is_vdw) 1968 use siestaxc, only: getxc 1969 1970 logical, intent(out) :: need_grads 1971 logical, intent(out) :: is_vdw 1972 1973 integer :: nf, nXCfunc 1974 character(len=20):: XCauth(10), XCfunc(10) 1975 1976 need_grads = .false. 1977 is_vdw = .false. 1978 call getXC( nXCfunc, XCfunc, XCauth ) 1979 do nf = 1,nXCfunc 1980 if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then 1981 need_grads = .true. 1982 endif 1983 if ( XCfunc(nf).eq.'VDW' .or. XCfunc(nf).eq.'vdw') then 1984 need_grads = .true. 1985 is_vdw = .true. 1986 endif 1987 enddo 1988 end subroutine xc_setup 1989 1990 end module m_dhscf 1991