ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/nanoRodBuilder/nanorodBuilder.cpp
(Generate patch)

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

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines