1 //  $Id: mmdb_symop.h $
2 //  =================================================================
3 //
4 //   CCP4 Coordinate Library: support of coordinate-related
5 //   functionality in protein crystallography applications.
6 //
7 //   Copyright (C) Eugene Krissinel 2000-2013.
8 //
9 //    This library is free software: you can redistribute it and/or
10 //    modify it under the terms of the GNU Lesser General Public
11 //    License version 3, modified in accordance with the provisions
12 //    of the license to address the requirements of UK law.
13 //
14 //    You should have received a copy of the modified GNU Lesser
15 //    General Public License along with this library. If not, copies
16 //    may be downloaded from http://www.ccp4.ac.uk/ccp4license.php
17 //
18 //    This program is distributed in the hope that it will be useful,
19 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
20 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21 //    GNU Lesser General Public License for more details.
22 //
23 //  =================================================================
24 //
25 //    12.09.13   <--  Date of Last Modification.
26 //                   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
27 //  -----------------------------------------------------------------
28 //
29 //  **** Module  :   MMDB_SymOp <implementation>
30 //       ~~~~~~~~~
31 //  **** Project :   MacroMolecular Data Base (MMDB)
32 //       ~~~~~~~~~
33 //
34 //  **** Classes :   mmdb::SymOp  ( symmetry operators )
35 //       ~~~~~~~~~   mmdb::SymOps ( container of symmetry operators )
36 //
37 //   (C) E. Krissinel 2000-2013
38 //
39 //  =================================================================
40 //
41 
42 #include <string.h>
43 #include <stdlib.h>
44 #include <math.h>
45 
46 #include "mmdb_symop.h"
47 
48 namespace mmdb  {
49 
50   //  ====================  SymOp  ========================
51 
SymOp()52   SymOp::SymOp() : io::Stream()  {
53     InitSymOp();
54   }
55 
SymOp(io::RPStream Object)56   SymOp::SymOp ( io::RPStream Object ) : io::Stream(Object)  {
57     InitSymOp();
58   }
59 
~SymOp()60   SymOp::~SymOp()  {
61     FreeMemory();
62   }
63 
InitSymOp()64   void SymOp::InitSymOp()  {
65   int i,j;
66     XYZOp = NULL;
67     for (i=0;i<4;i++)  {
68       for (j=0;j<4;j++)
69         T[i][j] = 0.0;
70       T[i][i] = 1.0;
71     }
72   }
73 
FreeMemory()74   void SymOp::FreeMemory()  {
75     if (XYZOp)  delete[] XYZOp;
76     XYZOp = NULL;
77   }
78 
SetSymOp(cpstr XYZOperation)79   int  SymOp::SetSymOp ( cpstr XYZOperation )  {
80   int  i,j;
81 
82     CreateCopy ( XYZOp,XYZOperation );
83     DelSpaces  ( XYZOp );
84 
85     for (i=0;i<4;i++)
86       for (j=0;j<4;j++)
87         T[i][j] = 0.0;
88 
89     i = GetOperation ( 0 );
90     if (!i) i = GetOperation ( 1 );
91     if (!i) i = GetOperation ( 2 );
92     T[3][3] = 1.0;
93 
94     return i;
95 
96   }
97 
GetSymOp()98   pstr SymOp::GetSymOp()  {
99     if (XYZOp)  return XYZOp;
100           else  return pstr("");
101   }
102 
103 
GetOperation(int n)104   int  SymOp::GetOperation ( int n )  {
105   char     L[100];
106   pstr     p1,p2;
107   int      len;
108   realtype V;
109 
110     p1 = XYZOp;
111     p2 = FirstOccurence ( p1,',' );
112     if (!p2)  return SYMOP_WrongSyntax;
113     if (n>0)  {
114       p1 = p2+1;
115       p2 = FirstOccurence ( p1,',' );
116       if (!p2)  return SYMOP_WrongSyntax;
117     }
118     if (n>1)  {
119       p1 = p2+1;
120       p2 = NULL;
121     }
122 
123     if (p2)  *p2 = char(0);
124     strcpy ( L,p1 );
125     if (p2)  *p2 = ',';
126 
127     DelSpaces ( L );
128     if (!L[0])  return SYMOP_WrongSyntax;
129     UpperCase ( L );
130 
131     len = strlen ( L );
132     T[n][0] = 0.0;
133     if (L[0]=='X')  {
134       T[n][0] += 1.0;
135       L[0] = ' ';
136     }
137     do  {
138       p1 = strstr ( L,"+X" );
139       if (p1)  {
140         T[n][0] += 1.0;
141         strncpy ( p1,"  ",2 );
142       }
143     } while (p1);
144     do  {
145       p1 = strstr ( L,"-X" );
146       if (p1)  {
147         T[n][0] -= 1.0;
148         strncpy ( p1,"  ",2 );
149       }
150     } while (p1);
151 
152     T[n][1] = 0.0;
153     if (L[0]=='Y')  {
154       T[n][1] += 1.0;
155       L[0] = ' ';
156     }
157     do  {
158       p1 = strstr ( L,"+Y" );
159       if (p1)  {
160         T[n][1] += 1.0;
161         strncpy ( p1,"  ",2 );
162       }
163     } while (p1);
164     do  {
165       p1 = strstr ( L,"-Y" );
166       if (p1)  {
167         T[n][1] -= 1.0;
168         strncpy ( p1,"  ",2 );
169       }
170     } while (p1);
171 
172     T[n][2] = 0.0;
173     if (L[0]=='Z')  {
174       T[n][2] += 1.0;
175       L[0] = ' ';
176     }
177     do  {
178       p1 = strstr ( L,"+Z" );
179       if (p1)  {
180         T[n][2] += 1.0;
181         strncpy ( p1,"  ",2 );
182       }
183     } while (p1);
184     do  {
185       p1 = strstr ( L,"-Z" );
186       if (p1)  {
187         T[n][2] -= 1.0;
188         strncpy ( p1,"  ",2 );
189       }
190     } while (p1);
191 
192     DelSpaces ( L );
193     if ((int)strlen(L)>=len)  return SYMOP_NotAnOperation;
194 
195     // translational part
196     p1 = L;
197     T[n][3] = strtod ( p1,&p2 );
198     if (*p2=='/')  {
199       p1 = p2+1;
200       V  = strtod ( p1,&p2 );
201       if (V==0.0)  return SYMOP_ZeroDenominator;
202       T[n][3] /= V;
203     }
204 
205     return SYMOP_Ok;
206 
207   }
208 
MakeSign(pstr S,realtype V,realtype & AV)209   void  MakeSign ( pstr S, realtype V, realtype & AV )  {
210   int      l;
211     if (V>0.0)  {
212       l = strlen ( S );
213       if (l>0)  {
214         if (S[l-1]!=',')  {
215           strcat ( S,"+" );
216         }
217       }
218       AV = V;
219     } else if (V<0.0)  {
220       strcat ( S,"-" );
221       AV = -V;
222     } else  {
223       AV = V;
224       return;
225     }
226   }
227 
228 
229   #define  __eps  1.0e-5
230 
GenTranslation(pstr S,realtype V)231   void  GenTranslation ( pstr S, realtype V )  {
232   realtype AV,nAV;
233   char     N[50];
234   int      n,d;
235 
236     if (fabs(V)<=__eps)  return;
237     MakeSign ( S,V,AV );
238 
239     d = 0;
240     n = -1;
241     while ((d<=20) && (n<0))  {
242       d++;
243       nAV = AV*d;
244       n   = mround(nAV);
245       if (fabs(nAV-n)>__eps)  n = -1;
246     }
247 
248     if (d<=1)      sprintf ( N,"%i"    ,n   );
249     else if (n>=0) sprintf ( N,"%i/%i" ,n,d );
250               else sprintf ( N,"%-.10g",AV  );
251     strcat ( S,N );
252 
253   }
254 
GenTransformation(pstr S,realtype V,pstr Axis)255   void  GenTransformation ( pstr S, realtype V, pstr Axis ) {
256   realtype AV,nAV;
257   char     N[50];
258   int      n,d;
259 
260     if (fabs(V)<=__eps)  return;
261     MakeSign ( S,V,AV );
262 
263     if (fabs(AV-1.0)>__eps)  {
264 
265       d = 0;
266       n = -1;
267       while ((d<=20) && (n<0))  {
268         d++;
269         nAV = AV*d;
270         n   = mround(nAV);
271         if (fabs(nAV-n)>__eps)  n = -1;
272       }
273 
274       if (n>=0)  sprintf ( N,"%i/%i*",n,d );
275            else  sprintf ( N,"%-.10g*",AV );
276       strcat ( S,N );
277 
278     }
279 
280     strcat ( S,Axis );
281 
282   }
283 
284 
285   /*
286   void  GenTranslation ( pstr S, realtype V )  {
287   realtype AV,fAV;
288   int      n,d;
289   char     N[50];
290 
291     if (V==0.0)  return;
292     MakeSign ( S,V,AV );
293 
294     n = mround(floor(AV+0.00000001));
295     fAV = AV-n;
296 
297     if (fabs(fAV-0.5)<=__eps)                 { n += 1;  d = 2; }
298     else if (fabs(fAV-0.25)<=__eps)           { n += 1;  d = 4; }
299     else if (fabs(fAV-0.75)<=__eps)           { n += 3;  d = 4; }
300     else if (fabs(fAV-0.33333333333)<=__eps)  { n += 1;  d = 3; }
301     else if (fabs(fAV-0.66666666666)<=__eps)  { n += 2;  d = 3; }
302     else if (fabs(fAV-0.16666666666)<=__eps)  { n += 1;  d = 6; }
303     else if (fabs(fAV-0.83333333333)<=__eps)  { n += 5;  d = 6; }
304                                         else  d = 1;
305 
306     N[0] = char(0);
307     if (d>1)       sprintf  ( N,"%i/%i",n,d );
308     else if (n>0)  sprintf  ( N,"%i",n );
309              else  ParamStr ( N,pstr(""),AV );
310     strcat ( S,N );
311 
312   }
313 
314   void  GenTransformation ( pstr S, realtype V, pstr Axis ) {
315   realtype AV;
316 
317     if (V==0.0)  return;
318     MakeSign ( S,V,AV );
319 
320     if (fabs(AV-0.5)<=__eps)                 strcat   ( S,"1/2*" );
321     else if (fabs(AV-0.25)<=__eps)           strcat   ( S,"1/4*" );
322     else if (fabs(AV-0.75)<=__eps)           strcat   ( S,"3/4*" );
323     else if (fabs(AV-0.33333333333)<=__eps)  strcat   ( S,"1/3*" );
324     else if (fabs(AV-0.66666666666)<=__eps)  strcat   ( S,"2/3*" );
325     else if (fabs(AV-0.16666666666)<=__eps)  strcat   ( S,"1/6*" );
326     else if (fabs(AV-0.83333333333)<=__eps)  strcat   ( S,"5/6*" );
327     else if (fabs(AV-1.0)>__eps)             ParamStr ( S,pstr(""),AV,
328                                                         10,pstr("*") );
329 
330     strcat ( S,Axis );
331 
332   }
333 
334   */
335 
CompileOpTitle(pstr S)336   bool  SymOp::CompileOpTitle ( pstr S )  {
337     return CompileOpTitle ( S,T,true );
338   }
339 
CompileOpTitle(pstr S,mat44 symMat,bool compare)340   bool  SymOp::CompileOpTitle ( pstr S, mat44 symMat,
341                                     bool compare )  {
342     S[0] = char(0);
343     GenTransformation ( S,symMat[0][0],pstr("X") );
344     GenTransformation ( S,symMat[0][1],pstr("Y") );
345     GenTransformation ( S,symMat[0][2],pstr("Z") );
346     GenTranslation    ( S,symMat[0][3]           );
347     strcat ( S,"," );
348     GenTransformation ( S,symMat[1][0],pstr("X") );
349     GenTransformation ( S,symMat[1][1],pstr("Y") );
350     GenTransformation ( S,symMat[1][2],pstr("Z") );
351     GenTranslation    ( S,symMat[1][3]           );
352     strcat ( S,"," );
353     GenTransformation ( S,symMat[2][0],pstr("X") );
354     GenTransformation ( S,symMat[2][1],pstr("Y") );
355     GenTransformation ( S,symMat[2][2],pstr("Z") );
356     GenTranslation    ( S,symMat[2][3]           );
357     DelSpaces ( S );
358     if ((!compare) || (!strcmp(S,XYZOp)))  return true;
359     else  {
360       S[0] = char(0);
361       GenTranslation    ( S,symMat[0][3]           );
362       GenTransformation ( S,symMat[0][0],pstr("X") );
363       GenTransformation ( S,symMat[0][1],pstr("Y") );
364       GenTransformation ( S,symMat[0][2],pstr("Z") );
365       strcat ( S,"," );
366       GenTranslation    ( S,symMat[1][3]           );
367       GenTransformation ( S,symMat[1][0],pstr("X") );
368       GenTransformation ( S,symMat[1][1],pstr("Y") );
369       GenTransformation ( S,symMat[1][2],pstr("Z") );
370       strcat ( S,"," );
371       GenTranslation    ( S,symMat[2][3]           );
372       GenTransformation ( S,symMat[2][0],pstr("X") );
373       GenTransformation ( S,symMat[2][1],pstr("Y") );
374       GenTransformation ( S,symMat[2][2],pstr("Z") );
375       DelSpaces ( S );
376       if (!strcmp(S,XYZOp))  return true;
377     }
378     return false;
379   }
380 
Transform(realtype & x,realtype & y,realtype & z)381   void  SymOp::Transform ( realtype & x, realtype & y, realtype & z )  {
382   realtype x1,y1,z1;
383     x1 = T[0][0]*x + T[0][1]*y + T[0][2]*z + T[0][3];
384     y1 = T[1][0]*x + T[1][1]*y + T[1][2]*z + T[1][3];
385     z1 = T[2][0]*x + T[2][1]*y + T[2][2]*z + T[2][3];
386     x = x1;
387     y = y1;
388     z = z1;
389   }
390 
GetTMatrix(mat44 & TMatrix)391   void  SymOp::GetTMatrix ( mat44 & TMatrix )  {
392   // copies T to TMatrix
393   int i,j;
394     for (i=0;i<4;i++)
395       for (j=0;j<4;j++)
396         TMatrix[i][j] = T[i][j];
397   }
398 
SetTMatrix(mat44 & TMatrix)399   void  SymOp::SetTMatrix ( mat44 & TMatrix )  {
400   // copies TMatrix to T
401   int i,j;
402     for (i=0;i<4;i++)
403       for (j=0;j<4;j++)
404         T[i][j] = TMatrix[i][j];
405   }
406 
407 
Print()408   void  SymOp::Print()  {
409   int i;
410     printf ( "  operation '%s'\n",XYZOp );
411     for (i=0;i<4;i++)
412       printf ( "               %10.3g %10.3g %10.3g  %10.3g\n",
413                T[i][0],T[i][1],T[i][2],T[i][3] );
414   }
415 
Copy(PSymOp SymOp)416   void  SymOp::Copy ( PSymOp SymOp )  {
417   int i,j;
418     CreateCopy ( XYZOp,SymOp->XYZOp );
419     for (i=0;i<4;i++)
420       for (j=0;j<4;j++)
421         T[i][j] = SymOp->T[i][j];
422   }
423 
write(io::RFile f)424   void  SymOp::write ( io::RFile f )  {
425   int  i,j;
426   byte Version=1;
427     f.WriteByte   ( &Version );
428     f.CreateWrite ( XYZOp    );
429     for (i=0;i<4;i++)
430       for (j=0;j<4;j++)
431         f.WriteReal ( &(T[i][j]) );
432   }
433 
read(io::RFile f)434   void  SymOp::read ( io::RFile f )  {
435   int  i,j;
436   byte Version;
437     f.ReadByte   ( &Version );
438     f.CreateRead ( XYZOp    );
439     for (i=0;i<4;i++)
440       for (j=0;j<4;j++)
441         f.ReadReal ( &(T[i][j]) );
442   }
443 
MakeStreamFunctions(SymOp)444   MakeStreamFunctions(SymOp)
445 
446 
447 
448   //  ====================  SymOps  ========================
449 
450   SymOps::SymOps() : io::Stream()  {
451     InitSymOps();
452   }
453 
SymOps(io::RPStream Object)454   SymOps::SymOps ( io::RPStream Object ) : io::Stream(Object)  {
455     InitSymOps();
456   }
457 
~SymOps()458   SymOps::~SymOps()  {
459     FreeMemory();
460   }
461 
InitSymOps()462   void SymOps::InitSymOps()  {
463     SpGroup = NULL;
464     Nops    = 0;
465     symOp   = NULL;
466   }
467 
FreeMemory()468   void SymOps::FreeMemory()  {
469   int i;
470     if (SpGroup)  delete[] SpGroup;
471     SpGroup = NULL;
472     if (symOp)  {
473       for (i=0;i<Nops;i++)
474         if (symOp[i])  delete symOp[i];
475       delete[] symOp;
476       symOp = NULL;
477     }
478     Nops = 0;
479   }
480 
481   #define  symop_file  cpstr("symop.lib")
482 
SetGroupSymopLib(cpstr SpaceGroup,cpstr symop_lib)483   int  SymOps::SetGroupSymopLib ( cpstr SpaceGroup,
484                                   cpstr symop_lib )  {
485   char     S[500];
486   char     G[100];
487   pstr     p;
488   io::File f;
489   int      i,RC;
490 
491     FreeMemory();
492 
493     CreateCopy ( SpGroup,SpaceGroup );
494 
495     if (!symop_lib)          p = pstr(symop_file);
496     else if (!symop_lib[0])  p = pstr(symop_file);
497                        else  p = pstr(symop_lib);
498     f.assign ( p,true );
499     if (!f.reset(true))  {
500       p = getenv ( "SYMOP" );
501       if (p)
502         strcpy ( S,p );
503       else  {
504         p = getenv ( "CLIBD" );
505         if (p)  {
506           strcpy ( S,p );
507           if (S[strlen(S)-1]!='/')  strcat ( S,"/" );
508           strcat ( S,"symop.lib" );
509         } else
510           strcpy ( S,"symop.lib" );
511       }
512       f.assign ( S,true );
513       if (!f.reset(true))  return SYMOP_NoLibFile;
514     }
515 
516     strcpy ( G," '"    );
517     strcat ( G,SpGroup );
518     strcat ( G,"'"     );
519     S[0] = char(0);
520     while (!f.FileEnd() && (!strstr(S,G)))
521       f.ReadLine ( S,sizeof(S) );
522     if (f.FileEnd())  {
523       f.shut();
524       return SYMOP_UnknownSpaceGroup;
525     }
526 
527     p = S;
528     while (*p==' ')  p++;
529     p = FirstOccurence ( p,' ' );
530     if (p)  Nops = mround(strtod(p,NULL));
531     if (Nops<=0)  return SYMOP_NoSymOps;
532 
533     symOp = new PSymOp[Nops];
534     RC    = SYMOP_Ok;
535     for (i=0;(i<Nops) && (!RC);i++)  {
536       f.ReadLine ( S,sizeof(S) );
537       symOp[i] = new SymOp();
538       RC = symOp[i]->SetSymOp ( S );
539     }
540 
541     f.shut();
542 
543     return RC;
544 
545   }
546 
547 
548   #define  syminfo_file  cpstr("syminfo.lib")
549 
SetGroup(cpstr SpaceGroup,cpstr syminfo_lib)550   int  SymOps::SetGroup ( cpstr SpaceGroup,
551                            cpstr syminfo_lib )  {
552   io::File f;
553   pstr     p;
554   psvector lines,lines1;
555   char     S[500];
556   char     G[100];
557   char     O[100];
558   mat44    T1,T2,T3;
559   int      i,j,k,l,m,n,RC;
560   int      nlines,npops,ncops;
561 
562     FreeMemory();
563 
564     npops = 0;
565     ncops = 0;
566 
567     CreateCopy ( SpGroup,SpaceGroup );
568 
569     if (!syminfo_lib)          p = pstr(syminfo_file);
570     else if (!syminfo_lib[0])  p = pstr(syminfo_file);
571                          else  p = pstr(syminfo_lib);
572     f.assign ( p,true );
573     if (!f.reset(true))  {
574       p = getenv ( "SYMINFO" );
575       if (p)
576         strcpy ( S,p );
577       else  {
578         p = getenv ( "CLIBD" );
579         if (p)  {
580           strcpy ( S,p );
581           if (S[strlen(S)-1]!='/')  strcat ( S,"/" );
582           strcat ( S,"syminfo.lib" );
583         } else
584           strcpy ( S,"syminfo.lib" );
585       }
586       f.assign ( S,true );
587       if (!f.reset(true))  return SYMOP_NoLibFile;
588     }
589 
590 
591     if (strncasecmp(SpGroup,"Hall:",5))  {
592       // normal space group symbol on input
593       strcpy ( G," '"    );
594       strcat ( G,SpGroup );
595       strcat ( G,"'"     );
596       S[0] = char(0);
597       while (!f.FileEnd() && !(strstr(S,G) && (strstr(S,"symbol xHM") ||
598               strstr(S,"symbol old"))))
599         f.ReadLine ( S,sizeof(S) );
600     } else  {
601       // hall descriptor on input
602       strcpy ( G," ' " );
603       p = &(SpGroup[5]);
604       while (*p==' ')  p++;
605       strcat ( G,p     );
606       strcat ( G,"'"   );
607       S[0] = char(0);
608       while (!f.FileEnd() && !(strstr(S,G) && strstr(S,"symbol Hall")))
609         f.ReadLine ( S,sizeof(S) );
610     }
611     if (f.FileEnd())  {
612       f.shut();
613       return SYMOP_UnknownSpaceGroup;
614     }
615 
616     // found spacegroup, move to symop lines
617     while (!f.FileEnd() && (!strstr(S,"symop")))
618       f.ReadLine ( S,sizeof(S) );
619 
620     nlines = 256;
621     GetVectorMemory ( lines,nlines,0 );
622     for (i=0;i<nlines;i++)
623       lines[i] = NULL;
624     n = 0;
625     CreateCopy ( lines[n++],S );
626 
627     // count primitive operators
628     while (!f.FileEnd() && (strstr(S,"symop"))) {
629       npops++;
630       f.ReadLine ( S,sizeof(S) );
631       if (n>=nlines)  {
632         nlines += + 256;
633         GetVectorMemory ( lines1,nlines,0 );
634         for (i=0;i<n;i++)
635           lines1[i] = lines[i];
636         for (i=n;i<nlines;i++)
637           lines1[i] = NULL;
638         FreeVectorMemory ( lines,0 );
639         lines = lines1;
640       }
641       CreateCopy ( lines[n++],S );
642     }
643 
644     // count centering operators
645     while (!f.FileEnd() && (strstr(S,"cenop"))) {
646       ncops++;
647       f.ReadLine ( S,sizeof(S) );
648       if (n>=nlines)  {
649         nlines += + 256;
650         GetVectorMemory ( lines1,nlines,0 );
651         for (i=0;i<n;i++)
652           lines1[i] = lines[i];
653         for (i=n;i<nlines;i++)
654           lines1[i] = NULL;
655         FreeVectorMemory ( lines,0 );
656         lines = lines1;
657       }
658       CreateCopy ( lines[n++],S );
659     }
660 
661     Nops = npops*ncops;
662     symOp = new PSymOp[Nops];
663     RC    = SYMOP_Ok;
664 
665     n = 0;  // start second pass here
666 
667     // read primitive operators
668     for (i=0;(i<npops) && (!RC);i++)  {
669       symOp[i] = new SymOp();
670       RC = symOp[i]->SetSymOp ( lines[n++]+6 );
671     }
672 
673     // loop over non-trivial centering operators, and for each loop
674     // over primtive operators
675     for (i=1;(i<ncops) && (!RC);i++)  {
676       n++;  // this also skips the identity operator
677       for (j=0;(j<npops) && (!RC);j++)  {
678         symOp[i*npops+j] = new SymOp();
679         RC = symOp[i*npops+j]->SetSymOp ( lines[n]+6 );
680         symOp[i*npops+j]->GetTMatrix(T1);
681         symOp[j]->GetTMatrix(T2);
682         for (k=0;k<4;k++)
683           for (l=0;l<4;l++) {
684             T3[k][l] = 0.0;
685             for (m=0;m<4;m++)
686               T3[k][l] += T1[k][m]*T2[m][l];
687           }
688         for (k=0;k<3;k++)                  // kdc fix
689           T3[k][3] -= floor ( T3[k][3] );  // kdc fix
690         symOp[i*npops+j]->CompileOpTitle ( O,T3,false );
691         symOp[i*npops+j]->SetSymOp ( O );
692       }
693     }
694 
695     f.shut();
696 
697     for (i=0;i<nlines;i++)
698       if (lines[i])  delete[] lines[i];
699     FreeVectorMemory ( lines,0 );
700 
701     return RC;
702 
703   }
704 
705   /*
706   int  SymOps::SetGroup ( cpstr SpaceGroup,
707                            cpstr syminfo_lib )  {
708   CFile  f;
709   pstr   p;
710   char   S[500];
711   char   G[100];
712   char   O[100];
713   mat44  T1,T2,T3;
714   int    i,j,k,l,m,RC;
715   int    npops,ncops;
716   long   symop_start;
717   //fpos_t symop_start;
718 
719     FreeMemory();
720 
721     npops = 0;
722     ncops = 0;
723 
724     CreateCopy ( SpGroup,SpaceGroup );
725 
726     if (!syminfo_lib)          p = pstr(syminfo_file);
727     else if (!syminfo_lib[0])  p = pstr(syminfo_file);
728                          else  p = pstr(syminfo_lib);
729     f.assign ( p,true );
730     if (!f.reset(true))  {
731       p = getenv ( "SYMINFO" );
732       if (p)
733         strcpy ( S,p );
734       else  {
735         p = getenv ( "CLIBD" );
736         if (p)  {
737           strcpy ( S,p );
738           if (S[strlen(S)-1]!='/')  strcat ( S,"/" );
739           strcat ( S,"syminfo.lib" );
740         } else
741           strcpy ( S,"syminfo.lib" );
742       }
743       f.assign ( S,true );
744       if (!f.reset(true))  return SYMOP_NoLibFile;
745     }
746 
747     if (strncasecmp(SpGroup,"Hall:",5))  {
748       // normal space group symbol on input
749       strcpy ( G," '"    );
750       strcat ( G,SpGroup );
751       strcat ( G,"'"     );
752       S[0] = char(0);
753       while (!f.FileEnd() && !(strstr(S,G) && (strstr(S,"symbol xHM") ||
754               strstr(S,"symbol old"))))
755         f.ReadLine ( S,sizeof(S) );
756     } else  {
757       // hall descriptor on input
758       strcpy ( G," ' " );
759       p = &(SpGroup[5]);
760       while (*p==' ')  p++;
761       strcat ( G,p     );
762       strcat ( G,"'"   );
763       S[0] = char(0);
764       while (!f.FileEnd() && !(strstr(S,G) && strstr(S,"symbol Hall")))
765         f.ReadLine ( S,sizeof(S) );
766     }
767     if (f.FileEnd())  {
768       f.shut();
769       return SYMOP_UnknownSpaceGroup;
770     }
771 
772     // found spacegroup, move to symop lines
773     while (!f.FileEnd() && (!strstr(S,"symop"))) {
774       symop_start = f.Position();
775   //#    fgetpos ( f.GetHandle(),&symop_start );
776       f.ReadLine ( S,sizeof(S) );
777     }
778     // count primitive operators
779     while (!f.FileEnd() && (strstr(S,"symop"))) {
780       npops++;
781       f.ReadLine ( S,sizeof(S) );
782     }
783     // count centering operators
784     while (!f.FileEnd() && (strstr(S,"cenop"))) {
785       ncops++;
786       f.ReadLine ( S,sizeof(S) );
787     }
788     Nops = npops*ncops;
789     f.seek(symop_start);
790   //#  fsetpos ( f.GetHandle(),&symop_start );
791     SymOp = new PSymOp[Nops];
792     RC    = SYMOP_Ok;
793 
794     // read primitive operators
795     for (i=0;(i<npops) && (!RC);i++)  {
796       f.ReadLine ( S,sizeof(S) );
797       SymOp[i] = new SymOp();
798       RC = SymOp[i]->SetSymOp ( S+6 );
799     }
800 
801     // skip identity centering operator
802     f.ReadLine ( S,sizeof(S) );
803     // loop over non-trivial centering operators, and for each loop
804     // over primtive operators
805     for (i=1;(i<ncops) && (!RC);i++)  {
806       f.ReadLine ( S,sizeof(S) );
807       for (j=0;(j<npops) && (!RC);j++)  {
808         SymOp[i*npops+j] = new SymOp();
809         RC = SymOp[i*npops+j]->SetSymOp ( S+6 );
810 
811         SymOp[i*npops+j]->GetTMatrix(T1);
812         SymOp[j]->GetTMatrix(T2);
813         for (k=0;k<4;k++)
814           for (l=0;l<4;l++) {
815             T3[k][l] = 0.0;
816             for (m=0;m<4;m++)
817               T3[k][l] += T1[k][m]*T2[m][l];
818           }
819         for (k=0;k<3;k++)                  // kdc fix
820           T3[k][3] -= floor ( T3[k][3] );  // kdc fix
821         SymOp[i*npops+j]->CompileOpTitle(O,T3,false);
822         SymOp[i*npops+j]->SetSymOp (O);
823       }
824     }
825 
826     f.shut();
827 
828     return RC;
829 
830   }
831   */
832 
Reset()833   void SymOps::Reset()  {
834   // removes all symmetry operations
835     FreeMemory();
836   }
837 
AddSymOp(cpstr XYZOperation)838   int  SymOps::AddSymOp ( cpstr XYZOperation )  {
839   // adds a symmetry operation
840   PPSymOp symOp1;
841   int      i;
842     symOp1 = new PSymOp[Nops+1];
843     for (i=0;i<Nops;i++)
844       symOp1[i] = symOp[i];
845     if (symOp) delete[] symOp;
846     symOp = symOp1;
847     i = Nops;
848     symOp[i] = new SymOp();
849     Nops++;
850     return symOp[i]->SetSymOp ( XYZOperation );
851   }
852 
PutGroupName(cpstr SpGroupName)853   void SymOps::PutGroupName ( cpstr SpGroupName )  {
854     CreateCopy ( SpGroup,SpGroupName );
855   }
856 
857 
GetNofSymOps()858   int  SymOps::GetNofSymOps()  {
859   //  GetNofSymOps()  returns Nops -- the number of symmetry operations
860     return Nops;
861   }
862 
GetSymOp(int Nop)863   pstr SymOps::GetSymOp ( int Nop )  {
864     if ((0<=Nop) && (Nop<Nops))  return symOp[Nop]->GetSymOp();
865                            else  return pstr("");
866   }
867 
Transform(realtype & x,realtype & y,realtype & z,int Nop)868   int  SymOps::Transform ( realtype & x, realtype & y, realtype & z,
869                             int Nop )  {
870   //  Transform(..) transforms the coordinates according to the
871   // symmetry operation Nop. The return code is non-zero if
872   // Nop is a wrong operation number (must range from 0 to Nops-1).
873     if ((Nop<0) || (Nop>=Nops))  return 1;
874     if (symOp[Nop])  {
875       symOp[Nop]->Transform ( x,y,z );
876       return 0;
877     } else
878       return 2;
879   }
880 
GetTMatrix(mat44 & TMatrix,int Nop)881   int  SymOps::GetTMatrix ( mat44 & TMatrix, int Nop )  {
882   //  GetTMatrix(..) returns the coordinate transformation matrix
883   // for the symmetry operation Nop. The return code is non-zero if
884   // Nop is a wrong operation number (must range from 0 to Nops-1).
885     if ((Nop<0) || (Nop>=Nops))  return 1;
886     if (symOp[Nop])  {
887       symOp[Nop]->GetTMatrix ( TMatrix );
888       return 0;
889     } else
890       return 2;
891   }
892 
Print()893   void  SymOps::Print()  {
894   int  i;
895   char S[200];
896     printf ( "  SPACE GROUP  '%s'\n",SpGroup );
897     for (i=0;i<Nops;i++)  {
898       symOp[i]->Print();
899       if (symOp[i]->CompileOpTitle(S))
900             printf ( " CHECK STATUS: Ok\n" );
901       else  printf ( " CHECK STATUS: Generated '%s'\n",S );
902     }
903   }
904 
905 
Copy(PSymOps SymOps)906   void  SymOps::Copy ( PSymOps SymOps )  {
907   int i;
908     FreeMemory();
909     CreateCopy ( SpGroup,SymOps->SpGroup );
910     Nops = SymOps->Nops;
911     if (Nops>0)  {
912       symOp = new PSymOp[Nops];
913       for (i=0;i<Nops;i++) {
914         symOp[i] = new SymOp();
915         symOp[i]->Copy ( SymOps->symOp[i] );
916       }
917     }
918   }
919 
920 
write(io::RFile f)921   void  SymOps::write ( io::RFile f )  {
922   int  i;
923   byte Version=1;
924     f.WriteByte   ( &Version );
925     f.CreateWrite ( SpGroup  );
926     f.WriteInt    ( &Nops    );
927     for (i=0;i<Nops;i++)
928       StreamWrite ( f,symOp[i] );
929   }
930 
read(io::RFile f)931   void  SymOps::read ( io::RFile f )  {
932   int  i;
933   byte Version;
934     FreeMemory();
935     f.ReadByte   ( &Version );
936     f.CreateRead ( SpGroup  );
937     f.ReadInt    ( &Nops    );
938     if (Nops>0)  {
939       symOp = new PSymOp[Nops];
940       for (i=0;i<Nops;i++) {
941         symOp[i] = NULL;
942         StreamRead ( f,symOp[i] );
943       }
944     }
945   }
946 
947 
948   MakeStreamFunctions(SymOps)
949 
950 }  // namespace mmdb
951 
952 
953 /*
954 
955 void TestSymOps()  {
956 pstr     p,p1;
957 int      RC;
958 char     S[500];
959 SymOps  SymOps;
960 CFile    f;
961 
962   p = getenv ( "PDB_ROOT" );
963   if (p)  {
964     strcpy ( S,p );
965     strcat ( S,"/lib/" );
966   } else
967     S[0] = char(0);
968   strcat   ( S,"symop.lib" );
969   f.assign ( S,true );
970   if (!f.reset())  {
971     printf ( " +++++ No symop.lib file found.\n" );
972     return;
973   }
974 
975   while (!f.FileEnd())  {
976     f.ReadLine ( S,sizeof(S) );
977     if (S[0] && (S[0]!=' '))  {
978       p = FirstOccurence ( S,'\'' );
979       if (p)  {
980         p++;
981         p1 = FirstOccurence ( p,'\'' );
982         if (!p1)  p = NULL;
983       }
984       if (!p)  {
985         printf ( " +++++ Strange line in symop.lib:\n"
986                  "%s\n",S );
987         return;
988       }
989       *p1 = char(0);
990       RC = SymOps.SetGroup ( p );
991       printf ( " =========================================================\n"
992                "  RC=%i\n",RC );
993       SymOps.Print();
994     }
995   }
996 
997   return;
998 
999 }
1000 
1001 */
1002