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

Diff Legend

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