Octopus
magnetic.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module magnetic_oct_m
23 use comm_oct_m
24 use debug_oct_m
27 use global_oct_m
28 use grid_oct_m
29 use ions_oct_m
32 use mesh_oct_m
34 use mpi_oct_m
43 use unit_oct_m
45
46 implicit none
47
48 private
49 public :: &
58
59contains
60
61 ! ---------------------------------------------------------
62 subroutine magnetic_density(mesh, std, rho, md)
63 class(mesh_t), intent(in) :: mesh
64 type(states_elec_dim_t), intent(in) :: std
65 real(real64), intent(in) :: rho(:,:)
66 real(real64), intent(out) :: md(:,:)
67
68 push_sub(magnetic_density)
69
70 select case (std%ispin)
71 case (unpolarized)
72 md = m_zero
73
74 case (spin_polarized)
75 md = m_zero
76 md(1:mesh%np, 3) = rho(1:mesh%np, 1) - rho(1:mesh%np, 2)
77
78 case (spinors)
79 md(1:mesh%np, 1) = m_two*rho(1:mesh%np, 3)
80 md(1:mesh%np, 2) = -m_two*rho(1:mesh%np, 4)
81 md(1:mesh%np, 3) = rho(1:mesh%np, 1) - rho(1:mesh%np, 2)
82 end select
83
84 pop_sub(magnetic_density)
85 end subroutine magnetic_density
86
87
88 ! ---------------------------------------------------------
89 subroutine magnetic_moment(mesh, st, rho, mm)
90 class(mesh_t), intent(in) :: mesh
91 type(states_elec_t), intent(in) :: st
92 real(real64), intent(in) :: rho(:,:)
93 real(real64), intent(out) :: mm(3)
94
95 real(real64), allocatable :: md(:,:)
96
97 push_sub(states_elec_magnetic_moment)
98
99 safe_allocate(md(1:mesh%np, 1:3))
100 call magnetic_density(mesh, st%d, rho, md)
101
102 mm(1) = dmf_integrate(mesh, md(:, 1), reduce = .false.)
103 mm(2) = dmf_integrate(mesh, md(:, 2), reduce = .false.)
104 mm(3) = dmf_integrate(mesh, md(:, 3), reduce = .false.)
105
106 if (mesh%parallel_in_domains) then
107 call mesh%allreduce(mm)
108 end if
109
110 safe_deallocate_a(md)
111
112 pop_sub(states_elec_magnetic_moment)
113 end subroutine magnetic_moment
115
116 ! ---------------------------------------------------------
117 subroutine write_magnetic_moments(mesh, st, ions, boundaries, lmm_r, iunit, namespace)
118 class(mesh_t), intent(in) :: mesh
119 type(states_elec_t), intent(in) :: st
120 type(ions_t), intent(in) :: ions
121 type(boundaries_t), intent(in) :: boundaries
122 real(real64), intent(in) :: lmm_r
123 integer, optional, intent(in) :: iunit
124 type(namespace_t), optional, intent(in) :: namespace
125
126 integer :: ia
127 real(real64) :: mm(max(mesh%box%dim, 3))
128 real(real64), allocatable :: lmm(:,:)
129
130 push_sub(write_magnetic_moments)
131
132 call magnetic_moment(mesh, st, st%rho, mm)
133 safe_allocate(lmm(1:max(mesh%box%dim, 3), 1:ions%natoms))
134 call magnetic_local_moments(mesh, st, ions, boundaries, st%rho, lmm_r, lmm)
135
136 message(1) = 'Total Magnetic Moment:'
137 call messages_info(1, iunit=iunit, namespace=namespace)
138 if (st%d%ispin == spin_polarized) then ! collinear spin
139 write(message(1), '(a,f10.6)') ' mz = ', mm(3)
140 call messages_info(1, iunit=iunit, namespace=namespace)
141 else if (st%d%ispin == spinors) then ! non-collinear
142 write(message(1), '(1x,3(a,f10.6,3x))') 'mx = ', mm(1),'my = ', mm(2),'mz = ', mm(3)
143 call messages_info(1, iunit=iunit, namespace=namespace)
144 end if
145
146 write(message(1), '(a,a,a,f7.3,a)') 'Local Magnetic Moments (sphere radius [', &
147 trim(units_abbrev(units_out%length)),'] = ', units_from_atomic(units_out%length, lmm_r), '):'
148 call messages_info(1, iunit=iunit, namespace=namespace)
149 if (st%d%ispin == spin_polarized) then ! collinear spin
150 write(message(1),'(a,6x,14x,a)') ' Ion','mz'
151 call messages_info(1, iunit=iunit, namespace=namespace)
152 do ia = 1, ions%natoms
153 write(message(1),'(i4,a10,f15.6)') ia, trim(ions%atom(ia)%species%get_label()), lmm(3, ia)
154 call messages_info(1, iunit=iunit, namespace=namespace)
155 end do
156 else if (st%d%ispin == spinors) then ! non-collinear
157 write(message(1),'(a,8x,13x,a,13x,a,13x,a)') ' Ion','mx','my','mz'
158 call messages_info(1, iunit=iunit, namespace=namespace)
159 do ia = 1, ions%natoms
160 write(message(1),'(i4,a10,9f15.6)') ia, trim(ions%atom(ia)%species%get_label()), lmm(:, ia)
161 call messages_info(1, iunit=iunit, namespace=namespace)
162 end do
163 end if
164
165 safe_deallocate_a(lmm)
166
168 end subroutine write_magnetic_moments
169
170 ! ---------------------------------------------------------
171 subroutine magnetic_local_moments(mesh, st, ions, boundaries, rho, rr, lmm)
172 class(mesh_t), intent(in) :: mesh
173 type(states_elec_t), intent(in) :: st
174 type(ions_t), intent(in) :: ions
175 type(boundaries_t), intent(in) :: boundaries
176 real(real64), intent(in) :: rho(:,:)
177 real(real64), intent(in) :: rr
178 real(real64), intent(out) :: lmm(max(mesh%box%dim, 3), ions%natoms)
179
180 integer :: ia, idir, is
181 real(real64), allocatable :: md(:, :)
182 type(submesh_t) :: sphere
183 real(real64) :: cosqr, sinqr
184 complex(real64), allocatable :: phase_spiral(:)
185
186 push_sub(magnetic_local_moments)
187
188 safe_allocate(md(1:mesh%np, 1:max(mesh%box%dim, 3)))
189
190 call magnetic_density(mesh, st%d, rho, md)
191 lmm = m_zero
192 do ia = 1, ions%natoms
193 call submesh_init(sphere, ions%space, mesh, ions%latt, ions%pos(:, ia), rr)
194
195 if (boundaries%spiral) then
196 safe_allocate(phase_spiral(1:sphere%np))
197 do is = 1, sphere%np
198 phase_spiral(is) = exp(+m_zi*sum((sphere%rel_x(:,is) + sphere%center - mesh%x(sphere%map(is),:)) &
199 *boundaries%spiral_q(1:mesh%box%dim)))
200 end do
201
202 if (mesh%box%dim>= 3) then
203 lmm(1,ia) = m_zero
204 lmm(2,ia) = m_zero
205
206 do is = 1, sphere%np
207 !There is a factor of 1/2 in phase_spiral
208 cosqr = real(phase_spiral(is), real64)
209 sinqr = aimag(phase_spiral(is))
210 lmm(1,ia) = lmm(1,ia)+md(sphere%map(is),1)*cosqr - md(sphere%map(is),2)*sinqr
211 lmm(2,ia) = lmm(2,ia)+md(sphere%map(is),1)*sinqr + md(sphere%map(is),2)*cosqr
212 end do
213 lmm(1,ia) = lmm(1,ia)*mesh%volume_element
214 lmm(2,ia) = lmm(2,ia)*mesh%volume_element
215 lmm(3,ia) = dsm_integrate_frommesh(mesh, sphere, md(1:mesh%np,3), reduce = .false.)
216 else
217 assert(.not. boundaries%spiral)
218 end if
219
220 safe_deallocate_a(phase_spiral)
221 else
222 do idir = 1, max(mesh%box%dim, 3)
223 lmm(idir, ia) = dsm_integrate_frommesh(mesh, sphere, md(1:mesh%np,idir), reduce = .false.)
224 end do
225 end if
226
227 call submesh_end(sphere)
228 end do
229
230 if (mesh%parallel_in_domains) then
231 call mesh%allreduce(lmm)
232 end if
233
234 safe_deallocate_a(md)
235
237 end subroutine magnetic_local_moments
238
239 ! ---------------------------------------------------------
240 subroutine magnetic_total_magnetization(mesh, st, qq, trans_mag)
241 class(mesh_t), intent(in) :: mesh
242 type(states_elec_t), intent(in) :: st
243 real(real64), intent(in) :: qq(:)
244 complex(real64), intent(out) :: trans_mag(6)
245
246 integer :: ip
247 complex(real64), allocatable :: tmp(:,:)
248 real(real64), allocatable :: md(:, :)
249 complex(real64) :: expqr
250
252
253 call profiling_in("TOTAL_MAGNETIZATION")
254
255 safe_allocate(tmp(1:mesh%np, 1:6))
256 safe_allocate(md(1:mesh%np, 1:max(mesh%box%dim, 3)))
257
258 call magnetic_density(mesh, st%d, st%rho, md)
259 do ip = 1, mesh%np
260 expqr = exp(-m_zi*sum(mesh%x(ip, 1:mesh%box%dim)*qq(1:mesh%box%dim)))
261 tmp(ip,1) = expqr*md(ip,1)
262 tmp(ip,2) = expqr*md(ip,2)
263 tmp(ip,3) = expqr*md(ip,3)
264 tmp(ip,4) = conjg(expqr)*md(ip,1)
265 tmp(ip,5) = conjg(expqr)*md(ip,2)
266 tmp(ip,6) = conjg(expqr)*md(ip,3)
267 end do
268
269 do ip = 1, 6
270 trans_mag(ip) = zmf_integrate(mesh, tmp(:,ip), reduce=.false.)
271 end do
272 if (mesh%parallel_in_domains) then
273 call mesh%allreduce(trans_mag)
274 end if
275
276
277 safe_deallocate_a(md)
278 safe_deallocate_a(tmp)
279
280 call profiling_out("TOTAL_MAGNETIZATION")
281
283 end subroutine magnetic_total_magnetization
284
285 ! ---------------------------------------------------------
289 subroutine magnetic_induced(namespace, gr, st, psolver, kpoints, a_ind, b_ind)
290 type(namespace_t), intent(in) :: namespace
291 type(grid_t), intent(in) :: gr
292 type(states_elec_t), intent(inout) :: st
293 type(poisson_t), intent(in) :: psolver
294 type(kpoints_t), intent(in) :: kpoints
295 real(real64), contiguous, intent(out) :: a_ind(:, :)
296 real(real64), contiguous, intent(out) :: b_ind(:, :)
297
299
300 integer :: idir
301 real(real64), allocatable :: jj(:, :, :)
302
303 push_sub(magnetic_induced)
304
305 ! If the states are real, we should never have reached here, but
306 ! just in case we return zero.
307 if (states_are_real(st)) then
308 a_ind = m_zero
309 b_ind = m_zero
310 pop_sub(magnetic_induced)
311 return
312 end if
313
314 safe_allocate(jj(1:gr%np_part, 1:gr%der%dim, 1:st%d%nspin))
315 call states_elec_calc_quantities(gr, st, kpoints, .false., paramagnetic_current = jj)
316
317 !We sum the current for up and down, valid for collinear and noncollinear spins
318 if (st%d%nspin > 1) then
319 do idir = 1, gr%der%dim
320 jj(:, idir, 1) = jj(:, idir, 1) + jj(:, idir, 2)
321 end do
322 end if
323
324 a_ind = m_zero
325 do idir = 1, gr%der%dim
326 call dpoisson_solve(psolver, namespace, a_ind(:, idir), jj(:, idir, 1))
327 end do
328 ! This minus sign is introduced here because the current that has been used
329 ! before is the "number-current density", and not the "charge-current density",
330 ! and therefore there is a minus sign missing (electrons are negative charges...)
331 a_ind(1:gr%np, 1:gr%der%dim) = - a_ind(1:gr%np, 1:gr%der%dim) / p_c
332
333 call dderivatives_curl(gr%der, a_ind, b_ind)
334
335 safe_deallocate_a(jj)
336 pop_sub(magnetic_induced)
337 end subroutine magnetic_induced
338
339 subroutine write_total_xc_torque(iunit, mesh, vxc, st)
340 integer, intent(in) :: iunit
341 class(mesh_t), intent(in) :: mesh
342 real(real64), intent(in) :: vxc(:,:)
343 type(states_elec_t), intent(in) :: st
344
345 real(real64), allocatable :: torque(:,:)
346 real(real64) :: tt(3)
347
348 push_sub(write_total_xc_torque)
349
350 safe_allocate(torque(1:mesh%np, 1:3))
351
352 call calc_xc_torque(mesh, vxc, st, torque)
353
354 tt(1) = dmf_integrate(mesh, torque(:,1))
355 tt(2) = dmf_integrate(mesh, torque(:,2))
356 tt(3) = dmf_integrate(mesh, torque(:,3))
357
358 if (mpi_grp_is_root(mpi_world)) then
359 write(iunit, '(a)') 'Total xc torque:'
360 write(iunit, '(1x,3(a,es10.3,3x))') 'Tx = ', tt(1),'Ty = ', tt(2),'Tz = ', tt(3)
361 end if
362
363 safe_deallocate_a(torque)
364
365 pop_sub(write_total_xc_torque)
366 end subroutine write_total_xc_torque
367
368 ! ---------------------------------------------------------
369 subroutine calc_xc_torque(mesh, vxc, st, torque)
370 class(mesh_t), intent(in) :: mesh
371 real(real64), intent(in) :: vxc(:,:)
372 type(states_elec_t), intent(in) :: st
373 real(real64), intent(inout) :: torque(:,:)
374
375 real(real64) :: mag(3), bxc(3)
376 integer :: ip
377
378 push_sub(calc_xc_torque)
379
380 assert(st%d%ispin == spinors)
381
382 do ip = 1, mesh%np
383 mag(1) = m_two * st%rho(ip, 3)
384 mag(2) = -m_two * st%rho(ip, 4)
385 mag(3) = st%rho(ip, 1) - st%rho(ip, 2)
386 bxc(1) = -m_two * vxc(ip, 3)
387 bxc(2) = m_two * vxc(ip, 4)
388 bxc(3) = -(vxc(ip, 1) - vxc(ip, 2))
389 torque(ip, 1) = mag(2) * bxc(3) - mag(3) * bxc(2)
390 torque(ip, 2) = mag(3) * bxc(1) - mag(1) * bxc(3)
391 torque(ip, 3) = mag(1) * bxc(2) - mag(2) * bxc(1)
392 end do
393
394 pop_sub(calc_xc_torque)
395 end subroutine calc_xc_torque
396
397
398
399end module magnetic_oct_m
400
401!! Local Variables:
402!! mode: f90
403!! coding: utf-8
404!! End:
double exp(double __x) __attribute__((__nothrow__
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:122
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_curl(der, ff, op_ff, ghost_update, set_bc)
apply the curl operator to a vector of mesh functions
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
Definition: global.F90:224
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public magnetic_local_moments(mesh, st, ions, boundaries, rho, rr, lmm)
Definition: magnetic.F90:265
subroutine, public write_magnetic_moments(mesh, st, ions, boundaries, lmm_r, iunit, namespace)
Definition: magnetic.F90:211
subroutine, public calc_xc_torque(mesh, vxc, st, torque)
Definition: magnetic.F90:463
subroutine, public magnetic_moment(mesh, st, rho, mm)
Definition: magnetic.F90:183
subroutine, public magnetic_total_magnetization(mesh, st, qq, trans_mag)
Definition: magnetic.F90:334
subroutine, public write_total_xc_torque(iunit, mesh, vxc, st)
Definition: magnetic.F90:433
subroutine, public magnetic_density(mesh, std, rho, md)
Definition: magnetic.F90:156
subroutine, public magnetic_induced(namespace, gr, st, psolver, kpoints, a_ind, b_ind)
This subroutine receives as input a current, and produces as an output the vector potential that it i...
Definition: magnetic.F90:383
This module defines various routines, operating on mesh functions.
complex(real64) function, public zmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
real(real64) function, public dmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
Definition: mpi.F90:430
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:266
subroutine, public dpoisson_solve(this, namespace, pot, rho, all_nodes, kernel, reset)
Calculates the Poisson equation. Given the density returns the corresponding potential.
Definition: poisson.F90:892
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_calc_quantities(gr, st, kpoints, nlcc, kinetic_energy_density, paramagnetic_current, density_gradient, density_laplacian, gi_kinetic_energy_density, st_end)
calculated selected quantities
real(real64) function, public dsm_integrate_frommesh(mesh, sm, ff, reduce)
Definition: submesh.F90:1205
subroutine, public submesh_end(this)
Definition: submesh.F90:735
subroutine, public submesh_init(this, space, mesh, latt, center, rc)
Definition: submesh.F90:280
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:223
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
This class contains information about the boundary conditions.
Definition: boundaries.F90:157
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:168
Describes mesh distribution to nodes.
Definition: mesh.F90:186
The states_elec_t class contains all electronic wave functions.
A submesh is a type of mesh, used for the projectors in the pseudopotentials It contains points on a ...
Definition: submesh.F90:175