Branch data 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
+ + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
|
