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

Comparing trunk/OOPSE-3.0/src/applications/nanoRodBuilder/nanorodBuilder.cpp (file contents):
Revision 2181 by tim, Tue Apr 12 21:58:09 2005 UTC vs.
Revision 2196 by chuckv, Thu Apr 14 16:10:57 2005 UTC

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

Diff Legend

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