diff --git a/src/DMMinModule.f90 b/src/DMMinModule.f90 index ac67afd30..10205443c 100644 --- a/src/DMMinModule.f90 +++ b/src/DMMinModule.f90 @@ -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 @@ -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) @@ -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 @@ -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) @@ -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 @@ -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 diff --git a/src/force_module.f90 b/src/force_module.f90 index 6077503fb..6c432db0c 100644 --- a/src/force_module.f90 +++ b/src/force_module.f90 @@ -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, & @@ -999,7 +1001,7 @@ 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, & @@ -1007,7 +1009,7 @@ subroutine pulay_force(p_force, KE_force, fixed_potential, vary_mu, & 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, & @@ -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 @@ -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 diff --git a/src/global_module.f90 b/src/global_module.f90 index 3b521481e..a4eef0650 100644 --- a/src/global_module.f90 +++ b/src/global_module.f90 @@ -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 !!*** diff --git a/testsuite/test_002_bulk_Si_1proc_OrderN/Conquest_out.ref b/testsuite/test_002_bulk_Si_1proc_OrderN/Conquest_out.ref index d097422a4..db5fb4165 100644 --- a/testsuite/test_002_bulk_Si_1proc_OrderN/Conquest_out.ref +++ b/testsuite/test_002_bulk_Si_1proc_OrderN/Conquest_out.ref @@ -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 @@ -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 @@ -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 @@ -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