1
2 /*
3 A* -------------------------------------------------------------------
4 B* This file contains source code for the PyMOL computer program
5 C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific.
6 D* -------------------------------------------------------------------
7 E* It is unlawful to modify or remove this copyright notice.
8 F* -------------------------------------------------------------------
9 G* Please see the accompanying LICENSE file for further information.
10 H* -------------------------------------------------------------------
11 I* Additional authors of this source file include:
12 -*
13 -*
14 -*
15 Z* -------------------------------------------------------------------
16 */
17 #include"os_python.h"
18
19 #include"os_predef.h"
20 #include"os_std.h"
21 #include"os_gl.h"
22 #include"Base.h"
23 #include"OOMac.h"
24 #include"CGO.h"
25
26 #include"Map.h"
27
28 #include"Shaker.h"
29 #include"vla.h"
30 #include<memory>
31
CShaker(PyMOLGlobals * G)32 CShaker::CShaker(PyMOLGlobals * G)
33 {
34 this->G = G;
35 this->DistCon = pymol::vla<ShakerDistCon>(1000);
36 this->PyraCon = pymol::vla<ShakerPyraCon>(1000);
37 this->PlanCon = pymol::vla<ShakerPlanCon>(1000);
38 this->TorsCon = pymol::vla<ShakerTorsCon>(1000);
39 this->LineCon = pymol::vla<ShakerLineCon>(100);
40 this->NDistCon = 0;
41 this->NPyraCon = 0;
42 this->NPlanCon = 0;
43 this->NLineCon = 0;
44 this->NTorsCon = 0;
45 }
46
ShakerReset(CShaker * I)47 void ShakerReset(CShaker * I)
48 {
49 I->NDistCon = 0;
50 I->NPyraCon = 0;
51 I->NPlanCon = 0;
52 I->NLineCon = 0;
53 I->NTorsCon = 0;
54 }
55
ShakerAddDistCon(CShaker * I,int atom0,int atom1,float target,int type,float wt)56 void ShakerAddDistCon(CShaker * I, int atom0, int atom1, float target, int type, float wt)
57 {
58 ShakerDistCon *sdc;
59
60 VLACheck(I->DistCon, ShakerDistCon, I->NDistCon);
61 sdc = I->DistCon + I->NDistCon;
62 sdc->at0 = atom0;
63 sdc->at1 = atom1;
64 sdc->targ = target;
65 sdc->type = type;
66 sdc->weight = wt;
67 I->NDistCon++;
68 }
69
ShakerGetPyra(float * targ2,const float * v0,const float * v1,const float * v2,const float * v3)70 float ShakerGetPyra(float *targ2,
71 const float *v0, const float *v1, const float *v2, const float *v3)
72 {
73 float d0[3], cp[3], d2[3], d3[3];
74 float av[3], t0[3];
75
76 add3f(v1, v2, av);
77 subtract3f(v2, v1, d2);
78 add3f(v3, av, av);
79 subtract3f(v3, v1, d3);
80 subtract3f(av, v0, t0);
81 cross_product3f(d2, d3, cp);
82 scale3f(av, 0.33333333F, av);
83 normalize3f(cp);
84 subtract3f(av, v0, d0);
85
86 (*targ2) = length3f(d0);
87 return (dot_product3f(d0, cp));
88 }
89
ShakerDoPyra(float targ1,float targ2,const float * v0,const float * v1,const float * v2,const float * v3,float * p0,float * p1,float * p2,float * p3,float wt,float inv_wt)90 float ShakerDoPyra(float targ1, float targ2,
91 const float *v0, const float *v1, const float *v2, const float *v3,
92 float *p0, float *p1, float *p2, float *p3, float wt, float inv_wt)
93 {
94 float d0[3], cp[3], d2[3], d3[3];
95 float av[3], t0[3], push[3];
96
97 float cur, dev, sc, result1, result2 = 0.0F;
98
99 add3f(v1, v2, av);
100 subtract3f(v2, v1, d2);
101 add3f(v3, av, av);
102 subtract3f(v3, v1, d3);
103 subtract3f(av, v0, t0);
104 cross_product3f(d2, d3, cp);
105 scale3f(av, 0.33333333F, av);
106 normalize3f(cp);
107 subtract3f(av, v0, d0);
108
109 cur = dot_product3f(d0, cp);
110 dev = cur - targ1;
111 result1 = (float) fabs(dev);
112 if(result1 > R_SMALL8) {
113 sc = wt * dev;
114 if((cur * targ1) < 0.0) /* inverted */
115 sc = sc * inv_wt; /* inversion fixing weight */
116 scale3f(cp, sc, push);
117 add3f(push, p0, p0);
118 scale3f(push, 0.333333F, push);
119 subtract3f(p1, push, p1);
120 subtract3f(p2, push, p2);
121 subtract3f(p3, push, p3);
122 }
123
124 if((targ2 >= 0.0F) && ((cur * targ1 > 0.0) || (fabs(targ1) < 0.1))) {
125 /* so long as we're not inverted...
126 also make sure v0 is the right distance from the average point */
127 cur = length3f(d0);
128 normalize3f(d0);
129 dev = cur - targ2;
130 result2 = (float) fabs(dev);
131 if(result2 > R_SMALL4) {
132 sc = wt * dev * 2.0F;
133 scale3f(d0, sc, push);
134 add3f(push, p0, p0);
135 scale3f(push, 0.333333F, push);
136 subtract3f(p1, push, p1);
137 subtract3f(p2, push, p2);
138 subtract3f(p3, push, p3);
139 }
140 }
141
142 return result1 + result2;
143 }
144
ShakerDoLine(const float * v0,const float * v1,const float * v2,float * p0,float * p1,float * p2,float wt)145 float ShakerDoLine(const float *v0, const float *v1, const float *v2,
146 float *p0, float *p1, float *p2, float wt)
147 {
148 /* v0-v1-v2 */
149
150 float d0[3], d1[3], cp[3], d2[3], d3[3], d4[3], push[3];
151 float dev, sc, lcp, result;
152
153 subtract3f(v2, v1, d2);
154 subtract3f(v0, v1, d1);
155 normalize3f(d2);
156 normalize23f(d1, d0);
157
158 cross_product3f(d2, d0, cp);
159 lcp = (float) length3f(cp);
160 if(lcp > R_SMALL4) {
161 lcp = 1.0F / lcp;
162 scale3f(cp, lcp, cp); /* axis 0 */
163
164 subtract3f(v2, v0, d3);
165 normalize3f(d3); /* axis 1 */
166
167 cross_product3f(cp, d3, d4);
168 normalize3f(d4); /* displacement direction */
169
170 dev = dot_product3f(d1, d4); /* current deviation */
171
172 if((result = (float) fabs(dev)) > R_SMALL8) {
173 sc = wt * dev;
174 scale3f(d4, sc, push);
175 add3f(push, p1, p1);
176 scale3f(push, 0.5F, push);
177 subtract3f(p0, push, p0);
178 subtract3f(p2, push, p2);
179 } else {
180 result = 0.0;
181 }
182 } else
183 result = 0.0;
184 return result;
185
186 }
187
ShakerDoPlan(const float * v0,const float * v1,const float * v2,const float * v3,float * p0,float * p1,float * p2,float * p3,float target,int fixed,float wt)188 float ShakerDoPlan(
189 const float *v0,
190 const float *v1,
191 const float *v2,
192 const float *v3,
193 float *p0, float *p1, float *p2, float *p3,
194 float target, int fixed, float wt)
195 {
196
197 float result;
198
199 float d01[3], d12[3], d23[3], d03[3], cp0[3], cp1[3], dp, sc, dev, d0[3], push[3];
200 double s01, s12, s23, s03;
201
202 subtract3f(v0, v1, d01);
203 subtract3f(v1, v2, d12);
204 subtract3f(v2, v3, d23);
205 subtract3f(v0, v3, d03);
206
207 s03 = lengthsq3f(d03);
208 s01 = lengthsq3f(d01);
209 s12 = lengthsq3f(d12);
210 s23 = lengthsq3f(d23);
211
212 if((s03 < s01) || (s03 < s12) || (s03 < s23))
213 return 0.0F;
214
215 cross_product3f(d01, d12, cp0);
216 cross_product3f(d12, d23, cp1);
217
218 normalize3f(cp0);
219 normalize3f(cp1);
220
221 dp = dot_product3f(cp0, cp1);
222
223 result = (dev = 1.0F - (float) fabs(dp));
224
225 if(dev > R_SMALL4) {
226
227 /*
228 add3f(cp0,cp1,d0);
229 normalize3f(d0);
230
231 cross_product3f(cp0,d12,pos);
232 dp2 = dot_product3f(cp1,pos);
233 */
234
235 if(fixed && (dp * target < 0.0F)) {
236
237 /* fixed & backwards... */
238
239 if(dp < 0.0F) {
240 sc = -wt * dev * 0.5F;
241 } else {
242 sc = wt * dev * 0.5F;
243 }
244 sc *= 0.02F; /* weaken considerably to allow resolution of
245 inconsistencies (folded rings, etc.) */
246
247 } else if(dp > 0) {
248 sc = -wt * dev * 0.5F;
249 } else {
250 sc = wt * dev * 0.5F;
251 }
252
253 if(fixed && (fixed < 7)) {
254 /* in small rings, ramp up the planarity factor */
255 sc *= 8;
256 } else {
257 sc *= 0.2F;
258 }
259
260 /* pair-wise nudges */
261
262 subtract3f(v0, v3, d0);
263 normalize3f(d0);
264 scale3f(d0, sc, push);
265 add3f(push, p0, p0);
266 subtract3f(p3, push, p3);
267
268 subtract3f(v1, v2, d0);
269 normalize3f(d0);
270 scale3f(d0, sc, push);
271 add3f(push, p1, p1);
272 subtract3f(p2, push, p2);
273
274 sc = -sc;
275 subtract3f(v0, v2, d0);
276 normalize3f(d0);
277 scale3f(d0, sc, push);
278 add3f(push, p0, p0);
279 subtract3f(p2, push, p2);
280
281 subtract3f(v1, v3, d0);
282 normalize3f(d0);
283 scale3f(d0, sc, push);
284 add3f(push, p1, p1);
285 subtract3f(p3, push, p3);
286
287 } else {
288 result = 0.0;
289 }
290 return result;
291
292 }
293
ShakerAddPyraCon(CShaker * I,int atom0,int atom1,int atom2,int atom3,float targ1,float targ2)294 void ShakerAddPyraCon(CShaker * I, int atom0, int atom1, int atom2, int atom3,
295 float targ1, float targ2)
296 {
297 ShakerPyraCon *spc;
298
299 VLACheck(I->PyraCon, ShakerPyraCon, I->NPyraCon);
300 spc = I->PyraCon + I->NPyraCon;
301 spc->at0 = atom0;
302 spc->at1 = atom1;
303 spc->at2 = atom2;
304 spc->at3 = atom3;
305 spc->targ1 = targ1;
306 spc->targ2 = targ2;
307 I->NPyraCon++;
308 }
309
ShakerAddTorsCon(CShaker * I,int atom0,int atom1,int atom2,int atom3,int type)310 void ShakerAddTorsCon(CShaker * I, int atom0, int atom1, int atom2, int atom3, int type)
311 {
312 ShakerTorsCon *stc;
313
314 VLACheck(I->TorsCon, ShakerTorsCon, I->NTorsCon);
315 stc = I->TorsCon + I->NTorsCon;
316 stc->at0 = atom0;
317 stc->at1 = atom1;
318 stc->at2 = atom2;
319 stc->at3 = atom3;
320 stc->type = type;
321 I->NTorsCon++;
322
323 }
324
ShakerAddPlanCon(CShaker * I,int atom0,int atom1,int atom2,int atom3,float target,int fixed)325 void ShakerAddPlanCon(CShaker * I, int atom0, int atom1, int atom2, int atom3,
326 float target, int fixed)
327 {
328 ShakerPlanCon *spc;
329
330 VLACheck(I->PlanCon, ShakerPlanCon, I->NPlanCon);
331 spc = I->PlanCon + I->NPlanCon;
332 spc->at0 = atom0;
333 spc->at1 = atom1;
334 spc->at2 = atom2;
335 spc->at3 = atom3;
336 spc->fixed = fixed;
337 spc->target = target;
338 I->NPlanCon++;
339 }
340
ShakerAddLineCon(CShaker * I,int atom0,int atom1,int atom2)341 void ShakerAddLineCon(CShaker * I, int atom0, int atom1, int atom2)
342 {
343 ShakerLineCon *slc;
344
345 VLACheck(I->LineCon, ShakerLineCon, I->NLineCon);
346 slc = I->LineCon + I->NLineCon;
347 slc->at0 = atom0;
348 slc->at1 = atom1;
349 slc->at2 = atom2;
350 I->NLineCon++;
351 }
352
353