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 1075 by gezelter, Tue Oct 17 17:51:52 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 85 | Line 85 | int main(int argc, char *argv []) {
85    double latticeConstant;
86    std::vector<double> lc;
87  
88 <  double particleRadius;
88 >  RealType particleRadius;
89  
90    Mat3x3d hmat;
91    std::vector<Vector3d> latticePos;
# Line 104 | 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"
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();
# Line 115 | 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 126 | 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.", (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();
197 >  int nSites = sites.size() - vacancyTargets.size();
198  
199    std::vector<Component*> components = simParams->getComponents();
200    std::vector<RealType> molFractions;
# Line 140 | Line 204 | int main(int argc, char *argv []) {
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 "
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();
# Line 168 | Line 232 | int main(int argc, char *argv []) {
232        painCave.isFatal = 1;
233        simError();
234      }
235 <  } else if ((int)args_info.ShellRadius_given) {
236 <    if ((int)args_info.ShellRadius_given == nComponents) {
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]);
238 >        shellRadii.push_back(args_info.shellRadius_arg[i]);
239        }
240 <    } else if ((int)args_info.ShellRadius_given == nComponents-1) {
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]);
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 shell radii "
247 <              "for all of the components in the <MetaData> block.");
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    } else {
252 <    sprintf(painCave.errMsg, "You have a multi-component <MetaData> block, but have not "
253 <            "specified either molFraction or ShellRadius arguments.");
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    }
# Line 253 | Line 317 | int main(int argc, char *argv []) {
317        }
318      }
319    }
320 <          
321 <  vector<int> ids;
258 <  for (int i = 0; i < sites.size(); i++) ids.push_back(i);
259 <  /* Random particle is the default case*/
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  
270 –    Vector3d myLoc;
271 –    RealType myR;
338      RealType smallestSoFar;
339      int myComponent = -1;
340      nMol.clear();
# Line 277 | Line 343 | int main(int argc, char *argv []) {
343      for (int i = 0; i < sites.size(); i++) {
344        myLoc = sites[i];
345        myR = myLoc.length();
346 <      smallestSoFar = particleRadius;
347 <    
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;
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 <      componentFromSite[i] = myComponent;
360 <      nMol[myComponent]++;
292 <    }
293 <  }  
359 >    }      
360 >  }
361    
362    outputFileName = args_info.output_arg;
363 <
297 <  
363 >  
364    //creat new .md file on fly which corrects the number of molecule    
365    createMdFile(inputFileName, outputFileName, nMol);
366    
# Line 308 | Line 374 | int main(int argc, char *argv []) {
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 <
385 <    if (args_info.ShellRadius_given) {
384 >    
385 >    if (!args_info.molFraction_given) {
386        for (int n = 0; n < sites.size(); n++) {
387 <        if (componentFromSite[n] == i) {
388 <          mol = NewInfo->getMoleculeByGlobalIndex(l);
389 <          locator->placeMol(sites[n], orientations[n], mol);
390 <          l++;
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 {
# Line 329 | Line 399 | int main(int argc, char *argv []) {
399          l++;
400        }
401      }
402 <  }
402 >  }
403    
404    //fill Hmat
405    hmat(0, 0)=  10.0*particleRadius;
# Line 352 | 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 ");
425 >    sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
426      painCave.isFatal = 1;
427      simError();
428    }
# Line 364 | Line 434 | int main(int argc, char *argv []) {
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();
# Line 382 | Line 452 | void createMdFile(const std::string&oldMdFileName,
452    //create new .md file based on old .md file
453    oldMdFile.open(oldMdFileName.c_str());
454    newMdFile.open(newMdFileName.c_str());
385 –  
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 403 | Line 472 | void createMdFile(const std::string&oldMdFileName,
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