Branch data 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 : : !!***
|
