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