Skip to content

Commit

Permalink
xTB/DFTB stress tensor with efield (#1405)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Mar 9, 2021
1 parent 7743848 commit 41c7a35
Show file tree
Hide file tree
Showing 7 changed files with 328 additions and 19 deletions.
108 changes: 90 additions & 18 deletions src/efield_tb_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,8 +46,11 @@ MODULE efield_tb_methods
init_efield_matrices
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE sap_kind_types, ONLY: release_sap_int,&
sap_int_type
USE virial_methods, ONLY: virial_pair_force
USE virial_types, ONLY: virial_type
USE xtb_coulomb, ONLY: xtb_dsint_list
#include "./base/base_uses.f90"

IMPLICIT NONE
Expand Down Expand Up @@ -278,18 +281,20 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo

COMPLEX(KIND=dp) :: zdeta
COMPLEX(KIND=dp), DIMENSION(3) :: zi(3)
INTEGER :: atom_a, atom_b, blk, handle, ia, iatom, &
ic, icol, idir, ikind, irow, is, &
ispin, jatom, jkind, natom, nimg, nspin
INTEGER :: atom_a, atom_b, blk, handle, ia, iac, &
iatom, ic, icol, idir, ikind, irow, &
is, ispin, jatom, jkind, natom, nimg, &
nkind, nspin
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
INTEGER, DIMENSION(3) :: cellind
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: found, use_virial
REAL(KIND=dp) :: charge, dd, fdir
REAL(KIND=dp) :: charge, dd, dr, fdir, fi
REAL(KIND=dp), DIMENSION(3) :: fieldpol, fij, forcea, fpolvec, kvec, &
qi, rab, ria, rib
qi, rab, ria, rib, rij
REAL(KIND=dp), DIMENSION(3, 3) :: hmat
REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, ks_block, p_block, s_block
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_para_env_type), POINTER :: para_env
Expand All @@ -304,6 +309,7 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
TYPE(virial_type), POINTER :: virial

CALL timeset(routineN, handle)
Expand All @@ -318,10 +324,6 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
! calculate stress only if forces requested also
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
use_virial = use_virial .AND. calculate_forces
! disable stress calculation
IF (use_virial) THEN
CPABORT("Stress tensor for periodic E-field not implemented")
END IF

fieldpol = dft_control%period_efield%polarisation
fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
Expand Down Expand Up @@ -374,7 +376,8 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
ria = particle_set(ia)%r
ria = pbc(ria, cell)
forcea(1:3) = fieldpol(1:3)*charge
CALL virial_pair_force(virial%pv_virial, 1.0_dp, forcea, ria)
CALL virial_pair_force(virial%pv_virial, -0.5_dp, forcea, ria)
CALL virial_pair_force(virial%pv_virial, -0.5_dp, ria, forcea)
END IF
END DO
ELSE
Expand Down Expand Up @@ -432,9 +435,76 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fdir*fij(1:3)
END IF

ENDDO
CALL dbcsr_iterator_stop(iter)
!
! stress tensor for Gamma point needs to recalculate overlap integral derivatives
!
IF (use_virial) THEN
! derivative overlap integral (non collapsed)
NULLIFY (sap_int)
IF (dft_control%qs_control%dftb) THEN
CPABORT("DFTB stress tensor for periodic efield not implemented")
ELSEIF (dft_control%qs_control%xtb) THEN
CALL xtb_dsint_list(qs_env, sap_int)
ELSE
CPABORT("TB method unknown")
END IF
!
CALL get_qs_env(qs_env, nkind=nkind)
DO ikind = 1, nkind
DO jkind = 1, nkind
iac = ikind + nkind*(jkind - 1)
IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
DO ia = 1, sap_int(iac)%nalist
IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
iatom = sap_int(iac)%alist(ia)%aatom
DO ic = 1, sap_int(iac)%alist(ia)%nclist
jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
rij = sap_int(iac)%alist(ia)%clist(ic)%rac
dr = SQRT(SUM(rij(:)**2))
IF (dr > 1.e-6_dp) THEN
dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
icol = MAX(iatom, jatom)
irow = MIN(iatom, jatom)
IF (irow == iatom) rij = -rij
fdir = 0.0_dp
ria = particle_set(irow)%r
rib = particle_set(icol)%r
DO idir = 1, 3
kvec(:) = twopi*cell%h_inv(idir, :)
dd = SUM(kvec(:)*ria(:))
zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
dd = SUM(kvec(:)*rib(:))
zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
END DO
fi = 1.0_dp
IF (iatom == jatom) fi = 0.5_dp
DO ispin = 1, nspin
NULLIFY (p_block)
CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
row=irow, col=icol, block=p_block, found=found)
CPASSERT(found)
fij = 0.0_dp
DO idir = 1, 3
IF (irow == iatom) THEN
fij(idir) = SUM(p_block*dsint(:, :, idir))
ELSE
fij(idir) = SUM(TRANSPOSE(p_block)*dsint(:, :, idir))
END IF
END DO
IF (irow == iatom) fij = -fij
CALL virial_pair_force(virial%pv_virial, fi, fdir*fij(1:3), rij)
END DO
END IF
END DO
END DO
END DO
END DO
CALL release_sap_int(sap_int)
END IF
ELSE
CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
Expand Down Expand Up @@ -489,8 +559,15 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
IF (irow == iatom) fij = -fij
force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fdir*fij(1:3)
force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) + fdir*fij(1:3)
IF (use_virial) THEN
dr = SQRT(SUM(rab(:)**2))
IF (dr > 1.e-6_dp) THEN
fi = 1.0_dp
IF (iatom == jatom) fi = 0.5_dp
CALL virial_pair_force(virial%pv_virial, fi, -fdir*fij(1:3), rab)
END IF
END IF
END IF

END DO
CALL neighbor_list_iterator_release(nl_iterator)
END IF
Expand Down Expand Up @@ -578,7 +655,7 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
use_virial = use_virial .AND. calculate_forces
! disable stress calculation
IF (use_virial) THEN
CPABORT("Stress tensor for periodic E-field not implemented")
CPABORT("Stress tensor for periodic D-field not implemented")
END IF

dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
Expand Down Expand Up @@ -671,11 +748,6 @@ SUBROUTINE dfield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
iatom = atom_of_kind(ia)
ikind = kind_of(ia)
force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + di(:)*charge
IF (use_virial) THEN
ria = particle_set(ia)%r
ria = pbc(ria, cell)
CALL virial_pair_force(virial%pv_virial, 1.0_dp, di*charge, ria)
END IF
END DO
END IF
END IF
Expand Down
2 changes: 1 addition & 1 deletion src/xtb_coulomb.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ MODULE xtb_coulomb

CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_coulomb'

PUBLIC :: build_xtb_coulomb, gamma_rab_sr
PUBLIC :: build_xtb_coulomb, gamma_rab_sr, xtb_dsint_list

CONTAINS

Expand Down
4 changes: 4 additions & 0 deletions tests/xTB/regtest-5/TEST_FILES
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,8 @@ AdeThyvdW.inp 1 1.0E-12 -58.89501875
ice.inp 1 1.0E-10 -46.31849039679943
ice2.inp 1 1.0E-10 -46.30633640965694
Ru_geo.inp 1 1.0E-10 -7.01246011345805
ef_stress1.inp 0
ef_stress2.inp 0
ef_stress3.inp 0
ef_stress4.inp 0
#EOF
59 changes: 59 additions & 0 deletions tests/xTB/regtest-5/ef_stress1.inp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
&GLOBAL
PRINT_LEVEL low
RUN_TYPE DEBUG
&END GLOBAL

&DEBUG
DEBUG_FORCES yes
DEBUG_STRESS_TENSOR yes
EPS_NO_ERROR_CHECK 1.0E-4
DX 1.E-3
STOP_ON_MISMATCH yes
&END DEBUG

&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR analytical
&PRINT
&STRESS_TENSOR ON
COMPONENTS
&END STRESS_TENSOR
&END
&DFT
&KPOINTS
SCHEME GAMMA
&END
&QS
EPS_DEFAULT 1.0E-14
EXTRAPOLATION USE_PREV_P
METHOD xTB
&xTB
DO_EWALD T
&END xTB
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 200
&MIXING
METHOD DIRECT_P_MIXING
ALPHA 0.25
&END
&END SCF
&PERIODIC_EFIELD
POLARISATION 1 1 1
INTENSITY 0.005
&END
&END DFT
&SUBSYS
&CELL
ABC 4.0000 4.0000 4.0000
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

56 changes: 56 additions & 0 deletions tests/xTB/regtest-5/ef_stress2.inp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
&GLOBAL
PRINT_LEVEL low
RUN_TYPE DEBUG
&END GLOBAL

&DEBUG
DEBUG_FORCES yes
DEBUG_STRESS_TENSOR yes
EPS_NO_ERROR_CHECK 1.0E-4
DX 1.E-3
STOP_ON_MISMATCH yes
&END DEBUG

&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR analytical
&PRINT
&STRESS_TENSOR ON
COMPONENTS
&END STRESS_TENSOR
&END
&DFT
&QS
EPS_DEFAULT 1.0E-14
EXTRAPOLATION USE_PREV_P
METHOD xTB
&xTB
DO_EWALD T
&END xTB
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 100
&MIXING
METHOD DIRECT_P_MIXING
ALPHA 0.15
&END
&END SCF
&PERIODIC_EFIELD
POLARISATION 1 1 1
INTENSITY 0.005
&END
&END DFT
&SUBSYS
&CELL
ABC 4.0000 4.0000 4.0000
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

61 changes: 61 additions & 0 deletions tests/xTB/regtest-5/ef_stress3.inp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
&GLOBAL
PRINT_LEVEL low
RUN_TYPE DEBUG
&END GLOBAL

&DEBUG
DEBUG_FORCES yes
DEBUG_STRESS_TENSOR yes
EPS_NO_ERROR_CHECK 1.0E-4
DX 1.E-3
STOP_ON_MISMATCH yes
&END DEBUG

&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR analytical
&PRINT
&STRESS_TENSOR ON
COMPONENTS
&END STRESS_TENSOR
&END
&DFT
LSD
&KPOINTS
SCHEME GAMMA
&END
&QS
EPS_DEFAULT 1.0E-14
EXTRAPOLATION USE_PREV_P
METHOD xTB
&xTB
DO_EWALD T
CHECK_ATOMIC_CHARGES Off
&END xTB
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 400
&MIXING
METHOD DIRECT_P_MIXING
ALPHA 0.25
&END
&END SCF
&PERIODIC_EFIELD
POLARISATION 1 1 1
INTENSITY 0.005
&END
&END DFT
&SUBSYS
&CELL
ABC 4.0000 4.0000 4.0000
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

0 comments on commit 41c7a35

Please sign in to comment.