Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2022 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_get_int_perturbation
8 : USE m_juDFT
9 : USE m_fft3d
10 : USE m_types
11 :
12 : IMPLICIT NONE
13 :
14 : CONTAINS
15 :
16 0 : SUBROUTINE get_int_local_perturbation(sym, stars, atoms, sphhar, &
17 : & input, den, den1, den1im, starsq)
18 :
19 : USE m_constants
20 :
21 : TYPE(t_input), INTENT(IN) :: input
22 : TYPE(t_sym), INTENT(IN) :: sym
23 : TYPE(t_stars), INTENT(IN) :: stars, starsq
24 : TYPE(t_sphhar), INTENT(IN) :: sphhar
25 : TYPE(t_atoms), INTENT(IN) :: atoms
26 : TYPE(t_potden), INTENT(IN) :: den
27 : TYPE(t_potden), INTENT(INOUT) :: den1, den1im
28 :
29 : INTEGER :: iden, jspin, ifft3
30 : INTEGER :: ityp, iri, imesh
31 : REAL :: rho_11, rho_22, m
32 : REAL :: rhotot, rho_up, rho_down, theta, phi
33 : COMPLEX :: m1, mx1, my1, mz1, n1, t1, p1, rho1_up, rho1_down
34 :
35 0 : REAL, ALLOCATABLE :: ris(:,:), ris_real(:,:), ris_imag(:,:)
36 : REAL, ALLOCATABLE :: fftwork(:)
37 :
38 0 : ifft3 = 27*stars%mx1*stars%mx2*stars%mx3
39 :
40 : !TODO: Make sure the indices for rho1 are 1,2,3,4 == n1,mx1,my1,mz1
41 0 : ALLOCATE (ris(ifft3,2),fftwork(ifft3))
42 0 : ALLOCATE (ris_real(ifft3,4),ris_imag(ifft3,4))
43 :
44 0 : ALLOCATE(den1%phi_pw(ifft3),den1%theta_pw(ifft3))
45 0 : ALLOCATE(den1im%phi_pw(ifft3),den1im%theta_pw(ifft3))
46 :
47 0 : DO iden = 1, 2
48 0 : CALL fft3d(ris(:,iden), fftwork, den%pw(:,iden), stars, +1)
49 0 : CALL fft3d(ris_real(:,iden), ris_imag(:,iden), den1%pw(:,iden), starsq, +1)
50 : END DO
51 :
52 0 : CALL fft3d(ris_real(:,3), ris_imag(:,3), den1%pw(:,3), starsq, +1)
53 0 : CALL fft3d(ris_real(:,4), ris_imag(:,4), den1%pw(:,4), starsq, +1)
54 :
55 0 : DO imesh = 1, ifft3
56 : ! Get real space density matrix elements
57 0 : rho_11 = ris(imesh,1)
58 0 : rho_22 = ris(imesh,2)
59 :
60 : ! Calculate unperturbed magnetization density
61 0 : m = rho_11 - rho_22
62 :
63 : ! Calculate perturbed total and magnetization density
64 0 : n1 = ris_real(imesh,1) + Imagunit * ris_imag(imesh,1)
65 0 : mx1 = ris_real(imesh,2) + Imagunit * ris_imag(imesh,2)
66 0 : my1 = ris_real(imesh,3) + Imagunit * ris_imag(imesh,3)
67 0 : mz1 = ris_real(imesh,4) + Imagunit * ris_imag(imesh,4)
68 :
69 0 : theta = den%theta_pw(imesh)
70 0 : phi = den%phi_pw(imesh)
71 :
72 : ! Calculate the perturbed absolute value of the magnetization density
73 0 : m1 = cmplx(0.0, 0.0)
74 0 : m1 = m1 + mx1 * SIN(theta) * COS(phi)
75 0 : m1 = m1 + my1 * SIN(theta) * SIN(phi)
76 0 : m1 = m1 + mz1 * COS(theta)
77 :
78 : ! Calculate the perturbed angles
79 0 : t1 = cmplx(0.0, 0.0)
80 0 : t1 = t1 + mx1 * COS(theta) * COS(phi) / m
81 0 : t1 = t1 + my1 * COS(theta) * SIN(phi) / m
82 0 : t1 = t1 - mz1 * SIN(theta) / m
83 :
84 0 : p1 = cmplx(0.0, 0.0)
85 0 : p1 = p1 - mx1 * SIN(theta) * SIN(phi) / m
86 0 : p1 = p1 + my1 * SIN(theta) * COS(phi) / m
87 :
88 0 : rho1_up = (n1 + m1)/2
89 0 : rho1_down = (n1 - m1)/2
90 :
91 0 : ris_real(imesh,1) = REAL(rho1_up )
92 0 : ris_real(imesh,2) = REAL(rho1_down)
93 0 : ris_imag(imesh,1) = AIMAG(rho1_up )
94 0 : ris_imag(imesh,2) = AIMAG(rho1_down)
95 0 : den1%theta_pw(imesh) = REAL(t1)
96 0 : den1%phi_pw(imesh) = REAL(p1)
97 0 : den1im%theta_pw(imesh) = AIMAG(t1)
98 0 : den1im%phi_pw(imesh) = AIMAG(p1)
99 : END DO
100 :
101 : ! Fourier tranform the up- and down density perturbations back to reciprocal space:
102 0 : DO jspin = 1, input%jspins
103 0 : CALL fft3d(ris_real(:,jspin),ris_imag(:,jspin),den1%pw(:,jspin),starsq,-1)
104 : END DO
105 0 : END SUBROUTINE get_int_local_perturbation
106 :
107 0 : SUBROUTINE get_int_global_perturbation(stars,atoms,sym,input,den,den1,den1im,vTot,vTot1,starsq)
108 : TYPE(t_input), INTENT(IN) :: input
109 : TYPE(t_sym), INTENT(IN) :: sym
110 : TYPE(t_stars), INTENT(IN) :: stars, starsq
111 : TYPE(t_atoms), INTENT(IN) :: atoms
112 : TYPE(t_potden), INTENT(IN) :: den, den1, den1im, vTot
113 : TYPE(t_potden), INTENT(INOUT) :: vTot1
114 :
115 : INTEGER :: imeshpt, ipot, jspin
116 : INTEGER :: ifft3, i
117 : REAL :: theta, phi, v11, v22
118 :
119 0 : REAL, ALLOCATABLE :: vis(:,:), vis_re(:,:), vis_im(:,:), vis2_re(:,:), vis2_im(:,:), v1re(:,:), v1im(:,:)
120 : REAL, ALLOCATABLE :: fftwork(:)
121 :
122 : COMPLEX :: a11, a22, a21, a12, av11, av22, av21, av12
123 : COMPLEX :: t1, p1, v21, v12, v1up, v1down, v1, b1, v1mat11, v1mat22, v1mat21, v1mat12
124 :
125 0 : ifft3 = 27*stars%mx1*stars%mx2*stars%mx3 !TODO: What if starsq/=stars in that regard?
126 0 : IF (ifft3.NE.SIZE(den%theta_pw)) CALL judft_error("Wrong size of angles")
127 :
128 0 : ALLOCATE (vis(ifft3,4),vis_re(ifft3,4),vis_im(ifft3,4),fftwork(ifft3),vis2_re(ifft3,4),vis2_im(ifft3,4),v1re(ifft3,4),v1im(ifft3,4))
129 :
130 0 : DO jspin = 1, input%jspins
131 0 : CALL fft3d(vis(:, jspin), fftwork, vTot%pw(:, jspin), stars, +1)
132 0 : CALL fft3d(v1re(:, jspin), v1im(:, jspin), vTot1%pw(:, jspin), starsq, +1)
133 : END DO
134 :
135 0 : CALL fft3d(vis(:, 3), vis(:, 4), vTot%pw(:, 3), stars, +1)
136 :
137 0 : DO imeshpt = 1, ifft3
138 0 : v11 = vis(imeshpt, 1)
139 0 : v22 = vis(imeshpt, 2)
140 0 : v21 = vis(imeshpt, 3) + ImagUnit * vis(imeshpt, 4)
141 0 : v12 = vis(imeshpt, 3) - ImagUnit * vis(imeshpt, 4)
142 :
143 0 : theta = den%theta_pw(imeshpt)
144 0 : phi = den%phi_pw(imeshpt)
145 :
146 0 : v1up = v1re(imeshpt,1) + ImagUnit * v1im(imeshpt,1)
147 0 : v1down = v1re(imeshpt,2) + ImagUnit * v1im(imeshpt,2)
148 :
149 0 : t1 = den1%theta_pw(imeshpt) + ImagUnit * den1im%theta_pw(imeshpt)
150 0 : p1 = den1%phi_pw(imeshpt) + ImagUnit * den1im%phi_pw(imeshpt)
151 :
152 0 : v1 = (v1up + v1down) / 2.0
153 0 : b1 = (v1up - v1down) / 2.0
154 :
155 0 : a11 = -ImagUnit * p1 / 2.0
156 0 : a22 = ImagUnit * p1 / 2.0
157 0 : a21 = EXP( ImagUnit*phi) * t1 / 2.0
158 0 : a12 = -EXP(-ImagUnit*phi) * t1 / 2.0
159 :
160 0 : v1mat11 = v1 + b1 * COS(theta) !11
161 0 : v1mat22 = v1 - b1 * COS(theta) !22
162 0 : v1mat21 = b1 * SIN(theta)*EXP( Imagunit*phi) !21
163 0 : v1mat12 = b1 * SIN(theta)*EXP(-Imagunit*phi) !12
164 :
165 0 : av11 = a11 * v11 + a12 * v21 !11
166 0 : av22 = a21 * v12 + a22 * v22 !22
167 0 : av21 = a21 * v11 + a22 * v21 !21
168 0 : av12 = a11 * v12 + a12 * v22 !12
169 :
170 0 : v1mat11 = v1mat11 + av11 + CONJG(av11)
171 0 : v1mat22 = v1mat22 + av22 + CONJG(av22)
172 0 : v1mat21 = v1mat21 + av21 + CONJG(av12)
173 0 : v1mat12 = v1mat12 + av12 + CONJG(av21)
174 :
175 0 : vis_re(imeshpt, 1) = REAL(v1mat11)
176 0 : vis_re(imeshpt, 2) = REAL(v1mat22)
177 0 : vis_re(imeshpt, 3) = REAL(v1mat21)
178 0 : vis_re(imeshpt, 4) = REAL(v1mat12)
179 :
180 0 : vis2_re(imeshpt, 1) = REAL(v1mat11 * stars%ufft(imeshpt-1) + v11 * starsq%ufft(imeshpt-1))
181 0 : vis2_re(imeshpt, 2) = REAL(v1mat22 * stars%ufft(imeshpt-1) + v22 * starsq%ufft(imeshpt-1))
182 0 : vis2_re(imeshpt, 3) = REAL(v1mat21 * stars%ufft(imeshpt-1) + v21 * starsq%ufft(imeshpt-1))
183 0 : vis2_re(imeshpt, 4) = REAL(v1mat12 * stars%ufft(imeshpt-1) + v12 * starsq%ufft(imeshpt-1))
184 :
185 0 : vis_im(imeshpt, 1) = AIMAG(v1mat11)
186 0 : vis_im(imeshpt, 2) = AIMAG(v1mat22)
187 0 : vis_im(imeshpt, 3) = AIMAG(v1mat21)
188 0 : vis_im(imeshpt, 4) = AIMAG(v1mat12)
189 :
190 0 : vis2_im(imeshpt, 1) = AIMAG(v1mat11 * stars%ufft(imeshpt-1) + v11 * starsq%ufft(imeshpt-1))
191 0 : vis2_im(imeshpt, 2) = AIMAG(v1mat22 * stars%ufft(imeshpt-1) + v22 * starsq%ufft(imeshpt-1))
192 0 : vis2_im(imeshpt, 3) = AIMAG(v1mat21 * stars%ufft(imeshpt-1) + v21 * starsq%ufft(imeshpt-1))
193 0 : vis2_im(imeshpt, 4) = AIMAG(v1mat12 * stars%ufft(imeshpt-1) + v12 * starsq%ufft(imeshpt-1))
194 :
195 : END DO
196 :
197 0 : DO ipot = 1, 4
198 0 : CALL fft3d(vis_re(:, ipot), vis_im(:, ipot), vTot1%pw(1, ipot), starsq, -1)
199 0 : CALL fft3d(vis2_re(:, ipot), vis2_im(:, ipot), vTot1%pw_w(1, ipot), starsq, -1)
200 : END DO
201 :
202 0 : END SUBROUTINE get_int_global_perturbation
203 : END MODULE m_get_int_perturbation
|
