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.
#	Content
1	/*
2	* 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	* 1. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
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.
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	*
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	* [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>
44	#include <fstream>
45	#include <string>
46
47	#include "applications/dump2Xyz/Dump2XYZCmd.h"
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	#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
70	gengetopt_args_info args_info;
71	string dumpFileName;
72	string xyzFileName;
73
74	bool printVel(false);
75	bool printFrc(false);
76	bool printVec(false);
77	bool printChrg(false);
78
79	//parse the command line option
80	if (cmdline_parser (argc, argv, &args_info) != 0) {
81	exit(1) ;
82	}
83
84	//get the dumpfile name and meta-data file name
85	if (args_info.input_given){
86	dumpFileName = args_info.input_arg;
87	} else {
88	cerr << "Does not have input file name" << endl;
89	exit(1);
90	}
91
92	if (args_info.output_given){
93	xyzFileName = args_info.output_arg;
94	} else {
95	xyzFileName = dumpFileName;
96	xyzFileName = xyzFileName.substr(0, xyzFileName.rfind(".")) + ".xyz";
97	}
98
99	//parse md file and set up the system
100	SimCreator creator;
101	SimInfo* info = creator.createSim(dumpFileName, false);
102	ForceManager* forceMan = new ForceManager(info);
103
104	//create visitor list
105	CompositeVisitor* compositeVisitor = new CompositeVisitor();
106
107	//create RigidBody Visitor
108	if(args_info.rigidbody_flag){
109	RBCOMVisitor* rbCOMVisitor = new RBCOMVisitor(info);
110	compositeVisitor->addVisitor(rbCOMVisitor, 900);
111	}
112
113	//create SSD atom visitor
114	SSDAtomVisitor* ssdVisitor = new SSDAtomVisitor(info);
115	compositeVisitor->addVisitor(ssdVisitor, 800);
116
117	//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	//create default atom visitor
126	DefaultAtomVisitor* defaultAtomVisitor = new DefaultAtomVisitor(info);
127	compositeVisitor->addVisitor(defaultAtomVisitor, 700);
128
129	// 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	}
137
138	if (args_info.basetype_flag) {
139	AtomNameVisitor* atomNameVisitor = new AtomNameVisitor(info);
140	compositeVisitor->addVisitor(atomNameVisitor, 550);
141	cout << compositeVisitor->toString();
142	}
143
144	//create ZconsVisitor
145	if(args_info.zconstraint_flag){
146
147	ZConsVisitor* zconsVisitor = new ZConsVisitor(info);
148
149	if(zconsVisitor->haveZconsMol()) {
150	compositeVisitor->addVisitor(zconsVisitor, 500);
151	} else {
152	delete zconsVisitor;
153	}
154	}
155
156	//create wrapping visitor
157
158	//if(args_info.periodicBox_flag){
159	// WrappingVisitor* wrappingVisitor = new WrappingVisitor(info);
160	// compositeVisitor->addVisitor(wrappingVisitor, 400);
161	//}
162
163	//create replicate visitor
164	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	compositeVisitor->addVisitor(replicateVisitor, 300);
173	}
174
175
176	//create rotation visitor
177	if (args_info.refsele_given&& args_info.originsele_given) {
178	compositeVisitor->addVisitor(new LipidTransVisitor(info,
179	args_info.originsele_arg,
180	args_info.refsele_arg),
181	250);
182	} else if (args_info.refsele_given \|\| args_info.originsele_given) {
183	cerr << "Both of --refsele and --originsele should appear by pair"
184	<< endl;
185	exit(1);
186	}
187
188	//create xyzVisitor
189	XYZVisitor* xyzVisitor;
190
191	if (args_info.selection_given) {
192	xyzVisitor = new XYZVisitor(info, args_info.selection_arg);
193	} else {
194	xyzVisitor = new XYZVisitor(info);
195	}
196
197	if(args_info.velocities_flag){
198	printVel = true;
199	xyzVisitor->doVelocities(printVel);
200	}
201	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	compositeVisitor->addVisitor(xyzVisitor, 200);
215
216	//create prepareVisitor
217	PrepareVisitor* prepareVisitor = new PrepareVisitor();
218
219	//open dump file
220	DumpReader* dumpReader = new DumpReader(info, dumpFileName);
221	int nframes = dumpReader->getNFrames();
222
223	ofstream xyzStream(xyzFileName.c_str());
224
225	SimInfo::MoleculeIterator miter;
226	Molecule::IntegrableObjectIterator iiter;
227	Molecule::RigidBodyIterator rbIter;
228	Molecule* mol;
229	StuntDouble* integrableObject;
230	RigidBody* rb;
231	Vector3d molCom;
232	Vector3d newMolCom;
233	Vector3d displacement;
234	Mat3x3d hmat;
235	Snapshot* currentSnapshot;
236
237	for (int i = 0; i < nframes; i += args_info.frame_arg){
238	dumpReader->readFrame(i);
239
240	if (printFrc) forceMan->calcForces();
241
242	//wrapping the molecule
243	if(args_info.periodicBox_flag) {
244	currentSnapshot = info->getSnapshotManager()->getCurrentSnapshot();
245	for (mol = info->beginMolecule(miter); mol != NULL;
246	mol = info->nextMolecule(miter)) {
247	molCom = mol->getCom();
248	newMolCom = molCom;
249	currentSnapshot->wrapVector(newMolCom);
250	displacement = newMolCom - molCom;
251	for (integrableObject = mol->beginIntegrableObject(iiter);
252	integrableObject != NULL;
253	integrableObject = mol->nextIntegrableObject(iiter)) {
254	integrableObject->setPos(integrableObject->getPos() + displacement);
255	}
256	}
257	}
258	//update atoms of rigidbody
259	for (mol = info->beginMolecule(miter); mol != NULL;
260	mol = info->nextMolecule(miter)) {
261
262	//change the positions of atoms which belong to the rigidbodies
263	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
264	rb = mol->nextRigidBody(rbIter)) {
265	rb->updateAtoms();
266	if (printVel) rb->updateAtomVel();
267	}
268	}
269
270	//prepare visit
271	for (mol = info->beginMolecule(miter); mol != NULL;
272	mol = info->nextMolecule(miter)) {
273	for (integrableObject = mol->beginIntegrableObject(iiter);
274	integrableObject != NULL;
275	integrableObject = mol->nextIntegrableObject(iiter)) {
276	integrableObject->accept(prepareVisitor);
277	}
278	}
279
280	//update visitor
281	compositeVisitor->update();
282
283
284	//visit stuntdouble
285	for (mol = info->beginMolecule(miter); mol != NULL;
286	mol = info->nextMolecule(miter)) {
287	for (integrableObject = mol->beginIntegrableObject(iiter);
288	integrableObject != NULL;
289	integrableObject = mol->nextIntegrableObject(iiter)) {
290	integrableObject->accept(compositeVisitor);
291	}
292	}
293
294	xyzVisitor->writeFrame(xyzStream);
295	xyzVisitor->clear();
296
297	}//end for (int i = 0; i < nframes; i += args_info.frame_arg)
298
299	xyzStream.close();
300	delete prepareVisitor;
301	delete compositeVisitor;
302	delete info;
303	}