| 35 |
|
* |
| 36 |
|
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
| 37 |
|
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
| 38 |
< |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
| 38 |
> |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
| 39 |
|
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
| 40 |
|
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 45 |
|
#include "io/DumpReader.hpp" |
| 46 |
|
#include "primitives/Molecule.hpp" |
| 47 |
|
#include "utils/NumericConstant.hpp" |
| 48 |
+ |
#include "types/MultipoleAdapter.hpp" |
| 49 |
|
namespace OpenMD { |
| 50 |
|
|
| 51 |
< |
|
| 52 |
< |
RippleOP::RippleOP(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2) |
| 53 |
< |
: StaticAnalyser(info, filename), |
| 54 |
< |
selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), evaluator2_(info), |
| 55 |
< |
seleMan1_(info), seleMan2_(info){ |
| 51 |
> |
|
| 52 |
> |
RippleOP::RippleOP(SimInfo* info, const std::string& filename, |
| 53 |
> |
const std::string& sele1, const std::string& sele2) |
| 54 |
> |
: StaticAnalyser(info, filename), |
| 55 |
> |
selectionScript1_(sele1), selectionScript2_(sele2), |
| 56 |
> |
seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info) { |
| 57 |
|
|
| 58 |
|
setOutputName(getPrefix(filename) + ".rp2"); |
| 59 |
|
|
| 125 |
|
int nTail=0; |
| 126 |
|
RealType sumZ = 0.0; |
| 127 |
|
|
| 128 |
< |
for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
| 128 |
> |
for (mol = info_->beginMolecule(mi); mol != NULL; |
| 129 |
> |
mol = info_->nextMolecule(mi)) { |
| 130 |
|
//change the positions of atoms which belong to the rigidbodies |
| 131 |
< |
for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
| 131 |
> |
for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
| 132 |
> |
rb = mol->nextRigidBody(rbIter)) { |
| 133 |
|
rb->updateAtoms(); |
| 134 |
|
} |
| 135 |
|
} |
| 136 |
|
|
| 137 |
< |
for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; sd3 = seleMan2_.nextSelected(i1)) { |
| 137 |
> |
for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; |
| 138 |
> |
sd3 = seleMan2_.nextSelected(i1)) { |
| 139 |
|
Vector3d pos1 = sd3->getPos(); |
| 140 |
|
if (usePeriodicBoundaryConditions_) |
| 141 |
|
currentSnapshot_->wrapVector(pos1); |
| 142 |
|
sd3->setPos(pos1); |
| 143 |
|
} |
| 144 |
|
|
| 145 |
< |
for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; sd3 = seleMan2_.nextSelected(i1)) { |
| 145 |
> |
for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; |
| 146 |
> |
sd3 = seleMan2_.nextSelected(i1)) { |
| 147 |
|
Vector3d pos1 = sd3->getPos(); |
| 148 |
|
sumZ += pos1.z(); |
| 149 |
|
} |
| 150 |
|
RealType avgZ = sumZ / (RealType) nMolecules; |
| 151 |
|
|
| 152 |
< |
Mat3x3d orderTensorHeadUpper(0.0), orderTensorTail(0.0), orderTensorHeadLower(0.0); |
| 153 |
< |
// for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) { |
| 154 |
< |
for (j1 = seleMan1_.beginSelected(i1); j1 != NULL; j1 = seleMan1_.nextSelected(i1)) { |
| 152 |
> |
Mat3x3d orderTensorHeadUpper(0.0); |
| 153 |
> |
Mat3x3d orderTensorTail(0.0); |
| 154 |
> |
Mat3x3d orderTensorHeadLower(0.0); |
| 155 |
> |
for (j1 = seleMan1_.beginSelected(i1); j1 != NULL; |
| 156 |
> |
j1 = seleMan1_.nextSelected(i1)) { |
| 157 |
|
Vector3d pos = j1->getPos(); |
| 158 |
|
if (usePeriodicBoundaryConditions_) |
| 159 |
|
currentSnapshot_->wrapVector(pos); |
| 160 |
|
Vector3d vecHeadUpper; |
| 161 |
|
if (pos.z() >= avgZ){ |
| 162 |
< |
vecHeadUpper = j1->getElectroFrame().getColumn(2); |
| 162 |
> |
AtomType* atype1 = static_cast<Atom*>(j1)->getAtomType(); |
| 163 |
> |
MultipoleAdapter ma1 = MultipoleAdapter(atype1); |
| 164 |
> |
if (ma1.isDipole()) |
| 165 |
> |
vecHeadUpper = j1->getDipole(); |
| 166 |
> |
else |
| 167 |
> |
vecHeadUpper = j1->getA().transpose()*V3Z; |
| 168 |
|
nUpper++; |
| 169 |
|
} |
| 170 |
|
Vector3d vecHeadLower; |
| 171 |
|
if (pos.z() <= avgZ){ |
| 172 |
< |
vecHeadLower = j1->getElectroFrame().getColumn(2); |
| 172 |
> |
AtomType* atype1 = static_cast<Atom*>(j1)->getAtomType(); |
| 173 |
> |
MultipoleAdapter ma1 = MultipoleAdapter(atype1); |
| 174 |
> |
if (ma1.isDipole()) |
| 175 |
> |
vecHeadLower = j1->getDipole(); |
| 176 |
> |
else |
| 177 |
> |
vecHeadLower = j1->getA().transpose() * V3Z; |
| 178 |
|
nLower++; |
| 179 |
|
} |
| 180 |
|
orderTensorHeadUpper +=outProduct(vecHeadUpper, vecHeadUpper); |
| 181 |
|
orderTensorHeadLower +=outProduct(vecHeadLower, vecHeadLower); |
| 182 |
|
} |
| 183 |
< |
for (j2 = seleMan2_.beginSelected(i1); j2 != NULL; j2 = seleMan2_.nextSelected(i1)) { |
| 183 |
> |
for (j2 = seleMan2_.beginSelected(i1); j2 != NULL; |
| 184 |
> |
j2 = seleMan2_.nextSelected(i1)) { |
| 185 |
|
// The lab frame vector corresponding to the body-fixed |
| 186 |
|
// z-axis is simply the second column of A.transpose() |
| 187 |
|
// or, identically, the second row of A itself. |
| 200 |
|
orderTensorTail -= (RealType)(1.0/3.0) * Mat3x3d::identity(); |
| 201 |
|
|
| 202 |
|
Vector3d eigenvaluesHeadUpper, eigenvaluesHeadLower, eigenvaluesTail; |
| 203 |
< |
Mat3x3d eigenvectorsHeadUpper, eigenvectorsHeadLower, eigenvectorsTail; |
| 204 |
< |
Mat3x3d::diagonalize(orderTensorHeadUpper, eigenvaluesHeadUpper, eigenvectorsHeadUpper); |
| 205 |
< |
Mat3x3d::diagonalize(orderTensorHeadLower, eigenvaluesHeadLower, eigenvectorsHeadLower); |
| 203 |
> |
Mat3x3d eigenvectorsHeadUpper, eigenvectorsHeadLower, eigenvectorsTail; |
| 204 |
> |
Mat3x3d::diagonalize(orderTensorHeadUpper, eigenvaluesHeadUpper, |
| 205 |
> |
eigenvectorsHeadUpper); |
| 206 |
> |
Mat3x3d::diagonalize(orderTensorHeadLower, eigenvaluesHeadLower, |
| 207 |
> |
eigenvectorsHeadLower); |
| 208 |
|
Mat3x3d::diagonalize(orderTensorTail, eigenvaluesTail, eigenvectorsTail); |
| 209 |
|
|
| 210 |
< |
int whichUpper, whichLower, whichTail; |
| 210 |
> |
int whichUpper(-1), whichLower(-1), whichTail(-1); |
| 211 |
|
RealType maxEvalUpper = 0.0; |
| 212 |
|
RealType maxEvalLower = 0.0; |
| 213 |
|
RealType maxEvalTail = 0.0; |
| 233 |
|
} |
| 234 |
|
RealType p2Tail = 1.5 * maxEvalTail; |
| 235 |
|
|
| 236 |
< |
//the eigen vector is already normalized in SquareMatrix3::diagonalize |
| 236 |
> |
//the eigenvector is already normalized in SquareMatrix3::diagonalize |
| 237 |
|
Vector3d directorHeadUpper = eigenvectorsHeadUpper.getColumn(whichUpper); |
| 238 |
|
if (directorHeadUpper[0] < 0) { |
| 239 |
|
directorHeadUpper.negate(); |
