ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/dump2Xyz/Dump2XYZ.cpp
(Generate patch)

Comparing trunk/src/applications/dump2Xyz/Dump2XYZ.cpp (file contents):
Revision 250 by tim, Wed Jan 12 23:21:42 2005 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

# Line 1 | Line 1
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
# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
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, 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   */
42  
43   #include <iostream>
# Line 47 | Line 48
48   #include "brains/Register.hpp"
49   #include "brains/SimCreator.hpp"
50   #include "brains/SimInfo.hpp"
51 + #include "brains/ForceManager.hpp"
52   #include "io/DumpReader.hpp"
53   #include "utils/simError.h"
54   #include "visitors/AtomVisitor.hpp"
55 + #include "visitors/ReplacementVisitor.hpp"
56   #include "visitors/CompositeVisitor.hpp"
57   #include "visitors/RigidBodyVisitor.hpp"
58   #include "visitors/OtherVisitor.hpp"
59   #include "visitors/ZconsVisitor.hpp"
60 < using namespace oopse;
60 > #include "selection/SelectionEvaluator.hpp"
61 > #include "selection/SelectionManager.hpp"
62 > #include "visitors/LipidTransVisitor.hpp"
63 > #include "visitors/AtomNameVisitor.hpp"
64  
65 + using namespace OpenMD;
66 +
67 + using namespace std;
68   int main(int argc, char* argv[]){
69    
61  //register force fields
62  registerForceFields();
63  
70    gengetopt_args_info args_info;
71 <  std::string dumpFileName;
72 <  std::string mdFileName;
73 <  std::string xyzFileName;
71 >  string dumpFileName;
72 >  string xyzFileName;
73 >
74 >  bool printVel(false);
75 >  bool printFrc(false);
76 >  bool printVec(false);
77 >  bool printChrg(false);
78 >  bool printField(false);
79    
80    //parse the command line option
81    if (cmdline_parser (argc, argv, &args_info) != 0) {
# Line 75 | Line 86 | int main(int argc, char* argv[]){
86    if (args_info.input_given){
87      dumpFileName = args_info.input_arg;
88    } else {
89 <    std::cerr << "Does not have input file name" << std::endl;
89 >    cerr << "Does not have input file name" << endl;
90      exit(1);
91    }
92    
82  mdFileName = dumpFileName;
83  mdFileName = mdFileName.substr(0, mdFileName.rfind(".")) + ".md";
84  
93    if (args_info.output_given){
94      xyzFileName = args_info.output_arg;
95    } else {
# Line 91 | Line 99 | int main(int argc, char* argv[]){
99    
100    //parse md file and set up the system
101    SimCreator creator;
102 <  SimInfo* info = creator.createSim(mdFileName);
102 >  SimInfo* info = creator.createSim(dumpFileName, false);
103 >  ForceManager* forceMan = new ForceManager(info);
104    
105 <  
97 <  
98 <  //creat visitor list
105 >  //create visitor list
106    CompositeVisitor* compositeVisitor = new CompositeVisitor();
107    
108 <  //creat ignore visitor
102 <  if(args_info.ignore_given ||args_info.water_flag){
103 <    
104 <    IgnoreVisitor* ignoreVisitor = new IgnoreVisitor();
105 <    
106 <    for(int i = 0; i < args_info.ignore_given; i++)
107 <      ignoreVisitor->addIgnoreType(args_info.ignore_arg[i]);
108 <    
109 <    //ignore water
110 <    if(args_info.water_flag){
111 <      ignoreVisitor->addIgnoreType("SSD");
112 <      ignoreVisitor->addIgnoreType("SSD1");
113 <      ignoreVisitor->addIgnoreType("SSD_E");
114 <      ignoreVisitor->addIgnoreType("SSD_RF");
115 <      ignoreVisitor->addIgnoreType("TIP3P_RB_0");
116 <      ignoreVisitor->addIgnoreType("TIP4P_RB_0");
117 <      ignoreVisitor->addIgnoreType("TIP5P_RB_0");
118 <      ignoreVisitor->addIgnoreType("SPCE_RB_0");      
119 <      ignoreVisitor->addIgnoreType("DPD_RB_0");
120 <    }
121 <    
122 <    compositeVisitor->addVisitor(ignoreVisitor, 1000);
123 <  }
124 <  
125 <  //creat RigidBody Visitor
108 >  //create RigidBody Visitor
109    if(args_info.rigidbody_flag){
110      RBCOMVisitor* rbCOMVisitor = new RBCOMVisitor(info);
111      compositeVisitor->addVisitor(rbCOMVisitor, 900);
112    }
113    
114 <  //compositeVisitor->addVisitor(lipidVisitor, 900);
132 <  
133 <  //creat SSD atom visitor
114 >  //create SSD atom visitor
115    SSDAtomVisitor* ssdVisitor = new SSDAtomVisitor(info);
116    compositeVisitor->addVisitor(ssdVisitor, 800);
117 +
118 +  //create GBtail atom visitor
119 +  GBtailVisitor* gbtVisitor = new GBtailVisitor(info);
120 +  compositeVisitor->addVisitor(gbtVisitor, 790);
121 +
122 +  //create GBhead atom visitor
123 +  GBheadVisitor* gbhVisitor = new GBheadVisitor(info);
124 +  compositeVisitor->addVisitor(gbhVisitor, 789);
125    
126 <  LinearAtomVisitor* linearVisitor = new LinearAtomVisitor(info);
138 <  compositeVisitor->addVisitor(linearVisitor, 750);
139 <  
140 <  //creat default atom visitor
126 >  //create default atom visitor
127    DefaultAtomVisitor* defaultAtomVisitor = new DefaultAtomVisitor(info);
128    compositeVisitor->addVisitor(defaultAtomVisitor, 700);
129    
130 <  //creat waterType visitor
131 <  if(args_info.watertype_flag){
132 <    WaterTypeVisitor* waterTypeVisitor = new WaterTypeVisitor;
133 <    compositeVisitor->addVisitor(waterTypeVisitor, 600);
130 >  // if we gave the -w option, we want to skip the waters:
131 >  if (!args_info.water_given) {
132 >    //create waterType visitor
133 >    if(args_info.watertype_flag){
134 >      WaterTypeVisitor* waterTypeVisitor = new WaterTypeVisitor;
135 >      compositeVisitor->addVisitor(waterTypeVisitor, 600);
136 >    }
137 >  }
138 >  
139 >  if (args_info.basetype_flag) {
140 >    AtomNameVisitor* atomNameVisitor = new AtomNameVisitor(info);
141 >    compositeVisitor->addVisitor(atomNameVisitor, 550);    
142 >    cout << compositeVisitor->toString();
143    }
144    
145    //create ZconsVisitor
# Line 159 | Line 154 | int main(int argc, char* argv[]){
154      }
155    }
156    
157 <  //creat wrapping visitor
157 >  //create wrapping visitor
158    
159 <  if(args_info.periodicBox_flag){
160 <    WrappingVisitor* wrappingVisitor = new WrappingVisitor(info);
161 <    compositeVisitor->addVisitor(wrappingVisitor, 400);
159 >  //if(args_info.periodicBox_flag){
160 >  //  WrappingVisitor* wrappingVisitor = new WrappingVisitor(info);
161 >  //  compositeVisitor->addVisitor(wrappingVisitor, 400);
162 >  //}
163 >
164 >  //create replicate visitor
165 >  if(args_info.repeatX_given > 0 ||
166 >     args_info.repeatY_given > 0 ||
167 >     args_info.repeatZ_given > 0) {
168 >    Vector3i replicateOpt(args_info.repeatX_arg,
169 >                          args_info.repeatY_arg,
170 >                          args_info.repeatZ_arg);
171 >    ReplicateVisitor* replicateVisitor = new ReplicateVisitor(info,
172 >                                                              replicateOpt);
173 >    compositeVisitor->addVisitor(replicateVisitor, 300);
174    }
175 +
176 +
177 +  //create rotation visitor
178 +  if (args_info.refsele_given&& args_info.originsele_given) {
179 +    compositeVisitor->addVisitor(new LipidTransVisitor(info,
180 +                                                       args_info.originsele_arg,
181 +                                                       args_info.refsele_arg),
182 +                                 250);
183 +  } else if (args_info.refsele_given || args_info.originsele_given) {
184 +    cerr << "Both of --refsele and --originsele should appear by pair"
185 +         << endl;
186 +    exit(1);
187 +  }
188    
189 <  //creat replicate visitor
190 <  if(args_info.repeatX_given > 0 || args_info.repeatY_given > 0 ||args_info.repeatY_given > 0){
191 <    Vector3i replicateOpt(args_info.repeatX_arg, args_info.repeatY_arg, args_info.repeatZ_arg);
192 <    ReplicateVisitor* replicateVisitor = new ReplicateVisitor(info, replicateOpt);
193 <    compositeVisitor->addVisitor(replicateVisitor, 300);
189 >  //create xyzVisitor
190 >  XYZVisitor* xyzVisitor;
191 >
192 >  if (args_info.selection_given) {
193 >    xyzVisitor = new XYZVisitor(info, args_info.selection_arg);
194 >  } else {
195 >    xyzVisitor = new XYZVisitor(info);
196    }
197 +
198 +  if(args_info.velocities_flag){
199 +    printVel = true;
200 +    xyzVisitor->doVelocities(printVel);
201 +  }
202 +  if(args_info.forces_flag){
203 +    printFrc = true;
204 +    xyzVisitor->doForces(printFrc);
205 +  }
206 +  if(args_info.vectors_flag){
207 +    printVec = true;
208 +    xyzVisitor->doVectors(printVec);
209 +  }
210 +  if(args_info.charges_flag){
211 +    printChrg = true;
212 +    xyzVisitor->doCharges(printChrg);
213 +  }
214 +  if(args_info.efield_flag){
215 +    printField = true;
216 +    xyzVisitor->doElectricFields(printField);
217 +  }
218    
219 <  //creat xyzVisitor
177 <  XYZVisitor* xyzVisitor = new XYZVisitor(info);
178 <  compositeVisitor->addVisitor(xyzVisitor, 200);
219 >  compositeVisitor->addVisitor(xyzVisitor, 200);
220    
221 <  std::cout << compositeVisitor->toString();
181 <  
182 <  //creat prepareVisitor
221 >  //create prepareVisitor
222    PrepareVisitor* prepareVisitor = new PrepareVisitor();
223    
224    //open dump file
225    DumpReader* dumpReader = new DumpReader(info, dumpFileName);
226    int nframes = dumpReader->getNFrames();
227    
228 +  ofstream xyzStream(xyzFileName.c_str());
229    
190  std::ofstream xyzStream;    
191  xyzStream .open(xyzFileName.c_str());
192  
193  
230    SimInfo::MoleculeIterator miter;
231    Molecule::IntegrableObjectIterator  iiter;
232    Molecule::RigidBodyIterator rbIter;
233    Molecule* mol;
234 <  StuntDouble* integrableObject;
234 >  StuntDouble* sd;
235    RigidBody* rb;
236 <  
236 >  Vector3d molCom;
237 >  Vector3d newMolCom;
238 >  Vector3d displacement;
239 >  Mat3x3d hmat;
240 >  Snapshot* currentSnapshot;
241 >      
242    for (int i = 0; i < nframes; i += args_info.frame_arg){
243      dumpReader->readFrame(i);
244      
245 +    if (printFrc) forceMan->calcForces();
246 +    
247 +    //wrapping the molecule
248 +    if(args_info.periodicBox_flag) {
249 +      currentSnapshot = info->getSnapshotManager()->getCurrentSnapshot();    
250 +      for (mol = info->beginMolecule(miter); mol != NULL;
251 +           mol = info->nextMolecule(miter)) {
252 +        
253 +        molCom = mol->getCom();
254 +        newMolCom = molCom;
255 +        currentSnapshot->wrapVector(newMolCom);
256 +        displacement = newMolCom - molCom;
257 +
258 +        for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
259 +             sd = mol->nextIntegrableObject(iiter)) {  
260 +
261 +          sd->setPos(sd->getPos() + displacement);
262 +          
263 +        }
264 +      }    
265 +    }
266 +
267      //update atoms of rigidbody
268 <    for (mol = info->beginMolecule(miter); mol != NULL; mol = info->nextMolecule(miter)) {
268 >    for (mol = info->beginMolecule(miter); mol != NULL;
269 >         mol = info->nextMolecule(miter)) {
270        
271        //change the positions of atoms which belong to the rigidbodies
272 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
272 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
273 >           rb = mol->nextRigidBody(rbIter)) {
274 >
275          rb->updateAtoms();
276 +        if (printVel) rb->updateAtomVel();
277 +
278        }
279      }
280      
281      //prepare visit
282 <    for (mol = info->beginMolecule(miter); mol != NULL; mol = info->nextMolecule(miter)) {
283 <      for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
284 <           integrableObject = mol->nextIntegrableObject(iiter)) {
285 <        integrableObject->accept(prepareVisitor);
282 >    for (mol = info->beginMolecule(miter); mol != NULL;
283 >         mol = info->nextMolecule(miter)) {
284 >
285 >      for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
286 >           sd = mol->nextIntegrableObject(iiter)) {
287 >
288 >        sd->accept(prepareVisitor);
289 >
290        }
291      }
292      
293      //update visitor
294      compositeVisitor->update();
295 <    
295 >
296 >
297      //visit stuntdouble
298 <    for (mol = info->beginMolecule(miter); mol != NULL; mol = info->nextMolecule(miter)) {
299 <      for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
300 <           integrableObject = mol->nextIntegrableObject(iiter)) {
301 <        integrableObject->accept(compositeVisitor);
298 >    for (mol = info->beginMolecule(miter); mol != NULL;
299 >         mol = info->nextMolecule(miter)) {
300 >
301 >      for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
302 >           sd = mol->nextIntegrableObject(iiter)) {
303 >
304 >        sd->accept(compositeVisitor);
305 >
306        }
307      }
308      
# Line 233 | Line 310 | int main(int argc, char* argv[]){
310      xyzVisitor->clear();
311      
312    }//end for (int i = 0; i < nframes; i += args_info.frame_arg)
313 <  
313 >
314    xyzStream.close();
315 <  
239 <  
315 >  delete prepareVisitor;
316    delete compositeVisitor;
317    delete info;
242  
318   }

Comparing trunk/src/applications/dump2Xyz/Dump2XYZ.cpp (property svn:keywords):
Revision 250 by tim, Wed Jan 12 23:21:42 2005 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines