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 2164 by chuckv, Mon Apr 11 21:37:30 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"
55 < #include "applications/simpleBuilder/LatticeFactory.hpp"
56 < #include "applications/simpleBuilder/MoLocator.hpp"
57 < #include "applications/simpleBuilder/Lattice.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"
60   #include "brains/Register.hpp"
61   #include "brains/SimInfo.hpp"
62   #include "brains/SimCreator.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 <    
241 <    gengetopt_args_info args_info;
242 <    std::string latticeType;
243 <    std::string inputFileName;
244 <    std::string outPrefix;
245 <    std::string outMdFileName;
246 <    std::string outInitFileName;
247 <                std::string outGeomFileName;
248 <                
249 <                
250 <    BaseLattice *simpleLat;
251 <    int numMol;
252 <    double latticeConstant;
253 <    std::vector<double> lc;
254 <    double mass;
255 <    const double rhoConvertConst = 1.661;
256 <    double density;
90 <                double rodLength;
91 <                double rodDiameter;
92 <                
93 <                
94 <    int nx,
95 <    ny,
96 <    nz;
97 <    Mat3x3d hmat;
98 <    MoLocator *locator;
99 <    std::vector<Vector3d> latticePos;
100 <    std::vector<Vector3d> latticeOrt;
101 <    int numMolPerCell;
102 <    int curMolIndex;
103 <    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 <
114 <    if (!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)) {
115 <        std::cerr << latticeType << " is an invalid lattice type" << std::endl;
116 <        std::cerr << LatticeFactory::getInstance()->toString() << std::endl;
117 <        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  
130 –    //parse md file and set up the system
131 –    SimCreator oldCreator;
132 –    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
133 –                
134 –    if (oldInfo->getNMoleculeStamp()> 1) {
135 –        std::cerr << "can not build nanorod with more than one components"
136 –            << std::endl;
137 –        exit(1);
138 –    }
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);
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 <    if (simpleLat == NULL) {
323 <        std::cerr << "Error in creating lattice" << std::endl;
324 <        exit(1);
325 <    }
326 <
327 <    numMolPerCell = simpleLat->getNumSitesPerCell();
328 <
329 <    //calculate lattice constant (in Angstrom)
330 <    //latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
331 <    //                      1.0 / 3.0);
332 <
333 <                latticeConstant = args_info.latticeCnst_arg;
334 <                rodLength = args_info.length_arg;
335 <                rodDiameter = args_info.width_arg;
336 <                
337 <    //set lattice constant
338 <    lc.push_back(latticeConstant);
339 <    simpleLat->setLatticeConstant(lc);
340 <
341 <    
342 <    //determine the output file names  
343 <    if (args_info.output_given)
344 <        outInitFileName = args_info.output_arg;
345 <    else
346 <        outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
347 <    
348 <    //creat Molocator
349 <    locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
350 <
351 <                /*
352 <                Assume we are carving nanorod out of a cublic block of material and that
353 <                the shape the material will fit within that block....
354 <                The model in geometry builder assumes the long axis is in the y direction and the x-z plane is the
355 <                diameter of the particle.                
356 <                 */
357 <                // Number of Unit Cells in Length first
358 <                ny = (int)(rodLength/latticeConstant);
359 <                // Number of unit cells in Width
360 <                nx = (int)(rodDiameter/latticeConstant);
361 <                nz = (int)(rodDiameter/latticeConstant);
362 <                
363 <                
364 <                
190 <                // Create geometry for nanocrystal
191 <                //GeometryBuilder myGeometry(rodLength,rodDiameter);
192 <                
193 <                /*
194 <                 We have to build the system first to figure out how many molecules there are
195 <                 then create a md file and then actually build the system.
196 <                 */
197 <                
198 <                //place the molecules
199 <                
200 <    curMolIndex = 0;
201 <                
202 <    //get the orientation of the cell sites
203 <    //for the same type of molecule in same lattice, it will not change
204 <    latticeOrt = simpleLat->getLatticePointsOrt();
205 <        
206 <                
207 <                /*
208 <                 void BaseLattice::getLatticePointsPos(std::vector<Vector3d>& latticePos, int nx, int ny, int nz){
209 <                        
210 <                         latticePos.resize(nCellSites);
211 <                        
212 <                         for( int i=0;i < nCellSites;i++){
213 <                                
214 <                                 latticePos[i][0] = origin[0] + cellSitesPos[i][0] + cellLen[0] * (double(nx) - 0.5);
215 <                                 latticePos[i][1] = origin[1] + cellSitesPos[i][1] + cellLen[1] * (double(ny) - 0.5);
216 <                                 latticePos[i][2] = origin[2] + cellSitesPos[i][2] + cellLen[2] * (double(nz) - 0.5);    
217 <                         }
218 <                        
219 <                         */
220 <                
221 <                
222 <                
223 <                
224 <                numMol = 0;
225 <    for(int i = 0; i < nx; i++) {
226 <                        for(int j = 0; j < ny; j++) {
227 <                                for(int k = 0; k < nz; k++) {
228 <                                        //if (oldInfo->getNGlobalMolecules() != numMol) {
229 <                                        
230 <                                        
231 <                                        
232 <                                        //get the position of the cell sites
233 <                                        simpleLat->getLatticePointsPos(latticePos, i, j, k);
234 <                                        
235 <                                        for(int l = 0; l < numMolPerCell; l++) {
236 <                                                /*
237 <                                                if (myGeometry.isInsidePolyhedron(latticePos[l][0],latticePos[l][1],latticePos[l][2])){
238 <                                                        numMol++;
239 <                                                }
240 <                                                 */
241 <                                                numMol++;
242 <                                        }
243 <                                        }
244 <                                }
245 <                        }
246 <                
247 <                        
248 <                // needed for writing out new md file.
249 <                
250 <                outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
251 <                outMdFileName = outPrefix + ".md";
252 <                
253 <                //creat new .md file on fly which corrects the number of molecule    
254 <                createMdFile(inputFileName, outMdFileName, numMol);
255 <                
256 <                if (oldInfo != NULL)
257 <                        delete oldInfo;
258 <
259 <                        
260 <                        // We need to read in new siminfo object.      
261 <                        //parse md file and set up the system
262 <    //SimCreator NewCreator;
263 <                
264 <    SimInfo* NewInfo = oldCreator.createSim(outMdFileName, false);
265 <                
266 <        // This was so much fun the first time, lets do it again.
267 <                
268 <    Molecule* mol;
269 <    SimInfo::MoleculeIterator mi;
270 <    mol = NewInfo->beginMolecule(mi);
271 <                
272 <                for(int i = 0; i < nx; i++) {
273 <                        for(int j = 0; j < ny; j++) {
274 <                                for(int k = 0; k < nz; k++) {
275 <                                        
276 <                                        //get the position of the cell sites
277 <                                        simpleLat->getLatticePointsPos(latticePos, i, j, k);
278 <                                        
279 <                                        for(int l = 0; l < numMolPerCell; l++) {
280 <                                                if (mol != NULL) {
281 <                                                        /*
282 <                                                        if (myGeometry.isInsidePolyhedron(latticePos[l][0],latticePos[l][1],latticePos[l][2])){
283 <                                                                locator->placeMol(latticePos[l], latticeOrt[l], mol);
284 <                                                        }              
285 <                                                         */
286 <                                                        locator->placeMol(latticePos[l], latticeOrt[l], mol);
287 <                                                } else {
288 <                                                        std::cerr << std::endl;                    
289 <                                                }
290 <                                                mol = NewInfo->nextMolecule(mi);
291 <                                        }
292 <                                }
293 <                        }
294 <    }
295 <                
296 <                
297 <                
298 <    //fill Hmat
299 <    hmat(0, 0)= nx * latticeConstant;
300 <    hmat(0, 1) = 0.0;
301 <    hmat(0, 2) = 0.0;
302 <
303 <    hmat(1, 0) = 0.0;
304 <    hmat(1, 1) = ny * latticeConstant;
305 <    hmat(1, 2) = 0.0;
306 <
307 <    hmat(2, 0) = 0.0;
308 <    hmat(2, 1) = 0.0;
309 <    hmat(2, 2) = nz * latticeConstant;
310 <
311 <    //set Hmat
312 <    NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
313 <
314 <    
315 <    //create dumpwriter and write out the coordinates
316 <    NewInfo->setFinalConfigFileName(outInitFileName);
317 <    writer = new DumpWriter(NewInfo);
318 <
319 <    if (writer == NULL) {
320 <        std::cerr << "error in creating DumpWriter" << std::endl;
321 <        exit(1);
322 <    }
323 <
324 <    writer->writeDump();
325 <    std::cout << "new initial configuration file: " << outInitFileName
326 <        << " is generated." << std::endl;
327 <
328 <    //delete objects
329 <
330 <    //delete oldInfo and oldSimSetup
331 <                
332 <                                if (NewInfo != NULL)
333 <                        delete NewInfo;
334 <                
335 <    if (writer != NULL)
336 <                        delete writer;
337 <                
338 <    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
357 <        if (strstr(buffer, "nMol") != NULL) {
358 <            sprintf(buffer, "\tnMol = %i;", numMol);                            
359 <            newMdFile << buffer << std::endl;
360 <        } else
361 <            newMdFile << buffer << std::endl;
362 <
363 <        oldMdFile.getline(buffer, MAXLEN);
364 <    }
365 <
366 <    oldMdFile.close();
367 <    newMdFile.close();
390 >  }
391 >  
392 >  oldMdFile.close();
393 >  newMdFile.close();
394   }
395  

Diff Legend

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