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