1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Contains the setup for the calculation of properties by linear response 8!> by the application of second order density functional perturbation theory. 9!> The knowledge of the ground state energy, density and wavefunctions is assumed. 10!> Uses the self consistent approach. 11!> Properties that can be calculated : none 12!> \par History 13!> created 06-2005 [MI] 14!> \author MI 15! ************************************************************************************************** 16MODULE qs_linres_module 17 USE bibliography, ONLY: Weber2009,& 18 cite_reference 19 USE cp_control_types, ONLY: dft_control_type 20 USE cp_log_handling, ONLY: cp_get_default_logger,& 21 cp_logger_type 22 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 23 cp_print_key_unit_nr 24 USE dbcsr_api, ONLY: dbcsr_p_type 25 USE force_env_types, ONLY: force_env_get,& 26 force_env_type,& 27 use_qmmm,& 28 use_qs_force 29 USE input_constants, ONLY: lr_current,& 30 lr_none,& 31 ot_precond_full_all,& 32 ot_precond_full_kinetic,& 33 ot_precond_full_single,& 34 ot_precond_full_single_inverse,& 35 ot_precond_none,& 36 ot_precond_s_inverse 37 USE input_section_types, ONLY: section_vals_get,& 38 section_vals_get_subs_vals,& 39 section_vals_type,& 40 section_vals_val_get 41 USE qs_environment_types, ONLY: get_qs_env,& 42 qs_environment_type,& 43 set_qs_env 44 USE qs_linres_current, ONLY: current_build_chi,& 45 current_build_current 46 USE qs_linres_current_utils, ONLY: current_env_cleanup,& 47 current_env_init,& 48 current_response 49 USE qs_linres_epr_nablavks, ONLY: epr_nablavks 50 USE qs_linres_epr_ownutils, ONLY: epr_g_print,& 51 epr_g_so,& 52 epr_g_soo,& 53 epr_g_zke,& 54 epr_ind_magnetic_field 55 USE qs_linres_epr_utils, ONLY: epr_env_cleanup,& 56 epr_env_init 57 USE qs_linres_issc_utils, ONLY: issc_env_cleanup,& 58 issc_env_init,& 59 issc_issc,& 60 issc_print,& 61 issc_response 62 USE qs_linres_methods, ONLY: linres_localize 63 USE qs_linres_nmr_shift, ONLY: nmr_shift,& 64 nmr_shift_print 65 USE qs_linres_nmr_utils, ONLY: nmr_env_cleanup,& 66 nmr_env_init 67 USE qs_linres_op, ONLY: current_operators,& 68 issc_operators,& 69 polar_operators 70 USE qs_linres_polar_utils, ONLY: polar_env_init,& 71 polar_polar,& 72 polar_print,& 73 polar_response 74 USE qs_linres_types, ONLY: current_env_type,& 75 epr_env_type,& 76 issc_env_type,& 77 linres_control_create,& 78 linres_control_release,& 79 linres_control_type,& 80 nmr_env_type 81 USE qs_mo_methods, ONLY: calculate_density_matrix 82 USE qs_mo_types, ONLY: mo_set_p_type 83 USE qs_p_env_methods, ONLY: p_env_create,& 84 p_env_psi0_changed 85 USE qs_p_env_types, ONLY: p_env_release,& 86 qs_p_env_type 87 USE qs_rho_methods, ONLY: qs_rho_update_rho 88 USE qs_rho_types, ONLY: qs_rho_get,& 89 qs_rho_type 90#include "./base/base_uses.f90" 91 92 IMPLICIT NONE 93 94 PRIVATE 95 PUBLIC :: linres_calculation, linres_calculation_low 96 97 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_module' 98 99! ************************************************************************************************** 100 101CONTAINS 102 103! ************************************************************************************************** 104!> \brief Driver for the linear response calculatios 105!> \param force_env ... 106!> \par History 107!> 06.2005 created [MI] 108!> \author MI 109! ************************************************************************************************** 110 SUBROUTINE linres_calculation(force_env) 111 112 TYPE(force_env_type), POINTER :: force_env 113 114 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation', & 115 routineP = moduleN//':'//routineN 116 117 INTEGER :: handle 118 TYPE(qs_environment_type), POINTER :: qs_env 119 120 CALL timeset(routineN, handle) 121 122 NULLIFY (qs_env) 123 124 CPASSERT(ASSOCIATED(force_env)) 125 CPASSERT(force_env%ref_count > 0) 126 127 SELECT CASE (force_env%in_use) 128 CASE (use_qs_force) 129 CALL force_env_get(force_env, qs_env=qs_env) 130 CASE (use_qmmm) 131 qs_env => force_env%qmmm_env%qs_env 132 CASE DEFAULT 133 CPABORT("Does not recognize this force_env") 134 END SELECT 135 136 qs_env%linres_run = .TRUE. 137 138 CALL linres_calculation_low(qs_env) 139 140 CALL timestop(handle) 141 142 END SUBROUTINE linres_calculation 143 144! ************************************************************************************************** 145!> \brief Linear response can be called as run type or as post scf calculation 146!> Initialize the perturbation environment 147!> Define which properties is to be calculated 148!> Start up the optimization of the response density and wfn 149!> \param qs_env ... 150!> \par History 151!> 06.2005 created [MI] 152!> 02.2013 added polarizability section [SL] 153!> \author MI 154! ************************************************************************************************** 155 SUBROUTINE linres_calculation_low(qs_env) 156 157 TYPE(qs_environment_type), POINTER :: qs_env 158 159 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation_low', & 160 routineP = moduleN//':'//routineN 161 162 INTEGER :: handle, iounit 163 LOGICAL :: epr_present, issc_present, & 164 lr_calculation, nmr_present, & 165 polar_present 166 TYPE(cp_logger_type), POINTER :: logger 167 TYPE(dft_control_type), POINTER :: dft_control 168 TYPE(linres_control_type), POINTER :: linres_control 169 TYPE(qs_p_env_type), POINTER :: p_env 170 TYPE(section_vals_type), POINTER :: lr_section, prop_section 171 172 CALL timeset(routineN, handle) 173 174 lr_calculation = .FALSE. 175 nmr_present = .FALSE. 176 epr_present = .FALSE. 177 issc_present = .FALSE. 178 polar_present = .FALSE. 179 180 NULLIFY (dft_control, p_env, linres_control, logger, prop_section, lr_section) 181 logger => cp_get_default_logger() 182 183 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") 184 CALL section_vals_get(lr_section, explicit=lr_calculation) 185 186 IF (lr_calculation) THEN 187 CALL linres_init(lr_section, p_env, qs_env) 188 iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 189 extension=".linresLog") 190 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 191 linres_control=linres_control) 192 193 ! The type of perturbation has not been defined yet 194 linres_control%property = lr_none 195 196 ! We do NMR or EPR, then compute the current response 197 prop_section => section_vals_get_subs_vals(lr_section, "NMR") 198 CALL section_vals_get(prop_section, explicit=nmr_present) 199 prop_section => section_vals_get_subs_vals(lr_section, "EPR") 200 CALL section_vals_get(prop_section, explicit=epr_present) 201 202 IF (nmr_present .OR. epr_present) THEN 203 CALL nmr_epr_linres(linres_control, qs_env, p_env, dft_control, & 204 nmr_present, epr_present, iounit) 205 ENDIF 206 207 ! We do the indirect spin-spin coupling calculation 208 prop_section => section_vals_get_subs_vals(lr_section, "SPINSPIN") 209 CALL section_vals_get(prop_section, explicit=issc_present) 210 211 IF (issc_present) THEN 212 CALL issc_linres(linres_control, qs_env, p_env, dft_control) 213 ENDIF 214 215 ! We do the polarizability calculation 216 prop_section => section_vals_get_subs_vals(lr_section, "POLAR") 217 CALL section_vals_get(prop_section, explicit=polar_present) 218 IF (polar_present) THEN 219 CALL polar_linres(qs_env, p_env) 220 END IF 221 222 ! Other possible LR calculations can be introduced here 223 224 CALL p_env_release(p_env) 225 226 IF (iounit > 0) THEN 227 WRITE (UNIT=iounit, FMT="(/,T2,A,/,T25,A,/,T2,A,/)") & 228 REPEAT("=", 79), & 229 "ENDED LINRES CALCULATION", & 230 REPEAT("=", 79) 231 END IF 232 CALL cp_print_key_finished_output(iounit, logger, lr_section, & 233 "PRINT%PROGRAM_RUN_INFO") 234 END IF 235 236 CALL timestop(handle) 237 238 END SUBROUTINE linres_calculation_low 239 240! ************************************************************************************************** 241!> \brief Initialize some general settings like the p_env 242!> Localize the psi0 if required 243!> \param lr_section ... 244!> \param p_env ... 245!> \param qs_env ... 246!> \par History 247!> 06.2005 created [MI] 248!> \author MI 249!> \note 250!> - The localization should probably be always for all the occupied states 251! ************************************************************************************************** 252 SUBROUTINE linres_init(lr_section, p_env, qs_env) 253 254 TYPE(section_vals_type), POINTER :: lr_section 255 TYPE(qs_p_env_type), POINTER :: p_env 256 TYPE(qs_environment_type), POINTER :: qs_env 257 258 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_init', routineP = moduleN//':'//routineN 259 260 INTEGER :: iounit, ispin 261 LOGICAL :: do_it 262 TYPE(cp_logger_type), POINTER :: logger 263 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, rho_ao 264 TYPE(dft_control_type), POINTER :: dft_control 265 TYPE(linres_control_type), POINTER :: linres_control 266 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 267 TYPE(qs_rho_type), POINTER :: rho 268 TYPE(section_vals_type), POINTER :: loc_section 269 270 NULLIFY (logger) 271 logger => cp_get_default_logger() 272 iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 273 extension=".linresLog") 274 NULLIFY (dft_control, linres_control, loc_section, rho, mos, matrix_ks, rho_ao) 275 276 CPASSERT(.NOT. ASSOCIATED(p_env)) 277 278 CALL linres_control_create(linres_control) 279 CALL set_qs_env(qs_env=qs_env, linres_control=linres_control) 280 CALL linres_control_release(linres_control) 281 CALL get_qs_env(qs_env=qs_env, linres_control=linres_control, & 282 dft_control=dft_control, matrix_ks=matrix_ks, mos=mos, rho=rho) 283 CALL qs_rho_get(rho, rho_ao=rho_ao) 284 285 ! Localized Psi0 are required when the position operator has to be defined (nmr) 286 loc_section => section_vals_get_subs_vals(lr_section, "LOCALIZE") 287 CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", & 288 l_val=linres_control%localized_psi0) 289 IF (linres_control%localized_psi0) THEN 290 IF (iounit > 0) THEN 291 WRITE (UNIT=iounit, FMT="(/,T3,A,A)") & 292 "Localization of ground state orbitals", & 293 " before starting linear response calculation" 294 END IF 295 296 CALL linres_localize(qs_env, linres_control, dft_control%nspins) 297 298 DO ispin = 1, dft_control%nspins 299 CALL calculate_density_matrix(mos(ispin)%mo_set, rho_ao(ispin)%matrix) 300 ENDDO 301 ! ** update qs_env%rho 302 CALL qs_rho_update_rho(rho, qs_env=qs_env) 303 END IF 304 305 CALL section_vals_val_get(lr_section, "RESTART", l_val=linres_control%linres_restart) 306 CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter) 307 CALL section_vals_val_get(lr_section, "EPS", r_val=linres_control%eps) 308 CALL section_vals_val_get(lr_section, "EPS_FILTER", r_val=linres_control%eps_filter) 309 CALL section_vals_val_get(lr_section, "RESTART_EVERY", i_val=linres_control%restart_every) 310 CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type) 311 CALL section_vals_val_get(lr_section, "ENERGY_GAP", r_val=linres_control%energy_gap) 312 313 IF (iounit > 0) THEN 314 WRITE (UNIT=iounit, FMT="(/,T2,A,/,T25,A,/,T2,A,/)") & 315 REPEAT("=", 79), & 316 "START LINRES CALCULATION", & 317 REPEAT("=", 79) 318 319 WRITE (UNIT=iounit, FMT="(T2,A)") & 320 "LINRES| Properties to be calulated:" 321 CALL section_vals_val_get(lr_section, "NMR%_SECTION_PARAMETERS_", l_val=do_it) 322 IF (do_it) WRITE (UNIT=iounit, FMT="(T62,A)") "NMR Chemical Shift" 323 CALL section_vals_val_get(lr_section, "EPR%_SECTION_PARAMETERS_", l_val=do_it) 324 IF (do_it) WRITE (UNIT=iounit, FMT="(T68,A)") "EPR g Tensor" 325 CALL section_vals_val_get(lr_section, "SPINSPIN%_SECTION_PARAMETERS_", l_val=do_it) 326 IF (do_it) WRITE (UNIT=iounit, FMT="(T43,A)") "Indirect spin-spin coupling constants" 327 CALL section_vals_val_get(lr_section, "POLAR%_SECTION_PARAMETERS_", l_val=do_it) 328 IF (do_it) WRITE (UNIT=iounit, FMT="(T57,A)") "Electric Polarizability" 329 330 IF (linres_control%localized_psi0) WRITE (UNIT=iounit, FMT="(T2,A,T65,A)") & 331 "LINRES|", " LOCALIZED PSI0" 332 333 WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") & 334 "LINRES| Optimization algorithm", " Conjugate Gradients" 335 336 SELECT CASE (linres_control%preconditioner_type) 337 CASE (ot_precond_none) 338 WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") & 339 "LINRES| Preconditioner", " NONE" 340 CASE (ot_precond_full_single) 341 WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") & 342 "LINRES| Preconditioner", " FULL_SINGLE" 343 CASE (ot_precond_full_kinetic) 344 WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") & 345 "LINRES| Preconditioner", " FULL_KINETIC" 346 CASE (ot_precond_s_inverse) 347 WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") & 348 "LINRES| Preconditioner", " FULL_S_INVERSE" 349 CASE (ot_precond_full_single_inverse) 350 WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") & 351 "LINRES| Preconditioner", " FULL_SINGLE_INVERSE" 352 CASE (ot_precond_full_all) 353 WRITE (UNIT=iounit, FMT="(T2,A,T60,A)") & 354 "LINRES| Preconditioner", " FULL_ALL" 355 CASE DEFAULT 356 CPABORT("Preconditioner NYI") 357 END SELECT 358 359 WRITE (UNIT=iounit, FMT="(T2,A,T72,ES8.1)") & 360 "LINRES| EPS", linres_control%eps 361 WRITE (UNIT=iounit, FMT="(T2,A,T72,I8)") & 362 "LINRES| MAX_ITER", linres_control%max_iter 363 END IF 364 365 !------------------! 366 ! create the p_env ! 367 !------------------! 368 CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.TRUE., linres_control=linres_control) 369 370 ! update the m_epsilon matrix 371 CALL p_env_psi0_changed(p_env, qs_env) 372 373 p_env%os_valid = .FALSE. 374 p_env%new_preconditioner = .TRUE. 375 CALL cp_print_key_finished_output(iounit, logger, lr_section, & 376 "PRINT%PROGRAM_RUN_INFO") 377 378 END SUBROUTINE linres_init 379 380! ************************************************************************************************** 381!> \brief ... 382!> \param linres_control ... 383!> \param qs_env ... 384!> \param p_env ... 385!> \param dft_control ... 386!> \param nmr_present ... 387!> \param epr_present ... 388!> \param iounit ... 389! ************************************************************************************************** 390 SUBROUTINE nmr_epr_linres(linres_control, qs_env, p_env, dft_control, nmr_present, epr_present, iounit) 391 392 TYPE(linres_control_type), POINTER :: linres_control 393 TYPE(qs_environment_type), POINTER :: qs_env 394 TYPE(qs_p_env_type), POINTER :: p_env 395 TYPE(dft_control_type), POINTER :: dft_control 396 LOGICAL :: nmr_present, epr_present 397 INTEGER :: iounit 398 399 CHARACTER(LEN=*), PARAMETER :: routineN = 'nmr_epr_linres', routineP = moduleN//':'//routineN 400 401 INTEGER :: iB 402 LOGICAL :: do_qmmm 403 TYPE(current_env_type) :: current_env 404 TYPE(epr_env_type) :: epr_env 405 TYPE(nmr_env_type) :: nmr_env 406 407 linres_control%property = lr_current 408 409 CALL cite_reference(Weber2009) 410 411 IF (.NOT. linres_control%localized_psi0) THEN 412 CALL cp_abort(__LOCATION__, & 413 "Are you sure that you want to calculate the chemical "// & 414 "shift without localized psi0?") 415 CALL linres_localize(qs_env, linres_control, & 416 dft_control%nspins, centers_only=.TRUE.) 417 ENDIF 418 IF (dft_control%nspins /= 2 .AND. epr_present) THEN 419 CPABORT("LSD is needed to perform a g tensor calculation!") 420 ENDIF 421 ! 422 !Initialize the current environment 423 do_qmmm = .FALSE. 424 current_env%ref_count = 0 425 IF (qs_env%qmmm) do_qmmm = .TRUE. 426 current_env%do_qmmm = do_qmmm 427 !current_env%prop='nmr' 428 CALL current_env_init(current_env, qs_env) 429 CALL current_operators(current_env, qs_env) 430 CALL current_response(current_env, p_env, qs_env) 431 ! 432 IF (current_env%all_pert_op_done) THEN 433 !Initialize the nmr environment 434 IF (nmr_present) THEN 435 nmr_env%ref_count = 0 436 CALL nmr_env_init(nmr_env, qs_env) 437 ENDIF 438 ! 439 !Initialize the epr environment 440 IF (epr_present) THEN 441 epr_env%ref_count = 0 442 CALL epr_env_init(epr_env, qs_env) 443 CALL epr_g_zke(epr_env, qs_env) 444 CALL epr_nablavks(epr_env, qs_env) 445 ENDIF 446 ! 447 ! Build the rs_gauge if needed 448 !CALL current_set_gauge(current_env,qs_env) 449 ! 450 ! Loop over field direction 451 DO iB = 1, 3 452 ! 453 ! Build current response and succeptibility 454 CALL current_build_current(current_env, qs_env, iB) 455 CALL current_build_chi(current_env, qs_env, iB) 456 ! 457 ! Compute NMR shift 458 IF (nmr_present) THEN 459 CALL nmr_shift(nmr_env, current_env, qs_env, iB) 460 ENDIF 461 ! 462 ! Compute EPR 463 IF (epr_present) THEN 464 CALL epr_ind_magnetic_field(epr_env, current_env, qs_env, iB) 465 CALL epr_g_so(epr_env, current_env, qs_env, iB) 466 CALL epr_g_soo(epr_env, current_env, qs_env, iB) 467 ENDIF 468 ENDDO 469 ! 470 ! Finalized the nmr environment 471 IF (nmr_present) THEN 472 CALL nmr_shift_print(nmr_env, current_env, qs_env) 473 CALL nmr_env_cleanup(nmr_env) 474 ENDIF 475 ! 476 ! Finalized the epr environment 477 IF (epr_present) THEN 478 CALL epr_g_print(epr_env, qs_env) 479 CALL epr_env_cleanup(epr_env) 480 ENDIF 481 ! 482 ELSE 483 IF (iounit > 0) THEN 484 WRITE (iounit, "(T10,A,/T20,A,/)") & 485 "CURRENT: Not all responses to perturbation operators could be calculated.", & 486 " Hence: NO nmr and NO epr possible." 487 END IF 488 END IF 489 ! Finalized the current environment 490 CALL current_env_cleanup(current_env, qs_env) 491 492 END SUBROUTINE nmr_epr_linres 493 494! ************************************************************************************************** 495!> \brief ... 496!> \param linres_control ... 497!> \param qs_env ... 498!> \param p_env ... 499!> \param dft_control ... 500! ************************************************************************************************** 501 SUBROUTINE issc_linres(linres_control, qs_env, p_env, dft_control) 502 503 TYPE(linres_control_type), POINTER :: linres_control 504 TYPE(qs_environment_type), POINTER :: qs_env 505 TYPE(qs_p_env_type), POINTER :: p_env 506 TYPE(dft_control_type), POINTER :: dft_control 507 508 CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_linres', routineP = moduleN//':'//routineN 509 510 INTEGER :: iatom 511 LOGICAL :: do_qmmm 512 TYPE(current_env_type) :: current_env 513 TYPE(issc_env_type) :: issc_env 514 515 linres_control%property = lr_current 516 IF (.NOT. linres_control%localized_psi0) THEN 517 CALL cp_abort(__LOCATION__, & 518 "Are you sure that you want to calculate the chemical "// & 519 "shift without localized psi0?") 520 CALL linres_localize(qs_env, linres_control, & 521 dft_control%nspins, centers_only=.TRUE.) 522 ENDIF 523 ! 524 !Initialize the current environment 525 do_qmmm = .FALSE. 526 current_env%ref_count = 0 527 IF (qs_env%qmmm) do_qmmm = .TRUE. 528 current_env%do_qmmm = do_qmmm 529 !current_env%prop='issc' 530 !CALL current_env_init(current_env,qs_env) 531 !CALL current_response(current_env,p_env,qs_env) 532 ! 533 !Initialize the issc environment 534 issc_env%ref_count = 0 535 CALL issc_env_init(issc_env, qs_env) 536 ! 537 ! Loop over atoms 538 DO iatom = 1, issc_env%issc_natms 539 CALL issc_operators(issc_env, qs_env, iatom) 540 CALL issc_response(issc_env, p_env, qs_env) 541 CALL issc_issc(issc_env, qs_env, iatom) 542 ENDDO 543 ! 544 ! Finalized the issc environment 545 CALL issc_print(issc_env, qs_env) 546 CALL issc_env_cleanup(issc_env) 547 548 END SUBROUTINE issc_linres 549 550! ************************************************************************************************** 551!> \brief ... 552!> \param qs_env ... 553!> \param p_env ... 554!> \par History 555!> 06.2018 polar_env integrated into qs_env (MK) 556! ************************************************************************************************** 557 SUBROUTINE polar_linres(qs_env, p_env) 558 559 TYPE(qs_environment_type), POINTER :: qs_env 560 TYPE(qs_p_env_type), POINTER :: p_env 561 562 CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_linres', routineP = moduleN//':'//routineN 563 564 CALL polar_env_init(qs_env) 565 CALL polar_operators(qs_env) 566 CALL polar_response(p_env, qs_env) 567 CALL polar_polar(qs_env) 568 CALL polar_print(qs_env) 569 570 END SUBROUTINE polar_linres 571 572END MODULE qs_linres_module 573