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 1782 by gezelter, Wed Aug 22 02:28:28 2012 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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <cstdlib>
# Line 64 | Line 65 | using namespace std;
65   #include "utils/StringUtils.hpp"
66  
67   using namespace std;
68 < using namespace oopse;
68 > using namespace OpenMD;
69   void createMdFile(const std::string&oldMdFileName,
70                    const std::string&newMdFileName,
71                    std::vector<int> numMol);
72  
73   int main(int argc, char *argv []) {
74    
74  //register force fields
75  registerForceFields();
75    registerLattice();
76    
77    gengetopt_args_info args_info;
# Line 85 | Line 84 | int main(int argc, char *argv []) {
84    double latticeConstant;
85    std::vector<double> lc;
86  
87 <  double particleRadius;
87 >  RealType particleRadius;
88  
89    Mat3x3d hmat;
90    std::vector<Vector3d> latticePos;
# Line 104 | Line 103 | int main(int argc, char *argv []) {
103    if (args_info.inputs_num)
104      inputFileName = args_info.inputs[0];
105    else {
106 <    sprintf(painCave.errMsg, "No input .md file name was specified"
106 >    sprintf(painCave.errMsg, "No input .md file name was specified "
107              "on the command line");
108      painCave.isFatal = 1;
109      cmdline_parser_print_help();
# Line 115 | Line 114 | int main(int argc, char *argv []) {
114    SimCreator oldCreator;
115    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
116    
117 <  latticeConstant = args_info.latticeCnst_arg;
117 >  latticeConstant = args_info.latticeConstant_arg;
118    particleRadius = args_info.radius_arg;
119    Globals* simParams = oldInfo->getSimParams();
120    
# Line 127 | Line 126 | int main(int argc, char *argv []) {
126    vector<Vector3d> sites = nanoParticle.getSites();
127    vector<Vector3d> orientations = nanoParticle.getOrientations();
128  
130  std::cout <<"nSites: " << sites.size() << std::endl;
129  
130 +  std::vector<int> vacancyTargets;
131 +  vector<bool> isVacancy;
132 +  
133 +  Vector3d myLoc;
134 +  RealType myR;
135 +
136 +  for (int i = 0; i < sites.size(); i++)
137 +    isVacancy.push_back(false);
138 +
139 +  if (args_info.vacancyPercent_given) {
140 +    if (args_info.vacancyPercent_arg < 0.0 || args_info.vacancyPercent_arg > 100.0) {
141 +      sprintf(painCave.errMsg, "vacancyPercent was set to a non-sensical value.");
142 +      painCave.isFatal = 1;
143 +      simError();
144 +    } else {
145 +      RealType vF = args_info.vacancyPercent_arg / 100.0;
146 +      RealType vIR;
147 +      RealType vOR;
148 +      if (args_info.vacancyInnerRadius_given) {
149 +        vIR = args_info.vacancyInnerRadius_arg;
150 +      } else {
151 +        vIR = 0.0;
152 +      }
153 +      if (args_info.vacancyOuterRadius_given) {
154 +        vOR = args_info.vacancyOuterRadius_arg;
155 +      } else {
156 +        vOR = particleRadius;
157 +      }
158 +      if (vIR >= 0.0 && vOR <= particleRadius && vOR >= vIR) {
159 +        
160 +        for (int i = 0; i < sites.size(); i++) {
161 +          myLoc = sites[i];
162 +          myR = myLoc.length();
163 +          if (myR >= vIR && myR <= vOR) {
164 +            vacancyTargets.push_back(i);
165 +          }          
166 +        }
167 +        std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
168 +        
169 +        int nTargets = vacancyTargets.size();
170 +        vacancyTargets.resize((int)(vF * nTargets));
171 +        
172 +                  
173 +        sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
174 +                "\tsites between %lf and %lf.", (int) vacancyTargets.size(),
175 +                vIR, vOR);
176 +        painCave.isFatal = 0;
177 +        simError();
178 +
179 +        isVacancy.clear();
180 +        for (int i = 0; i < sites.size(); i++) {
181 +          bool vac = false;
182 +          for (int j = 0; j < vacancyTargets.size(); j++) {
183 +            if (i == vacancyTargets[j]) vac = true;
184 +          }
185 +          isVacancy.push_back(vac);
186 +        }
187 +              
188 +      } else {
189 +        sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
190 +                "\tinner or outer radii.  Check their values.");
191 +        painCave.isFatal = 1;
192 +        simError();
193 +      }
194 +    }
195 +  }
196 +
197    /* Get number of lattice sites */
198 <  int nSites = sites.size();
198 >  int nSites = sites.size() - vacancyTargets.size();
199  
200    std::vector<Component*> components = simParams->getComponents();
201    std::vector<RealType> molFractions;
# Line 140 | Line 205 | int main(int argc, char *argv []) {
205    std::map<int, int> componentFromSite;
206    nComponents = components.size();
207  
208 <  if (args_info.molFraction_given && args_info.ShellRadius_given) {
209 <    sprintf(painCave.errMsg, "Specify either molFraction or ShellRadius "
208 >  if (args_info.molFraction_given && args_info.shellRadius_given) {
209 >    sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
210              "arguments, but not both!");
211      painCave.isFatal = 1;
212      simError();
# Line 168 | Line 233 | int main(int argc, char *argv []) {
233        painCave.isFatal = 1;
234        simError();
235      }
236 <  } else if ((int)args_info.ShellRadius_given) {
237 <    if ((int)args_info.ShellRadius_given == nComponents) {
236 >  } else if ((int)args_info.shellRadius_given) {
237 >    if ((int)args_info.shellRadius_given == nComponents) {
238        for (int i = 0; i < nComponents; i++) {
239 <        shellRadii.push_back(args_info.ShellRadius_arg[i]);
239 >        shellRadii.push_back(args_info.shellRadius_arg[i]);
240        }
241 <    } else if ((int)args_info.ShellRadius_given == nComponents-1) {
241 >    } else if ((int)args_info.shellRadius_given == nComponents-1) {
242        for (int i = 0; i < nComponents-1; i++) {
243 <        shellRadii.push_back(args_info.ShellRadius_arg[i]);
243 >        shellRadii.push_back(args_info.shellRadius_arg[i]);
244        }
245        shellRadii.push_back(particleRadius);
246      } else {    
247 <      sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out the shell radii "
248 <              "for all of the components in the <MetaData> block.");
247 >      sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out the\n"
248 >              "\tshell radii for all of the components in the <MetaData> block.");
249        painCave.isFatal = 1;
250        simError();
251      }
252    } else {
253 <    sprintf(painCave.errMsg, "You have a multi-component <MetaData> block, but have not "
254 <            "specified either molFraction or ShellRadius arguments.");
253 >    sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
254 >            "\tbut have not specified either molFraction or shellRadius arguments.");
255      painCave.isFatal = 1;
256      simError();
257    }
# Line 253 | Line 318 | int main(int argc, char *argv []) {
318        }
319      }
320    }
321 <          
322 <  vector<int> ids;
258 <  for (int i = 0; i < sites.size(); i++) ids.push_back(i);
259 <  /* Random particle is the default case*/
321 >
322 >  vector<int> ids;          
323    if ((int)args_info.molFraction_given){
324      sprintf(painCave.errMsg, "Creating a randomized spherical nanoparticle.");
325      painCave.isFatal = 0;
326      simError();
327 +    /* Random particle is the default case*/
328 +
329 +    for (int i = 0; i < sites.size(); i++)
330 +      if (!isVacancy[i]) ids.push_back(i);
331 +    
332      std::random_shuffle(ids.begin(), ids.end());
333 +    
334    } else{
335      sprintf(painCave.errMsg, "Creating a core-shell spherical nanoparticle.");
336      painCave.isFatal = 0;
337      simError();
338  
270    Vector3d myLoc;
271    RealType myR;
339      RealType smallestSoFar;
340      int myComponent = -1;
341      nMol.clear();
# Line 277 | Line 344 | int main(int argc, char *argv []) {
344      for (int i = 0; i < sites.size(); i++) {
345        myLoc = sites[i];
346        myR = myLoc.length();
347 <      smallestSoFar = particleRadius;
348 <    
349 <      for (int j = 0; j < nComponents; j++) {
350 <        if (myR <= shellRadii[j]) {
351 <          if (shellRadii[j] <= smallestSoFar) {
352 <            smallestSoFar = shellRadii[j];
353 <            myComponent = j;
347 >      smallestSoFar = particleRadius;      
348 >      if (!isVacancy[i]) {
349 >        for (int j = 0; j < nComponents; j++) {
350 >          if (myR <= shellRadii[j]) {
351 >            if (shellRadii[j] <= smallestSoFar) {
352 >              smallestSoFar = shellRadii[j];
353 >              myComponent = j;
354 >            }
355            }
356          }
357 +        componentFromSite[i] = myComponent;
358 +        nMol[myComponent]++;
359        }
360 <      componentFromSite[i] = myComponent;
361 <      nMol[myComponent]++;
292 <    }
293 <  }  
360 >    }      
361 >  }
362    
363    outputFileName = args_info.output_arg;
364 <
297 <  
364 >  
365    //creat new .md file on fly which corrects the number of molecule    
366    createMdFile(inputFileName, outputFileName, nMol);
367    
# Line 308 | Line 375 | int main(int argc, char *argv []) {
375    Molecule* mol;
376    SimInfo::MoleculeIterator mi;
377    mol = NewInfo->beginMolecule(mi);
378 +
379    int l = 0;
380 +  int whichSite = 0;
381  
382    for (int i = 0; i < nComponents; i++){
383      locator = new MoLocator(NewInfo->getMoleculeStamp(i),
384                              NewInfo->getForceField());
385 <
386 <    if (args_info.ShellRadius_given) {
385 >    
386 >    if (!args_info.molFraction_given) {
387        for (int n = 0; n < sites.size(); n++) {
388 <        if (componentFromSite[n] == i) {
389 <          mol = NewInfo->getMoleculeByGlobalIndex(l);
390 <          locator->placeMol(sites[n], orientations[n], mol);
391 <          l++;
388 >        if (!isVacancy[n]) {
389 >          if (componentFromSite[n] == i) {
390 >            mol = NewInfo->getMoleculeByGlobalIndex(l);
391 >            locator->placeMol(sites[n], orientations[n], mol);
392 >            l++;
393 >          }
394          }
395        }
396      } else {
# Line 329 | Line 400 | int main(int argc, char *argv []) {
400          l++;
401        }
402      }
403 <  }
403 >  }
404    
405    //fill Hmat
406    hmat(0, 0)=  10.0*particleRadius;
# Line 352 | Line 423 | int main(int argc, char *argv []) {
423    writer = new DumpWriter(NewInfo, outputFileName);
424    
425    if (writer == NULL) {
426 <    sprintf(painCave.errMsg, "Error in creating dumpwrite object ");
426 >    sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
427      painCave.isFatal = 1;
428      simError();
429    }
# Line 364 | Line 435 | int main(int argc, char *argv []) {
435    delete writer;
436  
437    // cleanup a by calling sim error.....
438 <  sprintf(painCave.errMsg, "A new OOPSE MD file called \"%s\" has been "
438 >  sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
439            "generated.\n", outputFileName.c_str());
440    painCave.isFatal = 0;
441    simError();
# Line 382 | Line 453 | void createMdFile(const std::string&oldMdFileName,
453    //create new .md file based on old .md file
454    oldMdFile.open(oldMdFileName.c_str());
455    newMdFile.open(newMdFileName.c_str());
385  
456    oldMdFile.getline(buffer, MAXLEN);
457 <
457 >
458    int i = 0;
459    while (!oldMdFile.eof()) {
460 <    
460 >
461      //correct molecule number
462      if (strstr(buffer, "nMol") != NULL) {
463        if(i<nMol.size()){
464 <        sprintf(buffer, "\tnMol = %i;", nMol.at(i));                            
464 >        sprintf(buffer, "\tnMol = %i;", nMol.at(i));
465          newMdFile << buffer << std::endl;
466          i++;
467        }
# Line 403 | Line 473 | void createMdFile(const std::string&oldMdFileName,
473    
474    oldMdFile.close();
475    newMdFile.close();
476 +
477 +  if (i != nMol.size()) {
478 +    sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
479 +            "\tstatements in component blocks.  Make sure that all\n"
480 +            "\tcomponents in the template file have nMol=1");
481 +    painCave.isFatal = 1;
482 +    simError();
483 +  }
484 +    
485   }
486  

Comparing trunk/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp (property svn:keywords):
Revision 1075 by gezelter, Tue Oct 17 17:51:52 2006 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines