ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/eangle.c
Revision: 105
Committed: Fri Feb 20 21:01:21 2009 UTC (12 years, 4 months ago) by gilbertke
File size: 26517 byte(s)
Log Message:
add job_contro.h
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4
5 #include "utility.h"
6
7 //double acos(double);
8
9
10 EXTERN struct t_minim_values {
11 int iprint, ndc, nconst;
12 float dielc;
13 } minim_values;
14
15 void eangle(int nang,int **,int *use,double *x,double *y,double *z,int *angin,int *angtype,float *anat,float *acon,double *ebend);
16 void eangle1(int nang,int natom,int **,int *use,double *x,double *y,double *z,int *angin,int *angtype,float *anat,float *acon,double *ebend,double **dea);
17 void eangle2(int i,int nang,int **,int *angin, double *x,double *y,double *z,float *acon,float *anat,double **dea,float **hessx,float **hessy,float **hessz);
18 void eangle2a(int iatom,int nang,int **,int *angin,double *x,double *y,double *z,float *acon,float *anat,float **hessx,float **hessy,float **hessz);
19 void eangle2b(int i,int nang,int **,int *angin,double *x,double *y,double *z,float *acon,float *anat,double **dea);
20 // ======================
21 void eangle(int nang,int **i13,int *use,double *x,double *y,double *z,int *angin,int *angtype,float *anat,float *acon,double *ebend)
22 {
23 int i, ia, ib, ic, id;
24 double e, angle, dot, cosine;
25 double dt, dt2, dt3, dt4;
26 double xia, yia, zia;
27 double xib, yib, zib;
28 double xic, yic, zic;
29 double xab, yab, zab, rab2;
30 double xcb, ycb, zcb, rcb2;
31
32
33 *ebend = 0.0;
34 if (minim_values.iprint)
35 {
36 fprintf(pcmlogfile,"\nAngle Terms \n");
37 fprintf(pcmlogfile," At1 At2 At3 Angle Thet0 Tconst Ebend\n");
38 }
39
40 for (i=0; i < nang; i++)
41 {
42 ia = i13[i][0];
43 ib = i13[i][1];
44 ic = i13[i][2];
45 id = i13[i][3];
46 if (use[ia] || use[ib] || use[ic])
47 {
48 xia = x[ia];
49 yia = y[ia];
50 zia = z[ia];
51 xib = x[ib];
52 yib = y[ib];
53 zib = z[ib];
54 xic = x[ic];
55 yic = y[ic];
56 zic = z[ic];
57 xab = xia - xib;
58 yab = yia - yib;
59 zab = zia - zib;
60 rab2 = xab*xab + yab*yab + zab*zab;
61 xcb = xic - xib;
62 ycb = yic - yib;
63 zcb = zic - zib;
64 rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
65 if (rab2 != 0.0 && rcb2 != 0.00)
66 {
67 dot = xab*xcb + yab*ycb + zab*zcb;
68 cosine = dot / sqrt(rab2*rcb2);
69 if (cosine > 1.0)
70 cosine = 1.0;
71 if (cosine < -1.0)
72 cosine = -1.0;
73 angle = radian*acos(cosine);
74
75 dt = angle - anat[i];
76 dt2 = dt*dt;
77 dt3 = dt2*dt;
78 dt4 = dt2*dt2;
79 e = units.angunit*acon[i]*dt2
80 *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
81
82 *ebend += e;
83 if (minim_values.iprint)
84 {
85 fprintf(pcmlogfile,"Angle 1-3: (%-3d)- (%-3d)- (%-3d) %-8.3f %-8.3f %-8.4f = %-8.4f\n",
86 ia, ib, ic, angle, anat[i],acon[i], e);
87 }
88 }
89 }
90 }
91 }
92 // ===============================================
93 void eangle1(int nang,int natom,int **i13,int *use,double *x,double *y,double *z,int *angin,int *angtype,float *anat,float *acon,double *ebend,double **dea)
94 {
95 int i, ia, ib, ic, id;
96 double e, angle, dot, cosine;
97 double dt, dt2, dt3, dt4;
98 double deddt,terma,termc;
99 double xia, yia, zia;
100 double xib, yib, zib;
101 double xic, yic, zic;
102 double xab, yab, zab, rab2;
103 double xcb, ycb, zcb, rcb2;
104 double xp,yp,zp,rp;
105 double dedxia,dedyia,dedzia;
106 double dedxib,dedyib,dedzib;
107 double dedxic,dedyic,dedzic;
108 double dtemp;
109
110 *ebend = 0.0F;
111 for (i=0; i <= natom; i++)
112 {
113 dea[i][0] = 0.0;
114 dea[i][1] = 0.0;
115 dea[i][2] = 0.0;
116 }
117 for (i=0; i < nang; i++)
118 {
119 ia = i13[i][0];
120 ib = i13[i][1];
121 ic = i13[i][2];
122 id = i13[i][3];
123 if (use[ia] || use[ib] || use[ic])
124 {
125 xia = x[ia];
126 yia = y[ia];
127 zia = z[ia];
128 xib = x[ib];
129 yib = y[ib];
130 zib = z[ib];
131 xic = x[ic];
132 yic = y[ic];
133 zic = z[ic];
134 xab = xia - xib;
135 yab = yia - yib;
136 zab = zia - zib;
137 rab2 = xab*xab + yab*yab + zab*zab;
138 xcb = xic - xib;
139 ycb = yic - yib;
140 zcb = zic - zib;
141 rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
142 xp = ycb*zab - zcb*yab;
143 yp = zcb*xab - xcb*zab;
144 zp = xcb*yab - ycb*xab;
145 rp = sqrt(xp*xp + yp*yp + zp*zp);
146 if (rp != 0.0 )
147 {
148 dot = xab*xcb + yab*ycb + zab*zcb;
149 cosine = dot / sqrt(rab2*rcb2);
150 if (cosine > 1.0)
151 cosine = 1.0;
152 if (cosine < -1.0)
153 cosine = -1.0;
154 angle = radian*acos(cosine);
155 dt = angle - anat[i];
156 dt2 = dt*dt;
157 dt3 = dt2*dt;
158 dt4 = dt2*dt2;
159
160 e = units.angunit*acon[i]*dt2
161 *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
162 deddt = units.angunit*acon[i]* dt * radian
163 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
164 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
165
166
167 dtemp = dot*dot/(rab2*rcb2);
168 if (fabs(dtemp) >= 1.0)
169 dtemp = 0.9999*dtemp/(fabs(dtemp));
170 terma = sqrt(1.0-dtemp);
171 termc = sqrt(rab2*rcb2);
172 dedxia = -deddt*( xcb/termc - (xab*cosine)/rab2)/terma;
173 dedyia = -deddt*( ycb/termc - (yab*cosine)/rab2)/terma;
174 dedzia = -deddt*( zcb/termc - (zab*cosine)/rab2)/terma;
175
176 dedxic = -deddt*( xab/termc - (xcb*cosine)/rcb2)/terma;
177 dedyic = -deddt*( yab/termc - (ycb*cosine)/rcb2)/terma;
178 dedzic = -deddt*( zab/termc - (zcb*cosine)/rcb2)/terma;
179
180 dedxib = -deddt*( (-xab-xcb)/termc + (dot*(xcb*rab2 + xab*rcb2))/(rab2*rcb2*termc))/terma;
181 dedyib = -deddt*( (-yab-ycb)/termc + (dot*(ycb*rab2 + yab*rcb2))/(rab2*rcb2*termc))/terma;
182 dedzib = -deddt*( (-zab-zcb)/termc + (dot*(zcb*rab2 + zab*rcb2))/(rab2*rcb2*termc))/terma;
183
184
185 dea[ia][0] += dedxia;
186 dea[ia][1] += dedyia;
187 dea[ia][2] += dedzia;
188
189 dea[ib][0] += dedxib;
190 dea[ib][1] += dedyib;
191 dea[ib][2] += dedzib;
192
193 dea[ic][0] += dedxic;
194 dea[ic][1] += dedyic;
195 dea[ic][2] += dedzic;
196
197 *ebend += e;
198 }
199 }
200 }
201 }
202 // =========================
203 void eangle2(int i,int nang,int **i13,int *angin, double *x,double *y,double *z,float *acon,float *anat,double **dea,float **hessx,float **hessy,float **hessz)
204 {
205 eangle2a(i,nang,i13,angin,x,y,z,acon,anat,hessx,hessy,hessz);
206 }
207 // =================================================
208 void eangle2a(int iatom,int nang,int **i13,int *angin,double *x,double *y,double *z,float *acon,float *anat,float **hessx,float **hessy,float **hessz)
209 {
210 int i,ia,ib,ic;
211 double angle,dot,cosine;
212 double dt,dt2,dt3,dt4;
213 double xia,yia,zia;
214 double xib,yib,zib;
215 double xic,yic,zic;
216 double xab,yab,zab,rab;
217 double xcb,ycb,zcb,rcb;
218 double xp,yp,zp,rp,rp2;
219 double xrab,yrab,zrab,rab2;
220 double xrcb,yrcb,zrcb,rcb2;
221 double xabp,yabp,zabp;
222 double xcbp,ycbp,zcbp;
223 double deddt,d2eddt2,terma,termc;
224 double ddtdxia,ddtdyia,ddtdzia;
225 double ddtdxib,ddtdyib,ddtdzib;
226 double ddtdxic,ddtdyic,ddtdzic;
227 double dxiaxia,dxiayia,dxiazia;
228 double dxibxib,dxibyib,dxibzib;
229 double dxicxic,dxicyic,dxiczic;
230 double dyiayia,dyiazia,dziazia;
231 double dyibyib,dyibzib,dzibzib;
232 double dyicyic,dyiczic,dziczic;
233 double dxibxia,dxibyia,dxibzia;
234 double dyibxia,dyibyia,dyibzia;
235 double dzibxia,dzibyia,dzibzia;
236 double dxibxic,dxibyic,dxibzic;
237 double dyibxic,dyibyic,dyibzic;
238 double dzibxic,dzibyic,dzibzic;
239 double dxiaxic,dxiayic,dxiazic;
240 double dyiaxic,dyiayic,dyiazic;
241 double dziaxic,dziayic,dziazic;
242
243 for (i=0; i < nang; i++)
244 {
245 ia = i13[i][0];
246 ib = i13[i][1];
247 ic = i13[i][2];
248 if (iatom == ia || iatom == ib || iatom == ic)
249 {
250 xia = x[ia];
251 yia = y[ia];
252 zia = z[ia];
253 xib = x[ib];
254 yib = y[ib];
255 zib = z[ib];
256 xic = x[ic];
257 yic = y[ic];
258 zic = z[ic];
259
260 xab = xia - xib;
261 yab = yia - yib;
262 zab = zia - zib;
263 rab = sqrt(xab*xab + yab*yab + zab*zab);
264 xcb = xic - xib;
265 ycb = yic - yib;
266 zcb = zic - zib;
267 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
268 xp = ycb*zab - zcb*yab;
269 yp = zcb*xab - xcb*zab;
270 zp = xcb*yab - ycb*xab;
271 rp = sqrt(xp*xp + yp*yp + zp*zp);
272 if (rp != 0.0)
273 {
274 dot = xab*xcb + yab*ycb + zab*zcb;
275 cosine = dot / (rab*rcb);
276 if (cosine < -1.0)
277 cosine = -1.0;
278 if (cosine > 1.0)
279 cosine = 1.0;
280 angle = radian * acos(cosine);
281
282 dt = angle - anat[i];
283 dt2 = dt * dt;
284 dt3 = dt2 * dt;
285 dt4 = dt3 * dt;
286 deddt = units.angunit * acon[i] * dt
287 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
288 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
289 d2eddt2 = units.angunit * acon[i]
290 * (2.0 + 6.0*units.cang*dt + 12.0*units.qang*dt2
291 + 20.0*units.pang*dt3 + 30.0*units.sang*dt4);
292 terma = -radian/ (rab*rab*rp);
293 termc = radian / (rcb*rcb*rp);
294
295 ddtdxia = terma * (yab*zp-zab*yp);
296 ddtdyia = terma * (zab*xp-xab*zp);
297 ddtdzia = terma * (xab*yp-yab*xp);
298 ddtdxic = termc * (ycb*zp-zcb*yp);
299 ddtdyic = termc * (zcb*xp-xcb*zp);
300 ddtdzic = termc * (xcb*yp-ycb*xp);
301 ddtdxib = -ddtdxia - ddtdxic;
302 ddtdyib = -ddtdyia - ddtdyic;
303 ddtdzib = -ddtdzia - ddtdzic;
304 rab2 = 2.0 / (rab*rab);
305 xrab = xab * rab2;
306 yrab = yab * rab2;
307 zrab = zab * rab2;
308 rcb2 = 2.0 / (rcb*rcb);
309 xrcb = xcb * rcb2;
310 yrcb = ycb * rcb2;
311 zrcb = zcb * rcb2;
312 rp2 = 1.0 / (rp*rp);
313 xabp = (yab*zp-zab*yp) * rp2;
314 yabp = (zab*xp-xab*zp) * rp2;
315 zabp = (xab*yp-yab*xp) * rp2;
316 xcbp = (ycb*zp-zcb*yp) * rp2;
317 ycbp = (zcb*xp-xcb*zp) * rp2;
318 zcbp = (xcb*yp-ycb*xp) * rp2;
319
320 dxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab);
321 dxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab);
322 dxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab);
323 dyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab);
324 dyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab);
325 dziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab);
326 dxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb);
327 dxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb);
328 dxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb);
329 dyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb);
330 dyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb);
331 dziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb);
332 dxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp;
333 dxiayic = -terma*xab*yab - ddtdxia*yabp;
334 dxiazic = -terma*xab*zab - ddtdxia*zabp;
335 dyiaxic = -terma*xab*yab - ddtdyia*xabp;
336 dyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp;
337 dyiazic = -terma*yab*zab - ddtdyia*zabp;
338 dziaxic = -terma*xab*zab - ddtdzia*xabp;
339 dziayic = -terma*yab*zab - ddtdzia*yabp;
340 dziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp;
341
342 dxibxia = -dxiaxia - dxiaxic;
343 dxibyia = -dxiayia - dyiaxic;
344 dxibzia = -dxiazia - dziaxic;
345 dyibxia = -dxiayia - dxiayic;
346 dyibyia = -dyiayia - dyiayic;
347 dyibzia = -dyiazia - dziayic;
348 dzibxia = -dxiazia - dxiazic;
349 dzibyia = -dyiazia - dyiazic;
350 dzibzia = -dziazia - dziazic;
351 dxibxic = -dxicxic - dxiaxic;
352 dxibyic = -dxicyic - dxiayic;
353 dxibzic = -dxiczic - dxiazic;
354 dyibxic = -dxicyic - dyiaxic;
355 dyibyic = -dyicyic - dyiayic;
356 dyibzic = -dyiczic - dyiazic;
357 dzibxic = -dxiczic - dziaxic;
358 dzibyic = -dyiczic - dziayic;
359 dzibzic = -dziczic - dziazic;
360 dxibxib = -dxibxia - dxibxic;
361 dxibyib = -dxibyia - dxibyic;
362 dxibzib = -dxibzia - dxibzic;
363 dyibyib = -dyibyia - dyibyic;
364 dyibzib = -dyibzia - dyibzic;
365 dzibzib = -dzibzia - dzibzic;
366
367 if (ia == iatom)
368 {
369 hessx[ia][0] += deddt*dxiaxia + d2eddt2*ddtdxia*ddtdxia;
370 hessx[ia][1] += deddt*dxiayia + d2eddt2*ddtdxia*ddtdyia;
371 hessx[ia][2] += deddt*dxiazia + d2eddt2*ddtdxia*ddtdzia;
372 hessy[ia][0] += deddt*dxiayia + d2eddt2*ddtdyia*ddtdxia;
373 hessy[ia][1] += deddt*dyiayia + d2eddt2*ddtdyia*ddtdyia;
374 hessy[ia][2] += deddt*dyiazia + d2eddt2*ddtdyia*ddtdzia;
375 hessz[ia][0] += deddt*dxiazia + d2eddt2*ddtdzia*ddtdxia;
376 hessz[ia][1] += deddt*dyiazia + d2eddt2*ddtdzia*ddtdyia;
377 hessz[ia][2] += deddt*dziazia + d2eddt2*ddtdzia*ddtdzia;
378 hessx[ib][0] += deddt*dxibxia + d2eddt2*ddtdxia*ddtdxib;
379 hessx[ib][1] += deddt*dyibxia + d2eddt2*ddtdxia*ddtdyib;
380 hessx[ib][2] += deddt*dzibxia + d2eddt2*ddtdxia*ddtdzib;
381 hessy[ib][0] += deddt*dxibyia + d2eddt2*ddtdyia*ddtdxib;
382 hessy[ib][1] += deddt*dyibyia + d2eddt2*ddtdyia*ddtdyib;
383 hessy[ib][2] += deddt*dzibyia + d2eddt2*ddtdyia*ddtdzib;
384 hessz[ib][0] += deddt*dxibzia + d2eddt2*ddtdzia*ddtdxib;
385 hessz[ib][1] += deddt*dyibzia + d2eddt2*ddtdzia*ddtdyib;
386 hessz[ib][2] += deddt*dzibzia + d2eddt2*ddtdzia*ddtdzib;
387 hessx[ic][0] += deddt*dxiaxic + d2eddt2*ddtdxia*ddtdxic;
388 hessx[ic][1] += deddt*dxiayic + d2eddt2*ddtdxia*ddtdyic;
389 hessx[ic][2] += deddt*dxiazic + d2eddt2*ddtdxia*ddtdzic;
390 hessy[ic][0] += deddt*dyiaxic + d2eddt2*ddtdyia*ddtdxic;
391 hessy[ic][1] += deddt*dyiayic + d2eddt2*ddtdyia*ddtdyic;
392 hessy[ic][2] += deddt*dyiazic + d2eddt2*ddtdyia*ddtdzic;
393 hessz[ic][0] += deddt*dziaxic + d2eddt2*ddtdzia*ddtdxic;
394 hessz[ic][1] += deddt*dziayic + d2eddt2*ddtdzia*ddtdyic;
395 hessz[ic][2] += deddt*dziazic + d2eddt2*ddtdzia*ddtdzic;
396 } else if (ib == iatom)
397 {
398 hessx[ib][0] += deddt*dxibxib + d2eddt2*ddtdxib*ddtdxib;
399 hessx[ib][1] += deddt*dxibyib + d2eddt2*ddtdxib*ddtdyib;
400 hessx[ib][2] += deddt*dxibzib + d2eddt2*ddtdxib*ddtdzib;
401 hessy[ib][0] += deddt*dxibyib + d2eddt2*ddtdyib*ddtdxib;
402 hessy[ib][1] += deddt*dyibyib + d2eddt2*ddtdyib*ddtdyib;
403 hessy[ib][2] += deddt*dyibzib + d2eddt2*ddtdyib*ddtdzib;
404 hessz[ib][0] += deddt*dxibzib + d2eddt2*ddtdzib*ddtdxib;
405 hessz[ib][1] += deddt*dyibzib + d2eddt2*ddtdzib*ddtdyib;
406 hessz[ib][2] += deddt*dzibzib + d2eddt2*ddtdzib*ddtdzib;
407 hessx[ia][0] += deddt*dxibxia + d2eddt2*ddtdxib*ddtdxia;
408 hessx[ia][1] += deddt*dxibyia + d2eddt2*ddtdxib*ddtdyia;
409 hessx[ia][2] += deddt*dxibzia + d2eddt2*ddtdxib*ddtdzia;
410 hessy[ia][0] += deddt*dyibxia + d2eddt2*ddtdyib*ddtdxia;
411 hessy[ia][1] += deddt*dyibyia + d2eddt2*ddtdyib*ddtdyia;
412 hessy[ia][2] += deddt*dyibzia + d2eddt2*ddtdyib*ddtdzia;
413 hessz[ia][0] += deddt*dzibxia + d2eddt2*ddtdzib*ddtdxia;
414 hessz[ia][1] += deddt*dzibyia + d2eddt2*ddtdzib*ddtdyia;
415 hessz[ia][2] += deddt*dzibzia + d2eddt2*ddtdzib*ddtdzia;
416 hessx[ic][0] += deddt*dxibxic + d2eddt2*ddtdxib*ddtdxic;
417 hessx[ic][1] += deddt*dxibyic + d2eddt2*ddtdxib*ddtdyic;
418 hessx[ic][2] += deddt*dxibzic + d2eddt2*ddtdxib*ddtdzic;
419 hessy[ic][0] += deddt*dyibxic + d2eddt2*ddtdyib*ddtdxic;
420 hessy[ic][1] += deddt*dyibyic + d2eddt2*ddtdyib*ddtdyic;
421 hessy[ic][2] += deddt*dyibzic + d2eddt2*ddtdyib*ddtdzic;
422 hessz[ic][0] += deddt*dzibxic + d2eddt2*ddtdzib*ddtdxic;
423 hessz[ic][1] += deddt*dzibyic + d2eddt2*ddtdzib*ddtdyic;
424 hessz[ic][2] += deddt*dzibzic + d2eddt2*ddtdzib*ddtdzic;
425 }else if (ic == iatom)
426 {
427 hessx[ic][0] += deddt*dxicxic + d2eddt2*ddtdxic*ddtdxic;
428 hessx[ic][1] += deddt*dxicyic + d2eddt2*ddtdxic*ddtdyic;
429 hessx[ic][2] += deddt*dxiczic + d2eddt2*ddtdxic*ddtdzic;
430 hessy[ic][0] += deddt*dxicyic + d2eddt2*ddtdyic*ddtdxic;
431 hessy[ic][1] += deddt*dyicyic + d2eddt2*ddtdyic*ddtdyic;
432 hessy[ic][2] += deddt*dyiczic + d2eddt2*ddtdyic*ddtdzic;
433 hessz[ic][0] += deddt*dxiczic + d2eddt2*ddtdzic*ddtdxic;
434 hessz[ic][1] += deddt*dyiczic + d2eddt2*ddtdzic*ddtdyic;
435 hessz[ic][2] += deddt*dziczic + d2eddt2*ddtdzic*ddtdzic;
436 hessx[ib][0] += deddt*dxibxic + d2eddt2*ddtdxic*ddtdxib;
437 hessx[ib][1] += deddt*dyibxic + d2eddt2*ddtdxic*ddtdyib;
438 hessx[ib][2] += deddt*dzibxic + d2eddt2*ddtdxic*ddtdzib;
439 hessy[ib][0] += deddt*dxibyic + d2eddt2*ddtdyic*ddtdxib;
440 hessy[ib][1] += deddt*dyibyic + d2eddt2*ddtdyic*ddtdyib;
441 hessy[ib][2] += deddt*dzibyic + d2eddt2*ddtdyic*ddtdzib;
442 hessz[ib][0] += deddt*dxibzic + d2eddt2*ddtdzic*ddtdxib;
443 hessz[ib][1] += deddt*dyibzic + d2eddt2*ddtdzic*ddtdyib;
444 hessz[ib][2] += deddt*dzibzic + d2eddt2*ddtdzic*ddtdzib;
445 hessx[ia][0] += deddt*dxiaxic + d2eddt2*ddtdxic*ddtdxia;
446 hessx[ia][1] += deddt*dyiaxic + d2eddt2*ddtdxic*ddtdyia;
447 hessx[ia][2] += deddt*dziaxic + d2eddt2*ddtdxic*ddtdzia;
448 hessy[ia][0] += deddt*dxiayic + d2eddt2*ddtdyic*ddtdxia;
449 hessy[ia][1] += deddt*dyiayic + d2eddt2*ddtdyic*ddtdyia;
450 hessy[ia][2] += deddt*dziayic + d2eddt2*ddtdyic*ddtdzia;
451 hessz[ia][0] += deddt*dxiazic + d2eddt2*ddtdzic*ddtdxia;
452 hessz[ia][1] += deddt*dyiazic + d2eddt2*ddtdzic*ddtdyia;
453 hessz[ia][2] += deddt*dziazic + d2eddt2*ddtdzic*ddtdzia;
454 }
455 }
456 }
457 }
458 }
459 // ================================================================
460 void eangle2b(int i,int nang,int **i13,int *angin,double *x,double *y,double *z,float *acon,float *anat,double **dea)
461 {
462 int ia,ib,ic,id;
463 double angle,dot,cosine;
464 double dt,dt2,dt3,dt4;
465 double deddt,terma,termc,term;
466 double xia,yia,zia;
467 double xib,yib,zib;
468 double xic,yic,zic;
469 double xid,yid,zid;
470 double xad,yad,zad;
471 double xbd,ybd,zbd;
472 double xcd,ycd,zcd;
473 double xip,yip,zip;
474 double xap,yap,zap,rap2;
475 double xcp,ycp,zcp,rcp2;
476 double xt,yt,zt,rt2,ptrt2;
477 double xm,ym,zm,rm,delta,delta2;
478 double dedxia,dedyia,dedzia;
479 double dedxib,dedyib,dedzib;
480 double dedxic,dedyic,dedzic;
481 double dedxid,dedyid,dedzid;
482 double dedxip,dedyip,dedzip;
483 double dpdxia,dpdyia,dpdzia;
484 double dpdxic,dpdyic,dpdzic;
485
486 ia = i13[i][0];
487 ib = i13[i][1];
488 ic = i13[i][2];
489 id = i13[i][3];
490 xia = x[ia];
491 yia = y[ia];
492 zia = z[ia];
493 xib = x[ib];
494 yib = y[ib];
495 zib = z[ib];
496 xic = x[ic];
497 yic = y[ic];
498 zic = z[ic];
499 xid = x[id];
500 yid = y[id];
501 zid = z[id];
502
503 dea[ia][0] = 0.0;
504 dea[ia][1] = 0.0;
505 dea[ia][2] = 0.0;
506 dea[ib][0] = 0.0;
507 dea[ib][1] = 0.0;
508 dea[ib][2] = 0.0;
509 dea[ic][0] = 0.0;
510 dea[ic][1] = 0.0;
511 dea[ic][2] = 0.0;
512 dea[id][0] = 0.0;
513 dea[id][1] = 0.0;
514 dea[id][2] = 0.0;
515
516 xad = xia - xid;
517 yad = yia - yid;
518 zad = zia - zid;
519 xbd = xib - xid;
520 ybd = yib - yid;
521 zbd = zib - zid;
522 xcd = xic - xid;
523 ycd = yic - yid;
524 zcd = zic - zid;
525 xt = yad*zcd - zad*ycd;
526 yt = zad*xcd - xad*zcd;
527 zt = xad*ycd - yad*xcd;
528 rt2 = xt*xt + yt*yt + zt*zt;
529 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
530 xip = xib + xt*delta;
531 yip = yib + yt*delta;
532 zip = zib + zt*delta;
533 xap = xia - xip;
534 yap = yia - yip;
535 zap = zia - zip;
536 rap2 = xap*xap + yap*yap + zap*zap;
537 xcp = xic - xip;
538 ycp = yic - yip;
539 zcp = zic - zip;
540 rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
541 xm = ycp*zap - zcp*yap;
542 ym = zcp*xap - xcp*zap;
543 zm = xcp*yap - ycp*xap;
544 rm = sqrt(xm*xm + ym*ym + zm*zm);
545 if (rm != 0.0)
546 {
547 dot = xap*xcp + yap*ycp + zap*zcp;
548 cosine = dot / sqrt(rap2*rcp2);
549 if (cosine < -1.0)
550 cosine = -1.0;
551 if (cosine > 1.0)
552 cosine = 1.0;
553 angle = radian * acos(cosine);
554
555 dt = angle - anat[i];
556 dt2 = dt * dt;
557 dt3 = dt2 * dt;
558 dt4 = dt2 * dt2;
559 deddt = units.angunit * acon[i] * dt
560 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
561 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
562
563 deddt = deddt * radian;
564 terma = -deddt / (rap2*rm);
565 termc = deddt / (rcp2*rm);
566 dedxia = terma * (yap*zm-zap*ym);
567 dedyia = terma * (zap*xm-xap*zm);
568 dedzia = terma * (xap*ym-yap*xm);
569 dedxic = termc * (ycp*zm-zcp*ym);
570 dedyic = termc * (zcp*xm-xcp*zm);
571 dedzic = termc * (xcp*ym-ycp*xm);
572 dedxip = -dedxia - dedxic;
573 dedyip = -dedyia - dedyic;
574 dedzip = -dedzia - dedzic;
575
576 delta2 = 2.0 * delta;
577 ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
578 term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
579 dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
580 term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
581 dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
582 term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
583 dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
584 term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
585 dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
586 term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
587 dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
588 term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
589 dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
590 dedxia = dedxia + dpdxia;
591 dedyia = dedyia + dpdyia;
592 dedzia = dedzia + dpdzia;
593 dedxib = dedxip;
594 dedyib = dedyip;
595 dedzib = dedzip;
596 dedxic = dedxic + dpdxic;
597 dedyic = dedyic + dpdyic;
598 dedzic = dedzic + dpdzic;
599 dedxid = -dedxia - dedxib - dedxic;
600 dedyid = -dedyia - dedyib - dedyic;
601 dedzid = -dedzia - dedzib - dedzic;
602
603 dea[ia][0] = dedxia;
604 dea[ia][1] = dedyia;
605 dea[ia][2] = dedzia;
606 dea[ib][0] = dedxib;
607 dea[ib][1] = dedyib;
608 dea[ib][2] = dedzib;
609 dea[ic][0] = dedxic;
610 dea[ic][1] = dedyic;
611 dea[ic][2] = dedzic;
612 dea[id][0] = dedxid;
613 dea[id][1] = dedyid;
614 dea[id][2] = dedzid;
615 }
616 }
617