ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/eopbend.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 16965 byte(s)
Log Message:
test

Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "energies.h"
7 #include "angles.h"
8 #include "derivs.h"
9 #include "hess.h"
10 #include "utility.h"
11
12 void eopbend2b(int);
13
14 EXTERN struct t_units {
15 double bndunit, cbnd, qbnd;
16 double angunit, cang, qang, pang, sang, aaunit;
17 double stbnunit, ureyunit, torsunit, storunit, v14scale;
18 double aterm, bterm, cterm, dielec, chgscale;
19 } units;
20
21 EXTERN struct t_minim_values {
22 int iprint, ndc, nconst;
23 float dielc;
24 } minim_values;
25
26 void eopbend()
27 {
28 int i,ia,ib,ic,id;
29 double e,angle,dot,cosine;
30 double dt,dt2,dt3,dt4;
31 double xia,yia,zia;
32 double xib,yib,zib;
33 double xic,yic,zic;
34 double xid,yid,zid;
35 double xad,yad,zad;
36 double xbd,ybd,zbd,rbd2;
37 double xcd,ycd,zcd;
38 double xpd,ypd,zpd,rpd2;
39 double xt,yt,zt,rt2,delta;
40
41 energies.eopb = 0.0;
42
43 if (minim_values.iprint)
44 {
45 fprintf(pcmoutfile,"\nOut of Plane Bending Terms\n");
46 fprintf(pcmoutfile," At1 At2 At3 At4 Angle Oopconst Eoop\n");
47 }
48
49 for (i=0; i < angles.nang; i++)
50 {
51 if (angles.angin[i] == TRUE)
52 {
53 ia = angles.i13[i][0];
54 ib = angles.i13[i][1];
55 ic = angles.i13[i][2];
56 id = angles.i13[i][3];
57 if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use)
58 {
59 xia = atom[ia].x;
60 yia = atom[ia].y;
61 zia = atom[ia].z;
62 xib = atom[ib].x;
63 yib = atom[ib].y;
64 zib = atom[ib].z;
65 xic = atom[ic].x;
66 yic = atom[ic].y;
67 zic = atom[ic].z;
68 xid = atom[id].x;
69 yid = atom[id].y;
70 zid = atom[id].z;
71 xad = xia - xid;
72 yad = yia - yid;
73 zad = zia - zid;
74 xbd = xib - xid;
75 ybd = yib - yid;
76 zbd = zib - zid;
77 xcd = xic - xid;
78 ycd = yic - yid;
79 zcd = zic - zid;
80 xt = yad*zcd - zad*ycd;
81 yt = zad*xcd - xad*zcd;
82 zt = xad*ycd - yad*xcd;
83 rt2 = xt*xt + yt*yt + zt*zt;
84 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
85 xpd = xbd + xt*delta;
86 ypd = ybd + yt*delta;
87 zpd = zbd + zt*delta;
88 rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
89 rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
90 if (rbd2 != 0.0 && rpd2 != 0.0)
91 {
92 dot = xbd*xpd + ybd*ypd + zbd*zpd;
93 cosine = dot / sqrt(rbd2*rpd2);
94 if (cosine < -1.0)
95 cosine = -1.0;
96 if (cosine > 1.0)
97 cosine = 1.0;
98 angle = radian * acos(cosine);
99 dt = angle;
100 dt2 = dt * dt;
101 dt3 = dt2 * dt;
102 dt4 = dt2 * dt2;
103 e = units.angunit * angles.copb[i] * dt2
104 * (1.0+units.cang*dt+units.qang*dt2+units.pang*dt3+units.sang*dt4);
105 energies.eopb += e;
106 atom[ia].energy += e;
107 atom[ib].energy += e;
108 atom[ic].energy += e;
109 atom[id].energy += e;
110 if (minim_values.iprint)
111 fprintf(pcmoutfile,"Oop: %-4d- %-4d- %-4d- %-4d %-8.3f %-8.3f = %-8.4f\n",ia,ib,ic,id,
112 angle, angles.copb[i],e);
113 }
114
115 }
116 }
117 }
118 }
119
120 void eopbend1()
121 {
122 int i,ia,ib,ic,id;
123 double e,angle,dot,cosine;
124 double dt,dt2,dt3,dt4;
125 double deddt,termb,termp,term;
126 double xia,yia,zia;
127 double xib,yib,zib;
128 double xic,yic,zic;
129 double xid,yid,zid;
130 double xad,yad,zad;
131 double xbd,ybd,zbd,rbd2;
132 double xcd,ycd,zcd;
133 double xpd,ypd,zpd,rpd2;
134 double xt,yt,zt,rt2,ptrt2;
135 double xm,ym,zm,rm,delta,delta2;
136 double dedxia,dedyia,dedzia;
137 double dedxib,dedyib,dedzib;
138 double dedxic,dedyic,dedzic;
139 double dedxid,dedyid,dedzid;
140 double dedxip,dedyip,dedzip;
141 double dpdxia,dpdyia,dpdzia;
142 double dpdxic,dpdyic,dpdzic;
143
144 energies.eopb = 0.0;
145 for (i=0; i <= natom; i++)
146 {
147 deriv.deopb[i][0] = 0.0;
148 deriv.deopb[i][1] = 0.0;
149 deriv.deopb[i][2] = 0.0;
150 }
151
152 for (i=0; i < angles.nang; i++)
153 {
154 if (angles.angin[i] == TRUE)
155 {
156 ia = angles.i13[i][0];
157 ib = angles.i13[i][1];
158 ic = angles.i13[i][2];
159 id = angles.i13[i][3];
160 if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use)
161 {
162 xia = atom[ia].x;
163 yia = atom[ia].y;
164 zia = atom[ia].z;
165 xib = atom[ib].x;
166 yib = atom[ib].y;
167 zib = atom[ib].z;
168 xic = atom[ic].x;
169 yic = atom[ic].y;
170 zic = atom[ic].z;
171 xid = atom[id].x;
172 yid = atom[id].y;
173 zid = atom[id].z;
174 xad = xia - xid;
175 yad = yia - yid;
176 zad = zia - zid;
177 xbd = xib - xid;
178 ybd = yib - yid;
179 zbd = zib - zid;
180 xcd = xic - xid;
181 ycd = yic - yid;
182 zcd = zic - zid;
183 xt = yad*zcd - zad*ycd;
184 yt = zad*xcd - xad*zcd;
185 zt = xad*ycd - yad*xcd;
186 rt2 = xt*xt + yt*yt + zt*zt;
187 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
188 xpd = xbd + xt*delta;
189 ypd = ybd + yt*delta;
190 zpd = zbd + zt*delta;
191 rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
192 rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
193 xm = ypd*zbd - zpd*ybd;
194 ym = zpd*xbd - xpd*zbd;
195 zm = xpd*ybd - ypd*xbd;
196 rm = sqrt(xm*xm + ym*ym + zm*zm);
197 if (rm != 0.0)
198 {
199 dot = xbd*xpd + ybd*ypd + zbd*zpd;
200 cosine = dot / sqrt(rbd2*rpd2);
201 if (cosine < -1.0)
202 cosine = -1.0;
203 if (cosine > 1.0)
204 cosine = 1.0;
205 angle = radian * acos(cosine);
206 dt = angle;
207 dt2 = dt * dt;
208 dt3 = dt2 * dt;
209 dt4 = dt2 * dt2;
210 e = units.angunit * angles.copb[i] * dt2
211 * (1.0+units.cang*dt+units.qang*dt2+units.pang*dt3+units.sang*dt4);
212
213 deddt = units.angunit * angles.copb[i] * dt
214 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
215 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
216 deddt = deddt * radian;
217 termb = -deddt / (rbd2*rm);
218 termp = deddt / (rpd2*rm);
219 dedxib = termb * (ybd*zm-zbd*ym);
220 dedyib = termb * (zbd*xm-xbd*zm);
221 dedzib = termb * (xbd*ym-ybd*xm);
222 dedxip = termp * (ypd*zm-zpd*ym);
223 dedyip = termp * (zpd*xm-xpd*zm);
224 dedzip = termp * (xpd*ym-ypd*xm);
225
226 delta2 = 2.0 * delta;
227 ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
228 term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
229 dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
230 term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
231 dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
232 term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
233 dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
234 term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
235 dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
236 term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
237 dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
238 term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
239 dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
240
241 dedxia = dpdxia;
242 dedyia = dpdyia;
243 dedzia = dpdzia;
244 dedxic = dpdxic;
245 dedyic = dpdyic;
246 dedzic = dpdzic;
247 dedxid = -dedxia - dedxib - dedxic;
248 dedyid = -dedyia - dedyib - dedyic;
249 dedzid = -dedzia - dedzib - dedzic;
250
251 energies.eopb += e;
252 deriv.deopb[ia][0] += dedxia;
253 deriv.deopb[ia][1] += dedyia;
254 deriv.deopb[ia][2] += dedzia;
255
256 deriv.deopb[ib][0] += dedxib;
257 deriv.deopb[ib][1] += dedyib;
258 deriv.deopb[ib][2] += dedzib;
259
260 deriv.deopb[ic][0] += dedxic;
261 deriv.deopb[ic][1] += dedyic;
262 deriv.deopb[ic][2] += dedzic;
263
264 deriv.deopb[id][0] += dedxid;
265 deriv.deopb[id][1] += dedyid;
266 deriv.deopb[id][2] += dedzid;
267 virial.virx += xad*dedxia + xbd*dedxib + xcd*dedxic;
268 virial.viry += yad*dedyia + ybd*dedyib + ycd*dedyic;
269 virial.virz += zad*dedzia + zbd*dedzib + zcd*dedzic;
270 }
271
272 }
273 }
274 }
275 }
276
277 void eopbend2(int i)
278 {
279 int j,k;
280 int ia,ib,ic,id;
281 double old,eps;
282 double **d0; // d0[MAXATOM][3];
283
284 d0 = dmatrix(0,natom,0,3);
285 eps = 1.0e-7;
286
287 for (k=0; k < angles.nang; k++)
288 {
289 if (angles.angin[k] == TRUE)
290 {
291 ia = angles.i13[k][0];
292 ib = angles.i13[k][1];
293 ic = angles.i13[k][2];
294 id = angles.i13[k][3];
295 if (ia == i || ib == i || ic == i || id == i)
296 {
297 eopbend2b(k);
298 for (j = 0; j < 3; j++)
299 {
300 d0[ia][j] = deriv.deopb[ia][j];
301 d0[ib][j] = deriv.deopb[ib][j];
302 d0[ic][j] = deriv.deopb[ic][j];
303 d0[id][j] = deriv.deopb[id][j];
304 }
305 // numerical x
306 old = atom[i].x;
307 atom[i].x += eps;
308 eopbend2b(k);
309 atom[i].x = old;
310 for (j=0; j < 3; j++)
311 {
312 hess.hessx[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
313 hess.hessx[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
314 hess.hessx[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
315 hess.hessx[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
316 }
317 // numerical y
318 old = atom[i].y;
319 atom[i].y += eps;
320 eopbend2b(k);
321 atom[i].y = old;
322 for (j=0; j < 3; j++)
323 {
324 hess.hessy[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
325 hess.hessy[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
326 hess.hessy[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
327 hess.hessy[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
328 }
329 // numerical z
330 old = atom[i].z;
331 atom[i].z += eps;
332 eopbend2b(k);
333 atom[i].z = old;
334 for (j=0; j < 3; j++)
335 {
336 hess.hessz[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
337 hess.hessz[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
338 hess.hessz[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
339 hess.hessz[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
340 }
341 }
342 }
343 }
344 free_dmatrix( d0 ,0,natom,0,3);
345 }
346
347 void eopbend2b(int k)
348 {
349 int j,ia,ib,ic,id;
350 double angle,dot,cosine;
351 double dt,dt2,dt3,dt4;
352 double deddt,termb,termp,term;
353 double xia,yia,zia;
354 double xib,yib,zib;
355 double xic,yic,zic;
356 double xid,yid,zid;
357 double xad,yad,zad;
358 double xbd,ybd,zbd,rbd2;
359 double xcd,ycd,zcd;
360 double xpd,ypd,zpd,rpd2;
361 double xt,yt,zt,rt2,ptrt2;
362 double xm,ym,zm,rm,delta,delta2;
363 double dedxia,dedyia,dedzia;
364 double dedxib,dedyib,dedzib;
365 double dedxic,dedyic,dedzic;
366 double dedxid,dedyid,dedzid;
367 double dedxip,dedyip,dedzip;
368 double dpdxia,dpdyia,dpdzia;
369 double dpdxic,dpdyic,dpdzic;
370
371 ia = angles.i13[k][0];
372 ib = angles.i13[k][1];
373 ic = angles.i13[k][2];
374 id = angles.i13[k][3];
375
376 xia = atom[ia].x;
377 yia = atom[ia].y;
378 zia = atom[ia].z;
379 xib = atom[ib].x;
380 yib = atom[ib].y;
381 zib = atom[ib].z;
382 xic = atom[ic].x;
383 yic = atom[ic].y;
384 zic = atom[ic].z;
385 xid = atom[id].x;
386 yid = atom[id].y;
387 zid = atom[id].z;
388
389 for (j=0; j < 3; j++)
390 {
391 deriv.deopb[ia][j] = 0.0;
392 deriv.deopb[ib][j] = 0.0;
393 deriv.deopb[ic][j] = 0.0;
394 deriv.deopb[id][j] = 0.0;
395 }
396
397 xad = xia - xid;
398 yad = yia - yid;
399 zad = zia - zid;
400 xbd = xib - xid;
401 ybd = yib - yid;
402 zbd = zib - zid;
403 xcd = xic - xid;
404 ycd = yic - yid;
405 zcd = zic - zid;
406 xt = yad*zcd - zad*ycd;
407 yt = zad*xcd - xad*zcd;
408 zt = xad*ycd - yad*xcd;
409 rt2 = xt*xt + yt*yt + zt*zt;
410 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
411 xpd = xbd + xt*delta;
412 ypd = ybd + yt*delta;
413 zpd = zbd + zt*delta;
414 rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
415 rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
416 xm = ypd*zbd - zpd*ybd;
417 ym = zpd*xbd - xpd*zbd;
418 zm = xpd*ybd - ypd*xbd;
419 rm = sqrt(xm*xm + ym*ym + zm*zm);
420 if (rm != 0.0)
421 {
422 dot = xbd*xpd + ybd*ypd + zbd*zpd;
423 cosine = dot / sqrt(rbd2*rpd2);
424 if (cosine < -1.0)
425 cosine = -1.0;
426 if (cosine > 1.0)
427 cosine = 1.0;
428 angle = radian * acos(cosine);
429 dt = angle;
430 dt2 = dt * dt;
431 dt3 = dt2 * dt;
432 dt4 = dt2 * dt2;
433 deddt = units.angunit * angles.copb[k] * dt
434 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
435 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
436
437 deddt = deddt * radian;
438 termb = -deddt / (rbd2*rm);
439 termp = deddt / (rpd2*rm);
440 dedxib = termb * (ybd*zm-zbd*ym);
441 dedyib = termb * (zbd*xm-xbd*zm);
442 dedzib = termb * (xbd*ym-ybd*xm);
443 dedxip = termp * (ypd*zm-zpd*ym);
444 dedyip = termp * (zpd*xm-xpd*zm);
445 dedzip = termp * (xpd*ym-ypd*xm);
446
447 delta2 = 2.0 * delta;
448 ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
449 term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
450 dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
451 term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
452 dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
453 term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
454 dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
455 term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
456 dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
457 term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
458 dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
459 term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
460 dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
461
462 dedxia = dpdxia;
463 dedyia = dpdyia;
464 dedzia = dpdzia;
465 dedxic = dpdxic;
466 dedyic = dpdyic;
467 dedzic = dpdzic;
468 dedxid = -dedxia - dedxib - dedxic;
469 dedyid = -dedyia - dedyib - dedyic;
470 dedzid = -dedzia - dedzib - dedzic;
471
472 deriv.deopb[ia][0] = dedxia;
473 deriv.deopb[ia][1] = dedyia;
474 deriv.deopb[ia][2] = dedzia;
475 deriv.deopb[ib][0] = dedxib;
476 deriv.deopb[ib][1] = dedyib;
477 deriv.deopb[ib][2] = dedzib;
478 deriv.deopb[ic][0] = dedxic;
479 deriv.deopb[ic][1] = dedyic;
480 deriv.deopb[ic][2] = dedzic;
481 deriv.deopb[id][0] = dedxid;
482 deriv.deopb[id][1] = dedyid;
483 deriv.deopb[id][2] = dedzid;
484 }
485 }