ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp
Revision: 3046
Committed: Sat Oct 14 20:21:26 2006 UTC (17 years, 11 months ago) by gezelter
File size: 10588 byte(s)
Log Message:
fixed a few bugs

File Contents

# User Rev Content
1 chuckv 2352 /*
2     * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
6     * redistribute this software in source and binary code form, provided
7     * that the following conditions are met:
8     *
9     * 1. Acknowledgement of the program authors must be made in any
10     * publication of scientific results based in part on use of the
11     * program. An acceptable form of acknowledgement is citation of
12     * the article in which the program was described (Matthew
13     * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     * Parallel Simulation Engine for Molecular Dynamics,"
16     * J. Comput. Chem. 26, pp. 252-271 (2005))
17     *
18     * 2. Redistributions of source code must retain the above copyright
19     * notice, this list of conditions and the following disclaimer.
20     *
21     * 3. Redistributions in binary form must reproduce the above copyright
22     * notice, this list of conditions and the following disclaimer in the
23     * documentation and/or other materials provided with the
24     * distribution.
25     *
26     * This software is provided "AS IS," without a warranty of any
27     * kind. All express or implied conditions, representations and
28     * warranties, including any implied warranty of merchantability,
29     * fitness for a particular purpose or non-infringement, are hereby
30     * excluded. The University of Notre Dame and its licensors shall not
31     * be liable for any damages suffered by licensee as a result of
32     * using, modifying or distributing the software or its
33     * derivatives. In no event will the University of Notre Dame or its
34     * licensors be liable for any lost revenue, profit or data, or for
35     * direct, indirect, special, consequential, incidental or punitive
36     * damages, however caused and regardless of the theory of liability,
37     * arising out of the use of or inability to use software, even if the
38     * University of Notre Dame has been advised of the possibility of
39     * such damages.
40     */
41    
42     #include <cstdlib>
43     #include <cstdio>
44     #include <cstring>
45     #include <cmath>
46     #include <iostream>
47     #include <string>
48     #include <map>
49     #include <fstream>
50 chuckv 2657 #include <algorithm>
51 chuckv 2352
52     #include "config.h"
53 chuckv 2657 #include "shapedLatticeSpherical.hpp"
54 chuckv 2352 #include "nanoparticleBuilderCmd.h"
55     #include "lattice/LatticeFactory.hpp"
56     #include "utils/MoLocator.hpp"
57     #include "lattice/Lattice.hpp"
58     #include "brains/Register.hpp"
59     #include "brains/SimInfo.hpp"
60     #include "brains/SimCreator.hpp"
61     #include "io/DumpWriter.hpp"
62     #include "math/Vector3.hpp"
63     #include "math/SquareMatrix3.hpp"
64     #include "utils/StringUtils.hpp"
65    
66     using namespace std;
67     using namespace oopse;
68     void createMdFile(const std::string&oldMdFileName,
69     const std::string&newMdFileName,
70 chuckv 3044 std::vector<int> numMol);
71 chuckv 2352
72     int main(int argc, char *argv []) {
73    
74     //register force fields
75     registerForceFields();
76     registerLattice();
77    
78     gengetopt_args_info args_info;
79     std::string latticeType;
80     std::string inputFileName;
81 chuckv 3044 std::string outputFileName;
82 chuckv 2352
83     Lattice *simpleLat;
84 chuckv 2657 MoLocator* locator;
85     int* numMol;
86     int nComponents;
87 chuckv 2352 double latticeConstant;
88     std::vector<double> lc;
89 chuckv 3044 double mass;
90 chuckv 2352 const double rhoConvertConst = 1.661;
91     double density;
92 chuckv 2657 double particleRadius;
93 chuckv 2352
94     Mat3x3d hmat;
95     std::vector<Vector3d> latticePos;
96     std::vector<Vector3d> latticeOrt;
97     int numMolPerCell;
98     int nShells; /* Number of shells in nanoparticle*/
99 chuckv 3044
100 chuckv 2352 DumpWriter *writer;
101    
102 chuckv 2574 // Parse Command Line Arguments
103 chuckv 2352 if (cmdline_parser(argc, argv, &args_info) != 0)
104     exit(1);
105 gezelter 3046
106 chuckv 2352 /* get lattice type */
107 chuckv 3044 latticeType = "FCC";
108    
109 chuckv 2352 /* get input file name */
110     if (args_info.inputs_num)
111     inputFileName = args_info.inputs[0];
112     else {
113 chuckv 3044 sprintf(painCave.errMsg, "No input .md file name was specified"
114     "on the command line");
115     painCave.isFatal = 1;
116 chuckv 2352 cmdline_parser_print_help();
117 chuckv 3044 simError();
118 chuckv 2352 }
119    
120     /* parse md file and set up the system */
121     SimCreator oldCreator;
122     SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
123    
124     latticeConstant = args_info.latticeCnst_arg;
125     particleRadius = args_info.radius_arg;
126 chuckv 2657 Globals* simParams = oldInfo->getSimParams();
127 chuckv 2352
128 chuckv 2657 /* Create nanoparticle */
129 gezelter 3046 shapedLatticeSpherical nanoParticle(latticeConstant, latticeType,
130     particleRadius);
131 chuckv 2672
132 chuckv 2657 /* Build a lattice and get lattice points for this lattice constant */
133 gezelter 3046 vector<Vector3d> sites = nanoParticle.getSites();
134     vector<Vector3d> orientations = nanoParticle.getOrientations();
135 chuckv 3044
136     std::cout <<"nSites: " << sites.size() << std::endl;
137    
138 chuckv 2657 /* Get number of lattice sites */
139 chuckv 3044 int nSites = sites.size();
140    
141     std::vector<Component*> components = simParams->getComponents();
142     std::vector<RealType> molFractions;
143     std::vector<RealType> molecularMasses;
144     std::vector<int> nMol;
145     nComponents = components.size();
146    
147     if (nComponents == 1) {
148     molFractions.push_back(1.0);
149     } else {
150 gezelter 3046 if (args_info.molFraction_given == nComponents) {
151     for (int i = 0; i < nComponents; i++) {
152     molFractions.push_back(args_info.molFraction_arg[i]);
153     }
154 chuckv 3044 } else if (args_info.molFraction_given == nComponents-1) {
155 gezelter 3046 RealType remainingFraction = 1.0;
156     for (int i = 0; i < nComponents-1; i++) {
157     molFractions.push_back(args_info.molFraction_arg[i]);
158     remainingFraction -= molFractions[i];
159     }
160     molFractions.push_back(remainingFraction);
161     } else {
162     sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out molFractions "
163     "for all of the components in the <MetaData> block.");
164     painCave.isFatal = 1;
165     simError();
166     }
167     }
168    
169 chuckv 3044 RealType totalFraction = 0.0;
170 gezelter 3046
171 chuckv 3044 /* 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!");
177     painCave.isFatal = 1;
178     simError();
179     }
180     if (molFractions.at(i) > 1.0) {
181     sprintf(painCave.errMsg, "One of the requested molFractions was"
182     " greater than one!");
183     painCave.isFatal = 1;
184     simError();
185     }
186     totalFraction += molFractions.at(i);
187     }
188     if (abs(totalFraction - 1.0) > 1e-6) {
189     sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
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 chuckv 2352 /* Random particle is the default case*/
226     if (!args_info.ShellRadius_given){
227 chuckv 2737 /* do the iPod thing, Shuffle da vector */
228 chuckv 3044 std::random_shuffle(ids.begin(), ids.end());
229 chuckv 2352 } 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 chuckv 3044 sprintf(painCave.errMsg, "Number of .md components "
233     "does not match the number of shell radius specifications");
234     painCave.isFatal = 1;
235     simError();
236 chuckv 2657 }
237 chuckv 2352
238     }
239    
240 chuckv 2737
241 chuckv 2352
242 chuckv 2737
243 chuckv 3044 outputFileName = args_info.output_arg;
244    
245 chuckv 2737
246 chuckv 3044 //creat new .md file on fly which corrects the number of molecule
247     createMdFile(inputFileName, outputFileName, nMol);
248 chuckv 2352
249     if (oldInfo != NULL)
250     delete oldInfo;
251    
252    
253     // We need to read in new siminfo object.
254     //parse md file and set up the system
255     //SimCreator NewCreator;
256 chuckv 2737 SimCreator newCreator;
257 chuckv 3044 SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
258 chuckv 2352
259    
260 chuckv 2657 // Place molecules
261 chuckv 2352 Molecule* mol;
262     SimInfo::MoleculeIterator mi;
263     mol = NewInfo->beginMolecule(mi);
264 chuckv 2737 int l = 0;
265 chuckv 3044
266     for (int i = 0; i < nComponents; i++){
267     locator = new MoLocator(NewInfo->getMoleculeStamp(i),
268     NewInfo->getForceField());
269     for (int n = 0; n < nMol.at(i); n++) {
270     mol = NewInfo->getMoleculeByGlobalIndex(l);
271     locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
272     l++;
273     }
274 chuckv 2657 }
275 chuckv 2352
276 chuckv 3044
277    
278 chuckv 2352 //fill Hmat
279 chuckv 3044 hmat(0, 0)= 2.0*particleRadius;
280 chuckv 2352 hmat(0, 1) = 0.0;
281     hmat(0, 2) = 0.0;
282    
283     hmat(1, 0) = 0.0;
284 chuckv 3044 hmat(1, 1) = 2.0*particleRadius;
285 chuckv 2352 hmat(1, 2) = 0.0;
286    
287     hmat(2, 0) = 0.0;
288     hmat(2, 1) = 0.0;
289 chuckv 3044 hmat(2, 2) = 2.0*particleRadius;
290 chuckv 2352
291     //set Hmat
292     NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
293    
294    
295     //create dumpwriter and write out the coordinates
296 gezelter 3046 writer = new DumpWriter(NewInfo, outputFileName);
297 chuckv 2352
298     if (writer == NULL) {
299 chuckv 3044 sprintf(painCave.errMsg, "Error in creating dumpwrite object ");
300     painCave.isFatal = 1;
301     simError();
302 chuckv 2352 }
303    
304     writer->writeDump();
305 chuckv 3044
306     // deleting the writer will put the closing at the end of the dump file
307 gezelter 3046
308 chuckv 3044 delete writer;
309    
310     // cleanup a by calling sim error.....
311     sprintf(painCave.errMsg, "A new OOPSE MD file called \"%s\" has been "
312     "generated.\n", outputFileName.c_str());
313     painCave.isFatal = 0;
314     simError();
315 chuckv 2352 return 0;
316     }
317    
318     void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
319 chuckv 3044 std::vector<int> nMol) {
320 chuckv 2352 ifstream oldMdFile;
321     ofstream newMdFile;
322     const int MAXLEN = 65535;
323     char buffer[MAXLEN];
324    
325     //create new .md file based on old .md file
326     oldMdFile.open(oldMdFileName.c_str());
327     newMdFile.open(newMdFileName.c_str());
328    
329     oldMdFile.getline(buffer, MAXLEN);
330 chuckv 2657
331     int i = 0;
332 chuckv 2352 while (!oldMdFile.eof()) {
333    
334     //correct molecule number
335     if (strstr(buffer, "nMol") != NULL) {
336 chuckv 3044 if(i<nMol.size()){
337     sprintf(buffer, "\tnMol = %i;", nMol.at(i));
338 chuckv 2737 newMdFile << buffer << std::endl;
339     i++;
340     }
341 chuckv 2352 } else
342     newMdFile << buffer << std::endl;
343    
344     oldMdFile.getline(buffer, MAXLEN);
345     }
346    
347     oldMdFile.close();
348     newMdFile.close();
349     }
350