1!
2! Copyright (C) 2017 Quantum ESPRESSO Foundation
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!-----------------------------------------------------------------------------
9MODULE fft_rho
10  !-----------------------------------------------------------------------------
11  !
12  ! ... FFT and inverse FFT of rho on the dense grid
13  !
14  USE kinds,     ONLY : DP
15  USE fft_interfaces, ONLY: fwfft, invfft
16  USE control_flags,  ONLY: gamma_only
17  !
18  IMPLICIT NONE
19  PRIVATE
20  PUBLIC :: rho_r2g, rho_g2r
21  !
22  INTERFACE rho_g2r
23    MODULE PROCEDURE rho_g2r_1, rho_g2r_2, rho_g2r_sum_components
24  END INTERFACE
25  !
26CONTAINS
27  !
28  SUBROUTINE rho_r2g ( desc, rhor, rhog, v )
29    USE fft_types,              ONLY: fft_type_descriptor
30    USE fft_helper_subroutines, ONLY: fftx_threed2oned
31    !
32    TYPE(fft_type_descriptor), INTENT(in) :: desc
33    REAL(dp),    INTENT(in) :: rhor(:,:)
34    COMPLEX(dp), INTENT(OUT):: rhog(:,:)
35    REAL(dp),    OPTIONAL, INTENT(in) :: v(:)
36    !
37    INTEGER :: ir, ig, iss, isup, isdw
38    INTEGER :: nspin
39    COMPLEX(dp):: fp, fm
40    COMPLEX(dp), ALLOCATABLE :: psi(:)
41
42    nspin= SIZE (rhor, 2)
43
44    ALLOCATE( psi( desc%nnr ) )
45    IF( nspin == 1 ) THEN
46       iss=1
47       IF( PRESENT( v ) ) THEN
48          DO ir=1,desc%nnr
49             psi(ir)=CMPLX(rhor(ir,iss)+v(ir),0.0_dp,kind=dp)
50          END DO
51       ELSE
52          DO ir=1,desc%nnr
53             psi(ir)=CMPLX(rhor(ir,iss),0.0_dp,kind=dp)
54          END DO
55       END IF
56       CALL fwfft('Rho', psi, desc )
57       CALL fftx_threed2oned( desc, psi, rhog(:,iss) )
58    ELSE
59       IF ( gamma_only ) THEN
60          ! nspin/2 = 1 for LSDA, = 2 for noncolinear
61          DO iss=1,nspin/2
62             isup=1+(iss-1)*nspin/2 ! 1 for LSDA, 1 and 3 for noncolinear
63             isdw=2+(iss-1)*nspin/2 ! 2 for LSDA, 2 and 4 for noncolinear
64             IF( PRESENT( v ) ) THEN
65                DO ir=1,desc%nnr
66                    psi(ir)=CMPLX(rhor(ir,isup)+v(ir),rhor(ir,isdw)+v(ir),kind=dp)
67                END DO
68             ELSE
69                DO ir=1,desc%nnr
70                   psi(ir)=CMPLX(rhor(ir,isup),rhor(ir,isdw),kind=dp)
71                END DO
72             END IF
73             CALL fwfft('Rho', psi, desc )
74             CALL fftx_threed2oned( desc, psi, rhog(:,isup), rhog(:,isdw) )
75          END DO
76       ELSE
77          DO iss=1,nspin
78             IF( PRESENT( v ) ) THEN
79                DO ir=1,desc%nnr
80                    psi(ir)=CMPLX(rhor(ir,iss)+v(ir),0.0_dp,kind=dp)
81                END DO
82             ELSE
83                DO ir=1,desc%nnr
84                   psi(ir)=CMPLX(rhor(ir,iss),0.0_dp,kind=dp)
85                END DO
86             END IF
87             CALL fwfft('Rho', psi, desc )
88             CALL fftx_threed2oned( desc, psi, rhog(:,iss) )
89          END DO
90       END IF
91    ENDIF
92
93    DEALLOCATE( psi )
94
95  END SUBROUTINE rho_r2g
96  !
97  SUBROUTINE rho_g2r_1 ( desc, rhog, rhor )
98    USE fft_types,              ONLY: fft_type_descriptor
99    USE fft_helper_subroutines, ONLY: fftx_threed2oned, fftx_oned2threed
100    !
101    TYPE(fft_type_descriptor), INTENT(in) :: desc
102    COMPLEX(dp), INTENT(in ):: rhog(:)
103    REAL(dp),    INTENT(out):: rhor(:)
104    !
105    INTEGER :: ir
106    COMPLEX(dp), ALLOCATABLE :: psi(:)
107
108    ALLOCATE( psi( desc%nnr ) )
109    CALL fftx_oned2threed( desc, psi, rhog )
110    CALL invfft('Rho',psi, desc )
111!$omp parallel do
112    DO ir=1,desc%nnr
113       rhor(ir)=DBLE(psi(ir))
114    END DO
115!$omp end parallel do
116
117    DEALLOCATE( psi )
118
119  END SUBROUTINE rho_g2r_1
120  !
121  SUBROUTINE rho_g2r_2 ( desc, rhog, rhor )
122    USE fft_types,              ONLY: fft_type_descriptor
123    USE fft_helper_subroutines, ONLY: fftx_threed2oned, fftx_oned2threed
124    !
125    TYPE(fft_type_descriptor), INTENT(in) :: desc
126    COMPLEX(dp), INTENT(in ):: rhog(:,:)
127    REAL(dp),    INTENT(out):: rhor(:,:)
128    !
129    INTEGER :: ir, ig, iss, isup, isdw
130    INTEGER :: nspin
131    COMPLEX(dp), ALLOCATABLE :: psi(:)
132
133    nspin= SIZE (rhog, 2)
134
135    ALLOCATE( psi( desc%nnr ) )
136    IF ( gamma_only ) THEN
137       IF( nspin == 1 ) THEN
138          iss=1
139          CALL fftx_oned2threed( desc, psi, rhog(:,iss) )
140          CALL invfft('Rho',psi, desc )
141!$omp parallel do
142          DO ir=1,desc%nnr
143             rhor(ir,iss)=DBLE(psi(ir))
144          END DO
145!$omp end parallel do
146       ELSE
147          ! nspin/2 = 1 for LSDA, = 2 for noncolinear
148          DO iss=1,nspin/2
149             isup=1+(iss-1)*nspin/2 ! 1 for LSDA, 1 and 3 for noncolinear
150             isdw=2+(iss-1)*nspin/2 ! 2 for LSDA, 2 and 4 for noncolinear
151             CALL fftx_oned2threed( desc, psi, rhog(:,isup), rhog(:,isdw) )
152             CALL invfft('Rho',psi, desc )
153!$omp parallel do
154             DO ir=1,desc%nnr
155                rhor(ir,isup)= DBLE(psi(ir))
156                rhor(ir,isdw)=AIMAG(psi(ir))
157             END DO
158!$omp end parallel do
159          END DO
160       ENDIF
161       !
162    ELSE
163       !
164       DO iss=1, nspin
165          CALL fftx_oned2threed( desc, psi, rhog(:,iss) )
166          CALL invfft('Rho',psi, desc )
167!$omp parallel do
168          DO ir=1,desc%nnr
169             rhor(ir,iss)=DBLE(psi(ir))
170          END DO
171!$omp end parallel do
172       END DO
173    END IF
174
175    DEALLOCATE( psi )
176
177  END SUBROUTINE rho_g2r_2
178  !
179  SUBROUTINE rho_g2r_sum_components ( desc, rhog, rhor )
180    USE fft_types,              ONLY: fft_type_descriptor
181    USE fft_helper_subroutines, ONLY: fftx_threed2oned, fftx_oned2threed
182    !
183    TYPE(fft_type_descriptor), INTENT(in) :: desc
184    COMPLEX(dp), INTENT(in ):: rhog(:,:)
185    REAL(dp),    INTENT(out):: rhor(:)
186    !
187    INTEGER :: ir, ig, iss, isup, isdw
188    INTEGER :: nspin
189    COMPLEX(dp), ALLOCATABLE :: psi(:)
190
191    nspin= SIZE (rhog, 2)
192
193    ALLOCATE( psi( desc%nnr ) )
194    IF ( gamma_only ) THEN
195       IF( nspin == 1 ) THEN
196          iss=1
197          CALL fftx_oned2threed( desc, psi, rhog(:,iss) )
198          CALL invfft('Rho',psi, desc )
199!$omp parallel do
200          DO ir=1,desc%nnr
201             rhor(ir)=DBLE(psi(ir))
202          END DO
203!$omp end parallel do
204       ELSE IF ( nspin == 2) THEN
205          isup=1
206          isdw=2
207          CALL fftx_oned2threed( desc, psi, rhog(:,isup), rhog(:,isdw) )
208          CALL invfft('Rho',psi, desc )
209!$omp parallel do
210          DO ir=1,desc%nnr
211             rhor(ir)= DBLE(psi(ir))+AIMAG(psi(ir))
212          END DO
213!$omp end parallel do
214       ELSE
215          CALL errore('rho_g2r_sum_components','noncolinear case?',nspin)
216       ENDIF
217       !
218    ELSE
219       !
220       DO iss=1, nspin
221          CALL fftx_oned2threed( desc, psi, rhog(:,iss) )
222          CALL invfft('Rho',psi, desc )
223          IF( iss == 1 ) THEN
224!$omp parallel do
225             DO ir=1,desc%nnr
226                rhor(ir)=DBLE(psi(ir))
227             END DO
228!$omp end parallel do
229          ELSE
230!$omp parallel do
231             DO ir=1,desc%nnr
232                rhor(ir)=rhor(ir) + DBLE(psi(ir))
233             END DO
234!$omp end parallel do
235          END IF
236       END DO
237    END IF
238
239    DEALLOCATE( psi )
240
241  END SUBROUTINE rho_g2r_sum_components
242
243END MODULE fft_rho
244