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 2352 by chuckv, Tue Oct 11 21:57:22 2005 UTC vs.
Revision 3051 by gezelter, Tue Oct 17 17:51:52 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 +  double 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    
130  nShells = 0;
131  if (args_info.coreShellRadius_given){
132    nShells = args_info.coreShellRadius_given;
133  }
134  
135  nComponents = oldInfo->getNMoleculeStamp();
136  
137  /* Check to see if we have enough components to build that many shells. */
138  if (nShells){
139    if (oldInfo->getNMoleculeStamp() != nShells) {
140      std::cerr << "Not enough components present in MD file to build specified number of shells"
141      << std::endl;
142      exit(1);
143    }
144  }
145  
146  
147  //creat lattice
148  simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
149  
150  if (simpleLat == NULL) {
151    std::cerr << "Error in creating lattice" << std::endl;
152    exit(1);
153  }
154  
155  numMolPerCell = simpleLat->getNumSitesPerCell();
156  
157  /*calculate lattice constant (in Angstrom)
158  latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
159                        1.0 / 3.0);*/
160  
118    latticeConstant = args_info.latticeCnst_arg;
119    particleRadius = args_info.radius_arg;
120 <  particleDiameter = 2.0 * particleRadius;
120 >  Globals* simParams = oldInfo->getSimParams();
121    
122 <  /* set lattice constant */
123 <  lc.push_back(latticeConstant);
124 <  simpleLat->setLatticeConstant(lc);
122 >  /* Create nanoparticle */
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.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();
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 (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 <  /*determine the output file names*/
151 <  if (args_info.output_given)
152 <    outInitFileName = args_info.output_arg;
153 <  else
154 <    outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
155 <  
156 <  
157 <        
158 <  
159 <  
160 <  
161 <  /* create Molocators */
162 <  locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
163 <  
164 <  /* create a new spherical nanoparticle */
165 <  nanoparticle = new sphericalNanoparticle(particleRadius,latticeConstant);
166 <  /* Build a nanoparticle to see how many sites are there */
167 <  numSites = new int[nComponents]
168 <  nanoparticle.getNMol(numSites);
169 <  
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;
150 >  if (nComponents == 1) {
151 >    molFractions.push_back(1.0);    
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 <    int totComponents = 0;
172 <    for (int i = 0;i<nComponents-2;i++){ /* Figure out Percent for each component */
173 <      numMol[i] = int((double)numSites * args_info.molFraction_arg[i]);
174 <      totComponents += numMol[i];
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 <    numMol[nComponents-1] = numSites - totComponents;
188 <
189 <  } else{ /*Handle core-shell with multiple components.*/
190 <    std::cout << "Creating a core-shell nanoparticle." << std::endl;
191 <    if (nComponents != args_info.ShellRadius_given + 1){
192 <      std::cerr << "Number of components does not equal ShellRadius occurances." << std::endl;
193 <      exit 1;
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 >  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      
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 <   //get the orientation of the cell sites
271 <  //for the same type of molecule in same lattice, it will not change
272 <   latticeOrt = simpleLat->getLatticePointsOrt();
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    
295 +  outputFileName = args_info.output_arg;
296 +
297    
298 +  //creat new .md file on fly which corrects the number of molecule    
299 +  createMdFile(inputFileName, outputFileName, nMol);
300    
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  
301    if (oldInfo != NULL)
302      delete oldInfo;
303    
304 <  
305 <  // We need to read in new siminfo object.    
306 <  //parse md file and set up the system
307 <  //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 <  
304 >  SimCreator newCreator;
305 >  SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
306 >    
307 >  // Place molecules
308    Molecule* mol;
309    SimInfo::MoleculeIterator mi;
310    mol = NewInfo->beginMolecule(mi);
311 +  int l = 0;
312  
313 +  for (int i = 0; i < nComponents; i++){
314 +    locator = new MoLocator(NewInfo->getMoleculeStamp(i),
315 +                            NewInfo->getForceField());
316  
317 <  for(int i = -nx; i < nx; i++) {
318 <     for(int j = -ny; j < ny; j++) {
319 <        for(int k = -nz; k < nz; k++) {
320 <          
321 <           //get the position of the cell sites
322 <           simpleLat->getLatticePointsPos(latticePos, i, j, k);
254 <          
255 <           for(int l = 0; l < numMolPerCell; l++) {
256 < #ifdef HAVE_CGAL              
257 <              if (myGeometry->isInsidePolyhedron(latticePos[l][0],latticePos[l][1],latticePos[l][2])){
258 < #endif                              
259 <                 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 <           }
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 <  }
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 >  }
333    
273
274  
334    //fill Hmat
335 <  hmat(0, 0)= nx * latticeConstant;
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) = ny * latticeConstant;
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) = nz * latticeConstant;
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 <  NewInfo->setFinalConfigFileName(outInitFileName);
294 <  writer = new DumpWriter(NewInfo);
352 >  writer = new DumpWriter(NewInfo, outputFileName);
353    
354    if (writer == NULL) {
355 <    std::cerr << "error in creating DumpWriter" << std::endl;
356 <    exit(1);
355 >    sprintf(painCave.errMsg, "Error in creating dumpwrite object ");
356 >    painCave.isFatal = 1;
357 >    simError();
358    }
359    
360    writer->writeDump();
361 <  std::cout << "new initial configuration file: " << outInitFileName
362 <            << " is generated." << std::endl;
363 <  
364 <  //delete objects
365 <  
366 <  //delete oldInfo and oldSimSetup
367 <  
368 <  if (NewInfo != NULL)
369 <    delete NewInfo;
370 <  
312 <  if (writer != NULL)
313 <    delete writer;
314 <  delete simpleLat;    
315 <  cmdline_parser_free(&args_info);
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.....
367 >  sprintf(painCave.errMsg, "A new OOPSE MD file called \"%s\" has been "
368 >          "generated.\n", outputFileName.c_str());
369 >  painCave.isFatal = 0;
370 >  simError();
371    return 0;
372   }
373  
374 < void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
375 <                  int components,int &nummol) {
374 > void createMdFile(const std::string&oldMdFileName,
375 >                  const std::string&newMdFileName,
376 >                  std::vector<int> nMol) {
377    ifstream oldMdFile;
378    ofstream newMdFile;
379    const int MAXLEN = 65535;
# Line 328 | Line 384 | void createMdFile(const std::string&oldMdFileName, con
384    newMdFile.open(newMdFileName.c_str());
385    
386    oldMdFile.getline(buffer, MAXLEN);
387 <  
387 >
388 >  int i = 0;
389    while (!oldMdFile.eof()) {
390      
391      //correct molecule number
392      if (strstr(buffer, "nMol") != NULL) {
393 <      sprintf(buffer, "\tnMol = %i;", numMol);                          
394 <      newMdFile << buffer << std::endl;
393 >      if(i<nMol.size()){
394 >        sprintf(buffer, "\tnMol = %i;", nMol.at(i));                            
395 >        newMdFile << buffer << std::endl;
396 >        i++;
397 >      }
398      } else
399        newMdFile << buffer << std::endl;
400      

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines