Skip to content
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

Jump to bottom

Implement stresses without unfold_bz #511

Merged
merged 9 commits into from
Jul 30, 2021
Merged

Implement stresses without unfold_bz #511

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
03f6c31
Stresses without unfold_bz
jaemolihm Jul 30, 2021
d623830
Increase kgrid in stresses test
jaemolihm Jul 30, 2021
c23ddd8
Enable MPI for stresses test
jaemolihm Jul 30, 2021
d7180d4
Rename symmetrize_tensor to symmetrize_stresses
jaemolihm Jul 30, 2021
39aa683
review
jaemolihm Jul 30, 2021
d5460fb
runtests organization
jaemolihm Jul 30, 2021
3324ccb
Add doc on symmetrization
jaemolihm Jul 30, 2021
1ccda15
Fix stresses test for MPI
jaemolihm Jul 30, 2021
682ba70
Merge branch 'master' into stresses_sym
mfherbst Jul 30, 2021
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
4 changes: 2 additions & 2 deletions src/postprocess/stresses.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@ Compute the stresses (= 1/Vol dE/d(M*lattice), taken at M=I) of an obtained SCF
"""
@timing function compute_stresses(scfres)
# TODO optimize by only computing derivatives wrt 6 independent parameters
scfres = unfold_bz(scfres)
# compute the Hellmann-Feynman energy (with fixed ψ/occ/ρ)
function HF_energy(lattice)
T = eltype(lattice)
Expand All @@ -31,5 +30,6 @@ Compute the stresses (= 1/Vol dE/d(M*lattice), taken at M=I) of an obtained SCF
end
L = scfres.basis.model.lattice
Ω = scfres.basis.model.unit_cell_volume
ForwardDiff.gradient(M -> HF_energy((I+M) * L), zero(L)) / Ω
stresses_not_symmetrized = ForwardDiff.gradient(M -> HF_energy((I+M) * L), zero(L)) / Ω
symmetrize_tensor(stresses_not_symmetrized, scfres.basis.symmetries, L)
end
11 changes: 11 additions & 0 deletions src/symmetry.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -314,3 +314,14 @@ function unfold_bz(scfres)
new_scfres = (; basis=basis_unfolded, ψ, ham, eigenvalues, occupation)
merge(scfres, new_scfres)
end

# symmetrize a rank-2 tensor in reduced coordinates
function symmetrize_tensor(tensor, symmetries, lattice)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this file we generally follow the convention (system arguments, data arguments), i.e. it should be lattice, symmetries, tensor.

Also can you add ::Matrix or even better ::Mat3 to indicate this only works for matrices instead of general tensors.
In fact ::Mat3 would even allow to rename this function to symmetrize (same as the symmetrize we use for the densities) and add another method symmetrize(model, tensor::Mat3) and/or symmetrize(model, tensor::Mat3).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not too sure about this, for two reasons. The first is that this routine is specific to stresses-like tensors, ie that transform in a specific way, being the derivative of the energy wrt the lattice, which transforms in another way (I would imagine the term is that lattices are covariant, while stresses are contravariant? I'm hesitant here, this is definitely not my zone of confort). The second is that is that it's also specific to the fact that our stresses are in reduced coordinates. So I would call it symmetrize_stresses.

Unless it's canonical that "tensors" transform in this particular way? In which call calling it symmetrize_tensor is OK by me, since things are in reduced coordinates by default (we use _cart for non-reduced)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To clarify with "about this" you mean the renaming to symmetrize, but not the rest.

Yeah ok, I might have been a bit quick here. From the way stresses are formed I agree with your argument that they could be special (likely the only contravariant thing we deal with in the foreseeable time). I also never looked into this in detail, however.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To clarify with "about this" you mean the renaming to symmetrize, but not the rest.

I agree this takes Mat3 as inputs, but I'm not sure we should aim for a general symmetrize function that deduces what it should do by dispatch, since we generally don't use types a lot (eg a density is just an array). And if we keep explicit names there's no reason to type arguments. I do agree about the order though

Copy link
Member

@mfherbst mfherbst Jul 30, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

then we should rename to symmetrize_ρ or symmetrize_density?

Edit: I think we had this initially.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm, real is a bit generic... I think we already use "density" quite generically to denote "density or potential" at some other places in the code. I think it's better to stick with "density" even though it's more general, unless we can come up with another name to mean "real-space quantity"?

Copy link
Member

@mfherbst mfherbst Jul 30, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"symmetrize_densitylike" or "symmetrize_ρlike" ??

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

at that point it's simpler to call it density, no ? :-p

Copy link
Member

@mfherbst mfherbst Jul 30, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok. symmetrize_density it is.

@jaemolihm Sorry about that bikeshedding. If you want do the renaming, otherwise I'll do it in a followup.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll leave it as your followup. 😄

tensor_symmetrized = zero(tensor)
for (S, τ) in symmetries
S_reduced = inv(lattice) * S * lattice
tensor_symmetrized += S_reduced' * tensor * S_reduced
end
tensor_symmetrized /= length(symmetries)
tensor_symmetrized
end
2 changes: 1 addition & 1 deletion test/stresses.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ include("testcases.jl")
Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
atoms = [Si => silicon.positions]
model = model_DFT(lattice, atoms, [:lda_x, :lda_c_vwn]; symmetries)
kgrid = [2, 2, 1]
kgrid = [3, 3, 1]
mfherbst marked this conversation as resolved.
Show resolved Hide resolved
Ecut = 7
PlaneWaveBasis(model; Ecut, kgrid)
end
Expand Down