1
2
3 /*
4 A* -------------------------------------------------------------------
5 B* This file contains source code for the PyMOL computer program
6 C* Copyright (c) Schrodinger, LLC.
7 D* -------------------------------------------------------------------
8 E* It is unlawful to modify or remove this copyright notice.
9 F* -------------------------------------------------------------------
10 G* Please see the accompanying LICENSE file for further information.
11 H* -------------------------------------------------------------------
12 I* Additional authors of this source file include:
13 -*
14 -*
15 -*
16 Z* -------------------------------------------------------------------
17 */
18 #include"os_python.h"
19 #include"os_predef.h"
20 #include"os_std.h"
21 #include"os_gl.h"
22 #include"MemoryDebug.h"
23 #include"CGO.h"
24 #include"Err.h"
25 #include"Base.h"
26 #include"OOMac.h"
27 #include"Crystal.h"
28 #include"Feedback.h"
29 #include"Util.h"
30 #include"PConv.h"
31
CrystalAsPyList(CCrystal * I)32 PyObject *CrystalAsPyList(CCrystal * I)
33 {
34 PyObject *result = NULL;
35
36 if(I) {
37 result = PyList_New(2);
38 PyList_SetItem(result, 0, PConvFloatArrayToPyList(I->Dim, 3));
39 PyList_SetItem(result, 1, PConvFloatArrayToPyList(I->Angle, 3));
40 }
41 return (PConvAutoNone(result));
42 }
43
CrystalFromPyList(CCrystal * I,PyObject * list)44 int CrystalFromPyList(CCrystal * I, PyObject * list)
45 {
46 int ok = true, rok = true;
47 int ll = 0;
48 if(ok)
49 ok = (I != NULL);
50 if(ok)
51 ok = PyList_Check(list);
52 if(ok)
53 ll = PyList_Size(list);
54 rok = ok;
55 if(ok && (ll > 0))
56 ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list, 0), I->Dim, 3);
57 if(ok && (ll > 1))
58 ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list, 1), I->Angle, 3);
59 if(ok)
60 CrystalUpdate(I);
61
62 /* TO SUPPORT BACKWARDS COMPATIBILITY...
63 Always check ll when adding new PyList_GetItem's */
64
65 return (rok);
66 }
67
CCrystal(PyMOLGlobals * GParam)68 CCrystal::CCrystal(PyMOLGlobals* GParam) : G(GParam)
69 {
70 identity33f(RealToFrac);
71 identity33f(FracToReal);
72
73 }
74
75
CrystalUpdate(CCrystal * I)76 void CrystalUpdate(CCrystal * I)
77 {
78 float cabg[3]; /* Cosine of axis angle */
79 float sabg[3]; /* Singe of axis angle */
80 float cabgs[3];
81 float sabgs1;
82 int i;
83
84 /* if we just cleared out the memory, but didn't init
85 * then init the crystal and return */
86 if (!I->Dim[0] || !I->Dim[1] || !I->Dim[2]) {
87 *I = CCrystal(I->G);
88 return;
89 }
90
91 if (!I->Angle[0] || !I->Angle[1] || !I->Angle[2]) {
92 I->Angle[0] = I->Angle[1] = I->Angle[2] = 90.0F;
93 }
94
95 for(i = 0; i < 9; i++) {
96 I->RealToFrac[i] = 0.0;
97 I->FracToReal[i] = 0.0;
98 }
99
100 for(i = 0; i < 3; i++) {
101 cabg[i] = (float) cos(I->Angle[i] * PI / 180.0);
102 sabg[i] = (float) sin(I->Angle[i] * PI / 180.0);
103 }
104
105 cabgs[0] = (cabg[1] * cabg[2] - cabg[0]) / (sabg[1] * sabg[2]);
106 cabgs[1] = (cabg[2] * cabg[0] - cabg[1]) / (sabg[2] * sabg[0]);
107 cabgs[2] = (cabg[0] * cabg[1] - cabg[2]) / (sabg[0] * sabg[1]);
108
109 I->UnitCellVolume = (float) (I->Dim[0] * I->Dim[1] * I->Dim[2] *
110 sqrt1d(1.0 + (double) 2.0 * cabg[0] * cabg[1] * cabg[2] -
111 (double) (cabg[0] * cabg[0] +
112 (double) cabg[1] * cabg[1] +
113 (double) cabg[2] * cabg[2])));
114
115 I->RecipDim[0] = I->Dim[1] * I->Dim[2] * sabg[0] / I->UnitCellVolume;
116 I->RecipDim[1] = I->Dim[0] * I->Dim[2] * sabg[1] / I->UnitCellVolume;
117 I->RecipDim[2] = I->Dim[0] * I->Dim[1] * sabg[2] / I->UnitCellVolume;
118
119 sabgs1 = (float) sqrt1d(1.0 - cabgs[0] * cabgs[0]);
120
121 I->RealToFrac[0] = 1.0F / I->Dim[0];
122 I->RealToFrac[1] = -cabg[2] / (sabg[2] * I->Dim[0]);
123 I->RealToFrac[2] = -(cabg[2] * sabg[1] * cabgs[0] + cabg[1] * sabg[2]) /
124 (sabg[1] * sabgs1 * sabg[2] * I->Dim[0]);
125 I->RealToFrac[4] = 1.0F / (sabg[2] * I->Dim[1]);
126 I->RealToFrac[5] = cabgs[0] / (sabgs1 * sabg[2] * I->Dim[1]);
127 I->RealToFrac[8] = 1.0F / (sabg[1] * sabgs1 * I->Dim[2]);
128
129 I->FracToReal[0] = I->Dim[0];
130 I->FracToReal[1] = cabg[2] * I->Dim[1];
131 I->FracToReal[2] = cabg[1] * I->Dim[2];
132 I->FracToReal[4] = sabg[2] * I->Dim[1];
133 I->FracToReal[5] = -sabg[1] * cabgs[0] * I->Dim[2];
134 I->FracToReal[8] = sabg[1] * sabgs1 * I->Dim[2];
135
136 I->Norm[0] = (float) sqrt1d(I->RealToFrac[0] * I->RealToFrac[0] +
137 I->RealToFrac[1] * I->RealToFrac[1] +
138 I->RealToFrac[2] * I->RealToFrac[2]);
139 I->Norm[1] = (float) sqrt1d(I->RealToFrac[3] * I->RealToFrac[3] +
140 I->RealToFrac[4] * I->RealToFrac[4] +
141 I->RealToFrac[5] * I->RealToFrac[5]);
142 I->Norm[2] = (float) sqrt1d(I->RealToFrac[6] * I->RealToFrac[6] +
143 I->RealToFrac[7] * I->RealToFrac[7] +
144 I->RealToFrac[8] * I->RealToFrac[8]);
145
146 }
147
CrystalDump(CCrystal * I)148 void CrystalDump(CCrystal * I)
149 {
150 PyMOLGlobals *G = I->G;
151 int i;
152
153 PRINTF
154 " Crystal: Unit Cell %8.3f %8.3f %8.3f\n", I->Dim[0], I->Dim[1], I->Dim[2]
155 ENDF(G);
156 PRINTF
157 " Crystal: Alpha Beta Gamma %8.3f %8.3f %8.3f\n",
158 I->Angle[0], I->Angle[1], I->Angle[2]
159 ENDF(G);
160 PRINTF " Crystal: RealToFrac Matrix\n" ENDF(G);
161 for(i = 0; i < 3; i++) {
162 PRINTF " Crystal: %9.4f %9.4f %9.4f\n",
163 I->RealToFrac[i * 3], I->RealToFrac[i * 3 + 1], I->RealToFrac[i * 3 + 2]
164 ENDF(G);
165 }
166 PRINTF " Crystal: FracToReal Matrix\n" ENDF(G);
167 for(i = 0; i < 3; i++) {
168 PRINTF
169 " Crystal: %9.4f %9.4f %9.4f\n",
170 I->FracToReal[i * 3], I->FracToReal[i * 3 + 1], I->FracToReal[i * 3 + 2]
171 ENDF(G);
172 }
173 PRINTF " Crystal: Unit Cell Volume %8.0f.\n", I->UnitCellVolume ENDF(G);
174
175 }
176
177 static float unitCellVertices[][3] = { {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0}, // bottom 4 vertices
178 {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1} }; // top 4 vertices
179 static int unitCellLineIndices[] = { 0, 1, 1, 2, 2, 3, 3, 0, // bottom 4 lines
180 4, 5, 5, 6, 6, 7, 7, 4, // top 4 lines
181 0, 4, 1, 5, 2, 6, 3, 7 }; // 4 connector lines
182
CrystalGetUnitCellCGO(CCrystal * I)183 CGO *CrystalGetUnitCellCGO(CCrystal * I)
184 {
185 PyMOLGlobals *G = I->G;
186 float v[3];
187 CGO *cgo = NULL;
188 if(I) {
189 cgo = CGONew(G);
190 CGODisable(cgo, GL_LIGHTING);
191
192 float *vertexVals = CGODrawArrays(cgo, GL_LINES, CGO_VERTEX_ARRAY, 24);
193 for (int pl = 0; pl < 24 ; pl++){
194 transform33f3f(I->FracToReal, unitCellVertices[unitCellLineIndices[pl]], v);
195 copy3f(v, &vertexVals[pl*3]);
196 }
197
198 CGOEnable(cgo, GL_LIGHTING);
199 CGOStop(cgo);
200 }
201 return (cgo);
202 }
203