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 module m_forhar 9 10 implicit none 11 12 public :: forhar 13 private 14 15 CONTAINS 16 subroutine forhar( NTPL, NSPIN, NML, NTML, NTM, NPCC, 17 $ CELL, RHOATM, 18 & RHOPCC, VNA, DRHOOUT, VHARRIS1, VHARRIS2 ) 19 20C ********************************************************************** 21C Build the potentials needed for computing Harris forces: 22C (V_NA + V_Hartree(DeltaRho_in) - DV_xc(Rho_in)/Dn * (Rho_out-Rho_in)) 23C (V_NA + V_Hartree(DeltaRho_in) + V_xc(Rho_in) 24C **** BEHAVIOUR ******************************************************* 25C In the first SCF step, V_Hartree(DeltaRho_in) is zero, because 26C in that case, Rho_SCF(r) = Rho_atm(r) and therefore, DeltaRho(r) = 0 27C This calculation will be skipped. 28C NOTE: This is true only if the initial DM is built from atomic charges... 29C If Harris + Spin polarized in the first SCF step, then Vharris2 will 30C multiply to D Rho(Harris)/D R inside dfscf, and the change of 31C the harris density respect the displacement 32C on one atom will not depend on spin. We add in Vharris2 the 33C contributions of both spins. 34C Coded by J. Junquera 09/00 35C ********************************************************************** 36 37 use precision, only : dp, grid_p 38 use alloc, only : re_alloc, de_alloc 39 40 use parallel, only : nodes 41 use mesh, only : NSM, nsp, meshLim 42 use siestaXC, only : cellXC ! Finds xc energy and potential 43 44 use moreMeshSubs, only : setMeshDistr, distMeshData 45 use moreMeshSubs, only: UNIFORM, LINEAR 46 use moreMeshSubs, only: KEEP 47 48 use fdf, only: fdf_get 49 50 INTEGER :: NTPL, NML(3), NTML(3) 51 INTEGER, INTENT(IN) :: NSPIN, NPCC, NTM(3) 52 53 REAL(dp), INTENT(IN) :: CELL(3,3) 54 REAL(grid_p), INTENT(IN) :: VNA(NTPL), RHOATM(NTPL), 55 & RHOPCC(NTPL) 56 REAL(grid_p), INTENT(IN) :: DRHOOUT(NTPL,NSPIN) 57 REAL(grid_p), TARGET, INTENT(OUT) :: VHARRIS1(NTPL,NSPIN) 58 REAL(grid_p), INTENT(OUT) :: VHARRIS2(NTPL) 59 60 EXTERNAL bsc_cellxc 61 62! AG: Note: REAL*4 variables are really REAL(kind=grid_p) 63! 64C ***** INPUT ********************************************************** 65C INTEGER NTPL : Number of Mesh Total Points in unit cell 66C (including subpoints) locally. 67C INTEGER NSPIN : Spin polarizations 68C INTEGER NTM(3) : Number of mesh divisions of each cell 69C vector, including subgrid 70C INTEGER NPCC : Partial core corrections? (0=no,1=yes) 71C REAL*8 CELL(3,3) : Cell vectors 72C REAL*4 RHOATM(NTPL) : Harris density at mesh points 73C REAL*4 RHOPCC(NTPL) : Partial-core-correction density for xc 74C REAL*4 VNA(NTPL) : Sum of neutral atoms potentials 75C REAL*4 DRHOOUT(NTPL,NSPIN) : Charge density at the mesh points 76C in current step. 77C The charge density that enters in forhar 78C is Drho_out-Rhoatm. 79C ***** OUTPUT ********************************************************* 80C REAL*4 VHARRIS1(NTPL,NSPIN) : Vna + V_Hartree(DeltaRho_in) + V_xc(Rho_in) 81C REAL*4 VHARRIS2(NTPL) : Vna + V_Hartree(DeltaRho_in) + 82C - DV_xc(Rho_in)/DRho_in * (Rho_out-Rho_in) 83C If Harris forces are computed in the 84C first SCF step, it does not depend on spin 85C ***** INTERNAL VARIABLES ********************************************* 86C REAL*4 DVXDN(NTPL,NSPIN,NSPIN): Derivative of exchange-correlation 87C potential respect the charge density 88C ********************************************************************** 89 90C ---------------------------------------------------------------------- 91C Internal variables and arrays 92C ---------------------------------------------------------------------- 93 94 INTEGER IP, ISPIN, ISPIN2, myBox(2,3), NMPL 95 REAL(dp) EX, EC, DEX, DEC, STRESS(3,3) 96 97 real(grid_p) :: aux3(3,1) !! dummy arrays for cellxc 98 real(grid_p), pointer :: drhoin(:,:), drhoin_par(:,:), 99 & dvxcdn(:,:,:), dvxcdn_par(:,:,:), 100 & vharris1_par(:,:), fsrc(:), fdst(:) 101 logical :: use_bsc_cellxc 102 103 nullify( drhoin, dvxcdn ) 104 call re_alloc( drhoin, 1, ntpl, 1, nspin, 'drhoin', 'forhar' ) 105 call re_alloc( dvxcdn, 1, ntpl, 1, nspin, 1, nspin, 106 & 'dvxcdn', 'forhar' ) 107 108C ---------------------------------------------------------------------- 109C Initialize some variables 110C ---------------------------------------------------------------------- 111 VHARRIS1(:,:) = 0.0_grid_p 112 VHARRIS2(:) = 0.0_grid_p 113 DRHOIN(:,:) = 0.0_grid_p 114 DVXCDN(:,:,:) = 0.0_grid_p 115 116C ---------------------------------------------------------------------- 117C Compute exchange-correlation energy and potential and 118C their derivatives respect the input charge, that is, Harris charge 119C or the sum of atomic charges. 120C ---------------------------------------------------------------------- 121 122 ! All these arrays are in the UNIFORM distribution,in SEQUENTIAL form 123 DO ISPIN = 1, NSPIN 124 DRHOIN(1:NTPL,ISPIN) = RHOATM(1:NTPL)/NSPIN 125 IF (NPCC .EQ. 1) 126 . DRHOIN(1:NTPL,ISPIN) = DRHOIN(1:NTPL,ISPIN) + 127 . RHOPCC(1:NTPL)/NSPIN 128 ENDDO 129 130 ! Give the opportunity to use BSC's version 131 use_bsc_cellxc = fdf_get("XC.Use.BSC.Cellxc",.false.) 132 133 ! The input distribution is UNIFORM, but we need to work with the 134 ! "zero/not-zero rho" distribution (miscalled 'LINEAR') 135 if (nodes.gt.1) then 136 call setMeshDistr( LINEAR, nsm, nsp, nml, nmpl, ntml, ntpl ) 137 endif 138 139 nullify( drhoin_par, vharris1_par, dvxcdn_par ) 140 call re_alloc( drhoin_par, 1, ntpl, 1, nspin, 141 & 'drhoin_par', 'forhar' ) 142 call re_alloc( vharris1_par, 1, ntpl, 1, nspin, 143 & 'vharris1_par', 'forhar' ) 144 call re_alloc( dvxcdn_par, 1, ntpl, 1, nspin, 1, nspin, 145 & 'dvxcdn_par', 'forhar' ) 146 147 DO ISPIN = 1, NSPIN 148 fsrc => drhoin(:,ispin) 149 fdst => drhoin_par(:,ispin) 150 call distMeshData( UNIFORM, fsrc, LINEAR, fdst, KEEP ) 151 ENDDO 152 153 if (use_bsc_cellxc) then 154 155 STRESS(:,:) = 0.0_dp 156 CALL bsc_cellxc( 0, 1, CELL, NTML, NTML, NTPL, 0, AUX3, NSPIN, 157 & DRHOIN_PAR, EX, EC, DEX, DEC, VHARRIS1_PAR, 158 & DVXCDN_PAR, STRESS ) 159 160 else 161 162 myBox(1,:) = (meshLim(1,:)-1)*nsm + 1 163 myBox(2,:) = (meshLim(2,:)-1)*nsm + nsm 164 165 CALL CELLXC( 0, CELL, NTM, myBox(1,1), myBox(2,1), 166 . myBox(1,2), myBox(2,2), 167 . myBox(1,3), myBox(2,3), NSPIN, DRHOIN_par, 168 . EX, EC, DEX, DEC, STRESS, VHARRIS1_par, DVXCDN_par, 169 & keep_input_distribution = .true. ) 170 171 endif 172 173 if (nodes.gt.1) then 174! Everything back to UNIFORM, sequential 175 call setMeshDistr( UNIFORM, nsm, nsp, 176 & nml, nmpl, ntml, ntpl ) 177 endif 178 179 DO ISPIN = 1, NSPIN 180 fsrc => VHARRIS1_PAR(:,ISPIN) 181 fdst => VHARRIS1(:,ispin) 182 call distMeshData( LINEAR, fsrc, UNIFORM, fdst, KEEP ) 183 DO ISPIN2 = 1, NSPIN 184 fsrc => DVXCDN_PAR(:,ISPIN,ISPIN2) 185 fdst => DVXCDN(:,ISPIN,ISPIN2) 186 call distMeshData( LINEAR, fsrc, UNIFORM, fdst, KEEP ) 187 ENDDO 188 ENDDO 189 call de_alloc( dvxcdn_par, 'dvxcdn_par', 'forhar' ) 190 call de_alloc( vharris1_par, 'vharris1_par', 'forhar' ) 191 call de_alloc( drhoin_par, 'drhoin_par', 'forhar' ) 192 193 194 DO ISPIN = 1, NSPIN 195 DO IP = 1, NTPL 196 VHARRIS1(IP,ISPIN) = VHARRIS1(IP,ISPIN) + VNA(IP) 197 ENDDO 198 ENDDO 199 200C ---------------------------------------------------------------------- 201C Compute the product DV_xc(Rho_in)/DRho_in * (Rho_out - Rho_in). 202C Since the charge that enters into forhar is DRHOOUT = Rho_out-Rhoatm 203C no extra transformation on the charge density is needed. 204C ---------------------------------------------------------------------- 205 206 DO ISPIN = 1, NSPIN 207 DO ISPIN2 = 1, NSPIN 208 DO IP = 1, NTPL 209 VHARRIS2(IP) = VHARRIS2(IP) + 210 & DVXCDN(IP,ISPIN2,ISPIN) * DRHOOUT(IP,ISPIN2) 211 ENDDO 212 ENDDO 213 ENDDO 214 215C ---------------------------------------------------------------------- 216C Since V_Hartree(DeltaRho_in) = 0.0, we only add to vharris2 the neutral 217C atom potential (note sign change of above intermediate result) 218C ---------------------------------------------------------------------- 219 DO IP = 1, NTPL 220 VHARRIS2(IP) = VNA(IP) - VHARRIS2(IP) 221 ENDDO 222 223 call de_alloc( dvxcdn, 'dvxcdn', 'forhar') 224 call de_alloc( drhoin, 'drhoin', 'forhar') 225 end subroutine forhar 226 227 end module m_forhar 228