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