ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/GhostBend.cpp
Revision: 301
Committed: Mon Mar 10 14:46:17 2003 UTC (21 years, 6 months ago) by gezelter
File size: 2467 byte(s)
Log Message:
GhostBend addition

File Contents

# Content
1 #include "SRI.hpp"
2 #include "Atom.hpp"
3
4 #include <math.h>
5 #include <iostream>
6 #include <stdlib.h>
7
8 void GhostBend::set_atoms( Atom &a, Atom &b ){
9
10 c_p_a = &a;
11 c_p_b = &b;
12
13 c_potential_E = 0.0;
14 }
15
16
17 void GhostBend::calc_forces(){
18
19 double dx,dy,dz,gx,gy,gz,dx2,dy2,dz2,gx2,gy2,gz2;
20 double rij2, rkj2, riji2, rkji2, dot, denom, cosang, angl;
21
22 double sina2, sinai;
23
24 double comf2, comf3, comf4;
25 double dcsidx, dcsidy, dcsidz, dcskdx, dcskdy, dcskdz;
26 // double dcsjdx, dcsjdy, dcsjdz;
27 double dadxi, dadyi, dadzi;
28 double dadxk, dadyk, dadzk;//, dadxj, dadyj, dadzj;
29 double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj;
30 double u[3];
31
32 dx = c_p_a->getX() - c_p_b->getX();
33 dy = c_p_a->getY() - c_p_b->getY();
34 dz = c_p_a->getZ() - c_p_b->getZ();
35
36 c_p_b->getU(u);
37
38 gx = u[0];
39 gy = u[1];
40 gz = u[2];
41
42 dx2 = dx * dx;
43 dy2 = dy * dy;
44 dz2 = dz * dz;
45
46 gx2 = gx * gx;
47 gy2 = gy * gy;
48 gz2 = gz * gz;
49
50 rij2 = dx2 + dy2 + dz2;
51 rkj2 = gx2 + gy2 + gz2;
52
53 riji2 = 1.0 / rij2;
54 rkji2 = 1.0 / rkj2;
55
56 dot = dx * gx + dy * gy + dz * gz;
57 denom = sqrt((riji2 * rkji2));
58 cosang = dot * denom;
59
60 if(cosang > 1.0)cosang = 1.0;
61 if(cosang < -1.0) cosang = -1.0;
62
63 angl = acos(cosang);
64 angl = angl * 180.0 / M_PI;
65
66 sina2 = 1.0 - cosang*cosang;
67 if(fabs(sina2) < 1.0E-12 ) sina2 = 1.0E-12;
68 sinai = 1.0 / sqrt(sina2);
69
70 comf2 = cosang * riji2;
71 comf3 = cosang * rkji2;
72 comf4 = bend_force(angl);
73
74 dcsidx = gx*denom - comf2*dx;
75 dcsidy = gy*denom - comf2*dy;
76 dcsidz = gz*denom - comf2*dz;
77
78 dcskdx = dx*denom - comf3*gx;
79 dcskdy = dy*denom - comf3*gy;
80 dcskdz = dz*denom - comf3*gz;
81
82 // dcsjdx = -dcsidx - dcskdx;
83 // dcsjdy = -dcsidy - dcskdy;
84 // dcsjdz = -dcsidz - dcskdz;
85
86 dadxi = -sinai*dcsidx;
87 dadyi = -sinai*dcsidy;
88 dadzi = -sinai*dcsidz;
89
90 dadxk = -sinai*dcskdx;
91 dadyk = -sinai*dcskdy;
92 dadzk = -sinai*dcskdz;
93
94 // dadxj = -dadxi - dadxk;
95 // dadyj = -dadyi - dadyk;
96 // dadzj = -dadzi - dadzk;
97
98 daxi = comf4*dadxi;
99 dayi = comf4*dadyi;
100 dazi = comf4*dadzi;
101
102 daxk = comf4*dadxk;
103 dayk = comf4*dadyk;
104 dazk = comf4*dadzk;
105
106 daxj = -daxi - daxk;
107 dayj = -dayi - dayk;
108 dazj = -dazi - dazk;
109
110 c_p_a->addFx(daxi);
111 c_p_a->addFy(dayi);
112 c_p_a->addFz(dazi);
113
114 c_p_b->addFx(daxj + daxk);
115 c_p_b->addFy(dayj + dayk);
116 c_p_b->addFz(dazj + dazk);
117
118 c_p_b->addTx(gy*dazk - gz*dayk);
119 c_p_b->addTy(gz*daxk - gx*dazk);
120 c_p_b->addTz(gx*dayk - gy*daxk);
121
122 return;
123 }