--- trunk/src/applications/staticProps/RippleOP.cpp 2008/01/23 21:21:50 1213 +++ trunk/src/applications/staticProps/RippleOP.cpp 2015/03/07 21:41:51 2071 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,16 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include "applications/staticProps/RippleOP.hpp" @@ -44,13 +45,15 @@ #include "io/DumpReader.hpp" #include "primitives/Molecule.hpp" #include "utils/NumericConstant.hpp" -namespace oopse { +#include "types/MultipoleAdapter.hpp" +namespace OpenMD { - - RippleOP::RippleOP(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2) - : StaticAnalyser(info, filename), - selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), evaluator2_(info), - seleMan1_(info), seleMan2_(info){ + + RippleOP::RippleOP(SimInfo* info, const std::string& filename, + const std::string& sele1, const std::string& sele2) + : StaticAnalyser(info, filename), + selectionScript1_(sele1), selectionScript2_(sele2), + seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info) { setOutputName(getPrefix(filename) + ".rp2"); @@ -62,7 +65,7 @@ namespace oopse { }else { sprintf( painCave.errMsg, "--sele1 must be static selection\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } @@ -72,7 +75,7 @@ namespace oopse { }else { sprintf( painCave.errMsg, "--sele2 must be static selection\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } @@ -80,7 +83,7 @@ namespace oopse { if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) { sprintf( painCave.errMsg, "The number of selected Stuntdoubles are not the same in --sele1 and sele2\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); @@ -122,46 +125,63 @@ namespace oopse { int nTail=0; RealType sumZ = 0.0; - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { //change the positions of atoms which belong to the rigidbodies - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { rb->updateAtoms(); } } - for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; sd3 = seleMan2_.nextSelected(i1)) { + for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; + sd3 = seleMan2_.nextSelected(i1)) { Vector3d pos1 = sd3->getPos(); if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(pos1); sd3->setPos(pos1); } - for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; sd3 = seleMan2_.nextSelected(i1)) { + for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL; + sd3 = seleMan2_.nextSelected(i1)) { Vector3d pos1 = sd3->getPos(); sumZ += pos1.z(); } RealType avgZ = sumZ / (RealType) nMolecules; - Mat3x3d orderTensorHeadUpper(0.0), orderTensorTail(0.0), orderTensorHeadLower(0.0); - // for (std::vector >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) { - for (j1 = seleMan1_.beginSelected(i1); j1 != NULL; j1 = seleMan1_.nextSelected(i1)) { + Mat3x3d orderTensorHeadUpper(0.0); + Mat3x3d orderTensorTail(0.0); + Mat3x3d orderTensorHeadLower(0.0); + for (j1 = seleMan1_.beginSelected(i1); j1 != NULL; + j1 = seleMan1_.nextSelected(i1)) { Vector3d pos = j1->getPos(); if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(pos); Vector3d vecHeadUpper; if (pos.z() >= avgZ){ - vecHeadUpper = j1->getElectroFrame().getColumn(2); + AtomType* atype1 = static_cast(j1)->getAtomType(); + MultipoleAdapter ma1 = MultipoleAdapter(atype1); + if (ma1.isDipole()) + vecHeadUpper = j1->getDipole(); + else + vecHeadUpper = j1->getA().transpose()*V3Z; nUpper++; } Vector3d vecHeadLower; if (pos.z() <= avgZ){ - vecHeadLower = j1->getElectroFrame().getColumn(2); + AtomType* atype1 = static_cast(j1)->getAtomType(); + MultipoleAdapter ma1 = MultipoleAdapter(atype1); + if (ma1.isDipole()) + vecHeadLower = j1->getDipole(); + else + vecHeadLower = j1->getA().transpose() * V3Z; nLower++; } orderTensorHeadUpper +=outProduct(vecHeadUpper, vecHeadUpper); orderTensorHeadLower +=outProduct(vecHeadLower, vecHeadLower); } - for (j2 = seleMan2_.beginSelected(i1); j2 != NULL; j2 = seleMan2_.nextSelected(i1)) { + for (j2 = seleMan2_.beginSelected(i1); j2 != NULL; + j2 = seleMan2_.nextSelected(i1)) { // The lab frame vector corresponding to the body-fixed // z-axis is simply the second column of A.transpose() // or, identically, the second row of A itself. @@ -180,12 +200,14 @@ namespace oopse { orderTensorTail -= (RealType)(1.0/3.0) * Mat3x3d::identity(); Vector3d eigenvaluesHeadUpper, eigenvaluesHeadLower, eigenvaluesTail; - Mat3x3d eigenvectorsHeadUpper, eigenvectorsHeadLower, eigenvectorsTail; - Mat3x3d::diagonalize(orderTensorHeadUpper, eigenvaluesHeadUpper, eigenvectorsHeadUpper); - Mat3x3d::diagonalize(orderTensorHeadLower, eigenvaluesHeadLower, eigenvectorsHeadLower); + Mat3x3d eigenvectorsHeadUpper, eigenvectorsHeadLower, eigenvectorsTail; + Mat3x3d::diagonalize(orderTensorHeadUpper, eigenvaluesHeadUpper, + eigenvectorsHeadUpper); + Mat3x3d::diagonalize(orderTensorHeadLower, eigenvaluesHeadLower, + eigenvectorsHeadLower); Mat3x3d::diagonalize(orderTensorTail, eigenvaluesTail, eigenvectorsTail); - int whichUpper, whichLower, whichTail; + int whichUpper(-1), whichLower(-1), whichTail(-1); RealType maxEvalUpper = 0.0; RealType maxEvalLower = 0.0; RealType maxEvalTail = 0.0; @@ -211,7 +233,7 @@ namespace oopse { } RealType p2Tail = 1.5 * maxEvalTail; - //the eigen vector is already normalized in SquareMatrix3::diagonalize + //the eigenvector is already normalized in SquareMatrix3::diagonalize Vector3d directorHeadUpper = eigenvectorsHeadUpper.getColumn(whichUpper); if (directorHeadUpper[0] < 0) { directorHeadUpper.negate();