Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Contains the setup for the calculation of properties by linear response
10 : !> by the application of second order density functional perturbation theory.
11 : !> The knowledge of the ground state energy, density and wavefunctions is assumed.
12 : !> Uses the self consistent approach.
13 : !> Properties that can be calculated : none
14 : !> \par History
15 : !> created 06-2005 [MI]
16 : !> \author MI
17 : ! **************************************************************************************************
18 : MODULE qs_linres_module
19 : USE bibliography, ONLY: Ditler2021,&
20 : Ditler2022,&
21 : Weber2009,&
22 : cite_reference
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_type
27 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
28 : cp_print_key_unit_nr
29 : USE force_env_types, ONLY: force_env_get,&
30 : force_env_type,&
31 : use_qmmm,&
32 : use_qs_force
33 : USE input_constants, ONLY: lr_current,&
34 : lr_none,&
35 : ot_precond_full_all,&
36 : ot_precond_full_kinetic,&
37 : ot_precond_full_single,&
38 : ot_precond_full_single_inverse,&
39 : ot_precond_none,&
40 : ot_precond_s_inverse
41 : USE input_section_types, ONLY: section_vals_get,&
42 : section_vals_get_subs_vals,&
43 : section_vals_type,&
44 : section_vals_val_get
45 : USE kinds, ONLY: dp
46 : USE qs_dcdr, ONLY: apt_dR,&
47 : apt_dR_localization,&
48 : dcdr_build_op_dR,&
49 : dcdr_response_dR,&
50 : prepare_per_atom
51 : USE qs_dcdr_utils, ONLY: dcdr_env_cleanup,&
52 : dcdr_env_init,&
53 : dcdr_print
54 : USE qs_density_matrices, ONLY: calculate_density_matrix
55 : USE qs_environment_types, ONLY: get_qs_env,&
56 : qs_environment_type,&
57 : set_qs_env
58 : USE qs_linres_current, ONLY: current_build_chi,&
59 : current_build_current
60 : USE qs_linres_current_utils, ONLY: current_env_cleanup,&
61 : current_env_init,&
62 : current_response
63 : USE qs_linres_epr_nablavks, ONLY: epr_nablavks
64 : USE qs_linres_epr_ownutils, ONLY: epr_g_print,&
65 : epr_g_so,&
66 : epr_g_soo,&
67 : epr_g_zke,&
68 : epr_ind_magnetic_field
69 : USE qs_linres_epr_utils, ONLY: epr_env_cleanup,&
70 : epr_env_init
71 : USE qs_linres_issc_utils, ONLY: issc_env_cleanup,&
72 : issc_env_init,&
73 : issc_issc,&
74 : issc_print,&
75 : issc_response
76 : USE qs_linres_methods, ONLY: linres_localize
77 : USE qs_linres_nmr_shift, ONLY: nmr_shift,&
78 : nmr_shift_print
79 : USE qs_linres_nmr_utils, ONLY: nmr_env_cleanup,&
80 : nmr_env_init
81 : USE qs_linres_op, ONLY: current_operators,&
82 : issc_operators,&
83 : polar_operators,&
84 : polar_operators_local,&
85 : polar_operators_local_wannier
86 : USE qs_linres_polar_utils, ONLY: polar_env_init,&
87 : polar_polar,&
88 : polar_print,&
89 : polar_response
90 : USE qs_linres_types, ONLY: &
91 : current_env_type, dcdr_env_type, epr_env_type, get_polar_env, issc_env_type, &
92 : linres_control_type, nmr_env_type, polar_env_type, vcd_env_type
93 : USE qs_mfp, ONLY: mfp_aat,&
94 : mfp_build_operator_gauge_dependent,&
95 : mfp_build_operator_gauge_independent,&
96 : mfp_response
97 : USE qs_mo_types, ONLY: mo_set_type
98 : USE qs_p_env_methods, ONLY: p_env_create,&
99 : p_env_psi0_changed
100 : USE qs_p_env_types, ONLY: p_env_release,&
101 : qs_p_env_type
102 : USE qs_rho_methods, ONLY: qs_rho_update_rho
103 : USE qs_rho_types, ONLY: qs_rho_get,&
104 : qs_rho_type
105 : USE qs_vcd, ONLY: aat_dV,&
106 : apt_dV,&
107 : prepare_per_atom_vcd,&
108 : vcd_build_op_dV,&
109 : vcd_response_dV
110 : USE qs_vcd_utils, ONLY: vcd_env_cleanup,&
111 : vcd_env_init,&
112 : vcd_print
113 : USE xc, ONLY: xc_uses_kinetic_energy_density
114 : #include "./base/base_uses.f90"
115 :
116 : IMPLICIT NONE
117 :
118 : PRIVATE
119 : PUBLIC :: linres_calculation, linres_calculation_low
120 :
121 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_module'
122 : CHARACTER(len=*), PARAMETER, PRIVATE :: dcdr_meta_gga_error = &
123 : "Analytical DCDR is not implemented for functionals that depend on the kinetic energy density. "// &
124 : "Use PROPERTIES%LINRES%DCDR%APT_FD T to calculate APTs by finite differences."
125 :
126 : CONTAINS
127 : ! *****************************************************************************
128 : !> \brief Calculates the derivatives of the MO coefficients dC/dV^lambda_beta
129 : !> wrt to nuclear velocities. The derivative is indexed by `beta`, the
130 : !> electric dipole operator by `alpha`.
131 : !> Calculates the APT and AAT in velocity form
132 : !> P^lambda_alpha,beta = d< mu_alpha >/dV^lambda_beta
133 : !> M^lambda_alpha,beta = d< m_alpha >/dV^lambda_beta
134 : !> \param qs_env ...
135 : !> \param p_env ...
136 : !> \author Edward Ditler
137 : ! **************************************************************************************************
138 2 : SUBROUTINE vcd_linres(qs_env, p_env)
139 : TYPE(qs_environment_type), POINTER :: qs_env
140 : TYPE(qs_p_env_type) :: p_env
141 :
142 : INTEGER :: beta, i, latom
143 : LOGICAL :: mfp_is_done, mfp_repeat
144 60 : TYPE(vcd_env_type) :: vcd_env
145 :
146 2 : CALL cite_reference(Ditler2022)
147 :
148 : ! We need the position perturbation for the velocity perturbation operator
149 2 : CALL vcd_env_init(vcd_env, qs_env)
150 :
151 2 : mfp_repeat = vcd_env%distributed_origin
152 2 : mfp_is_done = .FALSE.
153 :
154 2 : qs_env%linres_control%linres_restart = .TRUE.
155 :
156 : ! Iterate over the list of atoms for which we want to calculate the APTs/AATs
157 : ! default is all atoms.
158 8 : DO latom = 1, SIZE(vcd_env%dcdr_env%list_of_atoms)
159 6 : vcd_env%dcdr_env%lambda = vcd_env%dcdr_env%list_of_atoms(latom)
160 :
161 6 : CALL prepare_per_atom(vcd_env%dcdr_env, qs_env)
162 6 : CALL prepare_per_atom_vcd(vcd_env, qs_env)
163 :
164 24 : DO beta = 1, 3 ! in every direction
165 :
166 18 : vcd_env%dcdr_env%beta = beta
167 18 : vcd_env%dcdr_env%deltaR(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) = 1._dp
168 :
169 : ! Since we do the heavy lifting anyways, we might also calculate the length form APTs here
170 18 : CALL dcdr_build_op_dR(vcd_env%dcdr_env, qs_env)
171 18 : CALL dcdr_response_dR(vcd_env%dcdr_env, p_env, qs_env)
172 18 : CALL apt_dR(qs_env, vcd_env%dcdr_env)
173 :
174 : ! And with the position perturbation ready, we can calculate the NVP
175 18 : CALL vcd_build_op_dV(vcd_env, qs_env)
176 18 : CALL vcd_response_dV(vcd_env, p_env, qs_env)
177 :
178 18 : CALL apt_dV(vcd_env, qs_env)
179 18 : CALL aat_dV(vcd_env, qs_env)
180 :
181 24 : IF (vcd_env%do_mfp) THEN
182 : ! Since we came so far, we might as well calculate the MFP AATs
183 : ! If we use a distributed origin we need to compute the MFP response again for each
184 : ! atom, because the reference point changes.
185 0 : IF (.NOT. mfp_is_done .OR. mfp_repeat) THEN
186 0 : DO i = 1, 3
187 0 : IF (vcd_env%origin_dependent_op_mfp) THEN
188 0 : CPWARN("Using the origin dependent MFP operator")
189 0 : CALL mfp_build_operator_gauge_dependent(vcd_env, qs_env, i)
190 : ELSE
191 0 : CALL mfp_build_operator_gauge_independent(vcd_env, qs_env, i)
192 : END IF
193 0 : CALL mfp_response(vcd_env, p_env, qs_env, i)
194 : END DO
195 : mfp_is_done = .TRUE.
196 : END IF
197 :
198 0 : CALL mfp_aat(vcd_env, qs_env)
199 : END IF
200 : END DO ! beta
201 :
202 : vcd_env%dcdr_env%apt_total_dcdr(:, :, vcd_env%dcdr_env%lambda) = &
203 : vcd_env%dcdr_env%apt_el_dcdr(:, :, vcd_env%dcdr_env%lambda) &
204 78 : + vcd_env%dcdr_env%apt_nuc_dcdr(:, :, vcd_env%dcdr_env%lambda)
205 :
206 : vcd_env%apt_total_nvpt(:, :, vcd_env%dcdr_env%lambda) = &
207 78 : vcd_env%apt_el_nvpt(:, :, vcd_env%dcdr_env%lambda) + vcd_env%apt_nuc_nvpt(:, :, vcd_env%dcdr_env%lambda)
208 :
209 8 : IF (vcd_env%do_mfp) THEN
210 0 : vcd_env%aat_atom_mfp(:, :, vcd_env%dcdr_env%lambda) = vcd_env%aat_atom_mfp(:, :, vcd_env%dcdr_env%lambda)*4._dp
211 : END IF
212 :
213 : END DO !lambda
214 :
215 2 : CALL vcd_print(vcd_env, qs_env)
216 2 : CALL vcd_env_cleanup(qs_env, vcd_env)
217 :
218 2 : END SUBROUTINE vcd_linres
219 :
220 : ! **************************************************************************************************
221 : !> \brief Calculates the derivatives of the MO coefficients dC/dR^lambda_beta
222 : !> wrt to nuclear coordinates. The derivative is index by `beta`, the
223 : !> electric dipole operator by `alpha`.
224 : !> Also calculates the APT
225 : !> P^lambda_alpha,beta = d< mu_alpha >/dR^lambda_beta
226 : !> and calculates the sum rules for the APT elements.
227 : !> \param qs_env ...
228 : !> \param p_env ...
229 : ! **************************************************************************************************
230 22 : SUBROUTINE dcdr_linres(qs_env, p_env)
231 : TYPE(qs_environment_type), POINTER :: qs_env
232 : TYPE(qs_p_env_type) :: p_env
233 :
234 : INTEGER :: beta, latom
235 308 : TYPE(dcdr_env_type) :: dcdr_env
236 : TYPE(polar_env_type), POINTER :: polar_env
237 :
238 22 : CALL cite_reference(Ditler2021)
239 22 : CALL dcdr_env_init(dcdr_env, qs_env)
240 :
241 22 : IF (.NOT. dcdr_env%z_matrix_method) THEN
242 :
243 72 : DO latom = 1, SIZE(dcdr_env%list_of_atoms)
244 54 : dcdr_env%lambda = dcdr_env%list_of_atoms(latom)
245 54 : CALL prepare_per_atom(dcdr_env, qs_env)
246 :
247 216 : DO beta = 1, 3 ! in every direction
248 162 : dcdr_env%beta = beta
249 162 : dcdr_env%deltaR(dcdr_env%beta, dcdr_env%lambda) = 1._dp
250 :
251 162 : CALL dcdr_build_op_dR(dcdr_env, qs_env)
252 162 : CALL dcdr_response_dR(dcdr_env, p_env, qs_env)
253 :
254 216 : IF (.NOT. dcdr_env%localized_psi0) THEN
255 126 : CALL apt_dR(qs_env, dcdr_env)
256 : ELSE IF (dcdr_env%localized_psi0) THEN
257 36 : CALL apt_dR_localization(qs_env, dcdr_env)
258 : END IF
259 :
260 : END DO !beta
261 :
262 : dcdr_env%apt_total_dcdr(:, :, dcdr_env%lambda) = &
263 720 : dcdr_env%apt_el_dcdr(:, :, dcdr_env%lambda) + dcdr_env%apt_nuc_dcdr(:, :, dcdr_env%lambda)
264 : END DO !lambda
265 :
266 : ELSE
267 :
268 4 : CALL polar_env_init(qs_env)
269 4 : CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
270 4 : CALL get_polar_env(polar_env=polar_env)
271 :
272 4 : IF (.NOT. dcdr_env%localized_psi0) THEN
273 4 : CALL polar_operators_local(qs_env)
274 : ELSE
275 0 : CALL polar_operators_local_wannier(qs_env, dcdr_env)
276 : END IF
277 :
278 4 : polar_env%do_periodic = .FALSE.
279 4 : CALL polar_response(p_env, qs_env)
280 :
281 16 : DO latom = 1, SIZE(dcdr_env%list_of_atoms)
282 12 : dcdr_env%lambda = dcdr_env%list_of_atoms(latom)
283 12 : CALL prepare_per_atom(dcdr_env, qs_env)
284 :
285 48 : DO beta = 1, 3 ! in every direction
286 36 : dcdr_env%beta = beta
287 36 : dcdr_env%deltaR(dcdr_env%beta, dcdr_env%lambda) = 1._dp
288 :
289 36 : CALL dcdr_build_op_dR(dcdr_env, qs_env)
290 48 : IF (.NOT. dcdr_env%localized_psi0) THEN
291 36 : CALL apt_dR(qs_env, dcdr_env)
292 : ELSE
293 0 : CALL apt_dR_localization(qs_env, dcdr_env)
294 : END IF
295 : END DO !beta
296 :
297 : dcdr_env%apt_total_dcdr(:, :, dcdr_env%lambda) = &
298 160 : dcdr_env%apt_el_dcdr(:, :, dcdr_env%lambda) + dcdr_env%apt_nuc_dcdr(:, :, dcdr_env%lambda)
299 : END DO !lambda
300 :
301 : END IF
302 :
303 22 : CALL dcdr_print(dcdr_env, qs_env)
304 22 : CALL dcdr_env_cleanup(qs_env, dcdr_env)
305 22 : END SUBROUTINE dcdr_linres
306 :
307 : ! **************************************************************************************************
308 : !> \brief Driver for the linear response calculatios
309 : !> \param force_env ...
310 : !> \par History
311 : !> 06.2005 created [MI]
312 : !> \author MI
313 : ! **************************************************************************************************
314 188 : SUBROUTINE linres_calculation(force_env)
315 :
316 : TYPE(force_env_type), POINTER :: force_env
317 :
318 : CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation'
319 :
320 : INTEGER :: handle
321 : TYPE(qs_environment_type), POINTER :: qs_env
322 :
323 188 : CALL timeset(routineN, handle)
324 :
325 188 : NULLIFY (qs_env)
326 :
327 188 : CPASSERT(ASSOCIATED(force_env))
328 188 : CPASSERT(force_env%ref_count > 0)
329 :
330 370 : SELECT CASE (force_env%in_use)
331 : CASE (use_qs_force)
332 182 : CALL force_env_get(force_env, qs_env=qs_env)
333 : CASE (use_qmmm)
334 6 : qs_env => force_env%qmmm_env%qs_env
335 : CASE DEFAULT
336 188 : CPABORT("Does not recognize this force_env")
337 : END SELECT
338 :
339 188 : qs_env%linres_run = .TRUE.
340 :
341 188 : CALL linres_calculation_low(qs_env)
342 :
343 188 : CALL timestop(handle)
344 :
345 188 : END SUBROUTINE linres_calculation
346 :
347 : ! **************************************************************************************************
348 : !> \brief Linear response can be called as run type or as post scf calculation
349 : !> Initialize the perturbation environment
350 : !> Define which properties is to be calculated
351 : !> Start up the optimization of the response density and wfn
352 : !> \param qs_env ...
353 : !> \par History
354 : !> 06.2005 created [MI]
355 : !> 02.2013 added polarizability section [SL]
356 : !> \author MI
357 : ! **************************************************************************************************
358 27581 : SUBROUTINE linres_calculation_low(qs_env)
359 :
360 : TYPE(qs_environment_type), POINTER :: qs_env
361 :
362 : CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation_low'
363 :
364 : INTEGER :: every_n_step, handle, iounit
365 : LOGICAL :: dcdr_present, do_apt_fd, epr_present, &
366 : issc_present, lr_calculation, &
367 : nmr_present, polar_present, vcd_present
368 : TYPE(cp_logger_type), POINTER :: logger
369 : TYPE(dft_control_type), POINTER :: dft_control
370 : TYPE(linres_control_type), POINTER :: linres_control
371 : TYPE(qs_p_env_type) :: p_env
372 : TYPE(section_vals_type), POINTER :: lr_section, prop_section, xc_fun_section
373 :
374 27581 : CALL timeset(routineN, handle)
375 :
376 : lr_calculation = .FALSE.
377 : nmr_present = .FALSE.
378 : epr_present = .FALSE.
379 : issc_present = .FALSE.
380 : polar_present = .FALSE.
381 : dcdr_present = .FALSE.
382 : do_apt_fd = .FALSE.
383 :
384 27581 : NULLIFY (dft_control, linres_control, logger, prop_section, lr_section, xc_fun_section)
385 :
386 27581 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
387 27581 : CALL section_vals_get(lr_section, explicit=lr_calculation)
388 :
389 27581 : CALL section_vals_val_get(lr_section, "DCDR%APT_FD", l_val=do_apt_fd)
390 27581 : IF (do_apt_fd) THEN
391 28 : CALL timestop(handle)
392 28 : RETURN
393 : END IF
394 :
395 27553 : logger => cp_get_default_logger()
396 :
397 27553 : CALL section_vals_val_get(lr_section, "EVERY_N_STEP", i_val=every_n_step)
398 :
399 27553 : IF (lr_calculation .AND. MODULO(qs_env%sim_step, every_n_step) == 0) THEN
400 336 : prop_section => section_vals_get_subs_vals(lr_section, "DCDR")
401 336 : CALL section_vals_get(prop_section, explicit=dcdr_present)
402 336 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
403 336 : xc_fun_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL")
404 336 : IF (dcdr_present .AND. xc_uses_kinetic_energy_density(xc_fun_section, dft_control%lsd)) THEN
405 0 : CPABORT(dcdr_meta_gga_error)
406 : END IF
407 :
408 336 : CALL linres_init(lr_section, p_env, qs_env)
409 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
410 336 : extension=".linresLog")
411 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
412 336 : linres_control=linres_control)
413 :
414 : ! The type of perturbation has not been defined yet
415 336 : linres_control%property = lr_none
416 :
417 : ! We do NMR or EPR, then compute the current response
418 336 : prop_section => section_vals_get_subs_vals(lr_section, "NMR")
419 336 : CALL section_vals_get(prop_section, explicit=nmr_present)
420 336 : prop_section => section_vals_get_subs_vals(lr_section, "EPR")
421 336 : CALL section_vals_get(prop_section, explicit=epr_present)
422 :
423 336 : IF (nmr_present .OR. epr_present) THEN
424 : CALL nmr_epr_linres(linres_control, qs_env, p_env, dft_control, &
425 174 : nmr_present, epr_present, iounit)
426 : END IF
427 :
428 : ! We do the indirect spin-spin coupling calculation
429 336 : prop_section => section_vals_get_subs_vals(lr_section, "SPINSPIN")
430 336 : CALL section_vals_get(prop_section, explicit=issc_present)
431 :
432 336 : IF (issc_present) THEN
433 12 : CALL issc_linres(linres_control, qs_env, p_env, dft_control)
434 : END IF
435 :
436 : ! We do the polarizability calculation
437 336 : prop_section => section_vals_get_subs_vals(lr_section, "POLAR")
438 336 : CALL section_vals_get(prop_section, explicit=polar_present)
439 336 : IF (polar_present) THEN
440 126 : CALL polar_linres(qs_env, p_env)
441 : END IF
442 :
443 : ! Nuclear Position Perturbation
444 336 : IF (dcdr_present) THEN
445 22 : CALL dcdr_linres(qs_env, p_env)
446 : END IF
447 :
448 : ! VCD
449 336 : prop_section => section_vals_get_subs_vals(lr_section, "VCD")
450 336 : CALL section_vals_get(prop_section, explicit=vcd_present)
451 :
452 336 : IF (vcd_present) THEN
453 2 : CALL vcd_linres(qs_env, p_env)
454 : END IF
455 :
456 : ! Other possible LR calculations can be introduced here
457 :
458 336 : CALL p_env_release(p_env)
459 :
460 336 : IF (iounit > 0) THEN
461 : WRITE (UNIT=iounit, FMT="(/,T2,A,/,T25,A,/,T2,A,/)") &
462 168 : REPEAT("=", 79), &
463 168 : "ENDED LINRES CALCULATION", &
464 336 : REPEAT("=", 79)
465 : END IF
466 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
467 336 : "PRINT%PROGRAM_RUN_INFO")
468 : END IF
469 :
470 27553 : CALL timestop(handle)
471 :
472 137905 : END SUBROUTINE linres_calculation_low
473 :
474 : ! **************************************************************************************************
475 : !> \brief Initialize some general settings like the p_env
476 : !> Localize the psi0 if required
477 : !> \param lr_section ...
478 : !> \param p_env ...
479 : !> \param qs_env ...
480 : !> \par History
481 : !> 06.2005 created [MI]
482 : !> \author MI
483 : !> \note
484 : !> - The localization should probably be always for all the occupied states
485 : ! **************************************************************************************************
486 2016 : SUBROUTINE linres_init(lr_section, p_env, qs_env)
487 :
488 : TYPE(section_vals_type), POINTER :: lr_section
489 : TYPE(qs_p_env_type), INTENT(OUT) :: p_env
490 : TYPE(qs_environment_type), POINTER :: qs_env
491 :
492 : INTEGER :: iounit, ispin
493 : LOGICAL :: do_it
494 : TYPE(cp_logger_type), POINTER :: logger
495 336 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, rho_ao
496 : TYPE(dft_control_type), POINTER :: dft_control
497 : TYPE(linres_control_type), POINTER :: linres_control
498 336 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
499 : TYPE(qs_rho_type), POINTER :: rho
500 : TYPE(section_vals_type), POINTER :: loc_section
501 :
502 336 : NULLIFY (logger)
503 336 : logger => cp_get_default_logger()
504 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
505 336 : extension=".linresLog")
506 336 : NULLIFY (dft_control, linres_control, loc_section, rho, mos, matrix_ks, rho_ao)
507 :
508 336 : ALLOCATE (linres_control)
509 336 : CALL set_qs_env(qs_env=qs_env, linres_control=linres_control)
510 : CALL get_qs_env(qs_env=qs_env, &
511 336 : dft_control=dft_control, matrix_ks=matrix_ks, mos=mos, rho=rho)
512 336 : CALL qs_rho_get(rho, rho_ao=rho_ao)
513 :
514 : ! Localized Psi0 are required when the position operator has to be defined (nmr)
515 336 : loc_section => section_vals_get_subs_vals(lr_section, "LOCALIZE")
516 : CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", &
517 336 : l_val=linres_control%localized_psi0)
518 336 : IF (linres_control%localized_psi0) THEN
519 190 : IF (iounit > 0) THEN
520 : WRITE (UNIT=iounit, FMT="(/,T3,A,A)") &
521 95 : "Localization of ground state orbitals", &
522 190 : " before starting linear response calculation"
523 : END IF
524 :
525 190 : CALL linres_localize(qs_env, linres_control, dft_control%nspins)
526 :
527 458 : DO ispin = 1, dft_control%nspins
528 458 : CALL calculate_density_matrix(mos(ispin), rho_ao(ispin)%matrix)
529 : END DO
530 : ! ** update qs_env%rho
531 190 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
532 : END IF
533 :
534 336 : CALL section_vals_val_get(lr_section, "RESTART", l_val=linres_control%linres_restart)
535 336 : CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter)
536 336 : CALL section_vals_val_get(lr_section, "EPS", r_val=linres_control%eps)
537 336 : CALL section_vals_val_get(lr_section, "EPS_FILTER", r_val=linres_control%eps_filter)
538 336 : CALL section_vals_val_get(lr_section, "RESTART_EVERY", i_val=linres_control%restart_every)
539 336 : CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
540 336 : CALL section_vals_val_get(lr_section, "ENERGY_GAP", r_val=linres_control%energy_gap)
541 :
542 336 : IF (iounit > 0) THEN
543 : WRITE (UNIT=iounit, FMT="(/,T2,A,/,T25,A,/,T2,A,/)") &
544 168 : REPEAT("=", 79), &
545 168 : "START LINRES CALCULATION", &
546 336 : REPEAT("=", 79)
547 :
548 : WRITE (UNIT=iounit, FMT="(T2,A)") &
549 168 : "LINRES| Properties to be calculated:"
550 168 : CALL section_vals_val_get(lr_section, "NMR%_SECTION_PARAMETERS_", l_val=do_it)
551 168 : IF (do_it) WRITE (UNIT=iounit, FMT="(T62,A)") "NMR Chemical Shift"
552 168 : CALL section_vals_val_get(lr_section, "EPR%_SECTION_PARAMETERS_", l_val=do_it)
553 168 : IF (do_it) WRITE (UNIT=iounit, FMT="(T68,A)") "EPR g Tensor"
554 168 : CALL section_vals_val_get(lr_section, "SPINSPIN%_SECTION_PARAMETERS_", l_val=do_it)
555 168 : IF (do_it) WRITE (UNIT=iounit, FMT="(T43,A)") "Indirect spin-spin coupling constants"
556 168 : CALL section_vals_val_get(lr_section, "POLAR%_SECTION_PARAMETERS_", l_val=do_it)
557 168 : IF (do_it) WRITE (UNIT=iounit, FMT="(T57,A)") "Electric Polarizability"
558 :
559 168 : IF (linres_control%localized_psi0) WRITE (UNIT=iounit, FMT="(T2,A,T65,A)") &
560 95 : "LINRES|", " LOCALIZED PSI0"
561 :
562 : WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
563 168 : "LINRES| Optimization algorithm", " Conjugate Gradients"
564 :
565 169 : SELECT CASE (linres_control%preconditioner_type)
566 : CASE (ot_precond_none)
567 : WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
568 1 : "LINRES| Preconditioner", " NONE"
569 : CASE (ot_precond_full_single)
570 : WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
571 2 : "LINRES| Preconditioner", " FULL_SINGLE"
572 : CASE (ot_precond_full_kinetic)
573 : WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
574 3 : "LINRES| Preconditioner", " FULL_KINETIC"
575 : CASE (ot_precond_s_inverse)
576 : WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
577 12 : "LINRES| Preconditioner", " FULL_S_INVERSE"
578 : CASE (ot_precond_full_single_inverse)
579 : WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
580 32 : "LINRES| Preconditioner", " FULL_SINGLE_INVERSE"
581 : CASE (ot_precond_full_all)
582 : WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") &
583 118 : "LINRES| Preconditioner", " FULL_ALL"
584 : CASE DEFAULT
585 168 : CPABORT("Preconditioner NYI")
586 : END SELECT
587 :
588 : WRITE (UNIT=iounit, FMT="(T2,A,T72,ES8.1)") &
589 168 : "LINRES| EPS", linres_control%eps
590 : WRITE (UNIT=iounit, FMT="(T2,A,T72,I8)") &
591 168 : "LINRES| MAX_ITER", linres_control%max_iter
592 : END IF
593 :
594 : !------------------!
595 : ! create the p_env !
596 : !------------------!
597 336 : CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.TRUE., linres_control=linres_control)
598 :
599 : ! update the m_epsilon matrix
600 336 : CALL p_env_psi0_changed(p_env, qs_env)
601 :
602 336 : p_env%new_preconditioner = .TRUE.
603 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
604 336 : "PRINT%PROGRAM_RUN_INFO")
605 :
606 336 : END SUBROUTINE linres_init
607 :
608 : ! **************************************************************************************************
609 : !> \brief ...
610 : !> \param linres_control ...
611 : !> \param qs_env ...
612 : !> \param p_env ...
613 : !> \param dft_control ...
614 : !> \param nmr_present ...
615 : !> \param epr_present ...
616 : !> \param iounit ...
617 : ! **************************************************************************************************
618 174 : SUBROUTINE nmr_epr_linres(linres_control, qs_env, p_env, dft_control, nmr_present, epr_present, iounit)
619 :
620 : TYPE(linres_control_type), POINTER :: linres_control
621 : TYPE(qs_environment_type), POINTER :: qs_env
622 : TYPE(qs_p_env_type) :: p_env
623 : TYPE(dft_control_type), POINTER :: dft_control
624 : LOGICAL :: nmr_present, epr_present
625 : INTEGER :: iounit
626 :
627 : INTEGER :: iB
628 : LOGICAL :: do_qmmm
629 : TYPE(current_env_type) :: current_env
630 : TYPE(epr_env_type) :: epr_env
631 : TYPE(nmr_env_type) :: nmr_env
632 :
633 174 : linres_control%property = lr_current
634 :
635 174 : CALL cite_reference(Weber2009)
636 :
637 174 : IF (.NOT. linres_control%localized_psi0) THEN
638 : CALL cp_abort(__LOCATION__, &
639 : "Are you sure that you want to calculate the chemical "// &
640 0 : "shift without localized psi0?")
641 : CALL linres_localize(qs_env, linres_control, &
642 0 : dft_control%nspins, centers_only=.TRUE.)
643 : END IF
644 174 : IF (dft_control%nspins /= 2 .AND. epr_present) THEN
645 0 : CPABORT("LSD is needed to perform a g tensor calculation!")
646 : END IF
647 : !
648 : !Initialize the current environment
649 174 : do_qmmm = .FALSE.
650 174 : IF (qs_env%qmmm) do_qmmm = .TRUE.
651 174 : current_env%do_qmmm = do_qmmm
652 : !current_env%prop='nmr'
653 174 : CALL current_env_init(current_env, qs_env)
654 174 : CALL current_operators(current_env, qs_env)
655 174 : CALL current_response(current_env, p_env, qs_env)
656 : !
657 174 : IF (current_env%all_pert_op_done) THEN
658 : !Initialize the nmr environment
659 174 : IF (nmr_present) THEN
660 160 : CALL nmr_env_init(nmr_env, qs_env)
661 : END IF
662 : !
663 : !Initialize the epr environment
664 174 : IF (epr_present) THEN
665 14 : CALL epr_env_init(epr_env, qs_env)
666 14 : CALL epr_g_zke(epr_env, qs_env)
667 14 : CALL epr_nablavks(epr_env, qs_env)
668 : END IF
669 : !
670 : ! Build the rs_gauge if needed
671 : !CALL current_set_gauge(current_env,qs_env)
672 : !
673 : ! Loop over field direction
674 696 : DO iB = 1, 3
675 : !
676 : ! Build current response and succeptibility
677 522 : CALL current_build_current(current_env, qs_env, iB)
678 522 : CALL current_build_chi(current_env, qs_env, iB)
679 : !
680 : ! Compute NMR shift
681 522 : IF (nmr_present) THEN
682 480 : CALL nmr_shift(nmr_env, current_env, qs_env, iB)
683 : END IF
684 : !
685 : ! Compute EPR
686 696 : IF (epr_present) THEN
687 42 : CALL epr_ind_magnetic_field(epr_env, current_env, qs_env, iB)
688 42 : CALL epr_g_so(epr_env, current_env, qs_env, iB)
689 42 : CALL epr_g_soo(epr_env, current_env, qs_env, iB)
690 : END IF
691 : END DO
692 : !
693 : ! Finalized the nmr environment
694 174 : IF (nmr_present) THEN
695 160 : CALL nmr_shift_print(nmr_env, current_env, qs_env)
696 160 : CALL nmr_env_cleanup(nmr_env)
697 : END IF
698 : !
699 : ! Finalized the epr environment
700 174 : IF (epr_present) THEN
701 14 : CALL epr_g_print(epr_env, qs_env)
702 14 : CALL epr_env_cleanup(epr_env)
703 : END IF
704 : !
705 : ELSE
706 0 : IF (iounit > 0) THEN
707 : WRITE (iounit, "(T10,A,/T20,A,/)") &
708 0 : "CURRENT: Not all responses to perturbation operators could be calculated.", &
709 0 : " Hence: NO nmr and NO epr possible."
710 : END IF
711 : END IF
712 : ! Finalized the current environment
713 174 : CALL current_env_cleanup(current_env)
714 :
715 12702 : END SUBROUTINE nmr_epr_linres
716 :
717 : ! **************************************************************************************************
718 : !> \brief ...
719 : !> \param linres_control ...
720 : !> \param qs_env ...
721 : !> \param p_env ...
722 : !> \param dft_control ...
723 : ! **************************************************************************************************
724 12 : SUBROUTINE issc_linres(linres_control, qs_env, p_env, dft_control)
725 :
726 : TYPE(linres_control_type), POINTER :: linres_control
727 : TYPE(qs_environment_type), POINTER :: qs_env
728 : TYPE(qs_p_env_type) :: p_env
729 : TYPE(dft_control_type), POINTER :: dft_control
730 :
731 : INTEGER :: iatom
732 : LOGICAL :: do_qmmm
733 : TYPE(current_env_type) :: current_env
734 : TYPE(issc_env_type) :: issc_env
735 :
736 12 : linres_control%property = lr_current
737 12 : IF (.NOT. linres_control%localized_psi0) THEN
738 : CALL cp_abort(__LOCATION__, &
739 : "Are you sure that you want to calculate the chemical "// &
740 0 : "shift without localized psi0?")
741 : CALL linres_localize(qs_env, linres_control, &
742 0 : dft_control%nspins, centers_only=.TRUE.)
743 : END IF
744 : !
745 : !Initialize the current environment
746 : do_qmmm = .FALSE.
747 : IF (qs_env%qmmm) do_qmmm = .TRUE.
748 12 : current_env%do_qmmm = do_qmmm
749 : !current_env%prop='issc'
750 : !CALL current_env_init(current_env,qs_env)
751 : !CALL current_response(current_env,p_env,qs_env)
752 : !
753 : !Initialize the issc environment
754 12 : CALL issc_env_init(issc_env, qs_env)
755 : !
756 : ! Loop over atoms
757 56 : DO iatom = 1, issc_env%issc_natms
758 44 : CALL issc_operators(issc_env, qs_env, iatom)
759 44 : CALL issc_response(issc_env, p_env, qs_env)
760 56 : CALL issc_issc(issc_env, qs_env, iatom)
761 : END DO
762 : !
763 : ! Finalized the issc environment
764 12 : CALL issc_print(issc_env, qs_env)
765 12 : CALL issc_env_cleanup(issc_env)
766 :
767 888 : END SUBROUTINE issc_linres
768 :
769 : ! **************************************************************************************************
770 : !> \brief ...
771 : !> \param qs_env ...
772 : !> \param p_env ...
773 : !> \par History
774 : !> 06.2018 polar_env integrated into qs_env (MK)
775 : ! **************************************************************************************************
776 126 : SUBROUTINE polar_linres(qs_env, p_env)
777 :
778 : TYPE(qs_environment_type), POINTER :: qs_env
779 : TYPE(qs_p_env_type) :: p_env
780 :
781 126 : CALL polar_env_init(qs_env)
782 126 : CALL polar_operators(qs_env)
783 126 : CALL polar_response(p_env, qs_env)
784 126 : CALL polar_polar(qs_env)
785 126 : CALL polar_print(qs_env)
786 :
787 126 : END SUBROUTINE polar_linres
788 :
789 : END MODULE qs_linres_module
|
