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