ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/eimptors.c
Revision: 119
Committed: Tue Jun 16 21:36:00 2009 UTC (11 years, 10 months ago) by wdelano
File size: 32040 byte(s)
Log Message:
printf( to fprintf(pcmlogfile,
Line File contents
1 #define EXTERN extern
2 #include "pcwin.h"
3
4 void eimptors(int nimptors,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in,double *eimptors);
5 void eimptors1(int nimptors,int natom,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in,
6 double *eimptors,double **deimprop);
7 void eimptors2(int iatom,int nimptors,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in,
8 float **hessx,float **hessy,float **hessz);
9
10
11 EXTERN struct t_minim_values {
12 int iprint, ndc, nconst;
13 float dielc;
14 } minim_values;
15
16
17 void eimptors(int nimptors,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in,double *eimptors)
18 {
19 int i,ia,ib,ic,id;
20 double e,cosine,cosine2,cosine3;
21 double v1,v2,v3,s1,s2,s3;
22 double phi1,phi2,phi3;
23 double xia,yia,zia;
24 double xib,yib,zib;
25 double xic,yic,zic;
26 double xid,yid,zid;
27 double xba,yba,zba;
28 double xcb,ycb,zcb;
29 double xdc,ydc,zdc;
30 double xt,yt,zt,rt2;
31 double xu,yu,zu,ru2,rtru;
32
33 *eimptors = 0.0;
34 if (minim_values.iprint)
35 {
36 fprintf(pcmlogfile,"\nAmber Torsion Terms\n");
37 fprintf(pcmlogfile," At1 At2 At3 At4 Types Angle V1 V2 V3 Etor\n");
38 }
39
40 for (i = 0; i < nimptors; i++)
41 {
42 ia = iiprop[i][1];
43 ib = iiprop[i][0];
44 ic = iiprop[i][2];
45 id = iiprop[i][3];
46 if (use[ia] || use[ib] || use[ic] || use[id] )
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 xid = x[id];
58 yid = y[id];
59 zid = z[id];
60 xba = xib - xia;
61 yba = yib - yia;
62 zba = zib - zia;
63 xcb = xic - xib;
64 ycb = yic - yib;
65 zcb = zic - zib;
66 xdc = xid - xic;
67 ydc = yid - yic;
68 zdc = zid - zic;
69 xt = yba*zcb - ycb*zba;
70 yt = zba*xcb - zcb*xba;
71 zt = xba*ycb - xcb*yba;
72 xu = ycb*zdc - ydc*zcb;
73 yu = zcb*xdc - zdc*xcb;
74 zu = xcb*ydc - xdc*ycb;
75 rt2 = xt*xt + yt*yt + zt*zt;
76 ru2 = xu*xu + yu*yu + zu*zu;
77 rtru = sqrt(rt2 * ru2);
78 if (rtru != 0.0)
79 {
80 cosine = (xt*xu + yt*yu + zt*zu) / rtru;
81 v1 = v1in[i];
82 v2 = v2in[i];
83 v3 = v3in[i];
84 s1 = ph1in[i];
85 s2 = ph2in[i];
86 s3 = ph3in[i];
87 cosine2 = cosine * cosine;
88 cosine3 = cosine2 * cosine;
89 phi1 = 1.0 + s1*cosine;
90 phi2 = 1.0 + s2*(2.0*cosine2 - 1.0);
91 phi3 = 1.0 + s3*(4.0*cosine3 - 3.0*cosine);
92 e = units.torsunit * (v1*phi1 + v2*phi2 + v3*phi3);
93 *eimptors += e;
94 if(minim_values.iprint)
95 {
96 fprintf(pcmlogfile,"ImpTor: (%-3d)- (%-3d)- (%-3d)- (%-3d) %d %d %d %d %8.3f %-8.3f %-8.3f %-8.3f = %-8.4f\n",
97 ia,ib,ic,id, type[ia], type[ib],type[ic],type[id], radian*acos(cosine),v1,v2,v3,e);
98 }
99 }
100 }
101 }
102 }
103 // =======================
104 void eimptors1(int nimptors,int natom,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in,
105 double *eimptors,double **deimprop)
106 {
107 int i,ia,ib,ic,id;
108 double e,dedphi,sine;
109 double v1,v2,v3,s1,s2,s3;
110 double cosine,cosine2,cosine3;
111 double phi1,phi2,phi3;
112 double xia,yia,zia;
113 double xib,yib,zib;
114 double xic,yic,zic;
115 double xid,yid,zid;
116 double xba,yba,zba;
117 double xcb,ycb,zcb,rcb;
118 double xdc,ydc,zdc;
119 double xca,yca,zca;
120 double xdb,ydb,zdb;
121 double xt,yt,zt,rt2;
122 double xu,yu,zu,ru2;
123 double xtu,ytu,ztu,rtru;
124 double dphi1,dphi2,dphi3;
125 double dphidxt,dphidyt,dphidzt;
126 double dphidxu,dphidyu,dphidzu;
127 double dphidxia,dphidyia,dphidzia;
128 double dphidxib,dphidyib,dphidzib;
129 double dphidxic,dphidyic,dphidzic;
130 double dphidxid,dphidyid,dphidzid;
131 double dedxia,dedyia,dedzia;
132 double dedxib,dedyib,dedzib;
133 double dedxic,dedyic,dedzic;
134 double dedxid,dedyid,dedzid;
135
136 *eimptors = 0.0;
137
138 for(i=1; i <= natom; i++ )
139 {
140 deimprop[i][0] = 0.0;
141 deimprop[i][1] = 0.0;
142 deimprop[i][2] = 0.0;
143 }
144 for (i = 0; i < nimptors; i++)
145 {
146 ia = iiprop[i][0];
147 ib = iiprop[i][1];
148 ic = iiprop[i][2];
149 id = iiprop[i][3];
150 if (use[ia] || use[ib] || use[ic] || use[id] )
151 {
152 xia = x[ia];
153 yia = y[ia];
154 zia = z[ia];
155 xib = x[ib];
156 yib = y[ib];
157 zib = z[ib];
158 xic = x[ic];
159 yic = y[ic];
160 zic = z[ic];
161 xid = x[id];
162 yid = y[id];
163 zid = z[id];
164 xba = xib - xia;
165 yba = yib - yia;
166 zba = zib - zia;
167 xcb = xic - xib;
168 ycb = yic - yib;
169 zcb = zic - zib;
170 xdc = xid - xic;
171 ydc = yid - yic;
172 zdc = zid - zic;
173 xt = yba*zcb - ycb*zba;
174 yt = zba*xcb - zcb*xba;
175 zt = xba*ycb - xcb*yba;
176 xu = ycb*zdc - ydc*zcb;
177 yu = zcb*xdc - zdc*xcb;
178 zu = xcb*ydc - xdc*ycb;
179 xtu = yt*zu - yu*zt;
180 ytu = zt*xu - zu*xt;
181 ztu = xt*yu - xu*yt;
182 rt2 = xt*xt + yt*yt + zt*zt;
183 ru2 = xu*xu + yu*yu + zu*zu;
184 rtru = sqrt(rt2 * ru2);
185 if (rtru != 0.0)
186 {
187 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
188 cosine = (xt*xu + yt*yu + zt*zu) / rtru;
189 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru);
190 v1 = v1in[i];
191 v2 = v2in[i];
192 v3 = v3in[i];
193 s1 = ph1in[i];
194 s2 = ph2in[i];
195 s3 = ph3in[i];
196 cosine2 = cosine * cosine;
197 cosine3 = cosine2 * cosine;
198 phi1 = 1.0 + s1*cosine;
199 phi2 = 1.0 + s2*(2.0*cosine2 - 1.0);
200 phi3 = 1.0 + s3*(4.0*cosine3 - 3.0*cosine);
201 dphi1 = s1;
202 dphi2 = s2 * (4.0 * cosine);
203 dphi3 = s3 * (12.0*cosine2 - 3.0);
204
205 e = units.torsunit * (v1*phi1 + v2*phi2 + v3*phi3);
206 *eimptors += e;
207 dedphi = -sine * units.torsunit * (v1*dphi1+v2*dphi2+v3*dphi3);
208
209 xca = xic - xia;
210 yca = yic - yia;
211 zca = zic - zia;
212 xdb = xid - xib;
213 ydb = yid - yib;
214 zdb = zid - zib;
215 dphidxt = (yt*zcb - ycb*zt) / (rt2*rcb);
216 dphidyt = (zt*xcb - zcb*xt) / (rt2*rcb);
217 dphidzt = (xt*ycb - xcb*yt) / (rt2*rcb);
218 dphidxu = -(yu*zcb - ycb*zu) / (ru2*rcb);
219 dphidyu = -(zu*xcb - zcb*xu) / (ru2*rcb);
220 dphidzu = -(xu*ycb - xcb*yu) / (ru2*rcb);
221
222 dphidxia = zcb*dphidyt - ycb*dphidzt;
223 dphidyia = xcb*dphidzt - zcb*dphidxt;
224 dphidzia = ycb*dphidxt - xcb*dphidyt;
225 dphidxib = yca*dphidzt - zca*dphidyt + zdc*dphidyu - ydc*dphidzu;
226 dphidyib = zca*dphidxt - xca*dphidzt + xdc*dphidzu - zdc*dphidxu;
227 dphidzib = xca*dphidyt - yca*dphidxt + ydc*dphidxu - xdc*dphidyu;
228 dphidxic = zba*dphidyt - yba*dphidzt + ydb*dphidzu - zdb*dphidyu;
229 dphidyic = xba*dphidzt - zba*dphidxt + zdb*dphidxu - xdb*dphidzu;
230 dphidzic = yba*dphidxt - xba*dphidyt + xdb*dphidyu - ydb*dphidxu;
231 dphidxid = zcb*dphidyu - ycb*dphidzu;
232 dphidyid = xcb*dphidzu - zcb*dphidxu;
233 dphidzid = ycb*dphidxu - xcb*dphidyu;
234 dedxia = dedphi * dphidxia;
235 dedyia = dedphi * dphidyia;
236 dedzia = dedphi * dphidzia;
237 dedxib = dedphi * dphidxib;
238 dedyib = dedphi * dphidyib;
239 dedzib = dedphi * dphidzib;
240 dedxic = dedphi * dphidxic;
241 dedyic = dedphi * dphidyic;
242 dedzic = dedphi * dphidzic;
243 dedxid = dedphi * dphidxid;
244 dedyid = dedphi * dphidyid;
245 dedzid = dedphi * dphidzid;
246
247 deimprop[ia][0] += dedxia;
248 deimprop[ia][1] += dedyia;
249 deimprop[ia][2] += dedzia;
250 deimprop[ib][0] += dedxib;
251 deimprop[ib][1] += dedyib;
252 deimprop[ib][2] += dedzib;
253 deimprop[ic][0] += dedxic;
254 deimprop[ic][1] += dedyic;
255 deimprop[ic][2] += dedzic;
256 deimprop[id][0] += dedxid;
257 deimprop[id][1] += dedyid;
258 deimprop[id][2] += dedzid;
259 }
260 }
261 }
262 }
263 // ====================================
264 void eimptors2(int iatom,int nimptors,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in,
265 float **hessx,float **hessy,float **hessz)
266 {
267 int i,ia,ib,ic,id;
268 double dedphi,d2edphi2,sine;
269 double v1,v2,v3,s1,s2,s3;
270 double cosine,cosine2,cosine3;
271 double xia,yia,zia,xib,yib,zib;
272 double xic,yic,zic,xid,yid,zid;
273 double xba,yba,zba,xcb,ycb,zcb;
274 double xdc,ydc,zdc;
275 double xca,yca,zca,xdb,ydb,zdb;
276 double xt,yt,zt,xu,yu,zu,xtu,ytu,ztu;
277 double rt2,ru2,rtru,rcb;
278 double dphi1,dphi2,dphi3;
279 double d2phi1,d2phi2,d2phi3;
280 double dphidxt,dphidyt,dphidzt;
281 double dphidxu,dphidyu,dphidzu;
282 double dphidxia,dphidyia,dphidzia;
283 double dphidxib,dphidyib,dphidzib;
284 double dphidxic,dphidyic,dphidzic;
285 double dphidxid,dphidyid,dphidzid;
286 double xycb2,xzcb2,yzcb2;
287 double rcbxt,rcbyt,rcbzt,rcbt2;
288 double rcbxu,rcbyu,rcbzu,rcbu2;
289 double dphidxibt,dphidyibt,dphidzibt;
290 double dphidxibu,dphidyibu,dphidzibu;
291 double dphidxict,dphidyict,dphidzict;
292 double dphidxicu,dphidyicu,dphidzicu;
293 double dxiaxia,dyiayia,dziazia,dxibxib,dyibyib,dzibzib;
294 double dxicxic,dyicyic,dziczic,dxidxid,dyidyid,dzidzid;
295 double dxiayia,dxiazia,dyiazia,dxibyib,dxibzib,dyibzib;
296 double dxicyic,dxiczic,dyiczic,dxidyid,dxidzid,dyidzid;
297 double dxiaxib,dxiayib,dxiazib,dyiaxib,dyiayib,dyiazib;
298 double dziaxib,dziayib,dziazib,dxiaxic,dxiayic,dxiazic;
299 double dyiaxic,dyiayic,dyiazic,dziaxic,dziayic,dziazic;
300 double dxiaxid,dxiayid,dxiazid,dyiaxid,dyiayid,dyiazid;
301 double dziaxid,dziayid,dziazid,dxibxic,dxibyic,dxibzic;
302 double dyibxic,dyibyic,dyibzic,dzibxic,dzibyic,dzibzic;
303 double dxibxid,dxibyid,dxibzid,dyibxid,dyibyid,dyibzid;
304 double dzibxid,dzibyid,dzibzid,dxicxid,dxicyid,dxiczid;
305 double dyicxid,dyicyid,dyiczid,dzicxid,dzicyid,dziczid;
306
307 for (i = 0; i < nimptors; i++)
308 {
309 ia = iiprop[i][0];
310 ib = iiprop[i][1];
311 ic = iiprop[i][2];
312 id = iiprop[i][3];
313 if (ia == iatom || ib == iatom || ic == iatom || id == iatom)
314 {
315 xia = x[ia];
316 yia = y[ia];
317 zia = z[ia];
318 xib = x[ib];
319 yib = y[ib];
320 zib = z[ib];
321 xic = x[ic];
322 yic = y[ic];
323 zic = z[ic];
324 xid = x[id];
325 yid = y[id];
326 zid = z[id];
327 xba = xib - xia;
328 yba = yib - yia;
329 zba = zib - zia;
330 xcb = xic - xib;
331 ycb = yic - yib;
332 zcb = zic - zib;
333 xdc = xid - xic;
334 ydc = yid - yic;
335 zdc = zid - zic;
336 xt = yba*zcb - ycb*zba;
337 yt = zba*xcb - zcb*xba;
338 zt = xba*ycb - xcb*yba;
339 xu = ycb*zdc - ydc*zcb;
340 yu = zcb*xdc - zdc*xcb;
341 zu = xcb*ydc - xdc*ycb;
342 xtu = yt*zu - yu*zt;
343 ytu = zt*xu - zu*xt;
344 ztu = xt*yu - xu*yt;
345 rt2 = xt*xt + yt*yt + zt*zt;
346 ru2 = xu*xu + yu*yu + zu*zu;
347 rtru = sqrt(rt2 * ru2);
348 if (rtru != 0.0)
349 {
350 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
351 cosine = (xt*xu + yt*yu + zt*zu) / rtru;
352 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru);
353 v1 = v1in[i];
354 v2 = v2in[i];
355 v3 = v3in[i];
356 s1 = ph1in[i];
357 s2 = ph2in[i];
358 s3 = ph3in[i];
359 cosine2 = cosine * cosine;
360 cosine3 = cosine2 * cosine;
361
362 dphi1 = s1;
363 dphi2 = s2 * (4.0 * cosine);
364 dphi3 = s3 * (12.0*cosine2 - 3.0);
365 d2phi1 = -s1 * cosine;
366 d2phi2 = -s2 * (8.0*cosine2 - 4.0);
367 d2phi3 = -s3 * (36.0*cosine3 - 27.0*cosine);
368
369 dedphi = -sine * units.torsunit * (v1*dphi1 + v2*dphi2 + v3*dphi3);
370 d2edphi2 = units.torsunit * (v1*d2phi1 + v2*d2phi2 + v3*d2phi3);
371
372 xca = xic - xia;
373 yca = yic - yia;
374 zca = zic - zia;
375 xdb = xid - xib;
376 ydb = yid - yib;
377 zdb = zid - zib;
378 dphidxt = (yt*zcb - ycb*zt) / (rt2*rcb);
379 dphidyt = (zt*xcb - zcb*xt) / (rt2*rcb);
380 dphidzt = (xt*ycb - xcb*yt) / (rt2*rcb);
381 dphidxu = -(yu*zcb - ycb*zu) / (ru2*rcb);
382 dphidyu = -(zu*xcb - zcb*xu) / (ru2*rcb);
383 dphidzu = -(xu*ycb - xcb*yu) / (ru2*rcb);
384 xycb2 = xcb*xcb + ycb*ycb;
385 xzcb2 = xcb*xcb + zcb*zcb;
386 yzcb2 = ycb*ycb + zcb*zcb;
387 rcbxt = -2.0 * rcb * dphidxt;
388 rcbyt = -2.0 * rcb * dphidyt;
389 rcbzt = -2.0 * rcb * dphidzt;
390 rcbt2 = rcb * rt2;
391 rcbxu = 2.0 * rcb * dphidxu;
392 rcbyu = 2.0 * rcb * dphidyu;
393 rcbzu = 2.0 * rcb * dphidzu;
394 rcbu2 = rcb * ru2;
395 dphidxibt = yca*dphidzt - zca*dphidyt;
396 dphidxibu = zdc*dphidyu - ydc*dphidzu;
397 dphidyibt = zca*dphidxt - xca*dphidzt;
398 dphidyibu = xdc*dphidzu - zdc*dphidxu;
399 dphidzibt = xca*dphidyt - yca*dphidxt;
400 dphidzibu = ydc*dphidxu - xdc*dphidyu;
401 dphidxict = zba*dphidyt - yba*dphidzt;
402 dphidxicu = ydb*dphidzu - zdb*dphidyu;
403 dphidyict = xba*dphidzt - zba*dphidxt;
404 dphidyicu = zdb*dphidxu - xdb*dphidzu;
405 dphidzict = yba*dphidxt - xba*dphidyt;
406 dphidzicu = xdb*dphidyu - ydb*dphidxu;
407
408 dphidxia = zcb*dphidyt - ycb*dphidzt;
409 dphidyia = xcb*dphidzt - zcb*dphidxt;
410 dphidzia = ycb*dphidxt - xcb*dphidyt;
411 dphidxib = dphidxibt + dphidxibu;
412 dphidyib = dphidyibt + dphidyibu;
413 dphidzib = dphidzibt + dphidzibu;
414 dphidxic = dphidxict + dphidxicu;
415 dphidyic = dphidyict + dphidyicu;
416 dphidzic = dphidzict + dphidzicu;
417 dphidxid = zcb*dphidyu - ycb*dphidzu;
418 dphidyid = xcb*dphidzu - zcb*dphidxu;
419 dphidzid = ycb*dphidxu - xcb*dphidyu;
420
421 dxiaxia = rcbxt*dphidxia;
422 dxiayia = rcbxt*dphidyia - zcb*rcb/rt2;
423 dxiazia = rcbxt*dphidzia + ycb*rcb/rt2;
424 dxiaxib = rcbxt*dphidxibt + xcb*(zca*ycb-yca*zcb)/rcbt2;
425 dxiayib = rcbxt*dphidyibt + dphidzt + (xca*zcb*xcb+zca*yzcb2)/rcbt2;
426 dxiazib = rcbxt*dphidzibt - dphidyt - (xca*ycb*xcb+yca*yzcb2)/rcbt2;
427 dxiaxic = rcbxt*dphidxict + xcb*xt/rcbt2;
428 dxiayic = rcbxt*dphidyict - dphidzt - (xba*zcb*xcb+zba*yzcb2)/rcbt2;
429 dxiazic = rcbxt*dphidzict + dphidyt + (xba*ycb*xcb+yba*yzcb2)/rcbt2;
430 dxiaxid = 0.0;
431 dxiayid = 0.0;
432 dxiazid = 0.0;
433 dyiayia = rcbyt*dphidyia;
434 dyiazia = rcbyt*dphidzia - xcb*rcb/rt2;
435 dyiaxib = rcbyt*dphidxibt - dphidzt - (yca*zcb*ycb+zca*xzcb2)/rcbt2;
436 dyiayib = rcbyt*dphidyibt + ycb*(xca*zcb-zca*xcb)/rcbt2;
437 dyiazib = rcbyt*dphidzibt + dphidxt + (yca*xcb*ycb+xca*xzcb2)/rcbt2;
438 dyiaxic = rcbyt*dphidxict + dphidzt + (yba*zcb*ycb+zba*xzcb2)/rcbt2;
439 dyiayic = rcbyt*dphidyict + ycb*yt/rcbt2;
440 dyiazic = rcbyt*dphidzict - dphidxt - (yba*xcb*ycb+xba*xzcb2)/rcbt2;
441 dyiaxid = 0.0;
442 dyiayid = 0.0;
443 dyiazid = 0.0;
444 dziazia = rcbzt*dphidzia;
445 dziaxib = rcbzt*dphidxibt + dphidyt + (zca*ycb*zcb+yca*xycb2)/rcbt2;
446 dziayib = rcbzt*dphidyibt - dphidxt - (zca*xcb*zcb+xca*xycb2)/rcbt2;
447 dziazib = rcbzt*dphidzibt + zcb*(yca*xcb-xca*ycb)/rcbt2;
448 dziaxic = rcbzt*dphidxict - dphidyt - (zba*ycb*zcb+yba*xycb2)/rcbt2;
449 dziayic = rcbzt*dphidyict + dphidxt + (zba*xcb*zcb+xba*xycb2)/rcbt2;
450 dziazic = rcbzt*dphidzict + zcb*zt/rcbt2;
451 dziaxid = 0.0;
452 dziayid = 0.0;
453 dziazid = 0.0;
454 dxibxic = -xcb*dphidxib/(rcb*rcb) - (yca*(zba*xcb+yt)-zca*(yba*xcb-zt))/rcbt2
455 - 2.0*(yt*zba-yba*zt)*dphidxibt/rt2 - (zdc*(ydb*xcb+zu)-ydc*(zdb*xcb-yu))/rcbu2
456 + 2.0*(yu*zdb-ydb*zu)*dphidxibu/ru2;
457 dxibyic = -ycb*dphidxib/(rcb*rcb) + dphidzt + dphidzu - (yca*(zba*ycb-xt)+zca*(xba*xcb+zcb*zba))/rcbt2
458 - 2.0*(zt*xba-zba*xt)*dphidxibt/rt2 + (zdc*(xdb*xcb+zcb*zdb)+ydc*(zdb*ycb+xu))/rcbu2
459 + 2.0*(zu*xdb-zdb*xu)*dphidxibu/ru2;
460 dxibxid = rcbxu*dphidxibu + xcb*xu/rcbu2;
461 dxibyid = rcbyu*dphidxibu - dphidzu - (ydc*zcb*ycb+zdc*xzcb2)/rcbu2;
462 dxibzid = rcbzu*dphidxibu + dphidyu + (zdc*ycb*zcb+ydc*xycb2)/rcbu2;
463 dyibzib = ycb*dphidzib/(rcb*rcb) - (xca*(xca*xcb+zcb*zca)+yca*(ycb*xca+zt))/rcbt2
464 - 2.0*(xt*zca-xca*zt)*dphidzibt/rt2 + (ydc*(xdc*ycb-zu)+xdc*(xdc*xcb+zcb*zdc))/rcbu2
465 + 2.0*(xu*zdc-xdc*zu)*dphidzibu/ru2;
466 dyibxic = -xcb*dphidyib/(rcb*rcb) - dphidzt - dphidzu + (xca*(zba*xcb+yt)+zca*(zba*zcb+ycb*yba))/rcbt2
467 - 2.0*(yt*zba-yba*zt)*dphidyibt/rt2 - (zdc*(zdb*zcb+ycb*ydb)+xdc*(zdb*xcb-yu))/rcbu2
468 + 2.0*(yu*zdb-ydb*zu)*dphidyibu/ru2;
469 dyibyic = -ycb*dphidyib/(rcb*rcb) - (zca*(xba*ycb+zt)-xca*(zba*ycb-xt))/rcbt2
470 - 2.0*(zt*xba-zba*xt)*dphidyibt/rt2 - (xdc*(zdb*ycb+xu)-zdc*(xdb*ycb-zu))/rcbu2
471 + 2.0*(zu*xdb-zdb*xu)*dphidyibu/ru2;
472 dyibxid = rcbxu*dphidyibu + dphidzu + (xdc*zcb*xcb+zdc*yzcb2)/rcbu2;
473 dyibyid = rcbyu*dphidyibu + ycb*yu/rcbu2;
474 dyibzid = rcbzu*dphidyibu - dphidxu - (zdc*xcb*zcb+xdc*xycb2)/rcbu2;
475 dzibxic = -xcb*dphidzib/(rcb*rcb) + dphidyt + dphidyu - (xca*(yba*xcb-zt)+yca*(zba*zcb+ycb*yba))/rcbt2
476 - 2.0*(yt*zba-yba*zt)*dphidzibt/rt2 + (ydc*(zdb*zcb+ycb*ydb)+xdc*(ydb*xcb+zu))/rcbu2
477 + 2.0*(yu*zdb-ydb*zu)*dphidzibu/ru2;
478 dzibzic = -zcb*dphidzib/(rcb*rcb) - (xca*(yba*zcb+xt)-yca*(xba*zcb-yt))/rcbt2
479 - 2.0*(xt*yba-xba*yt)*dphidzibt/rt2 - (ydc*(xdb*zcb+yu)-xdc*(ydb*zcb-xu))/rcbu2
480 + 2.0*(xu*ydb-xdb*yu)*dphidzibu/ru2;
481 dzibxid = rcbxu*dphidzibu - dphidyu - (xdc*ycb*xcb+ydc*yzcb2)/rcbu2;
482 dzibyid = rcbyu*dphidzibu + dphidxu + (ydc*xcb*ycb+xdc*xzcb2)/rcbu2;
483 dzibzid = rcbzu*dphidzibu + zcb*zu/rcbu2;
484 dxicxid = rcbxu*dphidxicu - xcb*(zdb*ycb-ydb*zcb)/rcbu2;
485 dxicyid = rcbyu*dphidxicu + dphidzu + (ydb*zcb*ycb+zdb*xzcb2)/rcbu2;
486 dxiczid = rcbzu*dphidxicu - dphidyu - (zdb*ycb*zcb+ydb*xycb2)/rcbu2;
487 dyicxid = rcbxu*dphidyicu - dphidzu - (xdb*zcb*xcb+zdb*yzcb2)/rcbu2;
488 dyicyid = rcbyu*dphidyicu - ycb*(xdb*zcb-zdb*xcb)/rcbu2;
489 dyiczid = rcbzu*dphidyicu + dphidxu + (zdb*xcb*zcb+xdb*xycb2)/rcbu2;
490 dzicxid = rcbxu*dphidzicu + dphidyu + (xdb*ycb*xcb+ydb*yzcb2)/rcbu2;
491 dzicyid = rcbyu*dphidzicu - dphidxu - (ydb*xcb*ycb+xdb*xzcb2)/rcbu2;
492 dziczid = rcbzu*dphidzicu - zcb*(ydb*xcb-xdb*ycb)/rcbu2;
493 dxidxid = rcbxu*dphidxid;
494 dxidyid = rcbxu*dphidyid + zcb*rcb/ru2;
495 dxidzid = rcbxu*dphidzid - ycb*rcb/ru2;
496 dyidyid = rcbyu*dphidyid;
497 dyidzid = rcbyu*dphidzid + xcb*rcb/ru2;
498 dzidzid = rcbzu*dphidzid;
499
500 dxibxib = -dxiaxib - dxibxic - dxibxid;
501 dxibyib = -dyiaxib - dxibyic - dxibyid;
502 dxibzib = -dxiazib - dzibxic - dzibxid;
503 dxibzic = -dziaxib - dxibzib - dxibzid;
504 dyibyib = -dyiayib - dyibyic - dyibyid;
505 dyibzic = -dziayib - dyibzib - dyibzid;
506 dzibzib = -dziazib - dzibzic - dzibzid;
507 dzibyic = -dyiazib - dyibzib - dzibyid;
508 dxicxic = -dxiaxic - dxibxic - dxicxid;
509 dxicyic = -dyiaxic - dyibxic - dxicyid;
510 dxiczic = -dziaxic - dzibxic - dxiczid;
511 dyicyic = -dyiayic - dyibyic - dyicyid;
512 dyiczic = -dziayic - dzibyic - dyiczid;
513 dziczic = -dziazic - dzibzic - dziczid;
514
515 if (iatom == ia)
516 {
517 hessx[ia][0] += dedphi*dxiaxia + d2edphi2*dphidxia*dphidxia;
518 hessy[ia][0] += dedphi*dxiayia + d2edphi2*dphidxia*dphidyia;
519 hessz[ia][0] += dedphi*dxiazia + d2edphi2*dphidxia*dphidzia;
520 hessx[ia][1] += dedphi*dxiayia + d2edphi2*dphidxia*dphidyia;
521 hessy[ia][1] += dedphi*dyiayia + d2edphi2*dphidyia*dphidyia;
522 hessz[ia][1] += dedphi*dyiazia + d2edphi2*dphidyia*dphidzia;
523 hessx[ia][2] += dedphi*dxiazia + d2edphi2*dphidxia*dphidzia;
524 hessy[ia][2] += dedphi*dyiazia + d2edphi2*dphidyia*dphidzia;
525 hessz[ia][2] += dedphi*dziazia + d2edphi2*dphidzia*dphidzia;
526 hessx[ib][0] += dedphi*dxiaxib + d2edphi2*dphidxia*dphidxib;
527 hessy[ib][0] += dedphi*dyiaxib + d2edphi2*dphidyia*dphidxib;
528 hessz[ib][0] += dedphi*dziaxib + d2edphi2*dphidzia*dphidxib;
529 hessx[ib][1] += dedphi*dxiayib + d2edphi2*dphidxia*dphidyib;
530 hessy[ib][1] += dedphi*dyiayib + d2edphi2*dphidyia*dphidyib;
531 hessz[ib][1] += dedphi*dziayib + d2edphi2*dphidzia*dphidyib;
532 hessx[ib][2] += dedphi*dxiazib + d2edphi2*dphidxia*dphidzib;
533 hessy[ib][2] += dedphi*dyiazib + d2edphi2*dphidyia*dphidzib;
534 hessz[ib][2] += dedphi*dziazib + d2edphi2*dphidzia*dphidzib;
535 hessx[ic][0] += dedphi*dxiaxic + d2edphi2*dphidxia*dphidxic;
536 hessy[ic][0] += dedphi*dyiaxic + d2edphi2*dphidyia*dphidxic;
537 hessz[ic][0] += dedphi*dziaxic + d2edphi2*dphidzia*dphidxic;
538 hessx[ic][1] += dedphi*dxiayic + d2edphi2*dphidxia*dphidyic;
539 hessy[ic][1] += dedphi*dyiayic + d2edphi2*dphidyia*dphidyic;
540 hessz[ic][1] += dedphi*dziayic + d2edphi2*dphidzia*dphidyic;
541 hessx[ic][2] += dedphi*dxiazic + d2edphi2*dphidxia*dphidzic;
542 hessy[ic][2] += dedphi*dyiazic + d2edphi2*dphidyia*dphidzic;
543 hessz[ic][2] += dedphi*dziazic + d2edphi2*dphidzia*dphidzic;
544 hessx[id][0] += dedphi*dxiaxid + d2edphi2*dphidxia*dphidxid;
545 hessy[id][0] += dedphi*dyiaxid + d2edphi2*dphidyia*dphidxid;
546 hessz[id][0] += dedphi*dziaxid + d2edphi2*dphidzia*dphidxid;
547 hessx[id][1] += dedphi*dxiayid + d2edphi2*dphidxia*dphidyid;
548 hessy[id][1] += dedphi*dyiayid + d2edphi2*dphidyia*dphidyid;
549 hessz[id][1] += dedphi*dziayid + d2edphi2*dphidzia*dphidyid;
550 hessx[id][2] += dedphi*dxiazid + d2edphi2*dphidxia*dphidzid;
551 hessy[id][2] += dedphi*dyiazid + d2edphi2*dphidyia*dphidzid;
552 hessz[id][2] += dedphi*dziazid + d2edphi2*dphidzia*dphidzid;
553 } else if (iatom == ib)
554 {
555 hessx[ib][0] += dedphi*dxibxib + d2edphi2*dphidxib*dphidxib;
556 hessy[ib][0] += dedphi*dxibyib + d2edphi2*dphidxib*dphidyib;
557 hessz[ib][0] += dedphi*dxibzib + d2edphi2*dphidxib*dphidzib;
558 hessx[ib][1] += dedphi*dxibyib + d2edphi2*dphidxib*dphidyib;
559 hessy[ib][1] += dedphi*dyibyib + d2edphi2*dphidyib*dphidyib;
560 hessz[ib][1] += dedphi*dyibzib + d2edphi2*dphidyib*dphidzib;
561 hessx[ib][2] += dedphi*dxibzib + d2edphi2*dphidxib*dphidzib;
562 hessy[ib][2] += dedphi*dyibzib + d2edphi2*dphidyib*dphidzib;
563 hessz[ib][2] += dedphi*dzibzib + d2edphi2*dphidzib*dphidzib;
564 hessx[ia][0] += dedphi*dxiaxib + d2edphi2*dphidxib*dphidxia;
565 hessy[ia][0] += dedphi*dxiayib + d2edphi2*dphidyib*dphidxia;
566 hessz[ia][0] += dedphi*dxiazib + d2edphi2*dphidzib*dphidxia;
567 hessx[ia][1] += dedphi*dyiaxib + d2edphi2*dphidxib*dphidyia;
568 hessy[ia][1] += dedphi*dyiayib + d2edphi2*dphidyib*dphidyia;
569 hessz[ia][1] += dedphi*dyiazib + d2edphi2*dphidzib*dphidyia;
570 hessx[ia][2] += dedphi*dziaxib + d2edphi2*dphidxib*dphidzia;
571 hessy[ia][2] += dedphi*dziayib + d2edphi2*dphidyib*dphidzia;
572 hessz[ia][2] += dedphi*dziazib + d2edphi2*dphidzib*dphidzia;
573 hessx[ic][0] += dedphi*dxibxic + d2edphi2*dphidxib*dphidxic;
574 hessy[ic][0] += dedphi*dyibxic + d2edphi2*dphidyib*dphidxic;
575 hessz[ic][0] += dedphi*dzibxic + d2edphi2*dphidzib*dphidxic;
576 hessx[ic][1] += dedphi*dxibyic + d2edphi2*dphidxib*dphidyic;
577 hessy[ic][1] += dedphi*dyibyic + d2edphi2*dphidyib*dphidyic;
578 hessz[ic][1] += dedphi*dzibyic + d2edphi2*dphidzib*dphidyic;
579 hessx[ic][2] += dedphi*dxibzic + d2edphi2*dphidxib*dphidzic;
580 hessy[ic][2] += dedphi*dyibzic + d2edphi2*dphidyib*dphidzic;
581 hessz[ic][2] += dedphi*dzibzic + d2edphi2*dphidzib*dphidzic;
582 hessx[id][0] += dedphi*dxibxid + d2edphi2*dphidxib*dphidxid;
583 hessy[id][0] += dedphi*dyibxid + d2edphi2*dphidyib*dphidxid;
584 hessz[id][0] += dedphi*dzibxid + d2edphi2*dphidzib*dphidxid;
585 hessx[id][1] += dedphi*dxibyid + d2edphi2*dphidxib*dphidyid;
586 hessy[id][1] += dedphi*dyibyid + d2edphi2*dphidyib*dphidyid;
587 hessz[id][1] += dedphi*dzibyid + d2edphi2*dphidzib*dphidyid;
588 hessx[id][2] += dedphi*dxibzid + d2edphi2*dphidxib*dphidzid;
589 hessy[id][2] += dedphi*dyibzid + d2edphi2*dphidyib*dphidzid;
590 hessz[id][2] += dedphi*dzibzid + d2edphi2*dphidzib*dphidzid;
591 } else if (iatom == ic)
592 {
593 hessx[ic][0] += dedphi*dxicxic + d2edphi2*dphidxic*dphidxic;
594 hessy[ic][0] += dedphi*dxicyic + d2edphi2*dphidxic*dphidyic;
595 hessz[ic][0] += dedphi*dxiczic + d2edphi2*dphidxic*dphidzic;
596 hessx[ic][1] += dedphi*dxicyic + d2edphi2*dphidxic*dphidyic;
597 hessy[ic][1] += dedphi*dyicyic + d2edphi2*dphidyic*dphidyic;
598 hessz[ic][1] += dedphi*dyiczic + d2edphi2*dphidyic*dphidzic;
599 hessx[ic][2] += dedphi*dxiczic + d2edphi2*dphidxic*dphidzic;
600 hessy[ic][2] += dedphi*dyiczic + d2edphi2*dphidyic*dphidzic;
601 hessz[ic][2] += dedphi*dziczic + d2edphi2*dphidzic*dphidzic;
602 hessx[ia][0] += dedphi*dxiaxic + d2edphi2*dphidxic*dphidxia;
603 hessy[ia][0] += dedphi*dxiayic + d2edphi2*dphidyic*dphidxia;
604 hessz[ia][0] += dedphi*dxiazic + d2edphi2*dphidzic*dphidxia;
605 hessx[ia][1] += dedphi*dyiaxic + d2edphi2*dphidxic*dphidyia;
606 hessy[ia][1] += dedphi*dyiayic + d2edphi2*dphidyic*dphidyia;
607 hessz[ia][1] += dedphi*dyiazic + d2edphi2*dphidzic*dphidyia;
608 hessx[ia][2] += dedphi*dziaxic + d2edphi2*dphidxic*dphidzia;
609 hessy[ia][2] += dedphi*dziayic + d2edphi2*dphidyic*dphidzia;
610 hessz[ia][2] += dedphi*dziazic + d2edphi2*dphidzic*dphidzia;
611 hessx[ib][0] += dedphi*dxibxic + d2edphi2*dphidxic*dphidxib;
612 hessy[ib][0] += dedphi*dxibyic + d2edphi2*dphidyic*dphidxib;
613 hessz[ib][0] += dedphi*dxibzic + d2edphi2*dphidzic*dphidxib;
614 hessx[ib][1] += dedphi*dyibxic + d2edphi2*dphidxic*dphidyib;
615 hessy[ib][1] += dedphi*dyibyic + d2edphi2*dphidyic*dphidyib;
616 hessz[ib][1] += dedphi*dyibzic + d2edphi2*dphidzic*dphidyib;
617 hessx[ib][2] += dedphi*dzibxic + d2edphi2*dphidxic*dphidzib;
618 hessy[ib][2] += dedphi*dzibyic + d2edphi2*dphidyic*dphidzib;
619 hessz[ib][2] += dedphi*dzibzic + d2edphi2*dphidzic*dphidzib;
620 hessx[id][0] += dedphi*dxicxid + d2edphi2*dphidxic*dphidxid;
621 hessy[id][0] += dedphi*dyicxid + d2edphi2*dphidyic*dphidxid;
622 hessz[id][0] += dedphi*dzicxid + d2edphi2*dphidzic*dphidxid;
623 hessx[id][1] += dedphi*dxicyid + d2edphi2*dphidxic*dphidyid;
624 hessy[id][1] += dedphi*dyicyid + d2edphi2*dphidyic*dphidyid;
625 hessz[id][1] += dedphi*dzicyid + d2edphi2*dphidzic*dphidyid;
626 hessx[id][2] += dedphi*dxiczid + d2edphi2*dphidxic*dphidzid;
627 hessy[id][2] += dedphi*dyiczid + d2edphi2*dphidyic*dphidzid;
628 hessz[id][2] += dedphi*dziczid + d2edphi2*dphidzic*dphidzid;
629 } else if (iatom == id)
630 {
631 hessx[id][0] += dedphi*dxidxid + d2edphi2*dphidxid*dphidxid;
632 hessy[id][0] += dedphi*dxidyid + d2edphi2*dphidxid*dphidyid;
633 hessz[id][0] += dedphi*dxidzid + d2edphi2*dphidxid*dphidzid;
634 hessx[id][1] += dedphi*dxidyid + d2edphi2*dphidxid*dphidyid;
635 hessy[id][1] += dedphi*dyidyid + d2edphi2*dphidyid*dphidyid;
636 hessz[id][1] += dedphi*dyidzid + d2edphi2*dphidyid*dphidzid;
637 hessx[id][2] += dedphi*dxidzid + d2edphi2*dphidxid*dphidzid;
638 hessy[id][2] += dedphi*dyidzid + d2edphi2*dphidyid*dphidzid;
639 hessz[id][2] += dedphi*dzidzid + d2edphi2*dphidzid*dphidzid;
640 hessx[ia][0] += dedphi*dxiaxid + d2edphi2*dphidxid*dphidxia;
641 hessy[ia][0] += dedphi*dxiayid + d2edphi2*dphidyid*dphidxia;
642 hessz[ia][0] += dedphi*dxiazid + d2edphi2*dphidzid*dphidxia;
643 hessx[ia][1] += dedphi*dyiaxid + d2edphi2*dphidxid*dphidyia;
644 hessy[ia][1] += dedphi*dyiayid + d2edphi2*dphidyid*dphidyia;
645 hessz[ia][1] += dedphi*dyiazid + d2edphi2*dphidzid*dphidyia;
646 hessx[ia][2] += dedphi*dziaxid + d2edphi2*dphidxid*dphidzia;
647 hessy[ia][2] += dedphi*dziayid + d2edphi2*dphidyid*dphidzia;
648 hessz[ia][2] += dedphi*dziazid + d2edphi2*dphidzid*dphidzia;
649 hessx[ib][0] += dedphi*dxibxid + d2edphi2*dphidxid*dphidxib;
650 hessy[ib][0] += dedphi*dxibyid + d2edphi2*dphidyid*dphidxib;
651 hessz[ib][0] += dedphi*dxibzid + d2edphi2*dphidzid*dphidxib;
652 hessx[ib][1] += dedphi*dyibxid + d2edphi2*dphidxid*dphidyib;
653 hessy[ib][1] += dedphi*dyibyid + d2edphi2*dphidyid*dphidyib;
654 hessz[ib][1] += dedphi*dyibzid + d2edphi2*dphidzid*dphidyib;
655 hessx[ib][2] += dedphi*dzibxid + d2edphi2*dphidxid*dphidzib;
656 hessy[ib][2] += dedphi*dzibyid + d2edphi2*dphidyid*dphidzib;
657 hessz[ib][2] += dedphi*dzibzid + d2edphi2*dphidzid*dphidzib;
658 hessx[ic][0] += dedphi*dxicxid + d2edphi2*dphidxid*dphidxic;
659 hessy[ic][0] += dedphi*dxicyid + d2edphi2*dphidyid*dphidxic;
660 hessz[ic][0] += dedphi*dxiczid + d2edphi2*dphidzid*dphidxic;
661 hessx[ic][1] += dedphi*dyicxid + d2edphi2*dphidxid*dphidyic;
662 hessy[ic][1] += dedphi*dyicyid + d2edphi2*dphidyid*dphidyic;
663 hessz[ic][1] += dedphi*dyiczid + d2edphi2*dphidzid*dphidyic;
664 hessx[ic][2] += dedphi*dzicxid + d2edphi2*dphidxid*dphidzic;
665 hessy[ic][2] += dedphi*dzicyid + d2edphi2*dphidyid*dphidzic;
666 hessz[ic][2] += dedphi*dziczid + d2edphi2*dphidzid*dphidzic;
667 }
668 }
669 }
670 }
671 }