Skip to content

Commit

Permalink
Towards periodic efield for TB
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Mar 19, 2019
1 parent 6fa337f commit 3214a37
Show file tree
Hide file tree
Showing 6 changed files with 189 additions and 116 deletions.
160 changes: 109 additions & 51 deletions src/efield_tb_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,9 +266,9 @@ 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, handle, ia, iatom, ic, &
icol, idir, ikind, irow, is, ispin, &
jatom, jkind, natom, nimg, nspin
INTEGER :: atom_a, atom_b, blk, handle, ia, iatom, &
ic, icol, idir, ikind, irow, is, &
ispin, jatom, jkind, natom, nimg, nspin
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
INTEGER, DIMENSION(3) :: cellind
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
Expand All @@ -281,6 +281,7 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_iterator_type) :: iter
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
Expand Down Expand Up @@ -323,9 +324,8 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo

zi(:) = CMPLX(1._dp, 0._dp, dp)
DO ia = 1, natom
charge = -mcharge(ia)
charge = mcharge(ia)
ria = particle_set(ia)%r
ria = pbc(ria, cell)
DO idir = 1, 3
kvec(:) = twopi*cell%h_inv(idir, :)
dd = SUM(kvec(:)*ria(:))
Expand Down Expand Up @@ -368,63 +368,121 @@ SUBROUTINE efield_tb_berry(qs_env, ks_matrix, rho, mcharge, energy, calculate_fo
END IF
END IF

CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
iatom=iatom, jatom=jatom, r=rab, cell=cellind)
IF (nimg == 1) THEN
! no k-points; all matrices have been transformed to periodic bsf
CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)

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

icol = MAX(iatom, jatom)
irow = MIN(iatom, jatom)
DO is = 1, nspin
NULLIFY (ks_block)
CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
row=irow, col=icol, block=ks_block, found=found)
CPASSERT(found)
ks_block = ks_block+0.5_dp*fdir*s_block
END DO
IF (calculate_forces) THEN
ikind = kind_of(irow)
jkind = kind_of(icol)
atom_a = atom_of_kind(irow)
atom_b = atom_of_kind(icol)
fij = 0.0_dp
DO ispin = 1, nspin
CALL dbcsr_get_block_p(matrix=matrix_p(ispin, 1)%matrix, &
row=irow, col=icol, BLOCK=p_block, found=found)
CPASSERT(found)
DO idir = 1, 3
CALL dbcsr_get_block_p(matrix=matrix_s(idir+1, 1)%matrix, &
row=irow, col=icol, BLOCK=ds_block, found=found)
CPASSERT(found)
fij(idir) = fij(idir)+SUM(p_block*ds_block)
END DO
END DO
fdir = SUM(ria(1:3)*fieldpol(1:3))
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)
fdir = SUM(rib(1:3)*fieldpol(1:3))
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)
ELSE
CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
iatom=iatom, jatom=jatom, r=rab, cell=cellind)

icol = MAX(iatom, jatom)
irow = MIN(iatom, jatom)

IF (nimg > 1) THEN
ic = cell_to_index(cellind(1), cellind(2), cellind(3))
CPASSERT(ic > 0)
ELSE
ic = 1
END IF

ria = particle_set(irow)%r
ria = pbc(ria, cell)
rib = particle_set(icol)%r
rib = pbc(rib, cell)
fdir = 0.5_dp*SUM(fieldpol(1:3)*(ria(1:3)+rib(1:3)))
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

NULLIFY (s_block)
CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
row=irow, col=icol, block=s_block, found=found)
CPASSERT(found)
DO is = 1, nspin
NULLIFY (ks_block)
CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
row=irow, col=icol, block=ks_block, found=found)
NULLIFY (s_block)
CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
row=irow, col=icol, block=s_block, found=found)
CPASSERT(found)
ks_block = ks_block+fdir*s_block
END DO
IF (calculate_forces) THEN
atom_a = atom_of_kind(iatom)
atom_b = atom_of_kind(jatom)
fij = 0.0_dp
DO ispin = 1, nspin
CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
row=irow, col=icol, BLOCK=p_block, found=found)
DO is = 1, nspin
NULLIFY (ks_block)
CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
row=irow, col=icol, block=ks_block, found=found)
CPASSERT(found)
DO idir = 1, 3
CALL dbcsr_get_block_p(matrix=matrix_s(idir+1, ic)%matrix, &
row=irow, col=icol, BLOCK=ds_block, found=found)
ks_block = ks_block+0.5_dp*fdir*s_block
END DO
IF (calculate_forces) THEN
atom_a = atom_of_kind(iatom)
atom_b = atom_of_kind(jatom)
fij = 0.0_dp
DO ispin = 1, nspin
CALL dbcsr_get_block_p(matrix=matrix_p(ispin, ic)%matrix, &
row=irow, col=icol, BLOCK=p_block, found=found)
CPASSERT(found)
fij(idir) = fij(idir)+SUM(p_block*ds_block)
DO idir = 1, 3
CALL dbcsr_get_block_p(matrix=matrix_s(idir+1, ic)%matrix, &
row=irow, col=icol, BLOCK=ds_block, found=found)
CPASSERT(found)
fij(idir) = fij(idir)+SUM(p_block*ds_block)
END DO
END DO
END DO
fdir = SUM(ria(1:3)*fieldpol(1:3))
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)
fdir = SUM(rib(1:3)*fieldpol(1:3))
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
IF (irow == iatom) fij = -fij
force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a)+0.5_dp*fdir*fij(1:3)
force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b)-0.5_dp*fdir*fij(1:3)
fdir = SUM(rib(1:3)*fieldpol(1:3))
force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a)+0.5_dp*fdir*fij(1:3)
force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b)-0.5_dp*fdir*fij(1:3)
END IF

END DO
CALL neighbor_list_iterator_release(nl_iterator)
END DO
CALL neighbor_list_iterator_release(nl_iterator)
END IF

IF (calculate_forces) THEN
DEALLOCATE (atom_of_kind, kind_of)
Expand Down
20 changes: 11 additions & 9 deletions src/qs_environment.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -916,15 +916,17 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
dft_control%smear = scf_control%smear%do_smear

! Periodic efield needs equal occupation and orbital gradients
IF (dft_control%apply_period_efield) THEN
CALL get_qs_env(qs_env=qs_env, requires_mo_derivs=orb_gradient)
IF (.NOT. orb_gradient) THEN
CALL cp_abort(__LOCATION__, "Periodic Efield needs orbital gradient and direct optimization."// &
" Use the OT optimization method.")
END IF
IF (dft_control%smear) THEN
CALL cp_abort(__LOCATION__, "Periodic Efield needs equal occupation numbers."// &
" Smearing option is not possible.")
IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)) THEN
IF (dft_control%apply_period_efield) THEN
CALL get_qs_env(qs_env=qs_env, requires_mo_derivs=orb_gradient)
IF (.NOT. orb_gradient) THEN
CALL cp_abort(__LOCATION__, "Periodic Efield needs orbital gradient and direct optimization."// &
" Use the OT optimization method.")
END IF
IF (dft_control%smear) THEN
CALL cp_abort(__LOCATION__, "Periodic Efield needs equal occupation numbers."// &
" Smearing option is not possible.")
END IF
END IF
END IF

Expand Down
2 changes: 1 addition & 1 deletion tests/xTB/regtest-1/AdeThyvdW.inp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
&GLOBAL
PROJECT ATvdW
RUN_TYPE DEBUG
PRINT_LEVEL DEBUG
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
&DFT
Expand Down
9 changes: 7 additions & 2 deletions tests/xTB/regtest-1/h2o_dimer.inp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
&GLOBAL
PROJECT WaterDimer
RUN_TYPE DEBUG
PRINT_LEVEL DEBUG
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
&DFT
Expand Down Expand Up @@ -40,5 +40,10 @@
&END SUBSYS
&END FORCE_EVAL
&DEBUG
DEBUG_FORCES .TRUE.
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
DX 0.0001
EPS_NO_ERROR_CHECK 0.000001
STOP_ON_MISMATCH T
&END DEBUG

28 changes: 16 additions & 12 deletions tests/xTB/regtest-2/H2O-field.inp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,14 +10,18 @@
&QS
METHOD xTB
&xTB
DO_EWALD F
DO_EWALD T
&END XTB
&END QS
&KPOINTS
SCHEME GAMMA
&END
&SCF
EPS_SCF 1.e-7
MAX_SCF 30
&OT
MINIMIZER DIIS
EPS_SCF 1.e-8
MAX_SCF 90
&OT OFF
MINIMIZER CG
LINESEARCH 3PNT
PRECONDITIONER FULL_S_INVERSE
&END OT
&END SCF
Expand All @@ -28,18 +32,18 @@
GMAX 75
&END EWALD
&END POISSON
#&PERIODIC_EFIELD
# POLARISATION -1 -1 -1
# INTENSITY 0.001
#&END
&EFIELD
POLARISATION 0 0.5 1
&PERIODIC_EFIELD
POLARISATION 0 0 1
INTENSITY 0.05
&END
#&EFIELD
# POLARISATION 0 0.5 1
# INTENSITY 0.05
#&END
&END DFT
&SUBSYS
&CELL
ABC 30.0 30.0 30.0
ABC 6.0 6.0 6.0
&END CELL
&COORD
O 0.000000 0.000000 -0.065587 H2O
Expand Down

0 comments on commit 3214a37

Please sign in to comment.