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