TABLE OF CONTENTS


ABINIT/m_slc_potential [ Modules ]

[ Top ] [ Modules ]

NAME

 m_slc_potential

FUNCTION

 This module contains the spin-lattice coupling, and the methods for
 calculating effective magnetic field (torque), force, and total_energy

COPYRIGHT

 Copyright (C) 2001-2022 ABINIT group (TO, hexu, NH)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

SOURCE

 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   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     call init_mpi_info(master, iam_master, my_rank, comm, nproc)
 94     call self%mpsspin%initialize(ntasks=nspin, master=master, comm=comm)
 95     call self%mpslatt%initialize(ntasks=natom, master=master, comm=comm)
 96     self%label="SLCPotential"
 97     self%has_spin=.True.
 98     self%has_displacement=.True.
 99     self%has_strain=.False.
100     self%is_null=.False.
101     self%nspin=nspin
102     self%natom=natom
103     
104     call xmpi_bcast(self%nspin, master, comm, ierr)
105     call xmpi_bcast(self%natom, master, comm, ierr)
106     ABI_MALLOC(self%ms, (self%nspin))
107   end subroutine initialize
108 
109   subroutine finalize(self)
110 
111     class(slc_potential_t), intent(inout):: self
112     ABI_SFREE(self%ms)
113 
114     if(self%has_bilin) then 
115       call self%liu_sc%finalize()
116       ABI_SFREE(self%luref)
117     endif
118     if(self%has_quadlin) then
119       call self%oiju_sc%finalize()
120       ABI_SFREE(self%ouref)
121     endif
122     if(self%has_linquad) then
123       call self%niuv_sc%finalize()
124       call self%nuvref%finalize()
125     endif 
126     if(self%has_biquad) then
127       call self%tijuv_sc%finalize()
128       call self%tuvij_sc%finalize()
129       call self%tuvref%finalize()
130     endif
131 
132     call self%mpsspin%finalize()
133     call self%mpslatt%finalize()
134 
135   end subroutine finalize
136 
137   !-------------------------------------------------------------------!
138   !set_params: which coupling terms are used
139   !-------------------------------------------------------------------!
140   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     coupling = params%slc_coupling
148 
149     call init_mpi_info(master, iam_master, my_rank, comm, nproc) 
150 
151     if(iam_master) then
152       if(coupling .ge. 1000) then
153         self%has_biquad=.True.
154         call xmpi_bcast(self%has_biquad, master, comm, ierr)
155         coupling=coupling - 1000
156       endif
157 
158       if(coupling .ge. 100) then
159         self%has_linquad=.True.
160         call xmpi_bcast(self%has_linquad, master, comm, ierr)
161         coupling=coupling - 100
162       endif
163 
164       if(coupling .ge. 10) then
165         self%has_quadlin=.True.
166         call xmpi_bcast(self%has_quadlin, master, comm, ierr)
167         coupling=coupling - 10
168       endif
169 
170       if(coupling .ge. 1) then   
171         self%has_bilin=.True.
172         call xmpi_bcast(self%has_bilin, master, comm, ierr)
173       endif
174     endif
175   end subroutine set_params
176 
177   !-------------------------------------------------------------------!
178   !set_supercell: use the same supercell for all terms
179   !               copy magnetic moments
180   !-------------------------------------------------------------------!
181   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     call init_mpi_info(master, iam_master, my_rank, comm, nproc) 
188 
189     self%supercell=>supercell
190     self%ms(:)=supercell%spin%ms(:)
191 
192     call xmpi_bcast(self%ms, master, comm, ierr)
193   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   subroutine calculate_ref(self)
201     class(slc_potential_t), intent(inout) :: self
202 
203     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     spref(:) = reshape(self%supercell%spin%Sref, (/ 3*self%nspin/))
210 
211     beta = 0.5_dp
212 
213     ABI_MALLOC(force, (3*self%natom))
214     if(self%has_bilin) then 
215       ABI_MALLOC(self%luref, (3*self%natom))
216       self%luref=0.0d0
217       force = 0.0d0
218       call self%liu_sc%vec_product2d(1, spref, 2, force)
219       self%luref(:) = - force(:)
220     endif
221     if(self%has_quadlin) then
222       ABI_MALLOC(self%ouref, (3*self%natom))     
223       force = 0.0d0
224       call self%oiju_sc%vec_product(1, spref, 2, spref, 3, force)
225       self%ouref(:) = - beta*force(:)
226     endif
227     ABI_SFREE(force)
228 
229     call init_mpi_info(master, iam_master, my_rank, comm, nproc) 
230 
231     if(self%has_linquad) then
232       if(iam_master) then
233         call self%nuvref%initialize(mshape=[self%natom*3, self%natom*3])
234       endif
235       call self%niuv_sc%mv1vec(spref, 1, self%nuvref)
236     endif
237 
238     if(self%has_biquad) then
239       if(iam_master) then
240          call self%tuvref%initialize(mshape=[self%natom*3, self%natom*3])
241       endif
242       call self%tijuv_sc%mv2vec(0.5*spref, spref, 1, 2, self%tuvref)
243     endif
244 
245   end subroutine calculate_ref
246 
247   !-----------------------------------------------------
248   ! add different coupling terms to supercell potential
249   ! TODO: test liu, niuv, tijuv
250   !-----------------------------------------------------
251   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     call init_mpi_info(master, iam_master, my_rank, comm, nproc) 
260     if(iam_master) then
261        call self%liu_sc%add_entry(ind=[i,u],val=val)
262     endif
263   end subroutine add_liu_term
264 
265 
266   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     call init_mpi_info(master, iam_master, my_rank, comm, nproc) 
275     if(iam_master) then
276        call self%niuv_sc%add_entry(ind=[i,u,v],val=val)
277     endif
278   end subroutine add_niuv_term
279 
280 
281 
282   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     call init_mpi_info(master, iam_master, my_rank, comm, nproc) 
291     if(iam_master) then
292        call self%oiju_sc%add_entry(ind=[i,j,u],val=val)
293     endif
294   end subroutine add_oiju_term
295 
296   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     call init_mpi_info(master, iam_master, my_rank, comm, nproc) 
305     if(iam_master) then
306        call self%tijuv_sc%add_entry(ind=[i,j,u,v],val=val)
307        call self%tuvij_sc%add_entry(ind=[u,v,i,j],val=val)
308     endif
309   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   subroutine calculate(self, displacement, strain, spin, lwf, &
318        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     real(dp) :: disp(1:3*self%natom), sp(1:3*self%nspin), spref(1:3*self%nspin) 
328     real(dp) :: f1(1:3*self%natom), b1(1:3*self%nspin)
329     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     sp(:) = reshape(spin, (/ 3*self%nspin /))
337     spref(:) = reshape(self%supercell%spin%Sref, (/ 3*self%nspin/))
338     disp(:) = reshape(displacement, (/ 3*self%natom /))
339 
340     beta = 0.5_dp
341 
342     ! Magnetic field
343     if(present(bfield)) then
344       bslc(:) = 0.0d0
345       if(self%has_bilin) then
346         b1(:) = 0.0d0
347         call self%liu_sc%vec_product2d(2, disp, 1, b1)
348         bslc(:) = bslc(:) + b1(:)
349       endif      
350       if(self%has_linquad) then
351         b1(:) = 0.0d0
352         call self%niuv_sc%vec_product(2, disp, 3, disp, 1, b1)
353         bslc(:) = bslc(:) + beta*b1(:)
354       endif      
355       if(self%has_quadlin) then
356         b1(:) = 0.0d0
357         call self%oiju_sc%vec_product(1, sp, 3, disp, 2, b1)
358         bslc(:) = bslc(:) + 2.0d0*beta*b1(:)
359       endif
360       if(self%has_biquad) then
361         b1(:) = 0.0d0
362         call self%tuvij_sc%vec_product4d(disp, disp, 3, sp, 4, b1)
363         bslc(:) = bslc(:) + beta*b1(:)
364       endif
365       btmp = reshape(bslc, (/ 3, self%nspin /))
366       do ii = 1, self%nspin
367         btmp(:,ii) = btmp(:,ii)/self%ms(ii)
368       enddo
369       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     if(present(force) .or. present(energy) .or. present(energy_table)) then
380       fslc(:) = 0.0d0
381       eslc = 0.0d0
382       if(self%has_bilin) then
383         eterm = 0.0d0
384         f1(:) = 0.0d0
385         call self%liu_sc%vec_product2d(1, sp, 2, f1)
386         fslc(:) = fslc(:) + f1(:)
387         eterm =  - dot_product(f1, disp)
388         ! add contributions from reference spin structure
389         fslc(:) = fslc(:) + self%luref(:)
390         eterm = eterm - dot_product(self%luref, disp)
391         if(present(energy_table)) then
392           label=trim(self%label)//'_Liu'
393           call energy_table%put(label, eterm)
394         endif
395         eslc = eslc + eterm
396       endif      
397       if(self%has_linquad) then
398         f1(:) = 0.0d0
399         eterm = 0.0d0
400         call self%niuv_sc%vec_product(1, sp, 2, disp, 3, f1)
401         fslc(:) = fslc(:) + 2.0d0*beta*f1(:)
402         eterm =  - beta*dot_product(f1, disp)
403         ! add contributions from reference spin structure
404         f1(:) = 0.0d0
405         call self%nuvref%vec_product2d(1, disp, 2, f1)
406         fslc(:)= fslc+2.0*beta*f1(:)
407         eterm = eterm + beta*dot_product(f1, disp)
408         if(present(energy_table)) then
409           label=trim(self%label)//'_Niuv'
410           call energy_table%put(label, eterm)
411         endif
412         eslc = eslc + eterm
413       endif      
414       if(self%has_quadlin) then
415         f1(:) = 0.0d0
416         eterm = 0.0d0
417         call self%oiju_sc%vec_product(1, sp, 2, sp, 3, f1)
418         fslc(:) = fslc(:) + beta*f1(:)
419         eterm = - beta*dot_product(f1, disp)
420         ! add contributions from reference spin structure
421         fslc(:) = fslc(:) + self%ouref(:)
422         eterm = eterm - dot_product(self%ouref, disp)
423         if(present(energy_table)) then
424           label=trim(self%label)//'_Oiju'
425           call energy_table%put(label, eterm)
426         endif
427         eslc = eslc + eterm
428       endif
429       if(self%has_biquad) then
430         f1(:) = 0.0d0
431         eterm = 0.0d0
432         call self%tijuv_sc%vec_product4d(sp, sp, 3, disp, 4, f1)
433         fslc(:) = fslc(:) + beta*f1(:)
434         eterm = - 0.5_dp*beta*dot_product(f1, disp)
435         ! add contributions from reference spin structure
436         f1(:) = 0.0d0
437         call self%tuvref%vec_product2d(1, disp, 2, f1)
438         fslc(:)= fslc+2.0*beta*f1(:)
439         eterm = eterm + beta*dot_product(f1, disp)
440         if(present(energy_table)) then
441           label=trim(self%label)//'_Tijuv'
442           call energy_table%put(label, eterm)
443         endif
444         eslc = eslc + eterm
445       endif
446     endif !energy or force
447 
448     if(present(force)) then
449       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     if(present(energy)) energy =  energy + eslc
458 
459 
460     ABI_UNUSED_A(strain)
461     ABI_UNUSED_A(lwf)
462     ABI_UNUSED_A(stress)
463     ABI_UNUSED_A(lwf_force)
464     
465   end subroutine calculate