Cell info in netcdf output

[bananenpampe]

Dear Developers,

I have encountered a strange issue when writing trajectories from a non-orthogonal box in the netcdf format.
I find spurious correlations in the rdf functions I calculate from my trajectories, which make me believe that the cell info is not read correctly. I could reproduce the issue using both OVITO and MDanalysis.
Please find attached the trajectory in the .extyz format written by GPUMD which I believe is correct and the .netcdf format. I had to rename the file extensions into .txt to be accepted from github, (dump.txt is the extxyz and movie_4.txt the netcdf)
dump.txt
movie_4.txt

[brucefan1983]

Thanks for reporting this issue. I think the implementation of the netcdf cell in GPUMD is incorrect. If I understand it correctly, the angles are defined to be between the planes according to the manual of netcdf:

https://ambermd.org/netcdf/nctraj.xhtml

But in our code, it was implemented as angles between the cell vectors (https://github.com/brucefan1983/GPUMD/blob/master/src/measure/dump_netcdf.cu):

    const double* t = box.cpu_h;
    double cosgamma, cosbeta, cosalpha;
    cell_lengths[0] = sqrt(t[0] * t[0] + t[3] * t[3] + t[6] * t[6]); // a-side
    cell_lengths[1] = sqrt(t[1] * t[1] + t[4] * t[4] + t[7] * t[7]); // b-side
    cell_lengths[2] = sqrt(t[2] * t[2] + t[5] * t[5] + t[8] * t[8]); // c-side

    cosgamma = (t[0] * t[1] + t[3] * t[4] + t[6] * t[7]) / (cell_lengths[0] * cell_lengths[1]);
    cosbeta = (t[0] * t[2] + t[3] * t[5] + t[6] * t[8]) / (cell_lengths[0] * cell_lengths[2]);
    cosalpha = (t[1] * t[2] + t[4] * t[5] + t[7] * t[8]) / (cell_lengths[1] * cell_lengths[2]);

    cell_angles[0] = acos(cosalpha) * 180.0 / PI;
    cell_angles[1] = acos(cosbeta) * 180.0 / PI;
    cell_angles[2] = acos(cosgamma) * 180.0 / PI;

However, I myself do not use netcdf at all, so I am not in a position to test and correct it. Could you try it out? Otherwise, I can try to get help from some collaborators to check it.

[bananenpampe]

Hello, thank you for the suggestions, yes, I have given it a look:
This is what I tried:

   const double* t = box.cpu_h;
    double cosgamma, cosbeta, cosalpha;
    double normal_ab[3];
    double normal_ac[3];
    double normal_bc[3];
    double mag_normal_ab, mag_normal_ac, mag_normal_bc;
    
    cell_lengths[0] = sqrt(t[0] * t[0] + t[3] * t[3] + t[6] * t[6]); // a-side
    cell_lengths[1] = sqrt(t[1] * t[1] + t[4] * t[4] + t[7] * t[7]); // b-side
    cell_lengths[2] = sqrt(t[2] * t[2] + t[5] * t[5] + t[8] * t[8]); // c-side

    // from the AMBER docs:
    // alpha defines the angle between the a-b and a-c planes
    // beta defines the angle between the a-b and b-c planes
    // gamma defines the angle between the a-c and b-c planes.
    
    normal_ab[0] = t[3] * t[7] - t[6] * t[4];
    normal_ab[1] = t[6] * t[1] - t[0] * t[7];
    normal_ab[2] = t[0] * t[4] - t[3] * t[1];

    normal_ac[0] = t[3] * t[8] - t[6] * t[5];
    normal_ac[1] = t[6] * t[2] - t[0] * t[8];
    normal_ac[2] = t[0] * t[5] - t[3] * t[2];

    normal_bc[0] = t[4] * t[8] - t[7] * t[5];
    normal_bc[1] = t[7] * t[2] - t[1] * t[8];
    normal_bc[2] = t[1] * t[5] - t[4] * t[2];

    mag_normal_bc = sqrt(normal_bc[0] * normal_bc[0] + normal_bc[1] * normal_bc[1] + normal_bc[2] * normal_bc[2]);
    mag_normal_ac = sqrt(normal_ac[0] * normal_ac[0] + normal_ac[1] * normal_ac[1] + normal_ac[2] * normal_ac[2]);
    mag_normal_ab = sqrt(normal_ab[0] * normal_ab[0] + normal_ab[1] * normal_ab[1] + normal_ab[2] * normal_ab[2]);

    cosgamma = (normal_bc[0] * normal_ac[0] + normal_bc[1] * normal_ac[1] + normal_bc[2] * normal_ac[2]) / (mag_normal_bc * mag_normal_ac);
    cosbeta = (normal_bc[0] * normal_ab[0] + normal_bc[1] * normal_ab[1] + normal_bc[2] * normal_ab[2]) / (mag_normal_bc * mag_normal_ab);
    cosalpha = (normal_ac[0] * normal_ab[0] + normal_ac[1] * normal_ab[1] + normal_ac[2] * normal_ab[2]) / (mag_normal_ac * mag_normal_ab);

    cell_angles[0] = acos(cosalpha) * 180.0 / PI;
    cell_angles[1] = acos(cosbeta) * 180.0 / PI;
    cell_angles[2] = acos(cosgamma) * 180.0 / PI;

Unfortunately, this does not fix the issue, is there maybe some convention in how the cell is defined that I am missing?

[brucefan1983]

OK, the major problem is to figure out what conventions AMBER NetCDF uses. The manual is not very clear to me. At least I don't understand why it uses such complicated angles. Is there any reference source code?

[brucefan1983]

Could you try this out:

const double* t = box.cpu_h;
	
const double a[3] = {t[0], t[3], t[6]};
const double b[3] = {t[1], t[4], t[7]};
const double c[3] = {t[2], t[5], t[8]};
	
const double axb = {a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]};
const double axc = {a[1]*c[2]-a[2]*c[1], a[2]*c[0]-a[0]*c[2], a[0]*c[1]-a[1]*c[0]};
const double bxc = {b[1]*c[2]-b[2]*c[1], b[2]*c[0]-b[0]*c[2], b[0]*c[1]-b[1]*c[0]};
	
const double len_axb = sqrt(axb[0]*axb[0] +  axb[1]*axb[1]  + axb[2]*axb[2]);
const double len_axc = sqrt(axc[0]*axc[0] +  axc[1]*axc[1]  + axc[2]*axc[2]);
const double len_bxc = sqrt(bxc[0]*bxc[0] +  bxc[1]*bxc[1]  + bxc[2]*bxc[2]);
	
const double cos_alpha = (axb[0] * axc[0] + axb[1] * axc[1] + axb[2] * axc[2]) / (len_axb * len_axc);
const double cos_beta = -(axb[0] * bxc[0] + axb[1] * bxc[1] + axb[2] * bxc[2]) / (len_axb * len_bxc);
const double cos_gamma = (axc[0] * bxc[0] + axc[1] * bxc[1] + axc[2] * bxc[2]) / (len_axc * len_bxc);
	
cell_lengths[0] = sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
cell_lengths[1] = sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]);
cell_lengths[2] = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
	
cell_angles[0] = acos(cos_alpha) * 180.0 / PI;
cell_angles[1] = acos(cos_beta) * 180.0 / PI;
cell_angles[2] = acos(cos_gamma) * 180.0 / PI;