ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Rattle.cpp
Revision: 1254
Committed: Wed Jun 9 16:16:33 2004 UTC (20 years, 1 month ago) by tim
File size: 5477 byte(s)
Log Message: 
1. adding some useful math classes(Mat3x3d, Vector3d, Quaternion, Euler3)
 these classes use anonymous union and struct to support
 double[3], double[3][3] and double[4]
2. adding roll constraint algorithm

File Contents

# Content
1 #include "Rattle.hpp"
2 #include <cmath>
3 #include "SimInfo.hpp"
4 #include "MatVec3.h"
5 ////////////////////////////////////////////////////////////////////////////////
6 //Implementation of DCShakeFunctor
7 ////////////////////////////////////////////////////////////////////////////////
8 int DCRattleAFunctor::operator()(ConstraintAtom* consAtom1, ConstraintAtom* consAtom2){
9 double posA[3];
10 double posB[3];
11 double oldPosA[3];
12 double oldPosB[3];
13 double velA[3];
14 double velB[3];
15 double pab[3];
16 double rab[3];
17 double rma, rmb;
18 double dx, dy, dz;
19 double rpab;
20 double rabsq, pabsq, rpabsq;
21 double diffsq;
22 double gab;
23 double dt;
24
25 dt = info->dt;
26
27 consAtom1->getPos(posA);
28 consAtom2->getPos(posB);
29
30
31 pab[0] = posA[0] - posB[0];
32 pab[1] = posA[1] - posB[1];
33 pab[2] = posA[2] - posB[2];
34
35 //periodic boundary condition
36
37 info->wrapVector(pab);
38
39 pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
40
41 rabsq = curPair->getBondLength2();
42 diffsq = rabsq - pabsq;
43
44 // the original rattle code from alan tidesley
45 if (fabs(diffsq) > (consTolerance * rabsq * 2)){
46
47 consAtom1->getOldPos(oldPosA);
48 consAtom2->getOldPos(oldPosB);
49
50 rab[0] = oldPosA[0] - oldPosB[0];
51 rab[1] = oldPosA[1] - oldPosB[1];
52 rab[2] = oldPosA[2] - oldPosB[2];
53
54 info->wrapVector(rab);
55
56 rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
57
58 rpabsq = rpab * rpab;
59
60
61 if (rpabsq < (rabsq * -diffsq)){
62 return consFail;
63 }
64
65 rma = 1.0 / consAtom1->getMass();
66 rmb = 1.0 / consAtom2->getMass();
67
68 gab = diffsq / (2.0 * (rma + rmb) * rpab);
69
70 dx = rab[0] * gab;
71 dy = rab[1] * gab;
72 dz = rab[2] * gab;
73
74 //set atom1's position
75 posA[0] += rma * dx;
76 posA[1] += rma * dy;
77 posA[2] += rma * dz;
78 consAtom1->setPos(posA);
79
80 //set atom2's position
81 posB[0] -= rmb * dx;
82 posB[1] -= rmb * dy;
83 posB[2] -= rmb * dz;
84 consAtom2->setPos(posB);
85
86 dx = dx / dt;
87 dy = dy / dt;
88 dz = dz / dt;
89
90 //set atom1's velocity
91 consAtom1->getVel(velA);
92 velA[0] += rma * dx;
93 velA[1] += rma * dy;
94 velA[2] += rma * dz;
95 consAtom1->setVel(velA);
96
97 //set atom2's velocity
98 consAtom2->getVel(velB);
99 velB[0] -= rmb * dx;
100 velB[1] -= rmb * dy;
101 velB[2] -= rmb * dz;
102 consAtom2->setVel(velB);
103
104 return consSuccess;
105 }
106 else
107 return consAlready;
108
109 }
110
111
112 int DCRattleAFunctor::operator()(ConstraintAtom* consAtom,ConstraintRigidBody* consRB){
113 return consElemHandlerFail;
114 }
115
116 int DCRattleAFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
117 //calculate lamda
118
119 //integrate
120
121 //
122 return consElemHandlerFail;
123 }
124 ////////////////////////////////////////////////////////////////////////////////
125 //Implementation of JCShakeFunctor
126 ////////////////////////////////////////////////////////////////////////////////
127 int JCRattleAFunctor::operator()(ConstraintAtom* consAtom1, ConstraintAtom* consAtom2){
128 return consElemHandlerFail;
129 }
130
131 int JCRattleAFunctor::operator()(ConstraintAtom* consAtom,ConstraintRigidBody* consRB){
132 return consElemHandlerFail;
133
134 }
135
136 int JCRattleAFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
137 //calculate lamda,
138 //integrate
139 //
140 return consElemHandlerFail;
141 }
142
143
144 ////////////////////////////////////////////////////////////////////////////////
145 //Implementation of DCRattleBFunctor
146 ////////////////////////////////////////////////////////////////////////////////
147 int DCRattleBFunctor::operator()(ConstraintAtom* consAtom1, ConstraintAtom* consAtom2){
148 double posA[3], posB[3];
149 double velA[3], velB[3];
150 double vxab, vyab, vzab;
151 double rab[3];
152 double rma, rmb;
153 double dx, dy, dz;
154 double rvab;
155 double gab;
156
157 consAtom1->getVel(velA);
158 consAtom2->getVel(velB);
159
160 vxab = velA[0] - velB[0];
161 vyab = velA[1] - velB[1];
162 vzab = velA[2] - velB[2];
163
164 consAtom1->getPos(posA);
165 consAtom2->getPos(posB);
166
167 rab[0] = posA[0] - posB[0];
168 rab[1] = posA[1] - posB[1];
169 rab[2] = posA[2] - posB[2];
170
171 info->wrapVector(rab);
172
173 rma = 1.0 / consAtom1->getMass();
174 rmb = 1.0 / consAtom2->getMass();
175
176 rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
177
178 gab = -rvab / ((rma + rmb) * curPair->getBondLength2());
179
180 if (fabs(gab) > consTolerance){
181 dx = rab[0] * gab;
182 dy = rab[1] * gab;
183 dz = rab[2] * gab;
184
185 velA[0] += rma * dx;
186 velA[1] += rma * dy;
187 velA[2] += rma * dz;
188
189 consAtom1->setVel(velA);
190
191 velB[0] -= rmb * dx;
192 velB[1] -= rmb * dy;
193 velB[2] -= rmb * dz;
194
195 consAtom2->setVel(velB);
196 }
197 return consSuccess;
198 }
199
200
201
202 int DCRattleBFunctor::operator()(ConstraintAtom* consAtom,ConstraintRigidBody* consRB){
203 return consElemHandlerFail;
204 }
205
206 int DCRattleBFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
207 return consElemHandlerFail;
208 }
209
210 ////////////////////////////////////////////////////////////////////////////////
211 //Implementation of JCRattleBFunctor
212 ////////////////////////////////////////////////////////////////////////////////
213 int JCRattleBFunctor::operator()(ConstraintAtom* consAtom1, ConstraintAtom* consAtom2){
214 return consElemHandlerFail;
215 }
216
217
218
219 int JCRattleBFunctor::operator()(ConstraintAtom* consAtom,ConstraintRigidBody* consRB){
220 return consElemHandlerFail;
221 }
222
223 int JCRattleBFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
224 return consElemHandlerFail;
225 }

Properties

Name Value
svn:executable *