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

Diff Legend

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