sghkl.c

/* $Id: sghkl.c 851 2002-02-09 05:39:17Z wdelano $ */

/* The source code contained in this file is            */
/* Copyright (C) 1994-2000 by Ralf W. Grosse-Kunstleve. */
/* Please see the LICENSE file for more information.    */

#ifdef _PYMOL_WIN32
#include"os_predef.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>


#undef SG_GLOBAL
#include "sglite.h"


#define HmulR(Hm, H, R) \
  (Hm[0]) = (R)[0] * (H[0]) + (R)[3] * (H[1]) + (R)[6] * (H[2]); \
  (Hm[1]) = (R)[1] * (H[0]) + (R)[4] * (H[1]) + (R)[7] * (H[2]); \
  (Hm[2]) = (R)[2] * (H[0]) + (R)[5] * (H[1]) + (R)[8] * (H[2])


static int Is000(const int H[3])
{
  int  i;
  rangei(3) if (H[i]) return 0;
  return 1;
}


static int AreSameMIx(const int H1[3], const int H2[3])
{
  int  i;
  rangei(3) if (H1[i] != H2[i]) return 0;
  return 1;
}


static int AreFriedelMates(const int H1[3], const int H2[3])
{
  int  i;
  rangei(3) if (H1[i] != -H2[i]) return 0;
  return 1;
}


int CmpEqMIx(const int H1[3], const int H2[3])
{
  int        i;
  const int  P[3] = { 2, 0, 1 };

  rangei(3) {
    if (H1[P[i]] >= 0 && H2[P[i]] <  0) return -1;
    if (H1[P[i]] <  0 && H2[P[i]] >= 0) return  1;
  }

  rangei(3) {
    if (abs(H1[P[i]]) < abs(H2[P[i]])) return -1;
    if (abs(H1[P[i]]) > abs(H2[P[i]])) return  1;
  }

  return 0;
}


int IsSysAbsMIx(const T_SgOps *SgOps, const int H[3], int *TH_Restriction)
{
  int        iSMx, iLTr, i;
  int        Hm[3];
  int        TH, InvTmT[3];
  const int  *T, *TS, *TR;

  /* systematically absent reflection: if HR == H and HT != 0 mod 1
     restricted phase: if HR == -H: phi(H) = pi*HT + n*pi
   */

  if (TH_Restriction) *TH_Restriction = -1;

  range1(iSMx, SgOps->nSMx)
  {
    HmulR(Hm, H, SgOps->SMx[iSMx].s.R);
             T = SgOps->SMx[iSMx].s.T;
    TS = NULL;
    TR = NULL;

    if      (AreSameMIx(H, Hm)) {
      TS = T;
      if (TH_Restriction && SgOps->fInv == 2) {
        rangei(3) InvTmT[i] = SgOps->InvT[i] - T[i];
        TR = InvTmT;
      }
    }
    else if (AreFriedelMates(H, Hm)) {
      if (TH_Restriction) TR = T;
      if (SgOps->fInv == 2) {
        rangei(3) InvTmT[i] = SgOps->InvT[i] - T[i];
        TS = InvTmT;
      }
    }

    if (TS) {
      range1(iLTr, SgOps->nLTr) {
        TH = 0;
        rangei(3) TH += (TS[i] + SgOps->LTr[iLTr].v[i]) * H[i];
        TH %= STBF;
        if (TH != 0) return 1;
      }
    }

    if (TR) {
      range1(iLTr, SgOps->nLTr) {
        TH = 0;
        rangei(3) TH += (TR[i] + SgOps->LTr[iLTr].v[i]) * H[i];
        TH %= STBF; if (TH < 0) TH += STBF;
        if      (*TH_Restriction < 0) *TH_Restriction = TH;
        else if (*TH_Restriction != TH) return 1;
      }
    }
  }

  return 0;
}


int IsCentricMIx(const T_SgOps *SgOps, const int H[3])
{
  int  iSMx, Hm[3];

  if (SgOps->fInv == 2) return 1;

  range1(iSMx, SgOps->nSMx) {
    HmulR(Hm, H, SgOps->SMx[iSMx].s.R);
    if (AreFriedelMates(H, Hm)) return 1;
  }

  return 0;
}


int GetPhaseRestriction(const T_SgOps *SgOps, const int H[3])
{
  int        TH, iSMx, Hm[3], i;
  const int  *TR;

  /* restricted phase: if HR == -H: phi(H) = pi*HT + n*pi
   */

  TH = -1;
  TR = NULL;

  if (SgOps->fInv == 2) {
    TR = SgOps->InvT;
  }
  else {
    range1(iSMx, SgOps->nSMx) {
      HmulR(Hm, H, SgOps->SMx[iSMx].s.R);
      if (AreFriedelMates(H, Hm)) {
        TR = SgOps->SMx[iSMx].s.T;
        break;
      }
    }
  }

  if (TR) {
    TH = 0;
    rangei(3) TH += TR[i] * H[i];
    TH %= STBF; if (TH < 0) TH += STBF;
  }

  return TH;
}


int EpsilonMIx(const T_SgOps *SgOps, const int H[3])
{
  int  Epsilon;
  int  iSMx, Hm[3];

  Epsilon = 0;

  range1(iSMx, SgOps->nSMx) {
    HmulR(Hm, H, SgOps->SMx[iSMx].s.R);
    if (AreSameMIx(H, Hm) || (SgOps->fInv == 2 && AreFriedelMates(H, Hm)))
      Epsilon++;
  }

  if (Epsilon == 0 || SgOps->nSMx % Epsilon != 0) return IE(-1);

  return Epsilon;
}


int MultMIx(const T_SgOps *SgOps, int FriedelSym, const int H[3])
{
  int  Centro;
  int  M, R;
  int  iSMx, Hm[3];

  if (Is000(H)) return 1;

  Centro = ((SgOps->fInv == 2) || FriedelSym);

  M = 0;
  R = 0;

  range1(iSMx, SgOps->nSMx) {
    HmulR(Hm, H, SgOps->SMx[iSMx].s.R);
    if      (AreSameMIx(H, Hm)) M++;
    else if (AreFriedelMates(H, Hm)) R++;
  }

  if (M == 0 || SgOps->nSMx % M != 0 || (R != 0 && R != M)) return IE(-1);

  M = SgOps->nSMx / M;

  if (Centro && R == 0)
    M *= 2;

  return M;
}


int BuildEqMIx(const T_SgOps *SgOps, int FriedelSym, const int H[3],
               T_EqMIx *EqMIx)
{
  int      iSMx, Hm[3], iEq, i;
  T_EqMIx  BufEqMIx[1];

  if (EqMIx == NULL)
      EqMIx = BufEqMIx;

  EqMIx->fInv = 1;
  if ((SgOps->fInv == 2 || FriedelSym) && ! Is000(H))
    EqMIx->fInv = 2;

  EqMIx->N = 0;

  range1(iSMx, SgOps->nSMx)
  {
    HmulR(Hm, H, SgOps->SMx[iSMx].s.R);

    range1(iEq, EqMIx->N) {
      if (AreSameMIx(Hm, EqMIx->H[iEq])) break;
      if (EqMIx->fInv != 2) continue;
      if (AreFriedelMates(Hm, EqMIx->H[iEq])) break;
    }

    if (iEq == EqMIx->N)
    {
      if (EqMIx->N >= sizeof EqMIx->TH / sizeof (*EqMIx->TH))
        return IE(0);

      MemCpy(EqMIx->H[iEq], Hm, 3);

          EqMIx->TH[iEq] = 0;
      rangei(3)
          EqMIx->TH[iEq] += SgOps->SMx[iSMx].s.T[i] * H[i];
          EqMIx->TH[iEq] %= STBF;
      if (EqMIx->TH[iEq] < 0)
          EqMIx->TH[iEq] += STBF;

      EqMIx->N++;
    }
  }

  if (SgOps->nSMx % EqMIx->N) return IE(0); /* another error trap */

  return EqMIx->fInv * EqMIx->N;
}


static int OneMxCutPRange(const int R[9])
{
  int  m, ir, ic, s;

  m = 0;

  range1(ic, 3) {
    s = 0;
    range1(ir, 3) s += abs(R[ir * 3 + ic]);
    if (m < s)
        m = s;
  }

  return m + 1;
}


static int SetCheckCutPRange(const T_SgOps *SgOps)
{
  int  Range, iSMx, m;

  Range = 0;

  range1(iSMx, SgOps->nSMx) {
                m = OneMxCutPRange(SgOps->SMx[iSMx].s.R);
    if (Range < m)
        Range = m;
  }

  return Range;
}


static int CheckCutParam(const T_SgOps *SgOps, int FriedelSym, int CutP[3],
                         int Range, int FullBlock)
{
  int      iBV, iEq, i;
  int      AdjRange[3], Step[3], H[3];
  T_EqMIx  EqMIx[1];

  rangei(3) AdjRange[i] = Range;

  range1(iBV, 3)
  {
    rangei(3) Step[i] = 1;
    if (FullBlock == 0) Step[iBV] = 2 * Range;

    for (H[0] = -AdjRange[0]; H[0] <= AdjRange[0]; H[0] += Step[0])
    for (H[1] = -AdjRange[1]; H[1] <= AdjRange[1]; H[1] += Step[1])
    for (H[2] = -AdjRange[2]; H[2] <= AdjRange[2]; H[2] += Step[2])
    {
      if (BuildEqMIx(SgOps, FriedelSym, H, EqMIx) < 1)
        return IE(-1);

      /* search for equivalent hkl in an active octant */
      range1(iEq, EqMIx->N) {
        rangei(3) if (CutP[i] == 0 && EqMIx->H[iEq][i] < 0) break;
        if (i == 3) break;
        if (EqMIx->fInv != 2) continue;
        rangei(3) if (CutP[i] == 0 && EqMIx->H[iEq][i] > 0) break;
        if (i == 3) break;
      }

      if (iEq == EqMIx->N) return 0; /* CutParam does not work */
    }

    if (FullBlock != 0) break;

    AdjRange[iBV]--;
  }

  return 1; /* CutParam works fine */
}


int GetCutParamMIx(const T_SgOps *SgOps, int FriedelSym, int CutP[3])
{
  int  Range, iTrial, status, i;

  static int ListTrialCutP[7][3] =
    {
      {  0,  0,  0 },
      {  0, -1,  0 },
      { -1,  0,  0 },
      {  0,  0, -1 },
      { -1, -1,  0 },
      {  0, -1, -1 },
      { -1,  0, -1 }
    };

  Range = SetCheckCutPRange(SgOps);
#ifdef JUNK
printf("Range=%d\n", Range);
#endif

  range1(iTrial, 7) {
        status = CheckCutParam(SgOps, FriedelSym, ListTrialCutP[iTrial],
                               Range, 0);
    if (status < 0) return IE(-1);
#ifdef JUNK
{
  int stat2;
  stat2 = CheckCutParam(SgOps, FriedelSym, ListTrialCutP[iTrial], Range + 2, 1);
  if (stat2 != status) {
    printf("CutP %d %d %d\n",
      ListTrialCutP[iTrial][0],
      ListTrialCutP[iTrial][1],
      ListTrialCutP[iTrial][2]);
    DumpSgOps(SgOps, stdout);
    return IE(-1);
  }
}
#endif
    if (status > 0) {
      MemCpy(CutP, ListTrialCutP[iTrial], 3);
      return 0;
    }
  }

  rangei(3) CutP[i] = -1;

  return 0;
}


static int IsInActiveArea(const int CutP[3], const int H[3])
{
  int  i;

  rangei(3) if (CutP[i] == 0 && H[i] < 0) return 0;
  return 1;
}


int GetMasterMIx(const T_EqMIx *EqMIx, const int CutP[3], int MasterH[3])
{
  int  iEq, iInv, i;
  int  HaveMaster, EqH[3];

  HaveMaster = 0;

  range1(iEq, EqMIx->N) {
    MemCpy(EqH, EqMIx->H[iEq], 3);
    range1(iInv, EqMIx->fInv) {
      if (iInv) rangei(3) EqH[i] *= -1;
      if (IsInActiveArea(CutP, EqH)) {
        if (HaveMaster == 0 || CmpEqMIx(MasterH, EqH) > 0) {
          MemCpy(MasterH, EqH, 3);
          HaveMaster = 1;
        }
      }
    }
  }

  if (HaveMaster == 0) return IE(-1);

  return 0;
}


int GetMasterMIx_and_MateID(const T_SgOps *SgOps,
                            const int CutP[3], const int MIx[3],
                            int MasterMIx[3], int *MateID)
{
  int      MateMIx[3], MateMasterMIx[3], i;
  T_EqMIx  EqMIx[1];

  if (BuildEqMIx(SgOps, 0, MIx, EqMIx) == 0) return IE(-1);
  if (GetMasterMIx(EqMIx, CutP, MasterMIx) != 0) return IE(-1);
  *MateID = 0;
  if (SgOps->fInv == 1) {
    rangei(3) MateMIx[i] = -MIx[i];
    if (BuildEqMIx(SgOps, 0, MateMIx, EqMIx) == 0) return IE(-1);
    if (GetMasterMIx(EqMIx, CutP, MateMasterMIx) != 0) return IE(-1);
    if (CmpEqMIx(MasterMIx, MateMasterMIx) > 0) {
      rangei(3) MasterMIx[i] = MateMasterMIx[i];
      *MateID = 1;
    }
  }

  return 0;
}