ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/applications/nanoRodBuilder/nanorodBuilder.cpp
Revision: 2491
Committed: Tue Dec 6 17:53:20 2005 UTC (18 years, 6 months ago) by chuckv
File size: 10717 byte(s)
Log Message:
Changes to build twinned crystal

File Contents

# User Rev Content
1 gezelter 2204 /*
2 chuckv 2164 * 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. 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
19     * notice, this list of conditions and the following disclaimer.
20     *
21     * 3. Redistributions in binary form must reproduce the above copyright
22     * notice, this list of conditions and the following disclaimer in the
23     * documentation and/or other materials provided with the
24     * distribution.
25     *
26     * This software is provided "AS IS," without a warranty of any
27     * kind. All express or implied conditions, representations and
28     * warranties, including any implied warranty of merchantability,
29     * fitness for a particular purpose or non-infringement, are hereby
30     * excluded. The University of Notre Dame and its licensors shall not
31     * be liable for any damages suffered by licensee as a result of
32     * using, modifying or distributing the software or its
33     * derivatives. In no event will the University of Notre Dame or its
34     * licensors be liable for any lost revenue, profit or data, or for
35     * direct, indirect, special, consequential, incidental or punitive
36     * damages, however caused and regardless of the theory of liability,
37     * arising out of the use of or inability to use software, even if the
38     * University of Notre Dame has been advised of the possibility of
39     * such damages.
40     */
41 gezelter 2190
42 chuckv 2164 #include <cstdlib>
43     #include <cstdio>
44     #include <cstring>
45     #include <cmath>
46     #include <iostream>
47     #include <string>
48     #include <map>
49     #include <fstream>
50    
51 gezelter 2190 #include "config.h"
52    
53 chuckv 2164 #include "nanorodBuilderCmd.h"
54 gezelter 2190 #ifdef HAVE_CGAL
55     #include "GeometryBuilder.hpp"
56     #endif
57 gezelter 2178 #include "lattice/LatticeFactory.hpp"
58     #include "utils/MoLocator.hpp"
59     #include "lattice/Lattice.hpp"
60 chuckv 2164 #include "brains/Register.hpp"
61     #include "brains/SimInfo.hpp"
62     #include "brains/SimCreator.hpp"
63     #include "io/DumpWriter.hpp"
64     #include "math/Vector3.hpp"
65     #include "math/SquareMatrix3.hpp"
66     #include "utils/StringUtils.hpp"
67    
68     using namespace std;
69     using namespace oopse;
70 gezelter 2190 void createMdFile(const std::string&oldMdFileName,
71     const std::string&newMdFileName,
72 chuckv 2164 int numMol);
73    
74     int main(int argc, char *argv []) {
75 gezelter 2190
76     //register force fields
77     registerForceFields();
78     registerLattice();
79    
80     gengetopt_args_info args_info;
81     std::string latticeType;
82     std::string inputFileName;
83     std::string outPrefix;
84     std::string outMdFileName;
85     std::string outInitFileName;
86     std::string outGeomFileName;
87    
88    
89     Lattice *simpleLat;
90     int numMol;
91     double latticeConstant;
92     std::vector<double> lc;
93     double mass;
94     const double rhoConvertConst = 1.661;
95     double density;
96     double rodLength;
97     double rodDiameter;
98    
99    
100     int nx, ny, nz;
101     Mat3x3d hmat;
102     MoLocator *locator;
103     std::vector<Vector3d> latticePos;
104     std::vector<Vector3d> latticeOrt;
105     int numMolPerCell;
106     DumpWriter *writer;
107    
108     // parse command line arguments
109     if (cmdline_parser(argc, argv, &args_info) != 0)
110     exit(1);
111    
112    
113 gezelter 2204 // Check for lib CGAL, if we don't have it, we should exit....
114 chuckv 2196
115     #ifndef HAVE_CGAL
116 gezelter 2204 std::cerr << "nanoRodBuilder requires libCGAL to function, please rebuild OOPSE with libCGAL."
117     << std::endl;
118     exit(1);
119 chuckv 2196 #endif
120    
121    
122    
123 gezelter 2190 //get lattice type
124     latticeType = UpperCase(args_info.latticetype_arg);
125 chuckv 2239
126 gezelter 2190 //get input file name
127     if (args_info.inputs_num)
128     inputFileName = args_info.inputs[0];
129     else {
130     std::cerr << "You must specify a input file name.\n" << std::endl;
131     cmdline_parser_print_help();
132     exit(1);
133     }
134    
135     //parse md file and set up the system
136     SimCreator oldCreator;
137     SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
138    
139     if (oldInfo->getNMoleculeStamp()> 1) {
140     std::cerr << "can not build nanorod with more than one components"
141     << std::endl;
142     exit(1);
143     }
144    
145     //get mass of molecule.
146    
147     mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
148    
149     //creat lattice
150     simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
151    
152     if (simpleLat == NULL) {
153     std::cerr << "Error in creating lattice" << std::endl;
154     exit(1);
155     }
156    
157     numMolPerCell = simpleLat->getNumSitesPerCell();
158    
159     //calculate lattice constant (in Angstrom)
160     //latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
161     // 1.0 / 3.0);
162    
163     latticeConstant = args_info.latticeCnst_arg;
164     rodLength = args_info.length_arg;
165     rodDiameter = args_info.width_arg;
166    
167     //set lattice constant
168     lc.push_back(latticeConstant);
169     simpleLat->setLatticeConstant(lc);
170    
171    
172     //determine the output file names
173     if (args_info.output_given)
174     outInitFileName = args_info.output_arg;
175     else
176     outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
177    
178 chuckv 2215
179    
180    
181    
182    
183 gezelter 2190 //creat Molocator
184     locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
185    
186     /*
187     Assume we are carving nanorod out of a cublic block of material and that
188     the shape the material will fit within that block....
189     The model in geometry builder assumes the long axis is in the y direction and the x-z plane is the
190     diameter of the particle.
191     */
192     // Number of Unit Cells in Length first
193     ny = (int)(rodLength/latticeConstant);
194     // Number of unit cells in Width
195     nx = (int)(rodDiameter/latticeConstant);
196     nz = (int)(rodDiameter/latticeConstant);
197    
198    
199    
200     // Create geometry for nanocrystal
201 chuckv 2196 #ifdef HAVE_CGAL
202 chuckv 2218 GeometryBuilder *myGeometry;
203     // GeometryBuilder myGeometry(rodLength,rodDiameter);
204     if (args_info.twinnedCrystal_flag){
205     myGeometry = new GeometryBuilder(rodLength,rodDiameter);
206     }
207     else{
208     myGeometry = new GeometryBuilder(rodLength,rodDiameter);
209     }
210    
211 chuckv 2215 if (args_info.genGeomview_given){
212     if (args_info.genGeomview_flag){
213     outGeomFileName = getPrefix(inputFileName.c_str()) + ".off";
214 chuckv 2218 myGeometry->dumpGeometry(outGeomFileName);
215 chuckv 2215 }
216     }
217 chuckv 2218
218 gezelter 2190 #endif
219    
220     /*
221     We have to build the system first to figure out how many molecules
222     there are then create a md file and then actually build the
223     system.
224     */
225    
226     //place the molecules
227    
228 chuckv 2218
229 gezelter 2190
230     //get the orientation of the cell sites
231     //for the same type of molecule in same lattice, it will not change
232     latticeOrt = simpleLat->getLatticePointsOrt();
233    
234    
235    
236     numMol = 0;
237 chuckv 2215 for(int i = -nx; i < nx; i++) {
238     for(int j = -ny; j < ny; j++) {
239     for(int k = -nz; k < nz; k++) {
240 gezelter 2190 //if (oldInfo->getNGlobalMolecules() != numMol) {
241 chuckv 2164
242 gezelter 2190
243    
244     //get the position of the cell sites
245     simpleLat->getLatticePointsPos(latticePos, i, j, k);
246    
247     for(int l = 0; l < numMolPerCell; l++) {
248 chuckv 2164
249 chuckv 2196 #ifdef HAVE_CGAL
250 chuckv 2218 if (myGeometry->isInsidePolyhedron(latticePos[l][0],latticePos[l][1],latticePos[l][2])){
251 gezelter 2190 numMol++;
252     }
253 chuckv 2218 #else
254     numMol++;
255 gezelter 2190 #endif
256     }
257     }
258 chuckv 2164 }
259 gezelter 2190 }
260 chuckv 2218
261 gezelter 2190
262     // needed for writing out new md file.
263    
264     outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
265     outMdFileName = outPrefix + ".md";
266    
267     //creat new .md file on fly which corrects the number of molecule
268     createMdFile(inputFileName, outMdFileName, numMol);
269    
270     if (oldInfo != NULL)
271     delete oldInfo;
272    
273    
274     // We need to read in new siminfo object.
275     //parse md file and set up the system
276     //SimCreator NewCreator;
277    
278     SimInfo* NewInfo = oldCreator.createSim(outMdFileName, false);
279    
280     // This was so much fun the first time, lets do it again.
281    
282     Molecule* mol;
283     SimInfo::MoleculeIterator mi;
284     mol = NewInfo->beginMolecule(mi);
285 chuckv 2218
286    
287 chuckv 2215 for(int i = -nx; i < nx; i++) {
288     for(int j = -ny; j < ny; j++) {
289     for(int k = -nz; k < nz; k++) {
290    
291     //get the position of the cell sites
292     simpleLat->getLatticePointsPos(latticePos, i, j, k);
293    
294     for(int l = 0; l < numMolPerCell; l++) {
295 chuckv 2218 #ifdef HAVE_CGAL
296     if (myGeometry->isInsidePolyhedron(latticePos[l][0],latticePos[l][1],latticePos[l][2])){
297     #endif
298 chuckv 2215 if (mol != NULL) {
299     locator->placeMol(latticePos[l], latticeOrt[l], mol);
300     } else {
301 chuckv 2218 std::cerr<<"Error in placing molecule " << std::endl;
302 chuckv 2215 }
303     mol = NewInfo->nextMolecule(mi);
304 chuckv 2218 #ifdef HAVE_CGAL
305 chuckv 2215 }
306 chuckv 2218 #endif
307 chuckv 2215 }
308 gezelter 2190 }
309 chuckv 2215 }
310 gezelter 2190 }
311    
312 chuckv 2218
313 gezelter 2190
314     //fill Hmat
315     hmat(0, 0)= nx * latticeConstant;
316     hmat(0, 1) = 0.0;
317     hmat(0, 2) = 0.0;
318    
319     hmat(1, 0) = 0.0;
320     hmat(1, 1) = ny * latticeConstant;
321     hmat(1, 2) = 0.0;
322    
323     hmat(2, 0) = 0.0;
324     hmat(2, 1) = 0.0;
325     hmat(2, 2) = nz * latticeConstant;
326    
327     //set Hmat
328     NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
329    
330    
331     //create dumpwriter and write out the coordinates
332     NewInfo->setFinalConfigFileName(outInitFileName);
333     writer = new DumpWriter(NewInfo);
334    
335     if (writer == NULL) {
336     std::cerr << "error in creating DumpWriter" << std::endl;
337     exit(1);
338 gezelter 2204 }
339 gezelter 2190
340 chuckv 2491 writer->writeEor();
341 gezelter 2190 std::cout << "new initial configuration file: " << outInitFileName
342     << " is generated." << std::endl;
343    
344     //delete objects
345    
346     //delete oldInfo and oldSimSetup
347    
348     if (NewInfo != NULL)
349     delete NewInfo;
350    
351     if (writer != NULL)
352     delete writer;
353     delete simpleLat;
354     cmdline_parser_free(&args_info);
355     return 0;
356 chuckv 2164 }
357    
358     void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
359     int numMol) {
360 gezelter 2190 ifstream oldMdFile;
361     ofstream newMdFile;
362     const int MAXLEN = 65535;
363     char buffer[MAXLEN];
364    
365     //create new .md file based on old .md file
366     oldMdFile.open(oldMdFileName.c_str());
367     newMdFile.open(newMdFileName.c_str());
368    
369     oldMdFile.getline(buffer, MAXLEN);
370    
371     while (!oldMdFile.eof()) {
372    
373     //correct molecule number
374     if (strstr(buffer, "nMol") != NULL) {
375     sprintf(buffer, "\tnMol = %i;", numMol);
376     newMdFile << buffer << std::endl;
377     } else
378     newMdFile << buffer << std::endl;
379    
380 chuckv 2164 oldMdFile.getline(buffer, MAXLEN);
381 gezelter 2190 }
382    
383     oldMdFile.close();
384     newMdFile.close();
385 chuckv 2164 }
386