/* $Id: sggen.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" void ResetLLTr(T_LTr *LLTr, int *nLLTr) { int i; rangei(3) LLTr[0].v[i] = 0; *nLLTr = 1; } static int AddLLTr(int LTBF, int mLLTr, T_LTr *LLTr, int *nLLTr, const int *NewLTr) { int NLTr[3], iLLTr, i; rangei(3) NLTr[i] = iModPositive(NewLTr[i], LTBF); for (iLLTr = 0; iLLTr < *nLLTr; iLLTr++, LLTr++) if (MemCmp(LLTr->v, NLTr, 3) == 0) return 0; if (*nLLTr >= mLLTr) return -1; MemCpy(LLTr->v, NLTr, 3); (*nLLTr)++; return 1; } int ExpLLTr(int LTBF, int mLLTr, T_LTr *LLTr, int *nLLTr, const int *NewLTr) { int TrialLTr[3], i; int iLLTr, jLLTr; const T_LTr *LLTri, *LLTrj; iLLTr = *nLLTr; LLTri = &LLTr[iLLTr]; jLLTr = 1; LLTrj = &LLTr[jLLTr]; for (;;) { if (NewLTr) { if (AddLLTr(LTBF, mLLTr, LLTr, nLLTr, NewLTr) < 0) return -1; } if (jLLTr > iLLTr) { iLLTr++; LLTri++; jLLTr = 1; LLTrj = &LLTr[jLLTr]; } if (iLLTr == *nLLTr) break; rangei(3) TrialLTr[i] = LLTrj->v[i] + LLTri->v[i]; NewLTr = TrialLTr; jLLTr++; LLTrj++; } return 0; } void ResetSgOps(T_SgOps *SgOps) { SgOps->NoExpand = 0; SgOps->nLSL = 1; SgOps->nSSL = 1; MemSet(SgOps->LTr, SgOps_mLTr, 0); ResetLLTr(SgOps->LTr, &SgOps->nLTr); SgOps->fInv = 1; IntSetZero(SgOps->InvT, 3); SgOps->nSMx = 1; MemSet(SgOps->SMx, SgOps_mSMx, 0); InitRTMx(&SgOps->SMx[0], 1); } static int AddSgLTr(T_SgOps *SgOps, const int *NewLTr) { int Stat; Stat = AddLLTr(STBF, SgOps_mLTr, SgOps->LTr, &SgOps->nLTr, NewLTr); if (Stat < 0) { SetSgError("Internal Error: SgOps_mLTr too small"); return -1; } return Stat; } static int AddLtrDueToInvT(T_SgOps *SgOps, const T_RTMx *LSMx) { int NewLTr[3], i; RotMx_t_Vector(NewLTr, LSMx->s.R, SgOps->InvT, 0); for (i = 0; i < 3; i++) NewLTr[i] += 2 * LSMx->s.T[i] - SgOps->InvT[i]; return AddSgLTr(SgOps, NewLTr); } static int AddSgInv(T_SgOps *SgOps, const int *InvT) { int NewLTr[3], i; int iLSMx; T_RTMx *LSMx; const int NNN[] = { 0, 0, 0}; if (! InvT) InvT = NNN; if (SgOps->fInv == 2) /* there is a centre of inversion already */ { for (i = 0; i < 3; i++) NewLTr[i] = SgOps->InvT[i] - InvT[i]; return AddSgLTr(SgOps, NewLTr); } for (i = 0; i < 3; i++) SgOps->InvT[i] = iModPositive(InvT[i], STBF); SgOps->fInv = 2; if (! SgOps->NoExpand) { LSMx = &SgOps->SMx[1]; for (iLSMx = 1; iLSMx < SgOps->nSMx; iLSMx++, LSMx++) if (AddLtrDueToInvT(SgOps, LSMx) < 0) return -1; } return 1; } static int AddSgSMx(T_SgOps *SgOps, const T_RTMx *NewSMx) { int mR[9], NewLTr[3], InvT[3], i; int iLSMx; T_RTMx *LSMx; for (i = 0; i < 9; i++) mR[i] = -NewSMx->s.R[i]; LSMx = SgOps->SMx; for (iLSMx = 0; iLSMx < SgOps->nSMx; iLSMx++, LSMx++) { if (MemCmp(LSMx->s.R, NewSMx->s.R, 9) == 0) { for (i = 0; i < 3; i++) NewLTr[i] = LSMx->s.T[i] - NewSMx->s.T[i]; return AddSgLTr(SgOps, NewLTr); } if (MemCmp(LSMx->s.R, mR, 9) == 0) { for (i = 0; i < 3; i++) InvT[i] = LSMx->s.T[i] + NewSMx->s.T[i]; return AddSgInv(SgOps, InvT); } } if (SgOps->nSMx >= SgOps_mSMx) { SetSgError("Error: Non-crystallographic rotation matrix encountered"); return -1; } MemCpy(LSMx->s.R, NewSMx->s.R, 9); for (i = 0; i < 3; i++) LSMx->s.T[i] = iModPositive(NewSMx->s.T[i], STBF); SgOps->nSMx++; if ( ! SgOps->NoExpand && SgOps->fInv == 2 && AddLtrDueToInvT(SgOps, LSMx) < 0) return -1; return 1; } static int DoMulSMxLTr(T_SgOps *SgOps, int iLSMx, int iLLTr, int OldOnly) { const T_RTMx *LSMxi; const T_LTr *LLTri; int NewLTr[3]; const int iLLTr0 = iLLTr; LSMxi = &SgOps->SMx[iLSMx]; for (; iLSMx < SgOps->nSMx; iLSMx++, LSMxi++) { LLTri = &SgOps->LTr[iLLTr0]; for (iLLTr = iLLTr0; iLLTr < (OldOnly ? SgOps->nLSL : SgOps->nLTr); iLLTr++, LLTri++) { RotMx_t_Vector(NewLTr, LSMxi->s.R, LLTri->v, 0); if (AddSgLTr(SgOps, NewLTr) < 0) return -1; } } return 0; } int ExpSgLTr(T_SgOps *SgOps, const int *NewLTr) { int TrialLTr[3], i; int iLLTr, jLLTr; const T_LTr *LLTri, *LLTrj; if (SgOps->NoExpand) { if (NewLTr) return AddSgLTr(SgOps, NewLTr); return 0; } if (DoMulSMxLTr(SgOps, SgOps->nSSL, 1, 1) < 0) return -1; SgOps->nSSL = SgOps->nSMx; iLLTr = SgOps->nLSL; LLTri = &SgOps->LTr[iLLTr]; jLLTr = 1; LLTrj = &SgOps->LTr[jLLTr]; for (;;) { if (NewLTr) { if (AddSgLTr(SgOps, NewLTr) < 0) return -1; } if (DoMulSMxLTr(SgOps, 1, SgOps->nLSL, 0) < 0) return -1; SgOps->nLSL = SgOps->nLTr; if (jLLTr > iLLTr) { iLLTr++; LLTri++; jLLTr = 1; LLTrj = &SgOps->LTr[jLLTr]; } if (iLLTr == SgOps->nLTr) break; for (i = 0; i < 3; i++) TrialLTr[i] = LLTrj->v[i] + LLTri->v[i]; NewLTr = TrialLTr; jLLTr++; LLTrj++; } return 0; } int ExpSgInv(T_SgOps *SgOps, const int *InvT) { if (AddSgInv(SgOps, InvT) < 0) return -1; return ExpSgLTr(SgOps, NULL); } int ExpSgSMx(T_SgOps *SgOps, const T_RTMx *NewSMx) { int iLSMx, jLSMx; const T_RTMx *LSMxi, *LSMxj; T_RTMx TrialSMx[1]; if (SgOps->NoExpand) { if (NewSMx) return AddSgSMx(SgOps, NewSMx); return 0; } iLSMx = SgOps->nSMx; LSMxi = &SgOps->SMx[iLSMx]; jLSMx = 1; LSMxj = &SgOps->SMx[jLSMx]; for (;;) { if (NewSMx && AddSgSMx(SgOps, NewSMx) < 0) return -1; if (jLSMx > iLSMx) { iLSMx++; LSMxi++; jLSMx = 1; LSMxj = &SgOps->SMx[jLSMx]; } if (iLSMx == SgOps->nSMx) break; SeitzMxMultiply(TrialSMx, LSMxj, LSMxi); NewSMx = TrialSMx; jLSMx++; LSMxj++; } return ExpSgLTr(SgOps, NULL); } int ExpSgRMx(T_SgOps *SgOps, const int NewRMx[9]) { T_RTMx SMx[1]; MemCpy(SMx->s.R, NewRMx, 9); IntSetZero(SMx->s.T, 3); return ExpSgSMx(SgOps, SMx); }