Jungfraujoch/common/CrystalLattice.cpp

// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include "JFJochMath.h"
#include <cmath>
#include "CrystalLattice.h"
#include "JFJochException.h"

#include "gemmi/symmetry.hpp"
#include "gemmi/unitcell.hpp"
#include "gemmi/cellred.hpp"

#define DEG_TO_RAD static_cast<float>(PI/180.0)

CrystalLattice::CrystalLattice(const UnitCell &cell) {
    vec[0] = {cell.a, 0, 0};
    vec[1] = {cell.b * cosf(cell.gamma * DEG_TO_RAD), cell.b * sinf(cell.gamma * DEG_TO_RAD), 0};
    float cx = cell.c * cosf(cell.beta * DEG_TO_RAD);
    float cy = cell.c
                * (cosf(cell.alpha * DEG_TO_RAD) - cosf(cell.beta * DEG_TO_RAD) * cosf(cell.gamma * DEG_TO_RAD))
                / sinf(cell.gamma * DEG_TO_RAD);
    vec[2] = {cx, cy, sqrtf(cell.c*cell.c-cx*cx-cy*cy)};

    FixHandedness();
}

CrystalLattice::CrystalLattice(const Coord &a, const Coord &b, const Coord &c) {
    vec[0] = a;
    vec[1] = b;
    vec[2] = c;

    FixHandedness();
}

const Coord &CrystalLattice::Vec0() const {
    return vec[0];
}

const Coord &CrystalLattice::Vec1() const {
    return vec[1];
}

const Coord &CrystalLattice::Vec2() const {
    return vec[2];
}

UnitCell CrystalLattice::GetUnitCell() const {
    UnitCell cell{};
    cell.a = vec[0].Length();
    cell.b = vec[1].Length();
    cell.c = vec[2].Length();
    cell.alpha = angle_deg(vec[1], vec[2]);
    cell.beta = angle_deg(vec[0], vec[2]);
    cell.gamma = angle_deg(vec[0], vec[1]);
    return cell;
}

std::vector<float> CrystalLattice::GetVector() const {
    std::vector<float> output(9);
    for (int i = 0; i < 3; i++) {
        output[3 * i + 0] = vec[i].x;
        output[3 * i + 1] = vec[i].y;
        output[3 * i + 2] = vec[i].z;
    }
    return output;
}

float CrystalLattice::CalcVolume() const {
    // Calculate the cell volume
    // V = a · (b × c)
    Coord cross_product = vec[1] % vec[2];
    return vec[0] * cross_product;
}

float CrystalLattice::VolumeFraction() const {
    const float denom = vec[0].Length() * vec[1].Length() * vec[2].Length();
    if (!(denom > 0.0f))
        return 0.0f;
    return std::fabs(CalcVolume()) / denom;
}

void CrystalLattice::Sort() {
    if (vec[0].Length() > vec[1].Length())
        std::swap(vec[0], vec[1]);
    if (vec[1].Length() > vec[2].Length())
        std::swap(vec[1], vec[2]);
    if (vec[0].Length() > vec[1].Length())
        std::swap(vec[0], vec[1]);
}

void CrystalLattice::FlipSign(size_t i1) {
    if (i1 >= 3)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
            "index out of range (0..2)");
    vec[i1] *= -1;
}

void CrystalLattice::FixHandedness() {
    if (CalcVolume() < 0)
        FlipSign(2);
}


// The reciprocal basis divides by the cell volume, so three coplanar rows make all three vectors
// infinite - and an infinite a* does not fail, it quietly predicts nothing and poisons every
// residual built from it. Say so instead.
void CrystalLattice::CheckHasReciprocal() const {
    if (VolumeFraction() < MIN_BASIS_VOLUME_FRACTION)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Crystal lattice is coplanar and has no reciprocal cell");
}

Coord CrystalLattice::Astar() const {
    CheckHasReciprocal();
    return (vec[1] % vec[2]) * (1.0f / CalcVolume());
}

Coord CrystalLattice::Bstar() const {
    CheckHasReciprocal();
    return (vec[2] % vec[0]) * (1.0f / CalcVolume());
}

Coord CrystalLattice::Cstar() const {
    CheckHasReciprocal();
    return (vec[0] % vec[1]) * (1.0f / CalcVolume());
}

CrystalLattice::CrystalLattice(float a, float b, float c, float alpha, float beta, float gamma)
: CrystalLattice(UnitCell{.a = a, .b = b, .c = c, .alpha = alpha, .beta = beta, .gamma = gamma}) {}

CrystalLattice::CrystalLattice(const std::vector<float> &input) {
    if (input.size() != 9)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,"Wrong size of crystal lattice vector");
    for (int i = 0; i < 3; i++) {
        vec[i].x = input[3 * i + 0];
        vec[i].y = input[3 * i + 1];
        vec[i].z = input[3 * i + 2];
    }
}

void CrystalLattice::ReorderABEqual() {
    double la = vec[0].Length();
    double lb = vec[1].Length();
    double lc = vec[2].Length();

    double dab = std::abs(la - lb);
    double dbc = std::abs(lb - lc);
    double dac = std::abs(la - lc);

    Coord a = vec[0];
    Coord b = vec[1];
    Coord c = vec[2];

    if (dbc < dab && dbc < dac) {
        // b≈c → [b, c, a]
        vec[0] = b;
        vec[1] = c;
        vec[2] = a;
    } else if (dac < dab && dac < dbc) {
        // a≈c → [a, c, b]
        vec[0] = a;
        vec[1] = c;
        vec[2] = b;
    } // else a≈b → keep [a, b, c]

    FixHandedness();
}

void CrystalLattice::ReorderMonoclinic() {
    // Enforce obtuse beta (>= 90°). Beta is the angle between a and c.
    // Flip signs of a and b simultaneously to keep handedness and lengths unchanged,
    // which maps beta -> 180° - beta.
    float beta_now = angle_deg(vec[0], vec[2]);
    if (beta_now < 90.0f) {
        vec[0] *= -1.0f; // a -> -a
        vec[1] *= -1.0f; // b -> -b (preserves cell volume sign)
        // beta becomes 180 - beta_now (> 90°)
    }
}

CrystalLattice CrystalLattice::Multiply(const RotMatrix &input) const {
    CrystalLattice l;
    l.vec[0] = input * vec[0];
    l.vec[1] = input * vec[1];
    l.vec[2] = input * vec[2];
    return l;
}

CrystalLattice CrystalLattice::Multiply(const gemmi::Mat33 &c2p) const {
    CrystalLattice l;
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            l.vec[i][j] = c2p[i][0] * vec[0][j] + c2p[i][1] * vec[1][j] + c2p[i][2] * vec[2][j];
        }
    }
    return l;
}

CrystalLattice CrystalLattice::FromPrimitive(char centering) const {
    if (centering == 'P')
        return *this;

    // Transposed for the same reason as ToPrimitive below.
    return Multiply(gemmi::rot_as_mat33(gemmi::centred_to_primitive(centering)).transpose().inverse());
}

CrystalLattice CrystalLattice::ToPrimitive(char centering) const {
    if (centering == 'P')
        return *this;

    // gemmi states the change of basis as an operator on COORDINATES, and Multiply combines BASIS
    // VECTORS, so the matrix has to be transposed - as it already is everywhere else a gemmi Op::Rot
    // reaches Multiply. A, B, C, I and F are symmetric, so the transpose is a no-op for them and the
    // omission never showed; R and H are not, and without it an R-centred lattice was handed back a
    // "primitive" cell that is not that lattice. Its VOLUME was right either way, which is what hid
    // this: a determinant does not change under transposition, and most callers only take the volume.
    return Multiply(gemmi::rot_as_mat33(gemmi::centred_to_primitive(centering)).transpose());
}

void CrystalLattice::Regularize(const gemmi::CrystalSystem &input) {
    switch (input) {
        case gemmi::CrystalSystem::Monoclinic:
            ReorderMonoclinic();
            break;
        case gemmi::CrystalSystem::Tetragonal:
        case gemmi::CrystalSystem::Hexagonal:
            ReorderABEqual();
            break;
        default:
            Sort();
            FixHandedness();
            break;
    }
}

std::vector<float> CrystalLattice::GetUBMatrix() const {
    const Coord astar = Astar();
    const Coord bstar = Bstar();
    const Coord cstar = Cstar();

    return {
        astar.x, bstar.x, cstar.x,
        astar.y, bstar.y, cstar.y,
        astar.z, bstar.z, cstar.z
    };
}

CrystalLattice CrystalLattice::NiggliReduce() const {
    UnitCell uc = GetUnitCell();
    gemmi::UnitCell g_uc(uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
    gemmi::GruberVector g_vec(g_uc, 'P', /*track_change_of_basis=*/true);
    g_vec.niggli_reduce();

    if (g_vec.change_of_basis)
        return Multiply(gemmi::rot_as_mat33(g_vec.change_of_basis->rot).transpose());
    return *this;
}