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

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

Diff Legend

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