-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdensity.f90
114 lines (90 loc) · 3.59 KB
/
density.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
module density
use types
use gvect
use fft
use pseudopot
use constants
implicit none
complex(dp),allocatable :: density_r(:,:,:) !density in real space
complex(dp),allocatable :: density_g(:,:,:) !density in reciprocal space
complex(dp),allocatable :: core_density_g(:,:,:), total_density_g(:,:,:)
real(dp),allocatable :: FillingFactor(:)
contains
subroutine init_density(Nx,Ny,Nz,numGVects, params)
integer, intent(in) :: Nx,Ny,Nz,numGVects
type(GthPotParams) :: params
real(dp) :: chi, two_pi_over_l
chi = params%chi
two_pi_over_l = 2 * pi / params%box_length
allocate(density_r(0:Nx-1,0:Ny-1,0:Nz-1))
allocate(density_g(0:Nx-1,0:Ny-1,0:Nz-1))
allocate(core_density_g(0:Nx-1,0:Ny-1,0:Nz-1))
allocate(total_density_g(0:Nx-1,0:Ny-1,0:Nz-1))
density_r = 0.0_dp
density_g = 0.0_dp
allocate(FillingFactor(numGVects))
FillingFactor = 0
core_density_g = - params%Zeff/params%omega * exp(-0.5 *(g_grid_norm * chi * two_pi_over_l)**2 )
!core_density_g = core_density_g * structure_factor
end subroutine init_density
function layoutKIndexForFft(Nx,Ny,Nz,num_orbitals,C) result(retval)
!-----------------------------------------------------------
integer :: num_orbitals, Nx,Ny,Nz
complex(dp) :: C(num_orbitals)
complex(dp) :: retval(0:Nx-1,0:Ny-1,0:Nz-1)
integer :: size , i
integer :: toI,toJ, toK
retval = 0.0 !init to zero
do i = 1, num_orbitals
toI = g_indexes(1,i)
if(toI < 0) toI = toI + Nx
toJ = g_indexes(2,i)
if(toJ < 0) toJ = toJ + Ny
toK = g_indexes(3,i)
if(toK < 0) toK = toK + Nz
!print *, toI,toJ,toK, C(i)
retval(toI,toJ,toK) = C(i)
! print *,g_grid(:,toI,toJ,toK), g_indexes(:,i)
enddo
! ok, checkeado
end function layoutKIndexForFft
!---------------------------------------------
! compute the system density
! in real and complex espace
! upon the coefficients
! of the expansion of the wave
! functions in plane waves in
! reciprocal space
! --------------------------------------------
subroutine compute_density(N,C,omega)
implicit none
integer,intent(in) :: N !number of G vectors
complex(dp) :: C(N,N) ! the wave functions expansion coeffs, C(:,j) for the jth one
real(dp) :: omega ! the unit cell volumen
real(dp),allocatable :: Nj(:,:,:)
! aux variables for calculation
real(dp) :: val
complex(dp), allocatable :: psi_r(:,:,:)
complex(dp), allocatable :: psi_k(:,:,:)
integer :: i,k,l,m
allocate(psi_r(0:Nx-1,0:Ny-1,0:Nz-1))
density_r = 0.0_dp
do i = 1,numGVects
psi_k = layoutKIndexForFft(Nx,Ny,Nz,numGVects,C(:,i))
call fft_forward_3d(Nx,Ny,Nz,psi_k,psi_r) ! exp(1)
psi_r = psi_r/sqrt(omega)
val = sum(abs(psi_r)**2)*omega
density_r = density_r + FillingFactor(i)*psi_r*conjg(psi_r)
enddo
call fft_backward_3d(Nx,Ny,Nz,density_r,density_g) ! exp(-1)
density_g = density_g/(Nx*Ny*Nz)
end subroutine compute_density
!--------------------------------------
! map to linear layout, unimplemented
subroutine to_linear_layout(Nx,Ny,Nz, in, out)
!----------------------------------------
integer, intent(in) :: Nx,Ny,Nz
complex(dp), intent(in) :: in(Nx,Ny,Nz)
complex(dp), intent(out) :: out(Nx,Ny,Nz)
end subroutine to_linear_layout
end module density