Skip to content

Commit

Permalink
Fixed CDFT force bug
Browse files Browse the repository at this point in the history
  • Loading branch information
@chrisahart @schulkov
chrisahart authored and schulkov committed Oct 10, 2024
1 parent 6733ffb commit ae45d9a
Showing 1 changed file with 16 additions and 12 deletions.
28 changes: 16 additions & 12 deletions src/qs_cdft_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -683,7 +683,6 @@ SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
IF (PRESENT(just_gradients)) my_just_gradients = just_gradients
IF (my_just_gradients) THEN
cdft_control%in_memory = .TRUE.
cdft_control%atomic_charges = .FALSE.
hirshfeld_control%print_density = .FALSE.
END IF

Expand Down Expand Up @@ -756,8 +755,8 @@ SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
CALL pw_set(cdft_control%group(igroup)%hw_rho_total_constraint, 1.0_dp)

IF (igroup == 1) THEN
CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total)
CALL pw_set(cdft_control%group(igroup)%hw_rho_total, 1.0_dp)
CALL auxbas_pw_pool%create_pw(cdft_control%group(1)%hw_rho_total)
CALL pw_set(cdft_control%group(1)%hw_rho_total, 1.0_dp)

IF (hirshfeld_control%print_density) THEN
DO iatom = 1, cdft_control%natoms
Expand Down Expand Up @@ -825,7 +824,6 @@ SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
pab=pab, o1=0, o2=0, & ! without map_consistent
prefactor=1.0_dp, cutoff=0.0_dp)
! IF (para_env%is_source()) PRINT *, "radius", radius
END IF

IF (igroup == 1) THEN
Expand Down Expand Up @@ -880,7 +878,7 @@ SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
DEALLOCATE (pab)

IF (igroup == 1) THEN
CALL transfer_rs2pw(rs_rho_all, cdft_control%group(igroup)%hw_rho_total)
CALL transfer_rs2pw(rs_rho_all, cdft_control%group(1)%hw_rho_total)
END IF

CALL transfer_rs2pw(rs_rho_constr, cdft_control%group(igroup)%hw_rho_total_constraint)
Expand All @@ -898,7 +896,7 @@ SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
CALL transfer_rs2pw(rs_single_charge(i), cdft_control%group(igroup)%hw_rho_atomic_charge(i))
CALL hfun_scale(cdft_control%charge(i)%array, &
cdft_control%group(igroup)%hw_rho_atomic_charge(i)%array, &
cdft_control%group(igroup)%hw_rho_total%array, divide=.TRUE., &
cdft_control%group(1)%hw_rho_total%array, divide=.TRUE., &
small=hirshfeld_control%eps_cutoff)
END DO
END IF
Expand Down Expand Up @@ -1080,18 +1078,27 @@ SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
IF (cdft_control%in_memory) THEN
DO igroup = 1, SIZE(cdft_control%group)
! Coefficients
coefficients(:) = 0.0_dp
is_constraint = .FALSE.
DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
is_constraint(cdft_control%group(igroup)%atoms(i)) = .TRUE.
END DO
DO k = lb_pw(3), ub_pw(3)
DO j = lb_pw(2), ub_pw(2)
DO i = lb_pw(1), ub_pw(1)
DO iatom = 1, natom
ra(:) = particle_set(iatom)%r
IF (cdft_control%group(igroup)%hw_rho_total%array(i, j, k) > hirshfeld_control%eps_cutoff) THEN
IF (cdft_control%group(1)%hw_rho_total%array(i, j, k) > hirshfeld_control%eps_cutoff) THEN
exp_eval = (coefficients(iatom) - &
cdft_control%group(igroup)%weight%array(i, j, k))/ &
cdft_control%group(igroup)%hw_rho_total%array(i, j, k)
cdft_control%group(1)%hw_rho_total%array(i, j, k)
r2 = (/i*dr_pw(1), j*dr_pw(2), k*dr_pw(3)/) + origin
r_pbc = pbc(ra, r2, cell)
Expand All @@ -1116,11 +1123,8 @@ SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
END DO
END DO
END DO
IF (igroup == 1) THEN
CALL auxbas_pw_pool%give_back_pw(cdft_control%group(1)%hw_rho_total)
END IF
END DO
CALL auxbas_pw_pool%give_back_pw(cdft_control%group(1)%hw_rho_total)
END IF
CALL rs_grid_release(rs_rho_all)
Expand Down

0 comments on commit ae45d9a

Please sign in to comment.