| 1 | < | /* | 
| 1 | > | /* | 
| 2 |  | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 3 |  | * | 
| 4 |  | * The University of Notre Dame grants you ("Licensee") a | 
| 43 |  | #include "utils/simError.h" | 
| 44 |  | namespace oopse { | 
| 45 |  |  | 
| 46 | < | DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType) | 
| 47 | < | : Atom(dAtomType){ | 
| 48 | < | objType_= otDAtom; | 
| 49 | < | if (dAtomType->isMultipole()) { | 
| 46 | > | DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType) | 
| 47 | > | : Atom(dAtomType){ | 
| 48 | > | objType_= otDAtom; | 
| 49 | > | if (dAtomType->isMultipole()) { | 
| 50 |  | electroBodyFrame_ = dAtomType->getElectroBodyFrame(); | 
| 51 | < | } | 
| 51 | > | } | 
| 52 |  |  | 
| 53 | < | //check if one of the diagonal inertia tensor of this directional atom  is zero | 
| 54 | < | int nLinearAxis = 0; | 
| 55 | < | Mat3x3d inertiaTensor = getI(); | 
| 56 | < | for (int i = 0; i < 3; i++) { | 
| 53 | > | //check if one of the diagonal inertia tensor of this directional atom  is zero | 
| 54 | > | int nLinearAxis = 0; | 
| 55 | > | Mat3x3d inertiaTensor = getI(); | 
| 56 | > | for (int i = 0; i < 3; i++) { | 
| 57 |  | if (fabs(inertiaTensor(i, i)) < oopse::epsilon) { | 
| 58 | < | linear_ = true; | 
| 59 | < | linearAxis_ = i; | 
| 60 | < | ++ nLinearAxis; | 
| 58 | > | linear_ = true; | 
| 59 | > | linearAxis_ = i; | 
| 60 | > | ++ nLinearAxis; | 
| 61 |  | } | 
| 62 | < | } | 
| 62 | > | } | 
| 63 |  |  | 
| 64 | < | if (nLinearAxis > 1) { | 
| 64 | > | if (nLinearAxis > 1) { | 
| 65 |  | sprintf( painCave.errMsg, | 
| 66 | < | "Directional Atom error.\n" | 
| 67 | < | "\tOOPSE found more than one axis in this directional atom with a vanishing \n" | 
| 68 | < | "\tmoment of inertia."); | 
| 69 | < | painCave.isFatal = 1; | 
| 66 | > | "Directional Atom warning.\n" | 
| 67 | > | "\tOOPSE found more than one axis in this directional atom with a vanishing \n" | 
| 68 | > | "\tmoment of inertia."); | 
| 69 | > | painCave.isFatal = 0; | 
| 70 |  | simError(); | 
| 71 | < | } | 
| 71 | > | } | 
| 72 |  |  | 
| 73 | < | } | 
| 73 | > | } | 
| 74 |  |  | 
| 75 | < | Mat3x3d DirectionalAtom::getI() { | 
| 75 | > | Mat3x3d DirectionalAtom::getI() { | 
| 76 |  | return static_cast<DirectionalAtomType*>(getAtomType())->getI(); | 
| 77 | < | } | 
| 77 | > | } | 
| 78 |  |  | 
| 79 | < | void DirectionalAtom::setPrevA(const RotMat3x3d& a) { | 
| 79 | > | void DirectionalAtom::setPrevA(const RotMat3x3d& a) { | 
| 80 |  | ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a; | 
| 81 |  | if (atomType_->isMultipole()) { | 
| 82 | < | ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; | 
| 82 | > | ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; | 
| 83 |  | } | 
| 84 | < | } | 
| 84 | > | } | 
| 85 |  |  | 
| 86 |  |  | 
| 87 | < | void DirectionalAtom::setA(const RotMat3x3d& a) { | 
| 87 | > | void DirectionalAtom::setA(const RotMat3x3d& a) { | 
| 88 |  | ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a; | 
| 89 |  |  | 
| 90 |  | if (atomType_->isMultipole()) { | 
| 91 | < | ((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; | 
| 91 | > | ((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; | 
| 92 |  | } | 
| 93 | < | } | 
| 93 | > | } | 
| 94 |  |  | 
| 95 | < | void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) { | 
| 95 | > | void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) { | 
| 96 |  | ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a; | 
| 97 |  |  | 
| 98 |  | if (atomType_->isMultipole()) { | 
| 99 | < | ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; | 
| 99 | > | ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; | 
| 100 |  | } | 
| 101 | < | } | 
| 101 | > | } | 
| 102 |  |  | 
| 103 | < | void DirectionalAtom::rotateBy(const RotMat3x3d& m) { | 
| 103 | > | void DirectionalAtom::rotateBy(const RotMat3x3d& m) { | 
| 104 |  | setA(m *getA()); | 
| 105 | < | } | 
| 105 | > | } | 
| 106 |  |  | 
| 107 | < | std::vector<double> DirectionalAtom::getGrad() { | 
| 107 | > | std::vector<double> DirectionalAtom::getGrad() { | 
| 108 |  | std::vector<double> grad(6, 0.0); | 
| 109 |  | Vector3d force; | 
| 110 |  | Vector3d torque; | 
| 144 |  |  | 
| 145 |  | //gradient is equal to -force | 
| 146 |  | for (int j = 0 ; j<3; j++) | 
| 147 | < | grad[j] = -force[j]; | 
| 147 | > | grad[j] = -force[j]; | 
| 148 |  |  | 
| 149 |  | for (int j = 0; j < 3; j++ ) { | 
| 150 |  |  | 
| 151 | < | grad[3] += torque[j]*ephi[j]; | 
| 152 | < | grad[4] += torque[j]*etheta[j]; | 
| 153 | < | grad[5] += torque[j]*epsi[j]; | 
| 151 | > | grad[3] -= torque[j]*ephi[j]; | 
| 152 | > | grad[4] -= torque[j]*etheta[j]; | 
| 153 | > | grad[5] -= torque[j]*epsi[j]; | 
| 154 |  |  | 
| 155 |  | } | 
| 156 |  |  | 
| 157 |  | return grad; | 
| 158 | < | } | 
| 158 | > | } | 
| 159 |  |  | 
| 160 | < | void DirectionalAtom::accept(BaseVisitor* v) { | 
| 160 | > | void DirectionalAtom::accept(BaseVisitor* v) { | 
| 161 |  | v->visit(this); | 
| 162 | < | } | 
| 162 | > | } | 
| 163 |  |  | 
| 164 |  | } | 
| 165 |  |  |