Skip to content

Commit

Permalink
xTB coulomb energies
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Mar 19, 2019
1 parent e611413 commit a3b9d68
Show file tree
Hide file tree
Showing 9 changed files with 1,019 additions and 118 deletions.
2 changes: 1 addition & 1 deletion src/mulliken.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -743,7 +743,7 @@ SUBROUTINE ao_charges_1(p_matrix, s_matrix, charges, iatom, para_env)
ELSEIF (iblock_col == iatom) THEN
DO i = 1, SIZE(p_block, 1)
DO j = 1, SIZE(p_block, 2)
charges(j) = charges(j)+mult*p_block(i, j)*s_block(i, j)
charges(j) = charges(j)+p_block(i, j)*s_block(i, j)
END DO
END DO
END IF
Expand Down
70 changes: 39 additions & 31 deletions src/qs_dftb3_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,17 +65,20 @@ MODULE qs_dftb3_methods
!> \param rho ...
!> \param mcharge ...
!> \param energy ...
!> \param xgamma ...
!> \param zeff ...
!> \param calculate_forces ...
!> \param just_energy ...
! **************************************************************************************************
SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, &
SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeff, &
calculate_forces, just_energy)

TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
TYPE(qs_rho_type), POINTER :: rho
REAL(dp), DIMENSION(:) :: mcharge
TYPE(qs_energy_type), POINTER :: energy
REAL(dp), DIMENSION(:) :: xgamma, zeff
LOGICAL, INTENT(in) :: calculate_forces, just_energy

CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb3_diagonal', &
Expand All @@ -88,7 +91,7 @@ SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, &
INTEGER, DIMENSION(:), POINTER :: atom_of_kind, kind_of
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: found, use_virial
REAL(KIND=dp) :: dr, eb3, eloc, fi, gmij, ua, ui, uj, zeff
REAL(KIND=dp) :: dr, eb3, eloc, fi, gmij, ua, ui, uj
REAL(KIND=dp), DIMENSION(3) :: fij, rij
REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
Expand All @@ -104,7 +107,6 @@ SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, &
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: n_list
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(virial_type), POINTER :: virial
Expand All @@ -119,15 +121,14 @@ SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, &
eb3 = 0.0_dp
CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
DO ikind = 1, SIZE(local_particles%n_el)
CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
CALL get_dftb_atom_param(dftb_kind, dudq=ua, zeff=zeff)
ua = xgamma(ikind)
DO ia = 1, local_particles%n_el(ikind)
iatom = local_particles%list(ikind)%array(ia)
eloc = -1.0_dp/6.0_dp*ua*mcharge(iatom)**3
eb3 = eb3+eloc
IF (atprop%energy) THEN
! we have to add the part not covered by 0.5*Tr(FP)
eloc = -0.5_dp*eloc-0.25_dp*ua*zeff*mcharge(iatom)**2
eloc = -0.5_dp*eloc-0.25_dp*ua*zeff(ikind)*mcharge(iatom)**2
atprop%atecoul(iatom) = atprop%atecoul(iatom)+eloc
END IF
END DO
Expand Down Expand Up @@ -220,20 +221,24 @@ SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, &
ic = cell_to_index(cellind(1), cellind(2), cellind(3))
CPASSERT(ic > 0)

ikind = kind_of(iatom)
atom_i = atom_of_kind(iatom)
CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
CALL get_dftb_atom_param(dftb_kind, dudq=ui)
jkind = kind_of(jatom)
atom_j = atom_of_kind(jatom)
CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind)
CALL get_dftb_atom_param(dftb_kind, dudq=uj)
!
gmij = -0.5_dp*(ui*mcharge(iatom)**2+uj*mcharge(jatom)**2)
!
NULLIFY (pblock)
CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
row=irow, col=icol, block=pblock, found=found)

ikind = kind_of(iatom)
atom_i = atom_of_kind(iatom)
ui = xgamma(ikind)
jkind = kind_of(jatom)
atom_j = atom_of_kind(jatom)
uj = xgamma(jkind)
!
gmij = -0.5_dp*(ui*mcharge(iatom)**2+uj*mcharge(jatom)**2)
!
NULLIFY (pblock)
CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
row=irow, col=icol, block=pblock, found=found)
CPASSERT(found)
DO i = 1, 3
NULLIFY (dsblock)
CALL dbcsr_get_block_p(matrix=matrix_s(1+i, ic)%matrix, &
row=irow, col=icol, block=dsblock, found=found)
CPASSERT(found)
DO i = 1, 3
NULLIFY (dsblock)
Expand Down Expand Up @@ -320,17 +325,20 @@ SUBROUTINE build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, &
ic = cell_to_index(cellind(1), cellind(2), cellind(3))
CPASSERT(ic > 0)

ikind = kind_of(iatom)
CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
CALL get_dftb_atom_param(dftb_kind, dudq=ui)
jkind = kind_of(jatom)
CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind)
CALL get_dftb_atom_param(dftb_kind, dudq=uj)
gmij = -0.5_dp*(ui*mcharge(iatom)**2+uj*mcharge(jatom)**2)
!
NULLIFY (sblock)
CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
row=irow, col=icol, block=sblock, found=found)
ikind = kind_of(iatom)
ui = xgamma(ikind)
jkind = kind_of(jatom)
uj = xgamma(jkind)
gmij = -0.5_dp*(ui*mcharge(iatom)**2+uj*mcharge(jatom)**2)
!
NULLIFY (sblock)
CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
row=irow, col=icol, block=sblock, found=found)
CPASSERT(found)
DO is = 1, SIZE(ks_matrix, 1)
NULLIFY (ksblock)
CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
row=irow, col=icol, block=ksblock, found=found)
CPASSERT(found)
DO is = 1, SIZE(ks_matrix, 1)
NULLIFY (ksblock)
Expand Down
26 changes: 22 additions & 4 deletions src/qs_dftb_coulomb.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,7 @@ SUBROUTINE build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb_coulomb', &
routineP = moduleN//':'//routineN

<<<<<<< HEAD
INTEGER :: atom_i, atom_j, blk, ewald_type, handle, i, ia, iatom, ic, icol, ikind, img, &
irow, is, jatom, jkind, natom, natorb_a, natorb_b, nimg, nmat
INTEGER, DIMENSION(3) :: cellind, periodic
Expand All @@ -114,6 +115,18 @@ SUBROUTINE build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
REAL(KIND=dp) :: alpha, ddr, deth, dgam, dr, drm, drp, &
fi, ga, gb, gmat, gmij, hb_para, zeff
REAL(KIND=dp), DIMENSION(0:3) :: eta_a, eta_b
=======
INTEGER :: atom_i, atom_j, blk, ewald_type, handle, &
i, ia, iatom, ic, icol, ikind, img, &
irow, is, jatom, jkind, natom, nimg, &
nkind, nmat
INTEGER, DIMENSION(3) :: cellind, periodic
INTEGER, DIMENSION(:), POINTER :: atom_of_kind, kind_of
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: do_ewald, found, use_virial
REAL(KIND=dp) :: alpha, deth, dr, fi, gmij, zeff
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma, zeffk
>>>>>>> xTB coulomb energies
REAL(KIND=dp), DIMENSION(3) :: fij, rij
REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, gmcharge, ksblock, pblock, &
sblock
Expand Down Expand Up @@ -302,12 +315,10 @@ SUBROUTINE build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
CALL get_qs_env(qs_env=qs_env, &
local_particles=local_particles)
DO ikind = 1, SIZE(local_particles%n_el)
CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
DO ia = 1, local_particles%n_el(ikind)
iatom = local_particles%list(ikind)%array(ia)
atprop%atecoul(iatom) = atprop%atecoul(iatom)+ &
0.5_dp*zeff*gmcharge(iatom, 1)
0.5_dp*mcharge(iatom)*gmcharge(iatom, 1)
END DO
END DO
END IF
Expand Down Expand Up @@ -463,9 +474,16 @@ SUBROUTINE build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
END IF

IF (dft_control%qs_control%dftb_control%dftb3_diagonal) THEN
CALL get_qs_env(qs_env, nkind=nkind)
ALLOCATE (zeffk(nkind), xgamma(nkind))
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
CALL get_dftb_atom_param(dftb_kind, dudq=xgamma(ikind), zeff=zeffk(ikind))
END DO
! Diagonal 3rd order correction (DFTB3)
CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, &
CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
calculate_forces, just_energy)
DEALLOCATE (zeffk, xgamma)
END IF

! QMMM
Expand Down
48 changes: 17 additions & 31 deletions src/qs_initial_guess.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,8 @@ MODULE qs_initial_guess
USE qs_wf_history_methods, ONLY: wfi_update
USE scf_control_types, ONLY: scf_control_type
USE util, ONLY: sort
USE xtb_types, ONLY: get_xtb_atom_param,&
xtb_atom_type
#include "./base/base_uses.f90"

IMPLICIT NONE
Expand Down Expand Up @@ -164,6 +166,7 @@ SUBROUTINE calculate_first_density_matrix(scf_env, qs_env)
TYPE(qs_rho_type), POINTER :: rho
TYPE(scf_control_type), POINTER :: scf_control
TYPE(section_vals_type), POINTER :: dft_section, input, subsys_section
TYPE(xtb_atom_type), POINTER :: xtb_kind

logger => cp_get_default_logger()
NULLIFY (atomic_kind, qs_kind, mo_coeff, sv, orb_basis_set, atomic_kind_set, &
Expand Down Expand Up @@ -1167,20 +1170,23 @@ SUBROUTINE calculate_mopac_dm(pmat, matrix_s, has_unit_metric, &
INTEGER :: atom_a, group, handle, iatom, ikind, &
iset, isgf, isgfa, ishell, ispin, la, &
maxl, maxll, nao, natom, ncount, nset, &
nsgf, z
nsgf, z, na
INTEGER, DIMENSION(5) :: occupation
INTEGER, DIMENSION(25) :: naox, laox
INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf
INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nshell
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, l, last_sgfa
LOGICAL :: has_pot
REAL(KIND=dp) :: maxocc, my_sum, nelec, paa, rscale, &
trps1, trps2, yy, zeff
trps1, trps2, yy, zeff, occ
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: econf, pdiag, sdiag
REAL(KIND=dp), DIMENSION(0:3) :: edftb
TYPE(all_potential_type), POINTER :: all_potential
TYPE(dbcsr_type), POINTER :: matrix_p
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(sgp_potential_type), POINTER :: sgp_potential
TYPE(xtb_atom_type), POINTER :: xtb_kind

CALL timeset(routineN, handle)

Expand Down Expand Up @@ -1238,6 +1244,9 @@ SUBROUTINE calculate_mopac_dm(pmat, matrix_s, has_unit_metric, &
lmax=maxll, occupation=edftb)
maxll = MIN(maxll, maxl)
econf(0:maxl) = edftb(0:maxl)
ELSEIF (dft_control%qs_control%xtb) THEN
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
CALL get_xtb_atom_param(xtb_kind, natorb=nsgf, nao=naox, lao=laox, occupation=occupation)
ELSEIF (has_pot) THEN
CALL get_atomic_kind(atomic_kind_set(ikind), z=z)
CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf, elec_conf=elec_conf, zeff=zeff)
Expand Down Expand Up @@ -1284,35 +1293,12 @@ SUBROUTINE calculate_mopac_dm(pmat, matrix_s, has_unit_metric, &
DO iatom = 1, natom
atom_a = atom_list(iatom)
isgfa = first_sgf(atom_a)
!qxtb
! DO la = 0, maxll
! SELECT CASE (la)
! CASE (0)
! pdiag(isgfa) = econf(0)
! CASE (1)
! pdiag(isgfa+1) = econf(1)/3._dp
! pdiag(isgfa+2) = econf(1)/3._dp
! pdiag(isgfa+3) = econf(1)/3._dp
! CASE (2)
! pdiag(isgfa+4) = econf(2)/5._dp
! pdiag(isgfa+5) = econf(2)/5._dp
! pdiag(isgfa+6) = econf(2)/5._dp
! pdiag(isgfa+7) = econf(2)/5._dp
! pdiag(isgfa+8) = econf(2)/5._dp
! CASE (3)
! pdiag(isgfa+9) = econf(3)/7._dp
! pdiag(isgfa+10) = econf(3)/7._dp
! pdiag(isgfa+11) = econf(3)/7._dp
! pdiag(isgfa+12) = econf(3)/7._dp
! pdiag(isgfa+13) = econf(3)/7._dp
! pdiag(isgfa+14) = econf(3)/7._dp
! pdiag(isgfa+15) = econf(3)/7._dp
! CASE DEFAULT
! CPABORT("")
! END SELECT
! END DO
pdiag(:) = 1._dp
!qxtb
DO isgf = 1, nsgf
na = naox(isgf)
la = laox(isgf)
occ = REAL(occupation(na), dp)/REAL(2*la+1,dp)
pdiag(isgfa+isgf-1) = occ
END DO
END DO
ELSEIF (dft_control%qs_control%semi_empirical) THEN
yy = REAL(dft_control%charge, KIND=dp)/REAL(nao, KIND=dp)
Expand Down
11 changes: 10 additions & 1 deletion src/qs_scf_output.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -311,7 +311,9 @@ SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelect

IF (output_unit > 0) THEN
CALL qs_rho_get(rho, tot_rho_r=tot_rho_r)
IF (.NOT. (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb)) THEN
IF (.NOT. (dft_control%qs_control%semi_empirical .OR. &
dft_control%qs_control%xtb .OR. &
dft_control%qs_control%dftb)) THEN
WRITE (UNIT=output_unit, FMT="(/,(T3,A,T41,2F20.10))") &
"Electronic density on regular grids: ", &
accurate_sum(tot_rho_r), &
Expand Down Expand Up @@ -365,6 +367,13 @@ SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelect
IF (energy%dftb3 /= 0.0_dp) &
WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
"DFTB3 3rd order energy: ", energy%dftb3
ELSEIF (dft_control%qs_control%xtb) THEN
WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
"Core Hamiltonian energy: ", energy%core, &
"Repulsive potential energy: ", energy%repulsive, &
"Electronic energy: ", energy%hartree, &
"DFTB3 3rd order energy: ", energy%dftb3, &
"Dispersion energy: ", energy%dispersion
ELSE
IF (dft_control%do_admm) THEN
exc_energy = energy%exc+energy%exc_aux_fit
Expand Down

0 comments on commit a3b9d68

Please sign in to comment.