applications/dump2Xyz/Dump2XYZ.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 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.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * 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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
#include <iostream>
#include <fstream>
#include <string>

#include "applications/dump2Xyz/Dump2XYZCmd.h"
#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"
#include "visitors/ZconsVisitor.hpp"
#include "selection/SelectionEvaluator.hpp"
#include "selection/SelectionManager.hpp"
#include "visitors/LipidTransVisitor.hpp"
#include "visitors/AtomNameVisitor.hpp"

using namespace OpenMD;

using namespace std;
int main(int argc, char* argv[]){
  
  gengetopt_args_info args_info;
  string dumpFileName;
  string xyzFileName;

  bool printVel(false);
  bool printFrc(false);
  bool printVec(false);
  bool printChrg(false);
  
  //parse the command line option
  if (cmdline_parser (argc, argv, &args_info) != 0) {
    exit(1) ;
  }
  
  //get the dumpfile name and meta-data file name
  if (args_info.input_given){
    dumpFileName = args_info.input_arg;
  } else {
    cerr << "Does not have input file name" << endl;
    exit(1);
  }
  
  if (args_info.output_given){
    xyzFileName = args_info.output_arg;
  } else {
    xyzFileName = dumpFileName;
    xyzFileName = xyzFileName.substr(0, xyzFileName.rfind(".")) + ".xyz";
  }
  
  //parse md file and set up the system
  SimCreator creator;
  SimInfo* info = creator.createSim(dumpFileName, false);
  ForceManager* forceMan = new ForceManager(info);
  
  //create visitor list
  CompositeVisitor* compositeVisitor = new CompositeVisitor();
  
  //create RigidBody Visitor
  if(args_info.rigidbody_flag){
    RBCOMVisitor* rbCOMVisitor = new RBCOMVisitor(info);
    compositeVisitor->addVisitor(rbCOMVisitor, 900);
  }
  
  //create SSD atom visitor
  SSDAtomVisitor* ssdVisitor = new SSDAtomVisitor(info);
  compositeVisitor->addVisitor(ssdVisitor, 800);
  
  //create GBtail atom visitor
  GBtailVisitor* gbtVisitor = new GBtailVisitor(info);
  compositeVisitor->addVisitor(gbtVisitor, 790);
  
  //create GBhead atom visitor
  GBheadVisitor* gbhVisitor = new GBheadVisitor(info);
  compositeVisitor->addVisitor(gbhVisitor, 789);
  
  //create default atom visitor
  DefaultAtomVisitor* defaultAtomVisitor = new DefaultAtomVisitor(info);
  compositeVisitor->addVisitor(defaultAtomVisitor, 700);
  
  // 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);    
    cout << compositeVisitor->toString();
  }
  
  //create ZconsVisitor
  if(args_info.zconstraint_flag){

    ZConsVisitor* zconsVisitor = new ZConsVisitor(info);
    
    if(zconsVisitor->haveZconsMol()) {
      compositeVisitor->addVisitor(zconsVisitor, 500);
    } else {
      delete zconsVisitor;
    }
  }
  
  //create wrapping visitor
  
  //if(args_info.periodicBox_flag){
  //  WrappingVisitor* wrappingVisitor = new WrappingVisitor(info);
  //  compositeVisitor->addVisitor(wrappingVisitor, 400);
  //}

  //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); 
  } else if (args_info.refsele_given || args_info.originsele_given) {
    cerr << "Both of --refsele and --originsele should appear by pair" 
              << endl;
    exit(1);
  }
  
  //create xyzVisitor
  XYZVisitor* xyzVisitor;

  if (args_info.selection_given) {
    xyzVisitor = new XYZVisitor(info, args_info.selection_arg);
  } else {
    xyzVisitor = new XYZVisitor(info);
  }

  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);
  }
  
  compositeVisitor->addVisitor(xyzVisitor, 200); 
  
  //create prepareVisitor
  PrepareVisitor* prepareVisitor = new PrepareVisitor();
  
  //open dump file
  DumpReader* dumpReader = new DumpReader(info, dumpFileName);
  int nframes = dumpReader->getNFrames();
  
  ofstream xyzStream(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();
    
    //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)) {
      
      //change the positions of atoms which belong to the rigidbodies
      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;
           integrableObject = mol->nextIntegrableObject(iiter)) {
        integrableObject->accept(prepareVisitor);
      }
    }
    
    //update visitor
    compositeVisitor->update();


    //visit stuntdouble
    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);
      }
    }
    
    xyzVisitor->writeFrame(xyzStream);
    xyzVisitor->clear();
    
  }//end for (int i = 0; i < nframes; i += args_info.frame_arg)
 
  xyzStream.close();
  delete prepareVisitor; 
  delete compositeVisitor;
  delete info;
}
Revision:	1725
Committed:	Sat May 26 18:13:43 2012 UTC (12 years, 11 months ago) by gezelter
File size:	10109 byte(s)
Log Message:	Individual ForceField classes have been removed (they were essentially all duplicates anyway). ForceField has moved to brains, and since only one force field is in play at any time, the ForceFieldFactory and Register methods have been removed.
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	246	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	gezelter	246	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
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	gezelter	1390	*
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, 24107 (2008).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	246	*/
42
43	tim	12	#include <iostream>
44			#include <fstream>
45			#include <string>
46
47			#include "applications/dump2Xyz/Dump2XYZCmd.h"
48	gezelter	246	#include "brains/Register.hpp"
49			#include "brains/SimCreator.hpp"
50			#include "brains/SimInfo.hpp"
51	gezelter	1456	#include "brains/ForceManager.hpp"
52	gezelter	246	#include "io/DumpReader.hpp"
53			#include "utils/simError.h"
54	tim	12	#include "visitors/AtomVisitor.hpp"
55	gezelter	1455	#include "visitors/ReplacementVisitor.hpp"
56	tim	12	#include "visitors/CompositeVisitor.hpp"
57			#include "visitors/RigidBodyVisitor.hpp"
58			#include "visitors/OtherVisitor.hpp"
59			#include "visitors/ZconsVisitor.hpp"
60	tim	291	#include "selection/SelectionEvaluator.hpp"
61	tim	324	#include "selection/SelectionManager.hpp"
62	tim	369	#include "visitors/LipidTransVisitor.hpp"
63	tim	543	#include "visitors/AtomNameVisitor.hpp"
64	tim	369
65	gezelter	1390	using namespace OpenMD;
66	tim	12
67	gezelter	1706	using namespace std;
68	tim	12	int main(int argc, char* argv[]){
69	gezelter	246
70	tim	12	gengetopt_args_info args_info;
71	gezelter	1706	string dumpFileName;
72			string xyzFileName;
73
74			bool printVel(false);
75			bool printFrc(false);
76			bool printVec(false);
77			bool printChrg(false);
78	tim	12
79	gezelter	246	//parse the command line option
80			if (cmdline_parser (argc, argv, &args_info) != 0) {
81			exit(1) ;
82			}
83	tim	12
84			//get the dumpfile name and meta-data file name
85			if (args_info.input_given){
86			dumpFileName = args_info.input_arg;
87	gezelter	246	} else {
88	gezelter	1706	cerr << "Does not have input file name" << endl;
89	tim	12	exit(1);
90			}
91	gezelter	246
92	tim	12	if (args_info.output_given){
93			xyzFileName = args_info.output_arg;
94	gezelter	246	} else {
95	tim	12	xyzFileName = dumpFileName;
96			xyzFileName = xyzFileName.substr(0, xyzFileName.rfind(".")) + ".xyz";
97			}
98	gezelter	246
99	tim	12	//parse md file and set up the system
100	gezelter	246	SimCreator creator;
101	gezelter	1026	SimInfo* info = creator.createSim(dumpFileName, false);
102	gezelter	1456	ForceManager* forceMan = new ForceManager(info);
103	tim	12
104	chrisfen	527	//create visitor list
105	gezelter	246	CompositeVisitor* compositeVisitor = new CompositeVisitor();
106	gezelter	1455
107	chrisfen	527	//create RigidBody Visitor
108	tim	12	if(args_info.rigidbody_flag){
109	gezelter	246	RBCOMVisitor* rbCOMVisitor = new RBCOMVisitor(info);
110	tim	12	compositeVisitor->addVisitor(rbCOMVisitor, 900);
111			}
112	gezelter	246
113	chrisfen	527	//create SSD atom visitor
114	gezelter	246	SSDAtomVisitor* ssdVisitor = new SSDAtomVisitor(info);
115	tim	12	compositeVisitor->addVisitor(ssdVisitor, 800);
116	gezelter	246
117	gezelter	1629	//create GBtail atom visitor
118			GBtailVisitor* gbtVisitor = new GBtailVisitor(info);
119			compositeVisitor->addVisitor(gbtVisitor, 790);
120
121			//create GBhead atom visitor
122			GBheadVisitor* gbhVisitor = new GBheadVisitor(info);
123			compositeVisitor->addVisitor(gbhVisitor, 789);
124
125	chrisfen	527	//create default atom visitor
126	gezelter	246	DefaultAtomVisitor* defaultAtomVisitor = new DefaultAtomVisitor(info);
127	tim	12	compositeVisitor->addVisitor(defaultAtomVisitor, 700);
128	gezelter	246
129	cli2	1303	// if we gave the -w option, we want to skip the waters:
130			if (!args_info.water_given) {
131			//create waterType visitor
132			if(args_info.watertype_flag){
133			WaterTypeVisitor* waterTypeVisitor = new WaterTypeVisitor;
134			compositeVisitor->addVisitor(waterTypeVisitor, 600);
135			}
136	gezelter	1455	}
137
138	tim	543	if (args_info.basetype_flag) {
139	gezelter	1455	AtomNameVisitor* atomNameVisitor = new AtomNameVisitor(info);
140			compositeVisitor->addVisitor(atomNameVisitor, 550);
141	gezelter	1706	cout << compositeVisitor->toString();
142	tim	543	}
143	gezelter	246
144	tim	12	//create ZconsVisitor
145			if(args_info.zconstraint_flag){
146	gezelter	246
147			ZConsVisitor* zconsVisitor = new ZConsVisitor(info);
148	tim	12
149	gezelter	246	if(zconsVisitor->haveZconsMol()) {
150	tim	12	compositeVisitor->addVisitor(zconsVisitor, 500);
151	gezelter	246	} else {
152	tim	12	delete zconsVisitor;
153	gezelter	246	}
154	tim	12	}
155	gezelter	246
156	chrisfen	527	//create wrapping visitor
157	gezelter	246
158	tim	968	//if(args_info.periodicBox_flag){
159			// WrappingVisitor* wrappingVisitor = new WrappingVisitor(info);
160			// compositeVisitor->addVisitor(wrappingVisitor, 400);
161			//}
162	tim	369
163	chrisfen	527	//create replicate visitor
164	gezelter	1456	if(args_info.repeatX_given > 0 \|\|
165			args_info.repeatY_given > 0 \|\|
166			args_info.repeatY_given > 0) {
167			Vector3i replicateOpt(args_info.repeatX_arg,
168			args_info.repeatY_arg,
169			args_info.repeatZ_arg);
170			ReplicateVisitor* replicateVisitor = new ReplicateVisitor(info,
171			replicateOpt);
172	tim	12	compositeVisitor->addVisitor(replicateVisitor, 300);
173			}
174	tim	369
175
176			//create rotation visitor
177			if (args_info.refsele_given&& args_info.originsele_given) {
178	gezelter	1456	compositeVisitor->addVisitor(new LipidTransVisitor(info,
179			args_info.originsele_arg,
180			args_info.refsele_arg),
181			250);
182	tim	369	} else if (args_info.refsele_given \|\| args_info.originsele_given) {
183	gezelter	1706	cerr << "Both of --refsele and --originsele should appear by pair"
184			<< endl;
185	tim	369	exit(1);
186			}
187	gezelter	1456
188	chrisfen	527	//create xyzVisitor
189	tim	413	XYZVisitor* xyzVisitor;
190	chuckv	1117
191	tim	413	if (args_info.selection_given) {
192			xyzVisitor = new XYZVisitor(info, args_info.selection_arg);
193			} else {
194			xyzVisitor = new XYZVisitor(info);
195			}
196	chuckv	1117
197	gezelter	1456	if(args_info.velocities_flag){
198			printVel = true;
199			xyzVisitor->doVelocities(printVel);
200	chuckv	1117	}
201	gezelter	1456	if(args_info.forces_flag){
202			printFrc = true;
203			xyzVisitor->doForces(printFrc);
204			}
205			if(args_info.vectors_flag){
206			printVec = true;
207			xyzVisitor->doVectors(printVec);
208			}
209			if(args_info.charges_flag){
210			printChrg = true;
211			xyzVisitor->doCharges(printChrg);
212			}
213
214	gezelter	1455	compositeVisitor->addVisitor(xyzVisitor, 200);
215	gezelter	246
216	chrisfen	527	//create prepareVisitor
217	gezelter	246	PrepareVisitor* prepareVisitor = new PrepareVisitor();
218
219	tim	12	//open dump file
220	gezelter	246	DumpReader* dumpReader = new DumpReader(info, dumpFileName);
221			int nframes = dumpReader->getNFrames();
222
223	gezelter	1706	ofstream xyzStream(xyzFileName.c_str());
224	tim	12
225	gezelter	246	SimInfo::MoleculeIterator miter;
226			Molecule::IntegrableObjectIterator iiter;
227			Molecule::RigidBodyIterator rbIter;
228			Molecule* mol;
229			StuntDouble* integrableObject;
230			RigidBody* rb;
231	tim	968	Vector3d molCom;
232			Vector3d newMolCom;
233			Vector3d displacement;
234			Mat3x3d hmat;
235			Snapshot* currentSnapshot;
236	tim	413
237	tim	12	for (int i = 0; i < nframes; i += args_info.frame_arg){
238	gezelter	246	dumpReader->readFrame(i);
239	gezelter	1456
240	gezelter	1464	if (printFrc) forceMan->calcForces();
241	gezelter	1456
242	tim	968	//wrapping the molecule
243			if(args_info.periodicBox_flag) {
244			currentSnapshot = info->getSnapshotManager()->getCurrentSnapshot();
245	gezelter	1456	for (mol = info->beginMolecule(miter); mol != NULL;
246			mol = info->nextMolecule(miter)) {
247	tim	968	molCom = mol->getCom();
248			newMolCom = molCom;
249			currentSnapshot->wrapVector(newMolCom);
250			displacement = newMolCom - molCom;
251	gezelter	1456	for (integrableObject = mol->beginIntegrableObject(iiter);
252			integrableObject != NULL;
253	tim	968	integrableObject = mol->nextIntegrableObject(iiter)) {
254			integrableObject->setPos(integrableObject->getPos() + displacement);
255			}
256			}
257			}
258	tim	12	//update atoms of rigidbody
259	gezelter	1456	for (mol = info->beginMolecule(miter); mol != NULL;
260			mol = info->nextMolecule(miter)) {
261	gezelter	246
262			//change the positions of atoms which belong to the rigidbodies
263	gezelter	1456	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
264			rb = mol->nextRigidBody(rbIter)) {
265	gezelter	246	rb->updateAtoms();
266	gezelter	1456	if (printVel) rb->updateAtomVel();
267	gezelter	246	}
268			}
269	tim	12
270			//prepare visit
271	gezelter	1456	for (mol = info->beginMolecule(miter); mol != NULL;
272			mol = info->nextMolecule(miter)) {
273			for (integrableObject = mol->beginIntegrableObject(iiter);
274			integrableObject != NULL;
275	gezelter	246	integrableObject = mol->nextIntegrableObject(iiter)) {
276			integrableObject->accept(prepareVisitor);
277			}
278			}
279
280	tim	12	//update visitor
281			compositeVisitor->update();
282	tim	291
283	tim	413
284	tim	12	//visit stuntdouble
285	gezelter	1456	for (mol = info->beginMolecule(miter); mol != NULL;
286			mol = info->nextMolecule(miter)) {
287			for (integrableObject = mol->beginIntegrableObject(iiter);
288			integrableObject != NULL;
289	gezelter	246	integrableObject = mol->nextIntegrableObject(iiter)) {
290			integrableObject->accept(compositeVisitor);
291			}
292	tim	12	}
293	gezelter	246
294	tim	12	xyzVisitor->writeFrame(xyzStream);
295			xyzVisitor->clear();
296
297			}//end for (int i = 0; i < nframes; i += args_info.frame_arg)
298	gezelter	1282
299	tim	12	xyzStream.close();
300	tim	466	delete prepareVisitor;
301	tim	12	delete compositeVisitor;
302			delete info;
303			}