| 160 |
|
|
| 161 |
|
currentSnapshot_->wrapVector(pab); |
| 162 |
|
|
| 163 |
< |
double pabsq = pab.lengthSquare(); |
| 163 |
> |
RealType pabsq = pab.lengthSquare(); |
| 164 |
|
|
| 165 |
< |
double rabsq = consPair->getConsDistSquare(); |
| 166 |
< |
double diffsq = rabsq - pabsq; |
| 165 |
> |
RealType rabsq = consPair->getConsDistSquare(); |
| 166 |
> |
RealType diffsq = rabsq - pabsq; |
| 167 |
|
|
| 168 |
|
// the original rattle code from alan tidesley |
| 169 |
|
if (fabs(diffsq) > (consTolerance_ * rabsq * 2)){ |
| 175 |
|
|
| 176 |
|
currentSnapshot_->wrapVector(rab); |
| 177 |
|
|
| 178 |
< |
double rpab = dot(rab, pab); |
| 179 |
< |
double rpabsq = rpab * rpab; |
| 178 |
> |
RealType rpab = dot(rab, pab); |
| 179 |
> |
RealType rpabsq = rpab * rpab; |
| 180 |
|
|
| 181 |
|
if (rpabsq < (rabsq * -diffsq)){ |
| 182 |
|
return consFail; |
| 183 |
|
} |
| 184 |
|
|
| 185 |
< |
double rma = 1.0 / consElem1->getMass(); |
| 186 |
< |
double rmb = 1.0 / consElem2->getMass(); |
| 185 |
> |
RealType rma = 1.0 / consElem1->getMass(); |
| 186 |
> |
RealType rmb = 1.0 / consElem2->getMass(); |
| 187 |
|
|
| 188 |
< |
double gab = diffsq / (2.0 * (rma + rmb) * rpab); |
| 188 |
> |
RealType gab = diffsq / (2.0 * (rma + rmb) * rpab); |
| 189 |
|
|
| 190 |
|
Vector3d delta = rab * gab; |
| 191 |
|
|
| 234 |
|
|
| 235 |
|
currentSnapshot_->wrapVector(rab); |
| 236 |
|
|
| 237 |
< |
double rma = 1.0 / consElem1->getMass(); |
| 238 |
< |
double rmb = 1.0 / consElem2->getMass(); |
| 237 |
> |
RealType rma = 1.0 / consElem1->getMass(); |
| 238 |
> |
RealType rmb = 1.0 / consElem2->getMass(); |
| 239 |
|
|
| 240 |
< |
double rvab = dot(rab, dv); |
| 240 |
> |
RealType rvab = dot(rab, dv); |
| 241 |
|
|
| 242 |
< |
double gab = -rvab / ((rma + rmb) * consPair->getConsDistSquare()); |
| 242 |
> |
RealType gab = -rvab / ((rma + rmb) * consPair->getConsDistSquare()); |
| 243 |
|
|
| 244 |
|
if (fabs(gab) > consTolerance_){ |
| 245 |
|
Vector3d delta = rab * gab; |
