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 1075 by gezelter, Tue Oct 17 17:51:52 2006 UTC

# 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  
95 –  const double rhoConvertConst = 1.661;
96 –  double density;
88    double particleRadius;
98 –  
99 –  
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    
134 –  
118    latticeConstant = args_info.latticeCnst_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();
127 >  vector<Vector3d> sites = nanoParticle.getSites();
128 >  vector<Vector3d> orientations = nanoParticle.getOrientations();
129  
130    std::cout <<"nSites: " << sites.size() << std::endl;
131  
132    /* Get number of lattice sites */
133    int nSites = sites.size();
154 –
134  
156 –
157 –
135    std::vector<Component*> components = simParams->getComponents();
136    std::vector<RealType> molFractions;
137 +  std::vector<RealType> shellRadii;
138    std::vector<RealType> molecularMasses;
139    std::vector<int> nMol;
140 +  std::map<int, int> componentFromSite;
141    nComponents = components.size();
163 –  
142  
143 <
144 < if (nComponents == 1) {
143 >  if (args_info.molFraction_given && args_info.ShellRadius_given) {
144 >    sprintf(painCave.errMsg, "Specify either molFraction or ShellRadius "
145 >            "arguments, but not both!");
146 >    painCave.isFatal = 1;
147 >    simError();
148 >  }
149 >  
150 >  if (nComponents == 1) {
151      molFractions.push_back(1.0);    
152 <  } else {
153 <    if (args_info.molFraction_given == nComponents) {
152 >    shellRadii.push_back(particleRadius);
153 >  } else if (args_info.molFraction_given) {
154 >    if ((int)args_info.molFraction_given == nComponents) {
155        for (int i = 0; i < nComponents; i++) {
156          molFractions.push_back(args_info.molFraction_arg[i]);
157        }
158 <    } else if (args_info.molFraction_given == nComponents-1) {
158 >    } else if ((int)args_info.molFraction_given == nComponents-1) {
159        RealType remainingFraction = 1.0;
160        for (int i = 0; i < nComponents-1; i++) {
161          molFractions.push_back(args_info.molFraction_arg[i]);
# Line 183 | Line 168 | int main(int argc, char *argv []) {
168        painCave.isFatal = 1;
169        simError();
170      }
171 +  } else if ((int)args_info.ShellRadius_given) {
172 +    if ((int)args_info.ShellRadius_given == nComponents) {
173 +      for (int i = 0; i < nComponents; i++) {
174 +        shellRadii.push_back(args_info.ShellRadius_arg[i]);
175 +      }
176 +    } else if ((int)args_info.ShellRadius_given == nComponents-1) {
177 +      for (int i = 0; i < nComponents-1; i++) {
178 +        shellRadii.push_back(args_info.ShellRadius_arg[i]);
179 +      }
180 +      shellRadii.push_back(particleRadius);
181 +    } else {    
182 +      sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out the shell radii "
183 +              "for all of the components in the <MetaData> block.");
184 +      painCave.isFatal = 1;
185 +      simError();
186 +    }
187 +  } else {
188 +    sprintf(painCave.errMsg, "You have a multi-component <MetaData> block, but have not "
189 +            "specified either molFraction or ShellRadius arguments.");
190 +    painCave.isFatal = 1;
191 +    simError();
192    }
193 <
194 < RealType totalFraction = 0.0;
195 <
196 < /* Do some simple sanity checking*/
197 <
198 <
199 <
200 <  for (int i = 0; i < nComponents; i++) {
201 <    if (molFractions.at(i) < 0.0) {
202 <      sprintf(painCave.errMsg, "One of the requested molFractions was"
203 <              " less than zero!");
193 >    
194 >  if (args_info.molFraction_given) {
195 >    RealType totalFraction = 0.0;
196 >    
197 >    /* Do some simple sanity checking*/
198 >    
199 >    for (int i = 0; i < nComponents; i++) {
200 >      if (molFractions.at(i) < 0.0) {
201 >        sprintf(painCave.errMsg, "One of the requested molFractions was"
202 >                " less than zero!");
203 >        painCave.isFatal = 1;
204 >        simError();
205 >      }
206 >      if (molFractions.at(i) > 1.0) {
207 >        sprintf(painCave.errMsg, "One of the requested molFractions was"
208 >                " greater than one!");
209 >        painCave.isFatal = 1;
210 >        simError();
211 >      }
212 >      totalFraction += molFractions.at(i);
213 >    }
214 >    if (abs(totalFraction - 1.0) > 1e-6) {
215 >      sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
216        painCave.isFatal = 1;
217        simError();
218      }
219 <    if (molFractions.at(i) > 1.0) {
220 <      sprintf(painCave.errMsg, "One of the requested molFractions was"
221 <              " greater than one!");
219 >    
220 >    int remaining = nSites;
221 >    for (int i=0; i < nComponents-1; i++) {    
222 >      nMol.push_back(int((RealType)nSites * molFractions.at(i)));
223 >      remaining -= nMol.at(i);
224 >    }
225 >    nMol.push_back(remaining);
226 >    
227 >    // recompute actual mol fractions and perform final sanity check:
228 >    
229 >    int totalMolecules = 0;
230 >    for (int i=0; i < nComponents; i++) {
231 >      molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
232 >      totalMolecules += nMol.at(i);
233 >    }
234 >    
235 >    if (totalMolecules != nSites) {
236 >      sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
237 >              "to the number of lattice sites!");
238        painCave.isFatal = 1;
239        simError();
240      }
241 <    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");
211 <    painCave.isFatal = 1;
212 <    simError();
213 <  }
214 <
215 <  int remaining = nSites;
216 <  for (int i=0; i < nComponents-1; i++) {    
217 <    nMol.push_back(int((RealType)nSites * molFractions.at(i)));
218 <    remaining -= nMol.at(i);
219 <  }
220 <  nMol.push_back(remaining);
241 >  } else {
242  
243 <
244 <
245 < // recompute actual mol fractions and perform final sanity check:
246 <
247 <  int totalMolecules = 0;
248 <  RealType totalMass = 0.0;
249 <  for (int i=0; i < nComponents; i++) {
250 <    molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
251 <    totalMolecules += nMol.at(i);
252 <    molecularMasses.push_back(getMolMass(oldInfo->getMoleculeStamp(i),
253 <                                         oldInfo->getForceField()));
254 <    totalMass += (RealType)(nMol.at(i)) * molecularMasses.at(i);
243 >    for (int i = 0; i < shellRadii.size(); i++) {
244 >      if (shellRadii.at(i) > particleRadius + 1e-6 ) {
245 >        sprintf(painCave.errMsg, "One of the shellRadius values exceeds the particle Radius.");
246 >        painCave.isFatal = 1;
247 >        simError();
248 >      }
249 >      if (shellRadii.at(i) <= 0.0 ) {
250 >        sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!");
251 >        painCave.isFatal = 1;
252 >        simError();
253 >      }
254 >    }
255    }
256 <  RealType avgMass = totalMass / (RealType) totalMolecules;
236 <
237 <  if (totalMolecules != nSites) {
238 <    sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
239 <            "to the number of lattice sites!");
240 <    painCave.isFatal = 1;
241 <    simError();
242 <  }
243 <    
244 <
245 <
246 <
256 >          
257    vector<int> ids;
258    for (int i = 0; i < sites.size(); i++) ids.push_back(i);
259    /* Random particle is the default case*/
260 <  if (!args_info.ShellRadius_given){
261 <    /* do the iPod thing, Shuffle da vector */
260 >  if ((int)args_info.molFraction_given){
261 >    sprintf(painCave.errMsg, "Creating a randomized spherical nanoparticle.");
262 >    painCave.isFatal = 0;
263 >    simError();
264      std::random_shuffle(ids.begin(), ids.end());
265 <  } else{ /*Handle core-shell with multiple components.*/
266 <    std::cout << "Creating a core-shell nanoparticle." << std::endl;
267 <    if (nComponents != args_info.ShellRadius_given + 1){
268 <      sprintf(painCave.errMsg, "Number of .md components "
257 <              "does not match the number of shell radius specifications");
258 <      painCave.isFatal = 1;
259 <      simError();
260 <    }  
261 <    
262 <  }
265 >  } else{
266 >    sprintf(painCave.errMsg, "Creating a core-shell spherical nanoparticle.");
267 >    painCave.isFatal = 0;
268 >    simError();
269  
270 +    Vector3d myLoc;
271 +    RealType myR;
272 +    RealType smallestSoFar;
273 +    int myComponent = -1;
274 +    nMol.clear();
275 +    nMol.resize(nComponents);
276 +
277 +    for (int i = 0; i < sites.size(); i++) {
278 +      myLoc = sites[i];
279 +      myR = myLoc.length();
280 +      smallestSoFar = particleRadius;
281 +    
282 +      for (int j = 0; j < nComponents; j++) {
283 +        if (myR <= shellRadii[j]) {
284 +          if (shellRadii[j] <= smallestSoFar) {
285 +            smallestSoFar = shellRadii[j];
286 +            myComponent = j;
287 +          }
288 +        }
289 +      }
290 +      componentFromSite[i] = myComponent;
291 +      nMol[myComponent]++;
292 +    }
293 +  }  
294    
265 –  
266 –  
295    outputFileName = args_info.output_arg;
296  
297    
298 <   //creat new .md file on fly which corrects the number of molecule    
298 >  //creat new .md file on fly which corrects the number of molecule    
299    createMdFile(inputFileName, outputFileName, nMol);
300    
301    if (oldInfo != NULL)
302      delete oldInfo;
303    
276 –  
277 –  // We need to read in new siminfo object.    
278 –  //parse md file and set up the system
279 –  //SimCreator NewCreator;
304    SimCreator newCreator;
305    SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
306 <  
283 <  
306 >    
307    // Place molecules
308    Molecule* mol;
309    SimInfo::MoleculeIterator mi;
# Line 290 | Line 313 | int main(int argc, char *argv []) {
313    for (int i = 0; i < nComponents; i++){
314      locator = new MoLocator(NewInfo->getMoleculeStamp(i),
315                              NewInfo->getForceField());
316 <    for (int n = 0; n < nMol.at(i); n++) {
317 <      mol = NewInfo->getMoleculeByGlobalIndex(l);
318 <      locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
319 <      l++;
316 >
317 >    if (args_info.ShellRadius_given) {
318 >      for (int n = 0; n < sites.size(); n++) {
319 >        if (componentFromSite[n] == i) {
320 >          mol = NewInfo->getMoleculeByGlobalIndex(l);
321 >          locator->placeMol(sites[n], orientations[n], mol);
322 >          l++;
323 >        }
324 >      }
325 >    } else {
326 >      for (int n = 0; n < nMol.at(i); n++) {
327 >        mol = NewInfo->getMoleculeByGlobalIndex(l);
328 >        locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
329 >        l++;
330 >      }
331      }
332 <  }
332 >  }
333    
300 –
301 –  
334    //fill Hmat
335 <  hmat(0, 0)=  2.0*particleRadius;
335 >  hmat(0, 0)=  10.0*particleRadius;
336    hmat(0, 1) = 0.0;
337    hmat(0, 2) = 0.0;
338    
339    hmat(1, 0) = 0.0;
340 <  hmat(1, 1) =  2.0*particleRadius;
340 >  hmat(1, 1) =  10.0*particleRadius;
341    hmat(1, 2) = 0.0;
342    
343    hmat(2, 0) = 0.0;
344    hmat(2, 1) = 0.0;
345 <  hmat(2, 2) =  2.0*particleRadius;
345 >  hmat(2, 2) =  10.0*particleRadius;
346    
347    //set Hmat
348    NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
349    
350    
351    //create dumpwriter and write out the coordinates
352 <  writer = new DumpWriter(NewInfo,outputFileName);
352 >  writer = new DumpWriter(NewInfo, outputFileName);
353    
354    if (writer == NULL) {
355      sprintf(painCave.errMsg, "Error in creating dumpwrite object ");
356 <      painCave.isFatal = 1;
357 <      simError();
356 >    painCave.isFatal = 1;
357 >    simError();
358    }
359    
360    writer->writeDump();
329 –  std::cout << "new initial configuration file: " << outInitFileName
330 –            << " is generated." << std::endl;
331 –
361  
362    // deleting the writer will put the closing at the end of the dump file
363 +
364    delete writer;
365  
366    // cleanup a by calling sim error.....
# Line 341 | Line 371 | int main(int argc, char *argv []) {
371    return 0;
372   }
373  
374 < void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
374 > void createMdFile(const std::string&oldMdFileName,
375 >                  const std::string&newMdFileName,
376                    std::vector<int> nMol) {
377    ifstream oldMdFile;
378    ofstream newMdFile;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines