Skip to content
Permalink
Browse files

Local energy code, local stress not yet working (#543)

* Local energy for GPW method

* Local stress (1. step)

* Towards local stress

* Towards local energy and stress fields

* Add stop for not yet implemented local stress

* Fix regtests
  • Loading branch information...
juerghutter committed Sep 6, 2019
1 parent 9891594 commit 690020138288c1ce1ce0a85e3532497b8c382769
@@ -84,7 +84,7 @@ MODULE bibliography
Stoychev2016, Futera2017, Bailey2006, Papior2017, Lehtola2018, &
Brieuc2016, Barca2018, Scheiber2018, Huang2011, Heaton_Burgess2007, &
Schuett2018, Holmberg2018, Togo2018, Staub2019, Grimme2017, Kondov2007, &
Ren2011, Ren2013
Ren2011, Ren2013, Cohen2000, Rogers2002, Filippetti2000

CONTAINS

@@ -4048,6 +4048,49 @@ SUBROUTINE add_all_references()
"ER"), &
DOI="10.1103/PhysRevB.88.035120")

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", &
"TI Theory and applications of the stress density", &
"SO PHYSICAL REVIEW B", &
"DT Article", &
"PY 2000", &
"VL 61", &
"BP 8433", &
"ER"), &
DOI="10.1103/PhysRevB.61.8433")

END SUBROUTINE add_all_references

END MODULE bibliography
@@ -104,6 +104,7 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
sphi_b, zeta, zetb
@@ -139,6 +140,10 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
END IF
END IF

IF (use_virial) THEN
pv_loc = virial%pv_virial
END IF

maxder = ncoset(nder)

CALL get_qs_kind_set(qs_kind_set, &
@@ -400,6 +405,9 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
END DO
END IF
END IF
IF (calculate_forces .AND. use_virial) THEN
virial%pv_ppl = virial%pv_virial-pv_loc
END IF

CALL timestop(handle)

@@ -118,6 +118,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
REAL(KIND=dp), DIMENSION(1:10) :: aloc, bloc
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc
REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
set_radius_a, set_radius_b
@@ -178,6 +179,10 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
END IF
END DO

IF (use_virial) THEN
pv_loc = virial%pv_virial
END IF

nthread = 1
!$ nthread = omp_get_max_threads()

@@ -492,6 +497,10 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
END IF
END IF

IF (calculate_forces .AND. use_virial) THEN
virial%pv_ppl = virial%pv_virial-pv_loc
END IF

CALL timestop(handle)

END SUBROUTINE build_core_ppl
@@ -120,6 +120,7 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work
REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
REAL(KIND=dp), DIMENSION(3) :: fa, fb, rab, rac, rbc
REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc
REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, h_nl, p_block, rpgfa, &
sphi_a, vprj_ppnl, zeta
@@ -165,6 +166,10 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
END IF
END IF

IF (use_virial) THEN
pv_loc = virial%pv_virial
END IF

maxder = ncoset(nder)

CALL get_qs_kind_set(qs_kind_set, &
@@ -532,6 +537,10 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
END IF
END IF

IF (calculate_forces .AND. use_virial) THEN
virial%pv_ppnl = virial%pv_virial-pv_loc
END IF

END IF !ppnl_present

CALL timestop(handle)
@@ -157,6 +157,7 @@ MODULE force_env_methods
qs_rho_type
USE restraint, ONLY: restraint_control
USE string_utilities, ONLY: compress
USE virial_methods, ONLY: write_stress_components
USE virial_types, ONLY: cp_virial,&
sym_virial,&
virial_create,&
@@ -381,6 +382,9 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%NDIGITS", &
i_val=ndigits)
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, ndigits, virial%pv_numer)
IF ((.NOT. virial%pv_numer) .AND. (force_env%in_use == use_qs_force)) THEN
CALL write_stress_components(virial, output_unit)
END IF
END IF
CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
"PRINT%STRESS_TENSOR")
@@ -1535,6 +1535,48 @@ SUBROUTINE create_print_dft_section(section)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

CALL cp_print_key_section_create(print_key, __LOCATION__, "LOCAL_ENERGY_CUBE", &
description="Controls the printing of cube files with the local "// &
" energy. It is valid only for QS with GPW/GAPW formalism ."// &
" Meta and hybrid functionals are not possible.", &
print_level=debug_print_level, filename="")
CALL keyword_create(keyword, __LOCATION__, name="stride", &
description="The stride (X,Y,Z) used to write the cube file "// &
"(larger values result in smaller cube files). You can provide 3 numbers (for X,Y,Z) or"// &
" 1 number valid for all components.", &
usage="STRIDE 2 2 2", n_var=-1, default_i_vals=(/2, 2, 2/), type_of_var=integer_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="APPEND", &
description="append the cube files when they already exist", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)

CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

CALL cp_print_key_section_create(print_key, __LOCATION__, "LOCAL_STRESS_CUBE", &
description="Controls the printing of cube files with the local "// &
" stress. It is valid only for QS with GPW/GAPW formalism ."// &
" Meta and hybrid functionals are not possible.", &
print_level=debug_print_level, filename="")
CALL keyword_create(keyword, __LOCATION__, name="stride", &
description="The stride (X,Y,Z) used to write the cube file "// &
"(larger values result in smaller cube files). You can provide 3 numbers (for X,Y,Z) or"// &
" 1 number valid for all components.", &
usage="STRIDE 2 2 2", n_var=-1, default_i_vals=(/2, 2, 2/), type_of_var=integer_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="APPEND", &
description="append the cube files when they already exist", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)

CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

CALL create_pdos_section(print_key)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
@@ -16,6 +16,7 @@
!> JGH (26-Feb-2003) : OpenMP enabled
!> JGH (17-Nov-2007) : Removed mass arrays
!> JGH (01-Dec-2007) : Removed and renamed routines
!> JGH (04-Jul-2019) : added pw_multiply routine
!> 03.2008 [tlaino] : Splitting pw_types into pw_types and pw_methods
!> \author apsi
! **************************************************************************************************
@@ -61,7 +62,7 @@ MODULE pw_methods

PUBLIC :: pw_zero, pw_structure_factor, pw_smoothing
PUBLIC :: pw_copy, pw_axpy, pw_transfer, pw_scale
PUBLIC :: pw_derive, pw_dr2, pw_write
PUBLIC :: pw_derive, pw_dr2, pw_write, pw_multiply
PUBLIC :: pw_integral_ab, pw_integral_a2b
PUBLIC :: pw_dr2_gg, pw_integrate_function
PUBLIC :: pw_set
@@ -764,7 +765,7 @@ SUBROUTINE pw_axpy(pw1, pw2, alpha)
IF (pw1%in_use == COMPLEXDATA1D .AND. &
pw2%in_use == COMPLEXDATA1D .AND. &
pw1%in_space == RECIPROCALSPACE .AND. &
pw1%in_space == RECIPROCALSPACE) THEN
pw2%in_space == RECIPROCALSPACE) THEN

ng1 = SIZE(pw1%cc)
ng2 = SIZE(pw2%cc)
@@ -886,6 +887,89 @@ SUBROUTINE pw_axpy(pw1, pw2, alpha)

END SUBROUTINE pw_axpy

! **************************************************************************************************
!> \brief pw_out = pw_out + alpha * pw1 * pw2
!> alpha defaults to 1
!> \param pw_out ...
!> \param pw1 ...
!> \param pw2 ...
!> \param alpha ...
!> \author JGH
! **************************************************************************************************
SUBROUTINE pw_multiply(pw_out, pw1, pw2, alpha)

TYPE(pw_type), INTENT(INOUT) :: pw_out
TYPE(pw_type), INTENT(IN) :: pw1, pw2
REAL(KIND=dp), INTENT(in), OPTIONAL :: alpha

CHARACTER(len=*), PARAMETER :: routineN = 'pw_multiply', routineP = moduleN//':'//routineN

INTEGER :: handle
REAL(KIND=dp) :: flop, my_alpha

CALL timeset(routineN, handle)

CPASSERT(pw_out%ref_count > 0)
CPASSERT(pw1%ref_count > 0)
CPASSERT(pw2%ref_count > 0)

my_alpha = 1.0_dp
IF (PRESENT(alpha)) my_alpha = alpha

IF (pw_out%pw_grid%id_nr == pw2%pw_grid%id_nr .AND. &
pw_out%pw_grid%id_nr == pw2%pw_grid%id_nr) THEN

IF (pw1%in_use == REALDATA1D .AND. pw2%in_use == REALDATA1D .AND. &
pw_out%in_use == REALDATA1D) THEN
IF (my_alpha == 1.0_dp) THEN
pw_out%cr = pw_out%cr+pw1%cr+pw2%cr
flop = REAL(2*SIZE(pw2%cr), KIND=dp)
ELSE
pw_out%cr = pw_out%cr+my_alpha*pw1%cr*pw2%cr
flop = REAL(3*SIZE(pw2%cr), KIND=dp)
END IF
ELSE IF (pw1%in_use == COMPLEXDATA1D .AND. &
pw2%in_use == COMPLEXDATA1D .AND. &
pw_out%in_use == COMPLEXDATA1D) THEN
IF (my_alpha == 1.0_dp) THEN
pw_out%cc = pw_out%cc+pw1%cc*pw2%cc
flop = REAL(3*SIZE(pw2%cc), KIND=dp)
ELSE
pw_out%cc = pw_out%cc+my_alpha*pw1%cc*pw2%cc
flop = REAL(6*SIZE(pw2%cc), KIND=dp)
END IF
ELSE IF (pw1%in_use == REALDATA3D .AND. pw2%in_use == REALDATA3D .AND. &
pw_out%in_use == REALDATA3D) THEN
IF (my_alpha == 1.0_dp) THEN
pw_out%cr3d = pw_out%cr3d+pw1%cr3d*pw2%cr3d
flop = REAL(2*SIZE(pw2%cr3d), KIND=dp)
ELSE
pw_out%cr3d = pw_out%cr3d+my_alpha*pw1%cr3d*pw2%cr3d
flop = REAL(3*SIZE(pw2%cr3d), KIND=dp)
END IF
ELSE IF (pw1%in_use == COMPLEXDATA3D .AND. &
pw2%in_use == COMPLEXDATA3D .AND. &
pw_out%in_use == COMPLEXDATA3D) THEN
IF (my_alpha == 1.0_dp) THEN
pw_out%cc3d = pw_out%cc3d+pw1%cc3d*pw2%cc3d
flop = REAL(3*SIZE(pw2%cc3d), KIND=dp)
ELSE
pw_out%cc3d = pw_out%cc3d+my_alpha*pw1%cc3d*pw2%cc3d
flop = REAL(6*SIZE(pw2%cc3d), KIND=dp)
END IF
ELSE
CPABORT("No suitable data field")
END IF

ELSE
CPABORT("Grids not compatible")
END IF

flop = flop*1.e-6_dp
CALL timestop(handle)

END SUBROUTINE pw_multiply

! **************************************************************************************************
!> \brief Gathers the pw vector from a 3d data field
!> \param pw ...
@@ -211,6 +211,7 @@ SUBROUTINE calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular
rab2, rootaab, zab
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, radius, zeff
REAL(KIND=dp), DIMENSION(3) :: deab, rab
REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(atprop_type), POINTER :: atprop
TYPE(neighbor_list_iterator_p_type), &
@@ -279,6 +280,7 @@ SUBROUTINE calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular
END DO

ecore_overlap = 0.0_dp
pv_loc = 0.0_dp

CALL neighbor_list_iterator_create(nl_iterator, sab_core)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
@@ -308,7 +310,7 @@ SUBROUTINE calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular
force(ikind)%core_overlap(:, atom_a) = force(ikind)%core_overlap(:, atom_a)+deab(:)
force(jkind)%core_overlap(:, atom_b) = force(jkind)%core_overlap(:, atom_b)-deab(:)
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, 1._dp, deab, rab)
CALL virial_pair_force(pv_loc, 1._dp, deab, rab)
END IF
END IF
END IF
@@ -319,6 +321,10 @@ SUBROUTINE calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular
IF (calculate_forces) THEN
DEALLOCATE (atom_of_kind)
END IF
IF (calculate_forces .AND. use_virial) THEN
virial%pv_virial = virial%pv_virial+pv_loc
virial%pv_hartree = virial%pv_hartree+pv_loc
END IF

CALL mp_sum(ecore_overlap, group)

@@ -879,7 +879,7 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rcpbc
REAL(KIND=dp), DIMENSION(3) :: fdij, fdik, ra, rab, rb, rb0, rbc, rc, &
rc0, rca, rij, rik, sab_max
REAL(KIND=dp), DIMENSION(3, 3) :: dvirial, pv_virial_thread
REAL(KIND=dp), DIMENSION(3, 3) :: dvirial, pv_loc, pv_virial_thread
REAL(KIND=dp), DIMENSION(:), POINTER :: atener
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: atstr
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
@@ -964,6 +964,9 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
IF (use_virial .AND. debugall) THEN
dvirial = virial%pv_virial
END IF
IF (use_virial) THEN
pv_loc = virial%pv_virial
END IF
ELSE IF ((dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
dispersion_env%pp_type == vdw_pairpot_dftd3bj) .AND. dispersion_env%doabc) THEN
ALLOCATE (atom_of_kind(natom), kind_of(natom))
@@ -1746,6 +1749,10 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
DEALLOCATE (atom2mol)
END IF

IF (calculate_forces .AND. use_virial) THEN
virial%pv_vdw = virial%pv_virial-pv_loc
END IF

IF (ASSOCIATED(dispersion_env%dftd_section)) THEN
CALL cp_print_key_finished_output(unit_nr, logger, dispersion_env%dftd_section, "PRINT_DFTD")
END IF

0 comments on commit 6900201

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