--- trunk/src/applications/dump2Xyz/Dump2XYZ.cpp 2005/03/09 17:30:29 413 +++ trunk/src/applications/dump2Xyz/Dump2XYZ.cpp 2010/06/25 17:40:24 1456 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -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,15 @@ * 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, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). */ #include @@ -47,9 +47,11 @@ #include "brains/Register.hpp" #include "brains/SimCreator.hpp" #include "brains/SimInfo.hpp" +#include "brains/ForceManager.hpp" #include "io/DumpReader.hpp" #include "utils/simError.h" #include "visitors/AtomVisitor.hpp" +#include "visitors/ReplacementVisitor.hpp" #include "visitors/CompositeVisitor.hpp" #include "visitors/RigidBodyVisitor.hpp" #include "visitors/OtherVisitor.hpp" @@ -57,8 +59,9 @@ #include "selection/SelectionEvaluator.hpp" #include "selection/SelectionManager.hpp" #include "visitors/LipidTransVisitor.hpp" +#include "visitors/AtomNameVisitor.hpp" -using namespace oopse; +using namespace OpenMD; int main(int argc, char* argv[]){ @@ -67,8 +70,11 @@ int main(int argc, char* argv[]){ gengetopt_args_info args_info; std::string dumpFileName; - std::string mdFileName; std::string xyzFileName; + bool printVel; + bool printFrc; + bool printVec; + bool printChrg; //parse the command line option if (cmdline_parser (argc, argv, &args_info) != 0) { @@ -83,9 +89,6 @@ int main(int argc, char* argv[]){ exit(1); } - mdFileName = dumpFileName; - mdFileName = mdFileName.substr(0, mdFileName.rfind(".")) + ".md"; - if (args_info.output_given){ xyzFileName = args_info.output_arg; } else { @@ -95,34 +98,39 @@ int main(int argc, char* argv[]){ //parse md file and set up the system SimCreator creator; - SimInfo* info = creator.createSim(mdFileName, false); + SimInfo* info = creator.createSim(dumpFileName, false); + ForceManager* forceMan = new ForceManager(info); - - - //creat visitor list + //create visitor list CompositeVisitor* compositeVisitor = new CompositeVisitor(); - - //creat RigidBody Visitor + + //create RigidBody Visitor if(args_info.rigidbody_flag){ RBCOMVisitor* rbCOMVisitor = new RBCOMVisitor(info); compositeVisitor->addVisitor(rbCOMVisitor, 900); } - //creat SSD atom visitor + //create SSD atom visitor SSDAtomVisitor* ssdVisitor = new SSDAtomVisitor(info); compositeVisitor->addVisitor(ssdVisitor, 800); - LinearAtomVisitor* linearVisitor = new LinearAtomVisitor(info); - compositeVisitor->addVisitor(linearVisitor, 750); - - //creat default atom visitor + //create default atom visitor DefaultAtomVisitor* defaultAtomVisitor = new DefaultAtomVisitor(info); compositeVisitor->addVisitor(defaultAtomVisitor, 700); - //creat waterType visitor - if(args_info.watertype_flag){ - WaterTypeVisitor* waterTypeVisitor = new WaterTypeVisitor; - compositeVisitor->addVisitor(waterTypeVisitor, 600); + // if we gave the -w option, we want to skip the waters: + if (!args_info.water_given) { + //create waterType visitor + if(args_info.watertype_flag){ + WaterTypeVisitor* waterTypeVisitor = new WaterTypeVisitor; + compositeVisitor->addVisitor(waterTypeVisitor, 600); + } + } + + if (args_info.basetype_flag) { + AtomNameVisitor* atomNameVisitor = new AtomNameVisitor(info); + compositeVisitor->addVisitor(atomNameVisitor, 550); + std::cout << compositeVisitor->toString(); } //create ZconsVisitor @@ -137,74 +145,125 @@ int main(int argc, char* argv[]){ } } - //creat wrapping visitor + //create wrapping visitor - if(args_info.periodicBox_flag){ - WrappingVisitor* wrappingVisitor = new WrappingVisitor(info); - compositeVisitor->addVisitor(wrappingVisitor, 400); - } + //if(args_info.periodicBox_flag){ + // WrappingVisitor* wrappingVisitor = new WrappingVisitor(info); + // compositeVisitor->addVisitor(wrappingVisitor, 400); + //} - //creat replicate visitor - if(args_info.repeatX_given > 0 || args_info.repeatY_given > 0 ||args_info.repeatY_given > 0){ - Vector3i replicateOpt(args_info.repeatX_arg, args_info.repeatY_arg, args_info.repeatZ_arg); - ReplicateVisitor* replicateVisitor = new ReplicateVisitor(info, replicateOpt); + //create replicate visitor + if(args_info.repeatX_given > 0 || + args_info.repeatY_given > 0 || + args_info.repeatY_given > 0) { + Vector3i replicateOpt(args_info.repeatX_arg, + args_info.repeatY_arg, + args_info.repeatZ_arg); + ReplicateVisitor* replicateVisitor = new ReplicateVisitor(info, + replicateOpt); compositeVisitor->addVisitor(replicateVisitor, 300); } //create rotation visitor if (args_info.refsele_given&& args_info.originsele_given) { - compositeVisitor->addVisitor(new LipidTransVisitor(info, args_info.originsele_arg, args_info.refsele_arg), 250); + compositeVisitor->addVisitor(new LipidTransVisitor(info, + args_info.originsele_arg, + args_info.refsele_arg), + 250); } else if (args_info.refsele_given || args_info.originsele_given) { - std::cerr << "Both of --refsele and --originsele should appear by pair" << std::endl; + std::cerr << "Both of --refsele and --originsele should appear by pair" + << std::endl; exit(1); } - - //creat xyzVisitor + + //create xyzVisitor XYZVisitor* xyzVisitor; + if (args_info.selection_given) { xyzVisitor = new XYZVisitor(info, args_info.selection_arg); } else { xyzVisitor = new XYZVisitor(info); } - compositeVisitor->addVisitor(xyzVisitor, 200); + + if(args_info.velocities_flag){ + printVel = true; + xyzVisitor->doVelocities(printVel); + } + if(args_info.forces_flag){ + printFrc = true; + xyzVisitor->doForces(printFrc); + } + if(args_info.vectors_flag){ + printVec = true; + xyzVisitor->doVectors(printVec); + } + if(args_info.charges_flag){ + printChrg = true; + xyzVisitor->doCharges(printChrg); + } - std::cout << compositeVisitor->toString(); + compositeVisitor->addVisitor(xyzVisitor, 200); - //creat prepareVisitor + //create prepareVisitor PrepareVisitor* prepareVisitor = new PrepareVisitor(); //open dump file DumpReader* dumpReader = new DumpReader(info, dumpFileName); int nframes = dumpReader->getNFrames(); + std::ofstream xyzStream(xyzFileName.c_str()); - std::ofstream xyzStream; - xyzStream .open(xyzFileName.c_str()); - - SimInfo::MoleculeIterator miter; Molecule::IntegrableObjectIterator iiter; Molecule::RigidBodyIterator rbIter; Molecule* mol; StuntDouble* integrableObject; RigidBody* rb; + Vector3d molCom; + Vector3d newMolCom; + Vector3d displacement; + Mat3x3d hmat; + Snapshot* currentSnapshot; for (int i = 0; i < nframes; i += args_info.frame_arg){ dumpReader->readFrame(i); + if (printFrc) forceMan->calcForces(true, false); + + //wrapping the molecule + if(args_info.periodicBox_flag) { + currentSnapshot = info->getSnapshotManager()->getCurrentSnapshot(); + for (mol = info->beginMolecule(miter); mol != NULL; + mol = info->nextMolecule(miter)) { + molCom = mol->getCom(); + newMolCom = molCom; + currentSnapshot->wrapVector(newMolCom); + displacement = newMolCom - molCom; + for (integrableObject = mol->beginIntegrableObject(iiter); + integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(iiter)) { + integrableObject->setPos(integrableObject->getPos() + displacement); + } + } + } //update atoms of rigidbody - for (mol = info->beginMolecule(miter); mol != NULL; mol = info->nextMolecule(miter)) { + for (mol = info->beginMolecule(miter); mol != NULL; + mol = info->nextMolecule(miter)) { //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(); + if (printVel) rb->updateAtomVel(); } } //prepare visit - for (mol = info->beginMolecule(miter); mol != NULL; mol = info->nextMolecule(miter)) { - for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL; + for (mol = info->beginMolecule(miter); mol != NULL; + mol = info->nextMolecule(miter)) { + for (integrableObject = mol->beginIntegrableObject(iiter); + integrableObject != NULL; integrableObject = mol->nextIntegrableObject(iiter)) { integrableObject->accept(prepareVisitor); } @@ -215,8 +274,10 @@ int main(int argc, char* argv[]){ //visit stuntdouble - for (mol = info->beginMolecule(miter); mol != NULL; mol = info->nextMolecule(miter)) { - for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL; + for (mol = info->beginMolecule(miter); mol != NULL; + mol = info->nextMolecule(miter)) { + for (integrableObject = mol->beginIntegrableObject(iiter); + integrableObject != NULL; integrableObject = mol->nextIntegrableObject(iiter)) { integrableObject->accept(compositeVisitor); } @@ -226,11 +287,9 @@ int main(int argc, char* argv[]){ xyzVisitor->clear(); }//end for (int i = 0; i < nframes; i += args_info.frame_arg) - + xyzStream.close(); - - + delete prepareVisitor; delete compositeVisitor; delete info; - }