Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 154 : 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 154 : 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 154 : 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 154 : CALL timeset(routineN, handle)
105 :
106 154 : NULLIFY (dft_control)
107 154 : NULLIFY (linres_control)
108 154 : NULLIFY (logger)
109 154 : NULLIFY (matrix_s)
110 154 : NULLIFY (mos)
111 154 : NULLIFY (polar_env)
112 154 : NULLIFY (lr_section, polar_section)
113 :
114 154 : logger => cp_get_default_logger()
115 154 : 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 154 : extension=".linresLog")
119 :
120 154 : IF (iounit > 0) THEN
121 77 : WRITE (iounit, "(/,(T2,A))") "POLAR| Starting polarizability calculation", &
122 154 : "POLAR| Initialization of the polar environment"
123 : END IF
124 :
125 : polar_section => section_vals_get_subs_vals(qs_env%input, &
126 154 : "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 154 : mos=mos)
134 :
135 : ! Create polar environment if needed
136 154 : IF (.NOT. ASSOCIATED(polar_env)) THEN
137 126 : ALLOCATE (polar_env)
138 126 : CALL set_qs_env(qs_env=qs_env, polar_env=polar_env)
139 : END IF
140 :
141 154 : nspins = dft_control%nspins
142 :
143 154 : CALL section_vals_val_get(polar_section, "DO_RAMAN", l_val=polar_env%do_raman)
144 154 : 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 154 : IF (.NOT. ASSOCIATED(polar_env%polar)) THEN
148 126 : ALLOCATE (polar_env%polar(3, 3))
149 1638 : polar_env%polar(:, :) = 0.0_dp
150 : END IF
151 154 : IF (.NOT. ASSOCIATED(polar_env%dBerry_psi0)) THEN
152 938 : 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 154 : IF (.NOT. ASSOCIATED(polar_env%psi1_dBerry)) THEN
162 938 : 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 322 : DO ispin = 1, nspins
172 168 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
173 168 : CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao)
174 168 : NULLIFY (tmp_fm_struct)
175 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
176 : ncol_global=m, &
177 168 : context=mo_coeff%matrix_struct%context)
178 672 : DO idir = 1, 3
179 504 : CALL cp_fm_create(polar_env%dBerry_psi0(idir, ispin), tmp_fm_struct)
180 672 : CALL cp_fm_create(polar_env%psi1_dBerry(idir, ispin), tmp_fm_struct)
181 : END DO
182 490 : CALL cp_fm_struct_release(tmp_fm_struct)
183 : END DO
184 :
185 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
186 154 : "PRINT%PROGRAM_RUN_INFO")
187 :
188 154 : CALL timestop(handle)
189 :
190 154 : 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 150 : 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 150 : REAL(dp), DIMENSION(:, :), POINTER :: polar, polar_tmp
208 : TYPE(cell_type), POINTER :: cell
209 150 : 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 150 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
213 : TYPE(polar_env_type), POINTER :: polar_env
214 :
215 150 : CALL timeset(routineN, handle)
216 :
217 150 : NULLIFY (cell, dft_control, polar, psi1_dBerry, logger)
218 150 : NULLIFY (mos, dBerry_psi0)
219 150 : logger => cp_get_default_logger()
220 150 : 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 150 : polar_env=polar_env)
227 :
228 150 : nspins = dft_control%nspins
229 :
230 : CALL get_polar_env(polar_env=polar_env, &
231 : do_raman=do_raman, &
232 150 : run_stopped=run_stopped)
233 :
234 150 : IF (.NOT. run_stopped .AND. do_raman) THEN
235 :
236 150 : 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 150 : psi1_dBerry=psi1_dBerry)
243 :
244 : ! Initialize
245 150 : ALLOCATE (polar_tmp(3, 3))
246 1950 : polar_tmp(:, :) = 0.0_dp
247 :
248 600 : DO i = 1, 3 ! directions of electric field
249 1950 : DO z = 1, 3 !dipole directions
250 3258 : DO ispin = 1, dft_control%nspins
251 : !SL compute trace
252 : ptmp = 0.0_dp
253 1458 : CALL cp_fm_trace(psi1_dBerry(i, ispin), dBerry_psi0(z, ispin), ptmp)
254 2808 : polar_tmp(i, z) = polar_tmp(i, z) - 2.0_dp*ptmp
255 : END DO
256 : END DO
257 : END DO !spin
258 :
259 150 : IF (do_periodic) THEN
260 1820 : polar(:, :) = MATMUL(MATMUL(cell%hmat, polar_tmp), TRANSPOSE(cell%hmat))/(twopi*twopi)
261 : ELSE
262 3380 : polar(:, :) = polar_tmp(:, :)
263 : END IF
264 : !SL evtl maxocc instead?
265 150 : IF (dft_control%nspins == 1) THEN
266 1794 : polar(:, :) = 2.0_dp*polar(:, :)
267 : END IF
268 :
269 150 : IF (ASSOCIATED(polar_tmp)) THEN
270 150 : DEALLOCATE (polar_tmp)
271 : END IF
272 :
273 : END IF ! do_raman
274 :
275 150 : CALL timestop(handle)
276 :
277 150 : 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 150 : 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 150 : 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 150 : 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 150 : para_env=para_env)
307 :
308 150 : logger => cp_get_default_logger()
309 150 : iounit = cp_logger_get_default_io_unit(logger)
310 :
311 150 : 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 150 : run_stopped=run_stopped)
317 :
318 150 : IF (.NOT. run_stopped .AND. do_raman) THEN
319 :
320 150 : description = "[POLAR]"
321 150 : CALL cp_results_erase(results, description=description)
322 150 : CALL put_results(results, description=description, values=polar(:, :))
323 :
324 150 : 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 144 : extension=".data", middle_name="raman", log_filename=.FALSE.)
329 144 : IF (unit_p > 0) THEN
330 72 : IF (unit_p /= iounit) THEN
331 72 : WRITE (unit_p, *)
332 72 : WRITE (unit_p, '(T10,A)') 'POLARIZABILITY TENSOR (atomic units):'
333 72 : WRITE (unit_p, '(T10,A,3F15.5)') "xx,yy,zz", polar(1, 1), polar(2, 2), polar(3, 3)
334 72 : WRITE (unit_p, '(T10,A,3F15.5)') "xy,xz,yz", polar(1, 2), polar(1, 3), polar(2, 3)
335 72 : WRITE (unit_p, '(T10,A,3F15.5)') "yx,zx,zy", polar(2, 1), polar(3, 1), polar(3, 2)
336 72 : WRITE (unit_p, '(T10,A)') 'POLARIZABILITY TENSOR (Angstrom^3):'
337 72 : WRITE (unit_p, '(T10,A,3F15.5)') "xx,yy,zz", polar(1, 1)*angstrom**3, &
338 144 : polar(2, 2)*angstrom**3, polar(3, 3)*angstrom**3
339 72 : WRITE (unit_p, '(T10,A,3F15.5)') "xy,xz,yz", polar(1, 2)*angstrom**3, &
340 144 : polar(1, 3)*angstrom**3, polar(2, 3)*angstrom**3
341 72 : WRITE (unit_p, '(T10,A,3F15.5)') "yx,zx,zy", polar(2, 1)*angstrom**3, &
342 144 : polar(3, 1)*angstrom**3, polar(3, 2)*angstrom**3
343 : CALL cp_print_key_finished_output(unit_p, logger, polar_section, &
344 72 : "PRINT%POLAR_MATRIX")
345 : END IF
346 : END IF
347 : END IF
348 150 : IF (iounit > 0) THEN
349 : WRITE (iounit, '(/,T2,A)') &
350 75 : 'POLAR| Polarizability tensor [a.u.]'
351 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
352 75 : 'POLAR| xx,yy,zz', polar(1, 1), polar(2, 2), polar(3, 3)
353 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
354 75 : 'POLAR| xy,xz,yz', polar(1, 2), polar(1, 3), polar(2, 3)
355 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
356 75 : 'POLAR| yx,zx,zy', polar(2, 1), polar(3, 1), polar(3, 2)
357 : WRITE (iounit, '(/,T2,A)') &
358 75 : 'POLAR| Polarizability tensor [ang^3]'
359 : WRITE (iounit, '(T2,A,T24,3(1X,F18.12))') &
360 75 : '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 75 : '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 75 : 'POLAR| yx,zx,zy', polar(2, 1)*angstrom**3, polar(3, 1)*angstrom**3, polar(3, 2)*angstrom**3
365 : END IF
366 : IF (iounit > 0) THEN
367 975 : WRITE (UNIT=iounit, FMT="(T2,A,T61,E20.12)") ' POLAR : CheckSum =', SUM(polar)
368 : END IF
369 : END IF
370 :
371 150 : END SUBROUTINE polar_print
372 :
373 : ! **************************************************************************************************
374 : !> \brief Calculate the polarisability tensor using response theory
375 : !> \param p_env ...
376 : !> \param qs_env ...
377 : !> \par History
378 : !> 06.2018 polar_env integrated into qs_env (MK)
379 : ! **************************************************************************************************
380 154 : SUBROUTINE polar_response(p_env, qs_env)
381 :
382 : TYPE(qs_p_env_type) :: p_env
383 : TYPE(qs_environment_type), POINTER :: qs_env
384 :
385 : CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_response'
386 :
387 : INTEGER :: handle, idir, iounit, ispin, nao, nmo, &
388 : nspins
389 : LOGICAL :: do_periodic, do_raman, should_stop
390 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
391 154 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h1_psi0, psi0_order, psi1
392 154 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry
393 : TYPE(cp_fm_type), POINTER :: mo_coeff
394 : TYPE(cp_logger_type), POINTER :: logger
395 : TYPE(dft_control_type), POINTER :: dft_control
396 : TYPE(linres_control_type), POINTER :: linres_control
397 154 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
398 : TYPE(mp_para_env_type), POINTER :: para_env
399 : TYPE(polar_env_type), POINTER :: polar_env
400 : TYPE(qs_matrix_pools_type), POINTER :: mpools
401 : TYPE(section_vals_type), POINTER :: lr_section, polar_section
402 :
403 154 : CALL timeset(routineN, handle)
404 :
405 154 : NULLIFY (dft_control, linres_control, lr_section, polar_section)
406 154 : NULLIFY (logger, mpools, mo_coeff, para_env)
407 154 : NULLIFY (tmp_fm_struct, psi1_dBerry, dBerry_psi0)
408 :
409 154 : logger => cp_get_default_logger()
410 154 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
411 : polar_section => section_vals_get_subs_vals(qs_env%input, &
412 154 : "PROPERTIES%LINRES%POLAR")
413 :
414 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
415 154 : extension=".linresLog")
416 154 : IF (iounit > 0) THEN
417 : WRITE (UNIT=iounit, FMT="(T2,A,/)") &
418 77 : "POLAR| Self consistent optimization of the response wavefunctions"
419 : END IF
420 :
421 : CALL get_qs_env(qs_env=qs_env, &
422 : dft_control=dft_control, &
423 : mpools=mpools, &
424 : linres_control=linres_control, &
425 : mos=mos, &
426 : polar_env=polar_env, &
427 154 : para_env=para_env)
428 :
429 154 : nspins = dft_control%nspins
430 :
431 154 : CALL get_polar_env(polar_env=polar_env, do_raman=do_raman, do_periodic=do_periodic)
432 :
433 : ! Allocate the vectors
434 630 : ALLOCATE (psi0_order(nspins))
435 798 : ALLOCATE (psi1(nspins), h1_psi0(nspins))
436 322 : DO ispin = 1, nspins
437 168 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
438 168 : psi0_order(ispin) = mo_coeff
439 168 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao)
440 168 : NULLIFY (tmp_fm_struct)
441 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
442 : ncol_global=nmo, &
443 168 : context=mo_coeff%matrix_struct%context)
444 168 : CALL cp_fm_create(psi1(ispin), tmp_fm_struct)
445 168 : CALL cp_fm_create(h1_psi0(ispin), tmp_fm_struct)
446 490 : CALL cp_fm_struct_release(tmp_fm_struct)
447 : END DO
448 :
449 154 : IF (do_raman) THEN
450 : CALL get_polar_env(polar_env=polar_env, &
451 : psi1_dBerry=psi1_dBerry, &
452 154 : dBerry_psi0=dBerry_psi0)
453 616 : DO idir = 1, 3
454 1120 : DO ispin = 1, nspins
455 966 : CALL cp_fm_set_all(psi1_dBerry(idir, ispin), 0.0_dp)
456 : END DO
457 : END DO
458 : ! Restart
459 154 : IF (linres_control%linres_restart) THEN
460 24 : DO idir = 1, 3
461 24 : CALL linres_read_restart(qs_env, lr_section, psi1_dBerry(idir, :), idir, "psi1_dBerry")
462 : END DO
463 : END IF
464 616 : loop_idir: DO idir = 1, 3
465 462 : IF (iounit > 0) THEN
466 231 : IF (do_periodic) THEN
467 : WRITE (iounit, "(/,T2,A)") &
468 30 : "POLAR| Response to the perturbation operator Berry phase_"//ACHAR(idir + 119)
469 : ELSE
470 : WRITE (iounit, "(/,T2,A)") &
471 201 : "POLAR| Response to the perturbation operator R_"//ACHAR(idir + 119)
472 : END IF
473 : END IF
474 : ! Do scf cycle to optimize psi1
475 966 : DO ispin = 1, nspins
476 504 : CALL cp_fm_to_fm(psi1_dBerry(idir, ispin), psi1(ispin))
477 966 : CALL cp_fm_to_fm(dBerry_psi0(idir, ispin), h1_psi0(ispin))
478 : END DO
479 : !
480 462 : linres_control%lr_triplet = .FALSE. ! we do singlet response
481 462 : linres_control%do_kernel = .TRUE. ! we do coupled response
482 462 : linres_control%converged = .FALSE.
483 462 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop)
484 :
485 : ! Copy the response
486 966 : DO ispin = 1, nspins
487 966 : CALL cp_fm_to_fm(psi1(ispin), psi1_dBerry(idir, ispin))
488 : END DO
489 : !
490 : ! Write the new result to the restart file
491 1078 : IF (linres_control%linres_restart) THEN
492 18 : CALL linres_write_restart(qs_env, lr_section, psi1_dBerry(idir, :), idir, "psi1_dBerry")
493 : END IF
494 : END DO loop_idir
495 : END IF ! do_raman
496 :
497 154 : CALL set_polar_env(polar_env, run_stopped=should_stop)
498 :
499 : ! Clean up
500 154 : CALL cp_fm_release(psi1)
501 154 : CALL cp_fm_release(h1_psi0)
502 :
503 154 : DEALLOCATE (psi0_order)
504 :
505 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
506 154 : "PRINT%PROGRAM_RUN_INFO")
507 :
508 154 : CALL timestop(handle)
509 :
510 308 : END SUBROUTINE polar_response
511 :
512 : ! **************************************************************************************************
513 : !> \brief Prints the polarisability tensor to a file during MD runs
514 : !> \param force_env ...
515 : !> \param motion_section ...
516 : !> \param itimes ...
517 : !> \param time ...
518 : !> \param pos ...
519 : !> \param act ...
520 : !> \par History
521 : !> 06.2018 Creation (MK)
522 : !> \author Matthias Krack (MK)
523 : ! **************************************************************************************************
524 3834 : SUBROUTINE write_polarisability_tensor(force_env, motion_section, itimes, time, pos, act)
525 :
526 : TYPE(force_env_type), POINTER :: force_env
527 : TYPE(section_vals_type), POINTER :: motion_section
528 : INTEGER, INTENT(IN) :: itimes
529 : REAL(KIND=dp), INTENT(IN) :: time
530 : CHARACTER(LEN=default_string_length), INTENT(IN), &
531 : OPTIONAL :: pos, act
532 :
533 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
534 : INTEGER :: iounit
535 : LOGICAL :: do_raman, new_file, run_stopped
536 3834 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: polar
537 : TYPE(cp_logger_type), POINTER :: logger
538 : TYPE(polar_env_type), POINTER :: polar_env
539 : TYPE(qs_environment_type), POINTER :: qs_env
540 :
541 3834 : NULLIFY (qs_env)
542 :
543 3834 : CALL force_env_get(force_env, qs_env=qs_env)
544 3834 : IF (ASSOCIATED(qs_env)) THEN
545 3834 : NULLIFY (polar_env)
546 3834 : CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
547 3834 : IF (ASSOCIATED(polar_env)) THEN
548 : CALL get_polar_env(polar_env=polar_env, &
549 : polar=polar, &
550 : do_raman=do_raman, &
551 6 : run_stopped=run_stopped)
552 6 : IF (.NOT. run_stopped .AND. do_raman) THEN
553 6 : NULLIFY (logger)
554 6 : logger => cp_get_default_logger()
555 6 : my_pos = "APPEND"
556 6 : my_act = "WRITE"
557 6 : IF (PRESENT(pos)) my_pos = pos
558 6 : IF (PRESENT(act)) my_act = act
559 : iounit = cp_print_key_unit_nr(logger, motion_section, "PRINT%POLAR_MATRIX", &
560 : extension=".polar", file_position=my_pos, &
561 : file_action=my_act, file_form="FORMATTED", &
562 6 : is_new_file=new_file)
563 : ELSE
564 0 : iounit = 0
565 : END IF
566 6 : IF (iounit > 0) THEN
567 3 : IF (new_file) THEN
568 : WRITE (UNIT=iounit, FMT='(A,9(11X,A2," [a.u.]"),6X,A)') &
569 1 : "# Step Time [fs]", "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"
570 : END IF
571 3 : WRITE (UNIT=iounit, FMT='(I8,F12.3,9(1X,F19.8))') itimes, time, &
572 3 : polar(1, 1), polar(1, 2), polar(1, 3), &
573 3 : polar(2, 1), polar(2, 2), polar(2, 3), &
574 6 : polar(3, 1), polar(3, 2), polar(3, 3)
575 3 : CALL m_flush(iounit)
576 3 : CALL cp_print_key_finished_output(iounit, logger, motion_section, "PRINT%POLAR_MATRIX")
577 : END IF
578 : END IF ! polar_env
579 : END IF ! qs_env
580 :
581 3834 : END SUBROUTINE write_polarisability_tensor
582 :
583 20 : END MODULE qs_linres_polar_utils
|
