Skip to content
Permalink
Browse files

Debug using finite field differences for dipole and polarizability. … (

#549)

* Local energy for GPW method

* Local stress (1. step)

* Towards local stress

* Add stop for not yet implemented local stress

* Fix regtests

* Debug using finite field differences for dipole and polarizability.
Enable IR intensities for semi-empirical methods.

* More regtests for DEBUG option

* Bug fix on call to get_results
  • Loading branch information...
juerghutter committed Sep 11, 2019
1 parent c5227fa commit 7863fcd744fd48424bb8d1d15d365a3009dec56c
@@ -4110,6 +4110,37 @@ SUBROUTINE add_all_references()
"ER"), &
DOI="10.1103/PhysRevB.65.224117")

CALL add_reference(key=Cohen2000, ISI_record=s2a( &
"AU Cohen, Morrel H.", &
" Frydel, Derek", &
" Burke, Kieron", &
" Engel, Eberhard", &
"AF Cohen, Morrel H.", &
" Frydel, Derek", &
" Burke, Kieron", &
" Engel, Eberhard", &
"TI Total energy density as an interpretative tool", &
"SO The Journal of Chemical Physics", &
"PY 2000", &
"VL 113", &
"BP 2990", &
"DI 10.1063/1.1286805", &
"ER"), &
DOI="10.1063/1.1286805")

CALL add_reference(key=Rogers2002, ISI_record=s2a( &
"AU Rogers, Christopher L.", &
" Rappe, Andrew M.", &
"TI Geometric formulation of quantum stress fields", &
"SO PHYSICAL REVIEW B", &
"DT Article", &
"PY 2002", &
"VL 65", &
"AR 224117", &
"DI 10.1103/PhysRevB.65.224117", &
"ER"), &
DOI="10.1103/PhysRevB.65.224117")

CALL add_reference(key=Filippetti2000, ISI_record=s2a( &
"AU Filippetti, Alessio", &
" Fiorentini, Vincenzo", &
@@ -23,11 +23,15 @@
! **************************************************************************************************
MODULE cp2k_debug
USE cell_types, ONLY: cell_type
USE cp_control_types, ONLY: dft_control_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE cp_result_methods, ONLY: get_results,&
test_for_result
USE cp_result_types, ONLY: cp_result_type
USE cp_subsys_types, ONLY: cp_subsys_get,&
cp_subsys_type
USE force_env_methods, ONLY: force_env_calc_energy_force,&
@@ -78,19 +82,28 @@ SUBROUTINE cp2k_debug_energy_and_forces(force_env)
routineP = moduleN//':'//routineN

CHARACTER(LEN=3*default_string_length) :: message
CHARACTER(LEN=default_string_length) :: description
INTEGER :: i, ip, iw, j, k, np, stress_tensor
LOGICAL :: check_failed, debug_forces, &
LOGICAL :: check_failed, debug_dipole, &
debug_forces, debug_polar, &
debug_stress_tensor, skip, &
stop_on_mismatch
REAL(KIND=dp) :: difference, dx, eps_no_error_check, &
std_value, sum_of_differences
REAL(KIND=dp) :: amplitude, dd, de, derr, difference, dx, &
eps_no_error_check, std_value, &
sum_of_differences
REAL(KIND=dp), DIMENSION(2) :: numer_energy
REAL(KIND=dp), DIMENSION(3) :: err, my_maxerr
REAL(KIND=dp), DIMENSION(3) :: dipole_moment, dipole_numer, err, &
my_maxerr, poldir
REAL(KIND=dp), DIMENSION(3, 2) :: dipn
REAL(KIND=dp), DIMENSION(3, 3) :: polar_analytic, polar_numeric
REAL(KIND=dp), DIMENSION(9) :: pvals
REAL(KIND=dp), DIMENSION(:, :), POINTER :: analyt_forces, numer_forces
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(cp_result_type), POINTER :: results
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(dft_control_type), POINTER :: dft_control
TYPE(particle_type), DIMENSION(:), POINTER :: particles
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(section_vals_type), POINTER :: root_section, subsys_section
@@ -109,8 +122,14 @@ SUBROUTINE cp2k_debug_energy_and_forces(force_env)
l_val=debug_stress_tensor)
CALL section_vals_val_get(root_section, "DEBUG%DEBUG_FORCES", &
l_val=debug_forces)
CALL section_vals_val_get(root_section, "DEBUG%DEBUG_DIPOLE", &
l_val=debug_dipole)
CALL section_vals_val_get(root_section, "DEBUG%DEBUG_POLARIZABILITY", &
l_val=debug_polar)
CALL section_vals_val_get(root_section, "DEBUG%DX", &
r_val=dx)
CALL section_vals_val_get(root_section, "DEBUG%DE", &
r_val=de)
dx = ABS(dx)
CALL section_vals_val_get(root_section, "DEBUG%EPS_NO_ERROR_CHECK", &
r_val=eps_no_error_check)
@@ -332,9 +351,142 @@ SUBROUTINE cp2k_debug_energy_and_forces(force_env)
IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
END IF
CALL cp_print_key_finished_output(iw, logger, root_section, &
"DEBUG%PROGRAM_RUN_INFO")

IF (debug_dipole) THEN
IF (force_env%in_use == use_qs_force) THEN
CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
poldir = (/0.0_dp, 0.0_dp, 1.0_dp/)
amplitude = 0.0_dp
CALL set_efield(dft_control, amplitude, poldir)
CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
CALL get_qs_env(force_env%qs_env, results=results)
description = "[DIPOLE]"
IF (test_for_result(results, description=description)) THEN
CALL get_results(results, description=description, values=dipole_moment)
ELSE
CALL cp_warn(__LOCATION__, "Debug of dipole moments needs DFT/PRINT/MOMENTS section enabled")
CPABORT("DEBUG")
END IF
amplitude = de
DO k = 1, 3
poldir = 0.0_dp
poldir(k) = 1.0_dp
DO j = 1, 2
poldir = -1.0_dp*poldir
CALL set_efield(dft_control, amplitude, poldir)
CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
CALL force_env_get(force_env, potential_energy=numer_energy(j))
END DO
dipole_numer(k) = 0.5_dp*(numer_energy(1)-numer_energy(2))/de
END DO
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(/,(T2,A))") &
"DEBUG|========================= DIPOLE MOMENTS ================================", &
"DEBUG| Coordinate D(numerical) D(analytical) Difference Error [%]"
DO k = 1, 3
dd = dipole_moment(k)-dipole_numer(k)
IF (ABS(dipole_moment(k)) > eps_no_error_check) THEN
derr = 100._dp*dd/dipole_moment(k)
WRITE (UNIT=iw, FMT="(T13,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
ACHAR(119+k), dipole_numer(k), dipole_moment(k), dd, derr
ELSE
WRITE (UNIT=iw, FMT="(T13,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
ACHAR(119+k), dipole_numer(k), dipole_moment(k), dd
END IF
END DO
WRITE (UNIT=iw, FMT="((T2,A))") &
"DEBUG|========================================================================="
WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") ' DIPOLE : CheckSum =', SUM(dipole_moment)
END IF
ELSE
CALL cp_warn(__LOCATION__, "Debug of dipole moments only for Quickstep code available")
END IF
END IF

IF (debug_polar) THEN
IF (force_env%in_use == use_qs_force) THEN
CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
poldir = (/0.0_dp, 0.0_dp, 1.0_dp/)
amplitude = 0.0_dp
CALL set_efield(dft_control, amplitude, poldir)
CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
CALL get_qs_env(force_env%qs_env, results=results)
description = "[POLAR]"
IF (test_for_result(results, description=description)) THEN
CALL get_results(results, description=description, values=pvals)
polar_analytic(1:3, 1:3) = RESHAPE(pvals(1:9), (/3, 3/))
ELSE
CALL cp_warn(__LOCATION__, "Debug of polarizabilities needs PROPERTIES/LINRES/POLAR section enabled")
CPABORT("DEBUG")
END IF
description = "[DIPOLE]"
IF (.NOT. test_for_result(results, description=description)) THEN
CALL cp_warn(__LOCATION__, "Debug of polarizabilities need DFT/PRINT/MOMENTS section enabled")
CPABORT("DEBUG")
END IF
amplitude = de
DO k = 1, 3
poldir = 0.0_dp
poldir(k) = 1.0_dp
DO j = 1, 2
poldir = -1.0_dp*poldir
CALL set_efield(dft_control, amplitude, poldir)
CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., linres=.TRUE.)
CALL get_results(results, description=description, values=dipn(1:3, j))
END DO
polar_numeric(1:3, k) = 0.5_dp*(dipn(1:3, 2)-dipn(1:3, 1))/de
END DO
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(/,(T2,A))") &
"DEBUG|========================= POLARIZABILITY ================================", &
"DEBUG| Coordinates P(numerical) P(analytical) Difference Error [%]"
DO k = 1, 3
DO j = 1, 3
dd = polar_analytic(k, j)-polar_numeric(k, j)
IF (ABS(polar_analytic(k, j)) > eps_no_error_check) THEN
derr = 100._dp*dd/polar_analytic(k, j)
WRITE (UNIT=iw, FMT="(T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
ACHAR(119+k), ACHAR(119+j), polar_numeric(k, j), polar_analytic(k, j), dd, derr
ELSE
WRITE (UNIT=iw, FMT="(T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
ACHAR(119+k), ACHAR(119+j), polar_numeric(k, j), polar_analytic(k, j), dd
END IF
END DO
END DO
WRITE (UNIT=iw, FMT="((T2,A))") &
"DEBUG|========================================================================="
WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") ' POLAR : CheckSum =', SUM(polar_analytic)
END IF
ELSE
CALL cp_warn(__LOCATION__, "Debug of polarizabilities only for Quickstep code available")
END IF
END IF

CALL cp_print_key_finished_output(iw, logger, root_section, "DEBUG%PROGRAM_RUN_INFO")

END SUBROUTINE cp2k_debug_energy_and_forces

! **************************************************************************************************
!> \brief ...
!> \param dft_control ...
!> \param amplitude ...
!> \param poldir ...
! **************************************************************************************************
SUBROUTINE set_efield(dft_control, amplitude, poldir)
TYPE(dft_control_type), POINTER :: dft_control
REAL(KIND=dp), INTENT(IN) :: amplitude
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: poldir

IF (dft_control%apply_efield) THEN
dft_control%efield_fields(1)%efield%strength = amplitude
dft_control%efield_fields(1)%efield%polarisation(1:3) = poldir(1:3)
ELSEIF (dft_control%apply_period_efield) THEN
dft_control%period_efield%strength = amplitude
dft_control%period_efield%polarisation(1:3) = poldir(1:3)
ELSE
CPABORT("No EFIELD definition available")
END IF

END SUBROUTINE set_efield

END MODULE cp2k_debug
@@ -1685,9 +1685,9 @@ SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmo
DO l = 0, nmom
SELECT CASE (l)
CASE (0)
WRITE (unit_number, "(T3,A,T32,3F14.8)") "Reference Point [Bohr]", rcc
WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
WRITE (unit_number, "(T3,A)") "Charges"
WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T68,F14.8)") &
WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
"Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
CASE (1)
IF (periodic) THEN
@@ -1701,7 +1701,7 @@ SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmo
ENDIF
dd = SQRT(SUM(rmom(2:4, 3)**2))*debye
WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
WRITE (unit_number, "(T5,3(A,A,F14.8,1X),T60,A,T68,F14.8)") &
WRITE (unit_number, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
(TRIM(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
CASE (2)
WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
@@ -1743,7 +1743,7 @@ SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmo
IF (nmom >= 1) THEN
dd = SQRT(SUM(mmom(1:3)**2))
WRITE (unit_number, "(T3,A)") "Magnetic Dipole Moment (only orbital contrib.) [a.u.]"
WRITE (unit_number, "(T5,3(A,A,F14.8,1X),T60,A,T68,F14.8)") &
WRITE (unit_number, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
(TRIM(rlab(i+1)), "=", mmom(i), i=1, 3), "Total=", dd
END IF
END IF
@@ -28,6 +28,9 @@ MODULE qs_scf_post_se
cp_print_key_should_output,&
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE cp_result_methods, ONLY: cp_results_erase,&
put_results
USE cp_result_types, ONLY: cp_result_type
USE dbcsr_api, ONLY: dbcsr_get_block_p,&
dbcsr_p_type
USE input_section_types, ONLY: section_get_ival,&
@@ -225,6 +228,7 @@ SUBROUTINE qs_scf_post_moments(input, logger, qs_env)
REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_result_type), POINTER :: results
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: basis_set
@@ -234,6 +238,7 @@ SUBROUTINE qs_scf_post_moments(input, logger, qs_env)
TYPE(section_vals_type), POINTER :: print_key
TYPE(semi_empirical_type), POINTER :: se_kind

NULLIFY (results)
print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MOMENTS")
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
! Dipole Moments
@@ -257,6 +262,7 @@ SUBROUTINE qs_scf_post_moments(input, logger, qs_env)
natom=natom, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
results=results, &
dft_control=dft_control)

CALL qs_rho_get(rho, rho_ao=matrix_p)
@@ -369,6 +375,9 @@ SUBROUTINE qs_scf_post_moments(input, logger, qs_env)
END DO
DEALLOCATE (mom)
END DO
CALL cp_results_erase(results=results, description=description)
CALL put_results(results=results, description=description, &
values=dipole(1:3))
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(1X,A,T48,3F11.6)') "DIPOLE "//TRIM(dipole_type)//"(A.U.)|", dipole
WRITE (unit_nr, '(1X,A,T48,3F11.6)') "DIPOLE "//TRIM(dipole_type)//"(Debye)|", dipole*debye
@@ -330,12 +330,32 @@ SUBROUTINE create_debug_section(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="DEBUG_DIPOLE", &
description="Activates the debugging of the dipole moment", &
usage="DEBUG_DIPOLE {LOGICAL}", default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="DEBUG_POLARIZABILITY", &
description="Activates the debugging of the polarizability", &
usage="DEBUG_POLARIZABILITY {LOGICAL}", default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="DX", &
description="Increment for the calculation of the numerical derivatives", &
usage="DX {REAL}", default_r_val=0.001_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="DE", &
description="Increment for the calculation of the numerical electric field derivatives", &
usage="DE {REAL}", default_r_val=0.0001_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="EPS_NO_ERROR_CHECK", &
description="The mismatch between the numerical and the "// &
"analytical value is not checked for analytical "// &

0 comments on commit 7863fcd

Please sign in to comment.
You can’t perform that action at this time.