phonopy and phono3py find different number of phonon displacements

Check the develop branch
phonopy 210509ab99, phono3py 19fc8c7, both should be latest develop branch

Describe the bug
different numbers of phonon displacement configurations are found by phonopy and phono3py

To Reproduce
run

python << EOF
import phonopy
import phono3py

cell = [[5.9479774024925023,    0.0000000000000000,    0.0000000000000000],
        [0.0000000000000000,    5.9479774024925023,    0.0000000000000000],
        [0.0000000000000000,    0.0000000000000000,    5.9479774024925023]]

nums = [55,82,35,35,35]
pos = [[2.9739887012462511,  2.9739887012462511,  2.9739887012462511],
       [0.0000000000000000,  0.0000000000000000,  0.0000000000000000],
       [2.9739887012462511,  0.0000000000000000,  0.0000000000000000],
       [0.0000000000000000,  2.9739887012462511,  0.0000000000000000],
       [0.0000000000000000,  0.0000000000000000,  2.9739887012462511]]

at_ph = phonopy.structure.atoms.PhonopyAtoms(numbers=nums, positions=pos, cell=cell)

ph2 = phonopy.Phonopy(at_ph, supercell_matrix=[2,2,2])
ph2.generate_displacements(distance=0.1)
print("number of phonon displacements found by phonopy", len(ph2.supercells_with_displacements))

ph3 = phono3py.Phono3py(at_ph, supercell_matrix=[2,2,2], phonon_supercell_matrix=[2,2,2])
ph3.generate_displacements(distance=0.1, cutoff_pair_distance=4.0)

print("number of phonon displacements found by phono3py", len(ph3.phonon_supercells_with_displacements))
EOF

I get

number of phonon displacements found by phonopy 3
number of phonon displacements found by phono3py 4

and I expected to get the same number from both codes.

Please have a look at the docstrings:

phono3py/phono3py/api_phono3py.py

Lines 1480 to 1507 in 19fc8c7

    
               def generate_displacements( 
        
                   self, 
        
                   distance=0.03, 
        
                   cutoff_pair_distance=None, 
        
                   is_plusminus="auto", 
        
                   is_diagonal=True, 
        
               ): 
        
                   """Generate displacement dataset in supercell for fc3. 
        
                   This systematically generates single and pair atomic displacements 
        
                   in supercells to calculate fc3 considering crystal symmetry. 
        
                   For fc3, two atoms are displaced for each configuration 
        
                   considering crystal symmetry. The first displacement is chosen 
        
                   in the perfect supercell, and the second displacement in the 
        
                   displaced supercell. The first displacements are taken along 
        
                   the basis vectors of the supercell. This is because the 
        
                   symmetry is expected to be less broken by the introduced first 
        
                   displacement, and as the result, the number of second 
        
                   displacements may become smaller than the case that the first 
        
                   atom is displaced not along the basis vectors. 
        
                   Note 
        
                   ---- 
        
                   When phonon_supercell_matrix is not given, fc2 is also 
        
                   computed from the same set of the displacements for fc3 and 
        
                   respective supercell forces. When phonon_supercell_matrix is 
        
                   set, the displacements in phonon_supercell are generated.

phono3py/phono3py/api_phono3py.py

Lines 1548 to 1561 in 19fc8c7

    
               def generate_fc2_displacements( 
        
                   self, distance=0.03, is_plusminus="auto", is_diagonal=False 
        
               ): 
        
                   """Generate displacement dataset in phonon supercell for fc2. 
        
                   This systematically generates single atomic displacements 
        
                   in supercells to calculate phonon_fc2 considering crystal symmetry. 
        
                   Note 
        
                   ---- 
        
                   is_diagonal=False is chosen as the default setting intentionally 
        
                   to be consistent to the first displacements of the fc3 pair 
        
                   displacemets in supercell.

I'm afraid that I don't see anything in the docstrings that explains the difference between the values returned by Phonopy and Phono3py for the fc2 (phonon) displacements. Is there a specific item that you think is particularly relevant?

I did just check what happens when I call Phono3py.generate_fc2_displacements(), with and without is_diagonal=True, and in every case it also returns 4 displacements, while Phonopy.generate_displacements() returns only 3.

Maybe I see your issue.

I didn't consider rerun generate_displacements (or generate_fc2_displacements). So with the same Phono3py instance, supercells_with_displacements (or phonon_supercells_with_displacements) is not updated. Please see the following scripts and runs. Is this your issue?

Case 1

import phonopy
import phono3py

cell = [
    [5.9479774024925023, 0.0000000000000000, 0.0000000000000000],
    [0.0000000000000000, 5.9479774024925023, 0.0000000000000000],
    [0.0000000000000000, 0.0000000000000000, 5.9479774024925023],
]

nums = [55, 82, 35, 35, 35]
pos = [
    [2.9739887012462511, 2.9739887012462511, 2.9739887012462511],
    [0.0000000000000000, 0.0000000000000000, 0.0000000000000000],
    [2.9739887012462511, 0.0000000000000000, 0.0000000000000000],
    [0.0000000000000000, 2.9739887012462511, 0.0000000000000000],
    [0.0000000000000000, 0.0000000000000000, 2.9739887012462511],
]

at_ph = phonopy.structure.atoms.PhonopyAtoms(numbers=nums, positions=pos, cell=cell)

print("is_diagonal=True:")
ph3 = phono3py.Phono3py(
    at_ph, supercell_matrix=[2, 2, 2], phonon_supercell_matrix=[2, 2, 2]
)
ph3.generate_displacements(distance=0.1, is_diagonal=True)
print(
    "number of displacements found by phono3py",
    len(ph3.supercells_with_displacements),
)
ph3.generate_fc2_displacements(distance=0.1, is_diagonal=True)
print(
    "number of phonon displacements found by phono3py",
    len(ph3.phonon_supercells_with_displacements),
)

print("is_diagonal=False:")
ph3 = phono3py.Phono3py(
    at_ph, supercell_matrix=[2, 2, 2], phonon_supercell_matrix=[2, 2, 2]
)
ph3.generate_displacements(distance=0.1, is_diagonal=False)
print(
    "number of displacements found by phono3py",
    len(ph3.supercells_with_displacements),
)
ph3.generate_fc2_displacements(distance=0.1, is_diagonal=False)
print(
    "number of phonon displacements found by phono3py",
    len(ph3.phonon_supercells_with_displacements),
)

% python generate-disps.py
is_diagonal=True:
number of displacements found by phono3py 220
number of phonon displacements found by phono3py 3
is_diagonal=False:
number of displacements found by phono3py 300
number of phonon displacements found by phono3py 4

Case 2

import phonopy
import phono3py

cell = [
    [5.9479774024925023, 0.0000000000000000, 0.0000000000000000],
    [0.0000000000000000, 5.9479774024925023, 0.0000000000000000],
    [0.0000000000000000, 0.0000000000000000, 5.9479774024925023],
]

nums = [55, 82, 35, 35, 35]
pos = [
    [2.9739887012462511, 2.9739887012462511, 2.9739887012462511],
    [0.0000000000000000, 0.0000000000000000, 0.0000000000000000],
    [2.9739887012462511, 0.0000000000000000, 0.0000000000000000],
    [0.0000000000000000, 2.9739887012462511, 0.0000000000000000],
    [0.0000000000000000, 0.0000000000000000, 2.9739887012462511],
]

at_ph = phonopy.structure.atoms.PhonopyAtoms(numbers=nums, positions=pos, cell=cell)

print("is_diagonal=True:")
ph3 = phono3py.Phono3py(
    at_ph, supercell_matrix=[2, 2, 2], phonon_supercell_matrix=[2, 2, 2]
)
ph3.generate_displacements(distance=0.1, is_diagonal=True)
print(
    "number of displacements found by phono3py",
    len(ph3.supercells_with_displacements),
)
ph3.generate_fc2_displacements(distance=0.1, is_diagonal=True)
print(
    "number of phonon displacements found by phono3py",
    len(ph3.phonon_supercells_with_displacements),
)

print("is_diagonal=False:")
ph3.generate_displacements(distance=0.1, is_diagonal=False)
print(
    "number of displacements found by phono3py",
    len(ph3.supercells_with_displacements),
)
ph3.generate_fc2_displacements(distance=0.1, is_diagonal=False)
print(
    "number of phonon displacements found by phono3py",
    len(ph3.phonon_supercells_with_displacements),
)

% python generate-disps.py
is_diagonal=True:
number of displacements found by phono3py 220
number of phonon displacements found by phono3py 3
is_diagonal=False:
number of displacements found by phono3py 220
number of phonon displacements found by phono3py 3

Thank you - it does appear that the sequence of operations is indeed the problem. However, I find the behavior even more surprising than your examples show. What looks like a simple access of the attribute ph3.phonon_supercells_with_displacements, after the call to generate_displacements(), makes it so the second call, to generate_fc2_displacements(), does nothing. However, if I just call generate_displacements(), access onlysupercells_with_displacements, and then call generate_fc2_displacements(), then I get the expect number of phonon displacements. The script below shows this behavior.

I think it would be helpful if the documentation was explicit about the allowed calling order, or even better, if the code detected that something has been run before and gave an error, rather than giving inconsistent numbers of displacement configurations depending on the order of calls.

import phonopy
import phono3py

cell = [
    [5.9479774024925023, 0.0000000000000000, 0.0000000000000000],
    [0.0000000000000000, 5.9479774024925023, 0.0000000000000000],
    [0.0000000000000000, 0.0000000000000000, 5.9479774024925023],
]

nums = [55, 82, 35, 35, 35]
pos = [
    [2.9739887012462511, 2.9739887012462511, 2.9739887012462511],
    [0.0000000000000000, 0.0000000000000000, 0.0000000000000000],
    [2.9739887012462511, 0.0000000000000000, 0.0000000000000000],
    [0.0000000000000000, 2.9739887012462511, 0.0000000000000000],
    [0.0000000000000000, 0.0000000000000000, 2.9739887012462511],
]

at_ph = phonopy.structure.atoms.PhonopyAtoms(numbers=nums, positions=pos, cell=cell)

####################################################################################################
ph3 = phono3py.Phono3py(
    at_ph, supercell_matrix=[2, 2, 2], phonon_supercell_matrix=[2, 2, 2]
)
ph3.generate_displacements(distance=0.1, cutoff_pair_distance=4.0, is_diagonal=True)
l1 = len(ph3.supercells_with_displacements)
l2 = len(ph3.phonon_supercells_with_displacements)
ph3.generate_fc2_displacements(distance=0.1, is_diagonal=True)
print("number of phonon displacements from")
print("    generate_displacements()")
print("    len(supercells_With_displacements")
print("    len(phonon_supercells_with_displacements)")
print("    generate_fc2_displacements()")
print("    len(phonon_supercells_with_displacements)")
print(len(ph3.phonon_supercells_with_displacements))

print("")
####################################################################################################
ph3 = phono3py.Phono3py(
    at_ph, supercell_matrix=[2, 2, 2], phonon_supercell_matrix=[2, 2, 2]
)
ph3.generate_displacements(distance=0.1, cutoff_pair_distance=4.0, is_diagonal=True)
l1 = len(ph3.supercells_with_displacements)
# l2 = len(ph3.phonon_supercells_with_displacements)
ph3.generate_fc2_displacements(distance=0.1, is_diagonal=True)
print("number of phonon displacements from")
print("    generate_displacements()")
print("    len(supercells_with_displacements)")
print("    generate_fc2_displacements()")
print("    len(phonon_supercells_with_displacements")
print(len(ph3.phonon_supercells_with_displacements))

Output:

number of phonon displacements from
    generate_displacements()
    len(supercells_With_displacements
    len(phonon_supercells_with_displacements)
    generate_fc2_displacements()
    len(phonon_supercells_with_displacements)
4

number of phonon displacements from
    generate_displacements()
    len(supercells_with_displacements)
    generate_fc2_displacements()
    len(phonon_supercells_with_displacements
3

Thanks your inputs. I fix it at PR #65. To run this, you need the develop branch of phonopy. After this fix, your last script results in

number of phonon displacements from
    generate_displacements()
    len(supercells_With_displacements
    len(phonon_supercells_with_displacements)
    generate_fc2_displacements()
    len(phonon_supercells_with_displacements)
3

number of phonon displacements from
    generate_displacements()
    len(supercells_with_displacements)
    generate_fc2_displacements()
    len(phonon_supercells_with_displacements
3

I close this issue because I think the issue is fixed. Please reopen if you need.

	def generate_displacements(
	self,
	distance=0.03,
	cutoff_pair_distance=None,
	is_plusminus="auto",
	is_diagonal=True,
	):
	"""Generate displacement dataset in supercell for fc3.

	This systematically generates single and pair atomic displacements
	in supercells to calculate fc3 considering crystal symmetry.

	For fc3, two atoms are displaced for each configuration
	considering crystal symmetry. The first displacement is chosen
	in the perfect supercell, and the second displacement in the
	displaced supercell. The first displacements are taken along
	the basis vectors of the supercell. This is because the
	symmetry is expected to be less broken by the introduced first
	displacement, and as the result, the number of second
	displacements may become smaller than the case that the first
	atom is displaced not along the basis vectors.

	Note
	----
	When phonon_supercell_matrix is not given, fc2 is also
	computed from the same set of the displacements for fc3 and
	respective supercell forces. When phonon_supercell_matrix is
	set, the displacements in phonon_supercell are generated.

	def generate_fc2_displacements(
	self, distance=0.03, is_plusminus="auto", is_diagonal=False
	):
	"""Generate displacement dataset in phonon supercell for fc2.

	This systematically generates single atomic displacements
	in supercells to calculate phonon_fc2 considering crystal symmetry.


	Note
	----
	is_diagonal=False is chosen as the default setting intentionally
	to be consistent to the first displacements of the fc3 pair
	displacemets in supercell.