ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/eobpw.c
Revision: 93
Committed: Sat Jan 17 23:24:55 2009 UTC (13 years, 5 months ago) by wdelano
File size: 137922 byte(s)
Log Message: 
branch created
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5     #include "energies.h"
6     #include "angles.h"
7     #include "derivs.h"
8     #include "hess.h"
9    
10    
11     EXTERN struct t_minim_values {
12     int iprint, ndc, nconst;
13     float dielc;
14     } minim_values;
15    
16     void eopbend_wilson(void);
17     void eopbend_wilson1(void);
18     void eopbend_wilson2(int);
19     double term1aa (double,double,double,double,double, double, double);
20    
21     void eopbend_wilson()
22     {
23     int i,ia,ib,ic,id;
24     double e,angle,dot,cosine;
25     double dt,dt2;
26     double xia,yia,zia;
27     double xib,yib,zib;
28     double xic,yic,zic;
29     double xid,yid,zid;
30     double xab,yab,zab,rab2,rab;
31     double xbc,ybc,zbc,rbc2,rbc;
32     double xbd,ybd,zbd,rbd2,rbd;
33     double sine;
34    
35     energies.eopb = 0.0;
36    
37     if (minim_values.iprint)
38     {
39 wdelano 58 fprintf(pcmlogfile,"\nWilson Out of Plane Bending Terms\n");
40     fprintf(pcmlogfile," At1 At2 At3 At4 Angle Oopconst Eoop\n");
41 tjod 3 }
42    
43     for (i=0; i < angles.nang; i++)
44     {
45     if (angles.angin[i] == TRUE)
46     {
47     ia = angles.i13[i][0];
48     ib = angles.i13[i][1];
49     ic = angles.i13[i][2];
50     id = angles.i13[i][3];
51     if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use)
52     {
53     xia = atom[ia].x;
54     yia = atom[ia].y;
55     zia = atom[ia].z;
56     xib = atom[ib].x;
57     yib = atom[ib].y;
58     zib = atom[ib].z;
59     xic = atom[ic].x;
60     yic = atom[ic].y;
61     zic = atom[ic].z;
62     xid = atom[id].x;
63     yid = atom[id].y;
64     zid = atom[id].z;
65    
66     xab = xia - xib;
67     yab = yia - yib;
68     zab = zia - zib;
69     rab2 = (xab*xab+yab*yab+zab*zab);
70     rab = sqrt(rab2);
71    
72     xbc = xic - xib;
73     ybc = yic - yib;
74     zbc = zic - zib;
75     rbc2 = xbc*xbc+ybc*ybc+zbc*zbc;
76     rbc = sqrt(rbc2);
77    
78     if (rab2 != 0.0 && rbc2 != 0.0)
79     {
80     dot = xab*xbc + yab*ybc + zab*zbc;
81     cosine = dot /(rab*rbc);
82     if (cosine < -1.0)
83     cosine = -1.0;
84     if (cosine > 1.0)
85     cosine = 1.0;
86     sine = sqrt(1.00-cosine*cosine);
87    
88     xbd = atom[id].x - atom[ib].x;
89     ybd = atom[id].y - atom[ib].y;
90     zbd = atom[id].z - atom[ib].z;
91     rbd2 = xbd*xbd+ybd*ybd+zbd*zbd;
92     rbd = sqrt(rbd2);
93    
94     cosine = ( xbd*(yab*zbc-zab*ybc) + ybd*(zab*xbc-xab*zbc) + zbd*(xab*ybc-yab*xbc))/((rab*rbc*rbd)*sine);
95    
96     angle = radian * asin(cosine);
97     dt = angle;
98     dt2 = dt * dt;
99     e = units.angunit * angles.copb[i] * dt2;
100     energies.eopb += e;
101     atom[ia].energy += e;
102     atom[ib].energy += e;
103     atom[ic].energy += e;
104     atom[id].energy += e;
105     if (minim_values.iprint)
106 wdelano 58 fprintf(pcmlogfile,"Oop: %-4d- %-4d- %-4d- %-4d %-8.3f %-8.3f = %-8.4f\n",ia,ib,ic,id,
107 tjod 3 angle, angles.copb[i],e);
108     }
109    
110     }
111     }
112     }
113     }
114    
115     void eopbend_wilson1()
116     {
117     int i,ia,ib,ic,id;
118     double e,angle,dot,cosine;
119     double dt,dt2;
120     double xia,yia,zia;
121     double xib,yib,zib;
122     double xic,yic,zic;
123     double xid,yid,zid;
124     double xab,yab,zab,rab2,rab;
125     double xbc,ybc,zbc,rbc2,rbc;
126     double xbd,ybd,zbd,rbd2,rbd;
127     double sine;
128     double deddt, terma, denom,termc;
129     double dedxia, dedxib, dedxic, dedxid;
130     double dedyia, dedyib, dedyic, dedyid;
131     double dedzia, dedzib, dedzic, dedzid;
132    
133     energies.eopb = 0.0;
134     for (i=0; i <= natom; i++)
135     {
136     deriv.deopb[i][0] = 0.0;
137     deriv.deopb[i][1] = 0.0;
138     deriv.deopb[i][2] = 0.0;
139     }
140    
141     for (i=0; i < angles.nang; i++)
142     {
143     if (angles.angin[i] == TRUE)
144     {
145     ia = angles.i13[i][0];
146     ib = angles.i13[i][1];
147     ic = angles.i13[i][2];
148     id = angles.i13[i][3];
149     if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use)
150     {
151     xia = atom[ia].x;
152     yia = atom[ia].y;
153     zia = atom[ia].z;
154     xib = atom[ib].x;
155     yib = atom[ib].y;
156     zib = atom[ib].z;
157     xic = atom[ic].x;
158     yic = atom[ic].y;
159     zic = atom[ic].z;
160     xid = atom[id].x;
161     yid = atom[id].y;
162     zid = atom[id].z;
163    
164     xab = xia - xib;
165     yab = yia - yib;
166     zab = zia - zib;
167     rab2 = (xab*xab+yab*yab+zab*zab);
168     rab = sqrt(rab2);
169    
170     xbc = xic - xib;
171     ybc = yic - yib;
172     zbc = zic - zib;
173     rbc2 = xbc*xbc+ybc*ybc+zbc*zbc;
174     rbc = sqrt(rbc2);
175    
176     if (rab2 != 0.0 && rbc2 != 0.0)
177     {
178     dot = xab*xbc + yab*ybc + zab*zbc;
179     cosine = dot /(rab*rbc);
180     if (cosine < -1.0)
181     cosine = -1.0;
182     if (cosine > 1.0)
183     cosine = 1.0;
184     sine = sqrt(1.00-cosine*cosine);
185    
186     xbd = atom[id].x - atom[ib].x;
187     ybd = atom[id].y - atom[ib].y;
188     zbd = atom[id].z - atom[ib].z;
189     rbd2 = xbd*xbd+ybd*ybd+zbd*zbd;
190     rbd = sqrt(rbd2);
191    
192     cosine = ( xbd*(yab*zbc-zab*ybc) + ybd*(zab*xbc-xab*zbc) + zbd*(xab*ybc-yab*xbc))/((rab*rbc*rbd)*sine);
193    
194     angle = radian * asin(cosine);
195     dt = angle;
196     dt2 = dt * dt;
197     e = units.angunit * angles.copb[i] * dt2;
198     energies.eopb += e;
199    
200     deddt = 2.0*units.angunit*angles.copb[i]*dt;
201     deddt *= radian;
202    
203     terma = xbd*(yab*zbc-zab*ybc) + ybd*(zab*xbc-xab*zbc) + zbd*(xab*ybc - yab*xbc);
204     termc = rab*rbc*rbd*sine;
205     denom = sqrt(1.00 - (terma*terma)/(rab2*rbd2*rbc2*sine*sine));
206    
207     dedxia = (-(terma*( (2.0*xab*dot*dot)/(rab2*rab2*rbc2) - (2.0*xbc*dot)/(rab2*rbc2))/(2.0*termc*sine*sine) )
208     - (xab*terma)/(rab2*termc) + (ybc*zbd-ybd*zbc)/termc ) / denom;
209     dedyia = (-(terma*( (2.0*yab*dot*dot)/(rab2*rab2*rbc2) - (2.0*ybc*dot)/(rab2*rbc2))/(2.0*termc*sine*sine) )
210     - (yab*terma)/(rab2*termc) + (zbc*xbd-zbd*xbc)/termc ) / denom;
211     dedzia = (-(terma*( (2.0*zab*dot*dot)/(rab2*rab2*rbc2) - (2.0*zbc*dot)/(rab2*rbc2))/(2.0*termc*sine*sine) )
212     - (zab*terma)/(rab2*termc) + (xbc*ybd-xbd*ybc)/termc ) / denom;
213    
214     dedxib = ( (xbd*terma)/(termc*rbd2) + (ybc*zab-yab*zbc+ybd*(zic-zia)+zbd*(yia-yic))/termc
215     +(xbc*terma)/(termc*rbc2) + (xab*terma)/(termc*rab2)
216     -( (-(2.0*xbc*dot*dot)/(rab2*rbc2*rbc2) - (2.0*(2.0*xib-xia-xic)*dot)/(rab2*rbc2)
217     -(2.0*xab*dot*dot)/(rab2*rab2*rbc2))*terma)/(2.0*termc*sine*sine) )/ denom;
218    
219     dedyib = ( (ybd*terma)/(termc*rbd2) + (xab*zbc-zab*xbc+xbd*(zab-zbc)+zbd*(xic-xia))/termc
220     +(ybc*terma)/(termc*rbc2) + (yab*terma)/(termc*rab2)
221     -( (-(2.0*ybc*dot*dot)/(rab2*rbc2*rbc2) - (2.0*(2.0*yib-yia-yic)*dot)/(rab2*rbc2)
222     -(2.0*yab*dot*dot)/(rab2*rab2*rbc2))*terma)/(2.0*termc*sine*sine) )/ denom;
223    
224     dedzib = ( (zbd*terma)/(termc*rbd2) + (xbc*yab-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic))/termc
225     +(zbc*terma)/(termc*rbc2) + (zab*terma)/(termc*rab2)
226     -( (-(2.0*zbc*dot*dot)/(rab2*rbc2*rbc2) - (2.0*(2.0*zib-zia-zic)*dot)/(rab2*rbc2)
227     -(2.0*zab*dot*dot)/(rab2*rab2*rbc2))*terma)/(2.0*termc*sine*sine) )/ denom;
228    
229     dedxic = (-(terma*( (2.0*xbc*dot*dot)/(rab2*rbc2*rbc2) - (2.0*xab*dot)/(rab2*rbc2) )/(2.0*termc*sine*sine) )
230     - (xbc*terma)/(rbc2*termc) + (zab*ybd-yab*zbd)/termc ) / denom;
231     dedyic = (-(terma*( (2.0*ybc*dot*dot)/(rab2*rbc2*rbc2) - (2.0*yab*dot)/(rab2*rbc2) )/(2.0*termc*sine*sine) )
232     - (ybc*terma)/(rbc2*termc) + (xab*zbd-zab*xbd)/termc ) / denom;
233     dedzic = (-(terma*( (2.0*zbc*dot*dot)/(rab2*rbc2*rbc2) - (2.0*zab*dot)/(rab2*rbc2) )/(2.0*termc*sine*sine) )
234     - (zbc*terma)/(rbc2*termc) + (yab*xbd-xab*ybd)/termc ) / denom;
235    
236     dedxid = ( (yab*zbc-zab*ybc)/termc - (xbd*terma)/(rbd2*termc) )/ denom;
237     dedyid = ( (zab*xbc-xab*zbc)/termc - (ybd*terma)/(rbd2*termc) )/ denom;
238     dedzid = ( (xab*ybc-yab*xbc)/termc - (zbd*terma)/(rbd2*termc) )/ denom;
239    
240     deriv.deopb[ia][0] += deddt*dedxia;
241     deriv.deopb[ia][1] += deddt*dedyia;
242     deriv.deopb[ia][2] += deddt*dedzia;
243    
244     deriv.deopb[ib][0] += deddt*dedxib;
245     deriv.deopb[ib][1] += deddt*dedyib;
246     deriv.deopb[ib][2] += deddt*dedzib;
247    
248     deriv.deopb[ic][0] += deddt*dedxic;
249     deriv.deopb[ic][1] += deddt*dedyic;
250     deriv.deopb[ic][2] += deddt*dedzic;
251    
252     deriv.deopb[id][0] += deddt*dedxid;
253     deriv.deopb[id][1] += deddt*dedyid;
254     deriv.deopb[id][2] += deddt*dedzid;
255    
256     }
257    
258     }
259     }
260     }
261     }
262    
263    
264     double term1aa(double rab,double rcross,double terma,double term1xa,double termc,
265     double termss2, double termcrab2)
266     {
267     double term1;
268    
269     term1 = rcross/termc - (term1xa*terma)/termss2 - (rab*terma)/(termcrab2);
270     return(term1);
271     }
272    
273     void eopbend_wilson2(int iatom)
274     {
275     int i,ia,ib,ic,id;
276     double e,angle,dot,cosine;
277     double dt,dt2;
278     double xia,yia,zia;
279     double xib,yib,zib;
280     double xic,yic,zic;
281     double xid,yid,zid;
282     double xab,yab,zab,rab2,rab;
283     double xbc,ybc,zbc,rbc2,rbc;
284     double xbd,ybd,zbd,rbd2,rbd;
285     double sine, rab2rbc2, term1, term2, term3, term4, term5, term6;
286     double term1xa, term1ya, term1za;
287     double term1xc, term1yc, term1zc;
288     double deddt, terma, denom,termc;
289     double dedxiaxia, dedxiaxib, dedxiaxic, dedxiaxid, dedxiayia, dedxiayib, dedxiayic, dedxiayid, dedxiazia, dedxiazib, dedxiazic,dedxiazid;
290     double dedxibxib, dedxibxic, dedxibxid, dedxibyia, dedxibyib, dedxibyic, dedxibyid, dedxibzia, dedxibzib, dedxibzic, dedxibzid;
291     double dedxicxic, dedxicxid, dedxicyia, dedxicyib, dedxicyic, dedxicyid, dedxiczia, dedxiczib, dedxiczic, dedxiczid;
292     double dedxidxid, dedxidyia, dedxidyib, dedxidyic, dedxidyid, dedxidzia, dedxidzib, dedxidzic, dedxidzid;
293     double dedyiayia, dedyiayib, dedyiayic, dedyiayid, dedyiazia, dedyiazib, dedyiazic, dedyiazid;
294     double dedyibyib, dedyibyic, dedyibyid, dedyibzia, dedyibzib, dedyibzic, dedyibzid;
295     double dedyicyic, dedyicyid, dedyiczia, dedyiczib, dedyiczic, dedyiczid;
296     double dedyidyid, dedyidzia, dedyidzib, dedyidzic, dedyidzid;
297     double dedziazia, dedziazib, dedziazic, dedziazid;
298     double dedzibzib, dedzibzic, dedzibzid;
299     double dedziczic, dedziczid;
300     double dedzidzid;
301     double termss, denom2,denom3,dot2,terma2,termss2;
302     double termc2, terms4, ddtdenom, ddt1denom, ddt2denom;
303     double xab2,yab2,zab2,xbc2,ybc2,zbc2,xbd2,ybd2,zbd2;
304     double yabzbc,yabxbd, ybczbd,ybdzab, xbdzbc, xabzbd, xabybc, xbcybd, xbczab;
305     double term1xaa, term1zca, termcrab2, termcrbc2, xiabc, yiabc, ziabc;
306     double rrab22rbc2, rrab2rbc22;
307     double zabybc, xabzbc, yabxbc;
308     double termcrbd22;
309    
310     for (i=0; i < angles.nang; i++)
311     {
312     if (angles.angin[i] == TRUE)
313     {
314     ia = angles.i13[i][0];
315     ib = angles.i13[i][1];
316     ic = angles.i13[i][2];
317     id = angles.i13[i][3];
318     if (iatom == ia || iatom == ib || iatom == ic || iatom == id)
319     {
320     xia = atom[ia].x;
321     yia = atom[ia].y;
322     zia = atom[ia].z;
323     xib = atom[ib].x;
324     yib = atom[ib].y;
325     zib = atom[ib].z;
326     xic = atom[ic].x;
327     yic = atom[ic].y;
328     zic = atom[ic].z;
329     xid = atom[id].x;
330     yid = atom[id].y;
331     zid = atom[id].z;
332    
333     xab = xia - xib;
334     yab = yia - yib;
335     zab = zia - zib;
336     rab2 = (xab*xab+yab*yab+zab*zab);
337     rab = sqrt(rab2);
338    
339     xbc = xic - xib;
340     ybc = yic - yib;
341     zbc = zic - zib;
342     rbc2 = xbc*xbc+ybc*ybc+zbc*zbc;
343     rbc = sqrt(rbc2);
344    
345     if (rab2 != 0.0 && rbc2 != 0.0)
346     {
347     dot = xab*xbc + yab*ybc + zab*zbc;
348     dot2 = dot*dot;
349     cosine = dot /(rab*rbc);
350     if (cosine < -1.0)
351     cosine = -1.0;
352     if (cosine > 1.0)
353     cosine = 1.0;
354     sine = sqrt(1.00-cosine*cosine);
355    
356     xbd = atom[id].x - atom[ib].x;
357     ybd = atom[id].y - atom[ib].y;
358     zbd = atom[id].z - atom[ib].z;
359     rbd2 = xbd*xbd+ybd*ybd+zbd*zbd;
360     rbd = sqrt(rbd2);
361    
362     xab2 = xab*xab;
363     yab2 = yab*yab;
364     zab2 = zab*zab;
365     xbc2 = xbc*xbc;
366     ybc2 = ybc*ybc;
367     zbc2 = zbc*zbc;
368     xbd2 = xbd*xbd;
369     ybd2 = ybd*ybd;
370     zbd2 = zbd*zbd;
371    
372     cosine = ( xbd*(yab*zbc-zab*ybc) + ybd*(zab*xbc-xab*zbc) + zbd*(xab*ybc-yab*xbc))/((rab*rbc*rbd)*sine);
373    
374     angle = asin(cosine);
375     dt = angle*radian;
376     dt2 = dt * dt;
377     e = units.angunit * angles.copb[i] * dt2;
378    
379     deddt = 2.0*units.angunit*angles.copb[i]*radian*radian;
380    
381     terma = xbd*(yab*zbc-zab*ybc) + ybd*(zab*xbc-xab*zbc) + zbd*(xab*ybc - yab*xbc);
382     terma2 = terma*terma;
383     termc = rab*rbc*rbd*sine;
384     termc2 = termc*termc;
385     denom = sqrt(1.00 - (terma2)/termc2);
386     denom2 = (1.00 - (terma2)/termc2);
387     denom3 = denom2*denom;
388     ddtdenom = deddt/denom2;
389     ddt1denom = (deddt*angle)/denom;
390     ddt2denom = (0.5*deddt*angle)/denom3;
391     rab2rbc2 = rab2*rbc2;
392     rrab22rbc2 = rab2*rab2rbc2;
393     rrab2rbc22 = rbc2*rab2rbc2;
394     termcrbd22 = termc*rbd2*rbd2;
395    
396     termss = termc*sine*sine;
397     terms4 = 4.0*termss*sine*sine;
398     termss2 = (2.0*termss);
399     term1xa = (2.0*xab*dot2)/(rrab22rbc2) - (2.0*xbc*dot)/rab2rbc2;
400     term1ya = (2.0*yab*dot2)/(rrab22rbc2) - (2.0*ybc*dot)/rab2rbc2;
401     term1za = (2.0*zab*dot2)/(rrab22rbc2) - (2.0*zbc*dot)/rab2rbc2;
402     term1xc = (2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xab*dot)/rab2rbc2;
403     term1yc = (2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yab*dot)/rab2rbc2;
404     term1zc = (2.0*zbc*dot2)/(rrab2rbc22) - (2.0*zab*dot)/rab2rbc2;
405     term1xaa = term1xa*terma;
406     term1zca = term1zc*terma;
407     termcrab2 = termc*rab2;
408     termcrbc2 = termc*rbc2;
409    
410     ybczbd = ybc*zbd-ybd*zbc;
411     ybdzab = ybd*zab-yab*zbd;
412     yabzbc = yab*zbc-ybc*zab;
413     yabxbd = yab*xbd-xab*ybd;
414     xbdzbc = xbd*zbc-xbc*zbd;
415     xabzbd = xab*zbd-zab*xbd;
416     xabybc = xab*ybc-yab*xbc;
417     xbcybd = xbc*ybd-ybc*xbd;
418     xbczab = xbc*zab-xab*zbc;
419     xabzbc = (xab*zbc-zab*xbc+xbd*(zia-zic)+zbd*(xic-xia));
420     yabxbc = (xbc*yab-ybc*xab+xbd*(ybc-yab)+ybd*(xia-xic));
421     zabybc = (ybc*zab-yab*zbc+ybd*(zic-zia)+zbd*(yia-yic));
422     xiabc = (2.0*xib-xia-xic);
423     yiabc = (2.0*yib-yia-yic);
424     ziabc = (2.0*zib-zia-zic);
425    
426     term1 = term1aa(xab, ybczbd, terma, term1xa, termc, termss2, termcrab2);
427     term3 = (8.0*xab*xbc*dot)/(rrab22rbc2)- (2.0*xbc2)/rab2rbc2
428     -(8.0*xab2*dot2)/(rab2*rrab22rbc2) + (2.0*dot2)/(rrab22rbc2);
429     term2 = - (term1xa*(ybczbd))/(termss)
430     - (2.0*xab*(ybczbd))/(termcrab2)
431     + (3.0*terma*term1xa*term1xa)/(terms4)
432     - (terma*term3)/termss2
433     + (xab*term1xaa)/(termss*rab2)
434     + (3.0*xab2*terma)/(rab2*termcrab2) - terma/(termcrab2);
435     term4 = (- (2.0*terma*(ybczbd))/(termc2)
436     + (term1xaa*terma)/(termc*termss)
437     + (2.0*xab*terma2)/(termc2*rab2) )
438     *( (ybczbd)/termc - (term1xaa)/termss2
439     - (xab*terma)/(termcrab2));
440     dedxiaxia = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
441    
442     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
443     term3 = (4.0*xab*xbc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*xbc2*dot)/(rrab2rbc22) - (2.0*xiabc*xbc)/rab2rbc2
444     +(4.0*xab*xiabc*dot)/(rrab22rbc2) - (4.0*xab*xbc*dot)/(rrab22rbc2) + (2.0*dot)/rab2rbc2
445     +(8.0*xab2*dot2)/(rrab22rbc2*rab2) - (2.0*dot2)/(rrab22rbc2);
446     term1 = ( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2))
447     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
448     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
449     term2 = (xbd*(ybczbd))/(termc*rbd2)
450     -(xbd*term1xaa)/(2.0*termss*rbd2)
451     -(xab*xbd*terma)/(rab2*rbd2*termc)
452     -(term1xa*zabybc)/termss2
453     -(xab*zabybc)/(termcrab2)
454     -(term5*(ybczbd))/termss2
455     +(xbc*(ybczbd))/(termcrbc2)
456     +(xab*(ybczbd))/(termcrab2)
457     +(3.0*term5*term1xaa)/(terms4)
458     -(term3*terma)/termss2
459     +(xab*term5*terma)/(2.0*termss*rab2)
460     -(xbc*term1xaa)/(2.0*termss*rbc2)
461     -(xab*term1xaa)/(2.0*termss*rab2)
462     -(xab*xbc*terma)/(termcrab2*rbc2)
463     -(3.0*xab2*terma)/(termcrab2*rab2)
464     +terma/(termcrab2);
465     term4 = (-(2.0*xbd*terma2)/(termc2*rbd2) - (2.0*zabybc*terma)/(termc2)
466     +(term5*terma2)/(termc*termss) - (2.0*xbc*terma2)/(termc2*rbc2) - (2.0*xab*terma2)/(termc2*rab2))
467     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
468     dedxiaxib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
469    
470     term5 = (4.0*xbc2*dot)/(rrab2rbc22) - (4.0*xab*xbc*dot2)/(rab2rbc2*rab2rbc2) - (2.0*xab*xbc)/rab2rbc2
471     +(4.0*xab2*dot)/(rrab22rbc2) - 2.0*dot/rab2rbc2;
472     term1 = ( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2))
473     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
474     term2 = -(term1xa*(ybdzab))/termss2
475     -(xab*(ybdzab))/(termcrab2)
476     -(term1xc*(ybczbd))/termss2
477     -(xbc*(ybczbd))/(rbc2*termc)
478     +(3.0*term1xa*term1xc*terma)/(terms4)
479     -(term5*terma)/termss2
480     +(xab*term1xc*terma)/(2.0*rab2*termss)
481     +(xbc*term1xaa)/(2.0*rbc2*termss)
482     +(xab*xbc*terma)/(rab2*rbc2*termc);
483     term4 = (-(2.0*(ybdzab)*terma)/(termc2) + (term1xc*terma2)/(termc*termss)
484     +(2.0*xbc*terma2)/(termc2))
485     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
486     dedxiaxic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
487    
488     term1 = (-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc)*((ybczbd)/termc
489     -(term1xaa)/termss2 - (xab*terma)/(termcrab2));
490     term2 = -(xbd*(ybczbd))/(termc*rbd2)
491     +(xbd*term1xaa)/(2.0*termss*rbd2)
492     +(xab*xbd*terma)/(rab2*rbd2*termc)
493     -((yabzbc)*term1xa)/termss2
494     -(xab*(yabzbc))/(termcrab2);
495     term4 = ( (2.0*xbd*terma2)/(termc2*rbd2) - (2.0*(yabzbc)*terma)/(termc2))
496     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
497     dedxiaxid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
498    
499     term5 = (4.0*xbc*yab*dot)/(rrab22rbc2) + (4.0*xab*ybc*dot)/(rrab22rbc2) - (8.0*xab*yab*dot2)/(rrab22rbc2)
500     -(2.0*xbc*ybc)/rab2rbc2;
501     term1 = ( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2))
502     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
503     term2 = -(term1xa*(xbdzbc))/termss2
504     -(xab*(xbdzbc))/(termcrab2)
505     -(term1ya*(ybczbd))/termss2
506     -(yab*(ybczbd))/(termcrab2)
507     +(3.0*term1xa*term1ya*terma)/(terms4)
508     -(term5*terma)/termss2
509     +(yab*term1xaa)/(2.0*rab2*termss)
510     +(xab*term1ya*terma)/(2.0*rab2*termss)
511     +(3.0*xab*yab*terma)/(rab2*termcrab2);
512     term4 = (-(2.0*(xbdzbc)*terma)/(termc2) + (term1ya*terma2)/(termc*termss)
513     +(2.0*yab*terma2)/(termc2*rab2) )
514     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
515     dedxiayia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
516    
517     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
518     term3 = (4.0*xab*ybc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*xbc*ybc*dot)/(rrab2rbc22) - (2.0*yiabc*xbc)/rab2rbc2
519     +(4.0*xab*yiabc*dot)/(rrab22rbc2) - (4.0*yab*xbc*dot)/(rrab22rbc2)
520     +(8.0*xab*yab*dot2)/(rrab22rbc2*rab2);
521     term1 = ( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2))
522     *( (ybd*terma)/(termc*rbd2) + (xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))/termc
523     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
524     term2 = (ybd*(ybczbd))/(termc*rbd2)
525     -(ybd*term1xaa)/(2.0*termss*rbd2)
526     -(xab*ybd*terma)/(rab2*rbd2*termc)
527     +(zbc-zbd)/termc
528     -(term1xa*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/termss2
529     -(xab*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/(termcrab2)
530     -(term5*(ybczbd))/termss2
531     +(ybc*(ybczbd))/(termcrbc2)
532     +(yab*(ybczbd))/(termcrab2)
533     +(3.0*term5*term1xaa)/(terms4)
534     -(term3*terma)/termss2
535     -(ybc*term1xaa)/(2.0*termss*rbc2)
536     -(yab*term1xaa)/(2.0*termss*rab2)
537     +(xab*term5*terma)/(2.0*termss*rab2)
538     -(xab*ybc*terma)/(termcrab2*rbc2)
539     -(3.0*xab*yab*terma)/(termcrab2*rab2);
540     term4 = (-(2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))*terma)/(termc2)
541     +(term5*terma2)/(termc*termss) - (2.0*ybc*terma2)/(termc2*rbc2) - (2.0*yab*terma2)/(termc2*rab2))
542     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
543     dedxiayib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
544    
545     term5 = (4.0*xbc*ybc*dot)/(rrab2rbc22) - (4.0*xab*ybc*dot2)/(rab2rbc2*rab2rbc2) - (2.0*xbc*yab)/rab2rbc2
546     +(4.0*xab*yab*dot)/(rrab22rbc2);
547     term1 = ( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2))
548     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
549     term2 = (zbd)/termc
550     -(term1xa*(xabzbd))/termss2
551     -(xab*(xabzbd))/(termcrab2)
552     -(term1yc*(ybczbd))/termss2
553     -(ybc*(ybczbd))/(rbc2*termc)
554     +(3.0*term1yc*term1xaa)/(terms4)
555     -(term5*terma)/termss2
556     +(xab*term1yc*terma)/(2.0*rab2*termss)
557     +(ybc*term1xaa)/(2.0*termss*rbc2)
558     +(xab*ybc*terma)/(termcrab2*rbc2);
559     term4 = (-(2.0*(xabzbd)*terma)/(termc2) + (term1yc*terma2)/(termc*termss)
560     +(2.0*ybc*terma2)/(termc2*rbc2) )
561     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
562     dedxiayic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
563    
564     term1 = (-(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*( (ybczbd)/termc
565     -(term1xaa)/termss2 - (xab*terma)/(termcrab2));
566     term2 = -(ybd*(ybczbd))/(termc*rbd2)
567     +(ybd*term1xaa)/(2.0*termss*rbd2)
568     +(xab*ybd*terma)/(rab2*rbd2*termc)
569     -((xbczab)*term1xa)/termss2
570     -(zbc)/termc
571     -(xab*(xbczab))/(termcrab2);
572     term4 = ( (2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xbczab)*terma)/(termc2) )
573     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
574     dedxiayid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
575    
576     term5 = (4.0*xab*zbc*dot)/(rrab22rbc2) + (4.0*xbc*zab*dot)/(rrab22rbc2) - (2.0*xbc*zbc)/rab2rbc2
577     -(8.0*xab*zab*dot2)/(rab2*rrab22rbc2);
578     term1 = ( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2))
579     *( (xbcybd)/termc - (term1za*terma)/termss2 - (zab*terma)/(termcrab2));
580     term2 = -(term1xa*(xbcybd))/termss2
581     -(xab*(xbcybd))/(termcrab2)
582     -(zab*(ybczbd))/(termcrab2)
583     -(term1za*(ybczbd))/termss2
584     +(zab*term1xaa)/(2.0*rab2*termss)
585     +(3.0*xab*zab*terma)/(rab2*termcrab2)
586     -(term5*terma)/termss2
587     +(3.0*term1xa*term1za*terma)/(terms4)
588     +(xab*term1za*terma)/(2.0*rab2*termss);
589     term4 = (-(2.0*(xbcybd)*terma)/(termc2) + (2.0*zab*terma2)/(termc2*rab2)
590     +(term1za*terma2)/(termc*termss))
591     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
592     dedxiazia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
593    
594     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
595     term3 = (4.0*xab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*xbc*zbc*dot)/(rrab2rbc22)
596     -(2.0*xbc*ziabc)/rab2rbc2 - (4.0*xbc*zab*dot)/(rrab22rbc2)
597     +(4.0*xab*ziabc*dot)/(rrab22rbc2) + (8.0*xab*zab*dot2)/(rrab22rbc2*rab2);
598     term1 = ( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2))
599     *( (zbd*terma)/(termc*rbd2) + (xbc*yab-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic))/termc
600     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term5*terma)/termss2);
601     term2 = (zbd*(ybczbd))/(termc*rbd2)
602     -(zbd*term1xaa)/(2.0*termss*rbd2)
603     -(xab*zbd*terma)/(rab2*rbd2*termc)
604     -(term1xa*(xbc*yab-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic)))/termss2
605     -(xab*(xbc*yab-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic)))/(termcrab2)
606     +(ybd-ybc)/termc
607     +(zbc*(ybczbd))/(termcrbc2)
608     +(zab*(ybczbd))/(termcrab2)
609     -(term5*(ybczbd))/termss2
610     -(zbc*term1xaa)/(2.0*termss*rbc2)
611     -(zab*term1xaa)/(2.0*termss*rab2)
612     -(xab*zbc*terma)/(termcrab2*rbc2)
613     -(3.0*xab*zab*terma)/(termcrab2*rab2)
614     -(term3*terma)/termss2
615     +(3.0*term5*term1xaa)/(terms4)
616     +(xab*term5*terma)/(2.0*termss*rab2);
617     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xbc*yab-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic))*terma)/(termc2)
618     +(term5*terma2)/(termc*termss) - (2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2))
619     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
620     dedxiazib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
621    
622     term5 = (4.0*xbc*zbc*dot)/(rrab2rbc22) - (4.0*xab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (2.0*xbc*zab)/rab2rbc2
623     + (4.0*xab*zab*dot)/(rrab22rbc2);
624     term1 = ( (yabxbd)/termc - (term1zca)/termss2 - (zbc*terma)/(termcrbc2))
625     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
626     term2 = -(term1xa*(yabxbd))/termss2
627     -(xab*(yabxbd))/(termcrab2)
628     -ybd/termc
629     -(term1zc*(ybczbd))/termss2
630     -(zbc*(ybczbd))/(termcrbc2)
631     +(3.0*term1zc*term1xaa)/(terms4)
632     -(term5*terma)/termss2
633     +(xab*term1zca)/(2.0*rab2*termss)
634     +(zbc*term1xaa)/(2.0*rbc2*termss)
635     +(xab*zbc*terma)/(rab2*rbc2*termc);
636     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
637     +(2.0*zbc*terma2)/(rbc2*termc2))
638     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
639     dedxiazic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
640    
641     term1 = (-(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*( (ybczbd)/termc
642     -(term1xaa)/termss2 - (xab*terma)/(termcrab2));
643     term2 = -(zbd*(ybczbd))/(rbd2*termc)
644     +(term1xa*zbd*terma)/(2.0*termss*rbd2)
645     +(xab*zbd*terma)/(rab2*rbd2*termc)
646     -(term1xa*(xabybc))/termss2
647     -(xab*(xabybc))/(termcrab2)
648     + ybc/termc;
649     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
650     *( (ybczbd)/termc - (term1xaa)/termss2 - (xab*terma)/(termcrab2));
651     dedxiazid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
652    
653     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
654     term3 = -(8.0*xbc2*dot2)/(rrab2rbc22*rbc2) - (8.0*xbc*dot*xiabc)/(rrab2rbc22)
655     -(8.0*xab*xbc*dot2)/(rab2rbc2*rab2rbc2) + (2.0*dot2)/(rrab2rbc22)
656     -(2.0*xiabc*xiabc)/rab2rbc2
657     -(8.0*xab*xiabc*dot)/(rrab22rbc2) - (8.0*xab2*dot2)/(rrab22rbc2*rab2)
658     -(4.0*dot)/rab2rbc2 + (2.0*dot2)/(rrab22rbc2);
659     term1 = (xbd*terma)/(termc*rbd2) + zabybc/termc
660     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2);
661     term2 = (3.0*xbd2*terma)/(termcrbd22)
662     +(2.0*xbd*zabybc)/(termc*rbd2)
663     -(xbd*term5*terma)/(termss*rbd2)
664     +(2.0*xbc*xbd*terma)/(termc*rbd2*rbc2)
665     +(2.0*xab*xbd*terma)/(termc*rbd2*rab2)
666     -terma/(termc*rbd2)
667     -(term5*zabybc)/(termc*sine)
668     +(2.0*xbc*zabybc)/(termcrbc2)
669     +(2.0*xab*zabybc)/(termcrab2)
670     +(3.0*term5*term5*terma)/(terms4)
671     -(term3*terma)/termss2
672     -(xbc*term5*terma)/(termss*rbc2)
673     -(xab*term5*terma)/(termss*rab2)
674     +(3.0*xbc2*terma)/(termcrbc2*rbc2)
675     +(2.0*xab*xbc*terma)/(termcrab2*rbc2)
676     -terma/(termcrbc2)
677     +(3.0*xab2*terma)/(termcrab2*rab2)
678     -terma/(termcrab2);
679     term4 = (-(2.0*xbd*terma2)/(termc2*rbd2) - (2.0*zabybc*terma)/(termc2)
680     +(term5*terma2)/(termc*termss) - (2.0*xbc*terma2)/(termc2*rbc2) - (2.0*xab*terma2)/(termc2*rab2))
681     *( (xbd*terma)/(termc*rbd2) + zabybc/termc - (term5*terma)/(termc*2.0*sine*sine)
682     +(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2)) ;
683     dedxibxib = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
684    
685     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
686     term3 = (8.0*xbc2*dot2)/(rrab2rbc22*rbc2) - (4.0*xab*xbc*dot)/(rrab2rbc22)
687     +(4.0*xbc*dot*xiabc)/(rrab2rbc22) + (4.0*xab*xbc*dot2)/(rab2rbc2*rab2rbc2)
688     -(2.0*dot2)/(rrab2rbc22) - (2.0*xab*xiabc)/rab2rbc2 - (4.0*xab2*dot)/(rrab22rbc2)
689     +(2.0*dot)/rab2rbc2;
690     term1 = ( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2))
691     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
692     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
693     term2 = (xbd*(ybdzab))/(termc*rbd2)
694     -(xbd*term1xc*terma)/(2.0*termss*rbd2)
695     -(xbc*xbd*terma)/(termc*rbd2*rbc2)
696     -(term5*(ybdzab))/(termc*2.0*sine*sine)
697     +(xbc*(ybdzab))/(termcrbc2)
698     +(xab*(ybdzab))/(termcrab2)
699     -(term1xc*zabybc)/(termc*2.0*sine*sine)
700     -(xbc*zabybc)/(termcrbc2)
701     +(3.0*term1xc*term5*terma)/(terms4)
702     -(term3*terma)/termss2
703     -(xbc*term1xc*terma)/(2.0*termss*rbc2)
704     -(xab*term1xc*terma)/(2.0*termss*rab2)
705     +(xbc*term5*terma)/(2.0*termss*rbc2)
706     -(3.0*xbc2*terma)/(termcrbc2*rbc2)
707     -(xab*xbc*terma)/(termcrab2*rbc2)
708     +terma/(termcrbc2);
709     term4 = (-(2.0*(ybdzab)*terma)/(termc2) + (term1xc*terma2)/(termc*termss)
710     +(2.0*xbc*terma2)/(termc2*rbc2))
711     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
712     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
713     dedxibxic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
714    
715     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
716     term1 = (-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc)*( (xbd*terma)/(termc*rbd2)
717     +zabybc/termc
718     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
719     term2 = -(3.0*xbd2*terma)/(termcrbd22)
720     +(xbd*(yabzbc))/(termc*rbd2)
721     -(xbd*zabybc)/(termc*rbd2)
722     +(xbd*term5*terma)/(2.0*termss*rbd2)
723     -(xbd*xbc*terma)/(termcrbc2*rbd2)
724     -(xab*xbd*terma)/(termcrab2*rbd2)
725     +terma/(termc*rbd2)
726     -(term5*(yabzbc))/termss2
727     +(xbc*(yabzbc))/(termcrbc2)
728     +(xab*(yabzbc))/(termcrab2);
729     term4 = ( (2.0*xbd*terma2)/(termc2*rbd2) - (2.0*(yabzbc)*terma)/(termc2))
730     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
731     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
732     dedxibxid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
733    
734     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
735     term3 = (8.0*xab*yab*dot2)/(rrab22rbc2*rab2) - (4.0*ybc*xab*dot)/(rrab22rbc2)
736     +(4.0*yab*dot*xiabc)/(rrab22rbc2) + (4.0*yab*xbc*dot2)/(rab2rbc2*rab2rbc2)
737     -(2.0*ybc*xiabc)/(rab2rbc2) - (4.0*xbc*ybc*dot)/(rrab2rbc22);
738     term1 = ( (xbd*terma)/(termc*rbd2) + (ybc*zab-yab*zbc+ybd*(zic-zia)+zbd*(yia-yic) )/termc
739     -(term5*terma)/termss2+(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2))
740     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
741     term2 = (xbd*(xbdzbc))/(termc*rbd2)
742     -(xbd*term1ya*terma)/(2.0*termss*rbd2)
743     -(xbd*yab*terma)/(termcrab2*rbd2)
744     +(zbd-zbc)/termc
745     -(term5*(xbdzbc))/termss2
746     +(xbc*(xbdzbc))/(termcrbc2)
747     +(xab*(xbdzbc))/(termcrab2)
748     -(term1ya*zabybc)/termss2
749     -(yab*zabybc)/(termcrab2)
750     +(3.0*term1ya*term5*terma)/(terms4)
751     -(term3*terma)/termss2
752     +(yab*term5*terma)/(2.0*termss*rab2)
753     -(xbc*term1ya*terma)/(2.0*termss*rbc2)
754     -(xab*term1ya*terma)/(2.0*termss*rab2)
755     -(3.0*xab*yab*terma)/(termcrab2*rab2)
756     -(yab*xbc*terma)/(termcrbc2*rab2);
757     term4 = (-(2.0*(xbdzbc)*terma)/(termc2) + (term1ya*terma2)/(termc*termss)
758     +(2.0*yab*terma2)/(termc2*rab2))
759     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
760     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
761     dedxibyia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
762    
763     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
764     term6 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
765     term3 = -(8.0*xbd*ybc*dot2)/(rrab2rbc22*rbc2) - (4.0*xbd*dot*(yiabc))/(rrab2rbc22)
766     -(4.0*ybc*dot*xiabc)/(rrab2rbc22) - (4.0*xbd*yab*dot2)/(rab2rbc2*rab2rbc2)
767     -(4.0*xab*ybc*dot2)/(rab2rbc2*rab2rbc2) - (2.0*xiabc*yiabc)/rab2rbc2
768     -(4.0*yab*dot*xiabc)/(rrab22rbc2) - (4.0*xab*yiabc*dot)/(rrab22rbc2)
769     -(8.0*xab*yab*dot2)/(rrab22rbc2*rab2);
770     term1 = ( (xbd*terma)/(termc*rbd2) + (ybc*zab-yab*zbc+ybd*(zic-zia)+zbd*(yia-yic) )/termc
771     -(term5*terma)/termss2+(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2))
772     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
773     -(term6*terma)/termss2+(ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
774     term2 = (3.0*xbd*ybd*terma)/(termcrbd22)
775     +(xbd*xabzbc)/(termc*rbd2)
776     +(ybd*zabybc)/(termc*rbd2)
777     -(ybd*term5*terma)/(2.0*termss*rbd2)
778     -(xbd*term6*terma)/(2.0*termss*rbd2)
779     +(xbd*ybc*terma)/(termcrbc2*rbd2)
780     +(xbd*ybd*terma)/(termcrbc2*rbd2)
781     +(xbd*yab*terma)/(termcrab2*rbd2)
782     +(xab*ybd*terma)/(termcrab2*rbd2)
783     -(term5*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/termss2
784     +(xbc*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/(termcrbc2)
785     +(xab*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/(termcrab2)
786     -(term6*zabybc)/termss2
787     +(ybc*zabybc)/(termcrbc2)
788     +(yab*zabybc)/(termcrab2)
789     +(3.0*term5*term6*terma)/(terms4)
790     -(term3*terma)/termss2
791     -(ybc*term5*terma)/(2.0*termss*rbc2)
792     -(yab*term5*terma)/(2.0*termss*rab2)
793     -(xbc*term6*terma)/(2.0*termss*rbc2)
794     -(xab*term6*terma)/(2.0*termss*rab2)
795     +(3.0*xbc*ybc*terma)/(termcrbc2*rbc2)
796     +(xbc*yab*terma)/(termcrab2*rbc2)
797     +(xab*ybc*terma)/(termcrab2*rbc2)
798     +(3.0*xab*yab*terma)/(termcrab2*rab2);
799     term4 = (-(2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))*terma)/(termc2)
800     +(term6*terma2)/(termc*termss) - (2.0*ybc*terma2)/(termc2*rbc2) - (2.0*yab*terma2)/(termc2*rab2))
801     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
802     -(term5*terma)/termss2 + (xbd*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
803     dedxibyib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
804    
805     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
806     term3 = (8.0*xbc*ybc*dot2)/(rrab2rbc22*rbc2) - (4.0*xbc*yab*dot)/(rrab2rbc22)
807     +(4.0*ybc*dot*xiabc)/(rrab2rbc22) + (4.0*xab*ybc*dot2)/(rab2rbc2*rab2rbc2)
808     -(2.0*yab*xiabc)/(rab2rbc2) - (4.0*xab*yab*dot)/(rrab22rbc2);
809     term1 = ( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2))
810     *( (xbd*terma)/(termc*rbd2) + (ybc*zab-yab*zbc+ybd*(zic-zia)+zbd*(yia-yic) )/termc
811     -(term5*terma)/termss2+(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2) );
812     term2 = (xbd*(xabzbd))/(termc*rbd2)
813     -(xbd*term1yc*terma)/(2.0*termss*rbd2)
814     -(xbd*ybc*terma)/(termcrbc2*rbd2)
815     +(zab-zbd)/termc
816     -(term5*(xabzbd))/termss2
817     +(xbc*(xabzbd))/(termcrbc2)
818     +(xab*(xabzbd))/(termcrab2)
819     -(term1yc*zabybc)/termss2
820     -(ybc*zabybc)/(termcrbc2)
821     +(3.0*term1yc*term5*terma)/(terms4)
822     -(term3*terma)/termss2
823     -(xbc*term1yc*terma)/(2.0*termss*rbc2)
824     -(xab*term1yc*terma)/(2.0*termss*rab2)
825     +(ybc*term5*terma)/(2.0*termss*rbc2)
826     -(3.0*xbc*ybc*terma)/(termcrbc2*rbc2)
827     -(xab*ybc*terma)/(termcrbc2*rab2);
828     term4 = (-(2.0*(xabzbd)*terma)/(termc2) + (term1yc*terma2)/(termc*termss)
829     +(2.0*ybc*terma2)/(termc2*rbc2))
830     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
831     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
832     dedxibyic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
833    
834     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
835     term1 = (-(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*( (xbd*terma)/(termc*rbd2)
836     +zabybc/termc - (term5*terma)/termss2
837     +(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
838     term2 = -(3.0*xbd*ybd*terma)/(termcrbd22)
839     +(xbd*(xbczab))/(termc*rbd2)
840     -(ybd*zabybc)/(termc*rbd2)
841     +(ybd*term5*terma)/(2.0*termss*rbd2)
842     -(xbc*ybd*terma)/(termcrbc2*rbd2)
843     -(xab*ybd*terma)/(termcrab2*rbd2)
844     -((xbczab)*term5)/termss2
845     +(xbc*(xbczab))/(termcrbc2)
846     +(zbc-zab)/termc
847     +(xab*(xbczab))/(termcrab2);
848     term4 = ( (2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xbczab)*terma)/(termc2))
849     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
850     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
851     dedxibyid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
852    
853     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
854     term3 = (8.0*xab*zab*dot2)/(rrab22rbc2*rab2) - (4.0*zbc*xab*dot)/(rrab22rbc2)
855     +(4.0*zab*dot*xiabc)/(rrab22rbc2) + (4.0*zab*xbc*dot2)/(rab2rbc2*rab2rbc2)
856     -(2.0*zbc*xiabc)/(rab2rbc2) - (4.0*xbc*zbc*dot)/(rrab2rbc22);
857     term1 = ( (xbd*terma)/(termc*rbd2) + zabybc/termc
858     -(term5*terma)/termss2+(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2))
859     *( (xbcybd)/termc - (term1za*terma)/termss2 - (zab*terma)/(termcrab2));
860     term2 = (xbd*(xbcybd))/(termc*rbd2)
861     -(xbd*zab*terma)/(termcrab2*rbd2)
862     -(xbd*term1za*terma)/(2.0*termss*rbd2)
863     -(term5*(xbcybd))/termss2
864     +(xbc*(xbcybd))/(termcrbc2)
865     +(xab*(xbcybd))/(termcrab2)
866     +(ybc-ybd)/termc
867     -(zab*zabybc)/(termcrab2)
868     -(term1za*zabybc)/termss2
869     +(zab*term5*terma)/(2.0*termss*rab2)
870     -(xbc*zab*terma)/(termcrbc2*rab2)
871     -(3.0*xab*zab*terma)/(termcrab2*rab2)
872     -(term3*terma)/termss2
873     +(3.0*term1za*term5*terma)/(terms4)
874     -(xbc*term1za*terma)/(2.0*termss*rbc2)
875     -(xab*term1za*terma)/(2.0*termss*rab2);
876     term4 = (-(2.0*(xbcybd)*terma)/(termc2) + (term1za*terma2)/(termc*termss)
877     +(2.0*zab*terma2)/(termc2*rab2))
878     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
879     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
880     dedxibzia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
881    
882     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
883     term6 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
884     term3 = -(8.0*xbc*zbc*dot2)/(rrab2rbc22*rbc2) - (4.0*xbc*dot*(ziabc))/(rrab2rbc22)
885     -(4.0*zbc*dot*xiabc)/(rrab2rbc22) - (4.0*xbc*zab*dot2)/(rab2rbc2*rab2rbc2)
886     -(4.0*xab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (2.0*xiabc*ziabc)/rab2rbc2
887     -(4.0*zab*dot*xiabc)/(rrab22rbc2) - (4.0*xab*ziabc*dot)/(rrab22rbc2)
888     -(8.0*xab*zab*dot2)/(rrab22rbc2*rab2);
889     term1 = ( (xbd*terma)/(termc*rbd2) + zabybc/termc
890     -(term5*terma)/termss2+(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2))
891     *( (zbd*terma)/(termc*rbd2) + (xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
892     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2)-(term6*terma)/termss2);
893     term2 = (3.0*xbd*zbd*terma)/(termcrbd22)
894     +(xbd*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termc*rbd2)
895     +(zbd*zabybc)/(termc*rbd2)
896     +(xbd*zbc*terma)/(termcrbc2*rbd2)
897     +(xbd*zab*terma)/(termcrab2*rbd2)
898     -(xbd*term6*terma)/(2.0*termss*rbd2)
899     -(zbd*term5*terma)/(2.0*termss*rbd2)
900     +(xbc*zbd*terma)/(termcrbc2*rbd2)
901     +(xab*zbd*terma)/(termcrab2*rbd2)
902     -(term5*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/termss2
903     +(xbc*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrbc2)
904     +(xab*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrab2)
905     +(zbc*zabybc)/(termcrbc2)
906     +(zab*zabybc)/(termcrab2)
907     -(term6*zabybc)/termss2
908     -(zbc*term5*terma)/(2.0*termss*rbc2)
909     -(zab*term5*terma)/(2.0*termss*rab2)
910     +(3.0*xbc*zbc*terma)/(termcrbc2*rbc2)
911     +(xbc*zab*terma)/(termcrab2*rbc2)
912     +(xab*zbc*terma)/(termcrab2*rbc2)
913     +(3.0*xab*zab*terma)/(termcrab2*rab2)
914     -(term3*terma)/termss2
915     +(3.0*term5*term6*terma)/(terms4)
916     -(xbc*term6*terma)/(2.0*termss*rbc2)
917     -(xab*term6*terma)/(2.0*termss*rab2);
918     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))*terma)/(termc2)
919     +(term6*terma2)/(termc*termss) - (2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2))
920     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
921     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
922     dedxibzib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
923    
924     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
925     term3 = (8.0*xbc*zbc*dot2)/(rrab2rbc22*rbc2) - (4.0*xbc*zab*dot)/(rrab2rbc22)
926     +(4.0*zbc*dot*xiabc)/(rrab2rbc22) + (4.0*xab*zbc*dot2)/(rab2rbc2*rab2rbc2)
927     -(2.0*zab*xiabc)/(rab2rbc2) - (4.0*xab*zab*dot)/(rrab22rbc2);
928     term1 = ( (yabxbd)/termc - (term1zca)/termss2 - (zbc*terma)/(termcrbc2))
929     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
930     -(term5*terma)/termss2+(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
931     term2 = (xbd*(yabxbd))/(termc*rbd2)
932     -(xbd*term1zca)/(2.0*termss*rbd2)
933     -(xbd*zbc*terma)/(termcrbc2*rbd2)
934     -(term5*(yabxbd))/termss2
935     +(xbc*(yabxbd))/(termcrbc2)
936     +(xab*(yabxbd))/(termcrab2)
937     +(ybd-yab)/termc
938     -(term1zc*zabybc)/termss2
939     -(zbc*zabybc)/(termcrbc2)
940     +(3.0*term1zc*term5*terma)/(terms4)
941     -(term3*terma)/termss2
942     -(xbc*term1zca)/(2.0*termss*rbc2)
943     -(xab*term1zca)/(2.0*termss*rab2)
944     +(zbc*term5*terma)/(2.0*termss*rbc2)
945     -(3.0*xbc*zbc*terma)/(termcrbc2*rbc2)
946     -(xab*zbc*terma)/(termcrbc2*rab2);
947     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
948     +(2.0*zbc*terma2)/(termc2*rbc2))
949     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
950     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
951     dedxibzic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
952    
953     term5 = -(2.0*xbc*dot2)/(rrab2rbc22) - (2.0*xiabc*dot)/rab2rbc2 - (2.0*xab*dot2)/(rrab22rbc2);
954     term1 = (-(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*( (xbd*terma)/(termc*rbd2)
955     +zabybc/termc - (term5*terma)/termss2
956     +(xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
957     term2 = -(3.0*xbd*zbd*terma)/(termcrbd22)
958     +(xbd*(xabybc))/(termc*rbd2)
959     -(zbd*zabybc)/(termc*rbd2)
960     +(zbd*term5*terma)/(2.0*termss*rbd2)
961     -(xbc*zbd*terma)/(termcrbc2*rbd2)
962     -(xab*zbd*terma)/(termcrab2*rbd2)
963     -((xabybc)*term5)/termss2
964     +(xbc*(xabybc))/(termcrbc2)
965     +(yab-ybc)/termc
966     +(xab*(xabybc))/(termcrab2);
967     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
968     *( (xbd*terma)/(termc*rbd2) + zabybc/termc
969     -(term5*terma)/termss2 + (xbc*terma)/(termcrbc2) + (xab*terma)/(termcrab2));
970     dedxibzid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
971    
972     term3 = (8.0*xab*xbc*dot)/(rrab2rbc22)- (2.0*xab2)/rab2rbc2 - (8.0*xbc2*dot2)/(rbc2*rrab2rbc22)
973     + (2.0*dot2)/(rrab2rbc22);
974     term1 = term1aa(xbc, ybdzab, terma, term1xc, termc, termss2, termcrbc2);
975     term2 = - (term1xc*(ybdzab))/(termss)
976     - (2.0*xbc*(ybdzab))/(termcrbc2)
977     + (3.0*terma*term1xc*term1xc)/(terms4)
978     - (terma*term3)/termss2
979     + (xbc*term1xc*terma)/(termss*rbc2)
980     + (3.0*xbc2*terma)/(rbc2*rbc2*termc)
981     - terma/(rbc2*termc);
982     term4 = (-(2.0*terma*(ybdzab))/(termc2)
983     +(term1xc*terma2)/(termc*termss)
984     +(2.0*xbc*terma2)/(rbc2*termc2) )
985     *((ybdzab)/termc - (term1xc*terma)/termss2
986     -(xbc*terma)/(rbc2*termc));
987     dedxicxic = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
988    
989    
990     term1 = (-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc)*( (ybdzab)/termc
991     -(term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2));
992     term2 = -(xbd*(ybdzab))/(termc*rbd2)
993     +(xbd*term1xc*terma)/(2.0*termss*rbd2)
994     +(xbc*xbd*terma)/(rbc2*rbd2*termc)
995     -((yabzbc)*term1xc)/termss2
996     -(xbc*(yabzbc))/(rbc2*termc);
997     term4 = ( (2.0*xbd*terma2)/(termc2*rbd2) - (2.0*(yabzbc)*terma)/(termc2) )
998     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(rbc2*termc));
999     dedxicxid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1000    
1001     term5 = (4.0*xbc*ybc*dot)/(rrab2rbc22) - (4.0*xbc*yab*dot2)/(rab2rbc2*rab2rbc2) - (2.0*xab*ybc)/(rab2rbc2)
1002     + (4.0*xab*yab*dot)/(rrab22rbc2);
1003     term1 = ( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2))
1004     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1005     term2 = -(zbd)/termc
1006     -(term1xc*(xbdzbc))/termss2
1007     -(xbc*(xbdzbc))/(rbc2*termc)
1008     -(term1ya*(ybdzab))/termss2
1009     -(yab*(ybdzab))/(termcrab2)
1010     +(3.0*term1xc*term1ya*terma)/(terms4)
1011     -(term5*terma)/termss2
1012     +(yab*term1xc*terma)/(2.0*rab2*termss)
1013     +(xbc*term1ya*terma)/(2.0*rbc2*termss)
1014     +(xbc*yab*terma)/(rab2*rbc2*termc);
1015     term4 = (-(2.0*(xbdzbc)*terma)/(termc2) + (term1ya*terma2)/(termc*termss)
1016     +(2.0*yab*terma2)/(termc2*rab2))
1017     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(rbc2*termc));
1018     dedxicyia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1019    
1020     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1021     term3 = (4.0*yab*xbc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*xab*ybc*dot)/(rrab2rbc22) - (2.0*yiabc*xab)/rab2rbc2
1022     +(4.0*xbc*yiabc*dot)/(rrab22rbc2) - (4.0*xab*yab*dot)/(rrab22rbc2)
1023     +(8.0*xbc*ybc*dot2)/(rrab2rbc22*rbc2);
1024     term1 = ( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2))
1025     *( (ybd*terma)/(termc*rbd2) + (xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))/termc
1026     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1027     term2 = (ybd*(ybdzab))/(termc*rbd2)
1028     -(ybd*term1xc*terma)/(2.0*termss*rbd2)
1029     -(xbc*ybd*terma)/(rbc2*rbd2*termc)
1030     +(zbd-zab)/termc
1031     -(term1xc*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/termss2
1032     -(xbc*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/(termcrbc2)
1033     -(term5*(ybdzab))/termss2
1034     +(ybc*(ybdzab))/(termcrbc2)
1035     +(yab*(ybdzab))/(termcrab2)
1036     +(3.0*term5*term1xc*terma)/(terms4)
1037     -(term3*terma)/termss2
1038     -(ybc*term1xc*terma)/(2.0*termss*rbc2)
1039     -(yab*term1xc*terma)/(2.0*termss*rab2)
1040     +(xbc*term5*terma)/(termss*rbc2)
1041     -(yab*xbc*terma)/(termcrab2*rbc2)
1042     -(3.0*xbc*ybc*terma)/(termcrbc2*rbc2);
1043     term4 = (-(2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))*terma)/(termc2)
1044     +(term5*terma2)/(termc*termss) - (2.0*ybc*terma2)/(termc2*rbc2) - (2.0*yab*terma2)/(termc2*rab2))
1045     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(rbc2*termc));
1046     dedxicyib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1047    
1048     term5 = (4.0*xbc*yab*dot)/(rrab2rbc22) - (8.0*xbc*ybc*dot2)/(rrab2rbc22*rbc2) + (4.0*xab*ybc*dot)/(rrab2rbc22)
1049     - (2.0*xab*yab)/(rab2rbc2);
1050     term1 = ( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2))
1051     *( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2));
1052     term2 = -(term1xc*(xabzbd))/termss2
1053     -(xbc*(xabzbd))/(termcrbc2)
1054     -(term1yc*(ybdzab))/termss2
1055     -(ybc*(ybdzab))/(termcrbc2)
1056     +(3.0*term1xc*term1yc*terma)/(terms4)
1057     -(term5*terma)/termss2
1058     +(ybc*term1xc*terma)/(2.0*rbc2*termss)
1059     +(xbc*term1yc*terma)/(2.0*termss*rbc2)
1060     +(3.0*xbc*ybc*terma)/(termcrbc2*rbc2);
1061     term4 = (-(2.0*(xabzbd)*terma)/(termc2) + (term1yc*terma2)/(termc*termss)
1062     +(2.0*ybc*terma2)/(termc2*rbc2))
1063     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2));
1064     dedxicyic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1065    
1066     term1 = (-(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*( (ybdzab)/termc
1067     -(term1yc*terma)/termss2 - (xbc*terma)/(termcrbc2));
1068     term2 = -(ybd*(ybdzab))/(termc*rbd2)
1069     +(ybd*term1xc*terma)/(2.0*termss*rbd2)
1070     +(xbc*ybd*terma)/(termcrbc2*rbd2)
1071     -((xbczab)*term1xc)/termss2
1072     -(xbc*(xbczab))/(termcrbc2)
1073     +zab/termc;
1074     term4 = ( (2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xbczab)*terma)/(termc2))
1075     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(rbc2*termc));
1076     dedxicyid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1077    
1078     term5 = (4.0*xbc*zbc*dot)/(rrab2rbc22) + (4.0*xab*zab*dot)/(rrab22rbc2) - (2.0*xab*zbc)/rab2rbc2
1079     -(4.0*xbc*zab*dot2)/(rab2rbc2*rab2rbc2);
1080     term1 = ( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2))
1081     *( (xbcybd)/termc - (term1za*terma)/termss2 - (zab*terma)/(termcrab2));
1082     term2 = -(term1xc*(xbcybd))/termss2
1083     -(xbc*(xbcybd))/(termcrbc2)
1084     +ybd/termc
1085     -(zab*(ybdzab))/(termcrab2)
1086     -(term1za*(ybdzab))/termss2
1087     +(zab*term1xc*terma)/(2.0*rab2*termss)
1088     -(term5*terma)/termss2
1089     +(xbc*zab*terma)/(rab2*rbc2*termc)
1090     +(3.0*term1xc*term1za*terma)/(terms4)
1091     +(xbc*term1za*terma)/(2.0*rbc2*termss);
1092     term4 = (-(2.0*(xbcybd)*terma)/(termc2) + (2.0*zab*terma2)/(termc2*rab2)
1093     +(term1za*terma2)/(termc*termss))
1094     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(rbc2*termc));
1095     dedxiczia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1096    
1097     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1098     term3 = (4.0*zab*xbc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*xab*zbc*dot)/(rrab2rbc22) - (2.0*ziabc*xab)/rab2rbc2
1099     +(4.0*xbc*ziabc*dot)/(rrab22rbc2) - (4.0*zab*xab*dot)/(rrab22rbc2)
1100     +(8.0*xbc*zbc*dot2)/(rrab22rbc2*rab2);
1101     term1 = ( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2))
1102     *( (zbd*terma)/(termc*rbd2) + (xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
1103     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term5*terma)/termss2);
1104     term2 = (zbd*(ybdzab))/(termc*rbd2)
1105     -(zbd*term1xc*terma)/(2.0*termss*rbd2)
1106     -(xbc*zbd*terma)/(rbc2*rbd2*termc)
1107     -(term1xc*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/termss2
1108     -(xbc*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrbc2)
1109     +(yab-ybd)/termc
1110     +(zbc*(ybdzab))/(termcrbc2)
1111     +(zab*(ybdzab))/(termcrab2)
1112     -(term5*(ybdzab))/termss2
1113     -(zbc*term1xc*terma)/(2.0*termss*rbc2)
1114     -(zab*term1xc*terma)/(2.0*termss*rab2)
1115     -(term3*terma)/termss2
1116     -(3.0*xbc*zbc*terma)/(termcrbc2*rbc2)
1117     -(zab*xbc*terma)/(termcrab2*rbc2)
1118     +(3.0*term5*term1xc*terma)/(terms4)
1119     +(xbc*term5*terma)/(termss*rbc2);
1120     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))*terma)/(termc2)
1121     +(term5*terma2)/(termc*termss) - (2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2))
1122     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(rbc2*termc));
1123     dedxiczib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1124    
1125     term5 = (4.0*xbc*zab*dot)/(rrab2rbc22) + (4.0*xab*zbc*dot)/(rrab2rbc22) - (2.0*xab*zab)/rab2rbc2
1126     - (8.0*xbc*zbc*dot2)/(rbc2*rrab2rbc22);
1127     term1 = ( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2))
1128     *( (yabxbd)/termc - (term1zca)/termss2 - (zbc*terma)/(termcrbc2));
1129     term2 = -(term1xc*(yabxbd))/termss2
1130     -(xbc*(yabxbd))/(rbc2*termc)
1131     -(term1zc*(ybdzab))/termss2
1132     -(zbc*(ybdzab))/(termcrbc2)
1133     +(3.0*term1zc*term1xc*terma)/(terms4)
1134     -(term5*terma)/(2.0*termss*rbc2)
1135     +(xbc*term1zca)/(2.0*rbc2*termss)
1136     +(zbc*term1xc*terma)/(2.0*rbc2*termss)
1137     +(3.0*xbc*zbc*terma)/(rbc2*rbc2*termc);
1138     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
1139     +(2.0*zbc*terma2)/(rbc2*termc2))
1140     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(rbc2*termc));
1141     dedxiczic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1142    
1143     term1 = (-(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*( (ybdzab)/termc
1144     -(term1xc*terma)/termss2 - (xbc*terma)/(termcrbc2));
1145     term2 = -(zbd*(ybdzab))/(rbd2*termc)
1146     +(term1xc*zbd*terma)/(2.0*termss*rbd2)
1147     +(xbc*zbd*terma)/(rbc2*rbd2*termc)
1148     -(term1xc*(xabybc))/termss2
1149     -(xbc*(xabybc))/(rbc2*termc)
1150     -yab/termc;
1151     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1152     *( (ybdzab)/termc - (term1xc*terma)/termss2 - (xbc*terma)/(rbc2*termc));
1153     dedxiczid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1154    
1155     term1 = -(xbd*terma)/(termc*rbd2) + (yabzbc)/termc;
1156     term2 = (3.0*xbd2*terma)/(termcrbd22) - (2.0*xbd*(yabzbc))/(termc*rbd2)
1157     -(terma)/(termc*rbd2);
1158     term4 = ((2.0*xbd*terma2)/(termc2*rbd2) - (2.0*(yabzbc)*terma)/(termc2))
1159     *(-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc );
1160     dedxidxid = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1161    
1162    
1163     term1 = (-(xbd*terma)/(termc*rbd2)+(yabzbc)/termc)*( (xbdzbc)/termc
1164     -(term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1165     term2 = -(xbd*(xbdzbc))/(termc*rbd2)
1166     +(xbd*term1ya*terma)/(2.0*termss*rbd2)
1167     +(xbd*yab*terma)/(termcrab2*rbd2)
1168     -(term1ya*(yabzbc))/termss2
1169     +(zbc)/termc
1170     -(yab*(yabzbc))/(termcrab2);
1171     term4 = (-(2.0*(xbdzbc)*terma)/(termc2) + (term1ya*terma2)/(termc*termss)
1172     +(2.0*yab*terma2)/(termc2*rab2))
1173     *( -(xbd*terma)/(termc*rbd2) + (yabzbc)/termc );
1174     dedxidyia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1175    
1176     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1177     term1 = (-(xbd*terma)/(termc*rbd2)+(yabzbc)/termc)*( (ybd*terma)/(termc*rbd2)
1178     +(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))/termc - (term5*terma)/termss2
1179     +(ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1180     term2 = -(3.0*xbd*ybd*terma)/(termcrbd22)
1181     +(ybd*(yabzbc))/(termc*rbd2)
1182     -(xbd*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/(termc*rbd2)
1183     +(xbd*term5*terma)/(2.0*termss*rbd2)
1184     -(xbd*ybc*terma)/(termcrbc2*rbd2)
1185     -(xbd*yab*terma)/(termcrab2*rbd2)
1186     -(term5*(yabzbc))/termss2
1187     +(ybc*(yabzbc))/(termcrbc2)
1188     +(zab-zbc)/termc
1189     +(yab*(yabzbc))/(termcrab2);
1190     term4 = (-(2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))*terma)/(termc2)
1191     +(term5*terma2)/(termc*termss) - (2.0*ybc*terma2)/(termc2*rbc2) - (2.0*yab*terma2)/(termc2*rab2))
1192     *(-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc );
1193     dedxidyib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1194    
1195     term1 = (-(xbd*terma)/(termc*rbd2)+(yabzbc)/termc)*( (xabzbd)/termc
1196     -(term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2));
1197     term2 = -(xbd*(xabzbd))/(termc*rbd2)
1198     +(xbd*term1yc*terma)/(2.0*termss*rbd2)
1199     +(xbd*ybc*terma)/(termcrbc2*rbd2)
1200     -(term1yc*(yabzbc))/termss2
1201     -(ybc*(yabzbc))/(rbc2*termc)
1202     -zab/termc;
1203     term4 = (-(2.0*(xabzbd)*terma)/(termc2) + (term1yc*terma2)/(termc*termss)
1204     +(2.0*ybc*terma2)/(rbc2*termc2))
1205     *(-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc );
1206     dedxidyic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1207    
1208     term1 = (-(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*( -(xbd*terma)/(termc*rbd2)
1209     +(yabzbc)/termc);
1210     term2 = (3.0*xbd*ybd*terma)/(termcrbd22)
1211     -(xbd*(xbczab))/(termc*rbd2)
1212     -(ybd*(yabzbc))/(termc*rbd2);
1213     term4 = ( (2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xbczab)*terma)/(termc2))
1214     *(-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc);
1215     dedxidyid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1216    
1217     term1 = (-(xbd*terma)/(termc*rbd2)+(yabzbc)/termc)*( (xbcybd)/termc
1218     -(zab*terma)/(termcrab2) - (term1za*terma)/termss2);
1219     term2 = -(xbd*(xbcybd))/(termc*rbd2)
1220     +(xbd*zab*terma)/(termcrab2*rbd2)
1221     +(xbd*term1za*terma)/(2.0*termss*rbd2)
1222     -(ybc)/termc
1223     -(zab*(yabzbc))/(termcrab2)
1224     -(term1za*(yabzbc))/termss2;
1225     term4 = (-(2.0*(xbd*ybc-ybd*xbc)*terma)/(termc2) + (term1za*terma2)/(termc*termss)
1226     +(2.0*zab*terma2)/(termc2*rab2))
1227     *(-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc );
1228     dedxidzia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1229    
1230     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1231     term1 = (-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc)*( (zbd*terma)/(termc*rbd2)
1232     +(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc + (zbc*terma)/(termcrbc2)
1233     +(zab*terma)/(termcrab2) - (term5*terma)/termss2);
1234     term2 = -(3.0*xbd*zbd*terma)/(termcrbd22)
1235     -(xbd*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termc*rbd2)
1236     +(zbd*(yabzbc))/(termc*rbd2)
1237     -(xbd*zbc*terma)/(termcrbc2*rbd2)
1238     -(xbd*zab*terma)/(termcrab2*rbd2)
1239     +(xbd*term5*terma)/(2.0*termss*rbd2)
1240     +(zbc*(yabzbc))/(termcrbc2)
1241     +(ybc-yab)/termc
1242     +(zab*(yabzbc))/(termcrab2)
1243     -(term5*(yabzbc))/termss2;
1244     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))*terma)/(termc2)
1245     +(term5*terma2)/(termc*termss) - (2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2))
1246     *(-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc );
1247     dedxidzib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1248    
1249     term1 = (-(xbd*terma)/(termc*rbd2)+(yabzbc)/termc)*( (yabxbd)/termc
1250     -(term1zca)/termss2 - (zbc*terma)/(termcrbc2));
1251     term2 = -(xbd*(yabxbd))/(termc*rbd2)
1252     +(xbd*term1zca)/(2.0*termss*rbd2)
1253     +(xbd*zbc*terma)/(termcrbc2*rbd2)
1254     -(term1zc*(yabzbc))/termss2
1255     -(zbc*(yabzbc))/(rbc2*termc)
1256     +yab/termc;
1257     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
1258     +(2.0*zbc*terma2)/(rbc2*termc2))
1259     *(-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc );
1260     dedxidzic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1261    
1262     term1 = (-(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*(-(xbd*terma)/(termc*rbd2)
1263     +(yabzbc)/termc);
1264     term2 = (3.0*xbd*zbd*terma)/(termcrbd22)
1265     -(xbd*(ybc*xab-yab*xbc))/(termc*rbd2)
1266     -(zbd*(yabzbc))/(termc*rbd2);
1267     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(ybc*xab-yab*xbc)*terma)/(termc2))
1268     *(-(xbd*terma)/(termc*rbd2) + (yabzbc)/termc );
1269     dedxidzid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1270    
1271     term3 = (8.0*yab*ybc*dot)/(rrab22rbc2)- (2.0*ybc2)/rab2rbc2 - (8.0*yab2*dot2)/(rab2*rrab22rbc2)
1272     + (2.0*dot2)/(rrab22rbc2);
1273     term1 = term1aa(yab, xbdzbc, terma, term1ya, termc, termss2, termcrab2);
1274     term2 = -(term1ya*(xbdzbc))/(termss)
1275     -(2.0*yab*(xbdzbc))/(termcrab2)
1276     +(3.0*terma*term1ya*term1ya)/(terms4)
1277     -(terma*term3)/termss2
1278     +(yab*term1ya*terma)/(termss*rab2)
1279     +(3.0*yab2*terma)/(rab2*termcrab2)
1280     -terma/(termcrab2);
1281     term4 = (-(2.0*terma*(xbdzbc))/(termc2)
1282     +(term1ya*terma2)/(termc*termss)
1283     +(2.0*yab*terma2)/(termc2*rab2) )
1284     *( (xbdzbc)/termc - (term1ya*terma)/termss2
1285     -(yab*terma)/(termcrab2));
1286     dedyiayia = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1287    
1288     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1289     term3 = (4.0*yab*ybc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*ybc2*dot2)/(rrab2rbc22) - (2.0*yiabc*ybc)/rab2rbc2
1290     +(4.0*yab*yiabc*dot)/(rrab22rbc2) - (4.0*yab*ybc*dot)/(rrab22rbc2) +(2.0*dot)/rab2rbc2
1291     +(8.0*yab2*dot2)/(rrab22rbc2*rab2) - (2.0*dot2)/(rrab2rbc22);
1292     term1 = ( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2))
1293     *( (ybd*terma)/(termc*rbd2) + (xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))/termc
1294     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1295     term2 = (ybd*(xbdzbc))/(termc*rbd2)
1296     -(ybd*term1ya*terma)/(2.0*termss*rbd2)
1297     -(yab*ybd*terma)/(rab2*rbd2*termc)
1298     -(term5*(xbdzbc))/termss2
1299     +(ybc*(xbdzbc))/(termcrbc2)
1300     +(yab*(xbdzbc))/(termcrab2)
1301     -(term1ya*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/termss2
1302     -(yab*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia)))/(termcrab2)
1303     +(3.0*term5*term1ya*terma)/(terms4)
1304     -(term3*terma)/termss2
1305     +(yab*term5*terma)/(2.0*termss*rab2)
1306     -(ybc*term1ya*terma)/(2.0*termss*rbc2)
1307     -(yab*term1ya*terma)/(2.0*termss*rab2)
1308     -(yab*ybc*terma)/(termcrab2*rbc2)
1309     -(3.0*yab2*terma)/(termcrab2*rab2)
1310     +terma/(termcrab2);
1311     term4 = (-(2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xab*zbc-xbc*zab+xbd*(zia-zic)+zbd*(xic-xia))*terma)/(termc2)
1312     +(term5*terma2)/(termc*termss) - (2.0*ybc*terma2)/(termc2*rbc2)
1313     -(2.0*yab*terma2)/(termc2*rab2))
1314     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1315     dedyiayib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1316    
1317     term5 = (4.0*ybc2*dot)/(rrab2rbc22) + (4.0*yab2*dot)/(rrab22rbc2) - (4.0*yab*ybc*dot2)/(rab2rbc2*rab2rbc2)
1318     - (2.0*yab*ybc)/rab2rbc2;
1319     term1 = ( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2))
1320     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1321     term2 = -(term1ya*(xabzbd))/termss2
1322     -(yab*(xabzbd))/(termcrab2)
1323     -(term1yc*(xbdzbc))/termss2
1324     -(ybc*(xbdzbc))/(rbc2*termc)
1325     +(3.0*term1yc*term1ya*terma)/(terms4)
1326     -(term5*terma)/termss2
1327     +(yab*term1yc*terma)/(2.0*termss*rab2)
1328     +(ybc*term1ya*terma)/(2.0*termss*rbc2)
1329     +(yab*ybc*terma)/(rab2*rbc2*termc);
1330     term4 = (-(2.0*(xabzbd)*terma)/(termc2) + (term1yc*terma2)/(termc*termss)
1331     +(2.0*ybc*terma2)/(rbc2*termc2))
1332     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1333     dedyiayic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1334    
1335     term1 = (-(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*( (xbdzbc)/termc
1336     -(term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1337     term2 = -(ybd*(xbdzbc))/(termc*rbd2)
1338     +(ybd*term1ya*terma)/(2.0*termss*rbd2)
1339     +(yab*ybd*terma)/(rab2*rbd2*termc)
1340     -(term1ya*(xbczab))/termss2
1341     -(yab*(xbczab))/(termcrab2);
1342     term4 = ( (2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xbczab)*terma)/(termc2))
1343     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1344     dedyiayid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1345    
1346     term5 = (4.0*ybc*zab*dot)/(rrab22rbc2) + (4.0*yab*zbc*dot)/(rrab22rbc2) - (8.0*yab*zab*dot2)/(rrab22rbc2*rab2)
1347     - (2.0*zbc*ybc)/rab2rbc2;
1348     term1 = ( (xbdzbc)/termc - (term1ya*terma)/termss2 -(yab*terma)/(termcrab2))
1349     *( (xbcybd)/termc - (zab*terma)/(termcrab2) - (term1za*terma)/termss2);
1350     term2 = -(term1ya*(ybd*xbc-xbd*ybc))/termss2
1351     -(yab*(ybd*xbc-xbd*ybc))/(termcrab2)
1352     -(zab*(xbdzbc))/(termcrab2)
1353     -(term1za*(xbdzbc))/termss2
1354     +(zab*term1ya*terma)/(2.0*termss*rab2)
1355     +(3.0*yab*zab*terma)/(termcrab2*rab2)
1356     -(term5*terma)/termss2
1357     +(3.0*term1ya*term1za*terma)/(terms4)
1358     +(yab*term1za*terma)/(2.0*termss*rab2);
1359     term4 = (-(2.0*(ybd*xbc-xbd*ybc)*terma)/(termc2) + (2.0*zab*terma2)/(termc2*rab2)
1360     +(term1za*terma2)/(termc*termss))
1361     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1362     dedyiazia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1363    
1364     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1365     term3 = (4.0*yab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*ybc*zbc*dot)/(rrab2rbc22) - (2.0*ziabc*ybc)/rab2rbc2
1366     +(4.0*yab*ziabc*dot)/(rrab22rbc2) - (4.0*zab*ybc*dot)/(rrab22rbc2)
1367     +(8.0*yab*zab*dot2)/(rrab22rbc2*rab2);
1368     term1 = ( (xbdzbc)/termc - (term1ya*terma)/termss2 -(yab*terma)/(termcrab2))
1369     *( (zbd*terma)/(termc*rbd2) + (xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
1370     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term5*terma)/termss2);
1371     term2 = (zbd*(xbdzbc))/(termc*rbd2)
1372     -(term1ya*zbd*terma)/(2.0*termss*rbd2)
1373     -(yab*zbd*terma)/(termcrab2*rbd2)
1374     -(term1ya*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/termss2
1375     -(yab*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrab2)
1376     +(xbc-xbd)/termc
1377     +(zbc*(xbdzbc))/(termcrbc2)
1378     +(zab*(xbdzbc))/(termcrab2)
1379     -(term5*(xbdzbc))/termss2
1380     -(zbc*term1ya*terma)/(2.0*termss*rbc2)
1381     -(zab*term1ya*terma)/(2.0*termss*rab2)
1382     -(yab*zbc*terma)/(termcrab2*rbc2)
1383     -(3.0*yab*zab*terma)/(termcrab2*rab2)
1384     -(term3*terma)/termss2
1385     +(3.0*term1ya*term5*terma)/(terms4)
1386     +(yab*term5*terma)/(2.0*termss*rab2);
1387     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))*terma)/(termc2)
1388     -(2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2)
1389     +(term5*terma2)/(termc*termss))
1390     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1391     dedyiazib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1392    
1393     term5 = (4.0*ybc*zbc*dot)/(rrab2rbc22) + (4.0*yab*zab*dot)/(rrab22rbc2) - (4.0*yab*zbc*dot2)/(rab2rbc2*rab2rbc2)
1394     - (2.0*zab*ybc)/rab2rbc2;
1395     term1 = ( (yabxbd)/termc - (term1zca)/termss2 - (zbc*terma)/(termcrbc2))
1396     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1397     term2 = -(term1ya*(yabxbd))/termss2
1398     -(yab*(yabxbd))/(termcrab2)
1399     + xbd/termc
1400     -(term1zc*(xbdzbc))/termss2
1401     -(zbc*(xbdzbc))/(rbc2*termc)
1402     +(3.0*term1zc*term1ya*terma)/(terms4)
1403     -(term5*terma)/termss2
1404     +(yab*term1zca)/(2.0*termss*rab2)
1405     +(zbc*term1ya*terma)/(2.0*termss*rbc2)
1406     +(yab*zbc*terma)/(rab2*rbc2*termc);
1407     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
1408     +(2.0*zbc*terma2)/(rbc2*termc2))
1409     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1410     dedyiazic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1411    
1412     term1 = (-(zbd*terma)/(termc*rbd2) + (xab*ybc-xbc*yab)/termc)*( (xbdzbc)/termc
1413     -(term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1414     term2 = -(zbd*(xbdzbc))/(termc*rbd2)
1415     +(zbd*term1ya*terma)/(2.0*termss*rbd2)
1416     +(yab*zbd*terma)/(rab2*rbd2*termc)
1417     -(term1ya*(xabybc))/termss2
1418     -(yab*(xabybc))/(termcrab2)
1419     - xbc/termc;
1420     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1421     *( (xbdzbc)/termc - (term1ya*terma)/termss2 - (yab*terma)/(termcrab2));
1422     dedyiazid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1423    
1424     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1425     term3 = -(8.0*ybc2*dot2)/(rrab2rbc22*rbc2)
1426     -(8.0*yiabc*ybc*dot)/(rrab2rbc22) + (2.0*dot2)/(rrab2rbc22)
1427     -(8.0*yab*ybc*dot2)/(rab2rbc2*rab2rbc2) - (8.0*yab*yiabc)/(rab2rbc2)
1428     -(4.0*dot)/rab2rbc2 - (8.0*yab2*dot2)/(rrab22rbc2) + (2.0*dot2)/(rrab22rbc2)
1429     -(2.0*yiabc*yiabc)/rab2rbc2;
1430     term1 = (ybd*terma)/(termc*rbd2) + xabzbc/termc
1431     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2);
1432     term2 = (3.0*ybd2*terma)/(termcrbd22)
1433     +(2.0*ybd*xabzbc)/(termc*rbd2)
1434     -(term5*ybd*terma)/(termss*rbd2)
1435     +(2.0*ybc*ybd*terma)/(termc*rbd2*rbc2)
1436     +(2.0*yab*ybd*terma)/(termc*rbd2*rab2)
1437     -terma/(termc*rbd2)
1438     -(term5*xabzbc)/(termss)
1439     +(2.0*ybc*xabzbc)/(termcrbc2)
1440     +(2.0*yab*xabzbc)/(termcrab2)
1441     +(3.0*term5*term5*terma)/(terms4)
1442     -(term3*terma)/termss2
1443     -(ybc*term5*terma)/(termss*rbc2)
1444     -(yab*term5*terma)/(termss*rab2)
1445     +(3.0*ybc2*terma)/(termcrbc2*rbc2)
1446     +(2.0*yab*ybc*terma)/(termcrab2*rbc2)
1447     +(3.0*yab2*terma)/(termcrab2*rab2)
1448     -terma/(termcrab2) - terma/(termcrbc2);
1449     term4 = ( -(2.0*ybd*terma2)/(termc2*rbd2)
1450     -(2.0*xabzbc*terma)/(termc2)
1451     +(term5*terma2)/(termc*termss) - (2.0*ybc*terma2)/(termc2*rbc2)
1452     -(2.0*yab*terma2)/(termc2*rab2))
1453     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1454     +(ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2) - (term5*terma)/termss2);
1455     dedyibyib = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1456    
1457     term3 = (8.0*ybc2*dot2)/(rrab2rbc22*rbc2) - (4.0*yab*ybc*dot)/(rrab2rbc22)
1458     +(4.0*ybc*dot*yiabc)/(rrab2rbc22) + (4.0*yab*ybc*dot2)/(rab2rbc2*rab2rbc2)
1459     -(2.0*dot2)/(rrab2rbc22) - (2.0*yab*yiabc)/rab2rbc2 - (4.0*yab2*dot)/(rrab22rbc2)
1460     +(2.0*dot)/rab2rbc2;
1461     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1462     term1 = ( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2))
1463     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1464     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1465     term2 = (ybd*(xabzbd))/(termc*rbd2)
1466     -(ybd*term1yc*terma)/(2.0*termss*rbd2)
1467     -(ybc*ybd*terma)/(termcrbc2*rbd2)
1468     -(term5*(xabzbd))/termss2
1469     +(ybc*(xabzbd))/(termcrbc2)
1470     +(yab*(xabzbd))/(termcrab2)
1471     -(term1yc*xabzbc)/termss2
1472     -(ybc*xabzbc)/(termcrbc2)
1473     +(3.0*term1yc*term5*terma)/(terms4)
1474     -(term3*terma)/termss2
1475     -(ybc*term1yc*terma)/(2.0*termss*rbc2)
1476     -(yab*term1yc*terma)/(2.0*termss*rab2)
1477     +(ybc*term5*terma)/(2.0*termss*rbc2)
1478     -(3.0*ybc2*terma)/(termcrbc2*rbc2)
1479     -(yab*ybc*terma)/(termcrbc2*rab2)
1480     +(terma)/(termcrbc2);
1481     term4 = (-(2.0*(xabzbd)*terma)/(termc2) + (term1yc*terma2)/(termc*termss)
1482     +(2.0*ybc*terma2)/(termc2*rbc2))
1483     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1484     +(ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2) - (term5*terma)/termss2);
1485     dedyibyic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1486    
1487     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1488     term1 = (-(ybd*terma)/(termc*rbd2) + (xbczab)/termc )*( (ybd*terma)/(termc*rbd2)
1489     +xabzbc/termc
1490     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1491     term2 = -(3.0*ybd2*terma)/(termcrbd22)
1492     +(ybd*(xbczab))/(termc*rbd2)
1493     -(ybd*xabzbc)/(termc*rbd2)
1494     +(ybd*term5*terma)/(2.0*termss*rbd2)
1495     -(ybc*ybd*terma)/(termcrbc2*rbd2)
1496     -(yab*ybd*terma)/(termcrab2*rbd2)
1497     +terma/(termc*rbd2)
1498     -(term5*(xbczab))/termss2
1499     +(ybc*(xbczab))/(termcrbc2)
1500     +(yab*(xbczab))/(termcrab2);
1501     term4 = ( (2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xbczab))/(termc2))
1502     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1503     +(ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2) - (term5*terma)/termss2);
1504     dedyibyid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1505    
1506     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1507     term3 = (4.0*zab*ybc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*zbc*ybc*dot)/(rrab2rbc22) - (2.0*yiabc*zbc)/rab2rbc2
1508     +(4.0*zab*yiabc*dot)/(rrab22rbc2) - (4.0*yab*zbc*dot)/(rrab22rbc2)
1509     +(8.0*zab*yab*dot2)/(rrab22rbc2*rab2);
1510     term1 = ( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1511     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2))
1512     *( (xbcybd)/termc - (zab*terma)/(termcrab2) - (term1za*terma)/termss2);
1513     term2 = (ybd*(xbcybd))/(termc*rbd2)
1514     -(ybd*zab*terma)/(termcrab2*rbd2)
1515     -(ybd*term1za*terma)/(2.0*termss*rbd2)
1516     -(term5*(xbcybd))/termss2
1517     +(ybc*(xbcybd))/(termcrbc2)
1518     +(yab*(xbcybd))/(termcrab2)
1519     +(xbd-xbc)/termc
1520     -(zab*xabzbc)/(termcrab2)
1521     -(term1za*xabzbc)/termss2
1522     +(zab*term5*terma)/(2.0*termss*rab2)
1523     -(ybc*zab*terma)/(termcrab2*rbc2)
1524     -(3.0*yab*zab*terma)/(termcrab2*rab2)
1525     -(term3*terma)/termss2
1526     +(3.0*term5*term1za*terma)/(terms4)
1527     -(ybc*term1za*terma)/(2.0*termss*rbc2)
1528     -(yab*term1za*terma)/(2.0*termss*rab2);
1529     term4 = (-(2.0*(xbcybd)*terma)/(termc2) + (2.0*zab*terma2)/(termc2*rab2)
1530     +(term1za*terma2)/(termc*termss))
1531     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1532     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1533     dedyibzia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1534    
1535     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1536     term6 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1537     term3 = -(8.0*ybc*zbc*dot2)/(rrab2rbc22*rbc2) - (4.0*ybc*ziabc*dot)/(rrab2rbc22)
1538     -(4.0*zbc*dot*yiabc)/(rrab2rbc22) -(4.0*ybc*zab*dot2)/(rab2rbc2*rab2rbc2)
1539     -(4.0*yab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (2.0*yiabc*ziabc)/rab2rbc2
1540     -(4.0*zab*dot*yiabc)/(rrab22rbc2) - (4.0*yab*dot*ziabc)/(rrab22rbc2)
1541     -(8.0*yab*zab*dot2)/(rrab22rbc2*rab2);
1542     term1 = ( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1543     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2))
1544     *( (zbd*terma)/(termc*rbd2) + (yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
1545     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term6*terma)/termss2);
1546     term2 = (3.0*ybd*zbd*terma)/(termcrbd22)
1547     +(ybd*(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termc*rbd2)
1548     +(zbd*xabzbc)/(termc*rbd2)
1549     +(ybd*zbc*terma)/(termc*rbd2*rbc2)
1550     +(ybd*zab*terma)/(termc*rbd2*rab2)
1551     -(ybd*term6*terma)/(2.0*termss*rbd2)
1552     -(zbd*term5*terma)/(2.0*termss*rbd2)
1553     +(ybc*zbd*terma)/(termcrbc2*rbd2)
1554     +(yab*zbd*terma)/(termcrab2*rbd2)
1555     -(term5*(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/termss2
1556     +(ybc*(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrbc2)
1557     +(yab*(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrab2)
1558     +(zbc*xabzbc)/(termcrbc2)
1559     +(zab*xabzbc)/(termcrab2)
1560     -(term6*xabzbc)/termss2
1561     -(zbc*term5*terma)/(2.0*termss*rbc2)
1562     -(zab*term5*terma)/(2.0*termss*rab2)
1563     +(3.0*ybc*zbc*terma)/(termcrbc2*rbc2)
1564     +(ybc*zab*terma)/(termcrbc2*rab2)
1565     +(yab*zbc*terma)/(termcrbc2*rab2)
1566     +(3.0*yab*zab*terma)/(termcrab2*rab2)
1567     -(term3*terma)/termss2
1568     +(3.0*term5*term6*terma)/(terms4)
1569     -(ybc*term6*terma)/(2.0*termss*rbc2)
1570     -(yab*term6*terma)/(2.0*termss*rab2);
1571     term4 = (-(2.0*zbd*terma2)/(termc2) - (2.0*(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))*terma)/(termc2)
1572     -(2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2)
1573     +(term6*terma2)/(termc*termss))
1574     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1575     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1576     dedyibzib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1577    
1578     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1579     term3 = (4.0*yab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*yab*zab*dot)/(rrab2rbc22) - (2.0*yiabc*zab)/rab2rbc2
1580     +(4.0*zbc*yiabc*dot)/(rrab2rbc22) - (4.0*zab*ybc*dot)/(rrab2rbc22)
1581     +(8.0*zbc*ybc*dot2)/(rrab2rbc22*rbc2);
1582     term1 = ( (xbd*yab-ybd*xab)/termc - (term1zca)/termss2 - (zbc*terma)/(termcrbc2))
1583     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1584     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1585     term2 = (ybd*(yabxbd))/(termc*rbd2)
1586     -(ybd*term1zca)/(2.0*termss*rbd2)
1587     -(ybd*zbc*terma)/(termcrbc2*rbd2)
1588     -(term5*(yabxbd))/termss2
1589     +(ybc*(yabxbd))/(termcrbc2)
1590     +(yab*(yabxbd))/(termcrab2)
1591     +(xab-xbd)/termc
1592     -(term1zc*xabzbc)/termss2
1593     -(zbc*xabzbc)/(termcrbc2)
1594     +(3.0*term5*term1zca)/(terms4)
1595     -(term3*terma)/termss2
1596     -(ybc*term1zca)/(2.0*termss*rbc2)
1597     -(yab*term1zca)/(2.0*termss*rab2)
1598     +(zbc*term5*terma)/(2.0*termss*rbc2)
1599     -(3.0*ybc*zbc*terma)/(termcrbc2*rbc2)
1600     -(zbc*yab*terma)/(termcrab2*rbc2);
1601     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (2.0*zbc*terma2)/(termc2*rbc2)
1602     +(term1zca*terma)/(termc*termss))
1603     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1604     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1605     dedyibzic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1606    
1607     term5 = -(2.0*ybc*dot2)/(rrab2rbc22) - (2.0*yiabc*dot)/rab2rbc2 - (2.0*yab*dot2)/(rrab22rbc2);
1608     term1 = (-(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*((ybd*terma)/(termc*rbd2)
1609     +xabzbc/termc - (term5*terma)/termss2
1610     +(ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1611     term2 = -(3.0*ybd*zbd*terma)/(termcrbd22)
1612     +(ybd*(xabybc))/(termc*rbd2)
1613     -(zbd*xabzbc)/(termc*rbd2)
1614     +(term5*zbd*terma)/(2.0*termss*rbd2)
1615     -(ybc*zbd*terma)/(termcrbc2*rbd2)
1616     -(yab*zbd*terma)/(termcrab2*rbd2)
1617     -(term5*(xabybc))/termss2
1618     +(ybc*(xabybc))/(termcrbc2)
1619     +(yab*(xabybc))/(termcrab2)
1620     +(xbc-xab)/termc;
1621     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1622     *( (ybd*terma)/(termc*rbd2) + xabzbc/termc
1623     -(term5*terma)/termss2 + (ybc*terma)/(termcrbc2) + (yab*terma)/(termcrab2));
1624     dedyibzid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1625    
1626     term3 = (8.0*yab*ybc*dot)/(rrab2rbc22)- (2.0*yab2)/rab2rbc2 - (8.0*ybc2*dot2)/(rbc2*rrab2rbc22)
1627     + (2.0*dot2)/(rrab2rbc22);
1628     term1 = term1aa(ybc, xabzbd, terma, term1yc, termc, termss2, termcrbc2);
1629     term2 = -(term1yc*(xabzbd))/(termss)
1630     -(2.0*ybc*(xabzbd))/(termcrbc2)
1631     +(3.0*terma*term1yc*term1yc)/(terms4)
1632     -(terma*term3)/termss2
1633     +(ybc*term1yc*terma)/(termss*rbc2)
1634     +(3.0*ybc2*terma)/(rbc2*rbc2*termc) - terma/(rbc2*termc);
1635     term4 = (-(2.0*terma*(xabzbd))/(termc2)
1636     +(term1yc*terma2)/(termc*termss)
1637     +(2.0*ybc*terma2)/(rbc2*termc2) )
1638     *( (xabzbd)/termc - (term1yc*terma)/termss2
1639     -(ybc*terma)/(rbc2*termc));
1640     dedyicyic = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1641    
1642     term1 = (-(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*( (xabzbd)/termc
1643     -(term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2));
1644     term2 = -(ybd*(xabzbd))/(termc*rbd2)
1645     +(ybd*term1yc*terma)/(2.0*termss*rbd2)
1646     +(ybc*ybd*terma)/(rbc2*rbd2*termc)
1647     -(term1yc*(xbczab))/termss2
1648     -(ybc*(xbczab))/(rbc2*termc);
1649     term4 = ( (2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xbczab)*terma)/(termc2))
1650     *( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrab2));
1651     dedyicyid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1652    
1653     term5 = (4.0*ybc*zbc*dot)/(rrab2rbc22) + (4.0*yab*zab*dot)/(rrab22rbc2) - (4.0*ybc*zab*dot2)/(rab2rbc2*rab2rbc2)
1654     - (2.0*zbc*yab)/rab2rbc2;
1655     term1 = ( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2))
1656     *( (xbcybd)/termc - (zab*terma)/(termcrbc2) - (term1za*terma)/termss2);
1657     term2 = -(term1yc*(ybd*xbc-xbd*ybc))/termss2
1658     -(ybc*(ybd*xbc-xbd*ybc))/(rbc2*termc)
1659     - xbd/(termc)
1660     -(zab*(xabzbd))/(termcrab2)
1661     -(term1za*(xabzbd))/termss2
1662     +(zab*term1yc*terma)/(2.0*termss*rab2)
1663     +(ybc*zab*terma)/(termcrab2*rab2)
1664     -(term5*terma)/termss2
1665     +(ybc*zab*terma)/(termcrab2*rbc2)
1666     +(3.0*term1yc*term1za*terma)/(terms4)
1667     +(ybc*term1za*terma)/(2.0*termss*rab2);
1668     term4 = (-(2.0*(ybd*xbc-xbd*ybc)*terma)/(termc2) + (2.0*zab*terma2)/(termc2*rab2)
1669     +(term1za*terma2)/(termc*termss))
1670     *( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(rbc2*termc));
1671     dedyiczia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1672    
1673     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1674     term3 = (8.0*ybc*zbc*dot2)/(rrab2rbc22*rbc2) + (4.0*ybc*ziabc*dot)/(rrab2rbc22)
1675     -(4.0*yab*zbc*dot)/(rrab2rbc22) + (4.0*zab*ybc*dot2)/(rab2rbc2*rab2rbc2)
1676     -(4.0*yab*zab*dot)/(rrab2rbc22) - (2.0*ziabc*yab)/rab2rbc2;
1677     term1 = ( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2))
1678     *( (zbd*terma)/(termc*rbd2) + (xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
1679     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term5*terma)/termss2);
1680     term2 = (zbd*(xabzbd))/(termc*rbd2)
1681     -(term1yc*zbd*terma)/(2.0*termss*rbd2)
1682     -(ybc*zbd*terma)/(termcrbc2*rbd2)
1683     -(term1yc*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/termss2
1684     -(ybc*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrab2)
1685     +(xbd-xab)/termc
1686     +(zbc*(xabzbd))/(termcrbc2)
1687     +(zab*(xabzbd))/(termcrab2)
1688     -(term5*(xabzbd))/termss2
1689     -(zbc*term1yc*terma)/(2.0*termss*rbc2)
1690     -(zab*term1yc*terma)/(2.0*termss*rab2)
1691     -(term3*terma)/termss2
1692     -(3.0*ybc*zbc*terma)/(termcrbc2*rbc2)
1693     -(ybc*zab*terma)/(termcrab2*rbc2)
1694     +(3.0*term1yc*term5*terma)/(terms4)
1695     +(ybc*term5*terma)/(2.0*termss*rab2);
1696     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))*terma)/(termc2)
1697     -(2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2)
1698     +(term5*terma2)/(termc*termss))
1699     *( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2));
1700     dedyiczib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1701    
1702     term5 = (4.0*ybc*zab*dot)/(rrab2rbc22) + (4.0*yab*zbc*dot)/(rrab2rbc22) - (8.0*ybc*zbc*dot2)/(rrab2rbc22*rbc2)
1703     - (2.0*zab*yab)/rab2rbc2;
1704     term1 = ( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2))
1705     *( (yabxbd)/termc - (term1zca)/termss2 - (zbc*terma)/(termcrbc2));
1706     term2 = -(term1yc*(yabxbd))/termss2
1707     -(ybc*(yabxbd))/(rbc2*termc)
1708     -(term1zc*(xabzbd))/termss2
1709     -(zbc*(xabzbd))/(rbc2*termc)
1710     +(3.0*term1zc*term1yc*terma)/(terms4)
1711     -(term5*terma)/termss2
1712     +(ybc*term1zca)/(2.0*termss*rbc2)
1713     +(zbc*term1yc*terma)/(2.0*termss*rbc2)
1714     +(3.0*ybc*zbc*terma)/(rbc2*rbc2*termc);
1715     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
1716     +(2.0*zbc*terma2)/(rbc2*termc2))
1717     *( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(rbc2*termc));
1718     dedyiczic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1719    
1720     term1 = ( -(zbd*terma)/(termc*rbd2) + (xab*ybc-xbc*yab)/termc)*( (xabzbd)/termc
1721     -(term1yc*terma)/termss2 - (ybc*terma)/(termcrbc2));
1722     term2 = -(zbd*(xabzbd))/(termc*rbd2)
1723     +(zbd*term1yc*terma)/(2.0*termss*rbd2)
1724     +(ybc*zbd*terma)/(rbc2*rbd2*termc)
1725     -(term1yc*(xabybc))/termss2
1726     -(ybc*(xabybc))/(rbc2*termc) + xab/termc;
1727     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1728     *( (xabzbd)/termc - (term1yc*terma)/termss2 - (ybc*terma)/(rbc2*termc));
1729     dedyiczid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1730    
1731     term1 = -(ybd*terma)/(termc*rbd2) + (xbczab)/termc;
1732     term2 = (3.0*ybd2*terma)/(termcrbd22) - (2.0*ybd*(xbczab))/(termc*rbd2)
1733     -(terma)/(termc*rbd2);
1734     term4 = ( (2.0*ybd*terma2)/(termc2*rbd2) - (2.0*(xbczab)*terma)/(termc2))
1735     *(-(ybd*terma)/(termc*rbd2) + (xbczab)/termc);
1736     dedyidyid = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1737    
1738     term1 = (-(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*((xbcybd)/termc
1739     -(zab*terma)/(termcrab2) - (term1za*terma)/termss2);
1740     term2 = -(ybd*(xbcybd))/(termc*rbd2) + (ybd*zab*terma)/(rab2*rbd2*termc)
1741     +(ybd*term1za*terma)/(2.0*termss*rbd2) + xbc/termc
1742     -(zab*(xbczab))/(termcrab2)
1743     -((xbczab)*term1za)/termss2;
1744     term4 = (-(2.0*(xbcybd)*terma)/(termc2) + (2.0*zab*terma2)/(termc2*rab2)
1745     +(term1za*terma2)/(termc*termss))
1746     *(-(ybd*terma)/(termc*rbd2) + (xbczab)/termc ) ;
1747     dedyidzia = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1748    
1749     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1750     term1 = ( -(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*( (zbd*terma)/(termc*rbd2)
1751     +(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc + (zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2)
1752     -(term5*terma)/termss2 );
1753     term2 = -(3.0*ybd*zbd*terma)/(termcrbd22)
1754     -(ybd*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termc*rbd2)
1755     +(zbd*(xbczab))/(termc*rbd2)
1756     -(ybd*zbc*terma)/(termcrbc2*rbd2)
1757     -(ybd*zab*terma)/(termcrab2*rbd2)
1758     +(ybd*term5*terma)/(2.0*termss*rbd2)
1759     +(zbc*(xbczab))/(termcrbc2)
1760     +(zab*(xbczab))/(termcrab2)
1761     +(xab-xbc)/termc
1762     -(term5*(xbczab))/termss2;
1763     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2)
1764     -(2.0*terma*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termc2)
1765     -(2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2)
1766     +(term5*terma2)/(termc*termss*sine))
1767     *(-(ybd*terma)/(termc*rbd2) + (xbczab)/termc );
1768     dedyidzib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1769    
1770     term1 = ( -(ybd*terma)/(termc*rbd2) + (xbczab)/termc)*((yabxbd)/termc
1771     -(term1zca)/termss2 - (zbc*terma)/(termcrbc2) );
1772     term2 = -(ybd*(yabxbd))/(termc*rbd2)
1773     +(ybd*term1zca)/(2.0*termss*rbd2)
1774     +(ybd*zbc*terma)/(rbc2*rbd2*termc)
1775     -((xbczab)*term1zc)/termss2
1776     -(zbc*(xbczab))/(termcrbc2)
1777     -xab/termc ;
1778     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
1779     +(2.0*zbc*terma2)/(termc2*rbc2))
1780     *(-(ybd*terma)/(termc*rbd2) + (xbczab)/termc );
1781     dedyidzic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1782    
1783     term1 = ( -(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*(-(ybd*terma)/(termc*rbd2)
1784     +(xbczab)/termc);
1785     term2 = (3.0*ybd*zbd*terma)/(termcrbd22) - (ybd*(xabybc))/(termc*rbd2)
1786     -(zbd*(xbczab))/(termc*rbd2);
1787     term4 = ((2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1788     *(-(ybd*terma)/(termc*rbd2) + (xbczab)/termc );
1789     dedyidzid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1790    
1791     term3 = (8.0*zab*zbc*dot)/(rrab22rbc2)- (2.0*zbc2)/rab2rbc2 - (8.0*zab2*dot2)/(rab2*rrab22rbc2)
1792     + (2.0*dot2)/(rrab22rbc2);
1793     term1 = term1aa(zab, xbcybd, terma, term1za, termc, termss2, termcrab2);
1794     term2 = -(2.0*zab*(xbcybd))/(termcrab2)
1795     -(term1za*(xbcybd))/(termss)
1796     -(terma)/(termcrab2)
1797     +(3.0*zab2*terma)/(rab2*termcrab2)
1798     +(zab*term1za*terma)/(termss*rab2)
1799     +(3.0*terma*term1za*term1za)/(terms4)
1800     -(term3*terma)/termss2;
1801     term4 = ( -(2.0*terma*(xbcybd))/(termc2)
1802     +(term1za*terma2)/(termc*termss)
1803     +(2.0*zab*terma2)/(termc2*rab2) )
1804     *( (xbcybd)/termc - (term1za*terma)/termss2
1805     -(zab*terma)/(termcrab2));
1806     dedziazia = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1807    
1808     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1809     term3 = (4.0*zab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (4.0*zbc2*dot)/(rrab2rbc22) - (2.0*ziabc*zbc)/rab2rbc2
1810     +(4.0*zab*ziabc*dot)/(rrab22rbc2) - (4.0*zab*zbc*dot)/(rrab22rbc2)
1811     +(8.0*zab2*dot2)/(rrab22rbc2*rab2) + (2.0*dot)/rab2rbc2 - (2.0*dot2)/(rrab22rbc2);
1812     term1 = ( (zbd*terma)/(termc*rbd2) + (xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
1813     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term5*terma)/termss2)
1814     *( (xbcybd)/termc - (zab*terma)/(termcrab2) - (term1za*terma)/termss2);
1815     term2 = (zbd*(xbcybd))/(termc*rbd2)
1816     -(zab*zbd*terma)/(rab2*rbd2*termc)
1817     -(zbd*term1za*terma)/(2.0*termss*rbd2)
1818     +(zbc*(xbcybd))/(termcrbc2)
1819     -(zab*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrab2)
1820     +(zab*(xbcybd))/(termcrab2)
1821     -(term5*(xbcybd))/termss2
1822     -(term1za*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/termss2
1823     -(zab*zbc*terma)/(termcrab2*rbc2)
1824     +(terma)/(termcrab2)
1825     -(3.0*zab2*terma)/(termcrab2*rab2)
1826     +(zab*term5*terma)/(2.0*termss*rab2)
1827     -(zbc*term1za*terma)/(2.0*termss*rbc2)
1828     -(zab*term1za*terma)/(2.0*termss*rab2)
1829     +(3.0*term5*term1za*terma)/(terms4)
1830     -(term3*terma)/termss2;
1831     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xbc*yab-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))*terma)/(termc2)
1832     +(term5*terma2)/(termc*termss) - (2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2))
1833     *( (xbcybd)/termc - (term1za*terma)/termss2
1834     -(zab*terma)/(termcrab2));
1835     dedziazib = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1836    
1837     term3 = (4.0*zbc2*dot)/(rab2rbc2*rbd2) - (4.0*zab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (2.0*zab*zbc)/rab2rbc2
1838     - (2.0*dot)/rab2rbc2 + (4.0*zab2*dot)/(rrab22rbc2);
1839     term1 = ( (yabxbd)/termc - (term1zca)/termss2 - (zbc*terma)/(termcrbc2))
1840     *( (xbcybd)/termc - (zab*terma)/(termcrab2) - (term1za*terma)/termss2);
1841     term2 = -((xbcybd)*term1zc)/termss2
1842     -(zbc*(xbcybd))/(termcrbc2)
1843     -(zab*(yabxbd))/(termcrab2)
1844     -(term1za*(yabxbd))/termss2
1845     +(zab*term1zca)/(2.0*termss*rab2)
1846     -(term3*terma)/termss2
1847     +(zab*zbc*terma)/(termcrab2*rbc2)
1848     +(3.0*term1zc*term1za*terma)/(terms4)
1849     +(zbc*term1za*terma)/(2.0*termss*rbc2);
1850     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
1851     +(2.0*zbc*terma2)/(termc2*rbc2))
1852     *( (xbcybd)/termc - (zab*terma)/(termcrab2) - (term1za*terma)/termss2);
1853     dedziazic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1854    
1855     term1 = (-(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*( (xbcybd)/termc
1856     -(zab*terma)/(termcrab2) - (term1za*terma)/termss2 );
1857     term2 = -(zbd*(xbcybd))/(termc*rbd2)
1858     +(zab*zbd*terma)/(termcrab2*rbd2)
1859     +(zbd*term1za*terma)/(2.0*termss*rbd2)
1860     -(zab*(xabybc))/(termcrab2)
1861     -((xabybc)*term1za)/termss2;
1862     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1863     *( (xbcybd)/termc - (zab*terma)/(termcrab2) - (term1za*terma)/termss2);
1864     dedziazid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1865    
1866     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1867     term3 = -(8.0*zbc2*dot2)/(rrab2rbc22*rbc2)
1868     -(8.0*ziabc*zbc*dot)/(rrab2rbc22) + (2.0*dot2)/(rrab2rbc22)
1869     -(8.0*zab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (8.0*zab*ziabc)/(rab2rbc2)
1870     -(4.0*dot)/rab2rbc2 - (8.0*zab2*dot2)/(rrab22rbc2) + (2.0*dot2)/(rrab22rbc2)
1871     -(2.0*ziabc*ziabc)/rab2rbc2;
1872     term1 = (zbd*terma)/(termc*rbd2) + yabxbc/termc
1873     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term5*terma)/termss2;
1874     term2 = (3.0*zbd2*terma)/(termcrbd22)
1875     +(2.0*zbd*(yab*xbc-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic)))/(termc*rbd2)
1876     -terma/(termc*rbd2)
1877     +(2.0*zbc*zbd*terma)/(termc*rbd2*rbc2)
1878     +(2.0*zab*zbd*terma)/(termc*rbd2*rab2)
1879     -(term5*zbd*terma)/(termss*rbd2)
1880     +(2.0*zbc*(yab*xbc-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic)))/(termcrbc2)
1881     +(2.0*zab*(yab*xbc-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic)))/(termcrab2)
1882     -(term5*(yab*xbc-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic)))/(termss)
1883     +(3.0*zbc2*terma)/(termcrbc2*rbc2)
1884     -terma/(termcrbc2)
1885     +(2.0*zab*zbc*terma)/(termcrab2*rbc2)
1886     -terma/(termcrab2)
1887     +(3.0*zab2*terma)/(termcrab2*rab2)
1888     -(term5*zbc*terma)/(termcrbc2*sine)
1889     -(zab*term5*terma)/(termss*rab2)
1890     +(3.0*term5*term5*terma)/(terms4)
1891     -(term3*terma)/termss2;
1892     term4 = (-(2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(yab*xbc-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic))*terma)/(termc2)
1893     -(2.0*zbc*terma2)/(termc2*rbc2) - (2.0*zab*terma2)/(termc2*rab2)
1894     +(term5*terma2)/(termc*termss))
1895     *( (zbd*terma)/(termc*rbd2) + (yab*xbc-xab*ybc+xbd*(ybc-yab)+ybd*(xia-xic))/termc
1896     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term5*terma));
1897     dedzibzib = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1898    
1899     term3 = (8.0*zbc2*dot2)/(rrab2rbc22*rbc2) - (4.0*zab*zbc*dot)/(rrab2rbc22)
1900     +(4.0*ziabc*zbc*dot)/(rrab2rbc22) - (2.0*dot2)/(rrab2rbc22)
1901     +(4.0*zab*zbc*dot2)/(rab2rbc2*rab2rbc2) - (2.0*zab*ziabc)/(rab2rbc2)
1902     +(2.0*dot)/rab2rbc2 - (4.0*zab2*dot)/(rrab22rbc2);
1903     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1904     term1 = ( (yabxbd)/termc - (term1zca)/termss2 - (zbc*terma)/(termcrbc2))
1905     *( (zbd*terma)/(termc*rbd2) + (yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
1906     + (zbc*terma)/(termcrbc2) +(zab*terma)/(termcrab2) - (term5*terma)/termss2 );
1907     term2 = (zbd*(yabxbd))/(termc*rbd2)
1908     -(zbd*term1zca)/(2.0*termss*rbd2)
1909     -(zbc*zbd*terma)/(termcrbc2*rbd2)
1910     -(term1zc*(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/termss2
1911     +(zbc*(yabxbd))/(termcrbc2)
1912     -(zbc*(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termcrbc2)
1913     +(zab*(yabxbd))/(termcrab2)
1914     -(term5*(yabxbd))/termss2
1915     -(zbc*term1zca)/(2.0*termss*rbc2)
1916     -(zab*term1zca)/(2.0*termss*rab2)
1917     -(term3*terma)/termss2
1918     -(3.0*zbc2*terma)/(termcrbc2*rbc2)
1919     +terma/(termcrbc2)
1920     -(zab*zbc*terma)/(termcrab2*rbc2)
1921     +(3.0*term1zc*term5*terma)/(terms4)
1922     +(zbc*term5*terma)/(2.0*termss*rbc2);
1923     term4 = (-(2.0*(yabxbd)*terma)/(termc2) + (term1zca*terma)/(termc*termss)
1924     +(2.0*zbc*terma2)/(termc2*rbc2))
1925     *( (zbd*terma)/(termc*rbd2) + (yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
1926     +(zbc*terma)/(termcrbc2) + (zab*terma)/(termcrab2) - (term5*terma)/termss2);
1927     dedzibzic = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1928    
1929     term5 = -(2.0*zbc*dot2)/(rrab2rbc22) - (2.0*ziabc*dot)/rab2rbc2 - (2.0*zab*dot2)/(rrab22rbc2);
1930     term1 = ( -(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*( (zbd*terma)/(termc*rbd2)
1931     +(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc + (zbc*terma)/(termcrbc2)
1932     +(zab*terma)/(termcrab2) - (term5*terma)/termss2 );
1933     term2 = -(3.0*zbd2*terma)/(termcrbd22)
1934     +(zbd*(xabybc))/(termc*rbd2)
1935     -(zbd*(yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic)))/(termc*rbd2)
1936     +terma/(termc*rbd2)
1937     -(zbc*zbd*terma)/(termcrbc2*rbd2)
1938     -(zab*zbd*terma)/(termcrab2*rbd2)
1939     +(term5*zbd*terma)/(2.0*termss*rbd2)
1940     +(zbc*(xabybc))/(termcrbc2)
1941     +(zab*(xabybc))/(termcrab2)
1942     -(term5*(xabybc))/termss2;
1943     term4 = ((2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1944     *((zbd*terma)/(termc*rbd2) + (yab*xbc-xab*ybc+xbd*(yic-yia)+ybd*(xia-xic))/termc
1945     +(zbc*terma)/(termcrbc2) +(zab*terma)/(termcrab2) - (term5*terma)/termss2);
1946     dedzibzid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1947    
1948     term1 = term1aa(zbc, yabxbd, terma, term1zc, termc, termss2, termcrbc2);
1949     term3 = (8.0*zab*zbc*dot)/(rrab2rbc22)- (2.0*zab2)/rab2rbc2 - (8.0*zbc2*dot2)/(rbc2*rrab2rbc22)
1950     + (2.0*dot2)/(rrab2rbc22);
1951     term2 = - (term1zc*(yabxbd))/(termss)
1952     - (2.0*zbc*(yabxbd))/(termcrbc2)
1953     + (3.0*terma*term1zc*term1zc)/(terms4)
1954     - (terma*term3)/termss2
1955     + (zbc*term1zca)/(termss*rbc2)
1956     + (3.0*zbc2*terma)/(rbc2*rbc2*termc)
1957     - terma/(rbc2*termc);
1958     term4 = ( - (2.0*terma*(yabxbd))/(termc2)
1959     + (term1zca*terma)/(termc*termss)
1960     + (2.0*zbc*terma2)/(rbc2*termc2) )
1961     *( (yabxbd)/termc - (term1zca)/termss2
1962     - (zbc*terma)/(rbc2*termc));
1963     dedziczic = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1964    
1965     term1 = (-(zbd*terma)/(termc*rbd2) + (xabybc)/termc)*( (yabxbd)/termc
1966     -(term1zca)/termss2 - (zbc*terma)/(termcrbc2));
1967     term2 = -(zbd*(yabxbd))/(rbd2*termc) + (term1zc*zbd*terma)/(2.0*termss*rbd2)
1968     +(zbc*zbd*terma)/(rbc2*rbd2*termc) - (term1zc*(xabybc))/termss2
1969     -(zbc*(xabybc))/(termcrbc2);
1970     term4 = ( (2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1971     *( (yabxbd)/termc - (term1zca)/termss2 - (zbc*terma)/(rbc2*termc));
1972     dedziczid = ddtdenom*term1 + ddt1denom*term2 - ddt2denom*term4;
1973    
1974     term1 = -(zbd*terma)/(termc*rbd2) + (xabybc)/termc;
1975     term2 = (3.0*zbd2*terma)/(termcrbd22) - (2.0*(xabybc)*zbd)/(termc*rbd2)
1976     -(terma)/(termc*rbd2);
1977     term4 = ((2.0*zbd*terma2)/(termc2*rbd2) - (2.0*(xabybc)*terma)/(termc2))
1978     *( -(zbd*terma)/(termc*rbd2) + (xabybc)/termc );
1979     dedzidzid = ddtdenom*term1*term1 + ddt1denom*term2 - ddt2denom*term4;
1980    
1981     if (iatom == ia)
1982     {
1983     hess.hessx[ia][0] += dedxiaxia ;
1984     hess.hessx[ia][1] += dedxiayia ;
1985     hess.hessx[ia][2] += dedxiazia ;
1986    
1987     hess.hessy[ia][0] += dedxiayia ;
1988     hess.hessy[ia][1] += dedyiayia ;
1989     hess.hessy[ia][2] += dedyiazia ;
1990    
1991     hess.hessz[ia][0] += dedxiazia ;
1992     hess.hessz[ia][1] += dedyiazia ;
1993     hess.hessz[ia][2] += dedziazia ;
1994    
1995     hess.hessx[ib][0] += dedxiaxib ;
1996     hess.hessx[ib][1] += dedxiayib ;
1997     hess.hessx[ib][2] += dedxiazib ;
1998    
1999     hess.hessy[ib][0] += dedxibyia ;
2000     hess.hessy[ib][1] += dedyiayib ;
2001     hess.hessy[ib][2] += dedyiazib ;
2002    
2003     hess.hessz[ib][0] += dedxibzia ;
2004     hess.hessz[ib][1] += dedyibzia ;
2005     hess.hessz[ib][2] += dedziazib ;
2006    
2007     hess.hessx[ic][0] += dedxiaxic ;
2008     hess.hessy[ic][0] += dedxicyia ;
2009     hess.hessz[ic][0] += dedxiczia ;
2010    
2011     hess.hessx[ic][1] += dedxiayic ;
2012     hess.hessy[ic][1] += dedyiayic ;
2013     hess.hessz[ic][1] += dedyiczia ;
2014    
2015     hess.hessx[ic][2] += dedxiazic ;
2016     hess.hessy[ic][2] += dedyiazic ;
2017     hess.hessz[ic][2] += dedziazic ;
2018    
2019     hess.hessx[id][0] += dedxiaxid ;
2020     hess.hessy[id][0] += dedxidyia ;
2021     hess.hessz[id][0] += dedxidzia ;
2022     hess.hessx[id][1] += dedxiayid ;
2023     hess.hessy[id][1] += dedyiayid ;
2024     hess.hessz[id][1] += dedyidzia ;
2025     hess.hessx[id][2] += dedxiazid ;
2026     hess.hessy[id][2] += dedyiazid ;
2027     hess.hessz[id][2] += dedziazid ;
2028     }else if (iatom == ib)
2029     {
2030     hess.hessx[ib][0] += dedxibxib ;
2031     hess.hessy[ib][0] += dedxibyib ;
2032     hess.hessz[ib][0] += dedxibzib ;
2033     hess.hessx[ib][1] += dedxibyib ;
2034     hess.hessy[ib][1] += dedyibyib ;
2035     hess.hessz[ib][1] += dedyibzib ;
2036     hess.hessx[ib][2] += dedxibzib ;
2037     hess.hessy[ib][2] += dedyibzib ;
2038     hess.hessz[ib][2] += dedzibzib ;
2039    
2040     hess.hessx[ia][0] += dedxiaxib ;
2041     hess.hessy[ia][0] += dedxiayib ;
2042     hess.hessz[ia][0] += dedxiazib ;
2043     hess.hessx[ia][1] += dedxibyia ;
2044     hess.hessy[ia][1] += dedyiayib ;
2045     hess.hessz[ia][1] += dedyiazib ;
2046     hess.hessx[ia][2] += dedxibzia ;
2047     hess.hessy[ia][2] += dedyibzia ;
2048     hess.hessz[ia][2] += dedziazib ;
2049    
2050     hess.hessx[ic][0] += dedxibxic ;
2051     hess.hessy[ic][0] += dedxicyib ;
2052     hess.hessz[ic][0] += dedxiczib ;
2053    
2054     hess.hessx[ic][1] += dedxibyic ;
2055     hess.hessy[ic][1] += dedyibyic ;
2056     hess.hessz[ic][1] += dedyiczib ;
2057    
2058     hess.hessx[ic][2] += dedxibzic ;
2059     hess.hessy[ic][2] += dedyibzic ;
2060     hess.hessz[ic][2] += dedzibzic ;
2061    
2062     hess.hessx[id][0] += dedxibxid ;
2063     hess.hessy[id][0] += dedxidyib ;
2064     hess.hessz[id][0] += dedxidzib ;
2065    
2066     hess.hessx[id][1] += dedxibyid ;
2067     hess.hessy[id][1] += dedyibyid ;
2068     hess.hessz[id][1] += dedyidzib ;
2069    
2070     hess.hessx[id][2] += dedxibzid ;
2071     hess.hessy[id][2] += dedyibzid ;
2072     hess.hessz[id][2] += dedzibzid ;
2073     }else if (iatom == ic)
2074     {
2075     hess.hessx[ic][0] += dedxicxic ;
2076     hess.hessy[ic][0] += dedxicyic ;
2077     hess.hessz[ic][0] += dedxiczic ;
2078     hess.hessx[ic][1] += dedxicyic ;
2079     hess.hessy[ic][1] += dedyicyic ;
2080     hess.hessz[ic][1] += dedyiczic ;
2081     hess.hessx[ic][2] += dedxiczic ;
2082     hess.hessy[ic][2] += dedyiczic ;
2083     hess.hessz[ic][2] += dedziczic ;
2084     hess.hessx[ia][0] += dedxiaxic ;
2085     hess.hessy[ia][0] += dedxiayic ;
2086     hess.hessz[ia][0] += dedxiazic ;
2087     hess.hessx[ia][1] += dedxicyia ;
2088     hess.hessy[ia][1] += dedyiayic ;
2089     hess.hessz[ia][1] += dedyiazic ;
2090     hess.hessx[ia][2] += dedxiczia ;
2091     hess.hessy[ia][2] += dedyiczia ;
2092     hess.hessz[ia][2] += dedziazic ;
2093     hess.hessx[ib][0] += dedxibxic ;
2094     hess.hessy[ib][0] += dedxibyic ;
2095     hess.hessz[ib][0] += dedxibzic ;
2096     hess.hessx[ib][1] += dedxicyib ;
2097     hess.hessy[ib][1] += dedyibyic ;
2098     hess.hessz[ib][1] += dedyibzic ;
2099     hess.hessx[ib][2] += dedxiczib ;
2100     hess.hessy[ib][2] += dedyiczib ;
2101     hess.hessz[ib][2] += dedzibzic ;
2102     hess.hessx[id][0] += dedxicxid ;
2103     hess.hessy[id][0] += dedxidyic ;
2104     hess.hessz[id][0] += dedxidzic ;
2105     hess.hessx[id][1] += dedxicyid ;
2106     hess.hessy[id][1] += dedyicyid ;
2107     hess.hessz[id][1] += dedyidzic ;
2108     hess.hessx[id][2] += dedxiczid ;
2109     hess.hessy[id][2] += dedyiczid ;
2110     hess.hessz[id][2] += dedziczid ;
2111     }else if (iatom == id)
2112     {
2113     hess.hessx[id][0] += dedxidxid ;
2114     hess.hessy[id][0] += dedxidyid ;
2115     hess.hessz[id][0] += dedxidzid ;
2116    
2117     hess.hessx[id][1] += dedxidyid ;
2118     hess.hessy[id][1] += dedyidyid ;
2119     hess.hessz[id][1] += dedyidzid ;
2120    
2121     hess.hessx[id][2] += dedxidzid ;
2122     hess.hessy[id][2] += dedyidzid ;
2123     hess.hessz[id][2] += dedzidzid ;
2124    
2125     hess.hessx[ia][0] += dedxiaxid ;
2126     hess.hessy[ia][0] += dedxiayid ;
2127     hess.hessz[ia][0] += dedxiazid ;
2128     hess.hessx[ia][1] += dedxidyia ;
2129     hess.hessy[ia][1] += dedyiayid ;
2130     hess.hessz[ia][1] += dedyiazid ;
2131     hess.hessx[ia][2] += dedxidzia ;
2132     hess.hessy[ia][2] += dedyidzia ;
2133     hess.hessz[ia][2] += dedziazid ;
2134    
2135     hess.hessx[ib][0] += dedxibxid ;
2136     hess.hessy[ib][0] += dedxibyid ;
2137     hess.hessz[ib][0] += dedxibzid ;
2138     hess.hessx[ib][1] += dedxidyib ;
2139     hess.hessy[ib][1] += dedyibyid ;
2140     hess.hessz[ib][1] += dedyibzid ;
2141     hess.hessx[ib][2] += dedxidzib ;
2142     hess.hessy[ib][2] += dedyidzib ;
2143     hess.hessz[ib][2] += dedzibzid ;
2144    
2145     hess.hessx[ic][0] += dedxicxid ;
2146     hess.hessy[ic][0] += dedxicyid ;
2147     hess.hessz[ic][0] += dedxiczid ;
2148     hess.hessx[ic][1] += dedxidyic ;
2149     hess.hessy[ic][1] += dedyicyid ;
2150     hess.hessz[ic][1] += dedyiczid ;
2151     hess.hessx[ic][2] += dedxidzic ;
2152     hess.hessy[ic][2] += dedyidzic ;
2153     hess.hessz[ic][2] += dedziczid ;
2154     }
2155    
2156    
2157     }
2158    
2159     }
2160     }
2161     }
2162     }