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

Comparing trunk/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp (file contents):
Revision 1069 by chuckv, Fri Oct 13 20:16:59 2006 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   #include <cstdlib>
# Line 64 | Line 64 | using namespace std;
64   #include "utils/StringUtils.hpp"
65  
66   using namespace std;
67 < using namespace oopse;
67 > using namespace OpenMD;
68   void createMdFile(const std::string&oldMdFileName,
69                    const std::string&newMdFileName,
70                    std::vector<int> numMol);
# 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;
81    std::string outputFileName;
83 –  std::string outInitFileName;
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;                                                                          
87  
88 <  const double rhoConvertConst = 1.661;
96 <  double density;
97 <  double particleRadius;
98 <  
99 <  
88 >  RealType particleRadius;
89  
90    Mat3x3d hmat;
91    std::vector<Vector3d> latticePos;
92    std::vector<Vector3d> latticeOrt;
104 –  int numMolPerCell;
105 –  int nShells; /* Number of shells in nanoparticle*/
93  
107 –  
94    DumpWriter *writer;
95    
96    // Parse Command Line Arguments
97    if (cmdline_parser(argc, argv, &args_info) != 0)
98      exit(1);
99 <  
114 <        
115 <        
99 >        
100    /* get lattice type */
101    latticeType = "FCC";
102  
# Line 120 | 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 131 | Line 115 | int main(int argc, char *argv []) {
115    SimCreator oldCreator;
116    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
117    
118 <  
135 <  latticeConstant = args_info.latticeCnst_arg;
118 >  latticeConstant = args_info.latticeConstant_arg;
119    particleRadius = args_info.radius_arg;
120    Globals* simParams = oldInfo->getSimParams();
121    
139 –      
140 –  /* create Molocators */
141 –  locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
142 –  
122    /* Create nanoparticle */
123 <  shapedLatticeSpherical nanoParticle(latticeConstant,latticeType,particleRadius);
123 >  shapedLatticeSpherical nanoParticle(latticeConstant, latticeType,
124 >                                      particleRadius);
125    
126    /* Build a lattice and get lattice points for this lattice constant */
127 <  vector<Vector3d> sites = nanoParticle.getPoints();
128 <  vector<Vector3d> orientations = nanoParticle.getPointsOrt();
129 <
130 <  std::cout <<"nSites: " << sites.size() << std::endl;
131 <
132 <  /* Get number of lattice sites */
133 <  int nSites = sites.size();
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 +  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.", (int) 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();
163 –  
206  
207 <
208 < 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) {
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 (args_info.molFraction_given == nComponents-1) {
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]);
# Line 183 | Line 232 | int main(int argc, char *argv []) {
232        painCave.isFatal = 1;
233        simError();
234      }
235 <  }
236 <
237 < RealType totalFraction = 0.0;
238 <
239 < /* Do some simple sanity checking*/
240 <
241 <
242 <
243 <  for (int i = 0; i < nComponents; i++) {
244 <    if (molFractions.at(i) < 0.0) {
245 <      sprintf(painCave.errMsg, "One of the requested molFractions was"
246 <              " less than zero!");
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 <    if (molFractions.at(i) > 1.0) {
252 <      sprintf(painCave.errMsg, "One of the requested molFractions was"
253 <              " greater than one!");
204 <      painCave.isFatal = 1;
205 <      simError();
206 <    }
207 <    totalFraction += molFractions.at(i);
208 <  }
209 <  if (abs(totalFraction - 1.0) > 1e-6) {
210 <    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    }
235 –  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 <  }
243 <    
326 >    /* Random particle is the default case*/
327  
328 <
329 <
330 <  vector<int> ids;
248 <  for (int i = 0; i < sites.size(); i++) ids.push_back(i);
249 <  /* Random particle is the default case*/
250 <  if (!args_info.ShellRadius_given){
251 <    /* do the iPod thing, Shuffle da vector */
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());
253 –  } else{ /*Handle core-shell with multiple components.*/
254 –    std::cout << "Creating a core-shell nanoparticle." << std::endl;
255 –    if (nComponents != args_info.ShellRadius_given + 1){
256 –      sprintf(painCave.errMsg, "Number of .md components "
257 –              "does not match the number of shell radius specifications");
258 –      painCave.isFatal = 1;
259 –      simError();
260 –    }  
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    
265 –  
266 –  
362    outputFileName = args_info.output_arg;
363 <
364 <  
270 <   //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    
276 –  
277 –  // We need to read in new siminfo object.    
278 –  //parse md file and set up the system
279 –  //SimCreator NewCreator;
370    SimCreator newCreator;
371    SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
372 <  
283 <  
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    
300 –
301 –  
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);
419    
420    
421    //create dumpwriter and write out the coordinates
422 <  writer = new DumpWriter(NewInfo,outputFileName);
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();
329 –  std::cout << "new initial configuration file: " << outInitFileName
330 –            << " is generated." << std::endl;
331 –
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 "
437 >  sprintf(painCave.errMsg, "A new OpenMD 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,
444 > void createMdFile(const std::string&oldMdFileName,
445 >                  const std::string&newMdFileName,
446                    std::vector<int> nMol) {
447    ifstream oldMdFile;
448    ofstream newMdFile;
# Line 351 | 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());
354 –  
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 372 | 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