Line data Source code
1 : !!****f* ABINIT/m_xc_noncoll
2 : !! NAME
3 : !! m_xc_noncoll
4 : !!
5 : !! FUNCTION
6 : !! This module provides several routines used in non-collinear XC routines
7 : !! (rotation of the magnetization in order to align it)
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2001-2026 ABINIT group (EB, MT, FR, SPr)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! INPUTS
16 : !!
17 : !! OUTPUT
18 : !!
19 : !! SOURCE
20 :
21 : #if defined HAVE_CONFIG_H
22 : #include "config.h"
23 : #endif
24 :
25 : #include "abi_common.h"
26 :
27 : MODULE m_xc_noncoll
28 :
29 : use defs_basis
30 : use m_abicore
31 : use m_errors
32 :
33 : implicit none
34 :
35 : private
36 :
37 : ! public procedures
38 : public :: rotate_mag ! Rotate a non-collinear density wrt a magnetization
39 : public :: rotate_back_mag ! Rotate back a collinear XC potential wrt a magnetization
40 : public :: rotate_back_mag_dfpt ! Rotate back a collinear 1st-order XC potential wrt a magnetization
41 : public :: test_rotations ! test whether methods in rotate_back_mag_dfpt give similar results
42 :
43 : !Tolerance on magnetization norm
44 : real(dp),parameter :: m_norm_min=tol8
45 :
46 : !Default rotation method for DFPT
47 : integer,parameter :: rotation_method_default=3
48 :
49 : CONTAINS
50 :
51 : !===========================================================
52 : !!***
53 :
54 : !!****t* m_xc_noncoll/rotate_mag
55 : !! NAME
56 : !! rotate_mag
57 : !!
58 : !! FUNCTION
59 : !! Project (rotate) a non-collinear density (stored as density+magn.)
60 : !! on a magnetization and give a collinear density (stored as [up,dn] or [up+dn,up]).
61 : !! Align both z-axis.
62 : !!
63 : !!
64 : !! INPUTS
65 : !! rho_in(vectsize,4)=input non-collinear density and magnetization (1st or 0th order)
66 : !! mag(vectsize,3)=gs magnetization used for projection (0th order magnetization)
67 : !! vectsize=size of vector fields
68 : !! [mag_norm_in(vectsize)]= --optional-- norm of mag(:) at each point of the grid
69 : !! [rho_out_format]= 1=rho_out is stored as [up,dn]
70 : !! 2=rho_out is stored as [up+dn,up]
71 : !! Default=1
72 : ! OUTPUT
73 : !! rho_out(vectsize,2)=output (projected, collinear) (1st order if rho_in is 1st order NC density matrix)
74 : !! [mag_norm_out(vectsize)]= --optional-- norm of mag(:) at each point of the grid
75 : !!
76 : !! Explicit formulae:
77 : !! rho_out_format=1
78 : !! rho_out(1) = half*( rho_in(1) + (mag,rho_in(2:4))/|mag|) // where (*,*) is scalar product
79 : !! rho_out(2) = half*( rho_in(1) - (mag,rho_in(2:4))/|mag|)
80 : !!
81 : !! rho_out_format=2
82 : !! rho_out(1) = rho_in(1)
83 : !! rho_out(2) = half*( rho_in(1) + (mag,rho_in(2:4))/|mag|)
84 : !!
85 : !!
86 : !! SOURCE
87 :
88 34772 : subroutine rotate_mag(rho_in,rho_out,mag,vectsize,cplex,&
89 : & mag_norm_in,mag_norm_out,rho_out_format) ! optional arguments
90 :
91 : !Arguments ------------------------------------
92 : !scalars
93 : integer,intent(in) :: vectsize
94 : integer,intent(in) :: cplex
95 : integer,intent(in),optional :: rho_out_format
96 : !arrays
97 : real(dp),intent(in) :: rho_in(cplex*vectsize,4),mag(vectsize,3)
98 : real(dp),intent(out) :: rho_out(cplex*vectsize,2)
99 : real(dp),intent(in),optional :: mag_norm_in(vectsize)
100 : real(dp),intent(out),optional :: mag_norm_out(vectsize)
101 :
102 : !Local variables-------------------------------
103 : !scalars
104 : integer :: ipt
105 : logical :: has_mag_norm,out_mag_norm
106 : real(dp) :: m_norm,mm,rho_up,rhoin_dot_mag
107 : real(dp) :: rhoin_dot_mag_re,rhoin_dot_mag_im
108 : !arrays
109 :
110 : ! *************************************************************************
111 :
112 : !DBG_ENTER("COLL")
113 :
114 34772 : has_mag_norm=present(mag_norm_in)
115 34772 : out_mag_norm=present(mag_norm_out)
116 :
117 34772 : if(cplex==1) then
118 38760325 : do ipt=1,vectsize
119 38726066 : if (has_mag_norm) then
120 0 : m_norm=mag_norm_in(ipt)
121 : else
122 38726066 : m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2)
123 : end if
124 :
125 : rhoin_dot_mag=rho_in(ipt,2)*mag(ipt,1)+rho_in(ipt,3)*mag(ipt,2) &
126 38726066 : & +rho_in(ipt,4)*mag(ipt,3)
127 :
128 38726066 : if(m_norm>m_norm_min)then
129 18094210 : mm=rhoin_dot_mag/m_norm
130 18094210 : rho_out(ipt,1)=half*(rho_in(ipt,1)+mm)
131 18094210 : rho_out(ipt,2)=half*(rho_in(ipt,1)-mm)
132 : else
133 20631856 : rho_out(ipt,1)=half*rho_in(ipt,1)
134 20631856 : rho_out(ipt,2)=half*rho_in(ipt,1)
135 : end if
136 :
137 38760325 : if (out_mag_norm) then
138 5895866 : if (m_norm >m_norm_min) mag_norm_out(ipt)=m_norm
139 5895866 : if (m_norm<=m_norm_min) mag_norm_out(ipt)=zero
140 : end if
141 :
142 : end do
143 :
144 : else ! cplex==2
145 :
146 2345473 : do ipt=1,vectsize
147 2344960 : if (has_mag_norm) then
148 0 : m_norm=mag_norm_in(ipt)
149 : else
150 2344960 : m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2)
151 : end if
152 :
153 : ! real part of m.m^(1)
154 : rhoin_dot_mag_re=rho_in(2*ipt-1,2)*mag(ipt,1)+rho_in(2*ipt-1,3)*mag(ipt,2) &
155 2344960 : & +rho_in(2*ipt-1,4)*mag(ipt,3)
156 : ! imaginary part of m.m^(1)
157 : rhoin_dot_mag_im=rho_in(2*ipt ,2)*mag(ipt,1)+rho_in(2*ipt ,3)*mag(ipt,2) &
158 2344960 : & +rho_in(2*ipt ,4)*mag(ipt,3)
159 2344960 : if(m_norm>m_norm_min)then
160 1684571 : mm=rhoin_dot_mag_re/m_norm
161 1684571 : rho_out(2*ipt-1,1)=half*(rho_in(2*ipt-1,1)+mm)
162 1684571 : rho_out(2*ipt-1,2)=half*(rho_in(2*ipt-1,1)-mm)
163 1684571 : mm=rhoin_dot_mag_im/m_norm
164 1684571 : rho_out(2*ipt ,1)=half*(rho_in(2*ipt ,1)+mm)
165 1684571 : rho_out(2*ipt ,2)=half*(rho_in(2*ipt ,1)-mm)
166 : else
167 660389 : rho_out(2*ipt-1,1)=half*rho_in(2*ipt-1,1)
168 : !MR: BUG FIXED HERE
169 : ! rho_out(2*ipt-1,2)=half*rho_in(2*ipt-1,2)
170 660389 : rho_out(2*ipt-1,2)=half*rho_in(2*ipt-1,1)
171 660389 : rho_out(2*ipt ,1)=half*rho_in(2*ipt ,1)
172 : !MR: BUG FIXED HERE
173 : ! rho_out(2*ipt ,2)=half*rho_in(2*ipt ,2)
174 660389 : rho_out(2*ipt ,2)=half*rho_in(2*ipt ,1)
175 : end if
176 :
177 2345473 : if (out_mag_norm) then
178 2344960 : if (m_norm >m_norm_min) mag_norm_out(ipt)=m_norm
179 2344960 : if (m_norm<=m_norm_min) mag_norm_out(ipt)=zero
180 : end if
181 :
182 : end do
183 :
184 : end if
185 :
186 34772 : if (present(rho_out_format)) then
187 1086 : if (rho_out_format==2) then
188 8009388 : do ipt=1,cplex*vectsize
189 8008302 : rho_up=rho_out(ipt,1)
190 8008302 : rho_out(ipt,1)=rho_up+rho_out(ipt,2)
191 8009388 : rho_out(ipt,2)=rho_up
192 : end do
193 : end if
194 : end if
195 :
196 : !DBG_EXIT("COLL")
197 :
198 34772 : end subroutine rotate_mag
199 : !!***
200 :
201 : !----------------------------------------------------------------------
202 :
203 : !!****t* m_xc_noncoll/rotate_back_mag
204 : !! NAME
205 : !! rotate_back_mag
206 : !!
207 : !! FUNCTION
208 : !! Rotate back a collinear XC potential (stored as up+dn) with respect to
209 : !! a magnetization and give a non-collinear XC potential
210 : !! (stored as up_up, dn_dn, Re[up_dn], Im[up_dn])
211 : !! Note: works only for cplex=1
212 : !!
213 : !! INPUTS
214 : !! vxc_in(vectsize,2)=input collinear XC potential
215 : !! mag(vectsize,3)=gs magnetization used for projection
216 : !! vectsize=size of vector fields
217 : !! [mag_norm_in(vectsize)]= --optional-- norm of mag(:) at each point of the grid
218 : !!
219 : !! OUTPUT
220 : !! vxc_out(vectsize,4)=output non-collinear XC potential
221 : !!
222 : !! SOURCE
223 :
224 71629 : subroutine rotate_back_mag(vxc_in,vxc_out,mag,vectsize,&
225 : & mag_norm_in) ! optional argument
226 :
227 : !Arguments ------------------------------------
228 : !scalars
229 : integer,intent(in) :: vectsize
230 : !arrays
231 : real(dp),intent(in) :: vxc_in(vectsize,2),mag(vectsize,3)
232 : real(dp),intent(out) :: vxc_out(vectsize,4)
233 : real(dp),intent(in),optional :: mag_norm_in(vectsize)
234 :
235 : !Local variables-------------------------------
236 : !scalars
237 : integer :: ipt
238 : logical :: has_mag_norm
239 : real(dp) :: dvdn,dvdz,m_norm
240 : !arrays
241 :
242 : ! *************************************************************************
243 :
244 : !DBG_ENTER("COLL")
245 :
246 71629 : has_mag_norm=present(mag_norm_in)
247 :
248 98954912 : do ipt=1,vectsize
249 98883283 : if (has_mag_norm) then
250 908847 : m_norm=mag_norm_in(ipt)
251 : else
252 97974436 : m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2)
253 : end if
254 :
255 98883283 : dvdn=half*(vxc_in(ipt,1)+vxc_in(ipt,2))
256 :
257 98954912 : if (m_norm>m_norm_min) then
258 62550572 : dvdz=half*(vxc_in(ipt,1)-vxc_in(ipt,2))/m_norm
259 62550572 : vxc_out(ipt,1)=dvdn+mag(ipt,3)*dvdz
260 62550572 : vxc_out(ipt,2)=dvdn-mag(ipt,3)*dvdz
261 62550572 : vxc_out(ipt,3)= mag(ipt,1)*dvdz
262 62550572 : vxc_out(ipt,4)=-mag(ipt,2)*dvdz
263 : else
264 108998133 : vxc_out(ipt,1:2)=dvdn
265 108998133 : vxc_out(ipt,3:4)=zero
266 : end if
267 : end do
268 :
269 : !DBG_EXIT("COLL")
270 :
271 71629 : end subroutine rotate_back_mag
272 : !!***
273 :
274 : !!****t* m_xc_noncoll/rotate_back_mag_dfpt
275 : !! NAME
276 : !! rotate_back_mag_dfpt
277 : !!
278 : !! FUNCTION
279 : !! Rotate back a 1st-order collinear XC potential (stored as up+dn) with respect to
280 : !! a magnetization and give a 1st-order non-collinear XC potential
281 : !! (stored as up_up, dn_dn, Re{up_dn}, Im{up_dn}).
282 : !!
283 : !! INPUTS
284 : !! mag(vectsize,3)=0-order magnetization used for projection
285 : !! rho1(vectsize,4)=1st-order non-collinear density and magnetization
286 : !! vxc(vectsize,4)=0-order non-collinear XC potential
287 : !! kxc(vectsize,nkxc)=0-order XC kernel (associated to vxc)
288 : !! vxc1_in(vectsize,2)=input 1st-order collinear XC potential
289 : !! vectsize=size of vector fields
290 : !! [mag_norm_in(vectsize)]= --optional-- norm of 0-order mag(:) at each point of the grid
291 : !! [rot_method]=Select method used to compute rotation matrix (1, 2 or 3)
292 : !! option=if 0, compute only the U^0 vxc^(1) U^0 part
293 : !! if 1, full first order xc potential
294 : !!
295 : !! NOTES
296 : !! cplex=1:
297 : !! V is stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
298 : !! N is stored as : n, m_x, m_y, m_z (real)
299 : !! cplex=2:
300 : !! V is stored as : V^11, V^22, V^12, i.V^21 (complex)
301 : !! N is stored as : n, m_x, m_y, mZ (complex)
302 : !!
303 : !! OUTPUT
304 : !! vxc1_out(vectsize,4)=output 1st-order non-collinear XC potential
305 : !!
306 : !! SOURCE
307 :
308 997 : subroutine rotate_back_mag_dfpt(option,vxc1_in,vxc1_out,vxc,kxc,rho1,mag,vectsize,cplex,&
309 : & mag_norm_in,rot_method) ! optional arguments
310 :
311 : !Arguments ------------------------------------
312 : !scalars
313 : integer,intent(in) :: vectsize
314 : integer,intent(in) :: cplex
315 : integer,intent(in) :: option
316 : integer,intent(in),optional :: rot_method
317 : !arrays
318 : real(dp),intent(in) :: kxc(:,:),mag(vectsize,3),vxc(vectsize,4)
319 : real(dp),intent(in) :: rho1(cplex*vectsize,4)
320 : real(dp),intent(in) :: vxc1_in(cplex*vectsize,2)
321 : real(dp),intent(in),optional :: mag_norm_in(vectsize)
322 : real(dp),intent(out) :: vxc1_out(cplex*vectsize,4)
323 :
324 : !Local variables-------------------------------
325 : !scalars
326 : integer :: ipt,rotation_method
327 : logical :: has_mag_norm
328 : logical :: small_angle
329 : real(dp) :: bxc_over_m,d1,d2,dvdn,dvdz,fact,m_dot_m1,m_norm
330 : ! d3 and d4 were computed below but never used, and caused test
331 : ! fail problems
332 : !real(dp) :: d3,d4
333 : real(dp) :: dvdn_re,dvdn_im,dvdz_re,dvdz_im
334 : complex(dp) :: rho_updn
335 : real(dp) :: mdirx,mdiry,mdirz,mxy,mx1,my1,mz1,wx,wy,wx1,wy1
336 : real(dp) :: theta0,theta1,theta1_re,theta1_im
337 : real(dp) :: wx1_re,wx1_im
338 : real(dp) :: wy1_re,wy1_im
339 : real(dp) :: mx1_re,mx1_im,my1_re,my1_im,mz1_re,mz1_im
340 : real(dp) :: m_dot_m1_re,m_dot_m1_im
341 : real(dp) :: bxc
342 : !arrays
343 : real(dp) :: vxc_diag(2),v21tmp(2)
344 : complex(dp) :: r1tmp(2,2),u0(2,2),u0_1(2,2),u0_1r1(2,2),u0v1(2,2)
345 : complex(dp) :: rho1_updn(2,2),v1tmp(2,2),vxc1tmp(2,2)
346 : complex(dp) :: rho1_offdiag(2)
347 : ! *************************************************************************
348 :
349 : !DBG_ENTER("COLL")
350 :
351 : !Optional arguments
352 997 : has_mag_norm=present(mag_norm_in)
353 997 : rotation_method=rotation_method_default
354 997 : if (present(rot_method)) rotation_method=rot_method
355 :
356 : !Check Kxc
357 2991 : if (size(kxc)>3*vectsize) then
358 0 : ABI_ERROR('Cannot use Kxc from GGA!')
359 : end if
360 :
361 997 : if((rotation_method==1.or.rotation_method==2).and.cplex==2) then
362 0 : ABI_ERROR('rotation_method=1 and 2 are not available for cplex=2 case! use ixcrot=3')
363 : endif
364 :
365 :
366 : select case (rotation_method)
367 :
368 : !----------------------------------------
369 : ! Taylor expansion of U rotation matrix
370 : !----------------------------------------
371 : case (1)
372 :
373 743583 : do ipt=1,vectsize
374 :
375 743445 : if (has_mag_norm) then
376 743445 : m_norm=mag_norm_in(ipt)
377 : else
378 0 : m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2)
379 : end if
380 :
381 743583 : if(m_norm>m_norm_min) then
382 :
383 : ! Define the U^(0) transformation matrix
384 701352 : rho_updn=(mag(ipt,1)+(zero,one)*mag(ipt,2))
385 701352 : d1=sqrt(( m_norm+mag(ipt,3))**2+abs(rho_updn)**2)
386 701352 : d2=sqrt((-m_norm+mag(ipt,3))**2+abs(rho_updn)**2)
387 : ! d3 and d4 are computed here but never used
388 : !d3=sqrt(( m_norm-mag(ipt,3))**2+abs(rho_updn)**2)
389 : !d4=sqrt(( m_norm+mag(ipt,3))**2-abs(rho_updn)**2)
390 701352 : u0(1,1)=( m_norm+mag(ipt,3))/d1 ! ( m + mz)/d1
391 701352 : u0(2,2)=rho_updn/d2 ! ( mx +imy)/d2
392 701352 : u0(1,2)=(-m_norm+mag(ipt,3))/d2 ! (-m + mz)/d2
393 701352 : u0(2,1)=rho_updn/d1 ! ( mx +imy)/d1
394 :
395 : ! Define the inverse of U^(0): U^(0)^-1
396 701352 : if (abs(rho_updn) > m_norm_min) then
397 674325 : u0_1(1,1)= half*d1/m_norm
398 674325 : u0_1(2,2)= half*d2*(m_norm+mag(ipt,3))/(m_norm*rho_updn)
399 674325 : u0_1(1,2)= half*d1*(m_norm-mag(ipt,3))/(m_norm*rho_updn)
400 674325 : u0_1(2,1)=-half*d2/m_norm
401 : else
402 27027 : u0 = zero
403 27027 : u0(1,1) = one
404 27027 : u0(2,2) = one
405 27027 : u0_1 = zero
406 27027 : u0_1(1,1) = one
407 27027 : u0_1(2,2) = one
408 : end if
409 :
410 : ! Diagonalize the GS Vxc^(0): U^(0)^-1 Vxc^(0) U^(0)
411 : ! (Remember the abinit notation for vxc!)
412 : vxc_diag(1)=half*(vxc(ipt,1)+vxc(ipt,2) &
413 : & -sqrt((vxc(ipt,1)-vxc(ipt,2))**2 &
414 701352 : & +four*(vxc(ipt,3)**2+vxc(ipt,4)**2)))
415 : vxc_diag(2)=half*(vxc(ipt,1)+vxc(ipt,2) &
416 : & +sqrt((vxc(ipt,1)-vxc(ipt,2))**2 &
417 701352 : & +four*(vxc(ipt,3)**2+vxc(ipt,4)**2)))
418 701352 : v1tmp(1,1)=cmplx(real(vxc1_in(ipt,1),kind=dp),zero)
419 701352 : v1tmp(2,2)=cmplx(real(vxc1_in(ipt,2),kind=dp),zero)
420 :
421 : !Transforming the rhor1 with U0
422 701352 : rho1_updn(1,1)=half*(rho1(ipt,1)+rho1(ipt,4))
423 701352 : rho1_updn(2,2)=half*(rho1(ipt,1)-rho1(ipt,4))
424 701352 : rho1_updn(1,2)=half*(rho1(ipt,2)-(zero,one)*rho1(ipt,3))
425 701352 : rho1_updn(2,1)=half*(rho1(ipt,2)+(zero,one)*rho1(ipt,3))
426 10520280 : u0_1r1=matmul(u0_1,rho1_updn)
427 10520280 : r1tmp=matmul(u0_1r1,u0)
428 701352 : rho1_offdiag(1)=r1tmp(1,2) ; rho1_offdiag(2)=r1tmp(2,1)
429 :
430 701352 : if (option==0) then ! for xccc alone
431 449412 : v1tmp(1,2)=cmplx(zero,zero)
432 449412 : v1tmp(2,1)=cmplx(zero,zero)
433 : else
434 251940 : v1tmp(1,2)=-(rho1_offdiag(1)/m_norm)*(vxc_diag(2)-vxc_diag(1))
435 251940 : v1tmp(2,1)= (rho1_offdiag(2)/m_norm)*(vxc_diag(1)-vxc_diag(2))
436 : endif
437 :
438 : !Rotate back the "diagonal" xc computing the term U^(0) Vxc1_^(1) U^(0)^-1
439 10520280 : u0v1=matmul(u0,v1tmp)
440 10520280 : vxc1tmp=matmul(u0v1,u0_1)
441 701352 : vxc1_out(ipt,1)=real(vxc1tmp(1,1),kind=dp)
442 701352 : vxc1_out(ipt,2)=real(vxc1tmp(2,2),kind=dp)
443 701352 : vxc1_out(ipt,3)=real( real(vxc1tmp(1,2)),kind=dp)
444 701352 : vxc1_out(ipt,4)=real(aimag(vxc1tmp(1,2)),kind=dp)
445 :
446 : else ! Magnetization is zero
447 42093 : dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half
448 42093 : mx1=rho1(ipt,2) ; my1=rho1(ipt,3) ; mz1=rho1(ipt,4)
449 : ! Compute Bxc/|m| from Kxc (zero limit)
450 42093 : bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2))
451 42093 : vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1
452 42093 : vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1
453 42093 : vxc1_out(ipt,3)= bxc_over_m*mx1
454 42093 : vxc1_out(ipt,4)=-bxc_over_m*my1
455 : end if
456 :
457 : end do ! ipt
458 :
459 : !----------------------------------------
460 : ! Analytical expression of U rotation matrix
461 : !----------------------------------------
462 : case (2)
463 : !Alternative method (explicitly calculated rotation matrices)
464 : !Vxc^(1) = phixc^(1).Id + // <= change of "electrostatic" XC potential (phixc^(1) is denoted dvdn)
465 : ! + bxc^(1)*( Udag^(0).sigma_z.U^(0) ) + // <= this part describes the change of XC magnetic field magnitude bxc^(1)
466 : ! + bxc^(0)*( Udag^(1).sigma_z.U^(0) + Udag^(0).sigma_z.U^(1) ) // <= remaining terms describe the cost of magnetization rotation
467 :
468 748 : select case(cplex)
469 : case(1)
470 674436 : do ipt=1,vectsize
471 :
472 674325 : if (has_mag_norm) then
473 674325 : m_norm=mag_norm_in(ipt)
474 : else
475 0 : m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2)
476 : end if
477 :
478 :
479 674325 : mx1 =rho1(ipt,2);
480 674325 : my1 =rho1(ipt,3);
481 674325 : mz1 =rho1(ipt,4)
482 :
483 674325 : dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half !phixc^(1)
484 674325 : dvdz=(vxc1_in(ipt,1)-vxc1_in(ipt,2))*half !bxc^(1)
485 :
486 674436 : if (m_norm>m_norm_min) then
487 :
488 674325 : mxy = dsqrt(mag(ipt,1)**2+mag(ipt,2)**2)
489 674325 : small_angle=(mxy/m_norm<tol8) !condition for sin(x)~x to be valid
490 : !even possible to set to tol6
491 674325 : mdirx=mag(ipt,1)/m_norm
492 674325 : mdiry=mag(ipt,2)/m_norm
493 674325 : mdirz=mag(ipt,3)/m_norm
494 :
495 : ! dvdn is phixc^(1) (density only part)
496 : ! dvdz is bxc^(1) (magnetization magnitude part)
497 :
498 : !U^(0)*.Vxc1.U^(0) part
499 674325 : vxc1_out(ipt,1)= dvdn+dvdz*mdirz
500 674325 : vxc1_out(ipt,2)= dvdn-dvdz*mdirz
501 674325 : vxc1_out(ipt,3)= dvdz*mdirx ! Real part
502 674325 : vxc1_out(ipt,4)=-dvdz*mdiry ! Imaginary part, minus sign comes from sigma_y
503 :
504 : !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1) part
505 :
506 : !bxc = dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !bxc^(0)
507 674325 : bxc = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3)*m_norm
508 674325 : if (.not.small_angle) then
509 674325 : wx = mag(ipt,2)/mxy
510 674325 : wy =-mag(ipt,1)/mxy
511 674325 : theta0 = dacos(mag(ipt,3)/m_norm)
512 :
513 674325 : theta1 = (mdirz*(mdirx*mx1+mdiry*my1))/mxy - mz1*mxy/m_norm**2
514 674325 : wx1 = (+mag(ipt,1)**2*my1 - mag(ipt,1)*mag(ipt,2)*mx1)/mxy**2/m_norm ! wx1 multiplied by sin(theta)=mxy/m_norm
515 674325 : wy1 = (-mag(ipt,2)**2*mx1 + mag(ipt,1)*mag(ipt,2)*my1)/mxy**2/m_norm ! wx1 multiplied by sin(theta)=mxy/m_norm
516 :
517 674325 : vxc1_out(ipt,1) = vxc1_out(ipt,1) - bxc*dsin(theta0)*theta1
518 674325 : vxc1_out(ipt,2) = vxc1_out(ipt,2) + bxc*dsin(theta0)*theta1
519 674325 : vxc1_out(ipt,3) = vxc1_out(ipt,3) - bxc*(wy1+dcos(theta0)*wy*theta1)
520 674325 : vxc1_out(ipt,4) = vxc1_out(ipt,4) - bxc*(wx1+dcos(theta0)*wx*theta1)
521 : else
522 : !zero order terms O(1)
523 0 : vxc1_out(ipt,3) = vxc1_out(ipt,3) + bxc*mx1/abs(mag(ipt,3))
524 0 : vxc1_out(ipt,4) = vxc1_out(ipt,4) - bxc*my1/abs(mag(ipt,3))
525 : !first order terms O(theta)
526 0 : fact = bxc/(mag(ipt,3)*abs(mag(ipt,3)))
527 0 : vxc1_out(ipt,1) = vxc1_out(ipt,1) - (mag(ipt,1)*mx1+mag(ipt,2)*my1)*fact
528 0 : vxc1_out(ipt,2) = vxc1_out(ipt,2) + (mag(ipt,1)*mx1+mag(ipt,2)*my1)*fact
529 0 : vxc1_out(ipt,3) = vxc1_out(ipt,3) - mag(ipt,1)*mz1*fact
530 0 : vxc1_out(ipt,4) = vxc1_out(ipt,4) + mag(ipt,2)*mz1*fact
531 : endif
532 :
533 : else ! Magnetization is zero
534 : ! Compute Bxc/|m| from Kxc (zero limit)
535 0 : bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2))
536 0 : vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1
537 0 : vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1
538 0 : vxc1_out(ipt,3)= bxc_over_m*mx1
539 0 : vxc1_out(ipt,4)=-bxc_over_m*my1
540 : end if
541 : end do ! ipt
542 :
543 : case(2) !cplex=2
544 :
545 111 : do ipt=1,vectsize
546 :
547 0 : if (has_mag_norm) then
548 0 : m_norm=mag_norm_in(ipt)
549 : else
550 0 : m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2)
551 : end if
552 :
553 0 : mx1_re= rho1(2*ipt-1,2); mx1_im= rho1(2*ipt,2)
554 0 : my1_re= rho1(2*ipt-1,3); my1_im= rho1(2*ipt,3)
555 0 : mz1_re= rho1(2*ipt-1,4); mz1_im= rho1(2*ipt,4)
556 :
557 0 : dvdn_re=(vxc1_in(2*ipt-1,1)+vxc1_in(2*ipt-1,2))*half
558 0 : dvdz_re=(vxc1_in(2*ipt-1,1)-vxc1_in(2*ipt-1,2))*half
559 0 : dvdn_im=(vxc1_in(2*ipt ,1)+vxc1_in(2*ipt ,2))*half
560 0 : dvdz_im=(vxc1_in(2*ipt ,1)-vxc1_in(2*ipt ,2))*half
561 :
562 0 : if (m_norm>m_norm_min) then
563 :
564 0 : mdirx=mag(ipt,1)/m_norm
565 0 : mdiry=mag(ipt,2)/m_norm
566 0 : mdirz=mag(ipt,3)/m_norm
567 :
568 0 : mxy = dsqrt(mag(ipt,1)**2+mag(ipt,2)**2)
569 0 : small_angle=(mxy/m_norm<tol8) !condition for sin(x)~x to be valid
570 : !
571 0 : mdirx=mag(ipt,1)/m_norm
572 0 : mdiry=mag(ipt,2)/m_norm
573 0 : mdirz=mag(ipt,3)/m_norm
574 :
575 : ! dvdn is phixc^(1) (density only part)
576 : ! dvdz is bxc^(1) (magnetization magnitude part)
577 :
578 : !U^(0)*.Vxc1.U^(0) part
579 0 : vxc1_out(2*ipt-1,1)= dvdn_re+dvdz_re*mdirz
580 0 : vxc1_out(2*ipt ,1)= dvdn_im+dvdz_im*mdirz
581 0 : vxc1_out(2*ipt-1,2)= dvdn_re-dvdz_re*mdirz
582 0 : vxc1_out(2*ipt ,2)= dvdn_im-dvdz_im*mdirz
583 : !NOTE: change of definition of the potential matrix components
584 : ! vxc1_out(:,3) = V_updn
585 : ! vxc1_out(:,4) = i.V_updn
586 0 : vxc1_out(2*ipt-1,3)= dvdz_re*mdirx + dvdz_im*mdiry !Re[ V^12]
587 0 : vxc1_out(2*ipt ,3)= dvdz_im*mdirx - dvdz_re*mdiry !Im[ V^12]
588 :
589 : !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1) part
590 : !bxc = dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !bxc^(0)
591 0 : bxc = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3)*m_norm
592 0 : if (.not.small_angle) then
593 0 : wx = mag(ipt,2)/mxy
594 0 : wy =-mag(ipt,1)/mxy
595 0 : theta0 = dacos(mag(ipt,3)/m_norm)
596 :
597 :
598 0 : theta1_re = (mdirz*(mdirx*mx1_re+mdiry*my1_re))/mxy - mz1_re*mxy/m_norm**2
599 0 : theta1_im = (mdirz*(mdirx*mx1_im+mdiry*my1_im))/mxy - mz1_im*mxy/m_norm**2
600 :
601 0 : wx1_re = (+mag(ipt,1)**2*my1_re - mag(ipt,1)*mag(ipt,2)*mx1_re)/mxy**2/m_norm ! wx1 multiplied by sin(theta)=mxy/m_norm
602 0 : wx1_im = (+mag(ipt,1)**2*my1_im - mag(ipt,1)*mag(ipt,2)*mx1_im)/mxy**2/m_norm
603 0 : wy1_re = (-mag(ipt,2)**2*mx1_re + mag(ipt,1)*mag(ipt,2)*my1_re)/mxy**2/m_norm ! wy1 multiplied by sin(theta)=mxy/m_norm
604 0 : wy1_im = (-mag(ipt,2)**2*mx1_im + mag(ipt,1)*mag(ipt,2)*my1_im)/mxy**2/m_norm
605 :
606 : !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1)
607 0 : vxc1_out(2*ipt-1,1) = vxc1_out(2*ipt-1,1) - bxc*dsin(theta0)*theta1_re
608 0 : vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) - bxc*dsin(theta0)*theta1_im
609 0 : vxc1_out(2*ipt-1,2) = vxc1_out(2*ipt-1,2) + bxc*dsin(theta0)*theta1_re
610 0 : vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + bxc*dsin(theta0)*theta1_im
611 : !cplex=1 part:
612 : !v12 += -(bxc)*(wy1+dcos(theta0)*wy*theta1)-
613 : ! -i.(bxc)*(wx1+dcos(theta0)*wx*theta1)
614 0 : vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - bxc*(wy1_re+dcos(theta0)*wy*theta1_re)
615 0 : vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*(wx1_im+dcos(theta0)*wx*theta1_im)
616 0 : vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*(wy1_im+dcos(theta0)*wy*theta1_im)
617 0 : vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*(wx1_re+dcos(theta0)*wx*theta1_re)
618 : else
619 : !small theta case:
620 : !zero order terms O(1)
621 0 : vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*mx1_re/abs(mag(ipt,3))
622 0 : vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*my1_im/abs(mag(ipt,3))
623 0 : vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc*mx1_im/abs(mag(ipt,3))
624 0 : vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*my1_re/abs(mag(ipt,3))
625 : !first order terms:
626 0 : fact = bxc/(mag(ipt,3)*abs(mag(ipt,3)))
627 0 : vxc1_out(2*ipt-1,1) = vxc1_out(2*ipt-1,1) - (mag(ipt,1)*mx1_re+mag(ipt,2)*my1_re)*fact
628 0 : vxc1_out(2*ipt-1,2) = vxc1_out(2*ipt-1,2) + (mag(ipt,1)*mx1_re+mag(ipt,2)*my1_re)*fact
629 0 : vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) - (mag(ipt,1)*mx1_im+mag(ipt,2)*my1_im)*fact
630 0 : vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + (mag(ipt,1)*mx1_im+mag(ipt,2)*my1_im)*fact
631 :
632 0 : vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - mag(ipt,1)*mz1_re*fact
633 0 : vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - mag(ipt,2)*mz1_im*fact
634 0 : vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + mag(ipt,2)*mz1_re*fact
635 0 : vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - mag(ipt,1)*mz1_im*fact
636 : endif
637 :
638 : else ! Magnetization is practically zero
639 : ! Compute Bxc/|m| from Kxc (zero limit)
640 0 : bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2))
641 0 : vxc1_out(2*ipt-1,1)= dvdn_re + bxc_over_m*mz1_re
642 0 : vxc1_out(2*ipt ,1)= dvdn_im + bxc_over_m*mz1_im
643 0 : vxc1_out(2*ipt-1,2)= dvdn_re - bxc_over_m*mz1_re
644 0 : vxc1_out(2*ipt ,2)= dvdn_im - bxc_over_m*mz1_im
645 0 : vxc1_out(2*ipt-1,3)= bxc_over_m*( mx1_re+my1_im)
646 0 : vxc1_out(2*ipt ,3)= bxc_over_m*(-my1_re+mx1_im)
647 : end if
648 : !finally reconstruct i.V^12 from V^12
649 0 : vxc1_out(2*ipt-1,4) = vxc1_out(2*ipt ,3) ! Re[i.V^21] = Im[V^12]
650 0 : vxc1_out(2*ipt ,4) = vxc1_out(2*ipt-1,3) ! Im[i.V^21] = Re[V^12]
651 :
652 : end do ! ipt
653 :
654 : end select
655 :
656 : !----------------------------------------
657 : ! Explicit derivative of the rotated XC functional
658 : !----------------------------------------
659 : case (3)
660 : ! Brute-force derivative of Vxc
661 : ! Explicit calculation of the rotated xc functional
662 : ! (derivatives of the analytical expression) (Eq. A)
663 : ! Vxc^(1) = phixc^(1).Id + // <= change of "electrostatic" XC potential (phixc^(1) is denoted dvdn)
664 : ! + bxc^(1)*(sigma,m^(0))/|m^(0)| + // <= this term is equivalent to ( Udag^(0).sigma_z.U^(0) ) term in rotation_method=2
665 : ! + bxc^(0)*(sigma,m^(1)))/|m^(0)| - // <= the last terms are equivalent to ( Udag^(1).sigma_z.U^(0) + Udag^(0).sigma_z.U^(1) )
666 : ! - bxc^(0)*(sigma,m^(0))*(m^(1),m^(0))/|m^(0)|**3
667 997 : select case(cplex)
668 : case(1)
669 :
670 992000 : do ipt=1,vectsize
671 :
672 991765 : if (has_mag_norm) then
673 991765 : m_norm=mag_norm_in(ipt)
674 : else
675 0 : m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2)
676 : end if
677 :
678 : ! dvdn is phixc^(1) (density only part)
679 : ! dvdz is bxc^(1) (magnetization magnitude part)
680 991765 : dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half
681 991765 : dvdz=(vxc1_in(ipt,1)-vxc1_in(ipt,2))*half
682 :
683 991765 : mx1=rho1(ipt,2) ; my1=rho1(ipt,3) ; mz1=rho1(ipt,4)
684 :
685 992000 : if(m_norm>m_norm_min) then
686 :
687 798449 : mdirx=mag(ipt,1)/m_norm; mdiry=mag(ipt,2)/m_norm; mdirz=mag(ipt,3)/m_norm
688 :
689 : !This part describes the change of the magnitude of the xc magnetic field
690 : !and the change of the scalar part of the xc electrostatic potential, 1st + 2nd term in Eq.A
691 : !phixc^(1).Id + bxc^(1) (sigma,m^(0))/|m^(0)|
692 798449 : vxc1_out(ipt,1)= dvdn+dvdz*mdirz
693 798449 : vxc1_out(ipt,2)= dvdn-dvdz*mdirz
694 798449 : vxc1_out(ipt,3)= dvdz*mdirx ! Real part
695 798449 : vxc1_out(ipt,4)=-dvdz*mdiry ! Imaginary part, minus sign comes from sigma_y
696 :
697 798449 : if (option/=0) then
698 : !Add remaining contributions coming from the change of magnetization direction
699 : !projection of m^(1) on gs magnetization direction
700 255944 : m_dot_m1=(mdirx*rho1(ipt,2)+mdiry*rho1(ipt,3)+mdirz*rho1(ipt,4))
701 :
702 255944 : bxc_over_m =-dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !this is bxc^(0)
703 255944 : bxc_over_m = bxc_over_m/m_norm
704 255944 : vxc1_out(ipt,1) = vxc1_out(ipt,1) + bxc_over_m*( mz1 - mdirz*m_dot_m1 ) !
705 255944 : vxc1_out(ipt,2) = vxc1_out(ipt,2) + bxc_over_m*(-mz1 + mdirz*m_dot_m1 ) !
706 255944 : vxc1_out(ipt,3) = vxc1_out(ipt,3) + bxc_over_m*( mx1 - mdirx*m_dot_m1 ) !
707 255944 : vxc1_out(ipt,4) = vxc1_out(ipt,4) + bxc_over_m*(-my1 + mdiry*m_dot_m1 ) !
708 : endif
709 :
710 : else
711 193316 : if (option/=0) then
712 : !Compute bxc^(0)/|m| from kxc (|m^(0)| -> zero limit)
713 29621 : bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2))
714 29621 : vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1
715 29621 : vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1
716 29621 : vxc1_out(ipt,3)= bxc_over_m*mx1
717 29621 : vxc1_out(ipt,4)=-bxc_over_m*my1
718 : else
719 163695 : vxc1_out(ipt,1)= dvdn
720 163695 : vxc1_out(ipt,2)= dvdn
721 163695 : vxc1_out(ipt,3)= zero
722 163695 : vxc1_out(ipt,4)= zero
723 : endif
724 : end if
725 :
726 : end do ! ipt
727 :
728 : case(2)
729 : !cplex=2 case
730 :
731 2345708 : do ipt=1,vectsize
732 :
733 2344960 : if (has_mag_norm) then
734 2344960 : m_norm=mag_norm_in(ipt)
735 : else
736 0 : m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2)
737 : end if
738 :
739 : ! see cplex=1 case for details
740 2344960 : dvdn_re=(vxc1_in(2*ipt-1,1)+vxc1_in(2*ipt-1,2))*half
741 2344960 : dvdn_im=(vxc1_in(2*ipt ,1)+vxc1_in(2*ipt ,2))*half
742 2344960 : dvdz_re=(vxc1_in(2*ipt-1,1)-vxc1_in(2*ipt-1,2))*half
743 2344960 : dvdz_im=(vxc1_in(2*ipt ,1)-vxc1_in(2*ipt ,2))*half
744 :
745 2344960 : mx1_re=rho1(2*ipt-1,2); mx1_im=rho1(2*ipt,2)
746 2344960 : my1_re=rho1(2*ipt-1,3); my1_im=rho1(2*ipt,3)
747 2344960 : mz1_re=rho1(2*ipt-1,4); mz1_im=rho1(2*ipt,4)
748 :
749 2344960 : if(m_norm>m_norm_min) then
750 :
751 1684571 : mdirx=mag(ipt,1)/m_norm; mdiry=mag(ipt,2)/m_norm; mdirz=mag(ipt,3)/m_norm
752 :
753 : !first two terms:
754 1684571 : vxc1_out(2*ipt-1,1)= dvdn_re+dvdz_re*mdirz
755 1684571 : vxc1_out(2*ipt ,1)= dvdn_im+dvdz_im*mdirz
756 1684571 : vxc1_out(2*ipt-1,2)= dvdn_re-dvdz_re*mdirz
757 1684571 : vxc1_out(2*ipt ,2)= dvdn_im-dvdz_im*mdirz
758 : !NOTE: change of definition of the potential matrix components
759 : ! vxc1_out(:,3) = V_updn
760 : ! vxc1_out(:,4) = i.V_dnup
761 :
762 : ! V^12 = dvdz*mx/|m| - i.dvdz*my/|m| = (Re[dvdz]*mx/|m| + Im[dvdz]*my/|m|) + i.(Im[dvdz]*mx/|m| - Re[dvdz]*my/|m|) => vxc1(:,3)
763 : ! V^21 = dvdz*mx/|m| + i.dvdz*my/|m| = (Re[dvdz]*mx/|m| - Im[dvdz]*my/|m|) + i.(Im[dvdz]*mx/|m| + Re[dvdz]*my/|m|)
764 1684571 : vxc1_out(2*ipt-1,3)= dvdz_re*mdirx + dvdz_im*mdiry !Re[V^12]
765 1684571 : vxc1_out(2*ipt ,3)= dvdz_im*mdirx - dvdz_re*mdiry !Im[V^12]
766 1684571 : vxc1_out(2*ipt-1,4)= dvdz_re*mdirx - dvdz_im*mdiry !Re[V^21]
767 1684571 : vxc1_out(2*ipt ,4)= dvdz_im*mdirx + dvdz_re*mdiry !Im[V^21]
768 1684571 : if (option/=0) then
769 :
770 : !remaining contributions:
771 571295 : m_dot_m1_re= mdirx*mx1_re + mdiry*my1_re + mdirz*mz1_re
772 571295 : m_dot_m1_im= mdirx*mx1_im + mdiry*my1_im + mdirz*mz1_im
773 :
774 571295 : bxc_over_m =-dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !this is bxc^(0)
775 571295 : bxc_over_m = bxc_over_m/m_norm
776 : !bxc_over_m = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3)
777 :
778 571295 : vxc1_out(2*ipt-1,1) = vxc1_out(2*ipt-1,1) + bxc_over_m*( mz1_re - mdirz*m_dot_m1_re ) ! Re[V^11]
779 571295 : vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) + bxc_over_m*( mz1_im - mdirz*m_dot_m1_im ) ! Im[V^11]
780 571295 : vxc1_out(2*ipt-1,2) = vxc1_out(2*ipt-1,2) + bxc_over_m*(-mz1_re + mdirz*m_dot_m1_re ) ! Re[V^22]
781 571295 : vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + bxc_over_m*(-mz1_im + mdirz*m_dot_m1_im ) ! Im[V^22]
782 :
783 : ! v12 += bxc_over_m*( (mx1 - mdirx*m_dot_m1 ) - i.( my1 - mdiry*m_dot_m1 ) ) <= see cplex=1
784 : ! Re[v12] += bxc_over_m*( (mx1_re - mdirx*m_dot_m1_re) + ( my1_im - mdiry*m_dot_m1_im) )
785 : ! Im[v12] += bxc_over_m*( (mx1_im - mdirx*m_dot_m1_im) + (-my1_re + mdiry*m_dot_m1_re) )
786 571295 : vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc_over_m*( mx1_re - mdirx*m_dot_m1_re ) ! Re[V^12]
787 571295 : vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc_over_m*( my1_im - mdiry*m_dot_m1_im ) ! Re[V^12]
788 571295 : vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc_over_m*( mx1_im - mdirx*m_dot_m1_im ) ! Im[V^12]
789 571295 : vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc_over_m*(-my1_re + mdiry*m_dot_m1_re ) ! Im[V^12]
790 :
791 : ! v21 += bxc_over_m*( (mx1 - mdirx*m_dot_m1 ) + i.( my1 - mdiry*m_dot_m1 ) )
792 : ! Re[v21] += bxc_over_m*( (mx1_re - mdirx*m_dot_m1_re) + (-my1_im + mdiry*m_dot_m1_im) )
793 : ! Im[v21] += bxc_over_m*( (mx1_im - mdirx*m_dot_m1_im) + ( my1_re - mdiry*m_dot_m1_re) )
794 : ! the 4th component is actually not v21, but rather i.v21, this will be adjusted later
795 571295 : vxc1_out(2*ipt-1,4) = vxc1_out(2*ipt-1,4) + bxc_over_m*( mx1_re - mdirx*m_dot_m1_re ) ! Re[V^21]
796 571295 : vxc1_out(2*ipt-1,4) = vxc1_out(2*ipt-1,4) + bxc_over_m*(-my1_im + mdiry*m_dot_m1_im ) ! Re[V^21]
797 571295 : vxc1_out(2*ipt ,4) = vxc1_out(2*ipt ,4) + bxc_over_m*( mx1_im - mdirx*m_dot_m1_im ) ! Im[V^21]
798 571295 : vxc1_out(2*ipt ,4) = vxc1_out(2*ipt ,4) + bxc_over_m*( my1_re - mdiry*m_dot_m1_re ) ! Im[V^21]
799 : endif
800 : else
801 660389 : if(option/=0) then
802 : !Compute Bxc/|m| from Kxc (|m^(0)| -> zero limit)
803 84065 : bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2))
804 84065 : vxc1_out(2*ipt-1,1)= dvdn_re + bxc_over_m*mz1_re
805 84065 : vxc1_out(2*ipt-1,2)= dvdn_re - bxc_over_m*mz1_re
806 84065 : vxc1_out(2*ipt ,1)= dvdn_im + bxc_over_m*mz1_im
807 84065 : vxc1_out(2*ipt ,2)= dvdn_im - bxc_over_m*mz1_im
808 :
809 84065 : vxc1_out(2*ipt-1,3)= bxc_over_m*(mx1_re+my1_im)
810 84065 : vxc1_out(2*ipt ,3)= bxc_over_m*(mx1_im-my1_re)
811 84065 : vxc1_out(2*ipt-1,4)= bxc_over_m*(mx1_re-my1_im)
812 84065 : vxc1_out(2*ipt ,4)= bxc_over_m*(mx1_im+my1_re)
813 : else
814 576324 : vxc1_out(2*ipt-1,1)= dvdn_re
815 576324 : vxc1_out(2*ipt-1,2)= dvdn_re
816 576324 : vxc1_out(2*ipt ,1)= dvdn_im
817 576324 : vxc1_out(2*ipt ,2)= dvdn_im
818 :
819 576324 : vxc1_out(2*ipt-1,3)= zero
820 576324 : vxc1_out(2*ipt ,3)= zero
821 576324 : vxc1_out(2*ipt-1,4)= zero
822 576324 : vxc1_out(2*ipt ,4)= zero
823 : endif
824 : end if
825 :
826 : !finally reconstruct i.V^21 from V^21
827 2344960 : v21tmp(1) = vxc1_out(2*ipt-1,4) !Re[V^21]
828 2344960 : v21tmp(2) = vxc1_out(2*ipt ,4) !Im[V^21]
829 :
830 2344960 : vxc1_out(2*ipt-1,4) =-v21tmp(2) ! Re[i.V^21]=-Im[V^21]
831 2345473 : vxc1_out(2*ipt ,4) = v21tmp(1) ! Im[i.V^21]= Re[V^21]
832 :
833 : end do ! ipt
834 :
835 : end select !cplex
836 :
837 : end select ! rotation_method
838 :
839 : !DBG_EXIT("COLL")
840 :
841 997 : end subroutine rotate_back_mag_dfpt
842 : !!***
843 :
844 :
845 : !!****f* ABINIT/m_xc_noncoll/test_rotations
846 : !! NAME
847 : !! test_rotations
848 : !!
849 : !! FUNCTION
850 : !! Test three different methods in rotate_back_mag_dfpt
851 : !!
852 : !! INPUTS
853 : !! option= types of tests to perform
854 : !! 0=> only quick tests
855 : !! 1=> quick and slow tests
856 : !! cplex = complex or real potential and first order magnetization
857 : !!
858 : !! OUTPUT
859 : !!
860 : !! SIDE EFFECTS
861 : !!
862 : !! NOTES
863 : !!
864 : !! For debug purposes
865 : !!
866 : !! SOURCE
867 :
868 0 : subroutine test_rotations(option,cplex)
869 :
870 : !Arguments ------------------------------------
871 : integer , intent(in) :: option
872 : integer , intent(in) :: cplex
873 :
874 : !Local variables-------------------------------
875 : real(dp) :: m0(1,3),vxc0(1,4),kxc(1,3)
876 0 : real(dp) :: n1(cplex,4),vxc1_in(cplex,4),vxc1_out(cplex,4)
877 0 : real(dp) :: delta_23(cplex,4) !,delta_12(cplex,4)
878 : real(dp) :: m0_norm,dvdn,dvdz,err23 !,wrong_comp!,err12
879 : real(dp) :: theta0,phi0,theta1,phi1,err,m1_norm
880 : integer :: dir0,dir1
881 : ! *************************************************************************
882 :
883 : DBG_ENTER("COLL")
884 :
885 : ! if (option/=1 .and. option/=2 ) then
886 : ! write(msg,'(3a,i0)')&
887 : !& 'The argument option should be 1 or 2,',ch10,&
888 : !& 'however, option=',option
889 : ! ABI_BUG(msg)
890 : ! end if
891 : !
892 : ! if (sizein<1) then
893 : ! write(msg,'(3a,i0)')&
894 : !& ' The argument sizein should be a positive number,',ch10,&
895 : !& ' however, sizein=',sizein
896 : ! ABI_ERROR(msg)
897 : ! end if
898 :
899 : DBG_EXIT("COLL")
900 :
901 : !write(*,*) 'VXC_NONCOLL TESTS================================================================'
902 0 : if (cplex==1) then
903 : !write(*,*) ' cplex=1------------------------------------------------------------------------'
904 :
905 : !write(*,*) ' TEST: simple m* orietnations, bxc^(1) part'
906 : dvdn=zero;dvdz=1.0!
907 : err23=zero
908 0 : do dir0=1,3
909 0 : m0=zero; n1=zero
910 0 : do dir1=2,4
911 0 : m0(1,dir0)=0.1
912 0 : m0_norm=sqrt(m0(1,1)**2+m0(1,2)**2+m0(1,3)**2)
913 0 : n1(1,dir1)=0.8 ! any number would do here
914 :
915 0 : vxc0=zero; ! no bxc^(0) part at all
916 :
917 0 : vxc1_in(1,1)= dvdn+dvdz
918 0 : vxc1_in(1,2)= dvdn-dvdz
919 :
920 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1)
921 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2)
922 0 : delta_23=vxc1_out
923 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3)
924 0 : delta_23=abs(delta_23-vxc1_out)
925 0 : err=max(delta_23(1,1),delta_23(1,2),delta_23(1,3),delta_23(1,4))
926 0 : if (err23<err) err23=err;
927 :
928 : enddo
929 : enddo
930 : !write(*,*) ' maximum mismatch between methods 2 and 3:',err23
931 :
932 : !write(*,*) ' TEST: simple m* orietnations, bxc^(0) part'
933 :
934 : err23=zero
935 : dvdn=zero;dvdz=1.0
936 0 : do dir0=1,3
937 0 : m0=zero; n1=zero
938 0 : do dir1=2,4
939 0 : m0(1,dir0)=0.1
940 0 : m0_norm=sqrt(m0(1,1)**2+m0(1,2)**2+m0(1,3)**2)
941 0 : n1(1,dir1)=0.8 ! =m^1, any number would do here
942 :
943 0 : vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm
944 0 : vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm
945 0 : vxc0(1,3) = dvdz*m0(1,1)/m0_norm
946 0 : vxc0(1,4) =-dvdz*m0(1,2)/m0_norm
947 :
948 0 : vxc1_in=zero !vxc^(1) collinear is zero
949 :
950 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1)
951 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2)
952 0 : delta_23=vxc1_out
953 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3)
954 0 : delta_23=abs(delta_23-vxc1_out)
955 0 : err=maxval(abs(delta_23(1,:)))
956 0 : if (err23<err) err23=err;
957 : enddo
958 : enddo
959 : !write(*,*) ' maximum mismatch between methods 2 and 3:',err23
960 :
961 : !write(*,*) ' TEST: general m0 orietnations, bxc^(0) part'
962 :
963 : theta0=zero
964 : err23=zero
965 0 : m0_norm=0.3
966 0 : do while(theta0<=pi)
967 : phi0=zero
968 0 : do while(phi0<=2*pi)
969 0 : m0(1,1)=m0_norm*sin(theta0)*cos(phi0)
970 0 : m0(1,2)=m0_norm*sin(theta0)*sin(phi0)
971 0 : m0(1,3)=m0_norm*cos(theta0)
972 :
973 0 : do dir1=2,4
974 0 : n1=zero
975 0 : n1(1,dir1)=0.8 ! =m^1, any number would do here
976 :
977 : !vxc0=zero; !
978 0 : vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm
979 0 : vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm
980 0 : vxc0(1,3) = dvdz*m0(1,1)/m0_norm
981 0 : vxc0(1,4) =-dvdz*m0(1,2)/m0_norm
982 :
983 0 : vxc1_in=zero
984 :
985 : !call rotate_back_mag_dfpt(vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1)
986 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2)
987 0 : delta_23=vxc1_out
988 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3)
989 0 : delta_23=abs(delta_23-vxc1_out)
990 0 : err=maxval(abs(delta_23(1,:)))
991 0 : if (err23<err) err23=err;
992 : enddo
993 0 : phi0=phi0+2*pi/100.0
994 : enddo
995 0 : theta0=theta0+pi/100.0
996 : enddo
997 : !write(*,*) ' maximum mismatch between methods 2 and 3:',err23
998 :
999 0 : if(option==2) then
1000 : !write(*,*) ' TEST: general m* orietnations, bxc^(0) part'
1001 : dvdn=zero;dvdz=1.0
1002 :
1003 : theta0=zero
1004 : err23=zero
1005 0 : m0_norm=0.3
1006 : m1_norm=10.5
1007 0 : do while(theta0<=pi) !loops on orientation of m^(0)
1008 : phi0=zero
1009 0 : do while(phi0<=2*pi)
1010 0 : m0(1,1)=m0_norm*sin(theta0)*cos(phi0)
1011 0 : m0(1,2)=m0_norm*sin(theta0)*sin(phi0)
1012 0 : m0(1,3)=m0_norm*cos(theta0)
1013 :
1014 0 : vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm
1015 0 : vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm
1016 0 : vxc0(1,3) = dvdz*m0(1,1)/m0_norm
1017 0 : vxc0(1,4) =-dvdz*m0(1,2)/m0_norm
1018 :
1019 0 : vxc1_in=zero
1020 :
1021 : theta1=zero
1022 0 : do while(theta1<=pi) !loops on orientation of m^(1)
1023 : phi1=zero
1024 0 : do while(phi1<=2*pi)
1025 0 : n1(1,1)=zero
1026 0 : n1(1,2)=m1_norm*sin(theta1)*cos(phi1)
1027 0 : n1(1,3)=m1_norm*sin(theta1)*sin(phi1)
1028 0 : n1(1,4)=m1_norm*cos(theta1)
1029 :
1030 : !vxc0=zero; !
1031 : !call rotate_back_mag_dfpt(vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1)
1032 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2)
1033 0 : delta_23=vxc1_out
1034 0 : call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3)
1035 0 : delta_23=abs(delta_23-vxc1_out)
1036 0 : err=maxval(abs(delta_23(1,:)))
1037 : if (err23<err) err23=err;
1038 0 : phi1=phi1+2*pi/100.0
1039 : enddo
1040 0 : theta1=theta1+pi/100.0
1041 : enddo
1042 :
1043 0 : phi0=phi0+2*pi/100.0
1044 : enddo
1045 0 : theta0=theta0+pi/100.0
1046 : enddo
1047 : !write(*,*) ' maximum mismatch between methods 2 and 3:',err23
1048 : endif
1049 :
1050 : !else !cplex=2
1051 :
1052 : endif
1053 :
1054 0 : end subroutine test_rotations
1055 : !!***
1056 :
1057 :
1058 : !----------------------------------------------------------------------
1059 :
1060 : END MODULE m_xc_noncoll
1061 : !!***
|
