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 g tensor calculation by dfpt
10 : !> Initialization of the epr_env, creation of the special neighbor lists
11 : !> Perturbation Hamiltonians by application of the p and rxp oprtators to psi0
12 : !> Write output
13 : !> Deallocate everything
14 : !> \note
15 : !> The psi0 should be localized
16 : !> the Sebastiani method works within the assumption that the orbitals are
17 : !> completely contained in the simulation box
18 : !> \par History
19 : !> created 07-2005 [MI]
20 : !> \author MI
21 : ! **************************************************************************************************
22 : MODULE qs_linres_epr_utils
23 : USE atomic_kind_types, ONLY: atomic_kind_type
24 : USE cell_types, ONLY: cell_type
25 : USE cp_control_types, ONLY: dft_control_type
26 : USE cp_log_handling, ONLY: cp_get_default_logger,&
27 : cp_logger_type
28 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29 : cp_print_key_unit_nr
30 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
31 : section_vals_type
32 : USE kinds, ONLY: dp
33 : USE mathconstants, ONLY: fourpi,&
34 : twopi
35 : USE particle_types, ONLY: particle_type
36 : USE physcon, ONLY: a_fine,&
37 : e_gfactor
38 : USE pw_env_types, ONLY: pw_env_get,&
39 : pw_env_type
40 : USE pw_pool_types, ONLY: pw_pool_type
41 : USE pw_types, ONLY: pw_c1d_gs_type,&
42 : pw_r3d_rs_type
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_kind_types, ONLY: qs_kind_type
46 : USE qs_linres_types, ONLY: deallocate_nablavks_atom_set,&
47 : epr_env_type,&
48 : init_nablavks_atom_set,&
49 : linres_control_type,&
50 : nablavks_atom_type,&
51 : set_epr_env
52 : USE qs_matrix_pools, ONLY: qs_matrix_pools_type
53 : USE qs_mo_types, ONLY: mo_set_type
54 : USE qs_rho_atom_types, ONLY: deallocate_rho_atom_set
55 : USE qs_rho_types, ONLY: qs_rho_clear,&
56 : qs_rho_create,&
57 : qs_rho_set
58 : USE scf_control_types, ONLY: scf_control_type
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 :
63 : PRIVATE
64 :
65 : PUBLIC :: epr_env_cleanup, epr_env_init
66 :
67 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_utils'
68 :
69 : CONTAINS
70 :
71 : ! **************************************************************************************************
72 : !> \brief Initialize the epr environment
73 : !> \param epr_env ...
74 : !> \param qs_env ...
75 : !> \par History
76 : !> 07.2006 created [MI]
77 : !> \author MI
78 : ! **************************************************************************************************
79 14 : SUBROUTINE epr_env_init(epr_env, qs_env)
80 : !
81 : TYPE(epr_env_type) :: epr_env
82 : TYPE(qs_environment_type), POINTER :: qs_env
83 :
84 : CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_env_init'
85 :
86 : INTEGER :: handle, i_B, idir, ispin, n_mo(2), nao, &
87 : natom, nmoloc, nspins, output_unit
88 : LOGICAL :: gapw
89 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
90 : TYPE(cell_type), POINTER :: cell
91 : TYPE(cp_logger_type), POINTER :: logger
92 : TYPE(dft_control_type), POINTER :: dft_control
93 : TYPE(linres_control_type), POINTER :: linres_control
94 14 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
95 14 : TYPE(nablavks_atom_type), DIMENSION(:), POINTER :: nablavks_atom_set
96 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
97 14 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
98 : TYPE(pw_env_type), POINTER :: pw_env
99 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
100 14 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
101 14 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
102 : TYPE(qs_matrix_pools_type), POINTER :: mpools
103 : TYPE(scf_control_type), POINTER :: scf_control
104 : TYPE(section_vals_type), POINTER :: lr_section
105 :
106 14 : CALL timeset(routineN, handle)
107 :
108 14 : NULLIFY (atomic_kind_set, qs_kind_set, cell, dft_control, linres_control, scf_control)
109 14 : NULLIFY (logger, mos, mpools, particle_set)
110 14 : NULLIFY (auxbas_pw_pool, pw_env)
111 14 : NULLIFY (nablavks_atom_set)
112 :
113 : n_mo(1:2) = 0
114 14 : nao = 0
115 14 : nmoloc = 0
116 :
117 14 : logger => cp_get_default_logger()
118 14 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
119 :
120 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
121 14 : extension=".linresLog")
122 :
123 14 : CALL epr_env_cleanup(epr_env)
124 :
125 14 : IF (output_unit > 0) THEN
126 7 : WRITE (output_unit, "(/,T20,A,/)") "*** Start EPR g tensor calculation ***"
127 7 : WRITE (output_unit, "(T10,A,/)") "Initialization of the EPR environment"
128 : END IF
129 :
130 : CALL get_qs_env(qs_env=qs_env, &
131 : atomic_kind_set=atomic_kind_set, &
132 : qs_kind_set=qs_kind_set, &
133 : cell=cell, &
134 : dft_control=dft_control, &
135 : linres_control=linres_control, &
136 : mos=mos, &
137 : mpools=mpools, &
138 : particle_set=particle_set, &
139 : pw_env=pw_env, &
140 14 : scf_control=scf_control)
141 : !
142 : ! Check if restat also psi0 should be restarted
143 : !IF(epr_env%restart_epr .AND. scf_control%density_guess/=restart_guess)THEN
144 : ! CPABORT("restart_epr requires density_guess=restart")
145 : !ENDIF
146 : !
147 : ! check that the psi0 are localized and you have all the centers
148 14 : CPASSERT(linres_control%localized_psi0)
149 14 : IF (output_unit > 0) THEN
150 : WRITE (output_unit, '(A)') &
151 7 : ' To get EPR parameters within PBC you need localized zero order orbitals '
152 : END IF
153 14 : gapw = dft_control%qs_control%gapw
154 14 : nspins = dft_control%nspins
155 14 : natom = SIZE(particle_set, 1)
156 : !
157 : ! Conversion factors
158 : ! Magical constant twopi/cell%deth just like in NMR shift (basically undo scale_fac in qs_linres_nmr_current.F)
159 14 : epr_env%g_free_factor = -1.0_dp*e_gfactor
160 14 : epr_env%g_zke_factor = e_gfactor*(a_fine)**2
161 14 : epr_env%g_so_factor = (a_fine)**2*(-1.0_dp*e_gfactor - 1.0_dp)/2.0_dp*twopi/cell%deth
162 14 : epr_env%g_so_factor_gapw = (a_fine)**2*(-1.0_dp*e_gfactor - 1.0_dp)/2.0_dp
163 : ! * 2 because B_ind = 2 * B_beta
164 14 : epr_env%g_soo_factor = 2.0_dp*fourpi*(a_fine)**2*twopi/cell%deth
165 : ! 2 * 2 * 1/4 * e^2 / m * a_0^2 * 2/3 * mu_0 / (omega * 1e-30 )
166 14 : epr_env%g_soo_chicorr_factor = 2.0/3.0_dp*fourpi*(a_fine)**2/cell%deth
167 : !
168 : ! If the current density on the grid needs to be stored
169 14 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
170 : !
171 : ! Initialize local current density if GAPW calculation
172 14 : IF (gapw) THEN
173 10 : CALL init_nablavks_atom_set(nablavks_atom_set, atomic_kind_set, qs_kind_set, nspins)
174 : CALL set_epr_env(epr_env=epr_env, &
175 10 : nablavks_atom_set=nablavks_atom_set)
176 : END IF
177 : !
178 : ! Bind
179 182 : ALLOCATE (epr_env%bind_set(3, 3))
180 56 : DO i_B = 1, 3
181 182 : DO idir = 1, 3
182 126 : NULLIFY (epr_env%bind_set(idir, i_B)%rho, rho_r, rho_g)
183 126 : ALLOCATE (epr_env%bind_set(idir, i_B)%rho)
184 126 : CALL qs_rho_create(epr_env%bind_set(idir, i_B)%rho)
185 378 : ALLOCATE (rho_r(1), rho_g(1))
186 126 : CALL auxbas_pw_pool%create_pw(rho_r(1))
187 126 : CALL auxbas_pw_pool%create_pw(rho_g(1))
188 168 : CALL qs_rho_set(epr_env%bind_set(idir, i_B)%rho, rho_r=rho_r, rho_g=rho_g)
189 : END DO
190 : END DO
191 :
192 : ! Nabla_V_ks
193 154 : ALLOCATE (epr_env%nablavks_set(3, dft_control%nspins))
194 56 : DO idir = 1, 3
195 140 : DO ispin = 1, nspins
196 84 : NULLIFY (epr_env%nablavks_set(idir, ispin)%rho, rho_r, rho_g)
197 84 : ALLOCATE (epr_env%nablavks_set(idir, ispin)%rho)
198 84 : CALL qs_rho_create(epr_env%nablavks_set(idir, ispin)%rho)
199 252 : ALLOCATE (rho_r(1), rho_g(1))
200 84 : CALL auxbas_pw_pool%create_pw(rho_r(1))
201 84 : CALL auxbas_pw_pool%create_pw(rho_g(1))
202 : CALL qs_rho_set(epr_env%nablavks_set(idir, ispin)%rho, &
203 126 : rho_r=rho_r, rho_g=rho_g)
204 : END DO
205 : END DO
206 :
207 : ! Initialize the g tensor components
208 14 : ALLOCATE (epr_env%g_total(3, 3))
209 14 : ALLOCATE (epr_env%g_so(3, 3))
210 14 : ALLOCATE (epr_env%g_soo(3, 3))
211 182 : epr_env%g_total = 0.0_dp
212 14 : epr_env%g_zke = 0.0_dp
213 182 : epr_env%g_so = 0.0_dp
214 182 : epr_env%g_soo = 0.0_dp
215 :
216 : CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
217 14 : & "PRINT%PROGRAM_RUN_INFO")
218 :
219 14 : CALL timestop(handle)
220 :
221 14 : END SUBROUTINE epr_env_init
222 :
223 : ! **************************************************************************************************
224 : !> \brief Deallocate the epr environment
225 : !> \param epr_env ...
226 : !> \par History
227 : !> 07.2005 created [MI]
228 : !> \author MI
229 : ! **************************************************************************************************
230 28 : SUBROUTINE epr_env_cleanup(epr_env)
231 :
232 : TYPE(epr_env_type) :: epr_env
233 :
234 : INTEGER :: i_B, idir, ispin
235 :
236 : ! nablavks_set
237 28 : IF (ASSOCIATED(epr_env%nablavks_set)) THEN
238 42 : DO ispin = 1, SIZE(epr_env%nablavks_set, 2)
239 126 : DO idir = 1, SIZE(epr_env%nablavks_set, 1)
240 84 : CALL qs_rho_clear(epr_env%nablavks_set(idir, ispin)%rho)
241 112 : DEALLOCATE (epr_env%nablavks_set(idir, ispin)%rho)
242 : END DO
243 : END DO
244 14 : DEALLOCATE (epr_env%nablavks_set)
245 : END IF
246 : ! nablavks_atom_set
247 28 : IF (ASSOCIATED(epr_env%nablavks_atom_set)) THEN
248 10 : CALL deallocate_nablavks_atom_set(epr_env%nablavks_atom_set)
249 : END IF
250 : ! vks_atom_set
251 28 : IF (ASSOCIATED(epr_env%vks_atom_set)) THEN
252 0 : CALL deallocate_rho_atom_set(epr_env%vks_atom_set)
253 : END IF
254 : ! bind_set
255 28 : IF (ASSOCIATED(epr_env%bind_set)) THEN
256 56 : DO i_B = 1, SIZE(epr_env%bind_set, 2)
257 182 : DO idir = 1, SIZE(epr_env%bind_set, 1)
258 126 : CALL qs_rho_clear(epr_env%bind_set(idir, i_B)%rho)
259 168 : DEALLOCATE (epr_env%bind_set(idir, i_B)%rho)
260 : END DO
261 : END DO
262 14 : DEALLOCATE (epr_env%bind_set)
263 : END IF
264 : ! bind_atom_set
265 28 : IF (ASSOCIATED(epr_env%bind_atom_set)) THEN
266 0 : DEALLOCATE (epr_env%bind_atom_set)
267 : END IF
268 : ! g_total
269 28 : IF (ASSOCIATED(epr_env%g_total)) THEN
270 14 : DEALLOCATE (epr_env%g_total)
271 : END IF
272 : ! g_so
273 28 : IF (ASSOCIATED(epr_env%g_so)) THEN
274 14 : DEALLOCATE (epr_env%g_so)
275 : END IF
276 : ! g_soo
277 28 : IF (ASSOCIATED(epr_env%g_soo)) THEN
278 14 : DEALLOCATE (epr_env%g_soo)
279 : END IF
280 :
281 28 : END SUBROUTINE epr_env_cleanup
282 :
283 : END MODULE qs_linres_epr_utils
|
