From 4a6d2ced8a408732a7b2f15c6ae47eae766304eb Mon Sep 17 00:00:00 2001 From: Matthias Krack Date: Sun, 1 Jul 2018 16:55:45 +0000 Subject: [PATCH] Integrate polar_env into qs_env; enable linres restart for polar; new regression test input added svn-origin-rev: 18498 --- src/input_cp2k_motion_print.F | 11 +- src/motion/md_energies.F | 9 +- src/qs_environment_types.F | 31 +- src/qs_linres_methods.F | 25 +- src/qs_linres_module.F | 39 +- src/qs_linres_op.F | 22 +- src/qs_linres_polar_utils.F | 547 +++++++++++++----------- src/qs_linres_types.F | 119 +++++- tests/QS/regtest-polar/H2O_md_polar.inp | 108 +++++ tests/QS/regtest-polar/TEST_FILES | 5 +- tests/TEST_TYPES | 2 +- 11 files changed, 590 insertions(+), 328 deletions(-) create mode 100644 tests/QS/regtest-polar/H2O_md_polar.inp diff --git a/src/input_cp2k_motion_print.F b/src/input_cp2k_motion_print.F index 219543cfcc..22ba0a4ca1 100644 --- a/src/input_cp2k_motion_print.F +++ b/src/input_cp2k_motion_print.F @@ -155,7 +155,7 @@ CONTAINS CALL section_release(print_key) CALL cp_print_key_section_create(print_key, "SHELL_FORCES", & - description="controls the output of the forces on shells when shell-model is used", & + description="Controls the output of the forces on shells when shell-model is used", & print_level=high_print_level, common_iter_levels=1, & filename="", unit_str="hartree*bohr^-1") CALL add_format_keyword(keyword, print_key, pos=.FALSE., & @@ -164,7 +164,7 @@ CONTAINS CALL section_release(print_key) CALL cp_print_key_section_create(print_key, "CORE_FORCES", & - description="controls the output of the forces on cores when shell-model is used", & + description="Controls the output of the forces on cores when shell-model is used", & print_level=high_print_level, common_iter_levels=1, & filename="", unit_str="hartree*bohr^-1") CALL add_format_keyword(keyword, print_key, pos=.FALSE., & @@ -188,6 +188,13 @@ CONTAINS CALL section_add_subsection(section, print_key) CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, "POLAR_MATRIX", & + description="Controls the output of the polarisability tensor during an MD run", & + print_level=low_print_level, common_iter_levels=1, & + filename="") + CALL section_add_subsection(section, print_key) + CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, "RESTART", & description="Controls the dumping of the restart file during runs. "// & "By default keeps a short history of three restarts. See also RESTART_HISTORY", & diff --git a/src/motion/md_energies.F b/src/motion/md_energies.F index e368ae30b3..136d8aec6c 100644 --- a/src/motion/md_energies.F +++ b/src/motion/md_energies.F @@ -71,6 +71,7 @@ MODULE md_energies femtoseconds,& kelvin USE qmmm_types, ONLY: qmmm_env_type + USE qs_linres_polar_utils, ONLY: write_polarisability_tensor USE reftraj_types, ONLY: reftraj_type USE simpar_types, ONLY: simpar_type USE thermal_region_types, ONLY: thermal_regions_type @@ -459,8 +460,12 @@ CONTAINS CALL print_barostat_status(barostat, simpar, my_pos, my_act, cell, itimes, time) ! Print Stress Tensor - CALL write_stress_tensor(virial, cell, motion_section, itimes, time*femtoseconds, & - my_pos, my_act) + CALL write_stress_tensor(virial, cell, motion_section, itimes, time*femtoseconds, my_pos, my_act) + + ! Print Polarisability Tensor + IF (ASSOCIATED(force_env%qs_env)) THEN + CALL write_polarisability_tensor(force_env, motion_section, itimes, time*femtoseconds, my_pos, my_act) + END IF ! Temperature per Kinds IF (simpar%temperature_per_kind) THEN diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index 44bfc85cc0..87b5cb81e5 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -105,7 +105,9 @@ MODULE qs_environment_types set_ks_env USE qs_linres_types, ONLY: linres_control_release,& linres_control_retain,& - linres_control_type + linres_control_type,& + polar_env_release,& + polar_env_type USE qs_local_rho_types, ONLY: get_local_rho,& local_rho_set_create,& local_rho_set_release,& @@ -209,6 +211,7 @@ MODULE qs_environment_types !> \par History !> 11.2002 added doc and attribute description [fawzi] !> 08.2004 renamed some of the very short names (s,c,k,h) for easier grepping +!> 06.2018 polar_env added (MK) !> \author Matthias Krack & fawzi ! ************************************************************************************************** @@ -295,6 +298,8 @@ MODULE qs_environment_types ! Embedding potential TYPE(pw_p_type), POINTER :: embed_pot TYPE(pw_p_type), POINTER :: spin_embed_pot + ! Polarizability tensor + TYPE(polar_env_type), POINTER :: polar_env END TYPE qs_environment_type ! ************************************************************************************************** @@ -470,6 +475,7 @@ CONTAINS !> \param variable_history ... !> \param embed_pot ... !> \param spin_embed_pot ... +!> \param polar_env ... !> \date 23.01.2002 !> \author MK !> \version 1.0 @@ -500,7 +506,7 @@ CONTAINS lri_env, lri_density, dispersion_env, vee, rho_external, external_vxc, mask, & mp2_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, & s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, & - gradient_history, variable_history, embed_pot, spin_embed_pot) + gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env) TYPE(qs_environment_type), POINTER :: qs_env TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, & POINTER :: atomic_kind_set @@ -631,6 +637,7 @@ CONTAINS TYPE(almo_scf_env_type), OPTIONAL, POINTER :: almo_scf_env REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: gradient_history, variable_history TYPE(pw_p_type), OPTIONAL, POINTER :: embed_pot, spin_embed_pot + TYPE(polar_env_type), OPTIONAL, POINTER :: polar_env CHARACTER(len=*), PARAMETER :: routineN = 'get_qs_env', routineP = moduleN//':'//routineN @@ -696,6 +703,9 @@ CONTAINS IF (PRESENT(embed_pot)) embed_pot => qs_env%embed_pot IF (PRESENT(spin_embed_pot)) spin_embed_pot => qs_env%spin_embed_pot + ! Polarisability tensor + IF (PRESENT(polar_env)) polar_env => qs_env%polar_env + IF (PRESENT(rho_atom_set)) & CALL get_local_rho(qs_env%local_rho_set, rho_atom_set=rho_atom_set) IF (PRESENT(rho0_atom_set)) & @@ -860,6 +870,9 @@ CONTAINS NULLIFY (qs_env%embed_pot) NULLIFY (qs_env%spin_embed_pot) + ! Polarisability tensor + NULLIFY (qs_env%polar_env) + NULLIFY (qs_env%ks_env) NULLIFY (qs_env%ks_qmmm_env) NULLIFY (qs_env%wf_history) @@ -1000,6 +1013,7 @@ CONTAINS !> \param variable_history ... !> \param embed_pot ... !> \param spin_embed_pot ... +!> \param polar_env ... !> \date 23.01.2002 !> \author MK !> \version 1.0 @@ -1018,7 +1032,7 @@ CONTAINS se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, & do_transport, transport_env, lri_env, lri_density, dispersion_env, mp2_env, kg_env, & kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, & - spin_embed_pot) + spin_embed_pot, polar_env) TYPE(qs_environment_type), POINTER :: qs_env TYPE(cell_type), OPTIONAL, POINTER :: super_cell @@ -1086,6 +1100,7 @@ CONTAINS TYPE(almo_scf_env_type), OPTIONAL, POINTER :: almo_scf_env REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: gradient_history, variable_history TYPE(pw_p_type), OPTIONAL, POINTER :: embed_pot, spin_embed_pot + TYPE(polar_env_type), OPTIONAL, POINTER :: polar_env CHARACTER(len=*), PARAMETER :: routineN = 'set_qs_env', routineP = moduleN//':'//routineN @@ -1162,6 +1177,9 @@ CONTAINS IF (PRESENT(embed_pot)) qs_env%embed_pot => embed_pot IF (PRESENT(spin_embed_pot)) qs_env%spin_embed_pot => spin_embed_pot + ! Polarisability tensor + IF (PRESENT(polar_env)) qs_env%polar_env => polar_env + IF (PRESENT(qs_charges)) THEN CALL qs_charges_retain(qs_charges) CALL qs_charges_release(qs_env%qs_charges) @@ -1256,6 +1274,7 @@ CONTAINS IF (PRESENT(dispersion_env)) qs_env%dispersion_env => dispersion_env IF (PRESENT(WannierCentres)) qs_env%WannierCentres => WannierCentres IF (PRESENT(kpoints)) CALL set_ks_env(qs_env%ks_env, kpoints=kpoints) + END SUBROUTINE set_qs_env ! ************************************************************************************************** @@ -1298,6 +1317,7 @@ CONTAINS !> \param qs_env the object to release !> \par History !> 12.2002 created [fawzi] +!> 06.2018 polar_env added (MK) !> \author Fawzi Mohamed ! ************************************************************************************************** SUBROUTINE qs_env_release(qs_env) @@ -1424,6 +1444,9 @@ CONTAINS ENDIF ENDIF + ! Polarisability tensor + CALL polar_env_release(qs_env%polar_env) + CALL qs_charges_release(qs_env%qs_charges) CALL qs_ks_release(qs_env%ks_env) CALL qs_ks_qmmm_release(qs_env%ks_qmmm_env) @@ -1486,7 +1509,9 @@ CONTAINS DEALLOCATE (qs_env) END IF END IF + NULLIFY (qs_env) + END SUBROUTINE qs_env_release END MODULE qs_environment_types diff --git a/src/qs_linres_methods.F b/src/qs_linres_methods.F index 8ee6cf45bc..4ef8dc4a12 100644 --- a/src/qs_linres_methods.F +++ b/src/qs_linres_methods.F @@ -275,9 +275,6 @@ CONTAINS TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos TYPE(section_vals_type), POINTER :: lr_section -! -! - CALL timeset(routineN, handle) NULLIFY (dft_control, linres_control, matrix_s, matrix_t, matrix_ks, para_env) @@ -296,7 +293,6 @@ CONTAINS para_env=para_env, & mos=mos) - ! nspins = dft_control%nspins CALL get_mo_set(mos(1)%mo_set, nao=nao) maxnmo = 0 @@ -591,7 +587,7 @@ CONTAINS tr_rz0(ispin) = tr_rz1(ispin) ENDDO - ! ** can we exit the SCF loop ? + ! Can we exit the SCF loop? CALL external_control(should_stop, "LINRES", target_time=qs_env%target_time, & start_time=qs_env%start_time) @@ -1510,9 +1506,8 @@ CONTAINS extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.) IF (output_unit > 0) THEN - WRITE (UNIT=output_unit, FMT="(/,T10,A,A,/)") & - "Linres response functions are written to restart file", & - TRIM(filename) + WRITE (UNIT=output_unit, FMT="(T2,A)") & + "LINRES| Writing response functions to the restart file <"//TRIM(ADJUSTL(filename))//">" END IF ! @@ -1639,11 +1634,11 @@ CONTAINS file_status="OLD", & unit_number=rst_unit) - IF (output_unit > 0) WRITE (output_unit, "(/,T20,A,A)") & - "Read response wavefunction from restart file ", TRIM(filename) + IF (output_unit > 0) WRITE (output_unit, "(T2,A)") & + "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">" ELSE - IF (output_unit > 0) WRITE (output_unit, "(/,T10,A)") & - "Restart file not available filename=<"//TRIM(filename)//'>' + IF (output_unit > 0) WRITE (output_unit, "(T2,A)") & + "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found" ENDIF ENDIF @@ -1706,9 +1701,9 @@ CONTAINS IF (ivec .EQ. ivec_tmp) EXIT ENDDO - IF (iostat .NE. 0) THEN - IF (output_unit > 0) WRITE (output_unit, "(/,T10,A,A,/)") & - "Restart file: didnt find ", TRIM(filename) + IF (iostat /= 0) THEN + IF (output_unit > 0) WRITE (output_unit, "(T2,A)") & + "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found" ENDIF DEALLOCATE (vecbuffer) diff --git a/src/qs_linres_module.F b/src/qs_linres_module.F index 7e262e910a..c5710eea30 100644 --- a/src/qs_linres_module.F +++ b/src/qs_linres_module.F @@ -67,8 +67,7 @@ MODULE qs_linres_module USE qs_linres_op, ONLY: current_operators,& issc_operators,& polar_operators - USE qs_linres_polar_utils, ONLY: polar_env_cleanup,& - polar_env_init,& + USE qs_linres_polar_utils, ONLY: polar_env_init,& polar_polar,& polar_print,& polar_response @@ -78,8 +77,7 @@ MODULE qs_linres_module linres_control_create,& linres_control_release,& linres_control_type,& - nmr_env_type,& - polar_env_type + nmr_env_type USE qs_mo_methods, ONLY: calculate_density_matrix USE qs_mo_types, ONLY: mo_set_p_type USE qs_p_env_methods, ONLY: p_env_create,& @@ -130,7 +128,7 @@ CONTAINS CASE (use_qmmm) qs_env => force_env%qmmm_env%qs_env CASE DEFAULT - CPABORT("Doesnt recognize this force_env.") + CPABORT("Does not recognize this force_env") END SELECT qs_env%linres_run = .TRUE. @@ -170,6 +168,7 @@ CONTAINS TYPE(section_vals_type), POINTER :: lr_section, prop_section CALL timeset(routineN, handle) + lr_calculation = .FALSE. nmr_present = .FALSE. epr_present = .FALSE. @@ -189,9 +188,9 @@ CONTAINS CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & linres_control=linres_control) - !The type of perturbation has not been defined yet + ! The type of perturbation has not been defined yet linres_control%property = lr_none - ! + ! We do NMR or EPR, then compute the current response prop_section => section_vals_get_subs_vals(lr_section, "NMR") CALL section_vals_get(prop_section, explicit=nmr_present) @@ -202,7 +201,7 @@ CONTAINS CALL nmr_epr_linres(linres_control, qs_env, p_env, dft_control, & nmr_present, epr_present, output_unit) ENDIF - ! + ! We do the indirect spin-spin coupling calculation prop_section => section_vals_get_subs_vals(lr_section, "SPINSPIN") CALL section_vals_get(prop_section, explicit=issc_present) @@ -211,11 +210,9 @@ CONTAINS CALL issc_linres(linres_control, qs_env, p_env, dft_control) ENDIF - ! ! We do the polarizability calculation prop_section => section_vals_get_subs_vals(lr_section, "POLAR") CALL section_vals_get(prop_section, explicit=polar_present) - IF (polar_present) THEN CALL polar_linres(qs_env, p_env) END IF @@ -326,8 +323,8 @@ CONTAINS "START LINRES CALCULATION", & REPEAT("=", 77) - WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") & - "Properties to be Calulated:" + WRITE (UNIT=output_unit, FMT="(/,T2,A,/)") & + "LINRES| Properties to be calulated:" CALL section_vals_val_get(lr_section, "NMR%_SECTION_PARAMETERS_", & l_val=do_it) IF (do_it) WRITE (UNIT=output_unit, FMT="(T45,A)") & @@ -558,6 +555,8 @@ CONTAINS !> \brief ... !> \param qs_env ... !> \param p_env ... +!> \par History +!> 06.2018 polar_env integrated into qs_env (MK) ! ************************************************************************************************** SUBROUTINE polar_linres(qs_env, p_env) @@ -566,18 +565,12 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_linres', routineP = moduleN//':'//routineN - TYPE(polar_env_type) :: polar_env - - polar_env%ref_count = 0 - CALL polar_env_init(polar_env, qs_env) - CALL polar_operators(polar_env, qs_env) - CALL polar_response(polar_env, p_env, qs_env) - CALL polar_polar(polar_env, qs_env) - ! - CALL polar_print(polar_env, qs_env) - CALL polar_env_cleanup(polar_env) + CALL polar_env_init(qs_env) + CALL polar_operators(qs_env) + CALL polar_response(p_env, qs_env) + CALL polar_polar(qs_env) + CALL polar_print(qs_env) END SUBROUTINE polar_linres END MODULE qs_linres_module - diff --git a/src/qs_linres_op.F b/src/qs_linres_op.F index c893efecf5..fe2ceded36 100644 --- a/src/qs_linres_op.F +++ b/src/qs_linres_op.F @@ -629,16 +629,15 @@ CONTAINS !> the ground state wave function (see paper Putrino et al., JCP, 13, 7102) for the AOs; !> afterwards multiply with the ground state MO coefficients !> -!> \param polar_env ... !> \param qs_env ... !> \par History !> 01.2013 created [SL] +!> 06.2018 polar_env integrated into qs_env (MK) !> \author SL ! ************************************************************************************************** - SUBROUTINE polar_operators(polar_env, qs_env) + SUBROUTINE polar_operators(qs_env) - TYPE(polar_env_type) :: polar_env TYPE(qs_environment_type), POINTER :: qs_env CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_operators', & @@ -664,12 +663,14 @@ CONTAINS TYPE(dbcsr_type), POINTER :: cosmat, sinmat TYPE(dft_control_type), POINTER :: dft_control TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos + TYPE(polar_env_type), POINTER :: polar_env TYPE(section_vals_type), POINTER :: lr_section CALL timeset(routineN, handle) NULLIFY (dBerry_psi0, inv_work, mo_derivs, & op_fm_set, opvec, sinmat, cosmat) + NULLIFY (polar_env) logger => cp_get_default_logger() lr_section => section_vals_get_subs_vals(qs_env%input, & @@ -682,6 +683,7 @@ CONTAINS cell=cell, & dft_control=dft_control, & para_env=para_env, & + polar_env=polar_env, & mos=mos, & matrix_s=matrix_s) @@ -701,9 +703,6 @@ CONTAINS ENDDO END DO - ! - ! - ! initialize all work matrices needed ALLOCATE (opvec(2, dft_control%nspins)) ALLOCATE (op_fm_set(2, dft_control%nspins)) @@ -728,11 +727,9 @@ CONTAINS CALL cp_cfm_create(inv_mat(i, ispin)%matrix, op_fm_set(1, ispin)%matrix%matrix_struct) CALL cp_fm_create(inv_work(2, i, ispin)%matrix, op_fm_set(2, ispin)%matrix%matrix_struct) CALL cp_fm_create(inv_work(1, i, ispin)%matrix, op_fm_set(1, ispin)%matrix%matrix_struct) - END DO END DO - ! NULLIFY (cosmat, sinmat) ALLOCATE (cosmat, sinmat) CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM') @@ -758,7 +755,7 @@ CONTAINS ENDDO - !second step invert C^T S_berry C + ! Second step invert C^T S_berry C zdet = one DO ispin = 1, dft_control%nspins CALL cp_cfm_get_info(eigrmat(ispin)%matrix, ncol_local=tmp_dim) @@ -770,9 +767,8 @@ CONTAINS CALL cp_cfm_set_all(inv_mat(i, ispin)%matrix, zero, one) CALL cp_cfm_solve(eigrmat(ispin)%matrix, inv_mat(i, ispin)%matrix, zdeta) END DO - ! - ! - !compute the derivative and add the result to mo_derivatives + + ! Compute the derivative and add the result to mo_derivatives DO ispin = 1, dft_control%nspins CALL cp_cfm_get_info(eigrmat(ispin)%matrix, ncol_local=tmp_dim) CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo, maxocc=maxocc) @@ -816,7 +812,7 @@ CONTAINS DEALLOCATE (eigrmat) DEALLOCATE (inv_mat) - ENDIF ! raman + ENDIF ! do_raman CALL cp_print_key_finished_output(output_unit, logger, lr_section, & "PRINT%PROGRAM_RUN_INFO") diff --git a/src/qs_linres_polar_utils.F b/src/qs_linres_polar_utils.F index 81b3be30af..d2d1ffebba 100644 --- a/src/qs_linres_polar_utils.F +++ b/src/qs_linres_polar_utils.F @@ -17,8 +17,7 @@ MODULE qs_linres_polar_utils cite_reference USE cell_types, ONLY: cell_type USE cp_control_types, ONLY: dft_control_type - USE cp_fm_basic_linalg, ONLY: cp_fm_frobenius_norm,& - cp_fm_trace + USE cp_fm_basic_linalg, ONLY: cp_fm_trace USE cp_fm_struct, ONLY: cp_fm_struct_create,& cp_fm_struct_release,& cp_fm_struct_type @@ -41,20 +40,27 @@ MODULE qs_linres_polar_utils put_results USE cp_result_types, ONLY: cp_result_type USE dbcsr_api, ONLY: dbcsr_p_type + USE force_env_types, ONLY: force_env_get,& + force_env_type USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get USE kinds, ONLY: default_string_length,& dp + USE machine, ONLY: m_flush USE mathconstants, ONLY: twopi USE physcon, ONLY: angstrom USE qs_environment_types, ONLY: get_qs_env,& - qs_environment_type - USE qs_linres_methods, ONLY: linres_solver + qs_environment_type,& + set_qs_env + USE qs_linres_methods, ONLY: linres_read_restart,& + linres_solver,& + linres_write_restart USE qs_linres_types, ONLY: get_polar_env,& linres_control_type,& polar_env_create,& - polar_env_type + polar_env_type,& + set_polar_env USE qs_matrix_pools, ONLY: qs_matrix_pools_type USE qs_mo_types, ONLY: get_mo_set,& mo_set_p_type @@ -64,7 +70,8 @@ MODULE qs_linres_polar_utils IMPLICIT NONE PRIVATE - PUBLIC :: polar_env_cleanup, polar_env_init, polar_response, polar_polar, polar_print + + PUBLIC :: polar_env_init, polar_polar, polar_print, polar_response, write_polarisability_tensor CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_polar_utils' @@ -72,161 +79,116 @@ CONTAINS ! ************************************************************************************************** !> \brief Initialize the polar environment -!> \param polar_env ... -!> \param p_env ... !> \param qs_env ... +!> \par History +!> 06.2018 polar_env integrated into qs_env (MK) ! ************************************************************************************************** - SUBROUTINE polar_response(polar_env, p_env, qs_env) + SUBROUTINE polar_env_init(qs_env) - TYPE(polar_env_type) :: polar_env - TYPE(qs_p_env_type), POINTER :: p_env TYPE(qs_environment_type), POINTER :: qs_env - CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_response', routineP = moduleN//':'//routineN + CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_env_init', routineP = moduleN//':'//routineN - INTEGER :: handle, idir, ispin, nao, nmo, nspins, & - output_unit - LOGICAL :: do_raman, should_stop - REAL(dp) :: chk, fro - TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: h1_psi0, psi0_order, psi1 - TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry + INTEGER :: handle, idir, ispin, m, nao, nmo, & + nspins, output_unit TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(cp_logger_type), POINTER :: logger - TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s TYPE(dft_control_type), POINTER :: dft_control TYPE(linres_control_type), POINTER :: linres_control TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos - TYPE(qs_matrix_pools_type), POINTER :: mpools + TYPE(polar_env_type), POINTER :: polar_env TYPE(section_vals_type), POINTER :: lr_section, polar_section CALL timeset(routineN, handle) - ! - NULLIFY (dft_control, linres_control, lr_section, polar_section) - NULLIFY (logger, mpools, psi1, h1_psi0, mo_coeff, para_env) - NULLIFY (tmp_fm_struct, psi1_dBerry, dBerry_psi0) + + NULLIFY (dft_control) + NULLIFY (linres_control) + NULLIFY (logger) + NULLIFY (matrix_s) + NULLIFY (mos) + NULLIFY (polar_env) + NULLIFY (lr_section, polar_section) logger => cp_get_default_logger() lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") - polar_section => section_vals_get_subs_vals(qs_env%input, & - "PROPERTIES%LINRES%POLAR") output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & extension=".linresLog") + IF (output_unit > 0) THEN - WRITE (UNIT=output_unit, FMT="(T10,A,/)") & - "*** Self consistent optimization of the response wavefunctions ***" + WRITE (output_unit, "(/,(T2,A))") "POLAR| Starting polarizability calculation", & + "POLAR| Initialization of the polar environment" ENDIF + polar_section => section_vals_get_subs_vals(qs_env%input, & + "PROPERTIES%LINRES%POLAR") + CALL get_qs_env(qs_env=qs_env, & + polar_env=polar_env, & dft_control=dft_control, & - mpools=mpools, & + matrix_s=matrix_s, & linres_control=linres_control, & - mos=mos, & - para_env=para_env) + mos=mos) + + ! Create polar environment if needed + IF (.NOT. ASSOCIATED(polar_env)) THEN + CALL polar_env_create(polar_env) + CALL set_qs_env(qs_env=qs_env, polar_env=polar_env) + ENDIF nspins = dft_control%nspins - CALL get_polar_env(polar_env=polar_env, & - do_raman=do_raman) - ! - ! allocate the vectors - ALLOCATE (psi0_order(nspins)) - ALLOCATE (psi1(nspins), h1_psi0(nspins)) - DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) - psi0_order(ispin)%matrix => mo_coeff - CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao) - NULLIFY (tmp_fm_struct, psi1(ispin)%matrix, h1_psi0(ispin)%matrix) - CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, & - ncol_global=nmo, & - context=mo_coeff%matrix_struct%context) - CALL cp_fm_create(psi1(ispin)%matrix, tmp_fm_struct) - CALL cp_fm_create(h1_psi0(ispin)%matrix, tmp_fm_struct) - CALL cp_fm_struct_release(tmp_fm_struct) - ENDDO - chk = 0.0_dp - ! - IF (do_raman) THEN - CALL get_polar_env(polar_env=polar_env, & - psi1_dBerry=psi1_dBerry, & - dBerry_psi0=dBerry_psi0) - DO idir = 1, 3 - IF (output_unit > 0) THEN - WRITE (output_unit, "(T10,A)") & - "Response to the perturbation operator Berry phase_"//ACHAR(idir+119) - ENDIF - DO ispin = 1, nspins - CALL cp_fm_set_all(psi1_dBerry(idir, ispin)%matrix, 0.0_dp) - ENDDO - ! - !Initial guess for psi1 - DO ispin = 1, nspins - CALL cp_fm_set_all(psi1(ispin)%matrix, 0.0_dp) - ENDDO - ! - !DO scf cycle to optimize psi1 - DO ispin = 1, nspins - CALL cp_fm_to_fm(dBerry_psi0(idir, ispin)%matrix, h1_psi0(ispin)%matrix) - ENDDO - ! - ! - linres_control%lr_triplet = .FALSE. ! we do singlet response - linres_control%do_kernel = .TRUE. ! we do coupled response - linres_control%converged = .FALSE. - CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, should_stop) - ! - ! - ! copy the response - DO ispin = 1, nspins - CALL cp_fm_to_fm(psi1(ispin)%matrix, psi1_dBerry(idir, ispin)%matrix) - CALL cp_fm_frobenius_norm(psi1(ispin)%matrix, fro) - chk = chk+fro - ENDDO - ! - IF (output_unit > 0) THEN - WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet" - ENDIF - ! - ! Write the result in the restart file - ENDDO ! idir - ENDIF !raman + CALL section_vals_val_get(polar_section, "DO_RAMAN", l_val=polar_env%do_raman) - ! - ! print the checksum - ! IF(output_unit>0) THEN - ! WRITE(output_unit,'(T2,A,E23.16)') 'polar| response: CheckSum =',chk - ! ENDIF - ! - ! - ! clean up + ! Allocate components of the polar environment if needed + IF (.NOT. ASSOCIATED(polar_env%polar)) THEN + ALLOCATE (polar_env%polar(3, 3)) + polar_env%polar(:, :) = 0.0_dp + ENDIF + IF (.NOT. ASSOCIATED(polar_env%dBerry_psi0)) ALLOCATE (polar_env%dBerry_psi0(3, nspins)) + IF (.NOT. ASSOCIATED(polar_env%mo_derivs)) ALLOCATE (polar_env%mo_derivs(3, nspins)) + IF (.NOT. ASSOCIATED(polar_env%psi1_dBerry)) ALLOCATE (polar_env%psi1_dBerry(3, nspins)) DO ispin = 1, nspins - CALL cp_fm_release(psi1(ispin)%matrix) - CALL cp_fm_release(h1_psi0(ispin)%matrix) - ENDDO - DEALLOCATE (psi1, h1_psi0, psi0_order) - ! + CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nao=nao, nmo=nmo) + CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao) + NULLIFY (tmp_fm_struct) + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, & + ncol_global=m, & + context=mo_coeff%matrix_struct%context) + DO idir = 1, 3 + NULLIFY (polar_env%dBerry_psi0(idir, ispin)%matrix, & + polar_env%mo_derivs(idir, ispin)%matrix, & + polar_env%psi1_dBerry(idir, ispin)%matrix) + CALL cp_fm_create(polar_env%dBerry_psi0(idir, ispin)%matrix, tmp_fm_struct) + CALL cp_fm_create(polar_env%mo_derivs(idir, ispin)%matrix, tmp_fm_struct) + CALL cp_fm_create(polar_env%psi1_dBerry(idir, ispin)%matrix, tmp_fm_struct) + ENDDO + CALL cp_fm_struct_release(tmp_fm_struct) + END DO + CALL cp_print_key_finished_output(output_unit, logger, lr_section, & "PRINT%PROGRAM_RUN_INFO") - ! + CALL timestop(handle) - ! - END SUBROUTINE polar_response + + END SUBROUTINE polar_env_init ! ************************************************************************************************** !> \brief ... -!> \param polar_env ... !> \param qs_env ... +!> \par History +!> 06.2018 polar_env integrated into qs_env (MK) ! ************************************************************************************************** - SUBROUTINE polar_polar(polar_env, qs_env) + SUBROUTINE polar_polar(qs_env) - TYPE(polar_env_type) :: polar_env TYPE(qs_environment_type), POINTER :: qs_env CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_polar', routineP = moduleN//':'//routineN INTEGER :: handle, i, ispin, nspins, output_unit, z - LOGICAL :: do_raman + LOGICAL :: do_raman, run_stopped REAL(dp) :: ptmp REAL(dp), DIMENSION(:, :), POINTER :: polar, polar_tmp TYPE(cell_type), POINTER :: cell @@ -234,6 +196,7 @@ CONTAINS TYPE(cp_logger_type), POINTER :: logger TYPE(dft_control_type), POINTER :: dft_control TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos + TYPE(polar_env_type), POINTER :: polar_env CALL timeset(routineN, handle) @@ -245,27 +208,28 @@ CONTAINS CALL get_qs_env(qs_env=qs_env, & cell=cell, & dft_control=dft_control, & - mos=mos) + mos=mos, & + polar_env=polar_env) nspins = dft_control%nspins CALL get_polar_env(polar_env=polar_env, & - do_raman=do_raman) -! - IF (do_raman) THEN + do_raman=do_raman, & + run_stopped=run_stopped) + + IF (.NOT. run_stopped .AND. do_raman) THEN CALL cite_reference(Luber2014) CALL get_polar_env(polar_env=polar_env, & - psi1_dBerry=psi1_dBerry, & dBerry_psi0=dBerry_psi0, & + mo_derivs=mo_derivs, & polar=polar, & - mo_derivs=mo_derivs) - ! - ! + psi1_dBerry=psi1_dBerry) + ! polar_section => section_vals_get_subs_vals(qs_env%input, & ! & "PROPERTIES%LINRES%POLAR") - ! + ! Initialize ALLOCATE (polar_tmp(3, 3)) polar_tmp(:, :) = 0.0_dp @@ -284,60 +248,68 @@ CONTAINS polar(:, :) = MATMUL(MATMUL(cell%hmat, polar_tmp), TRANSPOSE(cell%hmat)) !SL evtl maxocc instead? - IF (dft_control%nspins .EQ. 1) THEN + IF (dft_control%nspins == 1) THEN polar(:, :) = 2.0_dp*polar(:, :) END IF - ! IF (ASSOCIATED(polar_tmp)) THEN DEALLOCATE (polar_tmp) END IF - ENDIF !raman - ! + ENDIF ! do_raman + ! Finalize CALL timestop(handle) - ! + END SUBROUTINE polar_polar ! ************************************************************************************************** -!> \brief ... -!> \param polar_env ... +!> \brief Print information related to the polarisability tensor !> \param qs_env ... +!> \par History +!> 06.2018 polar_env integrated into qs_env (MK) ! ************************************************************************************************** - SUBROUTINE polar_print(polar_env, qs_env) - TYPE(polar_env_type) :: polar_env + SUBROUTINE polar_print(qs_env) + TYPE(qs_environment_type), POINTER :: qs_env CHARACTER(len=*), PARAMETER :: routineN = 'polar_print', routineP = moduleN//':'//routineN CHARACTER(LEN=default_string_length) :: description INTEGER :: output_unit, unit_p - LOGICAL :: do_raman + LOGICAL :: do_raman, run_stopped REAL(dp), DIMENSION(:, :), POINTER :: polar TYPE(cp_logger_type), POINTER :: logger TYPE(cp_para_env_type), POINTER :: para_env TYPE(cp_result_type), POINTER :: results TYPE(dft_control_type), POINTER :: dft_control + TYPE(polar_env_type), POINTER :: polar_env TYPE(section_vals_type), POINTER :: polar_section NULLIFY (logger, dft_control, para_env, results) - CALL get_qs_env(qs_env=qs_env, results=results, dft_control=dft_control, & + + CALL get_qs_env(qs_env=qs_env, & + dft_control=dft_control, & + polar_env=polar_env, & + results=results, & para_env=para_env) logger => cp_get_default_logger() output_unit = cp_logger_get_default_io_unit(logger) - polar_section => section_vals_get_subs_vals(qs_env%input, & - "PROPERTIES%LINRES%POLAR") + polar_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%POLAR") - CALL get_polar_env(polar_env=polar_env, polar=polar, do_raman=do_raman) + CALL get_polar_env(polar_env=polar_env, & + polar=polar, & + do_raman=do_raman, & + run_stopped=run_stopped) - description = "[POLAR]" - CALL cp_results_erase(results, description=description) - CALL put_results(results, description=description, values=polar(:, :)) + IF (.NOT. run_stopped .AND. do_raman) THEN + + description = "[POLAR]" + CALL cp_results_erase(results, description=description) + CALL put_results(results, description=description, values=polar(:, :)) - IF (do_raman) THEN IF (BTEST(cp_print_key_should_output(logger%iter_info, polar_section, & "PRINT%POLAR_MATRIX"), cp_p_file)) THEN @@ -364,166 +336,247 @@ CONTAINS ENDIF IF (output_unit > 0) THEN WRITE (output_unit, *) - WRITE (output_unit, '(T10,A)') 'POLARIZABILITY TENSOR (atomic units):' - WRITE (output_unit, '(T10,A,3F15.5)') "xx,yy,zz", polar(1, 1), polar(2, 2), polar(3, 3) - WRITE (output_unit, '(T10,A,3F15.5)') "xy,xz,yz", polar(1, 2), polar(1, 3), polar(2, 3) - WRITE (output_unit, '(T10,A,3F15.5)') "yx,zx,zy", polar(2, 1), polar(3, 1), polar(3, 2) - WRITE (output_unit, '(T10,A)') 'POLARIZABILITY TENSOR (Angstrom^3):' - WRITE (output_unit, '(T10,A,3F15.5)') "xx,yy,zz", polar(1, 1)*angstrom**3, & + WRITE (output_unit, '(T2,A)') 'POLAR| POLARIZABILITY TENSOR (atomic units):' + WRITE (output_unit, '(T2,A,3F15.5)') "POLAR| xx,yy,zz", polar(1, 1), polar(2, 2), polar(3, 3) + WRITE (output_unit, '(T2,A,3F15.5)') "POLAR| xy,xz,yz", polar(1, 2), polar(1, 3), polar(2, 3) + WRITE (output_unit, '(T2,A,3F15.5)') "POLAR| yx,zx,zy", polar(2, 1), polar(3, 1), polar(3, 2) + WRITE (output_unit, '(T2,A)') 'POLAR| POLARIZABILITY TENSOR (Angstrom^3):' + WRITE (output_unit, '(T2,A,3F15.5)') "POLAR| xx,yy,zz", polar(1, 1)*angstrom**3, & polar(2, 2)*angstrom**3, polar(3, 3)*angstrom**3 - WRITE (output_unit, '(T10,A,3F15.5)') "xy,xz,yz", polar(1, 2)*angstrom**3, & + WRITE (output_unit, '(T2,A,3F15.5)') "POLAR| xy,xz,yz", polar(1, 2)*angstrom**3, & polar(1, 3)*angstrom**3, polar(2, 3)*angstrom**3 - WRITE (output_unit, '(T10,A,3F15.5)') "yx,zx,zy", polar(2, 1)*angstrom**3, & + WRITE (output_unit, '(T2,A,3F15.5)') "POLAR| yx,zx,zy", polar(2, 1)*angstrom**3, & polar(3, 1)*angstrom**3, polar(3, 2)*angstrom**3 END IF END IF - ! - ! + END SUBROUTINE polar_print ! ************************************************************************************************** -!> \brief Initialize the polar environment -!> \param polar_env ... +!> \brief Calculate the polarisability tensor using response theory +!> \param p_env ... !> \param qs_env ... +!> \par History +!> 06.2018 polar_env integrated into qs_env (MK) ! ************************************************************************************************** - SUBROUTINE polar_env_init(polar_env, qs_env) - ! - TYPE(polar_env_type) :: polar_env + SUBROUTINE polar_response(p_env, qs_env) + + TYPE(qs_p_env_type), POINTER :: p_env TYPE(qs_environment_type), POINTER :: qs_env - CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_env_init', routineP = moduleN//':'//routineN + CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_response', routineP = moduleN//':'//routineN - INTEGER :: handle, idir, ispin, m, nao, nmo, & - nspins, output_unit + INTEGER :: handle, idir, ispin, nao, nmo, nspins, & + output_unit + LOGICAL :: do_raman, should_stop + TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: h1_psi0, psi0_order, psi1, psi1_ptr + TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(cp_logger_type), POINTER :: logger - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s + TYPE(cp_para_env_type), POINTER :: para_env TYPE(dft_control_type), POINTER :: dft_control TYPE(linres_control_type), POINTER :: linres_control TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos + TYPE(polar_env_type), POINTER :: polar_env + TYPE(qs_matrix_pools_type), POINTER :: mpools TYPE(section_vals_type), POINTER :: lr_section, polar_section -! +!MK REAL(dp) :: chk, fro CALL timeset(routineN, handle) - NULLIFY (linres_control) - NULLIFY (logger, polar_section) + NULLIFY (dft_control, linres_control, lr_section, polar_section) + NULLIFY (logger, mpools, psi1, h1_psi0, mo_coeff, para_env) + NULLIFY (tmp_fm_struct, psi1_dBerry, dBerry_psi0) logger => cp_get_default_logger() lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") - - output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & - extension=".linresLog") - - IF (polar_env%ref_count /= 0) THEN - CALL polar_env_cleanup(polar_env) - ENDIF - - IF (output_unit > 0) THEN - WRITE (output_unit, "(/,T20,A,/)") "*** Start polarizability calculation ***" - WRITE (output_unit, "(T10,A,/)") "Initialization of the polar environment" - ENDIF - - CALL polar_env_create(polar_env) - ! polar_section => section_vals_get_subs_vals(qs_env%input, & "PROPERTIES%LINRES%POLAR") + output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & + extension=".linresLog") + IF (output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT="(T2,A,/)") & + "POLAR| Self consistent optimization of the response wavefunctions" + ENDIF + CALL get_qs_env(qs_env=qs_env, & dft_control=dft_control, & - matrix_s=matrix_s, & + mpools=mpools, & linres_control=linres_control, & - mos=mos) - ! - ! + mos=mos, & + polar_env=polar_env, & + para_env=para_env) + nspins = dft_control%nspins - ! - CALL section_vals_val_get(polar_section, "DO_RAMAN", l_val=polar_env%do_raman) - ! - ! - ! - ALLOCATE (polar_env%polar(3, 3)) - polar_env%polar(:, :) = 0.0_dp - ! - ! allocation - ALLOCATE (polar_env%dBerry_psi0(3, nspins), & - polar_env%psi1_dBerry(3, nspins), polar_env%mo_derivs(3, nspins)) + + CALL get_polar_env(polar_env=polar_env, & + do_raman=do_raman) + + ! Allocate the vectors + ALLOCATE (psi0_order(nspins)) + ALLOCATE (psi1(nspins), h1_psi0(nspins)) DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nao=nao, nmo=nmo) - CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao) - - NULLIFY (tmp_fm_struct) + CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) + psi0_order(ispin)%matrix => mo_coeff + CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao) + NULLIFY (tmp_fm_struct, psi1(ispin)%matrix, h1_psi0(ispin)%matrix) CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, & - ncol_global=m, & + ncol_global=nmo, & context=mo_coeff%matrix_struct%context) - DO idir = 1, 3 - NULLIFY (polar_env%psi1_dBerry(idir, ispin)%matrix, polar_env%dBerry_psi0(idir, ispin)%matrix, & - polar_env%mo_derivs(idir, ispin)%matrix) - CALL cp_fm_create(polar_env%psi1_dBerry(idir, ispin)%matrix, tmp_fm_struct) - CALL cp_fm_create(polar_env%dBerry_psi0(idir, ispin)%matrix, tmp_fm_struct) - CALL cp_fm_create(polar_env%mo_derivs(idir, ispin)%matrix, tmp_fm_struct) - ENDDO + CALL cp_fm_create(psi1(ispin)%matrix, tmp_fm_struct) + CALL cp_fm_create(h1_psi0(ispin)%matrix, tmp_fm_struct) CALL cp_fm_struct_release(tmp_fm_struct) + ENDDO +!MK chk = 0.0_dp - ! - END DO + IF (do_raman) THEN + CALL get_polar_env(polar_env=polar_env, & + psi1_dBerry=psi1_dBerry, & + dBerry_psi0=dBerry_psi0) + ! Restart + IF (linres_control%linres_restart) THEN + DO idir = 1, 3 + psi1_ptr => psi1_dBerry(idir, :) + CALL linres_read_restart(qs_env, lr_section, psi1_ptr, idir, "psi1_dBerry") + ENDDO + ELSE + DO idir = 1, 3 + DO ispin = 1, nspins + CALL cp_fm_set_all(psi1_dBerry(idir, ispin)%matrix, 0.0_dp) + ENDDO + ENDDO + ENDIF + loop_idir: DO idir = 1, 3 + IF (output_unit > 0) THEN + WRITE (output_unit, "(/,T2,A)") & + "POLAR| Response to the perturbation operator Berry phase_"//ACHAR(idir+119) + ENDIF + ! Do scf cycle to optimize psi1 + DO ispin = 1, nspins + CALL cp_fm_to_fm(psi1_dBerry(idir, ispin)%matrix, psi1(ispin)%matrix) + CALL cp_fm_to_fm(dBerry_psi0(idir, ispin)%matrix, h1_psi0(ispin)%matrix) + ENDDO + ! + linres_control%lr_triplet = .FALSE. ! we do singlet response + linres_control%do_kernel = .TRUE. ! we do coupled response + linres_control%converged = .FALSE. + CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, should_stop) + + ! Copy the response + DO ispin = 1, nspins + CALL cp_fm_to_fm(psi1(ispin)%matrix, psi1_dBerry(idir, ispin)%matrix) +!MK CALL cp_fm_frobenius_norm(psi1(ispin)%matrix, fro) +!MK chk = chk+fro + ENDDO + ! +!MK IF (output_unit > 0) THEN +!MK WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet" +!MK ENDIF + ! + ! Write the new result to the restart file + IF (linres_control%linres_restart) THEN + psi1_ptr => psi1_dBerry(idir, :) + CALL linres_write_restart(qs_env, lr_section, psi1_ptr, idir, "psi1_dBerry") + ENDIF + ENDDO loop_idir + ENDIF ! do_raman + + CALL set_polar_env(polar_env, run_stopped=should_stop) + + ! print the checksum + ! IF(output_unit>0) THEN + ! WRITE(output_unit,'(T2,A,E23.16)') 'polar| response: CheckSum =',chk + ! ENDIF + + ! Clean up + DO ispin = 1, nspins + CALL cp_fm_release(psi1(ispin)%matrix) + CALL cp_fm_release(h1_psi0(ispin)%matrix) + ENDDO + DEALLOCATE (psi1, h1_psi0, psi0_order) CALL cp_print_key_finished_output(output_unit, logger, lr_section, & "PRINT%PROGRAM_RUN_INFO") CALL timestop(handle) - END SUBROUTINE polar_env_init + END SUBROUTINE polar_response ! ************************************************************************************************** -!> \brief Deallocate the polar environment -!> \param polar_env ... +!> \brief Prints the polarisability tensor to a file during MD runs +!> \param force_env ... +!> \param motion_section ... +!> \param itimes ... +!> \param time ... +!> \param pos ... +!> \param act ... !> \par History +!> 06.2018 Creation (MK) +!> \author Matthias Krack (MK) ! ************************************************************************************************** - SUBROUTINE polar_env_cleanup(polar_env) + SUBROUTINE write_polarisability_tensor(force_env, motion_section, itimes, time, pos, act) - TYPE(polar_env_type) :: polar_env + TYPE(force_env_type), POINTER :: force_env + TYPE(section_vals_type), POINTER :: motion_section + INTEGER, INTENT(IN) :: itimes + REAL(KIND=dp), INTENT(IN) :: time + CHARACTER(LEN=default_string_length), INTENT(IN), & + OPTIONAL :: pos, act - CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_env_cleanup', & + CHARACTER(LEN=*), PARAMETER :: routineN = 'write_polarisability_tensor', & routineP = moduleN//':'//routineN - INTEGER :: idir, ispin + CHARACTER(LEN=default_string_length) :: my_act, my_pos + INTEGER :: output_unit + LOGICAL :: do_raman, new_file, run_stopped + REAL(KIND=dp), DIMENSION(:, :), POINTER :: polar + TYPE(cp_logger_type), POINTER :: logger + TYPE(polar_env_type), POINTER :: polar_env + TYPE(qs_environment_type), POINTER :: qs_env - polar_env%ref_count = polar_env%ref_count-1 - IF (polar_env%ref_count == 0) THEN - IF (ASSOCIATED(polar_env%polar)) THEN - DEALLOCATE (polar_env%polar) - ENDIF - IF (ASSOCIATED(polar_env%dBerry_psi0)) THEN - DO idir = 1, SIZE(polar_env%dBerry_psi0, 1) - DO ispin = 1, SIZE(polar_env%dBerry_psi0, 2) - CALL cp_fm_release(polar_env%dBerry_psi0(idir, ispin)%matrix) - ENDDO - ENDDO - DEALLOCATE (polar_env%dBerry_psi0) - ENDIF - IF (ASSOCIATED(polar_env%mo_derivs)) THEN - DO idir = 1, SIZE(polar_env%mo_derivs, 1) - DO ispin = 1, SIZE(polar_env%mo_derivs, 2) - CALL cp_fm_release(polar_env%mo_derivs(idir, ispin)%matrix) - ENDDO - ENDDO - DEALLOCATE (polar_env%mo_derivs) - ENDIF - IF (ASSOCIATED(polar_env%psi1_dBerry)) THEN - DO idir = 1, SIZE(polar_env%psi1_dBerry, 1) - DO ispin = 1, SIZE(polar_env%psi1_dBerry, 2) - CALL cp_fm_release(polar_env%psi1_dBerry(idir, ispin)%matrix) - ENDDO - ENDDO - DEALLOCATE (polar_env%psi1_dBerry) - ENDIF - ! + NULLIFY (qs_env) - ENDIF ! ref count + CALL force_env_get(force_env, qs_env=qs_env) + IF (ASSOCIATED(qs_env)) THEN + NULLIFY (polar_env) + CALL get_qs_env(qs_env=qs_env, polar_env=polar_env) + IF (ASSOCIATED(polar_env)) THEN + CALL get_polar_env(polar_env=polar_env, & + polar=polar, & + do_raman=do_raman, & + run_stopped=run_stopped) + IF (.NOT. run_stopped .AND. do_raman) THEN + NULLIFY (logger) + logger => cp_get_default_logger() + my_pos = "APPEND" + my_act = "WRITE" + IF (PRESENT(pos)) my_pos = pos + IF (PRESENT(act)) my_act = act + output_unit = cp_print_key_unit_nr(logger, motion_section, "PRINT%POLAR_MATRIX", & + extension=".polar", file_position=my_pos, & + file_action=my_act, file_form="FORMATTED", & + is_new_file=new_file) + ELSE + output_unit = 0 + END IF + IF (output_unit > 0) THEN + IF (new_file) THEN + WRITE (UNIT=output_unit, FMT='(A,9(11X,A2," [a.u.]"),6X,A)') & + "# Step Time [fs]", "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz" + END IF + WRITE (UNIT=output_unit, FMT='(I8,F12.3,9(1X,F19.8))') itimes, time, & + polar(1, 1), polar(1, 2), polar(1, 3), & + polar(2, 1), polar(2, 2), polar(2, 3), & + polar(3, 1), polar(3, 2), polar(3, 3) + CALL m_flush(output_unit) + CALL cp_print_key_finished_output(output_unit, logger, motion_section, "PRINT%POLAR_MATRIX") + END IF + END IF ! polar_env + END IF ! qs_env - END SUBROUTINE polar_env_cleanup + END SUBROUTINE write_polarisability_tensor END MODULE qs_linres_polar_utils diff --git a/src/qs_linres_types.F b/src/qs_linres_types.F index ece0993be9..8c7f19fa29 100644 --- a/src/qs_linres_types.F +++ b/src/qs_linres_types.F @@ -15,7 +15,8 @@ MODULE qs_linres_types gto_basis_set_type USE cp_array_utils, ONLY: cp_2d_i_p_type,& cp_2d_r_p_type - USE cp_fm_types, ONLY: cp_fm_p_type + USE cp_fm_types, ONLY: cp_fm_p_type,& + cp_fm_release USE dbcsr_api, ONLY: dbcsr_p_type USE kinds, ONLY: dp USE qs_grid_atom, ONLY: grid_atom_type @@ -118,12 +119,11 @@ MODULE qs_linres_types END TYPE current_env_type !********************************************************************************************************* -! \param type for polarizability calculation using Berry operator - +! \param type for polarisability calculation using Berry operator TYPE polar_env_type - INTEGER :: ref_count - LOGICAL :: do_raman - REAL(dp), DIMENSION(:, :), POINTER :: polar + INTEGER :: ref_count + LOGICAL :: do_raman, run_stopped + REAL(dp), DIMENSION(:, :), POINTER :: polar TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: psi1_dBerry, dBerry_psi0, mo_derivs END TYPE polar_env_type ! ********************************************************************************************************* @@ -227,7 +227,7 @@ MODULE qs_linres_types set2zero_jrho_atom_rad, get_epr_env, & nmr_env_create, epr_env_create, current_env_create, get_issc_env, & set_current_env, issc_env_create, get_polar_env, & - polar_env_create + polar_env_create, polar_env_release, set_polar_env CONTAINS @@ -343,6 +343,7 @@ CONTAINS NULLIFY (current_env%psi0_order) END SUBROUTINE current_env_create + ! ************************************************************************************************** !> \brief ... !> \param nmr_env ... @@ -1243,7 +1244,7 @@ CONTAINS DO iat = 1, nat iatom = atom_list(iat) - !*** allocate the radial density for each LM,for each atom *** + ! Allocate the radial density for each LM,for each atom ALLOCATE (jrho1_atom_set(iatom)%jrho_vec_rad_h(3, nspins), & jrho1_atom_set(iatom)%jrho_vec_rad_s(3, nspins), & jrho1_atom_set(iatom)%jrho_h(nspins), & @@ -1381,20 +1382,23 @@ CONTAINS CALL timestop(handle) END SUBROUTINE init_nablavks_atom_set + ! ************************************************************************************************** !> \brief ... !> \param polar_env ... ! ************************************************************************************************** SUBROUTINE polar_env_create(polar_env) - TYPE(polar_env_type) :: polar_env + TYPE(polar_env_type), POINTER :: polar_env CHARACTER(len=*), PARAMETER :: routineN = 'polar_env_create', & routineP = moduleN//':'//routineN - CPASSERT(polar_env%ref_count == 0) + CPASSERT(.NOT. ASSOCIATED(polar_env)) + ALLOCATE (polar_env) polar_env%ref_count = 1 -! polar_env%do_raman=.FALSE. + polar_env%do_raman = .FALSE. + polar_env%run_stopped = .FALSE. NULLIFY (polar_env%polar) NULLIFY (polar_env%psi1_dBerry) NULLIFY (polar_env%dBerry_psi0) @@ -1410,45 +1414,120 @@ CONTAINS !> \param polar ... !> \param psi1_dBerry ... !> \param mo_derivs ... +!> \param run_stopped ... +!> \par History +!> 06.2018 polar_env integrated into qs_env (MK) ! ************************************************************************************************** SUBROUTINE get_polar_env(polar_env, do_raman, dBerry_psi0, polar, & - psi1_dBerry, mo_derivs) + psi1_dBerry, mo_derivs, run_stopped) - TYPE(polar_env_type) :: polar_env + TYPE(polar_env_type), POINTER :: polar_env LOGICAL, OPTIONAL :: do_raman TYPE(cp_fm_p_type), DIMENSION(:, :), OPTIONAL, & POINTER :: dBerry_psi0 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: polar TYPE(cp_fm_p_type), DIMENSION(:, :), OPTIONAL, & POINTER :: psi1_dBerry, mo_derivs + LOGICAL, OPTIONAL :: run_stopped CHARACTER(len=*), PARAMETER :: routineN = 'get_polar_env', routineP = moduleN//':'//routineN CPASSERT(polar_env%ref_count > 0) IF (PRESENT(polar)) polar => polar_env%polar - IF (PRESENT(psi1_dBerry)) psi1_dBerry => polar_env%psi1_dBerry + IF (PRESENT(do_raman)) do_raman = polar_env%do_raman IF (PRESENT(dBerry_psi0)) dBerry_psi0 => polar_env%dBerry_psi0 IF (PRESENT(mo_derivs)) mo_derivs => polar_env%mo_derivs - IF (PRESENT(do_raman)) do_raman = polar_env%do_raman - END SUBROUTINE get_polar_env + IF (PRESENT(psi1_dBerry)) psi1_dBerry => polar_env%psi1_dBerry + IF (PRESENT(run_stopped)) run_stopped = polar_env%run_stopped -! ************************************************************************************************** + END SUBROUTINE get_polar_env ! ************************************************************************************************** !> \brief ... !> \param polar_env ... +!> \param do_raman ... +!> \param dBerry_psi0 ... +!> \param polar ... +!> \param psi1_dBerry ... +!> \param mo_derivs ... +!> \param run_stopped ... ! ************************************************************************************************** - SUBROUTINE set_polar_env(polar_env) + SUBROUTINE set_polar_env(polar_env, do_raman, dBerry_psi0, polar, & + psi1_dBerry, mo_derivs, run_stopped) - TYPE(polar_env_type) :: polar_env + TYPE(polar_env_type), POINTER :: polar_env + LOGICAL, OPTIONAL :: do_raman + TYPE(cp_fm_p_type), DIMENSION(:, :), OPTIONAL, & + POINTER :: dBerry_psi0 + REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: polar + TYPE(cp_fm_p_type), DIMENSION(:, :), OPTIONAL, & + POINTER :: psi1_dBerry, mo_derivs + LOGICAL, OPTIONAL :: run_stopped CHARACTER(len=*), PARAMETER :: routineN = 'set_polar_env', routineP = moduleN//':'//routineN CPASSERT(polar_env%ref_count > 0) + IF (PRESENT(polar)) polar_env%polar => polar + IF (PRESENT(do_raman)) polar_env%do_raman = do_raman + IF (PRESENT(psi1_dBerry)) polar_env%psi1_dBerry => psi1_dBerry + IF (PRESENT(dBerry_psi0)) polar_env%dBerry_psi0 => dBerry_psi0 + IF (PRESENT(mo_derivs)) polar_env%mo_derivs => mo_derivs + IF (PRESENT(run_stopped)) polar_env%run_stopped = run_stopped + END SUBROUTINE set_polar_env + ! ************************************************************************************************** +!> \brief Deallocate the polar environment +!> \param polar_env ... +!> \par History +!> 06.2018 polar_env integrated into qs_env (MK) +! ************************************************************************************************** + SUBROUTINE polar_env_release(polar_env) + + TYPE(polar_env_type), POINTER :: polar_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_env_release', & + routineP = moduleN//':'//routineN + + INTEGER :: idir, ispin + + IF (ASSOCIATED(polar_env)) THEN + CPASSERT(polar_env%ref_count > 0) + polar_env%ref_count = polar_env%ref_count-1 + IF (polar_env%ref_count < 1) THEN + IF (ASSOCIATED(polar_env%polar)) THEN + DEALLOCATE (polar_env%polar) + ENDIF + IF (ASSOCIATED(polar_env%dBerry_psi0)) THEN + DO idir = 1, SIZE(polar_env%dBerry_psi0, 1) + DO ispin = 1, SIZE(polar_env%dBerry_psi0, 2) + CALL cp_fm_release(polar_env%dBerry_psi0(idir, ispin)%matrix) + ENDDO + ENDDO + DEALLOCATE (polar_env%dBerry_psi0) + ENDIF + IF (ASSOCIATED(polar_env%mo_derivs)) THEN + DO idir = 1, SIZE(polar_env%mo_derivs, 1) + DO ispin = 1, SIZE(polar_env%mo_derivs, 2) + CALL cp_fm_release(polar_env%mo_derivs(idir, ispin)%matrix) + ENDDO + ENDDO + DEALLOCATE (polar_env%mo_derivs) + ENDIF + IF (ASSOCIATED(polar_env%psi1_dBerry)) THEN + DO idir = 1, SIZE(polar_env%psi1_dBerry, 1) + DO ispin = 1, SIZE(polar_env%psi1_dBerry, 2) + CALL cp_fm_release(polar_env%psi1_dBerry(idir, ispin)%matrix) + ENDDO + ENDDO + DEALLOCATE (polar_env%psi1_dBerry) + ENDIF + ENDIF + NULLIFY (polar_env) + ENDIF + + END SUBROUTINE polar_env_release END MODULE qs_linres_types - diff --git a/tests/QS/regtest-polar/H2O_md_polar.inp b/tests/QS/regtest-polar/H2O_md_polar.inp new file mode 100644 index 0000000000..7c21759485 --- /dev/null +++ b/tests/QS/regtest-polar/H2O_md_polar.inp @@ -0,0 +1,108 @@ +&GLOBAL + PROJECT H2O_md_polar + PRINT_LEVEL low + RUN_TYPE md +&END GLOBAL + +&MOTION + &MD + ENSEMBLE nve + STEPS 2 + TEMPERATURE [K] 298.0 + TIMESTEP [fs] 0.50 + &END MD + &PRINT + &POLAR_MATRIX on + &EACH + MD 1 + &END EACH + &END POLAR_MATRIX + &RESTART on + BACKUP_COPIES 0 + &EACH + MD 1 + &END EACH + &END RESTART + &RESTART_HISTORY off + &END RESTART_HISTORY + &STRUCTURE_DATA + DISTANCE 1 2 + DISTANCE 1 3 + ANGLE 2 1 3 + &END STRUCTURE_DATA + &TRAJECTORY on + FORMAT xyz + &EACH + MD 1 + &END EACH + &END TRAJECTORY + &VELOCITIES off + &END VELOCITIES + &END PRINT +&END MOTION + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 320 + &END MGRID + &QS + EPS_DEFAULT 1.0E-12 + &END QS + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS atomic + MAX_SCF 30 + &PRINT + &RESTART + ADD_LAST NUMERIC + BACKUP_COPIES 0 + &END RESTART + &END PRINT + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &XC_GRID + XC_DERIV SPLINE2_SMOOTH + &END XC_GRID + &END XC + &END DFT + &PROPERTIES + &LINRES + EPS 1.0E-6 + MAX_ITER 30 + PRECONDITIONER full_single_inverse + RESTART on + &POLAR on + DO_RAMAN + &PRINT + &POLAR_MATRIX off + &END POLAR_MATRIX + &END PRINT + &END POLAR + &END LINRES + &END PROPERTIES + &SUBSYS + &CELL + ABC [bohr] 12.0 12.0 12.0 + &END CELL + &COORD + UNIT bohr + O 0.000000 0.000000 0.224953 + H 0.000000 1.451310 -0.899812 + H 0.000000 -1.451310 -0.899812 + &END COORD + &KIND H + BASIS_SET DZVP-MOLOPT-SR-GTH-q1 + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-MOLOPT-SR-GTH-q6 + POTENTIAL GTH-PBE-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-polar/TEST_FILES b/tests/QS/regtest-polar/TEST_FILES index dd0b4adce3..e449eb9ff9 100644 --- a/tests/QS/regtest-polar/TEST_FILES +++ b/tests/QS/regtest-polar/TEST_FILES @@ -3,6 +3,7 @@ # e.g. 0 means do not compare anything, running is enough # 1 compares the last total energy in the file # for details see cp2k/tools/do_regtest -h2o_LRraman.inp 60 2e-05 0.00613 -h2o_LRraman_noort.inp 60 2e-05 0.023300000000000001 +h2o_LRraman.inp 60 2e-05 0.00613 +h2o_LRraman_noort.inp 60 2e-05 0.02330 +H2O_md_polar.inp 60 1e-04 0.63555 #EOF diff --git a/tests/TEST_TYPES b/tests/TEST_TYPES index 27e7a215e8..1152b03e10 100644 --- a/tests/TEST_TYPES +++ b/tests/TEST_TYPES @@ -58,7 +58,7 @@ Final value of function !6 Emp2-RI = !3 Final checksums !3 GLBOPT| Lowest reported potential energy !7 -xx,yy,zz !2 +POLAR| xx,yy,zz !3 POWELL| Final value of function !6 Gibbs energy correction !6 \ \ C !2