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 3046 by gezelter, Sat Oct 14 20:21:26 2006 UTC vs.
Revision 3053 by gezelter, Wed Oct 18 19:35:07 2006 UTC

# Line 80 | Line 80 | int main(int argc, char *argv []) {
80    std::string inputFileName;
81    std::string outputFileName;
82  
83 –  Lattice *simpleLat;
83    MoLocator* locator;
85 –  int* numMol;
84    int nComponents;
85    double latticeConstant;
86    std::vector<double> lc;
89 –  double mass;                                                                          
90 –  const double rhoConvertConst = 1.661;
91 –  double density;
92 –  double particleRadius;
87  
88 +  RealType particleRadius;
89 +
90    Mat3x3d hmat;
91    std::vector<Vector3d> latticePos;
92    std::vector<Vector3d> latticeOrt;
97 –  int numMolPerCell;
98 –  int nShells; /* Number of shells in nanoparticle*/
93  
94    DumpWriter *writer;
95    
# Line 110 | Line 104 | int main(int argc, char *argv []) {
104    if (args_info.inputs_num)
105      inputFileName = args_info.inputs[0];
106    else {
107 <      sprintf(painCave.errMsg, "No input .md file name was specified"
108 <              "on the command line");
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      simError();
# Line 121 | Line 115 | int main(int argc, char *argv []) {
115    SimCreator oldCreator;
116    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
117    
118 <  latticeConstant = args_info.latticeCnst_arg;
118 >  latticeConstant = args_info.latticeConstant_arg;
119    particleRadius = args_info.radius_arg;
120    Globals* simParams = oldInfo->getSimParams();
121    
# Line 132 | Line 126 | int main(int argc, char *argv []) {
126    /* Build a lattice and get lattice points for this lattice constant */
127    vector<Vector3d> sites = nanoParticle.getSites();
128    vector<Vector3d> orientations = nanoParticle.getOrientations();
129 +  std::vector<int> vacancyTargets;
130 +  vector<bool> isVacancy;
131 +  
132 +  Vector3d myLoc;
133 +  RealType myR;
134 +
135 +  for (int i = 0; i < sites.size(); i++)
136 +    isVacancy.push_back(false);
137  
138 <  std::cout <<"nSites: " << sites.size() << std::endl;
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();
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 (nComponents == 1) {
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 >  if (nComponents == 1) {
215      molFractions.push_back(1.0);    
216 <  } else {
217 <   if (args_info.molFraction_given == nComponents) {
218 <     for (int i = 0; i < nComponents; i++) {
219 <       molFractions.push_back(args_info.molFraction_arg[i]);
220 <     }
221 <    } else if (args_info.molFraction_given == nComponents-1) {
222 <     RealType remainingFraction = 1.0;
223 <     for (int i = 0; i < nComponents-1; i++) {
224 <       molFractions.push_back(args_info.molFraction_arg[i]);
225 <       remainingFraction -= molFractions[i];
226 <     }
227 <     molFractions.push_back(remainingFraction);
228 <   } else {    
229 <     sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out molFractions "
230 <             "for all of the components in the <MetaData> block.");
231 <     painCave.isFatal = 1;
165 <     simError();
166 <   }
167 < }
168 <
169 < RealType totalFraction = 0.0;
170 <
171 < /* Do some simple sanity checking*/
172 <
173 <  for (int i = 0; i < nComponents; i++) {
174 <    if (molFractions.at(i) < 0.0) {
175 <      sprintf(painCave.errMsg, "One of the requested molFractions was"
176 <              " less than zero!");
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 <    if (molFractions.at(i) > 1.0) {
236 <      sprintf(painCave.errMsg, "One of the requested molFractions was"
237 <              " greater than one!");
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 <    totalFraction += molFractions.at(i);
252 <  }
253 <  if (abs(totalFraction - 1.0) > 1e-6) {
189 <    sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
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 <  int remaining = nSites;
259 <  for (int i=0; i < nComponents-1; i++) {    
260 <    nMol.push_back(int((RealType)nSites * molFractions.at(i)));
261 <    remaining -= nMol.at(i);
262 <  }
263 <  nMol.push_back(remaining);
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 <
308 <
309 < // recompute actual mol fractions and perform final sanity check:
310 <
311 <  int totalMolecules = 0;
312 <  RealType totalMass = 0.0;
313 <  for (int i=0; i < nComponents; i++) {
314 <    molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
315 <    totalMolecules += nMol.at(i);
316 <    molecularMasses.push_back(getMolMass(oldInfo->getMoleculeStamp(i),
317 <                                         oldInfo->getForceField()));
318 <    totalMass += (RealType)(nMol.at(i)) * molecularMasses.at(i);
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    }
214 –  RealType avgMass = totalMass / (RealType) totalMolecules;
320  
321 <  if (totalMolecules != nSites) {
322 <    sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
323 <            "to the number of lattice sites!");
324 <    painCave.isFatal = 1;
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 <  }
327 <    
328 <  vector<int> ids;
329 <  for (int i = 0; i < sites.size(); i++) ids.push_back(i);
330 <  /* Random particle is the default case*/
226 <  if (!args_info.ShellRadius_given){
227 <    /* do the iPod thing, Shuffle da vector */
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());
229 –  } else{ /*Handle core-shell with multiple components.*/
230 –    std::cout << "Creating a core-shell nanoparticle." << std::endl;
231 –    if (nComponents != args_info.ShellRadius_given + 1){
232 –      sprintf(painCave.errMsg, "Number of .md components "
233 –              "does not match the number of shell radius specifications");
234 –      painCave.isFatal = 1;
235 –      simError();
236 –    }  
332      
333 <  }
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    
241 –  
242 –  
362    outputFileName = args_info.output_arg;
363 <
364 <  
246 <   //creat new .md file on fly which corrects the number of molecule    
363 >  
364 >  //creat new .md file on fly which corrects the number of molecule    
365    createMdFile(inputFileName, outputFileName, nMol);
366    
367    if (oldInfo != NULL)
368      delete oldInfo;
369    
252 –  
253 –  // We need to read in new siminfo object.    
254 –  //parse md file and set up the system
255 –  //SimCreator NewCreator;
370    SimCreator newCreator;
371    SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
372 <  
259 <  
372 >    
373    // Place molecules
374    Molecule* mol;
375    SimInfo::MoleculeIterator mi;
376    mol = NewInfo->beginMolecule(mi);
377 +
378    int l = 0;
379 +  int whichSite = 0;
380  
381    for (int i = 0; i < nComponents; i++){
382      locator = new MoLocator(NewInfo->getMoleculeStamp(i),
383                              NewInfo->getForceField());
384 <    for (int n = 0; n < nMol.at(i); n++) {
385 <      mol = NewInfo->getMoleculeByGlobalIndex(l);
386 <      locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
387 <      l++;
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    
276 –
277 –  
404    //fill Hmat
405 <  hmat(0, 0)=  2.0*particleRadius;
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) =  2.0*particleRadius;
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) =  2.0*particleRadius;
415 >  hmat(2, 2) =  10.0*particleRadius;
416    
417    //set Hmat
418    NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
# Line 296 | Line 422 | int main(int argc, char *argv []) {
422    writer = new DumpWriter(NewInfo, outputFileName);
423    
424    if (writer == NULL) {
425 <    sprintf(painCave.errMsg, "Error in creating dumpwrite object ");
426 <      painCave.isFatal = 1;
427 <      simError();
425 >    sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
426 >    painCave.isFatal = 1;
427 >    simError();
428    }
429    
430    writer->writeDump();
# Line 315 | Line 441 | void createMdFile(const std::string&oldMdFileName, con
441    return 0;
442   }
443  
444 < void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
444 > void createMdFile(const std::string&oldMdFileName,
445 >                  const std::string&newMdFileName,
446                    std::vector<int> nMol) {
447    ifstream oldMdFile;
448    ofstream newMdFile;
# Line 325 | 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());
328 –  
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<nMol.size()){
463 <        sprintf(buffer, "\tnMol = %i;", nMol.at(i));                            
463 >        sprintf(buffer, "\tnMol = %i;", nMol.at(i));
464          newMdFile << buffer << std::endl;
465          i++;
466        }
# Line 346 | 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