From 5cedbc1402073fa762a85e0a81467d8e83aa561c Mon Sep 17 00:00:00 2001 From: David Bowler Date: Mon, 24 Jul 2023 16:10:30 +0100 Subject: [PATCH 1/3] Updates to non-SCF NA stress Fixing bugs for non-SCF stress with neutral atom (mainly a problem with the formulation rather than code) --- src/force_module.f90 | 73 +++++++++++++++++++++----------------------- 1 file changed, 34 insertions(+), 39 deletions(-) diff --git a/src/force_module.f90 b/src/force_module.f90 index c6df3aff..d86d312c 100644 --- a/src/force_module.f90 +++ b/src/force_module.f90 @@ -3414,6 +3414,8 @@ contains !! Bugfix: correct calculation of rsq implemented, variables tidied !! 2023/07/19 10:40 dave !! Added minus sign to non-SCF PCC force for consistency with non-SCF part + !! 2023/07/24 16:00 dave + !! Bug fix for non-SCF NA stress (remove Hartree stress terms) !! SOURCE !! subroutine get_nonSC_correction_force(HF_force, density_out, inode, & @@ -3428,7 +3430,7 @@ contains area_moveatoms, IPRINT_TIME_THRES3, & flag_pcc_global, nspin, spin_factor, & flag_full_stress, flag_stress, & - flag_atomic_stress + flag_atomic_stress, flag_neutral_atom use XC, only: get_xc_potential, & get_dxc_potential, & flag_is_GGA @@ -3442,7 +3444,7 @@ contains use atomic_density, only: atomic_density_table use pseudo_tm_info, only: pseudo use dimens, only: grid_point_volume, n_my_grid_points - use GenBlas, only: axpy + use GenBlas, only: axpy, dot use density_module, only: density, density_scale, density_pcc use hartree_module, only: hartree, Hartree_stress use potential_module, only: potential @@ -3490,10 +3492,9 @@ contains real(double), dimension(nspin) :: pot_here, pot_here_pcc ! only for GGA with P.C.C. real(double), allocatable, dimension(:) :: h_potential_in, & - wk_grid_total, & - density_out_GGA_total + wk_grid_total real(double), allocatable, dimension(:,:) :: wk_grid, & - density_out_GGA + density_out_GGA !****lat<$ @@ -3563,10 +3564,16 @@ contains jacobian = zero ! use density_total (input charge) WITHOUT a factor of half so that this term corrects the GPV ! found in the main force routine - do ipoint = 1,nsize - jacobian = jacobian + (density_out_total(ipoint) - density_total(ipoint))* & - (h_potential(ipoint) + pseudopotential(ipoint)) ! Also calculate the correction to GPV for local potential - end do + jacobian = jacobian + dot(nsize,density_out_total,1,pseudopotential,1) - & + dot(nsize,density_total,1,pseudopotential,1) + if(flag_neutral_atom) then + ! These should be zero so ensure that they are; also don't include h_potential in jacobian + Hartree_stress = zero + loc_stress = zero + else + jacobian = jacobian + dot(nsize,density_out_total,1,h_potential,1) - & + dot(nsize,density_total,1,h_potential,1) + end if ! Don't apply gsum to jacobian because nonSCF_stress will be summed (end of this routine) jacobian = jacobian*grid_point_volume ! loc_stress is \int n^{out} V^{PAD}_{Har} @@ -3586,13 +3593,11 @@ contains ! DeltaXC is added in the main force routine ! For PCC we will do this in the PCC force routine (easier) if (.NOT.flag_pcc_global) then - call get_xc_potential(density, dVxc_drho(:,:,1), & - potential(:,1), y_pcc, nsize) + call get_xc_potential(density, potential(:,:), & + dVxc_drho(:,1,1), h_energy, nsize) ! NB dVxc_drho is a dummy here jacobian = zero do spin = 1, nspin - do ipoint = 1,nsize - jacobian = jacobian + spin_factor*density_out(ipoint,spin)*dVxc_drho(ipoint,spin,1) - end do + jacobian = jacobian + spin_factor*dot(nsize,density_out(:,spin),1,potential(:,spin),1) end do jacobian = jacobian*grid_point_volume call gsum(jacobian) ! gsum as XC_stress isn't summed elsewhere @@ -3612,6 +3617,7 @@ contains call stop_timer(tmr_l_tmp2, TIME_ACCUMULATE_NO) ! for P.C.C. + call start_timer (tmr_l_tmp1, WITH_LEVEL) if (flag_pcc_global) then allocate(wk_grid_total(nsize), wk_grid(nsize,nspin), STAT=stat) wk_grid_total = zero @@ -3626,15 +3632,13 @@ contains end do ! only for GGA if (flag_is_GGA) then - allocate(density_out_GGA_total(nsize), density_out_GGA(nsize,nspin), STAT=stat) + allocate(density_out_GGA(nsize,nspin), STAT=stat) if (stat /= 0) call cq_abort ('Error allocating & &density_out_GGAs in get_nonSC_force ', stat) - call reg_alloc_mem(area_moveatoms, (nspin + 1) * nsize, type_dbl) - density_out_GGA_total = zero + call reg_alloc_mem(area_moveatoms, nspin * nsize, type_dbl) density_out_GGA = zero do spin = 1, nspin density_out_GGA(:,spin) = density_out(:,spin) + half * density_pcc(:) - density_out_GGA_total(:) = density_out_GGA_total(:) + spin_factor * density_out_GGA(:,spin) end do ! copy hartree potential allocate(h_potential_in(nsize), STAT=stat) @@ -3643,15 +3647,18 @@ contains &get_nonSC_force ', stat) call reg_alloc_mem(area_moveatoms, nsize, type_dbl) h_potential_in = h_potential - end if - end if - - call start_timer (tmr_l_tmp1, WITH_LEVEL) - if (flag_pcc_global) then - if(flag_is_GGA) then call get_dxc_potential(wk_grid, dVxc_drho, nsize, density_out_GGA) ! GGA with spin not implemented ! potential(:,1) = potential(:,1) + dVxc_drho(:,1,1) + deallocate(density_out_GGA, STAT=stat) + if (stat /= 0) call cq_abort('Error deallocating density_out_GGAs in & + &get_nonSC_force ', stat) + call reg_dealloc_mem(area_moveatoms, nspin * nsize, type_dbl) + ! make a copy of potential at this point + ! use wk_grid as a temporary storage + do spin = 1, nspin + wk_grid(:,spin) = potential(:,spin) + end do else call get_dxc_potential(wk_grid, dVxc_drho, nsize) end if @@ -3664,22 +3671,8 @@ contains call get_dxc_potential(density, dVxc_drho, nsize) end if end if - ! deallocating density_out_GGA: only for P.C.C. - if (flag_pcc_global) then - if (flag_is_GGA) then - deallocate(density_out_GGA_total, density_out_GGA, STAT=stat) - if (stat /= 0) call cq_abort('Error deallocating density_out_GGAs in & - &get_nonSC_force ', stat) - call reg_dealloc_mem(area_moveatoms, (nspin + 1) * nsize, type_dbl) - ! make a copy of potential at this point - ! use wk_grid as a temporary storage - do spin = 1, nspin - wk_grid(:,spin) = potential(:,spin) - end do - end if - end if !flag_pcc_global - ! for LDA + ! for LDA: find dn* dxc_potential if (.NOT.(flag_is_GGA)) then do spin = 1, nspin do spin_2 = 1, nspin @@ -3698,12 +3691,14 @@ contains ! Restart of the timer; level assigned here call start_timer(tmr_l_tmp2, WITH_LEVEL) + ! Calculate the Hartree potential from the output density h_potential = zero ! Preserve the Hartree stress we've calculated loc_stress = Hartree_stress call hartree(density_out_total, h_potential, nsize, h_energy) ! And restore Hartree_stress = loc_stress + ! And subtract from potential so that we have delta V_Ha do spin = 1, nspin call axpy(nsize, -one, h_potential, 1, potential(:,spin), 1) end do From 5de2765181114b5aa90a8a04b44b13d6b8a5e587 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Tue, 25 Jul 2023 15:10:50 +0100 Subject: [PATCH 2/3] Updates for non-SCF stress Mainly adding GGA XC contribution but also tidying and simplifying non-SCF force routine to enable checking of formulae. This appears to have a small error in the stress when using GGA with non-SCF (and probably with PCC) at the level of 0.15GPa which I'm unable to find at the moment --- src/XC_CQ_module.f90 | 1 + src/XC_LibXC_v4_module.f90 | 1 + src/XC_LibXC_v5_module.f90 | 1 + src/force_module.f90 | 116 +++++++++++-------------------------- 4 files changed, 36 insertions(+), 83 deletions(-) diff --git a/src/XC_CQ_module.f90 b/src/XC_CQ_module.f90 index 6c86948b..bf824feb 100644 --- a/src/XC_CQ_module.f90 +++ b/src/XC_CQ_module.f90 @@ -229,6 +229,7 @@ contains ! Local variables real(double) :: loc_x_energy, exx_tmp + XC_GGA_stress = zero select case(flag_functional_type) case (functional_lda_pz81) ! NOT SPIN POLARISED diff --git a/src/XC_LibXC_v4_module.f90 b/src/XC_LibXC_v4_module.f90 index 8cf471d2..7202924e 100644 --- a/src/XC_LibXC_v4_module.f90 +++ b/src/XC_LibXC_v4_module.f90 @@ -366,6 +366,7 @@ contains ! Local variables real(double) :: loc_x_energy, exx_tmp + XC_GGA_stress = zero if(flag_use_libxc) then call get_libxc_potential(density=density, size=size,& xc_potential =xc_potential, & diff --git a/src/XC_LibXC_v5_module.f90 b/src/XC_LibXC_v5_module.f90 index aefb76e9..7fff7752 100644 --- a/src/XC_LibXC_v5_module.f90 +++ b/src/XC_LibXC_v5_module.f90 @@ -370,6 +370,7 @@ contains ! Local variables real(double) :: loc_x_energy, exx_tmp + XC_GGA_stress = zero if(flag_use_libxc) then call get_libxc_potential(density=density, size=size,& xc_potential =xc_potential, & diff --git a/src/force_module.f90 b/src/force_module.f90 index d86d312c..ac012b9e 100644 --- a/src/force_module.f90 +++ b/src/force_module.f90 @@ -391,11 +391,11 @@ contains GPV_stress(dir1,dir1) = (hartree_energy_total_rho + & local_ps_energy - core_correction) ! core contains 1/V term end if + ! XC_GGA_stress is zero for LDA if(flag_self_consistent .or. flag_mix_L_SC_min) then - XC_stress(dir1,dir1) = xc_energy + & - spin_factor*XC_GGA_stress(dir1,dir1) + XC_stress(dir1,dir1) = xc_energy + spin_factor*XC_GGA_stress(dir1,dir1) else ! nonSCF XC found later, along with corrections to Hartree - XC_stress(dir1,dir1) = delta_E_xc !xc_energy + spin_factor*XC_GGA_stress(direction) + XC_stress(dir1,dir1) = delta_E_xc + spin_factor*XC_GGA_stress(dir1,dir1) end if end do end if @@ -3418,7 +3418,7 @@ contains !! Bug fix for non-SCF NA stress (remove Hartree stress terms) !! SOURCE !! - subroutine get_nonSC_correction_force(HF_force, density_out, inode, & + subroutine get_nonSC_correction_force(NSCforce, density_out, inode, & ionode, n_atoms, nsize) use datatypes @@ -3462,7 +3462,7 @@ contains integer :: n_atoms, nsize integer :: inode, ionode real(double), dimension(:,:) :: density_out - real(double), dimension(:,:) :: HF_force + real(double), dimension(:,:) :: NSCforce ! Local variables character(len=80) :: sub_name = "get_nonSC_correction_force" @@ -3485,16 +3485,11 @@ contains type(cq_timer) :: backtrace_timer real(double), dimension(3,3) :: loc_stress - real(double), dimension(:), allocatable :: h_potential, & - density_total, & - density_out_total + real(double), dimension(:), allocatable :: h_potential, density_total, density_out_total real(double), dimension(:,:,:), allocatable :: dVxc_drho real(double), dimension(nspin) :: pot_here, pot_here_pcc ! only for GGA with P.C.C. - real(double), allocatable, dimension(:) :: h_potential_in, & - wk_grid_total - real(double), allocatable, dimension(:,:) :: wk_grid, & - density_out_GGA + real(double), allocatable, dimension(:,:) :: wk_grid, density_out_GGA !****lat<$ @@ -3503,36 +3498,24 @@ contains ! Spin-polarised PBE non-SCF forces not implemented, so exit if necessary if ((nspin == 2) .and. flag_is_GGA) then ! Only true for CQ not LibXC call cq_warn(sub_name, "NonSCF forces not implemented for spin and GGA; these will be set to zero.") - HF_force = zero + NSCforce = zero return end if stat = 0 call start_timer(tmr_std_allocation) allocate(h_potential(nsize), STAT=stat) - if (stat /= 0) & - call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize) + if (stat /= 0) call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize) allocate(density_total(nsize), STAT=stat) - if (stat /= 0) & - call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize) + if (stat /= 0) call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize) allocate(density_out_total(nsize), STAT=stat) - if (stat /= 0) & - call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize) + if (stat /= 0) call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize) allocate(dVxc_drho(nsize,nspin,nspin), STAT=stat) - if (stat /= 0) & - call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize) + if (stat /= 0) call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize) call reg_alloc_mem(area_moveatoms, (3+nspin*nspin)*nsize, type_dbl) call stop_timer(tmr_std_allocation) -!****lat<$ - !print*, size(density_total, dim=1) - !print*, size(dVxc_drho, dim=1) - !print*, size(density_out_total,dim=1) - !print*, size(potential, dim=1) - !print*, size(density_out_total,dim=1) -!****lat>$ - - HF_force = zero + NSCforce = zero dcellx_block = rcellx / blocks%ngcellx dcelly_block = rcelly / blocks%ngcelly @@ -3619,34 +3602,25 @@ contains ! for P.C.C. call start_timer (tmr_l_tmp1, WITH_LEVEL) if (flag_pcc_global) then - allocate(wk_grid_total(nsize), wk_grid(nsize,nspin), STAT=stat) - wk_grid_total = zero + allocate(wk_grid(nsize,nspin), STAT=stat) wk_grid = zero if (stat /= 0) & call cq_abort('Error allocating wk_grids in & &get_nonSC_correction ', stat) - call reg_alloc_mem(area_moveatoms, (nspin + 1) * nsize, type_dbl) + call reg_alloc_mem(area_moveatoms, nspin * nsize, type_dbl) do spin = 1, nspin wk_grid(:,spin) = density(:,spin) + half * density_pcc(:) - wk_grid_total(:) = wk_grid_total(:) + spin_factor * wk_grid(:,spin) end do - ! only for GGA + ! Find dVxc_drho (different for LDA and GGA) if (flag_is_GGA) then allocate(density_out_GGA(nsize,nspin), STAT=stat) if (stat /= 0) call cq_abort ('Error allocating & &density_out_GGAs in get_nonSC_force ', stat) call reg_alloc_mem(area_moveatoms, nspin * nsize, type_dbl) - density_out_GGA = zero + density_out_GGA = zero do spin = 1, nspin density_out_GGA(:,spin) = density_out(:,spin) + half * density_pcc(:) end do - ! copy hartree potential - allocate(h_potential_in(nsize), STAT=stat) - if (stat /= 0) & - call cq_abort('Error allocating h_potential_in in & - &get_nonSC_force ', stat) - call reg_alloc_mem(area_moveatoms, nsize, type_dbl) - h_potential_in = h_potential call get_dxc_potential(wk_grid, dVxc_drho, nsize, density_out_GGA) ! GGA with spin not implemented ! potential(:,1) = potential(:,1) + dVxc_drho(:,1,1) @@ -3654,11 +3628,6 @@ contains if (stat /= 0) call cq_abort('Error deallocating density_out_GGAs in & &get_nonSC_force ', stat) call reg_dealloc_mem(area_moveatoms, nspin * nsize, type_dbl) - ! make a copy of potential at this point - ! use wk_grid as a temporary storage - do spin = 1, nspin - wk_grid(:,spin) = potential(:,spin) - end do else call get_dxc_potential(wk_grid, dVxc_drho, nsize) end if @@ -3693,10 +3662,8 @@ contains call start_timer(tmr_l_tmp2, WITH_LEVEL) ! Calculate the Hartree potential from the output density h_potential = zero - ! Preserve the Hartree stress we've calculated - loc_stress = Hartree_stress + loc_stress = Hartree_stress ! Preserve Hartree stress call hartree(density_out_total, h_potential, nsize, h_energy) - ! And restore Hartree_stress = loc_stress ! And subtract from potential so that we have delta V_Ha do spin = 1, nspin @@ -3804,8 +3771,8 @@ contains ! on what we are doing here. do spin = 1, nspin do dir1 = 1, 3 - HF_force(dir1,ig_atom) = & - HF_force(dir1,ig_atom) + spin_factor * & + NSCforce(dir1,ig_atom) = & + NSCforce(dir1,ig_atom) + spin_factor * & fr_1(dir1,spin) * pot_here(spin) if (flag_stress) then if (flag_full_stress) then @@ -3848,23 +3815,16 @@ contains ! for GGA potential = zero do spin = 1, nspin - do i = 1, n_my_grid_points - ! -delta n * dxc_potential - potential(i,spin) = wk_grid(i,spin) - h_potential_in(i) - end do + potential(:,spin) = dVxc_drho(:,spin,spin) end do else ! For LDA potential = zero do spin = 1, nspin do spin_2 = 1, nspin - do i = 1, n_my_grid_points - potential(i,spin) = & - potential(i,spin) + & - spin_factor * & - (density(i,spin_2) - density_out(i,spin_2)) * & - dVxc_drho(i,spin_2,spin) - end do + potential(:,spin) = potential(:,spin) + & + spin_factor * (density(:,spin_2) - density_out(:,spin_2)) * & + dVxc_drho(:,spin_2,spin) end do end do end if @@ -3972,8 +3932,8 @@ contains ! components) do spin = 1, nspin do dir1 = 1, 3 - HF_force(dir1,ig_atom) = & - HF_force(dir1,ig_atom) + spin_factor * & + NSCforce(dir1,ig_atom) = & + NSCforce(dir1,ig_atom) + spin_factor * & fr_pcc(dir1,spin) * pot_here_pcc(spin) if (flag_stress) then if (flag_full_stress) then @@ -4009,7 +3969,7 @@ contains end if ! (flag_pcc_global) call start_timer(tmr_l_tmp1, WITH_LEVEL) - call gsum(HF_force, 3, n_atoms) + call gsum(NSCforce, 3, n_atoms) if (flag_stress) call gsum(nonSCF_stress,3,3) call stop_print_timer(tmr_l_tmp1, "NSC force - Compilation", & IPRINT_TIME_THRES3) @@ -4017,18 +3977,11 @@ contains ! deallocating temporary arrays call start_timer(tmr_std_allocation) if (flag_pcc_global) then - deallocate(wk_grid_total, wk_grid, STAT=stat) + deallocate(wk_grid, STAT=stat) if (stat /= 0) & call cq_abort('Error deallocating wk_grid in & &get_nonSC_correction_force ', stat) call reg_dealloc_mem(area_moveatoms, (nspin + 1) * nsize, type_dbl) - if (flag_is_GGA) then - deallocate(h_potential_in, STAT=stat) - if (stat /= 0) & - call cq_abort('Error deallocating h_potential_in in & - &get_nonSC_correction_force ', stat) - call reg_dealloc_mem(area_moveatoms, nsize, type_dbl) - end if ! for GGA end if ! flag_pcc_global deallocate(h_potential, density_total, density_out_total, dVxc_drho, & STAT=stat) @@ -4169,15 +4122,14 @@ contains real(double), dimension(nspin) :: pot_here_pcc real(double), dimension(3) :: r_pcc, fr_pcc, r ! allocatable arrays - real(double), dimension(:), allocatable :: xc_epsilon, density_wk_tot + real(double), dimension(:), allocatable :: xc_epsilon real(double), dimension(:,:), allocatable :: xc_potential, density_wk type(cq_timer) :: backtrace_timer call start_backtrace(t=backtrace_timer,who='get_PCC_force',where=7,level=3,echo=.true.) - allocate(xc_epsilon(size), density_wk_tot(size), & - xc_potential(size,nspin), density_wk(size,nspin), STAT=stat) + allocate(xc_epsilon(size), xc_potential(size,nspin), density_wk(size,nspin), STAT=stat) if (stat /= 0) call cq_abort("get_pcc_force: Error alloc mem: ", size) - call reg_alloc_mem(area_moveatoms, (2+2*nspin)*size, type_dbl) + call reg_alloc_mem(area_moveatoms, (1+2*nspin)*size, type_dbl) ! initialise arrays pcc_force = zero @@ -4194,10 +4146,8 @@ contains call start_timer (tmr_l_tmp2) density_wk = zero - density_wk_tot = zero do spin = 1, nspin density_wk(:,spin) = density(:,spin) + half * density_pcc(:) - density_wk_tot(:) = density_wk_tot(:) + spin_factor * density_wk(:,spin) end do call get_xc_potential(density_wk, xc_potential, & @@ -4360,9 +4310,9 @@ contains call stop_print_timer(tmr_l_tmp1, "PCC force - Compilation", & IPRINT_TIME_THRES3) - deallocate(xc_epsilon, density_wk_tot, xc_potential, density_wk, STAT=stat) + deallocate(xc_epsilon, xc_potential, density_wk, STAT=stat) if (stat /= 0) call cq_abort("get_pcc_force: Error dealloc mem") - call reg_dealloc_mem(area_moveatoms, (2+2*nspin)*size, type_dbl) + call reg_dealloc_mem(area_moveatoms, (1+2*nspin)*size, type_dbl) call stop_backtrace(t=backtrace_timer,who='get_PCC_force',echo=.true.) return From 4081dfa00fba8a8da55adfa7fb7a66b342dcf3f6 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Wed, 26 Jul 2023 12:24:21 +0100 Subject: [PATCH 3/3] Further updates to stress These changes are actually for non-SCF GGA calculations, but fall under the general umbrella of stress bug fixes. Note that the XC_GGA_stress (which comes from changes in the XC energy related to cell-induced changes in the density gradient) are not quite correctly implemented. For non-SCF there should be a mix of input and output densities (and Exc and Vxc) used for this term, but this is very complex to implement, so for now we have an approximation: we mix the stress contribution found using input and output densities. The proportion can be set using an input flag (but should not really be changed from 0.5) --- docs/input_tags.rst | 10 ++++++ src/force_module.f90 | 69 +++++++++++++++++++++++++++++++++---- src/initial_read_module.f90 | 6 ++++ 3 files changed, 78 insertions(+), 7 deletions(-) diff --git a/docs/input_tags.rst b/docs/input_tags.rst index 5ed3c9b7..e14203d5 100644 --- a/docs/input_tags.rst +++ b/docs/input_tags.rst @@ -1403,6 +1403,16 @@ General.GapThreshold (*real*) General.only_Dispersion (*boolean*) Selects only DFT\_D2 calculation (no electronic structure etc) +General.MixXCGGAInOut (*real*) + For non-SCF calculations only, chooses how to mix the proportions of + GGA XC stress contribution (from the change of the electron density + gradient) found using input (0.0 gives pure input) and output (1.0 + gives pure output) densities. Note that this is an approximation but + varying the value significantly away from 0.5 will give inconsistency + between stress and energy. + + *default*: 0.5 + Go to :ref:`top `. .. _advanced_atomic_spec_tags: diff --git a/src/force_module.f90 b/src/force_module.f90 index ac012b9e..6c5499ea 100644 --- a/src/force_module.f90 +++ b/src/force_module.f90 @@ -106,6 +106,9 @@ module force_module PP_stress, GPV_stress, XC_stress, & nonSCF_stress, pcc_stress, NA_stress + ! How much we mix input and output densities for the calculations of XC_GGA_stress + ! during non-SCF simulations (defaults to half for averaging) + real(double) :: mix_input_output_XC_GGA_stress ! Useful parameters for selecting force calculations in NL part integer, parameter :: HF = 1 integer, parameter :: Pulay = 2 @@ -394,13 +397,12 @@ contains ! XC_GGA_stress is zero for LDA if(flag_self_consistent .or. flag_mix_L_SC_min) then XC_stress(dir1,dir1) = xc_energy + spin_factor*XC_GGA_stress(dir1,dir1) - else ! nonSCF XC found later, along with corrections to Hartree - XC_stress(dir1,dir1) = delta_E_xc + spin_factor*XC_GGA_stress(dir1,dir1) + else ! The rest of nonSCF XC found later (nonSC or PCC routines) + XC_stress(dir1,dir1) = delta_E_xc end if - end do + end do end if WhichPulay = BothPulay - ! matK->matKatomf backtransformation for contracted SFs if (atomf.ne.sf) then do ispin = 1, nspin @@ -761,7 +763,6 @@ contains volume = rcellx*rcelly*rcellz ! We need pressure in GPa, and only diagonal terms output scale = -HaBohr3ToGPa/volume - !call print_stress("Total pressure: ", stress*scale, 0) if(inode==ionode.AND.iprint_MD + min_layer>=1) & write(io_lun,'(4x,a,f15.8,a4/)') trim(prefix)//" Average pressure: ", & third*scale*(stress(1,1) + stress(2,2) + stress(3,3))," GPa" @@ -3455,6 +3456,7 @@ contains WITH_LEVEL, TIME_ACCUMULATE_NO, & TIME_ACCUMULATE_YES use pseudopotential_common, only: pseudopotential + use XC, ONLY: XC_GGA_stress implicit none @@ -3576,8 +3578,30 @@ contains ! DeltaXC is added in the main force routine ! For PCC we will do this in the PCC force routine (easier) if (.NOT.flag_pcc_global) then + ! Find XC_GGA_stress for density_out - we will average the input and output + ! as an approximation to the correct (hideously complex) stress so add factor of half + call get_xc_potential(density_out, potential(:,:), & + dVxc_drho(:,1,1), h_energy, nsize) ! NB dVxc_drho is a dummy here + if (flag_stress) then + XC_stress(1,1) = XC_stress(1,1) + spin_factor*XC_GGA_stress(1,1) * & + mix_input_output_XC_GGA_stress + XC_stress(2,2) = XC_stress(2,2) + spin_factor*XC_GGA_stress(2,2) * & + mix_input_output_XC_GGA_stress + XC_stress(3,3) = XC_stress(3,3) + spin_factor*XC_GGA_stress(3,3) * & + mix_input_output_XC_GGA_stress + end if + ! Find XC potential for input density call get_xc_potential(density, potential(:,:), & - dVxc_drho(:,1,1), h_energy, nsize) ! NB dVxc_drho is a dummy here + dVxc_drho(:,1,1), h_energy, nsize) ! NB dVxc_drho is a dummy here + ! And now add XC_GGA_stress for density_in scaled by half + if (flag_stress) then + XC_stress(1,1) = XC_stress(1,1) + spin_factor*XC_GGA_stress(1,1) * & + (one - mix_input_output_XC_GGA_stress) + XC_stress(2,2) = XC_stress(2,2) + spin_factor*XC_GGA_stress(2,2) * & + (one - mix_input_output_XC_GGA_stress) + XC_stress(3,3) = XC_stress(3,3) + spin_factor*XC_GGA_stress(3,3) * & + (one - mix_input_output_XC_GGA_stress) + end if jacobian = zero do spin = 1, nspin jacobian = jacobian + spin_factor*dot(nsize,density_out(:,spin),1,potential(:,spin),1) @@ -4093,6 +4117,7 @@ contains print_timer, stop_print_timer, & WITH_LEVEL, TIME_ACCUMULATE_NO, & TIME_ACCUMULATE_YES + use XC, ONLY: XC_GGA_stress implicit none @@ -4145,6 +4170,26 @@ contains dcellz_grid = dcellz_block / nz_in_block call start_timer (tmr_l_tmp2) + ! Find XC_GGA_stress using output density for non-SCF case + ! we will average the input and output as an approximation + ! to the correct (hideously complex) stress so add mixing factor + if((.not.flag_self_consistent).and.(.not.flag_mix_L_SC_min)) then + if(.NOT.present(density_out)) call cq_abort("Output density not passed to PCC force for nonSCF calculation") + if (flag_stress) then + density_wk = zero + do spin = 1, nspin + density_wk(:,spin) = density_out(:,spin) + half * density_pcc(:) + end do + call get_xc_potential(density_wk, xc_potential, & + xc_epsilon, xc_energy, size) + XC_stress(1,1) = XC_stress(1,1) + spin_factor*XC_GGA_stress(1,1) * & + mix_input_output_XC_GGA_stress + XC_stress(2,2) = XC_stress(2,2) + spin_factor*XC_GGA_stress(2,2) * & + mix_input_output_XC_GGA_stress + XC_stress(3,3) = XC_stress(3,3) + spin_factor*XC_GGA_stress(3,3) * & + mix_input_output_XC_GGA_stress + end if + end if density_wk = zero do spin = 1, nspin density_wk(:,spin) = density(:,spin) + half * density_pcc(:) @@ -4152,10 +4197,20 @@ contains call get_xc_potential(density_wk, xc_potential, & xc_epsilon, xc_energy, size) + ! And now add XC_GGA_stress for density_in scaled by half + if((.not.flag_self_consistent).and.(.not.flag_mix_L_SC_min)) then + if (flag_stress) then + XC_stress(1,1) = XC_stress(1,1) + spin_factor*XC_GGA_stress(1,1) * & + (one - mix_input_output_XC_GGA_stress) + XC_stress(2,2) = XC_stress(2,2) + spin_factor*XC_GGA_stress(2,2) * & + (one - mix_input_output_XC_GGA_stress) + XC_stress(3,3) = XC_stress(3,3) + spin_factor*XC_GGA_stress(3,3) * & + (one - mix_input_output_XC_GGA_stress) + end if + end if if(PRESENT(xc_energy_ret)) xc_energy_ret = xc_energy ! We do this here to re-use xc_potential - for non-PCC we do it in get_nonSC_correction_force if((.not.flag_self_consistent).and.(.not.flag_mix_L_SC_min)) then - if(.NOT.present(density_out)) call cq_abort("Output density not passed to PCC force for nonSCF calculation") if (flag_stress) then jacobian = zero do spin=1,nspin diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index d0f6405e..5bd7ef47 100644 --- a/src/initial_read_module.f90 +++ b/src/initial_read_module.f90 @@ -941,6 +941,7 @@ contains use biblio, only: flag_dump_bib !2019/12/27 tsuyoshi use density_module, only: method_UpdateChargeDensity,DensityMatrix,AtomicCharge + use force_module, only: mix_input_output_XC_GGA_stress implicit none @@ -1572,6 +1573,10 @@ contains sqnm_trust_step = fdf_double ('AtomMove.MaxSQNMStep',0.2_double ) LBFGS_history = fdf_integer('AtomMove.LBFGSHistory', 5 ) flag_opt_cell = fdf_boolean('AtomMove.OptCell', .false.) + ! At present (2023/07/26 just before v1.2 release) neutral atom is required for cell opt + if(flag_opt_cell.and.(.not.flag_neutral_atom)) & + call cq_abort("You must use neutral atom for cell optimisation") + ! This can be removed when ewald update is implemented flag_variable_cell = flag_opt_cell optcell_method = fdf_integer('AtomMove.OptCellMethod', 1) cell_constraint_flag = fdf_string(20,'AtomMove.OptCell.Constraint','none') @@ -1583,6 +1588,7 @@ contains flag_stress = fdf_boolean('AtomMove.CalcStress', .true.) flag_full_stress = fdf_boolean('AtomMove.FullStress', .false.) flag_atomic_stress = fdf_boolean('AtomMove.AtomicStress', .false.) + mix_input_output_XC_GGA_stress = fdf_double('General.MixXCGGAInOut',half) ! flag_vary_basis = fdf_boolean('minE.VaryBasis', .false.) if(.NOT.flag_vary_basis) then