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

Comparing trunk/OOPSE-4/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp (file contents):
Revision 2737 by chuckv, Tue Apr 25 22:54:55 2006 UTC vs.
Revision 3053 by gezelter, Wed Oct 18 19:35:07 2006 UTC

# Line 67 | Line 67 | void createMdFile(const std::string&oldMdFileName,
67   using namespace oopse;
68   void createMdFile(const std::string&oldMdFileName,
69                    const std::string&newMdFileName,
70 <                  int components,int* numMol);
70 >                  std::vector<int> numMol);
71  
72   int main(int argc, char *argv []) {
73    
# Line 78 | Line 78 | int main(int argc, char *argv []) {
78    gengetopt_args_info args_info;
79    std::string latticeType;
80    std::string inputFileName;
81 <  std::string outPrefix;
82 <  std::string outMdFileName;
83 <  std::string outInitFileName;
81 >  std::string outputFileName;
82  
85 –  
86 –  
87 –  Lattice *simpleLat;
83    MoLocator* locator;
89 –  int* numMol;
84    int nComponents;
85    double latticeConstant;
86    std::vector<double> lc;
93 –  double mass;                                                                                      
94 –  const double rhoConvertConst = 1.661;
95 –  double density;
96 –  double particleRadius;
97 –  
98 –  
87  
88 +  RealType particleRadius;
89 +
90    Mat3x3d hmat;
91    std::vector<Vector3d> latticePos;
92    std::vector<Vector3d> latticeOrt;
93 <  int numMolPerCell;
104 <  int nShells; /* Number of shells in nanoparticle*/
105 <  int numSites;
106 <  
93 >
94    DumpWriter *writer;
95    
96    // Parse Command Line Arguments
97    if (cmdline_parser(argc, argv, &args_info) != 0)
98      exit(1);
99 <  
113 <        
114 <        
99 >        
100    /* get lattice type */
101 <  latticeType = UpperCase(args_info.latticetype_arg);
102 <    
101 >  latticeType = "FCC";
102 >
103    /* get input file name */
104    if (args_info.inputs_num)
105      inputFileName = args_info.inputs[0];
106    else {
107 <    std::cerr << "You must specify a input file name.\n" << std::endl;
107 >    sprintf(painCave.errMsg, "No input .md file name was specified "
108 >            "on the command line");
109 >    painCave.isFatal = 1;
110      cmdline_parser_print_help();
111 <    exit(1);
111 >    simError();
112    }
113    
114    /* parse md file and set up the system */
115    SimCreator oldCreator;
116    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
117    
118 <  
132 <  /*calculate lattice constant (in Angstrom)
133 <    latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
134 <    1.0 / 3.0);*/
135 <  
136 <  latticeConstant = args_info.latticeCnst_arg;
118 >  latticeConstant = args_info.latticeConstant_arg;
119    particleRadius = args_info.radius_arg;
120    Globals* simParams = oldInfo->getSimParams();
121    
140 –  /* Find out how many different components in this simualtion */
141 –  nComponents =simParams->getNComponents();
142 –  
143 –  /*determine the output file names*/
144 –  if (args_info.output_given){
145 –    outInitFileName = args_info.output_arg;
146 –  }else{
147 –    outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
148 –  }
149 –  
150 –  std::cout <<"Before build shaped lattice. "<<std::endl;
151 –  
152 –  /* create Molocators */
153 –  locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
154 –  
122    /* Create nanoparticle */
123 <  shapedLatticeSpherical nanoParticle(latticeConstant,latticeType,particleRadius);
123 >  shapedLatticeSpherical nanoParticle(latticeConstant, latticeType,
124 >                                      particleRadius);
125    
158 –  std::cout <<"Before build getPoints. "<<std::endl;
126    /* Build a lattice and get lattice points for this lattice constant */
127 <  vector<Vector3d> nanoParticleSites = nanoParticle.getPoints();
127 >  vector<Vector3d> sites = nanoParticle.getSites();
128 >  vector<Vector3d> orientations = nanoParticle.getOrientations();
129 >  std::vector<int> vacancyTargets;
130 >  vector<bool> isVacancy;
131    
132 <  /* Get number of lattice sites */
133 <  numSites = nanoParticleSites.size();
132 >  Vector3d myLoc;
133 >  RealType myR;
134  
135 <  //std::cout <<"numSites are %d "<<numSites<<std::endl;
136 <  // std::cout <<"nComponents are %d "<<nComponents<<std::endl;
137 <  numMol = new int[nComponents];
138 <
135 >  for (int i = 0; i < sites.size(); i++)
136 >    isVacancy.push_back(false);
137 >
138 >  if (args_info.vacancyPercent_given) {
139 >    if (args_info.vacancyPercent_arg < 0.0 || args_info.vacancyPercent_arg > 100.0) {
140 >      sprintf(painCave.errMsg, "vacancyPercent was set to a non-sensical value.");
141 >      painCave.isFatal = 1;
142 >      simError();
143 >    } else {
144 >      RealType vF = args_info.vacancyPercent_arg / 100.0;
145 >      RealType vIR;
146 >      RealType vOR;
147 >      if (args_info.vacancyInnerRadius_given) {
148 >        vIR = args_info.vacancyInnerRadius_arg;
149 >      } else {
150 >        vIR = 0.0;
151 >      }
152 >      if (args_info.vacancyOuterRadius_given) {
153 >        vOR = args_info.vacancyOuterRadius_arg;
154 >      } else {
155 >        vOR = particleRadius;
156 >      }
157 >      if (vIR >= 0.0 && vOR <= particleRadius && vOR >= vIR) {
158 >        
159 >        for (int i = 0; i < sites.size(); i++) {
160 >          myLoc = sites[i];
161 >          myR = myLoc.length();
162 >          if (myR >= vIR && myR <= vOR) {
163 >            vacancyTargets.push_back(i);
164 >          }          
165 >        }
166 >        std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
167 >        
168 >        int nTargets = vacancyTargets.size();
169 >        vacancyTargets.resize((int)(vF * nTargets));
170 >        
171 >                  
172 >        sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
173 >                "\tsites between %lf and %lf.", vacancyTargets.size(),
174 >                vIR, vOR);
175 >        painCave.isFatal = 0;
176 >        simError();
177 >
178 >        isVacancy.clear();
179 >        for (int i = 0; i < sites.size(); i++) {
180 >          bool vac = false;
181 >          for (int j = 0; j < vacancyTargets.size(); j++) {
182 >            if (i == vacancyTargets[j]) vac = true;
183 >          }
184 >          isVacancy.push_back(vac);
185 >        }
186 >              
187 >      } else {
188 >        sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
189 >                "\tinner or outer radii.  Check their values.");
190 >        painCave.isFatal = 1;
191 >        simError();
192 >      }
193 >    }
194 >  }
195 >
196 >  /* Get number of lattice sites */
197 >  int nSites = sites.size() - vacancyTargets.size();
198 >
199 >  std::vector<Component*> components = simParams->getComponents();
200 >  std::vector<RealType> molFractions;
201 >  std::vector<RealType> shellRadii;
202 >  std::vector<RealType> molecularMasses;
203 >  std::vector<int> nMol;
204 >  std::map<int, int> componentFromSite;
205 >  nComponents = components.size();
206 >
207 >  if (args_info.molFraction_given && args_info.shellRadius_given) {
208 >    sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
209 >            "arguments, but not both!");
210 >    painCave.isFatal = 1;
211 >    simError();
212 >  }
213    
214 <  /* Random particle is the default case*/
215 <  if (!args_info.ShellRadius_given){
216 <    std::cout << "Creating a random nanoparticle" << std::endl;
217 <    /* Check to see if we have enough components */
218 <    if (nComponents != args_info.molFraction_given && nComponents != 1){
219 <      std::cerr << "Number of components does not equal molFraction occurances." << std::endl;
220 <      exit(1);
214 >  if (nComponents == 1) {
215 >    molFractions.push_back(1.0);    
216 >    shellRadii.push_back(particleRadius);
217 >  } else if (args_info.molFraction_given) {
218 >    if ((int)args_info.molFraction_given == nComponents) {
219 >      for (int i = 0; i < nComponents; i++) {
220 >        molFractions.push_back(args_info.molFraction_arg[i]);
221 >      }
222 >    } else if ((int)args_info.molFraction_given == nComponents-1) {
223 >      RealType remainingFraction = 1.0;
224 >      for (int i = 0; i < nComponents-1; i++) {
225 >        molFractions.push_back(args_info.molFraction_arg[i]);
226 >        remainingFraction -= molFractions[i];
227 >      }
228 >      molFractions.push_back(remainingFraction);
229 >    } else {    
230 >      sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out molFractions "
231 >              "for all of the components in the <MetaData> block.");
232 >      painCave.isFatal = 1;
233 >      simError();
234      }
235 <    /* Build the mole fractions and number of molecules of each type */  
236 <    int totComponents = 0;
237 <    for (int i = 0;i<nComponents-1;i++){ /* Figure out Percent for each component */
238 <      numMol[i] = int((double)numSites * args_info.molFraction_arg[i]);
239 <      std::cout<<numMol[i]<<std::endl;
240 <      totComponents += numMol[i];
235 >  } else if ((int)args_info.shellRadius_given) {
236 >    if ((int)args_info.shellRadius_given == nComponents) {
237 >      for (int i = 0; i < nComponents; i++) {
238 >        shellRadii.push_back(args_info.shellRadius_arg[i]);
239 >      }
240 >    } else if ((int)args_info.shellRadius_given == nComponents-1) {
241 >      for (int i = 0; i < nComponents-1; i++) {
242 >        shellRadii.push_back(args_info.shellRadius_arg[i]);
243 >      }
244 >      shellRadii.push_back(particleRadius);
245 >    } else {    
246 >      sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out the\n"
247 >              "\tshell radii for all of the components in the <MetaData> block.");
248 >      painCave.isFatal = 1;
249 >      simError();
250      }
251 <    numMol[nComponents-1] = numSites - totComponents;
252 <  
253 <    /* do the iPod thing, Shuffle da vector */
254 <    std::random_shuffle(nanoParticleSites.begin(), nanoParticleSites.end());
255 <  } else{ /*Handle core-shell with multiple components.*/
256 <    std::cout << "Creating a core-shell nanoparticle." << std::endl;
191 <    if (nComponents != args_info.ShellRadius_given + 1){
192 <      std::cerr << "Number of components does not equal ShellRadius occurances." << std::endl;
193 <      exit(1);
194 <    }  
251 >  } else {
252 >    sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
253 >            "\tbut have not specified either molFraction or shellRadius arguments.");
254 >    painCave.isFatal = 1;
255 >    simError();
256 >  }
257      
258 +  if (args_info.molFraction_given) {
259 +    RealType totalFraction = 0.0;
260 +    
261 +    /* Do some simple sanity checking*/
262 +    
263 +    for (int i = 0; i < nComponents; i++) {
264 +      if (molFractions.at(i) < 0.0) {
265 +        sprintf(painCave.errMsg, "One of the requested molFractions was"
266 +                " less than zero!");
267 +        painCave.isFatal = 1;
268 +        simError();
269 +      }
270 +      if (molFractions.at(i) > 1.0) {
271 +        sprintf(painCave.errMsg, "One of the requested molFractions was"
272 +                " greater than one!");
273 +        painCave.isFatal = 1;
274 +        simError();
275 +      }
276 +      totalFraction += molFractions.at(i);
277 +    }
278 +    if (abs(totalFraction - 1.0) > 1e-6) {
279 +      sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
280 +      painCave.isFatal = 1;
281 +      simError();
282 +    }
283 +    
284 +    int remaining = nSites;
285 +    for (int i=0; i < nComponents-1; i++) {    
286 +      nMol.push_back(int((RealType)nSites * molFractions.at(i)));
287 +      remaining -= nMol.at(i);
288 +    }
289 +    nMol.push_back(remaining);
290 +    
291 +    // recompute actual mol fractions and perform final sanity check:
292 +    
293 +    int totalMolecules = 0;
294 +    for (int i=0; i < nComponents; i++) {
295 +      molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
296 +      totalMolecules += nMol.at(i);
297 +    }
298 +    
299 +    if (totalMolecules != nSites) {
300 +      sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
301 +              "to the number of lattice sites!");
302 +      painCave.isFatal = 1;
303 +      simError();
304 +    }
305 +  } else {
306 +
307 +    for (int i = 0; i < shellRadii.size(); i++) {
308 +      if (shellRadii.at(i) > particleRadius + 1e-6 ) {
309 +        sprintf(painCave.errMsg, "One of the shellRadius values exceeds the particle Radius.");
310 +        painCave.isFatal = 1;
311 +        simError();
312 +      }
313 +      if (shellRadii.at(i) <= 0.0 ) {
314 +        sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!");
315 +        painCave.isFatal = 1;
316 +        simError();
317 +      }
318 +    }
319    }
320  
321 +  vector<int> ids;          
322 +  if ((int)args_info.molFraction_given){
323 +    sprintf(painCave.errMsg, "Creating a randomized spherical nanoparticle.");
324 +    painCave.isFatal = 0;
325 +    simError();
326 +    /* Random particle is the default case*/
327 +
328 +    for (int i = 0; i < sites.size(); i++)
329 +      if (!isVacancy[i]) ids.push_back(i);
330 +    
331 +    std::random_shuffle(ids.begin(), ids.end());
332 +    
333 +  } else{
334 +    sprintf(painCave.errMsg, "Creating a core-shell spherical nanoparticle.");
335 +    painCave.isFatal = 0;
336 +    simError();
337 +
338 +    RealType smallestSoFar;
339 +    int myComponent = -1;
340 +    nMol.clear();
341 +    nMol.resize(nComponents);
342 +
343 +    for (int i = 0; i < sites.size(); i++) {
344 +      myLoc = sites[i];
345 +      myR = myLoc.length();
346 +      smallestSoFar = particleRadius;      
347 +      if (!isVacancy[i]) {
348 +        for (int j = 0; j < nComponents; j++) {
349 +          if (myR <= shellRadii[j]) {
350 +            if (shellRadii[j] <= smallestSoFar) {
351 +              smallestSoFar = shellRadii[j];
352 +              myComponent = j;
353 +            }
354 +          }
355 +        }
356 +        componentFromSite[i] = myComponent;
357 +        nMol[myComponent]++;
358 +      }
359 +    }      
360 +  }
361    
362 <  //get the orientation of the cell sites
363 <  //for the same type of molecule in same lattice, it will not change
201 <  latticeOrt = nanoParticle.getPointsOrt();
202 <  std::cout<<"Orientational vector Size: "<< std::endl;
203 <  std::cout<<latticeOrt.size()<< std::endl;
204 <  
205 <  
206 <  
207 <  // needed for writing out new md file.
208 <  
209 <  outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
210 <  outMdFileName = outPrefix + ".md";
211 <  
362 >  outputFileName = args_info.output_arg;
363 >  
364    //creat new .md file on fly which corrects the number of molecule    
365 <  createMdFile(inputFileName, outMdFileName, nComponents,numMol);
365 >  createMdFile(inputFileName, outputFileName, nMol);
366    
367    if (oldInfo != NULL)
368      delete oldInfo;
369    
218 –  
219 –  // We need to read in new siminfo object.    
220 –  //parse md file and set up the system
221 –  //SimCreator NewCreator;
370    SimCreator newCreator;
371 <  SimInfo* NewInfo = newCreator.createSim(outMdFileName, false);
372 <  
225 <  
371 >  SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
372 >    
373    // Place molecules
374    Molecule* mol;
375    SimInfo::MoleculeIterator mi;
376    mol = NewInfo->beginMolecule(mi);
377 +
378    int l = 0;
379 <  for (mol = NewInfo->beginMolecule(mi); mol != NULL; mol = NewInfo->nextMolecule(mi)) {
232 <    locator->placeMol(nanoParticleSites[l], latticeOrt[l], mol);
233 <    l++;
234 <  }
235 <
379 >  int whichSite = 0;
380  
381 <
381 >  for (int i = 0; i < nComponents; i++){
382 >    locator = new MoLocator(NewInfo->getMoleculeStamp(i),
383 >                            NewInfo->getForceField());
384 >    
385 >    if (!args_info.molFraction_given) {
386 >      for (int n = 0; n < sites.size(); n++) {
387 >        if (!isVacancy[n]) {
388 >          if (componentFromSite[n] == i) {
389 >            mol = NewInfo->getMoleculeByGlobalIndex(l);
390 >            locator->placeMol(sites[n], orientations[n], mol);
391 >            l++;
392 >          }
393 >        }
394 >      }
395 >    } else {
396 >      for (int n = 0; n < nMol.at(i); n++) {
397 >        mol = NewInfo->getMoleculeByGlobalIndex(l);
398 >        locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
399 >        l++;
400 >      }
401 >    }
402 >  }
403    
404    //fill Hmat
405 <  hmat(0, 0)=  latticeConstant;
405 >  hmat(0, 0)=  10.0*particleRadius;
406    hmat(0, 1) = 0.0;
407    hmat(0, 2) = 0.0;
408    
409    hmat(1, 0) = 0.0;
410 <  hmat(1, 1) =  latticeConstant;
410 >  hmat(1, 1) =  10.0*particleRadius;
411    hmat(1, 2) = 0.0;
412    
413    hmat(2, 0) = 0.0;
414    hmat(2, 1) = 0.0;
415 <  hmat(2, 2) =  latticeConstant;
415 >  hmat(2, 2) =  10.0*particleRadius;
416    
417    //set Hmat
418    NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
419    
420    
421    //create dumpwriter and write out the coordinates
422 <  NewInfo->setFinalConfigFileName(outInitFileName);
258 <  writer = new DumpWriter(NewInfo);
422 >  writer = new DumpWriter(NewInfo, outputFileName);
423    
424    if (writer == NULL) {
425 <    std::cerr << "error in creating DumpWriter" << std::endl;
426 <    exit(1);
425 >    sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
426 >    painCave.isFatal = 1;
427 >    simError();
428    }
429    
430    writer->writeDump();
431 <  std::cout << "new initial configuration file: " << outInitFileName
432 <            << " is generated." << std::endl;
433 <  
434 <  //delete objects
435 <  
436 <  //delete oldInfo and oldSimSetup
437 <  
438 <  if (NewInfo != NULL)
439 <    delete NewInfo;
440 <  
276 <  if (writer != NULL)
277 <    delete writer;      
278 <  cmdline_parser_free(&args_info);
431 >
432 >  // deleting the writer will put the closing at the end of the dump file
433 >
434 >  delete writer;
435 >
436 >  // cleanup a by calling sim error.....
437 >  sprintf(painCave.errMsg, "A new OOPSE MD file called \"%s\" has been "
438 >          "generated.\n", outputFileName.c_str());
439 >  painCave.isFatal = 0;
440 >  simError();
441    return 0;
442   }
443  
444 < void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
445 <                  int components,int* numMol) {
444 > void createMdFile(const std::string&oldMdFileName,
445 >                  const std::string&newMdFileName,
446 >                  std::vector<int> nMol) {
447    ifstream oldMdFile;
448    ofstream newMdFile;
449    const int MAXLEN = 65535;
# Line 289 | Line 452 | void createMdFile(const std::string&oldMdFileName, con
452    //create new .md file based on old .md file
453    oldMdFile.open(oldMdFileName.c_str());
454    newMdFile.open(newMdFileName.c_str());
292 –  
455    oldMdFile.getline(buffer, MAXLEN);
456 <
456 >
457    int i = 0;
458    while (!oldMdFile.eof()) {
459 <    
459 >
460      //correct molecule number
461      if (strstr(buffer, "nMol") != NULL) {
462 <      if(i<components){
463 <        sprintf(buffer, "\tnMol = %i;", numMol[i]);                            
462 >      if(i<nMol.size()){
463 >        sprintf(buffer, "\tnMol = %i;", nMol.at(i));
464          newMdFile << buffer << std::endl;
465          i++;
466        }
# Line 310 | Line 472 | void createMdFile(const std::string&oldMdFileName, con
472    
473    oldMdFile.close();
474    newMdFile.close();
475 +
476 +  if (i != nMol.size()) {
477 +    sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
478 +            "\tstatements in component blocks.  Make sure that all\n"
479 +            "\tcomponents in the template file have nMol=1");
480 +    painCave.isFatal = 1;
481 +    simError();
482 +  }
483 +    
484   }
485  

Diff Legend

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