Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Does all kind of post scf calculations for GPW/GAPW
10 : !> \par History
11 : !> Taken out from qs_scf_post_gpw
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_elf_methods
15 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
16 : USE kinds, ONLY: dp
17 : USE mathconstants, ONLY: pi
18 : USE pw_env_types, ONLY: pw_env_get,&
19 : pw_env_type
20 : USE pw_methods, ONLY: pw_derive,&
21 : pw_transfer,&
22 : pw_zero
23 : USE pw_pool_types, ONLY: pw_pool_p_type,&
24 : pw_pool_type
25 : USE pw_types, ONLY: pw_c1d_gs_type,&
26 : pw_r3d_rs_type
27 : USE qs_collocate_density, ONLY: calculate_rho_elec
28 : USE qs_environment_types, ONLY: get_qs_env,&
29 : qs_environment_type
30 : USE qs_ks_types, ONLY: qs_ks_env_type
31 : USE qs_rho_types, ONLY: qs_rho_get,&
32 : qs_rho_type
33 : #include "./base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 : PRIVATE
37 :
38 : ! Global parameters
39 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_elf_methods'
40 :
41 : PUBLIC :: qs_elf_calc
42 :
43 : ! **************************************************************************************************
44 :
45 : CONTAINS
46 :
47 : ! **************************************************************************************************
48 : !> \brief ...
49 : !> \param qs_env ...
50 : !> \param elf_r ...
51 : !> \param rho_cutoff ...
52 : ! **************************************************************************************************
53 82 : SUBROUTINE qs_elf_calc(qs_env, elf_r, rho_cutoff)
54 :
55 : TYPE(qs_environment_type), POINTER :: qs_env
56 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: elf_r
57 : REAL(kind=dp), INTENT(IN) :: rho_cutoff
58 :
59 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_elf_calc'
60 : INTEGER, DIMENSION(3, 3), PARAMETER :: nd = RESHAPE([1, 0, 0, 0, 1, 0, 0, 0, 1], [3, 3])
61 : REAL(KIND=dp), PARAMETER :: ELFCUT = 0.0001_dp, &
62 : f18 = (1.0_dp/8.0_dp), &
63 : f23 = (2.0_dp/3.0_dp), &
64 : f53 = (5.0_dp/3.0_dp)
65 :
66 : INTEGER :: handle, i, idir, ispin, j, k, nspin
67 : INTEGER, DIMENSION(2, 3) :: bo
68 : LOGICAL :: deriv_pw, drho_r_valid, tau_r_valid
69 : REAL(kind=dp) :: cfermi, elf_kernel, norm_drho, rho_53, &
70 : udvol
71 82 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
72 82 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_struct_ao
73 : TYPE(pw_c1d_gs_type) :: tmp_g
74 : TYPE(pw_env_type), POINTER :: pw_env
75 82 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
76 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
77 328 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: drho_r
78 82 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_struct_r, tau_struct_r
79 82 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_struct_r
80 : TYPE(pw_r3d_rs_type), POINTER :: rho_r, tau_r
81 : TYPE(qs_ks_env_type), POINTER :: ks_env
82 : TYPE(qs_rho_type), POINTER :: rho_struct
83 :
84 82 : CALL timeset(routineN, handle)
85 :
86 82 : NULLIFY (rho_struct, rho_r, tau_r, pw_env, auxbas_pw_pool, pw_pools, ks_env)
87 82 : NULLIFY (rho_struct_ao, rho_struct_r, tau_struct_r, drho_struct_r)
88 :
89 82 : CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env, rho=rho_struct)
90 :
91 : CALL qs_rho_get(rho_struct, &
92 : rho_ao_kp=rho_struct_ao, &
93 : rho_r=rho_struct_r, &
94 : tau_r=tau_struct_r, &
95 : drho_r=drho_struct_r, &
96 : tau_r_valid=tau_r_valid, &
97 82 : drho_r_valid=drho_r_valid)
98 :
99 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
100 82 : pw_pools=pw_pools)
101 82 : nspin = SIZE(rho_struct_r)
102 820 : bo = rho_struct_r(1)%pw_grid%bounds_local
103 82 : cfermi = (3.0_dp/10.0_dp)*(pi*pi*3.0_dp)**f23
104 :
105 : ! In this case, we need a work matrix containing tau in g space
106 : ! We will not have further use for it, so we will need only one
107 82 : IF (.NOT. tau_r_valid) THEN
108 82 : ALLOCATE (tau_r)
109 82 : CALL auxbas_pw_pool%create_pw(tau_r)
110 : END IF
111 82 : IF (.NOT. tau_r_valid .OR. .NOT. drho_r_valid) THEN
112 82 : CALL auxbas_pw_pool%create_pw(tmp_g)
113 : END IF
114 82 : IF (.NOT. drho_r_valid) THEN
115 328 : DO idir = 1, 3
116 328 : CALL auxbas_pw_pool%create_pw(drho_r(idir))
117 : END DO
118 : END IF
119 :
120 168 : DO ispin = 1, nspin
121 86 : rho_r => rho_struct_r(ispin)
122 86 : IF (tau_r_valid) THEN
123 0 : tau_r => tau_struct_r(ispin)
124 : ELSE
125 86 : rho_ao => rho_struct_ao(ispin, :)
126 86 : CALL pw_zero(tau_r)
127 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
128 : rho=tau_r, &
129 : rho_gspace=tmp_g, &
130 : ks_env=ks_env, soft_valid=.FALSE., &
131 86 : compute_tau=.TRUE.)
132 : END IF
133 :
134 86 : IF (drho_r_valid) THEN
135 0 : drho_r(:) = drho_struct_r(:, ispin)
136 : ELSE
137 86 : deriv_pw = .FALSE.
138 : IF (deriv_pw) THEN
139 : udvol = 1.0_dp/rho_r%pw_grid%dvol
140 : DO idir = 1, 3
141 : CALL pw_transfer(rho_r, tmp_g)
142 : CALL pw_derive(tmp_g, nd(:, idir))
143 : CALL pw_transfer(tmp_g, drho_r(idir))
144 : END DO
145 :
146 : ELSE
147 344 : DO idir = 1, 3
148 258 : rho_ao => rho_struct_ao(ispin, :)
149 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
150 : rho=drho_r(idir), &
151 : rho_gspace=tmp_g, &
152 : ks_env=ks_env, soft_valid=.FALSE., &
153 344 : compute_tau=.FALSE., compute_grad=.TRUE., idir=idir)
154 :
155 : END DO
156 : END IF
157 : END IF
158 :
159 : ! Calculate elf_r
160 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(bo,elf_r, ispin, drho_r,rho_r, tau_r, cfermi, rho_cutoff)&
161 168 : !$OMP PRIVATE(k,j,i, norm_drho, rho_53, elf_kernel)
162 : DO k = bo(1, 3), bo(2, 3)
163 : DO j = bo(1, 2), bo(2, 2)
164 : DO i = bo(1, 1), bo(2, 1)
165 : norm_drho = drho_r(1)%array(i, j, k)**2 + &
166 : drho_r(2)%array(i, j, k)**2 + &
167 : drho_r(3)%array(i, j, k)**2
168 : norm_drho = norm_drho/MAX(rho_r%array(i, j, k), rho_cutoff)
169 : rho_53 = cfermi*MAX(rho_r%array(i, j, k), rho_cutoff)**f53
170 : elf_kernel = (tau_r%array(i, j, k) - f18*norm_drho) + 2.87E-5_dp
171 : elf_kernel = (elf_kernel/rho_53)**2
172 : elf_r(ispin)%array(i, j, k) = 1.0_dp/(1.0_dp + elf_kernel)
173 : IF (elf_r(ispin)%array(i, j, k) < ELFCUT) elf_r(ispin)%array(i, j, k) = 0.0_dp
174 : END DO
175 : END DO
176 : END DO
177 : END DO
178 :
179 82 : IF (.NOT. drho_r_valid) THEN
180 328 : DO idir = 1, 3
181 328 : CALL auxbas_pw_pool%give_back_pw(drho_r(idir))
182 : END DO
183 : END IF
184 82 : IF (.NOT. tau_r_valid) THEN
185 82 : CALL auxbas_pw_pool%give_back_pw(tau_r)
186 82 : DEALLOCATE (tau_r)
187 : END IF
188 82 : IF (.NOT. tau_r_valid .OR. .NOT. drho_r_valid) THEN
189 82 : CALL auxbas_pw_pool%give_back_pw(tmp_g)
190 : END IF
191 :
192 82 : CALL timestop(handle)
193 :
194 82 : END SUBROUTINE qs_elf_calc
195 :
196 : END MODULE qs_elf_methods
|
