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 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_fm_basic_linalg, ONLY: cp_fm_trace
23 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
24 : cp_fm_struct_release,&
25 : cp_fm_struct_type
26 : USE cp_fm_types, ONLY: cp_fm_create,&
27 : cp_fm_get_info,&
28 : cp_fm_release,&
29 : cp_fm_set_all,&
30 : cp_fm_to_fm,&
31 : cp_fm_type
32 : USE cp_log_handling, ONLY: cp_get_default_logger,&
33 : cp_logger_get_default_io_unit,&
34 : cp_logger_type
35 : USE cp_output_handling, ONLY: cp_p_file,&
36 : cp_print_key_finished_output,&
37 : cp_print_key_should_output,&
38 : cp_print_key_unit_nr
39 : USE cp_result_methods, ONLY: cp_results_erase,&
40 : put_results
41 : USE cp_result_types, ONLY: cp_result_type
42 : USE dbcsr_api, ONLY: dbcsr_p_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 108 : 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 108 : 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 108 : 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 108 : CALL timeset(routineN, handle)
105 :
106 108 : NULLIFY (dft_control)
107 108 : NULLIFY (linres_control)
108 108 : NULLIFY (logger)
109 108 : NULLIFY (matrix_s)
110 108 : NULLIFY (mos)
111 108 : NULLIFY (polar_env)
112 108 : NULLIFY (lr_section, polar_section)
113 :
114 108 : logger => cp_get_default_logger()
115 108 : 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 108 : extension=".linresLog")
119 :
120 108 : IF (iounit > 0) THEN
121 54 : WRITE (iounit, "(/,(T2,A))") "POLAR| Starting polarizability calculation", &
122 108 : "POLAR| Initialization of the polar environment"
123 : END IF
124 :
125 : polar_section => section_vals_get_subs_vals(qs_env%input, &
126 108 : "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 108 : mos=mos)
134 :
135 : ! Create polar environment if needed
136 108 : IF (.NOT. ASSOCIATED(polar_env)) THEN
137 80 : ALLOCATE (polar_env)
138 80 : CALL set_qs_env(qs_env=qs_env, polar_env=polar_env)
139 : END IF
140 :
141 108 : nspins = dft_control%nspins
142 :
143 108 : CALL section_vals_val_get(polar_section, "DO_RAMAN", l_val=polar_env%do_raman)
144 108 : 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 108 : IF (.NOT. ASSOCIATED(polar_env%polar)) THEN
148 80 : ALLOCATE (polar_env%polar(3, 3))
149 1040 : polar_env%polar(:, :) = 0.0_dp
150 : END IF
151 108 : IF (.NOT. ASSOCIATED(polar_env%dBerry_psi0)) THEN
152 584 : 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 108 : IF (.NOT. ASSOCIATED(polar_env%psi1_dBerry)) THEN
162 584 : 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 222 : DO ispin = 1, nspins
172 114 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
173 114 : CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao)
174 114 : NULLIFY (tmp_fm_struct)
175 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
176 : ncol_global=m, &
177 114 : context=mo_coeff%matrix_struct%context)
178 456 : DO idir = 1, 3
179 342 : CALL cp_fm_create(polar_env%dBerry_psi0(idir, ispin), tmp_fm_struct)
180 456 : CALL cp_fm_create(polar_env%psi1_dBerry(idir, ispin), tmp_fm_struct)
181 : END DO
182 336 : CALL cp_fm_struct_release(tmp_fm_struct)
183 : END DO
184 :
185 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
186 108 : "PRINT%PROGRAM_RUN_INFO")
187 :
188 108 : CALL timestop(handle)
189 :
190 108 : 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 108 : 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 108 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h1_psi0, psi0_order, psi1
389 108 : 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 108 : 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 108 : CALL timeset(routineN, handle)
401 :
402 108 : NULLIFY (dft_control, linres_control, lr_section, polar_section)
403 108 : NULLIFY (logger, mpools, mo_coeff, para_env)
404 108 : NULLIFY (tmp_fm_struct, psi1_dBerry, dBerry_psi0)
405 :
406 108 : logger => cp_get_default_logger()
407 108 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
408 : polar_section => section_vals_get_subs_vals(qs_env%input, &
409 108 : "PROPERTIES%LINRES%POLAR")
410 :
411 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
412 108 : extension=".linresLog")
413 108 : IF (iounit > 0) THEN
414 : WRITE (UNIT=iounit, FMT="(T2,A,/)") &
415 54 : "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 108 : para_env=para_env)
425 :
426 108 : nspins = dft_control%nspins
427 :
428 108 : CALL get_polar_env(polar_env=polar_env, do_raman=do_raman, do_periodic=do_periodic)
429 :
430 : ! Allocate the vectors
431 438 : ALLOCATE (psi0_order(nspins))
432 768 : ALLOCATE (psi1(nspins), h1_psi0(nspins))
433 222 : DO ispin = 1, nspins
434 114 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
435 114 : psi0_order(ispin) = mo_coeff
436 114 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao)
437 114 : NULLIFY (tmp_fm_struct)
438 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
439 : ncol_global=nmo, &
440 114 : context=mo_coeff%matrix_struct%context)
441 114 : CALL cp_fm_create(psi1(ispin), tmp_fm_struct)
442 114 : CALL cp_fm_create(h1_psi0(ispin), tmp_fm_struct)
443 336 : CALL cp_fm_struct_release(tmp_fm_struct)
444 : END DO
445 :
446 108 : IF (do_raman) THEN
447 : CALL get_polar_env(polar_env=polar_env, &
448 : psi1_dBerry=psi1_dBerry, &
449 108 : dBerry_psi0=dBerry_psi0)
450 : ! Restart
451 108 : IF (linres_control%linres_restart) THEN
452 24 : DO idir = 1, 3
453 24 : CALL linres_read_restart(qs_env, lr_section, psi1_dBerry(idir, :), idir, "psi1_dBerry")
454 : END DO
455 : ELSE
456 408 : DO idir = 1, 3
457 732 : DO ispin = 1, nspins
458 630 : CALL cp_fm_set_all(psi1_dBerry(idir, ispin), 0.0_dp)
459 : END DO
460 : END DO
461 : END IF
462 432 : loop_idir: DO idir = 1, 3
463 324 : IF (iounit > 0) THEN
464 162 : IF (do_periodic) THEN
465 : WRITE (iounit, "(/,T2,A)") &
466 27 : "POLAR| Response to the perturbation operator Berry phase_"//ACHAR(idir + 119)
467 : ELSE
468 : WRITE (iounit, "(/,T2,A)") &
469 135 : "POLAR| Response to the perturbation operator R_"//ACHAR(idir + 119)
470 : END IF
471 : END IF
472 : ! Do scf cycle to optimize psi1
473 666 : DO ispin = 1, nspins
474 342 : CALL cp_fm_to_fm(psi1_dBerry(idir, ispin), psi1(ispin))
475 666 : CALL cp_fm_to_fm(dBerry_psi0(idir, ispin), h1_psi0(ispin))
476 : END DO
477 : !
478 324 : linres_control%lr_triplet = .FALSE. ! we do singlet response
479 324 : linres_control%do_kernel = .TRUE. ! we do coupled response
480 324 : linres_control%converged = .FALSE.
481 324 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop)
482 :
483 : ! Copy the response
484 666 : DO ispin = 1, nspins
485 666 : CALL cp_fm_to_fm(psi1(ispin), psi1_dBerry(idir, ispin))
486 : END DO
487 : !
488 : ! Write the new result to the restart file
489 756 : IF (linres_control%linres_restart) THEN
490 18 : CALL linres_write_restart(qs_env, lr_section, psi1_dBerry(idir, :), idir, "psi1_dBerry")
491 : END IF
492 : END DO loop_idir
493 : END IF ! do_raman
494 :
495 108 : CALL set_polar_env(polar_env, run_stopped=should_stop)
496 :
497 : ! Clean up
498 108 : CALL cp_fm_release(psi1)
499 108 : CALL cp_fm_release(h1_psi0)
500 :
501 108 : DEALLOCATE (psi0_order)
502 :
503 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
504 108 : "PRINT%PROGRAM_RUN_INFO")
505 :
506 108 : CALL timestop(handle)
507 :
508 216 : END SUBROUTINE polar_response
509 :
510 : ! **************************************************************************************************
511 : !> \brief Prints the polarisability tensor to a file during MD runs
512 : !> \param force_env ...
513 : !> \param motion_section ...
514 : !> \param itimes ...
515 : !> \param time ...
516 : !> \param pos ...
517 : !> \param act ...
518 : !> \par History
519 : !> 06.2018 Creation (MK)
520 : !> \author Matthias Krack (MK)
521 : ! **************************************************************************************************
522 3632 : SUBROUTINE write_polarisability_tensor(force_env, motion_section, itimes, time, pos, act)
523 :
524 : TYPE(force_env_type), POINTER :: force_env
525 : TYPE(section_vals_type), POINTER :: motion_section
526 : INTEGER, INTENT(IN) :: itimes
527 : REAL(KIND=dp), INTENT(IN) :: time
528 : CHARACTER(LEN=default_string_length), INTENT(IN), &
529 : OPTIONAL :: pos, act
530 :
531 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
532 : INTEGER :: iounit
533 : LOGICAL :: do_raman, new_file, run_stopped
534 3632 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: polar
535 : TYPE(cp_logger_type), POINTER :: logger
536 : TYPE(polar_env_type), POINTER :: polar_env
537 : TYPE(qs_environment_type), POINTER :: qs_env
538 :
539 3632 : NULLIFY (qs_env)
540 :
541 3632 : CALL force_env_get(force_env, qs_env=qs_env)
542 3632 : IF (ASSOCIATED(qs_env)) THEN
543 3632 : NULLIFY (polar_env)
544 3632 : CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
545 3632 : IF (ASSOCIATED(polar_env)) THEN
546 : CALL get_polar_env(polar_env=polar_env, &
547 : polar=polar, &
548 : do_raman=do_raman, &
549 6 : run_stopped=run_stopped)
550 6 : IF (.NOT. run_stopped .AND. do_raman) THEN
551 6 : NULLIFY (logger)
552 6 : logger => cp_get_default_logger()
553 6 : my_pos = "APPEND"
554 6 : my_act = "WRITE"
555 6 : IF (PRESENT(pos)) my_pos = pos
556 6 : IF (PRESENT(act)) my_act = act
557 : iounit = cp_print_key_unit_nr(logger, motion_section, "PRINT%POLAR_MATRIX", &
558 : extension=".polar", file_position=my_pos, &
559 : file_action=my_act, file_form="FORMATTED", &
560 6 : is_new_file=new_file)
561 : ELSE
562 0 : iounit = 0
563 : END IF
564 6 : IF (iounit > 0) THEN
565 3 : IF (new_file) THEN
566 : WRITE (UNIT=iounit, FMT='(A,9(11X,A2," [a.u.]"),6X,A)') &
567 1 : "# Step Time [fs]", "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"
568 : END IF
569 3 : WRITE (UNIT=iounit, FMT='(I8,F12.3,9(1X,F19.8))') itimes, time, &
570 3 : polar(1, 1), polar(1, 2), polar(1, 3), &
571 3 : polar(2, 1), polar(2, 2), polar(2, 3), &
572 6 : polar(3, 1), polar(3, 2), polar(3, 3)
573 3 : CALL m_flush(iounit)
574 3 : CALL cp_print_key_finished_output(iounit, logger, motion_section, "PRINT%POLAR_MATRIX")
575 : END IF
576 : END IF ! polar_env
577 : END IF ! qs_env
578 :
579 3632 : END SUBROUTINE write_polarisability_tensor
580 :
581 18 : END MODULE qs_linres_polar_utils
|
