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 kick_oct_m 22 use iso_c_binding 23 use geometry_oct_m 24 use global_oct_m 25 use ion_dynamics_oct_m 26 use kpoints_oct_m 27 use loct_math_oct_m 28 use math_oct_m 29 use mesh_oct_m 30 use messages_oct_m 31 use namespace_oct_m 32 use parser_oct_m 33 use pcm_eom_oct_m 34 use pcm_oct_m 35 use poisson_oct_m 36 use profiling_oct_m 37 use species_oct_m 38 use simul_box_oct_m 39 use states_elec_oct_m 40 use states_elec_dim_oct_m 41 use symm_op_oct_m 42 use symmetries_oct_m 43 use unit_oct_m 44 use unit_system_oct_m 45 use write_iter_oct_m 46 47 implicit none 48 49 private 50 public :: & 51 kick_t, & 52 kick_init, & 53 kick_copy, & 54 kick_end, & 55 kick_read, & 56 kick_write, & 57 kick_apply, & 58 kick_function_get, & 59 kick_get_type 60 61 62 integer, public, parameter :: & 63 KICK_FUNCTION_DIPOLE = 0, & 64 KICK_FUNCTION_MULTIPOLE = 1, & 65 KICK_FUNCTION_USER_DEFINED = 2 66 67 integer, public, parameter :: & 68 KICK_DENSITY_MODE = 0, & 69 KICK_SPIN_MODE = 1, & 70 KICK_SPIN_DENSITY_MODE = 2, & 71 KICK_MAGNON_MODE = 3 72 73 integer, public, parameter :: & 74 QKICKMODE_NONE = 0, & 75 QKICKMODE_EXP = 1, & 76 QKICKMODE_COS = 2, & 77 QKICKMODE_SIN = 3, & 78 QKICKMODE_BESSEL = 4 79 80 81 type kick_t 82 ! Components are public by default 83 84 !> The time which the kick is applied (normally, this is zero) 85 FLOAT :: time 86 !> The strength, and strength "mode". 87 integer, private :: delta_strength_mode 88 FLOAT :: delta_strength 89 !> In case we use a normal dipole kick: 90 FLOAT :: pol(MAX_DIM, MAX_DIM) 91 integer :: pol_dir 92 integer :: pol_equiv_axes 93 FLOAT :: wprime(MAX_DIM) 94 FLOAT :: easy_axis(MAX_DIM) 95 !> In case we have a general multipolar kick, 96 !! the form of this "kick" will be (atomic units): 97 !! \f[ 98 !! V(\vec{r}) = sum_{i=1}^{n\_multipoles} 99 !! weight(i) * (e^2 / a_0^(l+1)) * r^l(i) * Y_{l(i),m(i)} (\vec{r}) 100 !! \f] 101 !! which has units of energy; if we include the time-dependence (delta function): 102 !! \f[ 103 !! V(\vec{r}) = sum_{i=1}^{n\_multipoles} 104 !! weight(i) * (\hbar / a_0^l) * r^l(i) * Y_{l(i),m(i)} (\vec{r}) * \delta(t) 105 !! \f] 106 integer :: n_multipoles 107 integer, pointer :: l(:), m(:) 108 FLOAT, pointer :: weight(:) 109 integer :: nqmult(1:MAX_DIM) 110 integer :: nqvec 111 FLOAT, allocatable :: qvector(:,:) 112 FLOAT :: trans_vec(MAX_DIM,2) 113 FLOAT :: qlength 114 integer :: qkick_mode 115 integer :: qbessel_l, qbessel_m 116 !> In case we use a general function 117 integer :: function_mode 118 character(len=200), private:: user_defined_function 119 end type kick_t 120 121contains 122 123 ! --------------------------------------------------------- 124 subroutine kick_init(kick, namespace, sb, nspin) 125 type(kick_t), intent(out) :: kick 126 type(namespace_t), intent(in) :: namespace 127 type(simul_box_t), intent(in) :: sb 128 integer, intent(in) :: nspin 129 130 type(block_t) :: blk 131 integer :: n_rows, irow, idir, iop, iq, iqx, iqy, iqz 132 FLOAT :: norm, dot 133 FLOAT :: qtemp(1:MAX_DIM) 134 integer :: dim, periodic_dim 135 136 PUSH_SUB(kick_init) 137 138 dim = sb%dim 139 periodic_dim = sb%periodic_dim 140 141 !%Variable TDDeltaKickTime 142 !%Type float 143 !%Default 0.0 144 !%Section Time-Dependent::Response 145 !%Description 146 !% The delta-perturbation that can be applied by making use of the <tt>TDDeltaStrength</tt> variable, 147 !% can be applied at the time given by this variable. Usually, this time is zero, since one wants 148 !% to apply the delta pertubation or "kick" at a system at equilibrium, and no other time-dependent 149 !% external potential is used. However, one may want to apply a kick on top of a laser field, 150 !% for example. 151 !%End 152 call parse_variable(namespace, 'TDDeltaKickTime', M_ZERO, kick%time, units_inp%time) 153 if(kick%time > M_ZERO) then 154 call messages_experimental('TDDeltaKickTime > 0') 155 end if 156 157 !%Variable TDDeltaStrength 158 !%Type float 159 !%Default 0 160 !%Section Time-Dependent::Response 161 !%Description 162 !% When no laser is applied, a delta (in time) perturbation with 163 !% strength <tt>TDDeltaStrength</tt> can be applied. This is used to 164 !% calculate, <i>e.g.</i>, the linear optical spectra. If the ions are 165 !% allowed to move, the kick will affect them also. 166 !% The electric field is <math>-(\hbar k / e) \delta(t)</math> for a dipole with 167 !% zero wavevector, where <i>k</i> = <tt>TDDeltaStrength</tt>, which causes 168 !% the wavefunctions instantaneously to acquire a phase <math>e^{ikx}</math>. 169 !% The unit is inverse length. 170 !%End 171 call parse_variable(namespace, 'TDDeltaStrength', M_ZERO, kick%delta_strength, units_inp%length**(-1)) 172 173 nullify(kick%l) 174 nullify(kick%m) 175 nullify(kick%weight) 176 kick%function_mode = KICK_FUNCTION_DIPOLE 177 178 if(abs(kick%delta_strength) <= M_EPSILON) then 179 kick%delta_strength_mode = 0 180 kick%pol_equiv_axes = 0 181 kick%pol(1:3, 1) = (/ M_ONE, M_ZERO, M_ZERO /) 182 kick%pol(1:3, 2) = (/ M_ZERO, M_ONE, M_ZERO /) 183 kick%pol(1:3, 3) = (/ M_ZERO, M_ZERO, M_ONE /) 184 kick%pol_dir = 0 185 kick%wprime = M_ZERO 186 kick%n_multipoles = 0 187 kick%qkick_mode = QKICKMODE_NONE 188 kick%easy_axis(1:MAX_DIM) = M_ZERO 189 POP_SUB(kick_init) 190 return 191 end if 192 193 !%Variable TDDeltaStrengthMode 194 !%Type integer 195 !%Default kick_density 196 !%Section Time-Dependent::Response 197 !%Description 198 !% When calculating the density response via real-time propagation, 199 !% one needs to perform an initial kick on the KS system, at 200 !% time zero. Depending on what kind of response property one wants to obtain, 201 !% this kick may be done in several modes. For use to calculate triplet excitations, 202 !% see MJT Oliveira, A Castro, MAL Marques, and A Rubio, <i>J. Nanoscience and Nanotechnology</i> <b>8</b>, 3392 (2008). 203 !%Option kick_density 0 204 !% The total density of the system is perturbed. This mode is appropriate for 205 !% electric dipole response, as for optical absorption. 206 !%Option kick_spin 1 207 !% The individual spin densities are perturbed oppositely. Note that this mode 208 !% is only possible if the run is done in spin-polarized mode, or with spinors. 209 !% This mode is appropriate for the paramagnetic dipole response, which can couple 210 !% to triplet excitations. 211 !%Option kick_spin_and_density 2 212 !% A combination of the two above. Note that this mode 213 !% is only possible if the run is done in spin-polarized mode, or with spinors. 214 !% This mode is intended for use with symmetries to obtain both of the responses 215 !% at once, at described in the reference above. 216 !%Option kick_magnon 3 217 !% Rotates the magnetization. Only works for spinors. 218 !% Can be used in a supercell or my making use of the generalized Bloch theorem. 219 !% In the later case (see <tt>SpiralBoundaryConditions</tt>) spin-orbit coupling cannot be used. 220 !%End 221 call parse_variable(namespace, 'TDDeltaStrengthMode', KICK_DENSITY_MODE, kick%delta_strength_mode) 222 select case (kick%delta_strength_mode) 223 case (KICK_DENSITY_MODE) 224 case (KICK_SPIN_MODE, KICK_SPIN_DENSITY_MODE) 225 case (KICK_MAGNON_MODE) 226 if(nspin /= SPINORS) call messages_input_error(namespace, 'TDDeltaStrengthMode') 227 case default 228 call messages_input_error(namespace, 'TDDeltaStrengthMode') 229 end select 230 231 if(parse_is_defined(namespace, 'TDDeltaUserDefined')) then 232 233 kick%function_mode = KICK_FUNCTION_USER_DEFINED 234 kick%n_multipoles = 0 235 236 !%Variable TDDeltaUserDefined 237 !%Type string 238 !%Section Time-Dependent::Response 239 !%Description 240 !% By default, the kick function will be a dipole. This will change if (1) the variable 241 !% <tt>TDDeltaUserDefined</tt> is present in the inp file, or (2) if the block <tt>TDKickFunction</tt> 242 !% is present in the <tt>inp</tt> file. If both are present in the <tt>inp</tt> file, the <tt>TDKickFunction</tt> 243 !% block will be ignored. The value of <tt>TDDeltaUserDefined</tt> should be a string describing 244 !% the function that is going to be used as delta perturbation. 245 !%End 246 call parse_variable(namespace, 'TDDeltaUserDefined', '0', kick%user_defined_function) 247 248 !%Variable TDKickFunction 249 !%Type block 250 !%Section Time-Dependent::Response 251 !%Description 252 !% If the block <tt>TDKickFunction</tt> is present in the input file, and the variable 253 !% <tt>TDDeltaUserDefined</tt> is not present in the input file, the kick function to 254 !% be applied at time zero of the time-propagation will not be a "dipole" function 255 !% (<i>i.e.</i> <math>\phi \rightarrow e^{ikx} \phi</math>, but a general multipole in the form <math>r^l Y_{lm}(r)</math>. 256 !% 257 !% Each line has three columns: integers <i>l</i> and <i>m</i> that defines the 258 !% multipole, and a weight. Any number of lines may be given, and the kick will be the sum of those 259 !% multipoles with the given weights. 260 !% 261 !% This feature allows calculation of quadrupole, octupole, etc., response functions. 262 !%End 263 else if(parse_block(namespace, 'TDKickFunction', blk) == 0) then 264 265 kick%function_mode = KICK_FUNCTION_MULTIPOLE 266 n_rows = parse_block_n(blk) 267 kick%n_multipoles = n_rows 268 SAFE_ALLOCATE( kick%l(1:n_rows)) 269 SAFE_ALLOCATE( kick%m(1:n_rows)) 270 SAFE_ALLOCATE(kick%weight(1:n_rows)) 271 do irow = 1, n_rows 272 call parse_block_integer(blk, irow - 1, 0, kick%l(irow)) 273 call parse_block_integer(blk, irow - 1, 1, kick%m(irow)) 274 call parse_block_float(blk, irow - 1, 2, kick%weight(irow)) 275 if( (kick%l(irow) < 0) .or. (abs(kick%m(irow)) > abs(kick%l(irow))) ) then 276 call messages_input_error(namespace, 'TDkickFunction') 277 end if 278 end do 279 280 else 281 282 kick%function_mode = KICK_FUNCTION_DIPOLE 283 kick%n_multipoles = 0 284 285 ! Find out how many equivalent axes we have... 286 !%Variable TDPolarizationEquivAxes 287 !%Type integer 288 !%Default 0 289 !%Section Time-Dependent::Response::Dipole 290 !%Description 291 !% Defines how many of the <tt>TDPolarization</tt> axes are equivalent. This information is stored in a file and then 292 !% used by <tt>oct-propagation_spectrum</tt> to rebuild the full polarizability tensor from just the 293 !% first <tt>TDPolarizationEquivAxes</tt> directions. This variable is also used by <tt>CalculationMode = vdw</tt>. 294 !%End 295 call parse_variable(namespace, 'TDPolarizationEquivAxes', 0, kick%pol_equiv_axes) 296 297 !%Variable TDPolarizationDirection 298 !%Type integer 299 !%Section Time-Dependent::Response::Dipole 300 !%Description 301 !% When a delta potential is included in a time-dependent run, this 302 !% variable defines in which direction the field will be applied 303 !% by selecting one of the lines of <tt>TDPolarization</tt>. In a 304 !% typical run (without using symmetry), the <tt>TDPolarization</tt> block 305 !% would contain the three Cartesian unit vectors (the default 306 !% value), and one would make 3 runs varying 307 !% <tt>TDPolarization</tt> from 1 to 3. 308 !% If one is using symmetry, <tt>TDPolarization</tt> should run only from 1 309 !% to <tt>TDPolarizationEquivAxes</tt>. 310 !%End 311 312 call parse_variable(namespace, 'TDPolarizationDirection', 0, kick%pol_dir) 313 314 if(kick%delta_strength_mode /= KICK_MAGNON_MODE) then 315 if(kick%pol_dir < 1 .or. kick%pol_dir > dim) call messages_input_error(namespace, 'TDPolarizationDirection') 316 end if 317 318 !%Variable TDPolarization 319 !%Type block 320 !%Section Time-Dependent::Response::Dipole 321 !%Description 322 !% The (real) polarization of the delta electric field. Normally 323 !% one needs three perpendicular polarization directions to calculate a 324 !% spectrum (unless symmetry is used). 325 !% The format of the block is: 326 !% 327 !% <tt>%TDPolarization 328 !% <br> pol1x | pol1y | pol1z 329 !% <br> pol2x | pol2y | pol2z 330 !% <br> pol3x | pol3y | pol3z 331 !% <br>%</tt> 332 !% 333 !% <tt>Octopus</tt> uses both this block and the variable 334 !% <tt>TDPolarizationDirection</tt> to determine the polarization 335 !% vector for the run. For example, if 336 !% <tt>TDPolarizationDirection=2</tt> the polarization <tt>(pol2x, 337 !% pol2y, pol2z)</tt> would be used. 338 !% These directions may not be in periodic directions. 339 !% 340 !% The default value for <tt>TDPolarization</tt> is the three 341 !% Cartesian unit vectors (1,0,0), (0,1,0), and (0,0,1). 342 !% 343 !% Note that the directions do not necessarily need to be perpendicular 344 !% when symmetries are used. 345 !% 346 !% WARNING: If you want to obtain the cross-section tensor, the 347 !% <tt>TDPolarization</tt> block must be exactly the same for the run in 348 !% each direction. The direction must be selected by the 349 !% <tt>TDPolarizationDirection</tt> variable. 350 !% 351 !%End 352 353 kick%pol(:, :) = M_ZERO 354 if(parse_block(namespace, 'TDPolarization', blk)==0) then 355 n_rows = parse_block_n(blk) 356 357 if(n_rows < dim) call messages_input_error(namespace, 'TDPolarization') 358 359 do irow = 1, n_rows 360 do idir = 1, 3 361 call parse_block_float(blk, irow - 1, idir - 1, kick%pol(idir, irow)) 362 end do 363 end do 364 if(n_rows < 3) kick%pol(1:3, 3) = (/ M_ZERO, M_ZERO, M_ONE /) 365 if(n_rows < 2) kick%pol(1:3, 2) = (/ M_ZERO, M_ONE, M_ZERO /) 366 call parse_block_end(blk) 367 else 368 ! FIXME: Here the symmetry of the system should be analyzed, and the polarization 369 ! basis built accordingly. 370 kick%pol(1:3, 1) = (/ M_ONE, M_ZERO, M_ZERO /) 371 kick%pol(1:3, 2) = (/ M_ZERO, M_ONE, M_ZERO /) 372 kick%pol(1:3, 3) = (/ M_ZERO, M_ZERO, M_ONE /) 373 end if 374 375 ! Normalize 376 do idir = 1, 3 377 kick%pol(1:3, idir) = kick%pol(1:3, idir) / sqrt(sum(kick%pol(1:3, idir)**2)) 378 end do 379 380 if(kick%delta_strength_mode /= KICK_MAGNON_MODE) then 381 if(any(abs(kick%pol(1:periodic_dim, :)) > M_EPSILON)) then 382 message(1) = "Kick cannot be applied in a periodic direction. Use GaugeVectorField instead." 383 call messages_fatal(1, namespace=namespace) 384 end if 385 end if 386 387 !%Variable TDPolarizationWprime 388 !%Type block 389 !%Section Time-Dependent::Response::Dipole 390 !%Description 391 !% This block is needed only when 392 !% <tt>TDPolarizationEquivAxes</tt> is set to 3. In such a case, 393 !% the three directions (<i>pol1</i>, <i>pol2</i>, and <i>pol3</i>) defined in 394 !% the <tt>TDPolarization</tt> block should be related by symmetry 395 !% operations. If <i>A</i> is the symmetry operation that takes you 396 !% from <i>pol1</i> to <i>pol2</i>, then <tt>TDPolarizationWprime</tt> 397 !% should be set to the direction defined by <i>A</i><math>^{-1}</math><i>pol3</i>. 398 !% For more information see MJT Oliveira 399 !% <i>et al.</i>, <i>J. Nanoscience and Nanotechnology</i> <b>8</b>, 400 !% 3392 (2008). 401 !%End 402 if(parse_block(namespace, 'TDPolarizationWprime', blk)==0) then 403 do idir = 1, 3 404 call parse_block_float(blk, 0, idir - 1, kick%wprime(idir)) 405 end do 406 kick%wprime(1:3) = kick%wprime(1:3) / sqrt(sum(kick%wprime(1:3)**2)) 407 call parse_block_end(blk) 408 else 409 kick%wprime(1:3) = (/ M_ZERO, M_ZERO, M_ONE /) 410 end if 411 end if 412 413 ! for non-dipole, it is more complicated to check whether it is actually in the periodic direction 414 if(periodic_dim > 0 .and. kick%delta_strength_mode /= KICK_MAGNON_MODE) then 415 message(1) = "Kicks cannot be applied correctly in periodic directions." 416 call messages_warning(1, namespace=namespace) 417 end if 418 419 !%Variable TDMomentumTransfer 420 !%Type block 421 !%Section Time-Dependent::Response 422 !%Description 423 !% Momentum-transfer vector for the calculation of the dynamic structure factor. 424 !% When this variable is set, a non-dipole field is applied, and an output file 425 !% <tt>ftchd</tt> is created (it contains the Fourier transform of the charge density 426 !% at each time). The type of the applied external field can be set by 427 !% an optional last number. Possible options are <tt>qexp</tt> (default), <tt>qcos</tt>, 428 !% <tt>qsin</tt>, or <tt>qcos+qsin</tt>. In the formulae below, 429 !% <math>\vec{q}</math> is the momentum-transfer vector. 430 !%Option qexp 1 431 !% External field is <math>e^{i \vec{q} \cdot \vec{r}}</math>. 432 !%Option qcos 2 433 !% External field is <math>\cos \left( i \vec{q} \cdot \vec{r} \right)</math>. 434 !%Option qsin 3 435 !% External field is <math>\sin \left( i \vec{q} \cdot \vec{r} \right)</math>. 436 !%Option qbessel 4 437 !% External field is <math>j_l \left( \vec{q} \cdot \vec{r} \right) Y_{lm} \left(\vec{r} \right)</math>. 438 !% In this case, the block has to include two extra values (<i>l</i> and <i>m</i>). 439 !%End 440 441 if(parse_block(namespace, 'TDMomentumTransfer', blk)==0) then 442 kick%nqvec = 1 443 SAFE_ALLOCATE(kick%qvector(1:MAX_DIM,1)) 444 do idir = 1, MAX_DIM 445 call parse_block_float(blk, 0, idir - 1, kick%qvector(idir,1)) 446 kick%qvector(idir,1) = units_to_atomic(unit_one / units_inp%length, kick%qvector(idir,1)) 447 end do 448 449 ! Read the calculation mode (exp, cos, sin, or bessel) 450 if(parse_block_cols(blk, 0) > MAX_DIM) then 451 452 call parse_block_integer(blk, 0, idir - 1, kick%qkick_mode) 453 454 ! Read l and m if bessel mode (j_l*Y_lm) is used 455 if(kick%qkick_mode == QKICKMODE_BESSEL .and. parse_block_cols(blk, 0) == MAX_DIM+3) then 456 call parse_block_integer(blk, 0, idir + 0, kick%qbessel_l) 457 call parse_block_integer(blk, 0, idir + 1, kick%qbessel_m) 458 else 459 kick%qbessel_l = 0 460 kick%qbessel_m = 0 461 end if 462 463 else 464 kick%qkick_mode = QKICKMODE_EXP 465 end if 466 467 call parse_block_end(blk) 468 469 if(sb%kpoints%use_symmetries) then 470 do iop = 1, symmetries_number(sb%symm) 471 if(iop == symmetries_identity_index(sb%symm)) cycle 472 if(.not. symm_op_invariant_cart(sb%symm%ops(iop), kick%qvector(:,1), CNST(1e-5))) then 473 message(1) = "The TDMomentumTransfer breaks (at least) one of the symmetries used to reduce the k-points." 474 message(2) = "Set SymmetryBreakDir equal to TDMomemtumTransfer." 475 call messages_fatal(2, namespace=namespace) 476 end if 477 end do 478 end if 479 480 else 481 kick%qkick_mode = QKICKMODE_NONE 482 kick%nqvec = 1 483 SAFE_ALLOCATE(kick%qvector(1:MAX_DIM,1)) 484 kick%qvector(:,1) = M_ZERO 485 end if 486 487 kick%qlength = sqrt(sum(kick%qvector(:,1)**2)) 488 489 if(kick%delta_strength_mode == KICK_MAGNON_MODE) then 490 !%Variable TDEasyAxis 491 !%Type block 492 !%Section Time-Dependent::Response::Dipole 493 !%Description 494 !% For magnon kicks only. 495 !% This variable defines the direction of the easy axis of the crystal. 496 !% The magnetization is kicked in the plane transverse to this vector 497 !%End 498 if(parse_block(namespace, 'TDEasyAxis', blk)==0) then 499 n_rows = parse_block_n(blk) 500 501 do idir = 1, 3 502 call parse_block_float(blk, 0, idir - 1, kick%easy_axis(idir)) 503 end do 504 norm = sqrt(sum(kick%easy_axis(1:3)**2)) 505 if(norm < CNST(1e-9)) then 506 message(1) = "TDEasyAxis norm is too small." 507 call messages_fatal(1, namespace=namespace) 508 end if 509 kick%easy_axis(1:3) = kick%easy_axis(1:3)/norm 510 call parse_block_end(blk) 511 else 512 message(1) = "For magnons, the variable TDEasyAxis must be defined." 513 call messages_fatal(1, namespace=namespace) 514 end if 515 516 !We first two vectors defining a basis in the transverse plane 517 !For this we take two vectors not align with the first one 518 !and we perform a Gram-Schmidt orthogonalization 519 kick%trans_vec(1,1) = -kick%easy_axis(2) 520 kick%trans_vec(2,1) = M_TWO*kick%easy_axis(3) 521 kick%trans_vec(3,1) = M_THREE*kick%easy_axis(1) 522 523 dot = sum(kick%easy_axis(1:3)*kick%trans_vec(1:3,1)) 524 kick%trans_vec(1:3,1) = kick%trans_vec(1:3,1) - dot*kick%easy_axis(1:3) 525 norm = sum(kick%trans_vec(1:3,1)**2) 526 kick%trans_vec(1:3,1) = kick%trans_vec(1:3,1)/sqrt(norm) 527 528 !To get a direct basis, the last vector is obtained by the cross product 529 kick%trans_vec(1,2) = kick%easy_axis(2) * kick%trans_vec(3,1) - kick%easy_axis(3) * kick%trans_vec(2,1) 530 kick%trans_vec(2,2) = kick%easy_axis(3) * kick%trans_vec(1,1) - kick%easy_axis(1) * kick%trans_vec(3,1) 531 kick%trans_vec(3,2) = kick%easy_axis(1) * kick%trans_vec(2,1) - kick%easy_axis(2) * kick%trans_vec(1,1) 532 533 !The perturbation direction is defined as 534 !cos(q.r)*uvec + sin(q.r)*vvec 535 536 537 if(parse_is_defined(namespace, 'TDMomentumTransfer') & 538 .and. parse_is_defined(namespace, 'TDMultipleMomentumTransfer')) then 539 message(1) = "TDMomentumTransfer and TDMultipleMomentumTransfer cannot be defined at the same time." 540 call messages_fatal(1, namespace=namespace) 541 end if 542 543 if(parse_is_defined(namespace, 'TDMultipleMomentumTransfer')) then 544 545 kick%qkick_mode = QKICKMODE_EXP 546 547 !%Variable TDMultipleMomentumTransfer 548 !%Type block 549 !%Section Time-Dependent::Response 550 !%Description 551 !% For magnon kicks only. 552 !% A simple way to specify momentum-transfer vectors for the calculation of 553 !% the magnetization dynamics. This variable should be used for a supercell. 554 !% For each reciprocal lattice vectors, the code will kick the original magnetization 555 !% using all the multiples of it. 556 !% The syntax reads: 557 !% 558 !% <tt>%TDMultipleMomentumTransfer 559 !% <br> N_x | N_y | N_z 560 !% <br>%</tt> 561 !% 562 !% and will include the (2N_x+1)*(2N_y+1)*(2N_z+1) multiples vectors of the reciprocal 563 !% lattice vectors of the current cell. 564 !%End 565 if(parse_block(namespace, 'TDMultipleMomentumTransfer', blk) /= 0) then 566 write(message(1),'(a)') 'Internal error while reading TDMultipleMomentumTransfer.' 567 call messages_fatal(1, namespace=namespace) 568 end if 569 570 do idir = 1, 3 571 call parse_block_integer(blk, 0, idir-1, kick%nqmult(idir)) 572 end do 573 574 call parse_block_end(blk) 575 576 577 kick%nqvec = (2*kick%nqmult(1)+1)*(2*kick%nqmult(2)+1)*(2*kick%nqmult(3)+1) 578 !qvector has been allocated by default to a null vector before 579 SAFE_DEALLOCATE_A(kick%qvector) 580 SAFE_ALLOCATE(kick%qvector(1:MAX_DIM, 1:kick%nqvec)) 581 iq = 0 582 do iqx = -kick%nqmult(1), kick%nqmult(1) 583 do iqy = -kick%nqmult(2), kick%nqmult(2) 584 do iqz = -kick%nqmult(3), kick%nqmult(3) 585 iq = iq + 1 586 qtemp(1:3) = (/iqx, iqy, iqz/) 587 call kpoints_to_absolute(sb%klattice, qtemp, kick%qvector(1:3, iq), 3) 588 589 !Checking symmetries for all G vectors 590 if(sb%kpoints%use_symmetries) then 591 do iop = 1, symmetries_number(sb%symm) 592 if(iop == symmetries_identity_index(sb%symm)) cycle 593 if(.not. symm_op_invariant_cart(sb%symm%ops(iop), kick%qvector(:,iq), CNST(1e-5))) then 594 message(1) = "The TDMultipleMomentumTransfer breaks (at least) one " & 595 // "of the symmetries used to reduce the k-points." 596 message(2) = "Set SymmetryBreakDir accordingly." 597 call messages_fatal(2, namespace=namespace) 598 end if 599 end do 600 end if 601 end do 602 end do 603 end do 604 605 end if 606 607 else 608 kick%easy_axis(1:MAX_DIM) = M_ZERO 609 end if 610 611 if(kick%delta_strength_mode == KICK_MAGNON_MODE .and. kick%qkick_mode /= QKICKMODE_EXP) then 612 message(1) = "For magnons, the kick mode must be exponential." 613 call messages_fatal(1, namespace=namespace) 614 end if 615 616 POP_SUB(kick_init) 617 end subroutine kick_init 618 619 620 ! --------------------------------------------------------- 621 subroutine kick_copy(kick_out, kick_in) 622 type(kick_t), intent(inout) :: kick_out 623 type(kick_t), intent(in) :: kick_in 624 625 PUSH_SUB(kick_copy) 626 627 kick_out%time = kick_in%time 628 629 kick_out%delta_strength_mode = kick_in%delta_strength_mode 630 kick_out%delta_strength = kick_in%delta_strength 631 632 !> In case we use a normal dipole kick: 633 kick_out%pol(1:MAX_DIM, 1:MAX_DIM) = kick_in%pol(1:MAX_DIM, 1:MAX_DIM) 634 kick_out%pol_dir = kick_in%pol_dir 635 kick_out%pol_equiv_axes = kick_in%pol_equiv_axes 636 kick_out%wprime(1:MAX_DIM) = kick_in%wprime(1:MAX_DIM) 637 638 !> In case we have a general multipolar kick, 639 kick_out%n_multipoles = kick_in%n_multipoles 640 if (kick_out%n_multipoles > 0) then 641 SAFE_ALLOCATE(kick_out%l(1:kick_out%n_multipoles)) 642 SAFE_ALLOCATE(kick_out%m(1:kick_out%n_multipoles)) 643 SAFE_ALLOCATE(kick_out%weight(1:kick_out%n_multipoles)) 644 kick_out%l = kick_in%l 645 kick_out%m = kick_in%m 646 kick_out%weight = kick_in%weight 647 end if 648 kick_out%nqvec = kick_in%nqvec 649 SAFE_ALLOCATE(kick_out%qvector(1:MAX_DIM, 1:kick_in%nqvec)) 650 kick_out%qvector(1:MAX_DIM, 1:kick_in%nqvec) = kick_in%qvector(1:MAX_DIM, 1:kick_in%nqvec) 651 kick_out%qlength = kick_in%qlength 652 kick_out%qkick_mode = kick_in%qkick_mode 653 kick_out%qbessel_l = kick_in%qbessel_l 654 kick_out%qbessel_m = kick_in%qbessel_m 655 656 !> In case we use a general function 657 kick_out%function_mode = kick_in%function_mode 658 kick_out%user_defined_function = kick_in%user_defined_function 659 660 kick_out%easy_axis(1:MAX_DIM) = kick_in%easy_axis(1:MAX_DIM) 661 662 POP_SUB(kick_copy) 663 end subroutine kick_copy 664 665 ! --------------------------------------------------------- 666 subroutine kick_end(kick) 667 type(kick_t), intent(inout) :: kick 668 669 PUSH_SUB(kick_end) 670 671 kick%delta_strength_mode = 0 672 kick%pol_equiv_axes = 0 673 kick%pol = M_ZERO 674 kick%pol_dir = 0 675 kick%wprime = M_ZERO 676 if (kick%n_multipoles > 0) then 677 SAFE_DEALLOCATE_P(kick%l) 678 SAFE_DEALLOCATE_P(kick%m) 679 SAFE_DEALLOCATE_P(kick%weight) 680 end if 681 kick%n_multipoles = 0 682 kick%qkick_mode = QKICKMODE_NONE 683 SAFE_DEALLOCATE_A(kick%qvector) 684 kick%easy_axis(1:MAX_DIM) = M_ZERO 685 686 POP_SUB(kick_end) 687 end subroutine kick_end 688 689 ! --------------------------------------------------------- 690 subroutine kick_read(kick, iunit, namespace) 691 type(kick_t), intent(inout) :: kick 692 integer, intent(in) :: iunit 693 type(namespace_t), intent(in) :: namespace 694 695 integer :: im, ierr 696 character(len=100) :: line 697 698 PUSH_SUB(kick_read) 699 700 kick%function_mode = -1 701 702 read(iunit, '(15x,i2)') kick%delta_strength_mode 703 read(iunit, '(15x,f18.12)') kick%delta_strength 704 read(iunit, '(a)') line 705 if(index(line,'defined') /= 0) then 706 kick%function_mode = KICK_FUNCTION_USER_DEFINED 707 ! "# User defined: " 708 read(line,'(16x,a)') kick%user_defined_function 709 elseif(index(line,'multipole') /= 0) then 710 kick%function_mode = KICK_FUNCTION_MULTIPOLE 711 ! "# N multipoles " 712 read(line, '(15x,i3)') kick%n_multipoles 713 SAFE_ALLOCATE( kick%l(1:kick%n_multipoles)) 714 SAFE_ALLOCATE( kick%m(1:kick%n_multipoles)) 715 SAFE_ALLOCATE(kick%weight(1:kick%n_multipoles)) 716 do im = 1, kick%n_multipoles 717 ! "# multipole " 718 read(iunit, '(15x,2i3,f18.12)') kick%l(im), kick%m(im), kick%weight(im) 719 end do 720 else 721 kick%function_mode = KICK_FUNCTION_DIPOLE 722 kick%n_multipoles = 0 723 nullify(kick%l) 724 nullify(kick%m) 725 backspace(iunit) 726 727 read(iunit, '(15x,3f18.12)') kick%pol(1:3, 1) 728 read(iunit, '(15x,3f18.12)') kick%pol(1:3, 2) 729 read(iunit, '(15x,3f18.12)') kick%pol(1:3, 3) 730 read(iunit, '(15x,i2)') kick%pol_dir 731 read(iunit, '(15x,i2)') kick%pol_equiv_axes 732 read(iunit, '(15x,3f18.12)') kick%wprime(1:3) 733 end if 734 if(kick%delta_strength_mode == KICK_MAGNON_MODE) then 735 read(iunit, '(15x,i3)') kick%nqvec 736 SAFE_ALLOCATE(kick%qvector(1:MAX_DIM, 1:kick%nqvec)) 737 do im = 1, kick%nqvec 738 read(iunit, '(15x,3f18.12)') kick%qvector(1:3, im) 739 end do 740 read(iunit, '(15x,3f18.12)') kick%easy_axis(1:3) 741 read(iunit, '(15x,3f18.12)') kick%trans_vec(1:3,1) 742 read(iunit, '(15x,3f18.12)') kick%trans_vec(1:3,2) 743 end if 744 read(iunit, '(15x,f18.12)', iostat = ierr) kick%time 745 746 if(ierr /= 0) then 747 kick%time = M_ZERO 748 backspace(iunit) 749 end if 750 751 if (kick%function_mode < 0) then 752 message(1) = "No kick could be read from file." 753 call messages_fatal(1, namespace=namespace) 754 end if 755 756 POP_SUB(kick_read) 757 end subroutine kick_read 758 759 760 ! --------------------------------------------------------- 761 subroutine kick_write(kick, iunit, out) 762 type(kick_t), intent(in) :: kick 763 integer, optional, intent(in) :: iunit 764 type(c_ptr), optional, intent(inout) :: out 765 766 integer :: im 767 character(len=120) :: aux 768 769 PUSH_SUB(kick_write) 770 771 if(present(iunit)) then 772 write(iunit, '(a15,i1)') '# kick mode ', kick%delta_strength_mode 773 write(iunit, '(a15,f18.12)') '# kick strength', kick%delta_strength 774 ! if this were to be read by humans, we would want units_from_atomic(units_out%length**(-1)) 775 if(kick%function_mode == KICK_FUNCTION_USER_DEFINED) then 776 write(iunit,'(a15,1x,a)') '# User defined:', trim(kick%user_defined_function) 777 elseif(kick%n_multipoles > 0) then 778 write(iunit, '(a15,i3)') '# N multipoles ', kick%n_multipoles 779 do im = 1, kick%n_multipoles 780 write(iunit, '(a15,2i3,f18.12)') '# multipole ', kick%l(im), kick%m(im), kick%weight(im) 781 end do 782 else 783 write(iunit, '(a15,3f18.12)') '# pol(1) ', kick%pol(1:3, 1) 784 write(iunit, '(a15,3f18.12)') '# pol(2) ', kick%pol(1:3, 2) 785 write(iunit, '(a15,3f18.12)') '# pol(3) ', kick%pol(1:3, 3) 786 write(iunit, '(a15,i1)') '# direction ', kick%pol_dir 787 write(iunit, '(a15,i1)') '# Equiv. axes ', kick%pol_equiv_axes 788 write(iunit, '(a15,3f18.12)') '# wprime ', kick%wprime(1:3) 789 end if 790 write(iunit, '(a15,f18.12)') "# kick time ", kick%time 791 792 else if(present(out)) then 793 write(aux, '(a15,i2)') '# kick mode ', kick%delta_strength_mode 794 call write_iter_string(out, aux) 795 call write_iter_nl(out) 796 write(aux, '(a15,f18.12)') '# kick strength', kick%delta_strength 797 call write_iter_string(out, aux) 798 call write_iter_nl(out) 799 if(kick%function_mode == KICK_FUNCTION_USER_DEFINED) then 800 write(aux,'(a15,1x,a)') '# User defined:', trim(kick%user_defined_function) 801 call write_iter_string(out, aux) 802 call write_iter_nl(out) 803 elseif(kick%n_multipoles > 0) then 804 write(aux, '(a15,i3)') '# N multipoles ', kick%n_multipoles 805 call write_iter_string(out, aux) 806 call write_iter_nl(out) 807 do im = 1, kick%n_multipoles 808 write(aux, '(a15,2i3,f18.12)') '# multipole ', kick%l(im), kick%m(im), kick%weight(im) 809 call write_iter_string(out, aux) 810 call write_iter_nl(out) 811 end do 812 else 813 write(aux, '(a15,3f18.12)') '# pol(1) ', kick%pol(1:3, 1) 814 call write_iter_string(out, aux) 815 call write_iter_nl(out) 816 write(aux, '(a15,3f18.12)') '# pol(2) ', kick%pol(1:3, 2) 817 call write_iter_string(out, aux) 818 call write_iter_nl(out) 819 write(aux, '(a15,3f18.12)') '# pol(3) ', kick%pol(1:3, 3) 820 call write_iter_string(out, aux) 821 call write_iter_nl(out) 822 write(aux, '(a15,i2)') '# direction ', kick%pol_dir 823 call write_iter_string(out, aux) 824 call write_iter_nl(out) 825 write(aux, '(a15,i2)') '# Equiv. axes ', kick%pol_equiv_axes 826 call write_iter_string(out, aux) 827 call write_iter_nl(out) 828 write(aux, '(a15,3f18.12)') '# wprime ', kick%wprime(1:3) 829 call write_iter_string(out, aux) 830 call write_iter_nl(out) 831 end if 832 if(present(out) .and. kick%delta_strength_mode == KICK_MAGNON_MODE) then 833 write(aux, '(a15,i3)') '# N q-vectors ', kick%nqvec 834 call write_iter_string(out, aux) 835 call write_iter_nl(out) 836 do im = 1, kick%nqvec 837 write(aux, '(a15,3f18.12)') '# q-vector ', kick%qvector(1:3, im) 838 call write_iter_string(out, aux) 839 call write_iter_nl(out) 840 end do 841 write(aux, '(a15,3f18.12)') '# Easy axis ', kick%easy_axis(1:3) 842 call write_iter_string(out, aux) 843 call write_iter_nl(out) 844 write(aux, '(a15,3f18.12)') '# Trans. dir 1 ', kick%trans_vec(1:3,1) 845 call write_iter_string(out, aux) 846 call write_iter_nl(out) 847 write(aux, '(a15,3f18.12)') '# Trans. dir 2 ', kick%trans_vec(1:3,2) 848 call write_iter_string(out, aux) 849 call write_iter_nl(out) 850 end if 851 write(aux, '(a15,f18.12)') "# kick time ", kick%time 852 call write_iter_string(out, aux) 853 call write_iter_nl(out) 854 855 end if 856 857 POP_SUB(kick_write) 858 end subroutine kick_write 859 860 861 ! --------------------------------------------------------- 862 ! 863 subroutine kick_function_get(mesh, kick, kick_function, iq, to_interpolate) 864 type(mesh_t), intent(in) :: mesh 865 type(kick_t), intent(in) :: kick 866 CMPLX, intent(out) :: kick_function(:) 867 integer, intent(in) :: iq 868 logical, optional, intent(in) :: to_interpolate 869 870 integer :: ip, im 871 FLOAT :: xx(MAX_DIM) 872 FLOAT :: rkick, ikick, rr, ylm 873 874 integer :: np 875 876 PUSH_SUB(kick_function_get) 877 878 np = mesh%np 879 if(present(to_interpolate)) then 880 if(to_interpolate) np = mesh%np_part 881 end if 882 883 if(abs(kick%qlength) > M_EPSILON .or. kick%delta_strength_mode == KICK_MAGNON_MODE) then ! q-vector is set 884 885 select case (kick%qkick_mode) 886 case (QKICKMODE_COS) 887 write(message(1), '(a,3F9.5,a)') 'Info: Using cos(q.r) field with q = (', kick%qvector(1:3, iq), ')' 888 case (QKICKMODE_SIN) 889 write(message(1), '(a,3F9.5,a)') 'Info: Using sin(q.r) field with q = (', kick%qvector(1:3, iq), ')' 890 case (QKICKMODE_SIN + QKICKMODE_COS) 891 write(message(1), '(a,3F9.5,a)') 'Info: Using sin(q.r)+cos(q.r) field with q = (', kick%qvector(1:3, iq), ')' 892 case (QKICKMODE_EXP) 893 write(message(1), '(a,3F9.5,a)') 'Info: Using exp(iq.r) field with q = (', kick%qvector(1:3, iq), ')' 894 case (QKICKMODE_BESSEL) 895 write(message(1), '(a,I2,a,I2,a,F9.5)') 'Info: Using j_l(qr)*Y_lm(r) field with (l,m)= (', & 896 kick%qbessel_l, ",", kick%qbessel_m,') and q = ', kick%qlength 897 case default 898 write(message(1), '(a,3F9.6,a)') 'Info: Unknown field type!' 899 end select 900 call messages_info(1) 901 902 kick_function = M_z0 903 do ip = 1, np 904 call mesh_r(mesh, ip, rr, coords = xx) 905 select case (kick%qkick_mode) 906 case (QKICKMODE_COS) 907 kick_function(ip) = kick_function(ip) + cos(sum(kick%qvector(:, iq) * xx(:))) 908 case (QKICKMODE_SIN) 909 kick_function(ip) = kick_function(ip) + sin(sum(kick%qvector(:, iq) * xx(:))) 910 case (QKICKMODE_SIN+QKICKMODE_COS) 911 kick_function(ip) = kick_function(ip) + sin(sum(kick%qvector(:, iq) * xx(:))) 912 case (QKICKMODE_EXP) 913 kick_function(ip) = kick_function(ip) + exp(M_zI * sum(kick%qvector(:, iq) * xx(:))) 914 case (QKICKMODE_BESSEL) 915 call grylmr(mesh%x(ip, 1), mesh%x(ip, 2), mesh%x(ip, 3), kick%qbessel_l, kick%qbessel_m, ylm) 916 kick_function(ip) = kick_function(ip) + loct_sph_bessel(kick%qbessel_l, kick%qlength*sqrt(sum(xx(:)**2)))*ylm 917 end select 918 end do 919 920 else 921 if(kick%function_mode == KICK_FUNCTION_USER_DEFINED) then 922 923 kick_function = M_z0 924 do ip = 1, np 925 call mesh_r(mesh, ip, rr, coords = xx) 926 rkick = M_ZERO; ikick = M_ZERO 927 call parse_expression(rkick, ikick, mesh%sb%dim, xx, rr, M_ZERO, trim(kick%user_defined_function)) 928 kick_function(ip) = rkick 929 end do 930 931 elseif(kick%n_multipoles > 0) then 932 933 kick_function = M_z0 934 do im = 1, kick%n_multipoles 935 do ip = 1, np 936 call mesh_r(mesh, ip, rr, coords = xx) 937 call loct_ylm(1, xx(1), xx(2), xx(3), kick%l(im), kick%m(im), ylm) 938 kick_function(ip) = kick_function(ip) + kick%weight(im) * (rr**kick%l(im)) * ylm 939 end do 940 end do 941 else 942 do ip = 1, np 943 kick_function(ip) = sum(mesh%x(ip, 1:mesh%sb%dim) * & 944 kick%pol(1:mesh%sb%dim, kick%pol_dir)) 945 end do 946 end if 947 end if 948 949 POP_SUB(kick_function_get) 950 end subroutine kick_function_get 951 952 953 ! --------------------------------------------------------- 954 ! 955 subroutine kick_pcm_function_get(mesh, kick, psolver, pcm, kick_pcm_function) 956 type(mesh_t), intent(in) :: mesh 957 type(kick_t), intent(in) :: kick 958 type(poisson_t), intent(in) :: psolver 959 type(pcm_t), intent(inout) :: pcm 960 CMPLX, intent(out) :: kick_pcm_function(:) 961 962 CMPLX, allocatable :: kick_function_interpolate(:) 963 FLOAT, allocatable :: kick_function_real(:) 964 965 PUSH_SUB(kick_pcm_function_get) 966 967 kick_pcm_function = M_ZERO 968 if ( pcm%localf ) then 969 SAFE_ALLOCATE(kick_function_interpolate(1:mesh%np_part)) 970 kick_function_interpolate = M_ZERO 971 call kick_function_get(mesh, kick, kick_function_interpolate, 1, to_interpolate = .true.) 972 SAFE_ALLOCATE(kick_function_real(1:mesh%np_part)) 973 kick_function_real = DREAL(kick_function_interpolate) 974 if ( pcm%kick_like ) then 975 ! computing kick-like polarization due to kick 976 call pcm_calc_pot_rs(pcm, mesh, psolver, kick = kick%delta_strength * kick_function_real, kick_time = .true.) 977 else if ( .not.pcm%kick_like .and. pcm%which_eps == PCM_DEBYE_MODEL) then 978 ! computing the kick-like part of polarization due to kick for Debye dielectric model 979 pcm%kick_like = .true. 980 call pcm_calc_pot_rs(pcm, mesh, psolver, kick = kick%delta_strength * kick_function_real, kick_time = .true.) 981 pcm%kick_like = .false. 982 else if ( .not.pcm%kick_like .and. pcm%which_eps == PCM_DRUDE_MODEL ) then 983 POP_SUB(kick_pcm_function_get) 984 return 985 end if 986 kick_pcm_function = pcm%v_kick_rs / kick%delta_strength 987 end if 988 989 POP_SUB(kick_pcm_function_get) 990 end subroutine kick_pcm_function_get 991 992 993 ! --------------------------------------------------------- 994 !> Applies the delta-function electric field \f$ E(t) = E_0 \Delta(t) \f$ 995 !! where \f$ E_0 = \frac{- k \hbar}{e} \f$ k = kick\%delta_strength. 996 subroutine kick_apply(mesh, st, ions, geo, kick, psolver, pcm) 997 type(mesh_t), intent(in) :: mesh 998 type(states_elec_t), intent(inout) :: st 999 type(ion_dynamics_t), intent(in) :: ions 1000 type(geometry_t), intent(inout) :: geo 1001 type(kick_t), intent(in) :: kick 1002 type(poisson_t), intent(in) :: psolver 1003 type(pcm_t), optional, intent(inout) :: pcm 1004 1005 integer :: iqn, ist, idim, ip, ispin, iatom 1006 CMPLX :: cc(2), kick_value 1007 CMPLX, allocatable :: kick_function(:), psi(:, :) 1008 1009 CMPLX, allocatable :: kick_pcm_function(:) 1010 integer :: ns, iq 1011 FLOAT :: uvec(MAX_DIM), vvec(MAX_DIM), Gvec(MAX_DIM,MAX_DIM) 1012 FLOAT :: xx(MAX_DIM), rr 1013 1014 PUSH_SUB(kick_apply) 1015 1016 ! The wavefunctions at time delta t read 1017 ! psi(delta t) = psi(t) exp(i k x) 1018 delta_strength: if(kick%delta_strength /= M_ZERO) then 1019 1020 SAFE_ALLOCATE(kick_function(1:mesh%np)) 1021 if(kick%delta_strength_mode /= KICK_MAGNON_MODE .or. kick%nqvec == 1) then 1022 call kick_function_get(mesh, kick, kick_function, 1) 1023 end if 1024 1025 ! PCM - computing polarization due to kick 1026 if( present(pcm) ) then 1027 SAFE_ALLOCATE(kick_pcm_function(1:mesh%np)) 1028 call kick_pcm_function_get(mesh, kick, psolver, pcm, kick_pcm_function) 1029 kick_function = kick_function + kick_pcm_function 1030 end if 1031 1032 write(message(1),'(a,f11.6)') 'Info: Applying delta kick: k = ', kick%delta_strength 1033 select case (kick%function_mode) 1034 case (KICK_FUNCTION_DIPOLE) 1035 message(2) = "Info: kick function: dipole." 1036 case (KICK_FUNCTION_MULTIPOLE) 1037 message(2) = "Info: kick function: multipoles." 1038 case (KICK_FUNCTION_USER_DEFINED) 1039 message(2) = "Info: kick function: user defined function." 1040 end select 1041 select case (kick%delta_strength_mode) 1042 case (KICK_DENSITY_MODE) 1043 message(3) = "Info: Delta kick mode: Density mode" 1044 case (KICK_SPIN_MODE) 1045 message(3) = "Info: Delta kick mode: Spin mode" 1046 case (KICK_SPIN_DENSITY_MODE) 1047 message(3) = "Info: Delta kick mode: Density + Spin modes" 1048 end select 1049 call messages_info(3) 1050 1051 ns = 1 1052 if(st%d%nspin == 2) ns = 2 1053 1054 SAFE_ALLOCATE(psi(1:mesh%np, 1:st%d%dim)) 1055 1056 if(kick%delta_strength_mode /= KICK_MAGNON_MODE) then 1057 1058 do iqn = st%d%kpt%start, st%d%kpt%end 1059 do ist = st%st_start, st%st_end 1060 call states_elec_get_state(st, mesh, ist, iqn, psi) 1061 1062 select case (kick%delta_strength_mode) 1063 case (KICK_DENSITY_MODE) 1064 do idim = 1, st%d%dim 1065 do ip = 1, mesh%np 1066 psi(ip, idim) = exp(M_zI*kick%delta_strength*kick_function(ip))*psi(ip, idim) 1067 end do 1068 end do 1069 1070 case (KICK_SPIN_MODE) 1071 ispin = states_elec_dim_get_spin_index(st%d, iqn) 1072 do ip = 1, mesh%np 1073 kick_value = M_zI*kick%delta_strength*kick_function(ip) 1074 1075 cc(1) = exp(kick_value) 1076 cc(2) = exp(-kick_value) 1077 1078 select case (st%d%ispin) 1079 case (SPIN_POLARIZED) 1080 psi(ip, 1) = cc(ispin)*psi(ip, 1) 1081 case (SPINORS) 1082 psi(ip, 1) = cc(1)*psi(ip, 1) 1083 psi(ip, 2) = cc(2)*psi(ip, 2) 1084 end select 1085 end do 1086 1087 case (KICK_SPIN_DENSITY_MODE) 1088 do ip = 1, mesh%np 1089 kick_value = M_zI*kick%delta_strength*kick_function(ip) 1090 cc(1) = exp(M_TWO*kick_value) 1091 1092 select case (st%d%ispin) 1093 case (SPIN_POLARIZED) 1094 if(is_spin_up(iqn)) then 1095 psi(ip, 1) = cc(1)*psi(ip, 1) 1096 end if 1097 case (SPINORS) 1098 psi(ip, 1) = cc(1)*psi(ip, 1) 1099 end select 1100 end do 1101 end select 1102 1103 call states_elec_set_state(st, mesh, ist, iqn, psi) 1104 1105 end do 1106 end do 1107 1108 else 1109 ASSERT(st%d%ispin==SPINORS) 1110 1111 if(kick%nqvec == 1) then 1112 !The perturbation direction is defined as 1113 !cos(q.r)*uvec + sin(q.r)*vvec 1114 uvec(1:3) = kick%trans_vec(1:3,1) 1115 vvec(1:3) = kick%trans_vec(1:3,2) 1116 1117 do iqn = st%d%kpt%start, st%d%kpt%end, ns 1118 do ist = st%st_start, st%st_end 1119 1120 call states_elec_get_state(st, mesh, ist, iqn, psi) 1121 1122 do ip = 1, mesh%np 1123 1124 cc(1) = psi(ip, 1) 1125 cc(2) = psi(ip, 2) 1126 1127 !First part: 1I*cos(\lambda) 1128 psi(ip, 1) = cos(kick%delta_strength)* cc(1) 1129 psi(ip, 2) = cos(kick%delta_strength)* cc(2) 1130 1131 !We now add -i sin(\lambda) u.\sigma 1132 ! (u_z u_x-i*u_y) (v_z v_x-i*v_y) 1133 ! =cos(q.r) ( ) + sin(q.r)( ) 1134 ! (u_x+i*u_y -u_z ) (v_x+i*v_y -v_z ) 1135 psi(ip, 1) = psi(ip, 1) -M_zI*sin(kick%delta_strength)*( TOFLOAT(kick_function(ip)) & 1136 * (uvec(3)*cc(1) + (uvec(1)-M_zI*uvec(2))*cc(2)) & 1137 + aimag(kick_function(ip)) * (vvec(3)*cc(1) + (vvec(1)-M_zI*vvec(2))*cc(2))) 1138 psi(ip, 2) = psi(ip, 2) -M_zI*sin(kick%delta_strength)*( TOFLOAT(kick_function(ip)) & 1139 * (-uvec(3)*cc(2) + (uvec(1)+M_zI*uvec(2))*cc(1)) & 1140 + aimag(kick_function(ip)) * (-vvec(3)*cc(2) + (vvec(1)+M_zI*vvec(2))*cc(1))) 1141 1142 end do 1143 1144 call states_elec_set_state(st, mesh, ist, iqn, psi) 1145 1146 end do 1147 end do 1148 1149 else ! Multi-q kick 1150 1151 call kpoints_to_absolute(mesh%sb%klattice, (/M_ONE,M_ZERO,M_ZERO/), Gvec(1:3, 1), 3) 1152 call kpoints_to_absolute(mesh%sb%klattice, (/M_ZERO,M_ONE,M_ZERO/), Gvec(1:3, 2), 3) 1153 call kpoints_to_absolute(mesh%sb%klattice, (/M_ZERO,M_ZERO,M_ONE/), Gvec(1:3, 3), 3) 1154 1155 kick_function = M_ONE 1156 do ip = 1, mesh%np 1157 call mesh_r(mesh, ip, rr, coords = xx) 1158 do iq = 1, 3 1159 if(kick%nqmult(iq) == 0) cycle 1160 if(abs(sin(M_HALF*sum(Gvec(1:3, iq) * xx(1:3)))) <= M_EPSILON) cycle 1161 1162 kick_function(ip) = kick_function(ip)*sin(M_HALF*(2*kick%nqmult(iq)+1) & 1163 *sum(Gvec(1:3, iq) * xx(1:3)))/sin(M_HALF*sum(Gvec(1:3, iq) * xx(1:3))) 1164 end do 1165 kick_function(ip) = kick_function(ip)*kick%delta_strength 1166 end do 1167 1168 do iqn = st%d%kpt%start, st%d%kpt%end, ns 1169 do ist = st%st_start, st%st_end 1170 1171 call states_elec_get_state(st, mesh, ist, iqn, psi) 1172 1173 do ip = 1, mesh%np 1174 1175 cc(1) = psi(ip, 1) 1176 cc(2) = psi(ip, 2) 1177 1178 ! (cos(F) + in_x sin(F) sin(F)(u_y (u_x-iu_y)/(1+u_z) - iu_z)) 1179 ! = ( ) 1180 ! (-sin(F)(u_y (u_x+iu_y)/(1+u_z)+iu_z) cos(F) - in_x sin(F) ) 1181 1182 psi(ip, 1) = (cos(kick_function(ip))+M_zI*kick%easy_axis(1)*sin(kick_function(ip)))*cc(1) & 1183 +sin(kick_function(ip))*(kick%easy_axis(2)*(kick%easy_axis(1) & 1184 -M_zI*kick%easy_axis(2))/(1+kick%easy_axis(3))-M_zI*kick%easy_axis(3))*cc(2) 1185 psi(ip, 2) =-sin(kick_function(ip))*(kick%easy_axis(2)*(kick%easy_axis(1) & 1186 +M_zI*kick%easy_axis(2))/(1+kick%easy_axis(3))+M_zI*kick%easy_axis(3))*cc(1) & 1187 + (cos(kick_function(ip))-m_zI*kick%easy_axis(1)*sin(kick_function(ip)))*cc(2) 1188 1189 end do 1190 1191 call states_elec_set_state(st, mesh, ist, iqn, psi) 1192 1193 end do 1194 end do 1195 1196 end if 1197 end if 1198 1199 SAFE_DEALLOCATE_A(psi) 1200 1201 ! The nuclear velocities will be changed by 1202 ! Delta v_z = ( Z*e*E_0 / M) = - ( Z*k*\hbar / M) 1203 ! where M and Z are the ionic mass and charge, respectively. 1204 if(ion_dynamics_ions_move(ions) .and. kick%delta_strength /= M_ZERO) then 1205 if(kick%delta_strength_mode /= KICK_MAGNON_MODE) then 1206 do iatom = 1, geo%natoms 1207 geo%atom(iatom)%v(1:mesh%sb%dim) = geo%atom(iatom)%v(1:mesh%sb%dim) + & 1208 kick%delta_strength * kick%pol(1:mesh%sb%dim, kick%pol_dir) * & 1209 P_PROTON_CHARGE * species_zval(geo%atom(iatom)%species) / & 1210 species_mass(geo%atom(iatom)%species) 1211 end do 1212 end if 1213 end if 1214 1215 SAFE_DEALLOCATE_A(kick_function) 1216 end if delta_strength 1217 1218 POP_SUB(kick_apply) 1219 end subroutine kick_apply 1220 1221 pure integer function kick_get_type(kick) result(kick_type) 1222 type(kick_t), intent(in) :: kick 1223 1224 kick_type = kick%delta_strength_mode 1225 1226 end function kick_get_type 1227 1228end module kick_oct_m 1229 1230!! Local Variables: 1231!! mode: f90 1232!! coding: utf-8 1233!! End: 1234