Line data Source code
1 : !!****m* ABINIT/m_slc_potential
2 : !! NAME
3 : !! m_slc_potential
4 : !!
5 : !! FUNCTION
6 : !! This module contains the spin-lattice coupling, and the methods for
7 : !! calculating effective magnetic field (torque), force, and total_energy
8 : !!
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2001-2026 ABINIT group (TO, hexu, NH)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 : #include "abi_common.h"
23 : module m_slc_potential
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_xmpi
28 :
29 : use m_hashtable_strval, only: hash_table_t
30 : use m_abstract_potential, only : abstract_potential_t
31 : !use m_dynamic_array, only: int2d_array_type
32 : use m_mpi_scheduler, only: mb_mpi_info_t, init_mpi_info, mpi_scheduler_t
33 : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
34 : use m_multibinit_dataset, only: multibinit_dtset_type
35 : use m_spmat_ndcoo, only: ndcoo_mat_t
36 :
37 : implicit none
38 : private
39 :
40 : type, public, extends(abstract_potential_t) :: slc_potential_t
41 : integer :: nspin=0
42 : integer :: natom=0
43 : logical :: has_bilin=.False. ! bilinear coupling term, i.e. liu
44 : logical :: has_linquad=.False. ! spin first then lattice, i.e. niuv
45 : logical :: has_quadlin=.False. ! spin first then lattice, i.e. oiju
46 : logical :: has_biquad=.False. ! biquadratic coupling term, i.e. tijuv
47 :
48 : type(ndcoo_mat_t) :: liu_sc ! parameter values bilin term
49 : type(ndcoo_mat_t) :: niuv_sc ! parameter values linquad term
50 : type(ndcoo_mat_t) :: oiju_sc ! parameter values quadlin term
51 : type(ndcoo_mat_t) :: tijuv_sc ! parameter values biquad term
52 : type(ndcoo_mat_t) :: tuvij_sc ! same as tijuv but sorted differently
53 :
54 : ! magnetic moments
55 : real(dp), allocatable :: ms(:)
56 :
57 : ! precalculated things for reference spin structure
58 : real(dp), allocatable :: luref(:) ! force from liu for reference spin structure
59 : real(dp), allocatable :: ouref(:) ! force from oiju for reference spin structure
60 : type(ndcoo_mat_t) :: nuvref ! matrix in u and v from niuv for reference spin structure
61 : type(ndcoo_mat_t) :: tuvref ! matrix in u and v from tijuv for reference spin structure
62 :
63 : ! mpi
64 : type(mb_mpi_info_t) :: mpiinfo
65 : type(mpi_scheduler_t) :: mpsspin
66 : type(mpi_scheduler_t) :: mpslatt
67 :
68 : CONTAINS
69 : procedure :: initialize
70 : procedure :: finalize
71 : procedure :: set_params
72 : procedure :: set_supercell
73 : procedure :: calculate_ref
74 : procedure :: add_liu_term
75 : procedure :: add_niuv_term
76 : procedure :: add_oiju_term
77 : procedure :: add_tijuv_term
78 : procedure :: calculate
79 : end type slc_potential_t
80 :
81 : contains
82 :
83 2 : subroutine initialize(self, nspin, natom)
84 : !Arguments ------------------------------------
85 : !scalars
86 : class(slc_potential_t), intent(inout) :: self
87 : integer, intent(in) :: nspin
88 : integer, intent(in) :: natom
89 :
90 : integer :: master, my_rank, comm, nproc, ierr
91 : logical :: iam_master
92 :
93 2 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
94 2 : call self%mpsspin%initialize(ntasks=nspin, master=master, comm=comm)
95 2 : call self%mpslatt%initialize(ntasks=natom, master=master, comm=comm)
96 2 : self%label="SLCPotential"
97 2 : self%has_spin=.True.
98 2 : self%has_displacement=.True.
99 2 : self%has_strain=.False.
100 2 : self%is_null=.False.
101 2 : self%nspin=nspin
102 2 : self%natom=natom
103 :
104 2 : call xmpi_bcast(self%nspin, master, comm, ierr)
105 2 : call xmpi_bcast(self%natom, master, comm, ierr)
106 6 : ABI_MALLOC(self%ms, (self%nspin))
107 2 : end subroutine initialize
108 :
109 2 : subroutine finalize(self)
110 :
111 : class(slc_potential_t), intent(inout):: self
112 2 : ABI_SFREE(self%ms)
113 :
114 2 : if(self%has_bilin) then
115 0 : call self%liu_sc%finalize()
116 0 : ABI_SFREE(self%luref)
117 : endif
118 2 : if(self%has_quadlin) then
119 2 : call self%oiju_sc%finalize()
120 2 : ABI_SFREE(self%ouref)
121 : endif
122 2 : if(self%has_linquad) then
123 0 : call self%niuv_sc%finalize()
124 0 : call self%nuvref%finalize()
125 : endif
126 2 : if(self%has_biquad) then
127 1 : call self%tijuv_sc%finalize()
128 1 : call self%tuvij_sc%finalize()
129 1 : call self%tuvref%finalize()
130 : endif
131 :
132 2 : call self%mpsspin%finalize()
133 2 : call self%mpslatt%finalize()
134 :
135 2 : end subroutine finalize
136 :
137 : !-------------------------------------------------------------------!
138 : !set_params: which coupling terms are used
139 : !-------------------------------------------------------------------!
140 4 : subroutine set_params(self, params)
141 : class(slc_potential_t), intent(inout) :: self
142 : type(multibinit_dtset_type), intent(inout) :: params
143 :
144 : integer :: master, my_rank, comm, nproc, ierr, coupling
145 : logical :: iam_master
146 :
147 4 : coupling = params%slc_coupling
148 :
149 4 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
150 :
151 4 : if(iam_master) then
152 4 : if(coupling .ge. 1000) then
153 2 : self%has_biquad=.True.
154 2 : call xmpi_bcast(self%has_biquad, master, comm, ierr)
155 2 : coupling=coupling - 1000
156 : endif
157 :
158 4 : if(coupling .ge. 100) then
159 0 : self%has_linquad=.True.
160 0 : call xmpi_bcast(self%has_linquad, master, comm, ierr)
161 0 : coupling=coupling - 100
162 : endif
163 :
164 4 : if(coupling .ge. 10) then
165 4 : self%has_quadlin=.True.
166 4 : call xmpi_bcast(self%has_quadlin, master, comm, ierr)
167 4 : coupling=coupling - 10
168 : endif
169 :
170 4 : if(coupling .ge. 1) then
171 0 : self%has_bilin=.True.
172 0 : call xmpi_bcast(self%has_bilin, master, comm, ierr)
173 : endif
174 : endif
175 4 : end subroutine set_params
176 :
177 : !-------------------------------------------------------------------!
178 : !set_supercell: use the same supercell for all terms
179 : ! copy magnetic moments
180 : !-------------------------------------------------------------------!
181 8 : subroutine set_supercell(self, supercell)
182 : class(slc_potential_t), intent(inout) :: self
183 : type(mbsupercell_t), target, intent(inout) :: supercell
184 : integer :: master, my_rank, comm, nproc, ierr
185 : logical :: iam_master
186 :
187 4 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
188 :
189 4 : self%supercell=>supercell
190 868 : self%ms(:)=supercell%spin%ms(:)
191 :
192 4 : call xmpi_bcast(self%ms, master, comm, ierr)
193 4 : end subroutine set_supercell
194 :
195 : !-------------------------------------------------------------------!
196 : !calculate ref: calculates terms necessary for the force and energy
197 : ! of the reference structure terms
198 : !-------------------------------------------------------------------!
199 :
200 2 : subroutine calculate_ref(self)
201 : class(slc_potential_t), intent(inout) :: self
202 :
203 4 : real(dp) :: spref(1:3*self%nspin), beta
204 : real(dp), allocatable :: force(:)
205 :
206 : integer :: master, my_rank, comm, nproc
207 : logical :: iam_master
208 :
209 4 : spref(:) = reshape(self%supercell%spin%Sref, (/ 3*self%nspin/))
210 :
211 2 : beta = 0.5_dp
212 :
213 6 : ABI_MALLOC(force, (3*self%natom))
214 2 : if(self%has_bilin) then
215 0 : ABI_MALLOC(self%luref, (3*self%natom))
216 0 : self%luref=0.0d0
217 0 : force = 0.0d0
218 0 : call self%liu_sc%vec_product2d(1, spref, 2, force)
219 0 : self%luref(:) = - force(:)
220 : endif
221 2 : if(self%has_quadlin) then
222 6 : ABI_MALLOC(self%ouref, (3*self%natom))
223 6482 : force = 0.0d0
224 2 : call self%oiju_sc%vec_product(1, spref, 2, spref, 3, force)
225 6482 : self%ouref(:) = - beta*force(:)
226 : endif
227 2 : ABI_SFREE(force)
228 :
229 2 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
230 :
231 2 : if(self%has_linquad) then
232 0 : if(iam_master) then
233 0 : call self%nuvref%initialize(mshape=[self%natom*3, self%natom*3])
234 : endif
235 0 : call self%niuv_sc%mv1vec(spref, 1, self%nuvref)
236 : endif
237 :
238 2 : if(self%has_biquad) then
239 1 : if(iam_master) then
240 3 : call self%tuvref%initialize(mshape=[self%natom*3, self%natom*3])
241 : endif
242 649 : call self%tijuv_sc%mv2vec(0.5*spref, spref, 1, 2, self%tuvref)
243 : endif
244 :
245 2 : end subroutine calculate_ref
246 :
247 : !-----------------------------------------------------
248 : ! add different coupling terms to supercell potential
249 : ! TODO: test liu, niuv, tijuv
250 : !-----------------------------------------------------
251 0 : subroutine add_liu_term(self, i, u, val)
252 : class(slc_potential_t), intent(inout) :: self
253 : integer, intent(in) :: i, u
254 : real(dp), intent(in) :: val
255 :
256 : integer :: master, my_rank, comm, nproc
257 : logical :: iam_master
258 :
259 0 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
260 0 : if(iam_master) then
261 0 : call self%liu_sc%add_entry(ind=[i,u],val=val)
262 : endif
263 0 : end subroutine add_liu_term
264 :
265 :
266 0 : subroutine add_niuv_term(self, i, u, v, val)
267 : class(slc_potential_t), intent(inout) :: self
268 : integer, intent(in) :: i, u, v
269 : real(dp), intent(in) :: val
270 :
271 : integer :: master, my_rank, comm, nproc
272 : logical :: iam_master
273 :
274 0 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
275 0 : if(iam_master) then
276 0 : call self%niuv_sc%add_entry(ind=[i,u,v],val=val)
277 : endif
278 0 : end subroutine add_niuv_term
279 :
280 :
281 :
282 1041984 : subroutine add_oiju_term(self, i,j,u, val)
283 : class(slc_potential_t), intent(inout) :: self
284 : integer, intent(in) :: i, j, u
285 : real(dp), intent(in) :: val
286 :
287 : integer :: master, my_rank, comm, nproc
288 : logical :: iam_master
289 :
290 1041984 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
291 1041984 : if(iam_master) then
292 4167936 : call self%oiju_sc%add_entry(ind=[i,j,u],val=val)
293 : endif
294 1041984 : end subroutine add_oiju_term
295 :
296 69984 : subroutine add_tijuv_term(self, i,j,u,v, val)
297 : class(slc_potential_t), intent(inout) :: self
298 : integer, intent(in) :: i, j, u, v
299 : real(dp), intent(in) :: val
300 :
301 : integer :: master, my_rank, comm, nproc
302 : logical :: iam_master
303 :
304 69984 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
305 69984 : if(iam_master) then
306 349920 : call self%tijuv_sc%add_entry(ind=[i,j,u,v],val=val)
307 349920 : call self%tuvij_sc%add_entry(ind=[u,v,i,j],val=val)
308 : endif
309 69984 : end subroutine add_tijuv_term
310 :
311 :
312 : !-----------------------------------------------------------------
313 : ! Calculate forces, magnetic fields and energy for coupling terms
314 : ! TODO: precalculate terms containing Sref?
315 : !-----------------------------------------------------------------
316 :
317 8008 : subroutine calculate(self, displacement, strain, spin, lwf, &
318 8008 : force, stress, bfield, lwf_force, energy, energy_table)
319 : class(slc_potential_t), intent(inout) :: self
320 : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
321 : real(dp), optional, intent(inout) :: force(:,:), stress(:,:), bfield(:,:), lwf_force(:), energy
322 : type(hash_table_t),optional, intent(inout) :: energy_table
323 :
324 : integer :: ii
325 : character(len=80) :: label
326 : real(dp) :: eslc, beta, eterm
327 16016 : real(dp) :: disp(1:3*self%natom), sp(1:3*self%nspin), spref(1:3*self%nspin)
328 16016 : real(dp) :: f1(1:3*self%natom), b1(1:3*self%nspin)
329 16016 : real(dp) :: btmp(3, self%nspin), bslc(1:3*self%nspin), fslc(1:3*self%natom)
330 :
331 : ABI_UNUSED(lwf)
332 : ABI_UNUSED(strain)
333 : ABI_UNUSED(lwf_force)
334 : ABI_UNUSED(stress)
335 :
336 16016 : sp(:) = reshape(spin, (/ 3*self%nspin /))
337 16016 : spref(:) = reshape(self%supercell%spin%Sref, (/ 3*self%nspin/))
338 16016 : disp(:) = reshape(displacement, (/ 3*self%natom /))
339 :
340 8008 : beta = 0.5_dp
341 :
342 : ! Magnetic field
343 8008 : if(present(bfield)) then
344 2598596 : bslc(:) = 0.0d0
345 4004 : if(self%has_bilin) then
346 0 : b1(:) = 0.0d0
347 0 : call self%liu_sc%vec_product2d(2, disp, 1, b1)
348 0 : bslc(:) = bslc(:) + b1(:)
349 : endif
350 4004 : if(self%has_linquad) then
351 0 : b1(:) = 0.0d0
352 0 : call self%niuv_sc%vec_product(2, disp, 3, disp, 1, b1)
353 0 : bslc(:) = bslc(:) + beta*b1(:)
354 : endif
355 4004 : if(self%has_quadlin) then
356 2598596 : b1(:) = 0.0d0
357 4004 : call self%oiju_sc%vec_product(1, sp, 3, disp, 2, b1)
358 2598596 : bslc(:) = bslc(:) + 2.0d0*beta*b1(:)
359 : endif
360 4004 : if(self%has_biquad) then
361 1299298 : b1(:) = 0.0d0
362 2002 : call self%tuvij_sc%vec_product4d(disp, disp, 3, sp, 4, b1)
363 1299298 : bslc(:) = bslc(:) + beta*b1(:)
364 : endif
365 12012 : btmp = reshape(bslc, (/ 3, self%nspin /))
366 868868 : do ii = 1, self%nspin
367 3463460 : btmp(:,ii) = btmp(:,ii)/self%ms(ii)
368 : enddo
369 3463460 : bfield(:,:) = bfield(:,:) + btmp(:,:)
370 :
371 : ! TESTING: write magnetic fields to a file
372 : !write(201,*) 'Magnetic fields are'
373 : !do ii = 1, self%nspin
374 : ! write(201,*) ii, btmp(:,ii)
375 : !enddo
376 : endif
377 :
378 : ! Force and energy
379 8008 : if(present(force) .or. present(energy) .or. present(energy_table)) then
380 25953928 : fslc(:) = 0.0d0
381 8008 : eslc = 0.0d0
382 8008 : if(self%has_bilin) then
383 : eterm = 0.0d0
384 0 : f1(:) = 0.0d0
385 0 : call self%liu_sc%vec_product2d(1, sp, 2, f1)
386 0 : fslc(:) = fslc(:) + f1(:)
387 0 : eterm = - dot_product(f1, disp)
388 : ! add contributions from reference spin structure
389 0 : fslc(:) = fslc(:) + self%luref(:)
390 0 : eterm = eterm - dot_product(self%luref, disp)
391 0 : if(present(energy_table)) then
392 0 : label=trim(self%label)//'_Liu'
393 0 : call energy_table%put(label, eterm)
394 : endif
395 0 : eslc = eslc + eterm
396 : endif
397 8008 : if(self%has_linquad) then
398 0 : f1(:) = 0.0d0
399 : eterm = 0.0d0
400 0 : call self%niuv_sc%vec_product(1, sp, 2, disp, 3, f1)
401 0 : fslc(:) = fslc(:) + 2.0d0*beta*f1(:)
402 0 : eterm = - beta*dot_product(f1, disp)
403 : ! add contributions from reference spin structure
404 0 : f1(:) = 0.0d0
405 0 : call self%nuvref%vec_product2d(1, disp, 2, f1)
406 0 : fslc(:)= fslc+2.0*beta*f1(:)
407 0 : eterm = eterm + beta*dot_product(f1, disp)
408 0 : if(present(energy_table)) then
409 0 : label=trim(self%label)//'_Niuv'
410 0 : call energy_table%put(label, eterm)
411 : endif
412 0 : eslc = eslc + eterm
413 : endif
414 8008 : if(self%has_quadlin) then
415 25953928 : f1(:) = 0.0d0
416 : eterm = 0.0d0
417 8008 : call self%oiju_sc%vec_product(1, sp, 2, sp, 3, f1)
418 25953928 : fslc(:) = fslc(:) + beta*f1(:)
419 25953928 : eterm = - beta*dot_product(f1, disp)
420 : ! add contributions from reference spin structure
421 25953928 : fslc(:) = fslc(:) + self%ouref(:)
422 25953928 : eterm = eterm - dot_product(self%ouref, disp)
423 8008 : if(present(energy_table)) then
424 8008 : label=trim(self%label)//'_Oiju'
425 8008 : call energy_table%put(label, eterm)
426 : endif
427 8008 : eslc = eslc + eterm
428 : endif
429 8008 : if(self%has_biquad) then
430 12976964 : f1(:) = 0.0d0
431 : eterm = 0.0d0
432 4004 : call self%tijuv_sc%vec_product4d(sp, sp, 3, disp, 4, f1)
433 12976964 : fslc(:) = fslc(:) + beta*f1(:)
434 12976964 : eterm = - 0.5_dp*beta*dot_product(f1, disp)
435 : ! add contributions from reference spin structure
436 12976964 : f1(:) = 0.0d0
437 4004 : call self%tuvref%vec_product2d(1, disp, 2, f1)
438 12976964 : fslc(:)= fslc+2.0*beta*f1(:)
439 12976964 : eterm = eterm + beta*dot_product(f1, disp)
440 4004 : if(present(energy_table)) then
441 4004 : label=trim(self%label)//'_Tijuv'
442 4004 : call energy_table%put(label, eterm)
443 : endif
444 4004 : eslc = eslc + eterm
445 : endif
446 : endif !energy or force
447 :
448 8008 : if(present(force)) then
449 17309292 : force(:,:) = force(:,:) + reshape(fslc, (/3, self%natom /))
450 : !TESTING write forces to file
451 : !write(200,*) 'Forces are'
452 : !do ii = 1, self%natom
453 : ! write(200,*) ii, force(:,ii)
454 : !enddo
455 : endif
456 :
457 8008 : if(present(energy)) energy = energy + eslc
458 :
459 :
460 8008 : ABI_UNUSED_A(strain)
461 8008 : ABI_UNUSED_A(lwf)
462 8008 : ABI_UNUSED_A(stress)
463 8008 : ABI_UNUSED_A(lwf_force)
464 :
465 8008 : end subroutine calculate
466 :
467 : !!***
468 6 : end module m_slc_potential
|
