Skip to content

Commit

Permalink
charge mixing
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Mar 19, 2019
1 parent 6e8ba42 commit 9e14bf3
Show file tree
Hide file tree
Showing 5 changed files with 319 additions and 8 deletions.
269 changes: 269 additions & 0 deletions src/qs_charge_mixing.F
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,269 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2019 CP2K developers group !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
MODULE qs_charge_mixing

USE cp_para_types, ONLY: cp_para_env_type
USE kinds, ONLY: dp
USE mathlib, ONLY: get_pseudo_inverse_svd
USE message_passing, ONLY: mp_sum
USE qs_density_mixing_types, ONLY: broyden_mixing_new_nr,&
broyden_mixing_nr,&
gspace_mixing_nr,&
mixing_storage_type,&
multisecant_mixing_nr,&
pulay_mixing_nr
USE qs_environment_types, ONLY: qs_environment_type
#include "./base/base_uses.f90"

IMPLICIT NONE

PRIVATE

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

PUBLIC :: charge_mixing

CONTAINS

! **************************************************************************************************
!> \brief Driver for the charge mixing, calls the proper routine given the requested method
!> \param qs_env ...
!> \param mixing_method ...
!> \param mixing_store ...
!> \param charges ...
!> \param para_env ...
!> \param iter_count ...
!> \par History
!> \author JGH
! **************************************************************************************************
SUBROUTINE charge_mixing(qs_env, mixing_method, mixing_store, charges, para_env, iter_count)
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER, INTENT(IN) :: mixing_method
TYPE(mixing_storage_type), POINTER :: mixing_store
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
TYPE(cp_para_env_type), POINTER :: para_env
INTEGER, INTENT(IN) :: iter_count

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

INTEGER :: handle, ia, ii, imin, inow, nbuffer, ns, &
nvec
REAL(dp) :: alpha

CALL timeset(routineN, handle)

IF (mixing_method >= gspace_mixing_nr) THEN
CPASSERT(ASSOCIATED(mixing_store))
mixing_store%ncall = mixing_store%ncall+1
ns = SIZE(charges, 2)
ns = MIN(ns, mixing_store%max_shell)
alpha = mixing_store%alpha
nbuffer = mixing_store%nbuffer
inow = MOD(mixing_store%ncall-1, nbuffer)+1
imin = inow-1
IF (imin == 0) imin = nbuffer
IF (mixing_store%ncall > nbuffer) THEN
nvec = nbuffer
ELSE
nvec = mixing_store%ncall-1
END IF
IF (mixing_store%ncall > 1) THEN
! store in/out charge difference
DO ia = 1, mixing_store%nat_local
ii = mixing_store%atlist(ia)
mixing_store%dacharge(ia, 1:ns, imin) = mixing_store%acharge(ia, 1:ns, imin)-charges(ii, 1:ns)
END DO
END IF
IF ((iter_count == 1) .OR. (iter_count+1 <= mixing_store%nskip_mixing)) THEN
! skip mixing
mixing_store%iter_method = "NoMix"
ELSEIF (((iter_count+1-mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
mixing_store%iter_method = "Mixing"
ELSEIF (mixing_method == gspace_mixing_nr) THEN
CPABORT("Kerker method not available for Charge Mixing")
ELSEIF (mixing_method == pulay_mixing_nr) THEN
CPABORT("Pulay method not available for Charge Mixing")
ELSEIF (mixing_method == broyden_mixing_nr) THEN
CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
mixing_store%iter_method = "Broy."
ELSEIF (mixing_method == broyden_mixing_new_nr) THEN
CPABORT("Broyden_mixing_new method not available for Charge Mixing")
ELSEIF (mixing_method == multisecant_mixing_nr) THEN
CPABORT("Multisecant_mixing method not available for Charge Mixing")
END IF

! store new 'input' charges
DO ia = 1, mixing_store%nat_local
ii = mixing_store%atlist(ia)
mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
END DO

END IF

CALL timestop(handle)

END SUBROUTINE charge_mixing

! **************************************************************************************************
!> \brief Simple charge mixing
!> \param mixing_store ...
!> \param charges ...
!> \param alpha ...
!> \param imin ...
!> \param ns ...
!> \param para_env ...
!> \author JGH
! **************************************************************************************************
SUBROUTINE mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
TYPE(mixing_storage_type), POINTER :: mixing_store
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
REAL(KIND=dp), INTENT(IN) :: alpha
INTEGER, INTENT(IN) :: imin, ns
TYPE(cp_para_env_type), POINTER :: para_env

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

INTEGER :: ia, ii

charges = 0.0_dp

DO ia = 1, mixing_store%nat_local
ii = mixing_store%atlist(ia)
charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin)-mixing_store%acharge(ia, 1:ns, imin)
END DO

CALL mp_sum(charges, para_env%group)

END SUBROUTINE mix_charges_only

! **************************************************************************************************
!> \brief Broyden charge mixing
!> \param mixing_store ...
!> \param charges ...
!> \param inow ...
!> \param nvec ...
!> \param ns ...
!> \param para_env ...
!> \author JGH
! **************************************************************************************************
SUBROUTINE broyden_mixing(mixing_store, charges, inow, nvec, ns, para_env)
TYPE(mixing_storage_type), POINTER :: mixing_store
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
INTEGER, INTENT(IN) :: inow, nvec, ns
TYPE(cp_para_env_type), POINTER :: para_env

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

INTEGER :: i, ia, ii, imin, j, nbuffer, nv
REAL(KIND=dp) :: alpha, maxw, minw, omega0, rskip, wdf, &
wfac
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now

CPASSERT(nvec > 1)

nbuffer = mixing_store%nbuffer
alpha = mixing_store%alpha
imin = inow-1
IF (imin == 0) imin = nvec
nv = nvec-1

! charge vectors
q_now => mixing_store%acharge(:, :, inow)
q_last => mixing_store%acharge(:, :, imin)
dq_now => mixing_store%dacharge(:, :, inow)
dq_last => mixing_store%dacharge(:, :, imin)

IF (nvec == nbuffer) THEN
! reshuffel Broyden storage n->n-1
DO i = 1, nv-1
mixing_store%wbroy(i) = mixing_store%wbroy(i+1)
mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i+1)
mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i+1)
END DO
DO i = 1, nv-1
DO j = 1, nv-1
mixing_store%abroy(i, j) = mixing_store%abroy(i+1, j+1)
END DO
END DO
END IF

omega0 = 0.01_dp
minw = 1.0_dp
maxw = 100000.0_dp
wfac = 0.01_dp

mixing_store%wbroy(nv) = SUM(dq_now(:, :)**2)
CALL mp_sum(mixing_store%wbroy(nv), para_env%group)
mixing_store%wbroy(nv) = SQRT(mixing_store%wbroy(nv))
IF (mixing_store%wbroy(nv) > (wfac/maxw)) THEN
mixing_store%wbroy(nv) = wfac/mixing_store%wbroy(nv)
ELSE
mixing_store%wbroy(nv) = maxw
ENDIF
IF (mixing_store%wbroy(nv) < minw) mixing_store%wbroy(nv) = minw

! dfbroy
mixing_store%dfbroy(:, :, nv) = dq_now(:, :)-dq_last(:, :)
wdf = SUM(mixing_store%dfbroy(:, :, nv)**2)
CALL mp_sum(wdf, para_env%group)
wdf = 1.0_dp/SQRT(wdf)
mixing_store%dfbroy(:, :, nv) = wdf*mixing_store%dfbroy(:, :, nv)

! abroy matrix
DO i = 1, nv
wfac = SUM(mixing_store%dfbroy(:, :, i)*mixing_store%dfbroy(:, :, nv))
CALL mp_sum(wfac, para_env%group)
mixing_store%abroy(i, nv) = wfac
mixing_store%abroy(nv, i) = wfac
END DO

! broyden matrices
ALLOCATE (beta(nv, nv), cvec(nv), gammab(nv))
ALLOCATE (amat(nv, nv))
DO i = 1, nv
wfac = SUM(mixing_store%dfbroy(:, :, i)*dq_now(:, :))
CALL mp_sum(wfac, para_env%group)
cvec(i) = mixing_store%wbroy(i)*wfac
END DO

DO i = 1, nv
DO j = 1, nv
beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
END DO
beta(i, i) = beta(i, i)+omega0*omega0
END DO

rskip = 1.e-12_dp
CALL get_pseudo_inverse_svd(beta, amat, rskip)
gammab = MATMUL(cvec, amat)

! build ubroy
mixing_store%ubroy(:, :, nv) = alpha*mixing_store%dfbroy(:, :, nv)+wdf*(q_now(:, :)-q_last(:, :))

charges = 0.0_dp
DO ia = 1, mixing_store%nat_local
ii = mixing_store%atlist(ia)
charges(ii, 1:ns) = q_now(ia, 1:ns)+alpha*dq_now(ia, 1:ns)
END DO
DO i = 1, nv
DO ia = 1, mixing_store%nat_local
ii = mixing_store%atlist(ia)
charges(ii, 1:ns) = charges(ii, 1:ns)-mixing_store%wbroy(i)*gammab(i)*mixing_store%ubroy(ia, 1:ns, i)
END DO
END DO
CALL mp_sum(charges, para_env%group)

DEALLOCATE (beta, cvec, gammab)
DEALLOCATE (amat)

END SUBROUTINE broyden_mixing

END MODULE qs_charge_mixing
31 changes: 28 additions & 3 deletions src/qs_density_mixing_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,9 +75,14 @@ MODULE qs_density_mixing_types
REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm_res_buffer, pulay_matrix
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: fmat, gmat, smat
!
REAL(KIND=dp) :: nat_local, max_shell
REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: acharge
REAL(KIND=dp), DIMENSION(:), POINTER :: atlist
INTEGER :: nat_local, max_shell
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: acharge
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dacharge
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dfbroy
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ubroy
REAL(KIND=dp), DIMENSION(:, :), POINTER :: abroy
REAL(KIND=dp), DIMENSION(:), POINTER :: wbroy
INTEGER, DIMENSION(:), POINTER :: atlist
!
TYPE(cp_1d_z_p_type), DIMENSION(:), POINTER :: last_res, rhoin, rhoin_old
TYPE(cp_1d_z_p_type), DIMENSION(:, :), POINTER :: delta_res, u_vec, z_vec
Expand Down Expand Up @@ -141,6 +146,11 @@ SUBROUTINE mixing_storage_create(mixing_store, mixing_section, mixing_method, ec
NULLIFY (mixing_store%gmat)
NULLIFY (mixing_store%smat)
NULLIFY (mixing_store%acharge)
NULLIFY (mixing_store%dacharge)
NULLIFY (mixing_store%dfbroy)
NULLIFY (mixing_store%ubroy)
NULLIFY (mixing_store%abroy)
NULLIFY (mixing_store%wbroy)
NULLIFY (mixing_store%atlist)
NULLIFY (mixing_store%last_res)
NULLIFY (mixing_store%rhoin)
Expand Down Expand Up @@ -397,6 +407,21 @@ SUBROUTINE mixing_storage_release(mixing_store)
IF (ASSOCIATED(mixing_store%acharge)) THEN
DEALLOCATE (mixing_store%acharge)
END IF
IF (ASSOCIATED(mixing_store%dacharge)) THEN
DEALLOCATE (mixing_store%dacharge)
END IF
IF (ASSOCIATED(mixing_store%dfbroy)) THEN
DEALLOCATE (mixing_store%dfbroy)
END IF
IF (ASSOCIATED(mixing_store%ubroy)) THEN
DEALLOCATE (mixing_store%ubroy)
END IF
IF (ASSOCIATED(mixing_store%abroy)) THEN
DEALLOCATE (mixing_store%abroy)
END IF
IF (ASSOCIATED(mixing_store%wbroy)) THEN
DEALLOCATE (mixing_store%wbroy)
END IF
IF (ASSOCIATED(mixing_store%atlist)) THEN
DEALLOCATE (mixing_store%atlist)
END IF
Expand Down
17 changes: 13 additions & 4 deletions src/qs_mixing_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,17 +209,26 @@ SUBROUTINE mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mi
END DO
END DO
IF (ASSOCIATED(mixing_store%acharge)) DEALLOCATE (mixing_store%acharge)
ALLOCATE (mixing_store%acharge(na, max_shell, nbuffer, nspins))
ALLOCATE (mixing_store%acharge(na, max_shell, nbuffer))
IF (ASSOCIATED(mixing_store%dacharge)) DEALLOCATE (mixing_store%dacharge)
ALLOCATE (mixing_store%dacharge(na, max_shell, nbuffer))
END IF
IF (mixing_method == pulay_mixing_nr) THEN
IF (ASSOCIATED(mixing_store%pulay_matrix)) DEALLOCATE (mixing_store%pulay_matrix)
ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer))
ELSEIF (mixing_method == broyden_mixing_nr) THEN
!
IF (ASSOCIATED(mixing_store%abroy)) DEALLOCATE (mixing_store%abroy)
ALLOCATE (mixing_store%abroy(nbuffer, nbuffer))
IF (ASSOCIATED(mixing_store%wbroy)) DEALLOCATE (mixing_store%wbroy)
ALLOCATE (mixing_store%wbroy(nbuffer))
IF (ASSOCIATED(mixing_store%dfbroy)) DEALLOCATE (mixing_store%dfbroy)
ALLOCATE (mixing_store%dfbroy(na, max_shell, nbuffer))
IF (ASSOCIATED(mixing_store%ubroy)) DEALLOCATE (mixing_store%ubroy)
ALLOCATE (mixing_store%ubroy(na, max_shell, nbuffer))
ELSEIF (mixing_method == broyden_mixing_new_nr) THEN
!
CPABORT("broyden_mixing_new not available")
ELSEIF (mixing_method == multisecant_mixing_nr) THEN
!
CPABORT("multisecant_mixing not available")
END IF
ELSE
! *** allocate buffer for gspace mixing ***
Expand Down
8 changes: 8 additions & 0 deletions src/xtb_matrices.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ MODULE xtb_matrices
USE mulliken, ONLY: ao_charges
USE orbital_pointers, ONLY: ncoset
USE particle_types, ONLY: particle_type
USE qs_charge_mixing, ONLY: charge_mixing
USE qs_dispersion_pairpot, ONLY: d3_cnumber,&
d3_cnumber_bcast,&
dcnum_type
Expand Down Expand Up @@ -76,6 +77,7 @@ MODULE xtb_matrices
USE qs_overlap, ONLY: create_sab_matrix
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_scf_types, ONLY: qs_scf_env_type
USE virial_methods, ONLY: virial_pair_force
USE virial_types, ONLY: virial_type
USE xtb_coulomb, ONLY: build_xtb_coulomb
Expand Down Expand Up @@ -949,6 +951,7 @@ SUBROUTINE build_xtb_ks_matrix(qs_env, calculate_forces, just_energy)
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(qs_scf_env_type), POINTER :: scf_env
TYPE(section_vals_type), POINTER :: scf_section
TYPE(xtb_atom_type), POINTER :: xtb_kind

Expand Down Expand Up @@ -1019,6 +1022,11 @@ SUBROUTINE build_xtb_ks_matrix(qs_env, calculate_forces, just_energy)
END DO
DEALLOCATE (aocg)

! charge mixing
CALL get_qs_env(qs_env=qs_env, scf_env=scf_env)
CALL charge_mixing(qs_env, scf_env%mixing_method, scf_env%mixing_store, &
charges, para_env, scf_env%iter_count)

DO iatom = 1, natom
mcharge(iatom) = SUM(charges(iatom, :))
END DO
Expand Down
2 changes: 1 addition & 1 deletion src/xtb_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
! 3) Periodic displacement field (debugging)
! 4) Check for RTP and EMD
! 5) Wannier localization
! 6) Charge Mixing methods: Broyden/Pulay
! 6) Charge Mixing methods: Broyden/Pulay (more debugging needed, also add to DFTB)
! **************************************************************************************************
MODULE xtb_types

Expand Down

0 comments on commit 9e14bf3

Please sign in to comment.