ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/randomBuilder/randomBuilder.cpp
(Generate patch)

Comparing trunk/src/applications/randomBuilder/randomBuilder.cpp (file contents):
Revision 1064 by gezelter, Tue Oct 10 02:44:13 2006 UTC vs.
Revision 1065 by gezelter, Tue Oct 10 14:52:20 2006 UTC

# Line 42 | Line 42
42   *
43   *  Created by Charles F. Vardeman II on 10 Apr 2006.
44   *  @author  Charles F. Vardeman II
45 < *  @version $Id: randomBuilder.cpp,v 1.4 2006-10-10 02:44:13 gezelter Exp $
45 > *  @version $Id: randomBuilder.cpp,v 1.5 2006-10-10 14:52:20 gezelter Exp $
46   *
47   */
48  
# Line 73 | Line 73 | void createMdFile(const std::string&oldMdFileName,
73  
74   void createMdFile(const std::string&oldMdFileName,
75                    const std::string&newMdFileName,
76 <                  int components, int* numMol);
76 >                  std::vector<int> nMol);
77  
78   int main(int argc, char *argv []) {
79  
# Line 86 | Line 86 | int main(int argc, char *argv []) {
86    std::string inputFileName;
87    std::string outputFileName;
88    Lattice *simpleLat;
89 –  int* numMol;
89    RealType latticeConstant;
90    std::vector<RealType> lc;
92 –  RealType mass;
91    const RealType rhoConvertConst = 1.661;
92    RealType density;
93    int nx, ny, nz;
# Line 97 | Line 95 | int main(int argc, char *argv []) {
95    MoLocator *locator;
96    std::vector<Vector3d> latticePos;
97    std::vector<Vector3d> latticeOrt;
98 <  int numMolPerCell;
101 <  int curMolIndex;
98 >  int nMolPerCell;
99    DumpWriter *writer;
100  
101    // parse command line arguments
# Line 118 | Line 115 | int main(int argc, char *argv []) {
115      painCave.isFatal = 1;
116      simError();
117    }
118 +  nMolPerCell = simpleLat->getNumSitesPerCell();
119  
120    //get the number of unit cells in each direction:
121  
# Line 148 | Line 146 | int main(int argc, char *argv []) {
146      simError();
147    }
148  
149 +  int nSites = nMolPerCell * nx * ny * nz;
150 +
151    //get input file name
152    if (args_info.inputs_num)
153      inputFileName = args_info.inputs[0];
# Line 164 | Line 164 | int main(int argc, char *argv []) {
164    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
165    Globals* simParams = oldInfo->getSimParams();
166  
167 <  int nComponents =simParams->getNComponents();
168 <  if (oldInfo->getNMoleculeStamp() > 2) {
169 <    sprintf(painCave.errMsg, "randomBuilder can't yet build a system with "
170 <            "more than two components.");
171 <    painCave.isFatal = 1;
172 <    simError();
173 <  }
167 >  // Calculate lattice constant (in Angstroms)
168  
169 <  //get mass of molecule.
169 >  std::vector<Component*> components = simParams->getComponents();
170 >  std::vector<RealType> molFractions;
171 >  std::vector<RealType> molecularMasses;
172 >  std::vector<int> nMol;
173 >  int nComponents = components.size();
174  
175 <  mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
175 >  if (nComponents == 1) {
176 >    molFractions.push_back(1.0);    
177 >  } else {
178 >    if (args_info.molFraction_given == nComponents) {
179 >      for (int i = 0; i < nComponents; i++) {
180 >        molFractions.push_back(args_info.molFraction_arg[i]);
181 >      }
182 >    } else if (args_info.molFraction_given == nComponents-1) {
183 >      RealType remainingFraction = 1.0;
184 >      for (int i = 0; i < nComponents-1; i++) {
185 >        molFractions.push_back(args_info.molFraction_arg[i]);
186 >        remainingFraction -= molFractions[i];
187 >      }
188 >      molFractions.push_back(remainingFraction);
189 >    } else {    
190 >      sprintf(painCave.errMsg, "randomBuilder can't figure out molFractions "
191 >              "for all of the components in the <MetaData> block.");
192 >      painCave.isFatal = 1;
193 >      simError();
194 >    }
195 >  }
196  
197 <  // Create the lattice
197 >  // do some sanity checking:
198 >  
199 >  RealType totalFraction = 0.0;
200  
201 <  simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
202 <
203 <  if (simpleLat == NULL) {
204 <    sprintf(painCave.errMsg, "Error in creating lattice.");
201 >  for (int i = 0; i < nComponents; i++) {
202 >    if (molFractions.at(i) < 0.0) {
203 >      sprintf(painCave.errMsg, "One of the requested molFractions was"
204 >              " less than zero!");
205 >      painCave.isFatal = 1;
206 >      simError();
207 >    }
208 >    if (molFractions.at(i) > 1.0) {
209 >      sprintf(painCave.errMsg, "One of the requested molFractions was"
210 >              " greater than one!");
211 >      painCave.isFatal = 1;
212 >      simError();
213 >    }
214 >    totalFraction += molFractions.at(i);
215 >  }
216 >  if (abs(totalFraction - 1.0) > 1e-6) {
217 >    sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
218      painCave.isFatal = 1;
219      simError();
220    }
221  
222 <  numMolPerCell = simpleLat->getNumSitesPerCell();
223 <
224 <  // Calculate lattice constant (in Angstroms)
222 >  int remaining = nSites;
223 >  for (int i=0; i < nComponents-1; i++) {    
224 >    nMol.push_back(int((RealType)nSites * molFractions.at(i)));
225 >    remaining -= nMol.at(i);
226 >  }
227 >  nMol.push_back(remaining);
228  
229 <  latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
194 <                        (RealType)(1.0 / 3.0));
229 >  // recompute actual mol fractions and perform final sanity check:
230  
231 <  // Set the lattice constant
231 >  int totalMolecules = 0;
232 >  RealType totalMass = 0.0;
233 >  for (int i=0; i < nComponents; i++) {
234 >    molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
235 >    totalMolecules += nMol.at(i);
236 >    molecularMasses.push_back(getMolMass(oldInfo->getMoleculeStamp(i),
237 >                                         oldInfo->getForceField()));
238 >    totalMass += (RealType)(nMol.at(i)) * molecularMasses.at(i);
239 >  }
240 >  RealType avgMass = totalMass / (RealType) totalMolecules;
241  
242 +  if (totalMolecules != nSites) {
243 +    sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
244 +            "to the number of lattice sites!");
245 +    painCave.isFatal = 1;
246 +    simError();
247 +  }
248 +    
249 +  latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
250 +                        (RealType)(1.0 / 3.0));
251 +  
252 +  // Set the lattice constant
253 +  
254    lc.push_back(latticeConstant);
255    simpleLat->setLatticeConstant(lc);
256 <
201 <  // Calculate the total number of molecules
202 <
203 <  int totMol = nx * ny * nz * numMolPerCell;
204 <
256 >  
257    // Calculate the lattice sites and fill the lattice vector.
258  
259    // Get the standard orientations of the cell sites
# Line 219 | Line 271 | int main(int argc, char *argv []) {
271  
272          simpleLat->getLatticePointsPos(latticePos, i, j, k);
273  
274 <        for(int l = 0; l < numMolPerCell; l++) {
274 >        for(int l = 0; l < nMolPerCell; l++) {
275            sites.push_back(latticePos[l]);
276            orientations.push_back(latticeOrt[l]);
277          }
278        }
279      }
280    }
229 –
230 –  int numSites = sites.size();
231 –
232 –  numMol = new int[nComponents];
233 –  if (nComponents != args_info.molFraction_given && nComponents != 1){
234 –    sprintf(painCave.errMsg, "There needs to be the same number of "
235 –            "molFraction arguments as there are components in the "
236 –            "<MetaData> block.");
237 –    painCave.isFatal = 1;
238 –    simError();
239 –  }
240 –  int totComponents = 0;
241 –  for (int i = 0;i<nComponents-1;i++){
242 –    numMol[i] = int((RealType)numSites * args_info.molFraction_arg[i]);
243 –    std::cout<<numMol[i]<<std::endl;
244 –    totComponents += numMol[i];
245 –  }
246 –  numMol[nComponents-1] = numSites - totComponents;
281    
282    outputFileName = args_info.output_arg;
283  
284    // create a new .md file on the fly which corrects the number of molecules
285  
286 <  createMdFile(inputFileName, outputFileName, nComponents, numMol);
286 >  createMdFile(inputFileName, outputFileName, nMol);
287  
288    if (oldInfo != NULL)
289      delete oldInfo;
# Line 280 | Line 314 | int main(int argc, char *argv []) {
314  
315    // place the molecules
316  
283 –  curMolIndex = 0;
284 –
317    // Randomize a vector of ints:
318  
319    vector<int> ids;
# Line 293 | Line 325 | int main(int argc, char *argv []) {
325    for (int i = 0; i < nComponents; i++){
326      locator = new MoLocator(newInfo->getMoleculeStamp(i),
327                              newInfo->getForceField());
328 <    for (int n = 0; n < numMol[i]; n++) {
328 >    for (int n = 0; n < nMol.at(i); n++) {
329        mol = newInfo->getMoleculeByGlobalIndex(l);
330        locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
331        l++;
# Line 317 | Line 349 | int main(int argc, char *argv []) {
349    delete writer;
350  
351    sprintf(painCave.errMsg, "A new OOPSE MD file called \"%s\" has been "
352 <          "generated.", outputFileName.c_str());
352 >          "generated.\n", outputFileName.c_str());
353    painCave.isFatal = 0;
354    simError();
355    return 0;
# Line 325 | Line 357 | void createMdFile(const std::string&oldMdFileName,
357  
358   void createMdFile(const std::string&oldMdFileName,
359                    const std::string&newMdFileName,
360 <                  int components, int* numMol) {
360 >                  std::vector<int> nMol) {
361    ifstream oldMdFile;
362    ofstream newMdFile;
363    const int MAXLEN = 65535;
# Line 343 | Line 375 | void createMdFile(const std::string&oldMdFileName,
375      
376      //correct molecule number
377      if (strstr(buffer, "nMol") != NULL) {
378 <      if(i<components){
379 <        sprintf(buffer, "\tnMol = %i;", numMol[i]);
378 >      if(i<nMol.size()){
379 >        sprintf(buffer, "\tnMol = %i;", nMol.at(i));
380          newMdFile << buffer << std::endl;
381          i++;
382        }

Diff Legend

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