TABLE OF CONTENTS
ABINIT/m_slc_potential [ 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