Index exclusion by default in ncjoin.py

[jotelha]

@pastewka is there are particular reason for always excluding the "id" field from joint trajectories by default,

https://github.com/IMTEK-Simulation/code-snippets/blob/5d8e6a2f2f4e29d43db9d99150683de2f382ffd8/imteksimcs/cli/NetCDF/ncjoin.py#L128 ?

That renders the atoms unidentifiable after joining, I would like to remove that line, do you object?

As an example, the structure of a single netcdf file generated by LAMMPS might look like this,

$ ncdump -c single.nc
netcdf single {
dimensions:
    frame = UNLIMITED ; // (226 currently)
    atom = 583294 ;
    cell_spatial = 3 ;
    cell_angular = 3 ;
    label = 10 ;
    spatial = 3 ;
    Voigt = 6 ;
variables:
    char spatial(spatial) ;
    char cell_spatial(spatial) ;
    char cell_angular(spatial, label) ;
    double time(frame) ;
        time:units = "femtosecond" ;
        time:scale_factor = 2. ;
    double cell_origin(frame, cell_spatial) ;
        cell_origin:units = "Angstrom" ;
        cell_origin:scale_factor = 1. ;
    double cell_lengths(frame, cell_spatial) ;
        cell_lengths:units = "Angstrom" ;
        cell_lengths:scale_factor = 1. ;
    double cell_angles(frame, cell_angular) ;
        cell_angles:units = "degree" ;
    int id(frame, atom) ;
    int mol(frame, atom) ;
    int type(frame, atom) ;
    float mass(frame, atom) ;
    float coordinates(frame, atom, spatial) ;
    float velocities(frame, atom, spatial) ;
    float forces(frame, atom, spatial) ;
    float f_storeAnteShakeForces(frame, atom, spatial) ;
    float f_storeAnteStatForces(frame, atom, spatial) ;
    float f_storeAnteFreezeForces(frame, atom, spatial) ;
    float f_storeUnconstrainedForces(frame, atom, spatial) ;
    float f_storeForcesAve(frame, atom, spatial) ;
    float f_storeAnteShakeForcesAve(frame, atom, spatial) ;
    float f_storeAnteStatForcesAve(frame, atom, spatial) ;
    float f_storeAnteFreezeForcesAve(frame, atom, spatial) ;
    float f_storeUnconstrainedForcesAve(frame, atom, spatial) ;
    float c_peratom_stress(frame, atom, Voigt) ;
    float f_peratom_stress_ave(frame, atom, Voigt) ;

// global attributes:
        :Conventions = "AMBER" ;
        :ConventionVersion = "1.0" ;
        :program = "LAMMPS" ;
        :programVersion = "29 Oct 2020" ;
data:

 spatial = "xyz" ;

 cell_spatial = "abc" ;
}

and the structure of multiple such files concatenated via ncjoin.py like this,

$ ncdump -c joint.nc
netcdf joint {
dimensions:
    spatial = 3 ;
    label = 10 ;
    frame = UNLIMITED ; // (1751 currently)
    cell_spatial = 3 ;
    cell_angular = 3 ;
    atom = 583294 ;
    Voigt = 6 ;
variables:
    char spatial(spatial) ;
    char cell_spatial(spatial) ;
    char cell_angular(spatial, label) ;
    double time(frame) ;
        time:units = "femtosecond" ;
        time:scale_factor = 2. ;
    double cell_origin(frame, cell_spatial) ;
        cell_origin:units = "Angstrom" ;
        cell_origin:scale_factor = 1. ;
    double cell_lengths(frame, cell_spatial) ;
        cell_lengths:units = "Angstrom" ;
        cell_lengths:scale_factor = 1. ;
    double cell_angles(frame, cell_angular) ;
        cell_angles:units = "degree" ;
    int mol(frame, atom) ;
    int type(frame, atom) ;
    float mass(frame, atom) ;
    float coordinates(frame, atom, spatial) ;
    float velocities(frame, atom, spatial) ;
    float forces(frame, atom, spatial) ;
    float f_storeAnteShakeForces(frame, atom, spatial) ;
    float f_storeAnteStatForces(frame, atom, spatial) ;
    float f_storeAnteFreezeForces(frame, atom, spatial) ;
    float f_storeUnconstrainedForces(frame, atom, spatial) ;
    float f_storeForcesAve(frame, atom, spatial) ;
    float f_storeAnteShakeForcesAve(frame, atom, spatial) ;
    float f_storeAnteStatForcesAve(frame, atom, spatial) ;
    float f_storeAnteFreezeForcesAve(frame, atom, spatial) ;
    float f_storeUnconstrainedForcesAve(frame, atom, spatial) ;
    float c_peratom_stress(frame, atom, Voigt) ;
    float f_peratom_stress_ave(frame, atom, Voigt) ;

// global attributes:
        :Conventions = "AMBER" ;
        :ConventionVersion = "1.0" ;
        :program = "LAMMPS" ;
        :programVersion = "29 Oct 2020" ;
data:

 spatial = "" ;

 cell_spatial = "" ;
}

with the only important difference at the id field,

$ diff single.out joint.out 
1c1
< netcdf single {
---
> netcdf joint {
3,4c3,5
<   frame = UNLIMITED ; // (226 currently)
<   atom = 583294 ;
---
>   spatial = 3 ;
>   label = 10 ;
>   frame = UNLIMITED ; // (1751 currently)
7,8c8
<   label = 10 ;
<   spatial = 3 ;
---
>   atom = 583294 ;
25d24
<   int id(frame, atom) ;
51c50
<  spatial = "xyz" ;
---
>  spatial = "" ;
53c52
<  cell_spatial = "abc" ;
---
>  cell_spatial = "" ;

I want them both to have the id field.

[pastewka]

ncjoin operates under the assumption that id is only necessary if the atoms are not sorted. ncjoin sorts the atoms, therefore you no longer require id. The id is then identical to the index of the atom.

This is of course only true if there are no gaps in the ids, which may happen if you remove atoms at some point during the simulation.

[jotelha]

I see. I have introduced those modifications, https://github.com/IMTEK-Simulation/code-snippets/pull/4. The default behavior of excluding id won't change.

[pastewka]

Sounds good...

[pastewka]

Fixed with #4