Skip to content

Implement correction to S-Pulay force and stress #409

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
davidbowler merged 2 commits into from
Mar 28, 2025
Merged

Implement correction to S-Pulay force and stress #409

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7caf918
Implement correction to S-Pulay force and stress
davidbowler Mar 20, 2025
7c425ff
Update reference output for linear scaling test
davidbowler Mar 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 5 additions & 2 deletions src/DMMinModule.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -827,7 +827,7 @@ subroutine lateDM(ndone, n_L_iterations, done, deltaE, vary_mu, &
flag_mix_L_SC_min, &
flag_fix_spin_population, nspin, &
spin_factor, flag_dump_L, &
flag_SpinDependentSF, min_layer
flag_SpinDependentSF, min_layer, mu_DMM
use timer_module, only: cq_timer,start_timer, &
stop_print_timer, WITH_LEVEL
use io_module, only: dump_matrix, return_prefix
Expand Down Expand Up @@ -877,6 +877,7 @@ subroutine lateDM(ndone, n_L_iterations, done, deltaE, vary_mu, &
if (ndone > n_L_iterations) &
call cq_abort('lateDM: too many L iterations', ndone, n_L_iterations)

mu_DMM = zero
do spin = 1, nspin
do i = 1, maxpulayDMM
mat_Lstore(i,spin) = allocate_temp_matrix(Lrange,0)
Expand All @@ -888,7 +889,6 @@ subroutine lateDM(ndone, n_L_iterations, done, deltaE, vary_mu, &
matSphi(spin) = allocate_temp_matrix(Lrange,0)
mat_temp(spin) = allocate_temp_matrix(TLrange,0)
end do

! Update the charge density if flag is set
min_layer = min_layer - 1
if (flag_mix_L_SC_min) then
Expand Down Expand Up @@ -931,6 +931,7 @@ subroutine lateDM(ndone, n_L_iterations, done, deltaE, vary_mu, &
call matrix_product(mat_temp(spin), matT(spin_SF), matSphi(spin), mult(TL_T_L))
e_dot_n(spin) = matrix_product_trace(matSM3(spin), matphi(spin))
n_dot_n(spin) = matrix_product_trace(matSphi(spin), matphi(spin))
mu_DMM(spin) = e_dot_n(spin) / n_dot_n(spin)
if (inode == ionode .and. iprint_DM + min_layer >= 3) then
write(io_lun, '(4x,a,i1,") ",f16.6)') &
trim(prefix)//" e.n (spin=", spin, e_dot_n(spin)
Expand Down Expand Up @@ -1057,6 +1058,7 @@ subroutine lateDM(ndone, n_L_iterations, done, deltaE, vary_mu, &
mult(TL_T_L))
e_dot_n(spin) = matrix_product_trace(matSM3(spin), matphi(spin))
n_dot_n(spin) = matrix_product_trace(matSphi(spin), matphi(spin))
mu_DMM(spin) = e_dot_n(spin) / n_dot_n(spin)
end do
if (flag_fix_spin_population) then
do spin = 1, nspin
Expand Down Expand Up @@ -1185,6 +1187,7 @@ subroutine lateDM(ndone, n_L_iterations, done, deltaE, vary_mu, &
matSphi(spin), mult(TL_T_L))
e_dot_n(spin) = matrix_product_trace(matSM3(spin), matphi(spin))
n_dot_n(spin) = matrix_product_trace(matSphi(spin), matphi(spin))
mu_DMM(spin) = e_dot_n(spin) / n_dot_n(spin)
end do
if (flag_fix_spin_population) then
do spin = 1, nspin
Expand Down
16 changes: 11 additions & 5 deletions src/force_module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -986,6 +986,8 @@ end subroutine force
!! Added calculations for off-diagonal elements of PP_stress and SP_stress
!! 2019/05/08 zamaan
!! Added atomic stress contributions
!! 2025/03/20 16:07 dave
!! Added electron number gradient contribution to S-Pulay force and stress
!! SOURCE
!!
subroutine pulay_force(p_force, KE_force, fixed_potential, vary_mu, &
Expand All @@ -999,15 +1001,15 @@ subroutine pulay_force(p_force, KE_force, fixed_potential, vary_mu, &
use matrix_module, only: matrix, matrix_halo
use matrix_data, only: mat, halo, blip_trans, Srange, aSa_range
use mult_module, only: LNV_matrix_multiply, &
matM12, &
matM12, matM4, &
allocate_temp_matrix, &
free_temp_matrix, &
return_matrix_value, &
matrix_pos, ltrans, &
scale_matrix_value, &
matKatomf, &
return_matrix_block_pos, &
matrix_scale, &
matrix_scale, matrix_sum, &
SF_to_AtomF_transform
use global_module, only: iprint_MD, WhichPulay, &
BothPulay, PhiPulay, &
Expand All @@ -1019,7 +1021,7 @@ subroutine pulay_force(p_force, KE_force, fixed_potential, vary_mu, &
id_glob, species_glob, &
flag_diagonalisation, &
flag_full_stress, flag_stress, &
flag_atomic_stress, min_layer
flag_atomic_stress, min_layer, mu_DMM
use set_bucket_module, only: rem_bucket, atomf_atomf_rem
use blip_grid_transform_module, only: blip_to_support_new, &
blip_to_grad_new
Expand Down Expand Up @@ -1162,8 +1164,12 @@ subroutine pulay_force(p_force, KE_force, fixed_potential, vary_mu, &
! If we're diagonalising, we've already build data_M12
if (.not. flag_diagonalisation) then
call LNV_matrix_multiply(electrons, energy_tmp, dontK, doM1, &
doM2, dontM3, dontM4, dontphi, dontE, &
mat_M12=matM12)
doM2, dontM3, doM4, dontphi, dontE, &
mat_M12=matM12, mat_M4=matM4)
! Electron number gradient contribution to S-Pulay force and stress
do spin=1,nspin
call matrix_sum(one,matM12(spin),-half*mu_DMM(spin),matM4(spin))
end do
end if
t1 = mtime()
t0 = t1
Expand Down
2 changes: 2 additions & 0 deletions src/global_module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -424,5 +424,7 @@ module global_module
integer :: i_pol_dir_st, i_pol_dir_end ! Either 1,1 or 1,3
integer, dimension(3) :: i_pol_dir ! Either n,0,0 or 1,2,3

! Density matrix Lagrange multiplier for correct electron number (needed for forces and stress)
real(double), dimension(2) :: mu_DMM ! Allow for spin
end module global_module
!!***
44 changes: 24 additions & 20 deletions testsuite/test_002_bulk_Si_1proc_OrderN/Conquest_out.ref
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
M.J.Gillan (Keele, UCL)
________________________________________________________________________

Simulation cell dimensions: 10.3600 a0 x 10.3600 a0 x 10.3600 a0
Simulation cell dimensions: 10.3600 a0 x 10.3600 a0 x 10.3600 a0

Atomic coordinates (a0)
Atom X Y Z Species
Expand All @@ -34,9 +34,9 @@
7 2.5900 7.7700 7.7700 1
8 7.7700 2.5900 7.7700 1

This job was run on 2023/05/23 at 15:11 +0100
Code was compiled on 2023/05/23 at 14:56 +0100
Version comment: Git Branch: master; tag, hash: v1.1
This job was run on 2025/03/21 at 14:30 +0000
Code was compiled on 2025/03/20 at 16:02 +0000
Version comment: Git Branch: bugfix-correct-force-stress; tag, hash: v1.4-39-g049ac20c

Job title:
Job to be run: static calculation
Expand All @@ -57,30 +57,34 @@
--------------------------------------------------------

The calculation will be performed on 1 process

The calculation will be performed on 1 thread

Using the default matrix multiplication kernel
Density Matrix range = 16.0000 a0

The functional used will be GGA PBE96


lateDM: reached residual of 0.687257E-07 after 29 iterations
lateDM: reached residual of 0.693388E-07 after 29 iterations

| Number of electrons = 32.000031
|* Harris-Foulkes energy = -33.569389509372144 Ha
|* Harris-Foulkes energy = -33.569389505626887 Ha

force: Forces on atoms (Ha/a0)
force: Atom X Y Z
force: 1 -0.0021719988 -0.0040582799 -0.0060871907
force: 2 0.0000826094 0.0001737196 -0.0005472291
force: 3 0.0000821377 -0.0003627707 0.0002616019
force: 4 -0.0001763387 0.0001709996 0.0002604078
force: 5 0.0026208445 0.0027153199 0.0028116541
force: 6 0.0001194484 0.0002209308 0.0003000775
force: 7 -0.0014819902 0.0018086628 0.0019152353
force: 8 0.0009217186 -0.0005998024 0.0011300341

force: Maximum force : 0.00608719(Ha/a0) on atom 1 in z direction
force: Force Residual: 0.00340652 Ha/a0
force: Total stress: -9.65909297 -9.66104273 -9.66437275 GPa
force: 1 -0.0021132141 -0.0039395925 -0.0059091849
force: 2 0.0000893280 0.0001876750 -0.0004909937
force: 3 0.0000887525 -0.0003253930 0.0002828883
force: 4 -0.0001578731 0.0001849757 0.0002817526
force: 5 0.0025536669 0.0026340003 0.0027161702
force: 6 0.0001055402 0.0001928153 0.0002577064
force: 7 -0.0014596011 0.0017440213 0.0018365365
force: 8 0.0008896896 -0.0006100493 0.0010697559

force: Maximum force : 0.00590918(Ha/a0) on atom 1 in z direction
force: Force Residual: 0.00330389 Ha/a0
force: Total stress: -6.71161888 -6.71375049 -6.71725660 GPa

BIBLIOGRAPHY: Please consider citing the following references in the conquest.bib file

Expand All @@ -90,6 +94,6 @@
Pseudopotentials: Hamann2013, Bowler2019
XC functional: Perdew1996

Max total mem use is 105.637 MB
Max total mem use is 105.403 MB

Total run time was: 26.736 seconds
Total run time was: 51.565 seconds