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 653 by chuckv, Tue Oct 11 21:57:22 2005 UTC vs.
Revision 1077 by gezelter, Wed Oct 18 19:35:07 2006 UTC

# Line 47 | Line 47
47   #include <string>
48   #include <map>
49   #include <fstream>
50 + #include <algorithm>
51  
52   #include "config.h"
53 <
53 > #include "shapedLatticeSpherical.hpp"
54   #include "nanoparticleBuilderCmd.h"
54 – #include "sphericalNanoparticle.hpp"
55   #include "lattice/LatticeFactory.hpp"
56   #include "utils/MoLocator.hpp"
57   #include "lattice/Lattice.hpp"
# 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 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  
83 <  
84 <  
87 <  Lattice *simpleLat;
88 <  int numMol;
83 >  MoLocator* locator;
84 >  int nComponents;
85    double latticeConstant;
86    std::vector<double> lc;
91 –  double mass;
92 –  const double rhoConvertConst = 1.661;
93 –  double density;
94 –  
95 –  
87  
88 +  RealType particleRadius;
89 +
90    Mat3x3d hmat;
98 –  MoLocator *locator;
99 –  sphericalNanoparticle *nanoparticle;
91    std::vector<Vector3d> latticePos;
92    std::vector<Vector3d> latticeOrt;
93 <  int numMolPerCell;
103 <  int nShells; /* Number of shells in nanoparticle*/
104 <  int numSites;
105 <  
93 >
94    DumpWriter *writer;
95    
96 <  // parse command line arguments
96 >  // Parse Command Line Arguments
97    if (cmdline_parser(argc, argv, &args_info) != 0)
98      exit(1);
99 <  
112 <        
113 <        
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 <  nShells = 0;
119 <  if (args_info.coreShellRadius_given){
120 <    nShells = args_info.coreShellRadius_given;
133 <  }
118 >  latticeConstant = args_info.latticeConstant_arg;
119 >  particleRadius = args_info.radius_arg;
120 >  Globals* simParams = oldInfo->getSimParams();
121    
122 <  nComponents = oldInfo->getNMoleculeStamp();
122 >  /* Create nanoparticle */
123 >  shapedLatticeSpherical nanoParticle(latticeConstant, latticeType,
124 >                                      particleRadius);
125    
126 <  /* Check to see if we have enough components to build that many shells. */
127 <  if (nShells){
128 <    if (oldInfo->getNMoleculeStamp() != nShells) {
129 <      std::cerr << "Not enough components present in MD file to build specified number of shells"
130 <      << std::endl;
131 <      exit(1);
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 >  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 <  
197 <  //creat lattice
198 <  simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
199 <  
200 <  if (simpleLat == NULL) {
201 <    std::cerr << "Error in creating lattice" << std::endl;
202 <    exit(1);
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 <  numMolPerCell = simpleLat->getNumSitesPerCell();
215 <  
216 <  /*calculate lattice constant (in Angstrom)
217 <  latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
218 <                        1.0 / 3.0);*/
219 <  
220 <  latticeConstant = args_info.latticeCnst_arg;
221 <  particleRadius = args_info.radius_arg;
222 <  particleDiameter = 2.0 * particleRadius;
223 <  
224 <  /* set lattice constant */
225 <  lc.push_back(latticeConstant);
226 <  simpleLat->setLatticeConstant(lc);
227 <  
228 <  
229 <  /*determine the output file names*/
230 <  if (args_info.output_given)
231 <    outInitFileName = args_info.output_arg;
232 <  else
233 <    outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
175 <  
176 <  
177 <        
178 <  
179 <  
180 <  
181 <  /* create Molocators */
182 <  locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
183 <  
184 <  /* create a new spherical nanoparticle */
185 <  nanoparticle = new sphericalNanoparticle(particleRadius,latticeConstant);
186 <  /* Build a nanoparticle to see how many sites are there */
187 <  numSites = new int[nComponents]
188 <  nanoparticle.getNMol(numSites);
189 <  
190 <  numMol = new int[nComponents];
191 <  /* Random particle is the default case*/
192 <  if (!args_info.ShellRadius_given){
193 <    std::cout << "Creating a random nanoparticle" << std::endl;
194 <    /* Check to see if we have enough components */
195 <    if (nComponents != args_info.molFraction_given + 1){
196 <      std::cerr << "Number of components does not equal molFraction occurances." << std::endl;
197 <      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 <    int totComponents = 0;
236 <    for (int i = 0;i<nComponents-2;i++){ /* Figure out Percent for each component */
237 <      numMol[i] = int((double)numSites * args_info.molFraction_arg[i]);
238 <      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 <  } 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 <      std::cerr << "Number of components does not equal ShellRadius occurances." << std::endl;
257 <      exit 1;
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 <   //get the orientation of the cell sites
322 <  //for the same type of molecule in same lattice, it will not change
323 <   latticeOrt = simpleLat->getLatticePointsOrt();
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 +  outputFileName = args_info.output_arg;
363 +  
364 +  //creat new .md file on fly which corrects the number of molecule    
365 +  createMdFile(inputFileName, outputFileName, nMol);
366    
222 –  
223 –  // needed for writing out new md file.
224 –  
225 –    outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
226 –    outMdFileName = outPrefix + ".md";
227 –  
228 –    //creat new .md file on fly which corrects the number of molecule    
229 –    createMdFile(inputFileName, outMdFileName, numcomponents,numMol);
230 –  
367    if (oldInfo != NULL)
368      delete oldInfo;
369    
370 <  
371 <  // We need to read in new siminfo object.    
372 <  //parse md file and set up the system
373 <  //SimCreator NewCreator;
238 <  
239 <  SimInfo* NewInfo = oldCreator.createSim(outMdFileName, false);
240 <  
241 <  // This was so much fun the first time, lets do it again.
242 <  
370 >  SimCreator newCreator;
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 +  int whichSite = 0;
380  
381 <  for(int i = -nx; i < nx; i++) {
382 <     for(int j = -ny; j < ny; j++) {
383 <        for(int k = -nz; k < nz; k++) {
384 <          
385 <           //get the position of the cell sites
386 <           simpleLat->getLatticePointsPos(latticePos, i, j, k);
387 <          
388 <           for(int l = 0; l < numMolPerCell; l++) {
389 < #ifdef HAVE_CGAL              
390 <              if (myGeometry->isInsidePolyhedron(latticePos[l][0],latticePos[l][1],latticePos[l][2])){
391 < #endif                              
392 <                 if (mol != NULL) {
260 <                    locator->placeMol(latticePos[l], latticeOrt[l], mol);
261 <                 } else {
262 <                    std::cerr<<"Error in placing molecule " << std::endl;                    
263 <                 }
264 <                 mol = NewInfo->nextMolecule(mi);
265 < #ifdef HAVE_CGAL                
266 <              }
267 < #endif              
268 <           }
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 <     }
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    
273 –
274 –  
404    //fill Hmat
405 <  hmat(0, 0)= nx * 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) = ny * 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) = nz * 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);
294 <  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 <  
312 <  if (writer != NULL)
313 <    delete writer;
314 <  delete simpleLat;    
315 <  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 326 | 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());
329 –  
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 <      sprintf(buffer, "\tnMol = %i;", numMol);                          
463 <      newMdFile << buffer << std::endl;
462 >      if(i<nMol.size()){
463 >        sprintf(buffer, "\tnMol = %i;", nMol.at(i));
464 >        newMdFile << buffer << std::endl;
465 >        i++;
466 >      }
467      } else
468        newMdFile << buffer << std::endl;
469      
# Line 343 | 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