[QE-users] Details of Automatic FFT Grid Size Calculation
Dyer, Brock
brdyer at ursinus.edu
Thu Mar 6 01:19:30 CET 2025
I believe I've been able to implement a method that can determine the dimensions of the FFT grid for most cases (I neglected to include the check for if the system is a NEC SX-6 supercomputer since they are more than 20 years old at this point).
Huge thanks to Stefano de Gironcoli for helping point me in the right direction!
Also, if I have made a mistake here, I would really appreciate it if someone could let me know so I can rectify my error!
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
from math import gcd, sqrt, ceil, log2
import numpy.typing as npt
import numpy as np
def allowed_fft_dimension(nr: int, is_essl: bool = False) -> bool:
"""
Determine if the FFT dimension is allowable
### Parameters
1. nr : int
- Possible FFT dimension to test
2. is_essl : bool, (default False)
- True only if QE was compiled with the Engineering Scientific Subroutine Library (ESSL)
"""
if is_essl == True and log2(nr).is_integer():
return False
factors = [2, 3, 5, 7, 11]
pwr = [0, 0, 0, 0, 0]
mr = nr
count = 0
while mr > 1 and count < 5:
for j in range(5):
if mr % factors[j] == 0:
count = 0
pwr[j] += 1
mr = mr / factors[j]
else:
count += 1
continue
if mr != 1:
False
elif (
is_essl == True
and not pwr[0] >= 1
and not pwr[1] <= 2
and not pwr[2] <= 1
and not pwr[3] <= 1
and not pwr[4] <= 1
and not (pwr[1] == 0 and (pwr[2] + pwr[3] + pwr[4]) <= 2)
or not ((pwr[1] != 0 and pwr[2] + pwr[3] + pwr[4]) <= 1)
):
return False
elif (
(pwr[3] or pwr[4]) != 0 or count >= 5
):
return False
else:
return True
def get_fft_dimensions(
v1: npt.ArrayLike,
v2: npt.ArrayLike,
v3: npt.ArrayLike,
ecutwfc: int,
ecutrho: int = None,
is_essl: bool = False,
) -> list[int]:
"""
Determine FFT dimension from a given unit cell
### Parameters
1. v1 : ArrayLike
- First dimension of unit cell
2. v2 : ArrayLike
- Second dimension of unit cell
3. v3 : ArrayLike
- Third dimension of unit cell
4. ecutwfc : int
- Kinetic energy cutoff (in Rydbergs)
5. ecutrho : int, (default None)
- Kinetic energy cutoff for charge density and potential (in Rydbergs)
- If not given, `ecutrho = ecutwfc * 4`
6. is_essl : bool, (default False)
- Set True if QE was compiled with the Engineering Scientific Subroutine Library (ESSL)
"""
if ecutrho is None:
ecutrho = 4 * ecutwfc
# Lattice Parameter (corner to corner distance in bohrs)
alat = np.sqrt(v1[0] ** 2 + v1[1] ** 2 + v1[2] ** 2)
# Cell parameters but in alat units
at: list[list[float]] = [
v1 / alat,
v2 / alat,
v3 / alat,
]
tpiba = (2.0 * np.pi) / alat
gcutm = ecutrho / (tpiba**2)
# Initial guesses at FFT grid sizes
nr1 = 2 * int(np.sqrt(gcutm) * np.sqrt(at[0][0] ** 2 + at[0][1] ** 2 + at[0][2] ** 2)) + 1
nr2 = 2 * int(np.sqrt(gcutm) * np.sqrt(at[1][0] ** 2 + at[1][1] ** 2 + at[1][2] ** 2)) + 1
nr3 = 2 * int(np.sqrt(gcutm) * np.sqrt(at[2][0] ** 2 + at[2][1] ** 2 + at[2][2] ** 2)) + 1
nrx = [nr1, nr2, nr3]
for i in range(3):
while not allowed_fft_dimension(nrx[i], is_essl):
nrx[i] += 1
if nrx[i] > 16385:
raise ValueError("FFT Dimension exceeds maximum allowed dimension (16385)")
return nrx
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
________________________________
From: users on behalf of Stefano de Gironcoli
Sent: Sunday, March 2, 2025 5:18 AM
To: users at lists.quantum-espresso.org
Subject: Re: [QE-users] Details of Automatic FFT Grid Size Calculation
Dear Dyer,
you are right... the factor of two you find is the difference between the radius and the diameter of a sphere.
the bit of code in FFTXlib/src/fft_type.f90 I pointed to earlier defines the ranges nr1,nr2,nr3 of a parallelepiped box that contains the G-sphere of sqrt(gcutm) radius.
these are then passed to grid_set( dfft, bg, gcutm, dfft%nr1, dfft%nr2, dfft%nr3 ) where the actual size is determined in three nested loops that range from - the estimated limit to + the estimated limit
DO k = -nr3, nr3
DO j = -nr2, nr2
DO i = -nr1, nr1
where the maximum ACTUAL absolute index values are set as nb(1),nb(2),nb(3), which should coincide with nr1,nr2,nr3 if they were computed properly, and finally set the grid dimensions as
nr1 = 2 * nb(1) + 1
nr2 = 2 * nb(2) + 1
nr3 = 2 * nb(3) + 1
stefano
On 02/03/25 06:04, Dyer, Brock wrote:
I've been tracing down all the variables that are required to generate a good FFT grid, and I seem to be off by a factor of 2 in the end. I can't quite figure out what the issue is, but my guess may be that while I am storing 'ecutwfc' in Rydbergs there is a conversion to Hartrees somewhere in the code that I haven't seen yet.
My current process looks like this (with some values from a run I had recently so I can compare):
ecutwfc = 100 Ry
Ecutrho = 400 Ry
# Unit cell dimensions, given in Bohrs
v1 = [44.09733757, 0.0, 0.0]
v2 = [0.0, 44.09733757, 0.0]
v3 = [0.0, 0.0, 44.09733757]
alat = sqrt(v1[0]**2 + v1[1]**2 + v1[2]**2)
tpiba = (2.0 * pi) / alat
gcutm = ecutrho / (tpiba**2)
at = [v1/alat, v2/alat, v3/alat]
nr1 = int(sqrt(gcutm) * sqrt(at[0][0]**2 + at[0][1]**2 + at[0][2]**2)) + 1
nr2 = int(sqrt(gcutm) * sqrt(at[1][0]**2 + at[1][1]**2 + at[1][2]**2)) + 1
nr3 = int(sqrt(gcutm) * sqrt(at[2][0]**2 + at[2][1]**2 + at[2][2]**2)) + 1
These last lines are where I've noticed the problem. From looking at the output of my run with the given cell sizes, I expect an FFT grid of 288x288x288, however if I were to run this code (and the code that checks if it's a good size) I'd get an FFT grid that is only half that. I'd love some advice on this if it is at all possible. I also can send some more formatted code if it would help (I decided to cut down the python code so it looked a bit more like the original f90 code).
________________________________
From: users on behalf of Stefano de Gironcoli
Sent: Thursday, February 27, 2025 2:50 PM
To: users at lists.quantum-espresso.org<mailto:users at lists.quantum-espresso.org>
Subject: Re: [QE-users] Details of Automatic FFT Grid Size Calculation
it's in SUBROUTINE realspace_grid_init in FFTXlib/src/file fft_types.f90
!
! ... calculate the size of the real-space dense grid for FFT
! ... first, an estimate of nr1,nr2,nr3, based on the max values
! ... of n_i indices in: G = i*b_1 + j*b_2 + k*b_3
! ... We use G*a_i = n_i => n_i .le. |Gmax||a_i|
!
dfft%nr1 = int ( sqrt (gcutm) * sqrt (at(1, 1)**2 + at(2, 1)**2 + at(3, 1)**2) ) + 1
dfft%nr2 = int ( sqrt (gcutm) * sqrt (at(1, 2)**2 + at(2, 2)**2 + at(3, 2)**2) ) + 1
dfft%nr3 = int ( sqrt (gcutm) * sqrt (at(1, 3)**2 + at(2, 3)**2 + at(3, 3)**2) ) + 1
stefano
On 27/02/25 20:11, Dyer, Brock wrote:
Hello all, I have been working quite a bit lately on automating my QE workflow, and as part of that workflow I check the automatically calculated FFT grid sizes for the level of theory that I have been using in order to improve my parallelization.
I have tried tracing down and reading the code that calculates the FFT grid sizes, however I cannot find/understand the actual code to calculate the grid sizes. My current understanding is that the initial parameters to calculate the grid size are 'ecutwfc' and/or 'ecutrho', and the unit cell size, and then there seems to be some more math, and perhaps at the end the final dimensions get calculated in 'fft_ggen.f90'.
What I am looking for ideally is a mathematical formula that includes all of the input parameters and operations required to calculate the FFT grid sizes so that I can implement it into my workflow and hopefully not have to run double calculations to properly parallelize.
Thanks in advance for the help!
Brock Dyer, Ursinus College Class of 2025
_______________________________________________________________________________
The Quantum ESPRESSO Foundation stands in solidarity with all civilians worldwide who are victims of terrorism, military aggression, and indiscriminate warfare.
--------------------------------------------------------------------------------
Quantum ESPRESSO is supported by MaX (www.max-centre.eu<https://linkprotect.cudasvc.com/url?a=http%3a%2f%2fwww.max-centre.eu&c=E,1,o6mRyDmIfQXWXGr_JikolRfWqjEMwbk8xt5Q6J2l_N7Fky0gnXgD_aQU9TdYIxh51SmHmmgN9S5AJxHAS6vpqDSxTijqfOUXWuQTbchy3TZnaqcHjGvJ&typo=1>)
users mailing list users at lists.quantum-espresso.org<mailto:users at lists.quantum-espresso.org>
https://lists.quantum-espresso.org/mailman/listinfo/users<https://linkprotect.cudasvc.com/url?a=https%3a%2f%2flists.quantum-espresso.org%2fmailman%2flistinfo%2fusers&c=E,1,MVtAKpGoN3VwtTvpYzhgehNBw3lc6Ccvz18zJQYi8KwO2z3w1lLip_vRNXU25kTl-Lr0eMQKhqjTUVWShMvELzUclPM0hjrp4RgEFfOWL7-pZnGUAGdhpqG_ZBdI&typo=1>
_______________________________________________________________________________
The Quantum ESPRESSO Foundation stands in solidarity with all civilians worldwide who are victims of terrorism, military aggression, and indiscriminate warfare.
--------------------------------------------------------------------------------
Quantum ESPRESSO is supported by MaX (www.max-centre.eu<https://linkprotect.cudasvc.com/url?a=http%3a%2f%2fwww.max-centre.eu&c=E,1,Fk3xckU46d4kJt2B1k4dAx45Mwql_NLS2Q4mFkaiKohZ4Mfa7G89x8TkM4SYAanmfYlCxa06qxMoQUAjwQA91gR92qp0jD9P8YocmrWMlVNnTg,,&typo=1>)
users mailing list users at lists.quantum-espresso.org<mailto:users at lists.quantum-espresso.org>
https://lists.quantum-espresso.org/mailman/listinfo/users<https://linkprotect.cudasvc.com/url?a=https%3a%2f%2flists.quantum-espresso.org%2fmailman%2flistinfo%2fusers&c=E,1,WMy7IPJMHV7ccyILg4tXoP0iKMNXMCWvyCTs45Pq2v6EoS2BX6kzobvtwlk9bVU73K1T2PctdQG4gNDHs7YrqhJAP06L3Iec_QxBdAJhkcf8I_sYUBgzm9x9&typo=1>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20250306/7dd11864/attachment.html>
More information about the users
mailing list