ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/primitives/GhostBend.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/primitives/GhostBend.cpp (file contents):
Revision 1691, Thu Oct 28 22:34:02 2004 UTC vs.
Revision 1692 by tim, Mon Nov 1 20:15:58 2004 UTC

# Line 1 | Line 1
1 <
2 < #include <math.h>
3 < #include <iostream>
4 < #include <stdlib.h>
5 <
6 < #include "utils/simError.h"
7 < #include "primitives/SRI.hpp"
8 < #include "primitives/Atom.hpp"
9 <
10 <
11 <
12 < GhostBend::GhostBend( Atom &a, Atom &b ){
13 <  
14 <  c_p_a = &a;
15 <  
16 <  if( !b.isDirectional() ){
17 <    
18 <    // if atom b is not directional, then bad things will happen
19 <    
20 <    sprintf( painCave.errMsg,
21 <             " Ghost Bend error: Atom # %d of type \"%s\" is not "
22 <             "directional.\n",
23 <             b.getIndex(),
24 <             b.getType() );
25 <    painCave.isFatal = 1;
26 <    simError();
27 <  }    
28 <
29 <  atomB = ( DirectionalAtom* ) &b;
30 <  
31 <  c_potential_E = 0.0;
32 <
33 < }
34 <
35 <
36 < void GhostBend::calc_forces(){
37 <  
38 <  double dx,dy,dz,gx,gy,gz,dx2,dy2,dz2,gx2,gy2,gz2;
39 <  double rij2, rkj2, riji2, rkji2, dot, denom, cosang, angl;
40 <  
41 <  double sina2, sinai;
42 <
43 <  double comf2, comf3, comf4;
44 <  double dcsidx, dcsidy, dcsidz, dcskdx, dcskdy, dcskdz;
45 <  // double dcsjdx, dcsjdy, dcsjdz;
46 <  double dadxi, dadyi, dadzi;
47 <  double dadxk, dadyk, dadzk;//, dadxj, dadyj, dadzj;
48 <  double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj;
49 <  double u[3];
50 <  
51 <  double aR[3], bR[3];
52 <  double aF[3], bF[3], bTrq[3];
53 <
54 <  c_p_a->getPos( aR );
55 <  atomB->getPos( bR );
56 <  
57 <
58 <  dx = aR[0] - bR[0];
59 <  dy = aR[1] - bR[1];
60 <  dz = aR[2] - bR[2];
61 <
62 <  atomB->getU(u);
63 <
64 <  gx = u[0];
65 <  gy = u[1];
66 <  gz = u[2];
67 <  
68 <  dx2 = dx * dx;
69 <  dy2 = dy * dy;
70 <  dz2 = dz * dz;
71 <
72 <  gx2 = gx * gx;
73 <  gy2 = gy * gy;
74 <  gz2 = gz * gz;
75 <  
76 <  rij2 = dx2 + dy2 + dz2;
77 <  rkj2 = gx2 + gy2 + gz2;
78 <  
79 <  riji2 = 1.0 / rij2;
80 <  rkji2 = 1.0 / rkj2;
81 <
82 <  dot = dx * gx + dy * gy + dz * gz;
83 <  denom = sqrt((riji2 * rkji2));
84 <  cosang = dot * denom;
85 <
86 <  if(cosang > 1.0)cosang = 1.0;
87 <  if(cosang < -1.0) cosang = -1.0;
88 <
89 <  angl = acos(cosang);
90 <  angl = angl * 180.0 / M_PI;
91 <
92 <  sina2 = 1.0 - cosang*cosang;
93 <  if(fabs(sina2) < 1.0E-12 ) sina2 = 1.0E-12;
94 <  sinai = 1.0 / sqrt(sina2);
95 <
96 <  comf2 = cosang * riji2;
97 <  comf3 = cosang * rkji2;
98 <  comf4 = bend_force(angl);
99 <
100 <  dcsidx = gx*denom - comf2*dx;
101 <  dcsidy = gy*denom - comf2*dy;
102 <  dcsidz = gz*denom - comf2*dz;
103 <  
104 <  dcskdx = dx*denom - comf3*gx;
105 <  dcskdy = dy*denom - comf3*gy;
106 <  dcskdz = dz*denom - comf3*gz;
107 <  
108 < //   dcsjdx = -dcsidx - dcskdx;
109 < //   dcsjdy = -dcsidy - dcskdy;
110 < //   dcsjdz = -dcsidz - dcskdz;
111 <
112 <  dadxi = -sinai*dcsidx;
113 <  dadyi = -sinai*dcsidy;
114 <  dadzi = -sinai*dcsidz;
115 <
116 <  dadxk = -sinai*dcskdx;
117 <  dadyk = -sinai*dcskdy;
118 <  dadzk = -sinai*dcskdz;
119 <
120 < //   dadxj = -dadxi - dadxk;
121 < //   dadyj = -dadyi - dadyk;
122 < //   dadzj = -dadzi - dadzk;
123 <
124 <  daxi = comf4*dadxi;
125 <  dayi = comf4*dadyi;
126 <  dazi = comf4*dadzi;
127 <
128 <  daxk = comf4*dadxk;
129 <  dayk = comf4*dadyk;
130 <  dazk = comf4*dadzk;
131 <  
132 <  daxj = -daxi - daxk;
133 <  dayj = -dayi - dayk;
134 <  dazj = -dazi - dazk;
135 <
136 <  aF[0] = daxi;
137 <  aF[1] = dayi;
138 <  aF[2] = dazi;
139 <
140 <  bF[0] = daxj + daxk;
141 <  bF[1] = dayj + dayk;
142 <  bF[2] = dazj + dazk;
143 <
144 <  bTrq[0] = gy*dazk - gz*dayk;
145 <  bTrq[1] = gz*daxk - gx*dazk;
146 <  bTrq[2] = gx*dayk - gy*daxk;
147 <  
148 <  
149 <  c_p_a->addFrc( aF );
150 <  atomB->addFrc( bF );
151 <  atomB->addTrq( bTrq );
152 <  
153 <  return;
154 < }
155 <
156 < void GhostBend::setConstants( double the_c1, double the_c2, double the_c3,
157 <                                  double the_Th0 ){
158 <  c1 = the_c1;
159 <  c2 = the_c2;
160 <  c3 = the_c3;
161 <  theta0 = the_Th0;
162 < }
163 <
164 <
165 < double GhostBend::bend_force( double theta ){
166 <
167 <  double dt, dt2;
168 <  double force;
169 <
170 <  dt = ( theta - theta0 ) * M_PI / 180.0;
171 <  dt2 = dt * dt;
172 <
173 <  c_potential_E = ( c1 * dt2 ) + ( c2 * dt ) + c3;
174 <  force = -( ( 2.0 * c1 * dt ) + c2 );
175 <  return force;
176 < }
1 >
2 > #include <math.h>
3 > #include <iostream>
4 > #include <stdlib.h>
5 >
6 > #include "utils/simError.h"
7 > #include "primitives/SRI.hpp"
8 > #include "primitives/Atom.hpp"
9 >
10 >
11 >
12 > GhostBend::GhostBend( Atom &a, Atom &b ){
13 >  
14 >  c_p_a = &a;
15 >  
16 >  if( !b.isDirectional() ){
17 >    
18 >    // if atom b is not directional, then bad things will happen
19 >    
20 >    sprintf( painCave.errMsg,
21 >             " Ghost Bend error: Atom # %d of type \"%s\" is not "
22 >             "directional.\n",
23 >             b.getIndex(),
24 >             b.getType() );
25 >    painCave.isFatal = 1;
26 >    simError();
27 >  }    
28 >
29 >  atomB = ( DirectionalAtom* ) &b;
30 >  
31 >  c_potential_E = 0.0;
32 >
33 > }
34 >
35 >
36 > void GhostBend::calc_forces(){
37 >  
38 >  double dx,dy,dz,gx,gy,gz,dx2,dy2,dz2,gx2,gy2,gz2;
39 >  double rij2, rkj2, riji2, rkji2, dot, denom, cosang, angl;
40 >  
41 >  double sina2, sinai;
42 >
43 >  double comf2, comf3, comf4;
44 >  double dcsidx, dcsidy, dcsidz, dcskdx, dcskdy, dcskdz;
45 >  // double dcsjdx, dcsjdy, dcsjdz;
46 >  double dadxi, dadyi, dadzi;
47 >  double dadxk, dadyk, dadzk;//, dadxj, dadyj, dadzj;
48 >  double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj;
49 >  Vector3d u;
50 >  
51 >  double aR[3], bR[3];
52 >  double aF[3], bF[3], bTrq[3];
53 >
54 >  aR = c_p_a->getPos();
55 >  bR = atomB->getPos();
56 >  
57 >
58 >  dx = aR[0] - bR[0];
59 >  dy = aR[1] - bR[1];
60 >  dz = aR[2] - bR[2];
61 >
62 >  u = atomB->getU();
63 >
64 >  gx = u[0];
65 >  gy = u[1];
66 >  gz = u[2];
67 >  
68 >  dx2 = dx * dx;
69 >  dy2 = dy * dy;
70 >  dz2 = dz * dz;
71 >
72 >  gx2 = gx * gx;
73 >  gy2 = gy * gy;
74 >  gz2 = gz * gz;
75 >  
76 >  rij2 = dx2 + dy2 + dz2;
77 >  rkj2 = gx2 + gy2 + gz2;
78 >  
79 >  riji2 = 1.0 / rij2;
80 >  rkji2 = 1.0 / rkj2;
81 >
82 >  dot = dx * gx + dy * gy + dz * gz;
83 >  denom = sqrt((riji2 * rkji2));
84 >  cosang = dot * denom;
85 >
86 >  if(cosang > 1.0)cosang = 1.0;
87 >  if(cosang < -1.0) cosang = -1.0;
88 >
89 >  angl = acos(cosang);
90 >  angl = angl * 180.0 / M_PI;
91 >
92 >  sina2 = 1.0 - cosang*cosang;
93 >  if(fabs(sina2) < 1.0E-12 ) sina2 = 1.0E-12;
94 >  sinai = 1.0 / sqrt(sina2);
95 >
96 >  comf2 = cosang * riji2;
97 >  comf3 = cosang * rkji2;
98 >  comf4 = bend_force(angl);
99 >
100 >  dcsidx = gx*denom - comf2*dx;
101 >  dcsidy = gy*denom - comf2*dy;
102 >  dcsidz = gz*denom - comf2*dz;
103 >  
104 >  dcskdx = dx*denom - comf3*gx;
105 >  dcskdy = dy*denom - comf3*gy;
106 >  dcskdz = dz*denom - comf3*gz;
107 >  
108 > //   dcsjdx = -dcsidx - dcskdx;
109 > //   dcsjdy = -dcsidy - dcskdy;
110 > //   dcsjdz = -dcsidz - dcskdz;
111 >
112 >  dadxi = -sinai*dcsidx;
113 >  dadyi = -sinai*dcsidy;
114 >  dadzi = -sinai*dcsidz;
115 >
116 >  dadxk = -sinai*dcskdx;
117 >  dadyk = -sinai*dcskdy;
118 >  dadzk = -sinai*dcskdz;
119 >
120 > //   dadxj = -dadxi - dadxk;
121 > //   dadyj = -dadyi - dadyk;
122 > //   dadzj = -dadzi - dadzk;
123 >
124 >  daxi = comf4*dadxi;
125 >  dayi = comf4*dadyi;
126 >  dazi = comf4*dadzi;
127 >
128 >  daxk = comf4*dadxk;
129 >  dayk = comf4*dadyk;
130 >  dazk = comf4*dadzk;
131 >  
132 >  daxj = -daxi - daxk;
133 >  dayj = -dayi - dayk;
134 >  dazj = -dazi - dazk;
135 >
136 >  aF[0] = daxi;
137 >  aF[1] = dayi;
138 >  aF[2] = dazi;
139 >
140 >  bF[0] = daxj + daxk;
141 >  bF[1] = dayj + dayk;
142 >  bF[2] = dazj + dazk;
143 >
144 >  bTrq[0] = gy*dazk - gz*dayk;
145 >  bTrq[1] = gz*daxk - gx*dazk;
146 >  bTrq[2] = gx*dayk - gy*daxk;
147 >  
148 >  
149 >  c_p_a->addFrc( aF );
150 >  atomB->addFrc( bF );
151 >  atomB->addTrq( bTrq );
152 >  
153 >  return;
154 > }
155 >
156 > void GhostBend::setConstants( double the_c1, double the_c2, double the_c3,
157 >                                  double the_Th0 ){
158 >  c1 = the_c1;
159 >  c2 = the_c2;
160 >  c3 = the_c3;
161 >  theta0 = the_Th0;
162 > }
163 >
164 >
165 > double GhostBend::bend_force( double theta ){
166 >
167 >  double dt, dt2;
168 >  double force;
169 >
170 >  dt = ( theta - theta0 ) * M_PI / 180.0;
171 >  dt2 = dt * dt;
172 >
173 >  c_potential_E = ( c1 * dt2 ) + ( c2 * dt ) + c3;
174 >  force = -( ( 2.0 * c1 * dt ) + c2 );
175 >  return force;
176 > }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines