comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Integrate polar_env into qs_env; enable linres restart for polar; new regression test input added

Browse source 
svn-origin-rev: 18498
This commit is contained in:
Matthias Krack 2018-07-01 16:55:45 +00:00
parent af840dcaf6
commit 4a6d2ced8a
11 changed files with 590 additions and 328 deletions
Download patch file Download diff file Expand all files Collapse all files

11
src/input_cp2k_motion_print.F
View file

@ -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", &

9
src/motion/md_energies.F
View file

@ -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

31
src/qs_environment_types.F
View file

@ -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

25
src/qs_linres_methods.F
View file

@ -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)

39
src/qs_linres_module.F
View file

@ -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

22
src/qs_linres_op.F
View file

@ -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")

547
src/qs_linres_polar_utils.F
View file

@ -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

119
src/qs_linres_types.F
View file

@ -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

108
tests/QS/regtest-polar/H2O_md_polar.inp Normal file
View file

@ -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

5
tests/QS/regtest-polar/TEST_FILES
View file

@ -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

2
tests/TEST_TYPES
View file

@ -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