Significant slowdown upon change of the wavelength

[reddocksw]

Hi,

I have been trying to obtain X-ray patterns from a molecular dynamics simulation: I want to generate the diffraction pattern for every snapshot obtained from the simulation, average the intensities and create one plot. Everything has been working pretty smoothly for the default wavelength, however, when I switched to MoKa1 (based on the experiment), the time needed for generating the pattern for a single snapshot increased from 9 seconds up to roughly 300 seconds (with ~240 seconds needed for pm.simulate), and I have to iterate over 10 thousands of the snapshots. I wonder if that behavior is expected or if there is something bad happening in the background. The structures are read from CIF files.

def main():
    cell = cell_reader(files)
    intensities = []
    X = np.arange(0, 45, 0.015)
    for i in range(1001, 11000 + 1):
        mol = read('ZIF_8_Zn.xyz', index=i)
        mol.cell = cell[i]
        mol.pbc = True
        mol.write(f'temp/temp_{i}.cif')
        MOF = xu.materials.Crystal.fromCIF(os.path.join("temp", f"temp_{i}.cif"))
        MOF_powder = xu.simpack.Powder(MOF, 1)
        pm = xu.simpack.PowderModel(MOF_powder, wl='MoKa1')
        Y = pm.simulate(X)
        sf = max(Y) / 100
        intensities.append(Y / sf)

    intensities = np.array(intensities)
    intensities_avg = intensities.mean(axis=0)
    plt.plot(X, intensities_avg)
    plt.savefig(f'test_New.png', bbox_inches='tight', dpi=300)

if __name__ == '__main__':
    multiprocessing.freeze_support()
    main()

[dkriegner]

I guess some increase is expected since with MoKa1 you get more peaks. But what you describe sounds too much.

Can you provide a full working example including the needed input files? If confidential feel free to drop me an email.

One thing I see right away is that you should limit the angular range for the PowderModel. I believe PowderModel supports a tt_cutoff parameter which it forwards to PowerDiffraction.

Another thing I notice is that you should try to reuse the PowderModel object. There is a lot of caching going on which could help massively reduce the execution time. Currently I am not sure how you can load the new structure. Am I right that the CIF file is always in spacegroup P1? Is the number of atoms always the same?

On Mon 29. 7. 2024 at 0:05, reddocksw @.***> wrote:

Hi,

I have been trying to obtain XRD patterns from a molecular dynamics simulation: I want to generate the diffraction pattern for every snapshot obtained from the simulation, average the intensities and create one plot. Everything has been working pretty smoothly for the default wavelength, however, when I switched to MoKa1 (based on the experiment), the time needed for generating the pattern for a single snapshot increased from 9 seconds up to roughly 300 seconds (with ~240 seconds needed for pm.simulate), and I have to iterate over 10 thousands of the snapshots. I wonder if that behavior is expected or if there is something bad happening in the background. The structures are read from CIF files.

def main(): cell = cell_reader(files) intensities = [] X = np.arange(0, 45, 0.015) for i in range(1001, 11000 + 1): mol = read('ZIF_8Zn.xyz', index=i) mol.cell = cell[i] mol.pbc = True mol.write(f'temp/temp{i}.cif') MOF = xu.materials.Crystal.fromCIF(os.path.join("temp", f"temp_{i}.cif")) MOF_powder = xu.simpack.Powder(MOF, 1) pm = xu.simpack.PowderModel(MOF_powder, wl='MoKa1') Y = pm.simulate(X) sf = max(Y) / 100 intensities.append(Y / sf)

intensities = np.array(intensities) intensities_avg = intensities.mean(axis=0) plt.plot(X, intensities_avg) plt.savefig(f'test_New.png', bbox_inches='tight', dpi=300)

if name == 'main': multiprocessing.freeze_support() main()``

— Reply to this email directly, view it on GitHub https://github.com/dkriegner/xrayutilities/issues/191, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACKZJFMSOPH6EZQO3J5ZX73ZOVTL3AVCNFSM6AAAAABLTEUQCKVHI2DSMVQWIX3LMV43ASLTON2WKOZSGQZTIMJTGMZDAMQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[reddocksw]

Hello,

yes, the CIF is always in the P1 spacegroup. The number of atoms is constant, however, every snapshot is unique, as the atomic positions and cell vectors can vary throughout the simulation. Because of that, I am not sure if the PowderModel can be reused, as the system is slightly different every time.

Using the tt_cutoff has efficiently reduced the time from 300 seconds to 30 seconds, but sometimes it hits 120 seconds or even more (for 3 of the enclosed files, the pm.simulate() took 25-30 seconds but 234 seconds for temp_10997.cif; but it seems this file only to be "heavy" for some reason). I enclose 4 snapshots in the CIF formats (as the whole trajectory is almost 500 MB) and the code.

import xrayutilities as xu
from xrayutilities.materials.cif import CIFFile
import matplotlib.pyplot as plt 
import matplotlib.ticker as ticker
import numpy as np
from xrayutilities import config
import multiprocessing
import os
import time

config.NTHREADS = 0

def main():
    cifs = ['temp_4551.cif', 'temp_5110.cif', 'temp_6626.cif', 'temp_10999.cif']
    intensities = []
    X = np.arange(0, 45, 0.015)
    for cif in cifs:
        MOF = xu.materials.Crystal.fromCIF(os.path.join('cifs', cif))
        MOF_powder = xu.simpack.Powder(MOF, 1)
        start = time.time()
        pm = xu.simpack.PowderModel(MOF_powder, wl='MoKa1', tt_cutoff=45.0)
        end = time.time()
        print(end - start)
        X = np.arange(0, 45, 0.015)
        start = time.time()
        Y = pm.simulate(X)
        end = time.time()
        print(end - start)
        sf = max(Y) / 100
        intensities.append(Y / sf)

    intensities = np.array(intensities)
    intensities_avg = intensities.mean(axis=0)

    plt.plot(X, intensities_avg)

    plt.xlim(left=0, right=45)
    plt.savefig(f'test_New.png', bbox_inches='tight', dpi=300)

if __name__ == '__main__':
    multiprocessing.freeze_support()
    main()

cifs.zip

Although I played with the nthreads yesterday and I observed no significant speed-up whatsoever, setting the nthreads to 0 now reduces the 30 seconds further to 10-15 seconds (145 seconds for the problematic CIF mentioned above). I am running on Linux (Intel Gold 6240 CPU @ 2.60GHz, 36 threads).

The time needed for the calculation (provided I limit the number of snapshots to average) is pretty much acceptable now. However, when I tried iterating over 50 snapshots, it stopped after 8 with a message RuntimeError: can't start new thread.

[dkriegner]

I need a bit more time to look into all what you sent, but thanks for providing the extra details.

I believe if the number of atoms is constant you should be able to reuse the powdermodel. I need to check how you can update the position of the atoms, but it is certainly possible. I still would want to understand if the difference between CuKa1 and MoKa1 is only the number of peaks or if there is also another issue. One more note: I believe its better you set the tt_cutoff a bit above your simulation limit. In this way you can still get the influence of tails of peaks just above the cutoff.

The issue with the threads is a problem of not closing the PowderModel. The documentation says one should call PowderModel.close() after the end of the use. Somehow the python garbage collector has issues with this.

[dkriegner]

I now ran your example myself. I made few optimizations to the script which now reuses the PowderModel object for the various structures from the CIF file. Using this I get calculation times around 8-12s for each CIF (on my not very powerful notebook). The times do not vary much between the different CIF files. My guess is that because of the large memory use and the fact you did not close old powdermodels you might have run out of memory at some point which caused the much slower execution time for some CIF files.

On my machine I also get somewhat faster execution using more threads. But the improvement is is not overwhelming using my mobile processor.

The code I would propose for you is roughly the following (note that I removed the main function which on Linux/Unix is not needed, you might need to add it back). I also did not change your averaging which I find strange. Not sure if it makes sense to normalize the intensities before averaging. what is the reasoning for this?

import os
import time

import matplotlib.pyplot as plt
import numpy as np
import xrayutilities as xu
from xrayutilities import config
from xrayutilities.materials.cif import CIFFile

config.NTHREADS = 0

cifs = ['temp_4551.cif', 'temp_5110.cif', 'temp_6626.cif', 'temp_10997.cif']
intensities = []
X = np.arange(0, 45, 0.015)
MOF = xu.materials.Crystal.fromCIF(os.path.join('.', cifs[0]))
MOF.name = 'MOF'
MOF_powder = xu.simpack.Powder(MOF, 1)
start = time.time()
pm = xu.simpack.PowderModel(MOF_powder, wl='MoKa1', tt_cutoff=50.0)
end = time.time()
print('pm0', end - start, len(pm.pdiff[0].data))
lmpar = pm.create_fitparameters()
for cif in cifs:
    MOF = xu.materials.Crystal.fromCIF(os.path.join('.', cif))
    # update PowderModel for info from new CIF file
    lmpar['phase_MOF_a'].value = MOF.a
    lmpar['phase_MOF_b'].value = MOF.b
    lmpar['phase_MOF_c'].value = MOF.c
    lmpar['phase_MOF_alpha'].value = MOF.alpha
    lmpar['phase_MOF_beta'].value = MOF.beta
    lmpar['phase_MOF_gamma'].value = MOF.gamma
    for i, atom in enumerate(MOF.lattice.base()):
        name = atom[0].name
        lmpar[f'phase_MOF_at{i}_{name}_1a_occupation'].value = atom[2]
        lmpar[f'phase_MOF_at{i}_{name}_1a_biso'].value = atom[3]
        lmpar[f'phase_MOF_at{i}_{name}_1a_0_pos'].value = atom[1][0]
        lmpar[f'phase_MOF_at{i}_{name}_1a_1_pos'].value = atom[1][1]
        lmpar[f'phase_MOF_at{i}_{name}_1a_2_pos'].value = atom[1][2]
    pm.set_lmfit_parameters(lmpar)
    start = time.time()
    Y = pm.simulate(X)
    end = time.time()
    print(cif, 'sim', end - start, len(pm.pdiff[0].data))
    sf = max(Y) / 100
    intensities.append(Y / sf)

intensities = np.array(intensities)
intensities_avg = intensities.mean(axis=0)
plt.xlim(left=0, right=45)
plt.savefig(f'test_New.png', bbox_inches='tight', dpi=300)

Please closely check that all is really updated which can change in the CIF files. I hope the code works correclty like this.

the original issue of the large increase of execution time when going from CuKa1 to MoKa1: I can reproduce it. Total execution time of the script is roughly 70s for MoKa1 but <10s for CuKa. There are certainly a lot less peaks to consider for Cu radiation. On a quick look it seems this is sufficient to explain the difference. Up to 50deg I find ~920 peaks for Cu but 9200 for Mo. Due to the low symmetry they can't be combined and are calculated individually. I am not sure if there is much which can be done in this case to speed up the calculation without loss of precision. Limiting the angular range obviously helps quite a bit if it is an option.

[reddocksw]

Hi there,

I have managed to generate the patterns using your previous tips. Reusing PowderModel now speeds things up a lot, I also think it can't go any faster at this point. In the end, instead of averaging all the snapshots from the simulation, I decided to reduce the number of the snapshots used to generate the X-ray patterns, so everything is done in roughly 1 hour.

Indeed, normalising the intensities before averaging does not make much sense, we can blame my ignorance, thanks for pointing this out!

Thank you for taking the time to help me!

[dkriegner]

I was also thinking if any of the code above could be moved to some convenience function inside the repository, but it seems to me your application is a very special one and I am currently not sure how a "update_powdermodel_with_new_cif" function should exactly look. It would likely only work for a powdermodel with a single powder and only in cases like yours where the number of atoms is not changing, which are both big limitations.

glad that the suggestions helped. I will close this since I believe it is not an open issue anymore. feel free to reopen if you disagree.