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 : MODULE m_step_function
7 : !! Contains revised subroutines to construct the step function \f$\Theta(r)\f$
8 : !! both on a fine real space grid as well as in \f$G<G_{max}\f$ space.
9 : USE m_juDFT
10 : USE m_types
11 :
12 : IMPLICIT NONE
13 :
14 : CONTAINS
15 0 : SUBROUTINE stepf_analytical(sym, stars, atoms, input, cell, fmpi, fftgrid, qvec, iDtype, iDir, iOrd, stepf_array)
16 : !! Construct the analytical representation of the step function on a big
17 : !! reciprocal grid.
18 :
19 : USE m_constants
20 :
21 : TYPE(t_sym), INTENT(IN) :: sym
22 : TYPE(t_stars), INTENT(IN) :: stars
23 : TYPE(t_atoms), INTENT(IN) :: atoms
24 : TYPE(t_input), INTENT(IN) :: input
25 : TYPE(t_cell), INTENT(IN) :: cell
26 : TYPE(t_mpi), INTENT(IN) :: fmpi
27 :
28 : TYPE(t_fftgrid), INTENT(OUT) :: fftgrid
29 :
30 : REAL, OPTIONAL, INTENT(IN) :: qvec(3)
31 : INTEGER, OPTIONAL, INTENT(IN) :: iDtype, iDir, iOrd
32 :
33 : COMPLEX, OPTIONAL, INTENT(INOUT) :: stepf_array(0:,:,:)
34 :
35 : INTEGER :: x, y, z, z2, z_min, z_max
36 : INTEGER :: layerDim, ifftd, gInd, n, na, nn, iDir_loc
37 : INTEGER :: n_lims(2), dir_lims(2)
38 : REAL :: th, inv_omtil, fJ
39 : REAL :: g2, absg, absgr, help, fp_omtil, radg, gExtx, gExty
40 : REAL :: r_c, r_phs
41 : COMPLEX :: c_c, c_phs
42 :
43 : REAL :: gInt(3), gExt(3)
44 :
45 0 : CALL timestart("New stepf")
46 0 : ifftd = 27*stars%mx1*stars%mx2*stars%mx3
47 0 : layerDim = 9*stars%mx1*stars%mx2
48 :
49 : ! Initialize fftgrid
50 0 : CALL fftgrid%init((/3*stars%mx1,3*stars%mx2,3*stars%mx3/))
51 :
52 0 : inv_omtil = 1.0/cell%omtil
53 0 : fp_omtil = -fpi_const*inv_omtil
54 :
55 0 : fftgrid%grid = CMPLX(0.0,0.0)
56 :
57 0 : IF (.NOT.PRESENT(qvec)) THEN
58 0 : DO n = 1, atoms%ntype
59 0 : fftgrid%grid(0) = fftgrid%grid(0) + atoms%neq(n)*atoms%volmts(n)
60 : END DO
61 0 : fftgrid%grid(0) = 1.0 - fftgrid%grid(0)*inv_omtil
62 : END IF
63 :
64 0 : z_min = 0
65 0 : z_max = 3*stars%mx3 - 1
66 :
67 0 : n_lims = [1, atoms%ntype]
68 0 : dir_lims = [1, 1]
69 :
70 0 : IF (PRESENT(qvec)) THEN
71 0 : n_lims(1) = MERGE(1, iDtype,PRESENT(stepf_array).AND.iOrd==1)
72 0 : n_lims(2) = MERGE(atoms%ntype,iDtype,PRESENT(stepf_array).AND.iOrd==1)
73 :
74 0 : dir_lims(1) = MERGE(1,iDir,PRESENT(stepf_array))
75 0 : dir_lims(2) = MERGE(3,iDir,PRESENT(stepf_array))
76 : END IF
77 :
78 0 : DO z = z_min, z_max
79 0 : gInt(3) = REAL(z)
80 0 : IF (2*z > 3*stars%mx3) gInt(3) = gInt(3) - 3.0*stars%mx3
81 0 : DO y = 0, 3*stars%mx2 - 1
82 0 : gInt(2) = REAL(y)
83 0 : IF (2*y > 3*stars%mx2) gInt(2) = gInt(2) - 3.0*stars%mx2
84 0 : x_loop: DO x = 0, 3*stars%mx1 - 1
85 : ! TODO: Maybe add the possibility to skip calculations for elements
86 : ! that are already available by inversion; problematic when
87 : ! parallelized. Also: parallelize!
88 0 : gInt(1) = REAL(x)
89 0 : IF (2*x > 3*stars%mx1) gInt(1) = gInt(1) - 3.0*stars%mx1
90 0 : gInd = x + 3*stars%mx1*y + layerDim*z
91 :
92 0 : IF (PRESENT(qvec)) THEN
93 0 : IF (ALL([x,y,z]==0).AND.norm2(qvec)<1e-9) CYCLE x_loop
94 : ELSE
95 0 : IF (ALL([x,y,z]==0)) CYCLE x_loop
96 : END IF
97 :
98 0 : gExt = MATMUL(TRANSPOSE(cell%bmat), gInt)
99 0 : IF (PRESENT(qvec)) gExt = gExt + MATMUL(TRANSPOSE(cell%bmat), qvec)
100 0 : g2 = DOT_PRODUCT(gExt, gExt)
101 0 : absg = SQRT(g2)
102 0 : help = fp_omtil/g2
103 :
104 0 : c_c = CMPLX(0.0, 0.0)
105 :
106 0 : DO n = n_lims(1), n_lims(2)
107 0 : DO iDir_loc = dir_lims(1), dir_lims(2)
108 0 : c_phs = CMPLX(0.0, 0.0)
109 0 : na = atoms%firstAtom(n) - 1
110 : !na = MERGE(atoms%firstAtom(n) - 1,iDtype,.NOT.PRESENT(qvec))
111 0 : DO nn = 1, atoms%neq(n)
112 0 : IF (.NOT.PRESENT(qvec)) THEN
113 0 : th = -tpi_const*DOT_PRODUCT(gInt, atoms%taual(:, na + nn))
114 0 : c_phs = c_phs + EXP(CMPLX(0, th))
115 : ELSE
116 0 : th = -tpi_const*DOT_PRODUCT(gInt+qvec, atoms%taual(:, na + nn))
117 0 : IF (iOrd==1) THEN
118 0 : c_phs = c_phs + (-ImagUnit*gExt(iDir_loc))*EXP(CMPLX(0, th))
119 : ELSE
120 0 : c_phs = c_phs + (-gExt(iDir)*gExt(iDir_loc))*EXP(CMPLX(0, th))
121 : END IF
122 : END IF
123 : END DO
124 0 : absgr = absg*atoms%rmt(n)
125 0 : c_c = c_c + atoms%rmt(n)*(SIN(absgr)/absgr - COS(absgr))*c_phs
126 :
127 0 : IF (PRESENT(stepf_array)) THEN
128 0 : stepf_array(gInd, n, iDir_loc) = help * c_c
129 0 : c_c = CMPLX(0.0, 0.0)
130 : END IF
131 : END DO
132 : END DO
133 :
134 0 : IF (.NOT.PRESENT(stepf_array)) THEN
135 0 : fftgrid%grid(gInd) = help * c_c
136 :
137 0 : IF (((2*z .EQ. 3*stars%mx3) .OR. (2*y .EQ. 3*stars%mx2)) .OR. (2*x .EQ. 3*stars%mx1)) THEN
138 0 : fftgrid%grid(gInd) = CMPLX(0.0,0.0)
139 : END IF
140 : ELSE
141 0 : IF (((2*z .EQ. 3*stars%mx3) .OR. (2*y .EQ. 3*stars%mx2)) .OR. (2*x .EQ. 3*stars%mx1)) THEN
142 0 : stepf_array(gInd, :, :) = CMPLX(0.0,0.0)
143 : END IF
144 : END IF
145 :
146 : END DO x_loop ! 0, 3*stars%mx1 - 1
147 : END DO ! 0, 3*stars%mx2 - 1
148 : END DO ! 0, 3*stars%mx3 - 1
149 :
150 : !IF (fmpi%irank == 0) THEN
151 : ! IF (input%film) THEN
152 : ! fftgrid%grid(0) = fftgrid%grid(0) + cell%vol*inv_omtil - 1.0
153 : ! DO z2 = 1, 3*stars%mx3 - 1
154 : ! gInt(3) = REAL(z2)
155 : ! IF (gInt(3) > 1.5*stars%mx3) gInt(3) = gInt(3) - 3.0 * stars%mx3
156 : ! th = cell%bmat(3, 3)*gInt(3)*cell%z1
157 : ! fftgrid%grid(z2*layerDim) = fftgrid%grid(z2*layerDim) + cell%vol*inv_omtil*SIN(th)/th
158 : ! END DO
159 : ! END IF
160 : !END IF
161 :
162 0 : CALL timestop("New stepf")
163 0 : END SUBROUTINE
164 :
165 0 : SUBROUTINE stepf_stars(stars,fftgrid,qvec)
166 :
167 : TYPE(t_stars), INTENT(INOUT) :: stars
168 : TYPE(t_fftgrid), INTENT(INOUT) :: fftgrid
169 : REAL, OPTIONAL, INTENT(IN) :: qvec(3)
170 :
171 : !#ifdef CPP_MPI
172 : ! INTEGER :: ierr
173 : ! CALL MPI_BCAST(stars%ng3, 1, MPI_INTEGER, 0, fmpi%mpi_comm, ierr)
174 : ! CALL MPI_BCAST(stars%ig, stars%ng3, MPI_INTEGER, 0, fmpi%mpi_comm, ierr)
175 : !#endif
176 :
177 : ! Save G coefficients to ustep
178 0 : CALL fftgrid%takeFieldFromGrid(stars, stars%ustep)
179 0 : stars%ustep = stars%ustep * 3 * stars%mx1 * 3 * stars%mx2 * 3 * stars%mx3
180 0 : CALL fftgrid%perform_fft(forward=.FALSE.)
181 :
182 0 : IF (.NOT.PRESENT(qvec)) THEN
183 0 : stars%ufft=fftgrid%grid
184 : ELSE
185 0 : stars%ufft1=fftgrid%grid
186 : END IF
187 :
188 0 : END SUBROUTINE
189 : END MODULE m_step_function
|
