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 Polarizability calculation by dfpt
10 : !> Initialization of the polar_env,
11 : !> Perturbation Hamiltonian by application of the Berry phase operator to psi0
12 : !> Write output
13 : !> Deallocate everything
14 : !> periodic Raman SL February 2013
15 : !> \note
16 : ! **************************************************************************************************
17 : MODULE qs_linres_polar_utils
18 : USE bibliography, ONLY: Luber2014,&
19 : cite_reference
20 : USE cell_types, ONLY: cell_type
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
23 : USE cp_fm_basic_linalg, ONLY: cp_fm_trace
24 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
25 : cp_fm_struct_release,&
26 : cp_fm_struct_type
27 : USE cp_fm_types, ONLY: cp_fm_create,&
28 : cp_fm_get_info,&
29 : cp_fm_release,&
30 : cp_fm_set_all,&
31 : cp_fm_to_fm,&
32 : cp_fm_type
33 : USE cp_log_handling, ONLY: cp_get_default_logger,&
34 : cp_logger_get_default_io_unit,&
35 : cp_logger_type
36 : USE cp_output_handling, ONLY: cp_p_file,&
37 : cp_print_key_finished_output,&
38 : cp_print_key_should_output,&
39 : cp_print_key_unit_nr
40 : USE cp_result_methods, ONLY: cp_results_erase,&
41 : put_results
42 : USE cp_result_types, ONLY: cp_result_type
43 : USE force_env_types, ONLY: force_env_get,&
44 : force_env_type
45 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
46 : section_vals_type,&
47 : section_vals_val_get
48 : USE kinds, ONLY: default_string_length,&
49 : dp
50 : USE machine, ONLY: m_flush
51 : USE mathconstants, ONLY: twopi
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE physcon, ONLY: angstrom
54 : USE qs_environment_types, ONLY: get_qs_env,&
55 : qs_environment_type,&
56 : set_qs_env
57 : USE qs_linres_methods, ONLY: linres_read_restart,&
58 : linres_solver,&
59 : linres_write_restart
60 : USE qs_linres_types, ONLY: get_polar_env,&
61 : linres_control_type,&
62 : polar_env_type,&
63 : set_polar_env
64 : USE qs_matrix_pools, ONLY: qs_matrix_pools_type
65 : USE qs_mo_types, ONLY: get_mo_set,&
66 : mo_set_type
67 : USE qs_p_env_types, ONLY: qs_p_env_type
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 :
74 : PUBLIC :: polar_env_init, polar_polar, polar_print, polar_response, write_polarisability_tensor
75 :
76 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_polar_utils'
77 :
78 : CONTAINS
79 :
80 : ! **************************************************************************************************
81 : !> \brief Initialize the polar environment
82 : !> \param qs_env ...
83 : !> \par History
84 : !> 06.2018 polar_env integrated into qs_env (MK)
85 : ! **************************************************************************************************
86 112 : SUBROUTINE polar_env_init(qs_env)
87 :
88 : TYPE(qs_environment_type), POINTER :: qs_env
89 :
90 : CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_env_init'
91 :
92 : INTEGER :: handle, idir, iounit, ispin, m, nao, &
93 : nmo, nspins
94 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
95 : TYPE(cp_fm_type), POINTER :: mo_coeff
96 : TYPE(cp_logger_type), POINTER :: logger
97 112 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
98 : TYPE(dft_control_type), POINTER :: dft_control
99 : TYPE(linres_control_type), POINTER :: linres_control
100 112 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
101 : TYPE(polar_env_type), POINTER :: polar_env
102 : TYPE(section_vals_type), POINTER :: lr_section, polar_section
103 :
104 112 : CALL timeset(routineN, handle)
105 :
106 112 : NULLIFY (dft_control)
107 112 : NULLIFY (linres_control)
108 112 : NULLIFY (logger)
109 112 : NULLIFY (matrix_s)
110 112 : NULLIFY (mos)
111 112 : NULLIFY (polar_env)
112 112 : NULLIFY (lr_section, polar_section)
113 :
114 112 : logger => cp_get_default_logger()
115 112 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
116 :
117 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
118 112 : extension=".linresLog")
119 :
120 112 : IF (iounit > 0) THEN
121 56 : WRITE (iounit, "(/,(T2,A))") "POLAR| Starting polarizability calculation", &
122 112 : "POLAR| Initialization of the polar environment"
123 : END IF
124 :
125 : polar_section => section_vals_get_subs_vals(qs_env%input, &
126 112 : "PROPERTIES%LINRES%POLAR")
127 :
128 : CALL get_qs_env(qs_env=qs_env, &
129 : polar_env=polar_env, &
130 : dft_control=dft_control, &
131 : matrix_s=matrix_s, &
132 : linres_control=linres_control, &
133 112 : mos=mos)
134 :
135 : ! Create polar environment if needed
136 112 : IF (.NOT. ASSOCIATED(polar_env)) THEN
137 84 : ALLOCATE (polar_env)
138 84 : CALL set_qs_env(qs_env=qs_env, polar_env=polar_env)
139 : END IF
140 :
141 112 : nspins = dft_control%nspins
142 :
143 112 : CALL section_vals_val_get(polar_section, "DO_RAMAN", l_val=polar_env%do_raman)
144 112 : CALL section_vals_val_get(polar_section, "PERIODIC_DIPOLE_OPERATOR", l_val=polar_env%do_periodic)
145 :
146 : ! Allocate components of the polar environment if needed
147 112 : IF (.NOT. ASSOCIATED(polar_env%polar)) THEN
148 84 : ALLOCATE (polar_env%polar(3, 3))
149 1092 : polar_env%polar(:, :) = 0.0_dp
150 : END IF
151 112 : IF (.NOT. ASSOCIATED(polar_env%dBerry_psi0)) THEN
152 620 : ALLOCATE (polar_env%dBerry_psi0(3, nspins))
153 : ELSE
154 : ! Remove previous matrices
155 56 : DO ispin = 1, nspins
156 140 : DO idir = 1, 3
157 112 : CALL cp_fm_release(polar_env%dBerry_psi0(idir, ispin))
158 : END DO
159 : END DO
160 : END IF
161 112 : IF (.NOT. ASSOCIATED(polar_env%psi1_dBerry)) THEN
162 620 : ALLOCATE (polar_env%psi1_dBerry(3, nspins))
163 : ELSE
164 : ! Remove previous matrices
165 56 : DO ispin = 1, nspins
166 140 : DO idir = 1, 3
167 112 : CALL cp_fm_release(polar_env%psi1_dBerry(idir, ispin))
168 : END DO
169 : END DO
170 : END IF
171 232 : DO ispin = 1, nspins
172 120 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
173 120 : CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao)
174 120 : NULLIFY (tmp_fm_struct)
175 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
176 : ncol_global=m, &
177 120 : context=mo_coeff%matrix_struct%context)
178 480 : DO idir = 1, 3
179 360 : CALL cp_fm_create(polar_env%dBerry_psi0(idir, ispin), tmp_fm_struct)
180 480 : CALL cp_fm_create(polar_env%psi1_dBerry(idir, ispin), tmp_fm_struct)
181 : END DO
182 352 : CALL cp_fm_struct_release(tmp_fm_struct)
183 : END DO
184 :
185 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
186 112 : "PRINT%PROGRAM_RUN_INFO")
187 :
188 112 : CALL timestop(handle)
189 :
190 112 : END SUBROUTINE polar_env_init
191 :
192 : ! **************************************************************************************************
193 : !> \brief ...
194 : !> \param qs_env ...
195 : !> \par History
196 : !> 06.2018 polar_env integrated into qs_env (MK)
197 : ! **************************************************************************************************
198 108 : SUBROUTINE polar_polar(qs_env)
199 :
200 : TYPE(qs_environment_type), POINTER :: qs_env
201 :
202 : CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_polar'
203 :
204 : INTEGER :: handle, i, iounit, ispin, nspins, z
205 : LOGICAL :: do_periodic, do_raman, run_stopped
206 : REAL(dp) :: ptmp
207 108 : REAL(dp), DIMENSION(:, :), POINTER :: polar, polar_tmp
208 : TYPE(cell_type), POINTER :: cell
209 108 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry
210 : TYPE(cp_logger_type), POINTER :: logger
211 : TYPE(dft_control_type), POINTER :: dft_control
212 108 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
213 : TYPE(polar_env_type), POINTER :: polar_env
214 :
215 108 : CALL timeset(routineN, handle)
216 :
217 108 : NULLIFY (cell, dft_control, polar, psi1_dBerry, logger)
218 108 : NULLIFY (mos, dBerry_psi0)
219 108 : logger => cp_get_default_logger()
220 108 : iounit = cp_logger_get_default_io_unit(logger)
221 :
222 : CALL get_qs_env(qs_env=qs_env, &
223 : cell=cell, &
224 : dft_control=dft_control, &
225 : mos=mos, &
226 108 : polar_env=polar_env)
227 :
228 108 : nspins = dft_control%nspins
229 :
230 : CALL get_polar_env(polar_env=polar_env, &
231 : do_raman=do_raman, &
232 108 : run_stopped=run_stopped)
233 :
234 108 : IF (.NOT. run_stopped .AND. do_raman) THEN
235 :
236 108 : CALL cite_reference(Luber2014)
237 :
238 : CALL get_polar_env(polar_env=polar_env, &
239 : do_periodic=do_periodic, &
240 : dBerry_psi0=dBerry_psi0, &
241 : polar=polar, &
242 108 : psi1_dBerry=psi1_dBerry)
243 :
244 : ! Initialize
245 108 : ALLOCATE (polar_tmp(3, 3))
246 1404 : polar_tmp(:, :) = 0.0_dp
247 :
248 432 : DO i = 1, 3 ! directions of electric field
249 1404 : DO z = 1, 3 !dipole directions
250 2322 : DO ispin = 1, dft_control%nspins
251 : !SL compute trace
252 : ptmp = 0.0_dp
253 1026 : CALL cp_fm_trace(psi1_dBerry(i, ispin), dBerry_psi0(z, ispin), ptmp)
254 1998 : polar_tmp(i, z) = polar_tmp(i, z) - 2.0_dp*ptmp
255 : END DO
256 : END DO
257 : END DO !spin
258 :
259 108 : IF (do_periodic) THEN
260 2124 : polar(:, :) = MATMUL(MATMUL(cell%hmat, polar_tmp), TRANSPOSE(cell%hmat))/(twopi*twopi)
261 : ELSE
262 2340 : polar(:, :) = polar_tmp(:, :)
263 : END IF
264 : !SL evtl maxocc instead?
265 108 : IF (dft_control%nspins == 1) THEN
266 1326 : polar(:, :) = 2.0_dp*polar(:, :)
267 : END IF
268 :
269 108 : IF (ASSOCIATED(polar_tmp)) THEN
270 108 : DEALLOCATE (polar_tmp)
271 : END IF
272 :
273 : END IF ! do_raman
274 :
275 108 : CALL timestop(handle)
276 :
277 108 : END SUBROUTINE polar_polar
278 :
279 : ! **************************************************************************************************
280 : !> \brief Print information related to the polarisability tensor
281 : !> \param qs_env ...
282 : !> \par History
283 : !> 06.2018 polar_env integrated into qs_env (MK)
284 : ! **************************************************************************************************
285 108 : SUBROUTINE polar_print(qs_env)
286 :
287 : TYPE(qs_environment_type), POINTER :: qs_env
288 :
289 : CHARACTER(LEN=default_string_length) :: description
290 : INTEGER :: iounit, unit_p
291 : LOGICAL :: do_raman, run_stopped
292 108 : REAL(dp), DIMENSION(:, :), POINTER :: polar
293 : TYPE(cp_logger_type), POINTER :: logger
294 : TYPE(cp_result_type), POINTER :: results
295 : TYPE(dft_control_type), POINTER :: dft_control
296 : TYPE(mp_para_env_type), POINTER :: para_env
297 : TYPE(polar_env_type), POINTER :: polar_env
298 : TYPE(section_vals_type), POINTER :: polar_section
299 :
300 108 : NULLIFY (logger, dft_control, para_env, results)
301 :
302 : CALL get_qs_env(qs_env=qs_env, &
303 : dft_control=dft_control, &
304 : polar_env=polar_env, &
305 : results=results, &
306 108 : para_env=para_env)
307 :
308 108 : logger => cp_get_default_logger()
309 108 : iounit = cp_logger_get_default_io_unit(logger)
310 :
311 108 : polar_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%POLAR")
312 :
313 : CALL get_polar_env(polar_env=polar_env, &
314 : polar=polar, &
315 : do_raman=do_raman, &
316 108 : run_stopped=run_stopped)
317 :
318 108 : IF (.NOT. run_stopped .AND. do_raman) THEN
319 :
320 108 : description = "[POLAR]"
321 108 : CALL cp_results_erase(results, description=description)
322 108 : CALL put_results(results, description=description, values=polar(:, :))
323 :
324 108 : IF (BTEST(cp_print_key_should_output(logger%iter_info, polar_section, &
325 : "PRINT%POLAR_MATRIX"), cp_p_file)) THEN
326 :
327 : unit_p = cp_print_key_unit_nr(logger, polar_section, "PRINT%POLAR_MATRIX", &
328 102 : extension=".data", middle_name="raman", log_filename=.FALSE.)
329 102 : IF (unit_p > 0) THEN
330 51 : IF (unit_p /= iounit) THEN
331 51 : WRITE (unit_p, *)
332 51 : WRITE (unit_p, '(T10,A)') 'POLARIZABILITY TENSOR (atomic units):'
333 51 : WRITE (unit_p, '(T10,A,3F15.5)') "xx,yy,zz", polar(1, 1), polar(2, 2), polar(3, 3)
334 51 : WRITE (unit_p, '(T10,A,3F15.5)') "xy,xz,yz", polar(1, 2), polar(1, 3), polar(2, 3)
335 51 : WRITE (unit_p, '(T10,A,3F15.5)') "yx,zx,zy", polar(2, 1), polar(3, 1), polar(3, 2)
336 51 : WRITE (unit_p, '(T10,A)') 'POLARIZABILITY TENSOR (Angstrom^3):'
337 51 : WRITE (unit_p, '(T10,A,3F15.5)') "xx,yy,zz", polar(1, 1)*angstrom**3, &
338 102 : polar(2, 2)*angstrom**3, polar(3, 3)*angstrom**3
339 51 : WRITE (unit_p, '(T10,A,3F15.5)') "xy,xz,yz", polar(1, 2)*angstrom**3, &
340 102 : polar(1, 3)*angstrom**3, polar(2, 3)*angstrom**3
341 51 : WRITE (unit_p, '(T10,A,3F15.5)') "yx,zx,zy", polar(2, 1)*angstrom**3, &
342 102 : polar(3, 1)*angstrom**3, polar(3, 2)*angstrom**3
343 : CALL cp_print_key_finished_output(unit_p, logger, polar_section, &
344 51 : "PRINT%POLAR_MATRIX")
345 : END IF
346 : END IF
347 : END IF
348 108 : IF (iounit > 0) THEN
349 : WRITE (iounit, '(/,T2,A)') &
350 54 : 'POLAR| Polarizability tensor [a.u.]'
351 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
352 54 : 'POLAR| xx,yy,zz', polar(1, 1), polar(2, 2), polar(3, 3)
353 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
354 54 : 'POLAR| xy,xz,yz', polar(1, 2), polar(1, 3), polar(2, 3)
355 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
356 54 : 'POLAR| yx,zx,zy', polar(2, 1), polar(3, 1), polar(3, 2)
357 : WRITE (iounit, '(/,T2,A)') &
358 54 : 'POLAR| Polarizability tensor [ang^3]'
359 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
360 54 : 'POLAR| xx,yy,zz', polar(1, 1)*angstrom**3, polar(2, 2)*angstrom**3, polar(3, 3)*angstrom**3
361 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
362 54 : 'POLAR| xy,xz,yz', polar(1, 2)*angstrom**3, polar(1, 3)*angstrom**3, polar(2, 3)*angstrom**3
363 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
364 54 : 'POLAR| yx,zx,zy', polar(2, 1)*angstrom**3, polar(3, 1)*angstrom**3, polar(3, 2)*angstrom**3
365 : END IF
366 : END IF
367 :
368 108 : END SUBROUTINE polar_print
369 :
370 : ! **************************************************************************************************
371 : !> \brief Calculate the polarisability tensor using response theory
372 : !> \param p_env ...
373 : !> \param qs_env ...
374 : !> \par History
375 : !> 06.2018 polar_env integrated into qs_env (MK)
376 : ! **************************************************************************************************
377 112 : SUBROUTINE polar_response(p_env, qs_env)
378 :
379 : TYPE(qs_p_env_type) :: p_env
380 : TYPE(qs_environment_type), POINTER :: qs_env
381 :
382 : CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_response'
383 :
384 : INTEGER :: handle, idir, iounit, ispin, nao, nmo, &
385 : nspins
386 : LOGICAL :: do_periodic, do_raman, should_stop
387 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
388 112 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h1_psi0, psi0_order, psi1
389 112 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry
390 : TYPE(cp_fm_type), POINTER :: mo_coeff
391 : TYPE(cp_logger_type), POINTER :: logger
392 : TYPE(dft_control_type), POINTER :: dft_control
393 : TYPE(linres_control_type), POINTER :: linres_control
394 112 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
395 : TYPE(mp_para_env_type), POINTER :: para_env
396 : TYPE(polar_env_type), POINTER :: polar_env
397 : TYPE(qs_matrix_pools_type), POINTER :: mpools
398 : TYPE(section_vals_type), POINTER :: lr_section, polar_section
399 :
400 112 : CALL timeset(routineN, handle)
401 :
402 112 : NULLIFY (dft_control, linres_control, lr_section, polar_section)
403 112 : NULLIFY (logger, mpools, mo_coeff, para_env)
404 112 : NULLIFY (tmp_fm_struct, psi1_dBerry, dBerry_psi0)
405 :
406 112 : logger => cp_get_default_logger()
407 112 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
408 : polar_section => section_vals_get_subs_vals(qs_env%input, &
409 112 : "PROPERTIES%LINRES%POLAR")
410 :
411 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
412 112 : extension=".linresLog")
413 112 : IF (iounit > 0) THEN
414 : WRITE (UNIT=iounit, FMT="(T2,A,/)") &
415 56 : "POLAR| Self consistent optimization of the response wavefunctions"
416 : END IF
417 :
418 : CALL get_qs_env(qs_env=qs_env, &
419 : dft_control=dft_control, &
420 : mpools=mpools, &
421 : linres_control=linres_control, &
422 : mos=mos, &
423 : polar_env=polar_env, &
424 112 : para_env=para_env)
425 :
426 112 : nspins = dft_control%nspins
427 :
428 112 : CALL get_polar_env(polar_env=polar_env, do_raman=do_raman, do_periodic=do_periodic)
429 :
430 : ! Allocate the vectors
431 456 : ALLOCATE (psi0_order(nspins))
432 576 : ALLOCATE (psi1(nspins), h1_psi0(nspins))
433 232 : DO ispin = 1, nspins
434 120 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
435 120 : psi0_order(ispin) = mo_coeff
436 120 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao)
437 120 : NULLIFY (tmp_fm_struct)
438 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
439 : ncol_global=nmo, &
440 120 : context=mo_coeff%matrix_struct%context)
441 120 : CALL cp_fm_create(psi1(ispin), tmp_fm_struct)
442 120 : CALL cp_fm_create(h1_psi0(ispin), tmp_fm_struct)
443 352 : CALL cp_fm_struct_release(tmp_fm_struct)
444 : END DO
445 :
446 112 : IF (do_raman) THEN
447 : CALL get_polar_env(polar_env=polar_env, &
448 : psi1_dBerry=psi1_dBerry, &
449 112 : dBerry_psi0=dBerry_psi0)
450 448 : DO idir = 1, 3
451 808 : DO ispin = 1, nspins
452 696 : CALL cp_fm_set_all(psi1_dBerry(idir, ispin), 0.0_dp)
453 : END DO
454 : END DO
455 : ! Restart
456 112 : IF (linres_control%linres_restart) THEN
457 24 : DO idir = 1, 3
458 24 : CALL linres_read_restart(qs_env, lr_section, psi1_dBerry(idir, :), idir, "psi1_dBerry")
459 : END DO
460 : END IF
461 448 : loop_idir: DO idir = 1, 3
462 336 : IF (iounit > 0) THEN
463 168 : IF (do_periodic) THEN
464 : WRITE (iounit, "(/,T2,A)") &
465 27 : "POLAR| Response to the perturbation operator Berry phase_"//ACHAR(idir + 119)
466 : ELSE
467 : WRITE (iounit, "(/,T2,A)") &
468 141 : "POLAR| Response to the perturbation operator R_"//ACHAR(idir + 119)
469 : END IF
470 : END IF
471 : ! Do scf cycle to optimize psi1
472 696 : DO ispin = 1, nspins
473 360 : CALL cp_fm_to_fm(psi1_dBerry(idir, ispin), psi1(ispin))
474 696 : CALL cp_fm_to_fm(dBerry_psi0(idir, ispin), h1_psi0(ispin))
475 : END DO
476 : !
477 336 : linres_control%lr_triplet = .FALSE. ! we do singlet response
478 336 : linres_control%do_kernel = .TRUE. ! we do coupled response
479 336 : linres_control%converged = .FALSE.
480 336 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop)
481 :
482 : ! Copy the response
483 696 : DO ispin = 1, nspins
484 696 : CALL cp_fm_to_fm(psi1(ispin), psi1_dBerry(idir, ispin))
485 : END DO
486 : !
487 : ! Write the new result to the restart file
488 784 : IF (linres_control%linres_restart) THEN
489 18 : CALL linres_write_restart(qs_env, lr_section, psi1_dBerry(idir, :), idir, "psi1_dBerry")
490 : END IF
491 : END DO loop_idir
492 : END IF ! do_raman
493 :
494 112 : CALL set_polar_env(polar_env, run_stopped=should_stop)
495 :
496 : ! Clean up
497 112 : CALL cp_fm_release(psi1)
498 112 : CALL cp_fm_release(h1_psi0)
499 :
500 112 : DEALLOCATE (psi0_order)
501 :
502 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
503 112 : "PRINT%PROGRAM_RUN_INFO")
504 :
505 112 : CALL timestop(handle)
506 :
507 224 : END SUBROUTINE polar_response
508 :
509 : ! **************************************************************************************************
510 : !> \brief Prints the polarisability tensor to a file during MD runs
511 : !> \param force_env ...
512 : !> \param motion_section ...
513 : !> \param itimes ...
514 : !> \param time ...
515 : !> \param pos ...
516 : !> \param act ...
517 : !> \par History
518 : !> 06.2018 Creation (MK)
519 : !> \author Matthias Krack (MK)
520 : ! **************************************************************************************************
521 3768 : SUBROUTINE write_polarisability_tensor(force_env, motion_section, itimes, time, pos, act)
522 :
523 : TYPE(force_env_type), POINTER :: force_env
524 : TYPE(section_vals_type), POINTER :: motion_section
525 : INTEGER, INTENT(IN) :: itimes
526 : REAL(KIND=dp), INTENT(IN) :: time
527 : CHARACTER(LEN=default_string_length), INTENT(IN), &
528 : OPTIONAL :: pos, act
529 :
530 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
531 : INTEGER :: iounit
532 : LOGICAL :: do_raman, new_file, run_stopped
533 3768 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: polar
534 : TYPE(cp_logger_type), POINTER :: logger
535 : TYPE(polar_env_type), POINTER :: polar_env
536 : TYPE(qs_environment_type), POINTER :: qs_env
537 :
538 3768 : NULLIFY (qs_env)
539 :
540 3768 : CALL force_env_get(force_env, qs_env=qs_env)
541 3768 : IF (ASSOCIATED(qs_env)) THEN
542 3768 : NULLIFY (polar_env)
543 3768 : CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
544 3768 : IF (ASSOCIATED(polar_env)) THEN
545 : CALL get_polar_env(polar_env=polar_env, &
546 : polar=polar, &
547 : do_raman=do_raman, &
548 6 : run_stopped=run_stopped)
549 6 : IF (.NOT. run_stopped .AND. do_raman) THEN
550 6 : NULLIFY (logger)
551 6 : logger => cp_get_default_logger()
552 6 : my_pos = "APPEND"
553 6 : my_act = "WRITE"
554 6 : IF (PRESENT(pos)) my_pos = pos
555 6 : IF (PRESENT(act)) my_act = act
556 : iounit = cp_print_key_unit_nr(logger, motion_section, "PRINT%POLAR_MATRIX", &
557 : extension=".polar", file_position=my_pos, &
558 : file_action=my_act, file_form="FORMATTED", &
559 6 : is_new_file=new_file)
560 : ELSE
561 0 : iounit = 0
562 : END IF
563 6 : IF (iounit > 0) THEN
564 3 : IF (new_file) THEN
565 : WRITE (UNIT=iounit, FMT='(A,9(11X,A2," [a.u.]"),6X,A)') &
566 1 : "# Step Time [fs]", "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"
567 : END IF
568 3 : WRITE (UNIT=iounit, FMT='(I8,F12.3,9(1X,F19.8))') itimes, time, &
569 3 : polar(1, 1), polar(1, 2), polar(1, 3), &
570 3 : polar(2, 1), polar(2, 2), polar(2, 3), &
571 6 : polar(3, 1), polar(3, 2), polar(3, 3)
572 3 : CALL m_flush(iounit)
573 3 : CALL cp_print_key_finished_output(iounit, logger, motion_section, "PRINT%POLAR_MATRIX")
574 : END IF
575 : END IF ! polar_env
576 : END IF ! qs_env
577 :
578 3768 : END SUBROUTINE write_polarisability_tensor
579 :
580 18 : END MODULE qs_linres_polar_utils
|
