Saving the LJ forcefield in numpy format with generateLJFF.py raises error

[LauriKurki]

While trying to do relaxed scans with PPSTM, I found that generating the LJ forcefield with python generateLJFF.py -i INPUT_FILE.xyz -f npy seems to be broken at the moment in the main branch, as there is an IndexError: string index out of range raised by line 390 in the attached snippet .

https://github.com/Probe-Particle/ProbeParticleModel/blob/34e6a3603bfe4759fa68a0641ae202c30f64e778/pyProbeParticle/GridUtils.py#L382-L390

This happens because head is a string here, looks like something is wrong either in the indexing or the definition of head? Ultimately the variable head is defined here as atomstring:

https://github.com/Probe-Particle/ProbeParticleModel/blob/34e6a3603bfe4759fa68a0641ae202c30f64e778/pyProbeParticle/HighLevel.py#L150

So the question is: does this function work properly or is there something wrong in the atom-writing function in GridUtils.py?

PS. This feature of saving atoms was not implemented in the old main branch of ProbeParticleModel and generating the force field using the master_backup branch works fine.

Best,
Lauri

[ondrejkrejci]

Yeah, that was my try to put atoms into the npy format, so it can be used in other codes. The reason for this is mainly the problem with cube hartree files from FHI-aims which almost never starts at (0,0,0) - we have an issue about it in PPSTM: Probe-Particle/PPSTM#4 . I guess the best thing would be to change the whole strategy for save_vec_field https://github.com/Probe-Particle/ProbeParticleModel/blob/34e6a3603bfe4759fa68a0641ae202c30f64e778/pyProbeParticle/GridUtils.py#L402, load_vec_field https://github.com/Probe-Particle/ProbeParticleModel/blob/34e6a3603bfe4759fa68a0641ae202c30f64e778/pyProbeParticle/GridUtils.py#L414, save_scal_field https://github.com/Probe-Particle/ProbeParticleModel/blob/34e6a3603bfe4759fa68a0641ae202c30f64e778/pyProbeParticle/GridUtils.py#L430 and 'load_scal_field' https://github.com/Probe-Particle/ProbeParticleModel/blob/34e6a3603bfe4759fa68a0641ae202c30f64e778/pyProbeParticle/GridUtils.py#L442 ; so it would always give us and saves np.array([e, x, y, z], dtype=np.float ) so the infromation about atoms is always propagated forward ...
@NikoOinonen @ProkopHapala @yakutovicha - some ideas about this one?

[LauriKurki]

Okay, saving the atomic positions makes sense, and I guess the variable head contains the positions in xsf format. Does the format here have to be xsf, we already have the atomic positions in an array atoms here

https://github.com/Probe-Particle/ProbeParticleModel/blob/34e6a3603bfe4759fa68a0641ae202c30f64e778/pyProbeParticle/HighLevel.py#L148-L150

Could that be forwarded to save_vec_field instead of atomstring?

[NikoOinonen]

If we are changing save_vec_field, I would also change to using .npz file format, which can hold multiple arrays unlike .npy. So you could change to following lines https://github.com/Probe-Particle/ProbeParticleModel/blob/34e6a3603bfe4759fa68a0641ae202c30f64e778/pyProbeParticle/GridUtils.py#L383-L386 to just be

np.savez(f'{fname}.npz', FFx=FF[:,:,:,0], FFy=FF[:,:,:,1], FFz=FF[:,:,:,2], lvec=lvec)

It would make less files, probably faster IO, and you could also add the atoms information as an additional array.

[ondrejkrejci]

I would use it also for save_scal_field as we always propagate the lattice vectors (and in this issue we want to add there the information about atoms).
Knowing this new information, then I would save the information about atoms as np.array(e,dtype=np.int),np.array([x,y,z],dtype=np.float). This can ease handling of atoms in our PPSTM code.
Now the question is:

• should we go with the lowest amount of changes and just add a new parser, how to go from xsf type head to atoms. Here I am not sure, if this also creates proper head, once someone is starting to create everything from XYZ or some other file ...
• or we should add an exyz object, that would be parsed during the loadXSF ( https://github.com/Probe-Particle/ProbeParticleModel/blob/main/pyProbeParticle/GridUtils.py#L235 ) and loadCUBE ( https://github.com/Probe-Particle/ProbeParticleModel/blob/main/pyProbeParticle/GridUtils.py#L271 ) procedure - this way, we would have the head and exyz together - it is doubling the information in some cases, but it would mean the lowest amount of changes in the code ...
• or finally: to remove the head with the exyz object, and have an XSF-head-creator function in the case, that we would store the fields as xsf - this way, we would not have any doubling of information ...

[NikoOinonen]

or we should add an exyz object

I'm starting to think that using the ASE Atoms object everywhere in the code would really be to our benefit. It accomplishes exactly this functionality of carrying around multiple pieces of information that can be retrieved as needed, and has flexibility in that we can freely add any attributes we want to it without having to add new arguments to functions etc..

[ondrejkrejci]

@NikoOinonen - good option; still we should save the object in between the cases if you run it through the scripts and using the -f npy flag. Optionally we can just use one of the ASE native files ( https://wiki.fysik.dtu.dk/ase/ase/io/io.html ) like the traj file for that.

[ondrejkrejci]

Update the npy part of gpahene-spline test. is not working, because of this. (The whole test is not working at all ... but that is for a different issue)

[yakutovicha]

We've just experienced that bug as well. Can this be fixed anytime soon? Is any help needed here?

[ondrejkrejci]

At fastest on Friday... Sent from Outlook for Android<https://aka.ms/AAb9ysg>

…

________________________________ From: Aliaksandr Yakutovich ***@***.***> Sent: Tuesday, February 14, 2023 6:07:21 PM To: Probe-Particle/ppafm ***@***.***> Cc: Krejci Ondrej ***@***.***>; Comment ***@***.***> Subject: Re: [Probe-Particle/ppafm] Saving the LJ forcefield in numpy format with generateLJFF.py raises error (Issue #53) We've just experienced that bug as well. Can this be fixed anytime soon? Is any help needed here? — Reply to this email directly, view it on GitHub<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FProbe-Particle%2Fppafm%2Fissues%2F53%23issuecomment-1429993787&data=05%7C01%7Ckrejcio1%40aaltofi.mail.onmicrosoft.com%7C22c5390774234395ffcb08db0ea59036%7Cae1a772440414462a6dc538cb199707e%7C1%7C0%7C638119876466635711%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=2iSBmZd8Hnxr7sHfbxCkGl01FufgxaD24tAfrkUII5M%3D&reserved=0>, or unsubscribe<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FADFZAFIGGBXVX7JJJG4WLVLWXOUTTANCNFSM6AAAAAASORB4HQ&data=05%7C01%7Ckrejcio1%40aaltofi.mail.onmicrosoft.com%7C22c5390774234395ffcb08db0ea59036%7Cae1a772440414462a6dc538cb199707e%7C1%7C0%7C638119876466635711%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=ncTqGP%2BuBlenvEspkmHHsUbx%2FJJj28OVveQ3JnCK2og%3D&reserved=0>. You are receiving this because you commented.Message ID: ***@***.***>

[ondrejkrejci]

This goes much deeper than I expected:

If we are running XSF then the XSF_HEAD should be propagated through with the information about atoms and original system lattice.
I wanted to just simply add a possibility for propating just [e,x,y,z] for NPY format, but then the information about original system lattice is immediately lost, which is problem for PPSTM.
I also do not want how to deal with the possibility of double propagation - XSF_HEAD is a string, while in NPY I wanted to propagate numpy array (which could be easily changed for ase.atoms). I can easily propage any of them in one variable and since the format switcher will automatically decide which of them will be used.
Another option is to use just one of them and then parse/construct the other in the case, that I would use non-default option.
@NikoOinonen @ProkopHapala @yakutovicha - some ideas/preferences?

[ondrejkrejci]

@LauriKurki and @yakutovicha :
Please download and maybe test this branch: iss53. positions and original lattice vector is propagated through in both XSF and NPY .
In NPY the information is in (e.g. PPpos) PPpos_atoms.npy and the information is in form of at_array = np.array((number_of_atoms+3,4),floats). E.g.:

array([[ 6.     ,  0.165  ,  0.165  ,  0.165  ],
       [ 6.     ,  2.2987 ,  1.39372,  0.16643],
       [ 6.     ,  4.4293 ,  2.62743,  0.16468],
       [ 6.     ,  6.56208,  3.85683,  0.16431],
       [ 6.     ,  8.69719,  5.08888,  0.16545],
       [ 6.     , 10.83174,  6.32186,  0.16487],
       [ 6.     ,  2.29863, -1.06668,  0.16445],
       [ 6.     ,  4.43016,  0.16415,  0.16499],
       [ 6.     ,  6.5639 ,  1.39513,  0.16469],
       [ 6.     ,  8.69598,  2.62341,  0.1651 ],
       [ 6.     , 10.83364,  3.85745,  0.16497],
       [ 6.     , 12.9634 ,  5.0896 ,  0.16451],
       [ 6.     ,  4.43024, -2.29713,  0.16443],
       [ 6.     ,  6.56097, -1.0666 ,  0.16647],
       [ 6.     ,  8.69598,  0.15848,  0.16406],
       [ 6.     , 10.82917,  1.39705,  0.16555],
       [ 6.     , 12.96875,  2.62522,  0.1658 ],
       [ 6.     , 15.09725,  3.86194,  0.16486],
       [ 6.     ,  6.56402, -3.52838,  0.16595],
       [ 6.     ,  8.69593, -2.29218,  0.16398],
       [ 7.     , 10.80717, -1.06699,  0.16492],
       [ 7.     , 12.9742 ,  0.18437,  0.16464],
       [ 6.     , 15.09093,  1.40025,  0.16571],
       [ 6.     , 17.22761,  2.62833,  0.16442],
       [ 6.     ,  8.6958 , -4.75743,  0.16541],
       [ 6.     , 10.82899, -3.53111,  0.16548],
       [ 7.     , 12.97393, -2.31803,  0.16501],
       [ 6.     , 15.09692, -1.06662,  0.16467],
       [ 6.     , 17.22548,  0.16756,  0.16569],
       [ 6.     , 19.36005,  1.39865,  0.16487],
       [ 6.     , 10.83305, -5.99168,  0.16521],
       [ 6.     , 12.96825, -4.75927,  0.16496],
       [ 6.     , 15.09046, -3.53355,  0.16496],
       [ 6.     , 17.22521, -2.30088,  0.16503],
       [ 6.     , 19.35976, -1.06689,  0.16474],
       [ 6.     , 21.49483,  0.16594,  0.16528],
       [ 6.     ,  1.58723,  0.16476,  0.16539],
       [ 6.     ,  3.71933,  1.39559,  0.16547],
       [ 6.     ,  5.85277,  2.62757,  0.16414],
       [ 6.     ,  7.98679,  3.85576,  0.1651 ],
       [ 6.     , 10.12101,  5.08769,  0.16497],
       [ 6.     , 12.25153,  6.32169,  0.16458],
       [ 6.     ,  3.71994, -1.06648,  0.16432],
       [ 6.     ,  5.85253,  0.16432,  0.16569],
       [ 6.     ,  7.98668,  1.39415,  0.16439],
       [ 6.     , 10.12349,  2.61973,  0.16553],
       [ 6.     , 12.25335,  3.85752,  0.16464],
       [ 6.     , 14.3857 ,  5.09067,  0.16427],
       [ 6.     ,  5.85264, -2.2973 ,  0.16596],
       [ 6.     ,  7.98472, -1.06678,  0.16484],
       [ 6.     , 10.12107,  0.15961,  0.16495],
       [ 6.     , 12.25535,  1.39194,  0.16502],
       [ 6.     , 14.38587,  2.62898,  0.16599],
       [ 6.     , 16.51764,  3.86028,  0.16375],
       [ 6.     ,  7.98658, -3.52774,  0.16529],
       [ 6.     , 10.12107, -2.2935 ,  0.16469],
       [ 6.     , 14.37972,  0.16591,  0.16442],
       [ 6.     , 16.51606,  1.39682,  0.16514],
       [ 6.     , 18.65096,  2.62803,  0.16531],
       [ 6.     , 10.12318, -4.75383,  0.16572],
       [ 6.     , 12.25503, -3.52576,  0.16528],
       [ 6.     , 14.37952, -2.29909,  0.16466],
       [ 6.     , 16.50905, -1.06663,  0.16546],
       [ 6.     , 18.64744,  0.16574,  0.16472],
       [ 6.     , 20.78469,  1.3993 ,  0.16505],
       [ 6.     , 12.25281, -5.99166,  0.16443],
       [ 6.     , 14.3853 , -4.76242,  0.16549],
       [ 6.     , 16.51556, -3.52994,  0.16493],
       [ 6.     , 18.64732, -2.29915,  0.16451],
       [ 6.     , 20.78381, -1.06709,  0.16526],
       [ 6.     , 22.91809,  0.16726,  0.165  ],
       [ 0.     , 12.798  , -7.3889 ,  0.     ],
       [ 0.     , 12.798  ,  7.3889 ,  0.     ],
       [ 0.     ,  0.     ,  0.     , 20.     ]])

so in the begging you have [ei,xi,yi,zi] for each atoms, while on last 3 lines you have your lattice vector, with 0 in the 1st column lattice_vector=at_array[-3:,1:]. Is it clear?

[ondrejkrejci]

Ok, so we are switching to the NPZ format here: https://numpy.org/doc/stable/reference/generated/numpy.savez.html