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

Comparing trunk/OOPSE-4/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp (file contents):
Revision 3050 by gezelter, Sat Oct 14 20:21:26 2006 UTC vs.
Revision 3051 by gezelter, Tue Oct 17 17:51:52 2006 UTC

# Line 80 | Line 80 | int main(int argc, char *argv []) {
80    std::string inputFileName;
81    std::string outputFileName;
82  
83 –  Lattice *simpleLat;
83    MoLocator* locator;
85 –  int* numMol;
84    int nComponents;
85    double latticeConstant;
86    std::vector<double> lc;
87 <  double mass;                                                                          
90 <  const double rhoConvertConst = 1.661;
91 <  double density;
87 >
88    double particleRadius;
89  
90    Mat3x3d hmat;
91    std::vector<Vector3d> latticePos;
92    std::vector<Vector3d> latticeOrt;
97 –  int numMolPerCell;
98 –  int nShells; /* Number of shells in nanoparticle*/
93  
94    DumpWriter *writer;
95    
# Line 110 | 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 140 | Line 134 | int main(int argc, char *argv []) {
134  
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();
142  
143 < 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) {
154 <     for (int i = 0; i < nComponents; i++) {
155 <       molFractions.push_back(args_info.molFraction_arg[i]);
156 <     }
157 <    } else if (args_info.molFraction_given == nComponents-1) {
158 <     RealType remainingFraction = 1.0;
159 <     for (int i = 0; i < nComponents-1; i++) {
160 <       molFractions.push_back(args_info.molFraction_arg[i]);
161 <       remainingFraction -= molFractions[i];
162 <     }
163 <     molFractions.push_back(remainingFraction);
164 <   } else {    
165 <     sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out molFractions "
166 <             "for all of the components in the <MetaData> block.");
167 <     painCave.isFatal = 1;
165 <     simError();
166 <   }
167 < }
168 <
169 < RealType totalFraction = 0.0;
170 <
171 < /* Do some simple sanity checking*/
172 <
173 <  for (int i = 0; i < nComponents; i++) {
174 <    if (molFractions.at(i) < 0.0) {
175 <      sprintf(painCave.errMsg, "One of the requested molFractions was"
176 <              " less than zero!");
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 ((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]);
162 >        remainingFraction -= molFractions[i];
163 >      }
164 >      molFractions.push_back(remainingFraction);
165 >    } else {    
166 >      sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out molFractions "
167 >              "for all of the components in the <MetaData> block.");
168        painCave.isFatal = 1;
169        simError();
170      }
171 <    if (molFractions.at(i) > 1.0) {
172 <      sprintf(painCave.errMsg, "One of the requested molFractions was"
173 <              " greater than one!");
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 <    totalFraction += molFractions.at(i);
188 <  }
189 <  if (abs(totalFraction - 1.0) > 1e-6) {
189 <    sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
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 <  int remaining = nSites;
195 <  for (int i=0; i < nComponents-1; i++) {    
196 <    nMol.push_back(int((RealType)nSites * molFractions.at(i)));
197 <    remaining -= nMol.at(i);
198 <  }
199 <  nMol.push_back(remaining);
200 <
201 <
202 <
203 < // recompute actual mol fractions and perform final sanity check:
204 <
205 <  int totalMolecules = 0;
206 <  RealType totalMass = 0.0;
207 <  for (int i=0; i < nComponents; i++) {
208 <    molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
209 <    totalMolecules += nMol.at(i);
210 <    molecularMasses.push_back(getMolMass(oldInfo->getMoleculeStamp(i),
211 <                                         oldInfo->getForceField()));
212 <    totalMass += (RealType)(nMol.at(i)) * molecularMasses.at(i);
213 <  }
214 <  RealType avgMass = totalMass / (RealType) totalMolecules;
215 <
216 <  if (totalMolecules != nSites) {
217 <    sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
218 <            "to the number of lattice sites!");
219 <    painCave.isFatal = 1;
220 <    simError();
221 <  }
222 <    
223 <  vector<int> ids;
224 <  for (int i = 0; i < sites.size(); i++) ids.push_back(i);
225 <  /* Random particle is the default case*/
226 <  if (!args_info.ShellRadius_given){
227 <    /* do the iPod thing, Shuffle da vector */
228 <    std::random_shuffle(ids.begin(), ids.end());
229 <  } else{ /*Handle core-shell with multiple components.*/
230 <    std::cout << "Creating a core-shell nanoparticle." << std::endl;
231 <    if (nComponents != args_info.ShellRadius_given + 1){
232 <      sprintf(painCave.errMsg, "Number of .md components "
233 <              "does not match the number of shell radius specifications");
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 <    }  
218 >    }
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 +  } else {
242 +
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 +          
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 ((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{
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    
241 –  
242 –  
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    
252 –  
253 –  // We need to read in new siminfo object.    
254 –  //parse md file and set up the system
255 –  //SimCreator NewCreator;
304    SimCreator newCreator;
305    SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
306 <  
259 <  
306 >    
307    // Place molecules
308    Molecule* mol;
309    SimInfo::MoleculeIterator mi;
# Line 266 | 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    
276 –
277 –  
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);
# Line 297 | Line 353 | int main(int argc, char *argv []) {
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();
# Line 315 | Line 371 | void createMdFile(const std::string&oldMdFileName, con
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