ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/nanoparticleBuilder/nanorodBuilder.cpp
Revision: 2071
Committed: Sat Mar 7 21:41:51 2015 UTC (10 years, 1 month ago) by gezelter
File size: 16458 byte(s)
Log Message:
Reducing the number of warnings when using g++ to compile.

File Contents

# Content
1 /*
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. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
35 *
36 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 * Created by Kelsey M. Stocker on 2/9/12.
42 * @author Kelsey M. Stocker
43 *
44 */
45
46 #include <cstdlib>
47 #include <cstdio>
48 #include <cstring>
49 #include <cmath>
50 #include <iostream>
51 #include <string>
52 #include <map>
53 #include <fstream>
54 #include <algorithm>
55
56 #include "config.h"
57 #include "shapedLatticeRod.hpp"
58 #include "shapedLatticeEllipsoid.hpp"
59 #include "nanorodBuilderCmd.h"
60 #include "lattice/LatticeFactory.hpp"
61 #include "utils/MoLocator.hpp"
62 #include "lattice/Lattice.hpp"
63 #include "brains/Register.hpp"
64 #include "brains/SimInfo.hpp"
65 #include "brains/SimCreator.hpp"
66 #include "io/DumpWriter.hpp"
67 #include "math/Vector3.hpp"
68 #include "math/SquareMatrix3.hpp"
69 #include "utils/StringUtils.hpp"
70
71 using namespace std;
72 using namespace OpenMD;
73 void createMdFile(const std::string&oldMdFileName,
74 const std::string&newMdFileName,
75 std::vector<int> numMol);
76
77 int main(int argc, char *argv []) {
78
79 registerLattice();
80
81 gengetopt_args_info args_info;
82 std::string latticeType;
83 std::string inputFileName;
84 std::string outputFileName;
85 MoLocator* locator;
86 int nComponents;
87 double latticeConstant;
88 RealType rodRadius;
89 RealType rodLength;
90 Mat3x3d hmat;
91 DumpWriter *writer;
92
93 // Parse Command Line Arguments
94 if (cmdline_parser(argc, argv, &args_info) != 0)
95 exit(1);
96
97 /* get lattice type */
98 latticeType = "FCC";
99
100 /* get input file name */
101 if (args_info.inputs_num)
102 inputFileName = args_info.inputs[0];
103 else {
104 sprintf(painCave.errMsg, "No input .md file name was specified "
105 "on the command line");
106 painCave.isFatal = 1;
107 cmdline_parser_print_help();
108 simError();
109 }
110
111 /* parse md file and set up the system */
112 SimCreator oldCreator;
113 SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
114
115 latticeConstant = args_info.latticeConstant_arg;
116 rodRadius = args_info.radius_arg;
117 rodLength = args_info.length_arg;
118 Globals* simParams = oldInfo->getSimParams();
119
120 vector<Vector3d> sites;
121 vector<Vector3d> orientations;
122
123 if (args_info.ellipsoid_flag) {
124 shapedLatticeEllipsoid nanoEllipsoid(latticeConstant, latticeType,
125 rodLength, rodRadius);
126 sites = nanoEllipsoid.getSites();
127 orientations = nanoEllipsoid.getOrientations();
128 } else {
129
130 /* Create nanorod */
131 shapedLatticeRod nanoRod(latticeConstant, latticeType,
132 rodRadius, rodLength);
133 /* Build a lattice and get lattice points for this lattice constant */
134 sites = nanoRod.getSites();
135 orientations = nanoRod.getOrientations();
136 }
137
138 std::vector<std::size_t> vacancyTargets;
139 vector<bool> isVacancy;
140
141 Vector3d myLoc;
142 RealType myR;
143
144 for (std::size_t i = 0; i < sites.size(); i++)
145 isVacancy.push_back(false);
146
147 // cerr << "checking vacancyPercent" << "\n";
148 if (args_info.vacancyPercent_given) {
149 // cerr << "vacancyPercent given" << "\n";
150 if (args_info.vacancyPercent_arg < 0.0 || args_info.vacancyPercent_arg > 100.0) {
151 sprintf(painCave.errMsg,
152 "vacancyPercent was set to a non-sensical value.");
153 painCave.isFatal = 1;
154 simError();
155 } else {
156 RealType vF = args_info.vacancyPercent_arg / 100.0;
157 // cerr << "vacancyPercent = " << vF << "\n";
158 RealType vIR;
159 RealType vOR;
160 if (args_info.vacancyInnerRadius_given) {
161 vIR = args_info.vacancyInnerRadius_arg;
162 } else {
163 vIR = 0.0;
164 }
165 if (args_info.vacancyOuterRadius_given) {
166 vOR = args_info.vacancyOuterRadius_arg;
167 } else {
168 vOR = rodRadius;
169 }
170 if (vIR >= 0.0 && vOR <= rodRadius && vOR >= vIR) {
171
172 for (std::size_t i = 0; i < sites.size(); i++) {
173 myLoc = sites[i];
174 myR = myLoc.length();
175 if (myR >= vIR && myR <= vOR) {
176 vacancyTargets.push_back(i);
177 }
178 }
179 std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
180
181 std::size_t nTargets = vacancyTargets.size();
182 vacancyTargets.resize((int)(vF * nTargets));
183
184
185 sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
186 "\tsites between %lf and %lf.", (int) vacancyTargets.size(),
187 vIR, vOR);
188 painCave.isFatal = 0;
189 painCave.severity = OPENMD_INFO;
190 simError();
191
192 isVacancy.clear();
193 for (std::size_t i = 0; i < sites.size(); i++) {
194 bool vac = false;
195 for (std::size_t j = 0; j < vacancyTargets.size(); j++) {
196 if (i == vacancyTargets[j]) vac = true;
197 }
198 isVacancy.push_back(vac);
199 }
200
201 } else {
202 sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
203 "\tinner or outer radii. Check their values.");
204 painCave.isFatal = 1;
205 simError();
206 }
207 }
208 }
209
210 /* Get number of lattice sites */
211 std::size_t nSites = sites.size() - vacancyTargets.size();
212
213 // cerr << "sites.size() = " << sites.size() << "\n";
214 // cerr << "nSites = " << nSites << "\n";
215 // cerr << "vacancyTargets = " << vacancyTargets.size() << "\n";
216
217 std::vector<Component*> components = simParams->getComponents();
218 std::vector<RealType> molFractions;
219 std::vector<RealType> shellRadii;
220 std::vector<int> nMol;
221 std::map<int, int> componentFromSite;
222 nComponents = components.size();
223 // cerr << "nComponents = " << nComponents << "\n";
224
225 if (args_info.molFraction_given && args_info.shellRadius_given) {
226 sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
227 "arguments, but not both!");
228 painCave.isFatal = 1;
229 simError();
230 }
231
232 if (nComponents == 1) {
233 molFractions.push_back(1.0);
234 shellRadii.push_back(rodRadius);
235 } else if (args_info.molFraction_given) {
236 if ((int)args_info.molFraction_given == nComponents) {
237 for (int i = 0; i < nComponents; i++) {
238 molFractions.push_back(args_info.molFraction_arg[i]);
239 }
240 } else if ((int)args_info.molFraction_given == nComponents-1) {
241 RealType remainingFraction = 1.0;
242 for (int i = 0; i < nComponents-1; i++) {
243 molFractions.push_back(args_info.molFraction_arg[i]);
244 remainingFraction -= molFractions[i];
245 }
246 molFractions.push_back(remainingFraction);
247 } else {
248 sprintf(painCave.errMsg, "nanorodBuilder can't figure out molFractions "
249 "for all of the components in the <MetaData> block.");
250 painCave.isFatal = 1;
251 simError();
252 }
253 } else if ((int)args_info.shellRadius_given) {
254 if ((int)args_info.shellRadius_given == nComponents) {
255 for (int i = 0; i < nComponents; i++) {
256 shellRadii.push_back(args_info.shellRadius_arg[i]);
257 }
258 } else if ((int)args_info.shellRadius_given == nComponents-1) {
259 for (int i = 0; i < nComponents-1; i++) {
260 shellRadii.push_back(args_info.shellRadius_arg[i]);
261 }
262 shellRadii.push_back(rodRadius);
263 } else {
264 sprintf(painCave.errMsg, "nanorodBuilder can't figure out the\n"
265 "\tshell radii for all of the components in the <MetaData> block.");
266 painCave.isFatal = 1;
267 simError();
268 }
269 } else {
270 sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
271 "\tbut have not specified either molFraction or shellRadius arguments.");
272 painCave.isFatal = 1;
273 simError();
274 }
275
276 if (args_info.molFraction_given) {
277 RealType totalFraction = 0.0;
278
279 /* Do some simple sanity checking*/
280
281 for (int i = 0; i < nComponents; i++) {
282 if (molFractions.at(i) < 0.0) {
283 sprintf(painCave.errMsg, "One of the requested molFractions was"
284 " less than zero!");
285 painCave.isFatal = 1;
286 simError();
287 }
288 if (molFractions.at(i) > 1.0) {
289 sprintf(painCave.errMsg, "One of the requested molFractions was"
290 " greater than one!");
291 painCave.isFatal = 1;
292 simError();
293 }
294 totalFraction += molFractions.at(i);
295 }
296 if (abs(totalFraction - 1.0) > 1e-6) {
297 sprintf(painCave.errMsg,
298 "The sum of molFractions was not close enough to 1.0");
299 painCave.isFatal = 1;
300 simError();
301 }
302
303 int remaining = nSites;
304 for (int i=0; i < nComponents-1; i++) {
305 nMol.push_back(int((RealType)nSites * molFractions.at(i)));
306 remaining -= nMol.at(i);
307 }
308 nMol.push_back(remaining);
309
310 // recompute actual mol fractions and perform final sanity check:
311
312 std::size_t totalMolecules = 0;
313 for (int i=0; i < nComponents; i++) {
314 molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
315 totalMolecules += nMol.at(i);
316 }
317 if (totalMolecules != nSites) {
318 sprintf(painCave.errMsg,
319 "Computed total number of molecules is not equal "
320 "to the number of lattice sites!");
321 painCave.isFatal = 1;
322 simError();
323 }
324 } else {
325
326 for (unsigned int i = 0; i < shellRadii.size(); i++) {
327 if (shellRadii.at(i) > rodRadius + 1e-6 ) {
328 sprintf(painCave.errMsg,
329 "One of the shellRadius values exceeds the rod Radius.");
330 painCave.isFatal = 1;
331 simError();
332 }
333 if (shellRadii.at(i) <= 0.0 ) {
334 sprintf(painCave.errMsg,
335 "One of the shellRadius values is smaller than zero!");
336 painCave.isFatal = 1;
337 simError();
338 }
339 }
340 }
341
342 vector<int> ids;
343 if ((int)args_info.molFraction_given){
344 // cerr << "molFraction given 2" << "\n";
345 sprintf(painCave.errMsg,
346 "Creating a randomized spherically-capped nanorod.");
347 painCave.isFatal = 0;
348 painCave.severity = OPENMD_INFO;
349 simError();
350 /* Random rod is the default case*/
351
352 for (std::size_t i = 0; i < sites.size(); i++)
353 if (!isVacancy[i]) ids.push_back(i);
354
355 std::random_shuffle(ids.begin(), ids.end());
356
357 } else{
358 sprintf(painCave.errMsg, "Creating an fcc nanorod.");
359 painCave.isFatal = 0;
360 painCave.severity = OPENMD_INFO;
361 simError();
362
363 // RealType smallestSoFar;
364 int myComponent = -1;
365 nMol.clear();
366 nMol.resize(nComponents);
367
368 // cerr << "shellRadii[0] " << shellRadii[0] << "\n";
369 // cerr << "rodRadius " << rodRadius << "\n";
370
371 for (unsigned int i = 0; i < sites.size(); i++) {
372 myLoc = sites[i];
373 myR = myLoc.length();
374 // smallestSoFar = rodRadius;
375 // cerr << "vac = " << isVacancy[i]<< "\n";
376
377 if (!isVacancy[i]) {
378
379
380 // for (int j = 0; j < nComponents; j++) {
381 // if (myR <= shellRadii[j]) {
382 // if (shellRadii[j] <= smallestSoFar) {
383 // smallestSoFar = shellRadii[j];
384 // myComponent = j;
385 // }
386 // }
387 // }
388 myComponent = 0;
389 componentFromSite[i] = myComponent;
390 nMol[myComponent]++;
391 // cerr << "nMol for myComp(" << myComponent<<") = " << nMol[myComponent] << "\n";
392 }
393 }
394 }
395 // cerr << "nMol = " << nMol.at(0) << "\n";
396
397 outputFileName = args_info.output_arg;
398
399 //creat new .md file on fly which corrects the number of molecule
400
401 createMdFile(inputFileName, outputFileName, nMol);
402
403 delete oldInfo;
404
405 SimCreator newCreator;
406 SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
407
408 // Place molecules
409 Molecule* mol;
410 SimInfo::MoleculeIterator mi;
411 mol = NewInfo->beginMolecule(mi);
412
413 int l = 0;
414
415 for (int i = 0; i < nComponents; i++){
416 locator = new MoLocator(NewInfo->getMoleculeStamp(i),
417 NewInfo->getForceField());
418
419 // cerr << "nMol = " << nMol.at(i) << "\n";
420 if (!args_info.molFraction_given) {
421 for (unsigned int n = 0; n < sites.size(); n++) {
422 if (!isVacancy[n]) {
423 if (componentFromSite[n] == i) {
424 mol = NewInfo->getMoleculeByGlobalIndex(l);
425 locator->placeMol(sites[n], orientations[n], mol);
426 l++;
427 }
428 }
429 }
430 } else {
431 for (int n = 0; n < nMol.at(i); n++) {
432 mol = NewInfo->getMoleculeByGlobalIndex(l);
433 locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
434 l++;
435 }
436 }
437 }
438
439 //fill Hmat
440 hmat(0, 0)= 10.0*rodRadius;
441 hmat(0, 1) = 0.0;
442 hmat(0, 2) = 0.0;
443
444 hmat(1, 0) = 0.0;
445 hmat(1, 1) = 10.0*rodRadius;
446 hmat(1, 2) = 0.0;
447
448 hmat(2, 0) = 0.0;
449 hmat(2, 1) = 0.0;
450 hmat(2, 2) = 5.0*rodLength + 2.0*rodRadius;
451
452 //set Hmat
453 NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
454
455
456 //create dumpwriter and write out the coordinates
457 writer = new DumpWriter(NewInfo, outputFileName);
458
459 if (writer == NULL) {
460 sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
461 painCave.isFatal = 1;
462 simError();
463 }
464
465 writer->writeDump();
466
467 // deleting the writer will put the closing at the end of the dump file
468
469 delete writer;
470
471 // cleanup a by calling sim error.....
472 sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
473 "generated.\n", outputFileName.c_str());
474 painCave.isFatal = 0;
475 painCave.severity = OPENMD_INFO;
476 simError();
477 return 0;
478 }
479
480 void createMdFile(const std::string&oldMdFileName,
481 const std::string&newMdFileName,
482 std::vector<int> nMol) {
483 ifstream oldMdFile;
484 ofstream newMdFile;
485 const int MAXLEN = 65535;
486 char buffer[MAXLEN];
487
488 //create new .md file based on old .md file
489 oldMdFile.open(oldMdFileName.c_str());
490 newMdFile.open(newMdFileName.c_str());
491 oldMdFile.getline(buffer, MAXLEN);
492
493 unsigned int i = 0;
494 while (!oldMdFile.eof()) {
495
496 //correct molecule number
497 if (strstr(buffer, "nMol") != NULL) {
498 if(i<nMol.size()){
499 sprintf(buffer, "\tnMol = %i;", nMol.at(i));
500 newMdFile << buffer << std::endl;
501 i++;
502 }
503 } else
504 newMdFile << buffer << std::endl;
505
506 oldMdFile.getline(buffer, MAXLEN);
507 }
508
509 oldMdFile.close();
510 newMdFile.close();
511
512 if (i != nMol.size()) {
513 sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
514 "\tstatements in component blocks. Make sure that all\n"
515 "\tcomponents in the template file have nMol=1");
516 painCave.isFatal = 1;
517 simError();
518 }
519
520 }
521