ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/applications/nanoparticleBuilder/nanorodBuilder.cpp
Revision: 1864
Committed: Fri Apr 12 21:21:45 2013 UTC (12 years ago) by gezelter
File size: 16334 byte(s)
Log Message:
Added an ellipsoidal shaped lattice

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
86 MoLocator* locator;
87 int nComponents;
88 double latticeConstant;
89 std::vector<double> lc;
90
91 RealType rodRadius;
92 RealType rodLength;
93
94 Mat3x3d hmat;
95 std::vector<Vector3d> latticePos;
96 std::vector<Vector3d> latticeOrt;
97
98 DumpWriter *writer;
99
100 // Parse Command Line Arguments
101 if (cmdline_parser(argc, argv, &args_info) != 0)
102 exit(1);
103
104 /* get lattice type */
105 latticeType = "FCC";
106
107 /* get input file name */
108 if (args_info.inputs_num)
109 inputFileName = args_info.inputs[0];
110 else {
111 sprintf(painCave.errMsg, "No input .md file name was specified "
112 "on the command line");
113 painCave.isFatal = 1;
114 cmdline_parser_print_help();
115 simError();
116 }
117
118 /* parse md file and set up the system */
119 SimCreator oldCreator;
120 SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
121
122 latticeConstant = args_info.latticeConstant_arg;
123 rodRadius = args_info.radius_arg;
124 rodLength = args_info.length_arg;
125 Globals* simParams = oldInfo->getSimParams();
126
127 vector<Vector3d> sites;
128 vector<Vector3d> orientations;
129
130 if (args_info.ellipsoid_flag) {
131 shapedLatticeEllipsoid nanoEllipsoid(latticeConstant, latticeType,
132 rodLength, rodRadius);
133 sites = nanoEllipsoid.getSites();
134 orientations = nanoEllipsoid.getOrientations();
135 } else {
136
137 /* Create nanorod */
138 shapedLatticeRod nanoRod(latticeConstant, latticeType,
139 rodRadius, rodLength);
140 /* Build a lattice and get lattice points for this lattice constant */
141 sites = nanoRod.getSites();
142 orientations = nanoRod.getOrientations();
143 }
144
145 std::vector<int> vacancyTargets;
146 vector<bool> isVacancy;
147
148 Vector3d myLoc;
149 RealType myR;
150
151 for (int i = 0; i < sites.size(); i++)
152 isVacancy.push_back(false);
153
154 // cerr << "checking vacancyPercent" << "\n";
155 if (args_info.vacancyPercent_given) {
156 // cerr << "vacancyPercent given" << "\n";
157 if (args_info.vacancyPercent_arg < 0.0 || args_info.vacancyPercent_arg > 100.0) {
158 sprintf(painCave.errMsg, "vacancyPercent was set to a non-sensical value.");
159 painCave.isFatal = 1;
160 simError();
161 } else {
162 RealType vF = args_info.vacancyPercent_arg / 100.0;
163 // cerr << "vacancyPercent = " << vF << "\n";
164 RealType vIR;
165 RealType vOR;
166 if (args_info.vacancyInnerRadius_given) {
167 vIR = args_info.vacancyInnerRadius_arg;
168 } else {
169 vIR = 0.0;
170 }
171 if (args_info.vacancyOuterRadius_given) {
172 vOR = args_info.vacancyOuterRadius_arg;
173 } else {
174 vOR = rodRadius;
175 }
176 if (vIR >= 0.0 && vOR <= rodRadius && vOR >= vIR) {
177
178 for (int i = 0; i < sites.size(); i++) {
179 myLoc = sites[i];
180 myR = myLoc.length();
181 if (myR >= vIR && myR <= vOR) {
182 vacancyTargets.push_back(i);
183 }
184 }
185 std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
186
187 int nTargets = vacancyTargets.size();
188 vacancyTargets.resize((int)(vF * nTargets));
189
190
191 sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
192 "\tsites between %lf and %lf.", (int) vacancyTargets.size(),
193 vIR, vOR);
194 painCave.isFatal = 0;
195 simError();
196
197 isVacancy.clear();
198 for (int i = 0; i < sites.size(); i++) {
199 bool vac = false;
200 for (int j = 0; j < vacancyTargets.size(); j++) {
201 if (i == vacancyTargets[j]) vac = true;
202 }
203 isVacancy.push_back(vac);
204 }
205
206 } else {
207 sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
208 "\tinner or outer radii. Check their values.");
209 painCave.isFatal = 1;
210 simError();
211 }
212 }
213 }
214
215 /* Get number of lattice sites */
216 int nSites = sites.size() - vacancyTargets.size();
217
218 // cerr << "sites.size() = " << sites.size() << "\n";
219 // cerr << "nSites = " << nSites << "\n";
220 // cerr << "vacancyTargets = " << vacancyTargets.size() << "\n";
221
222 std::vector<Component*> components = simParams->getComponents();
223 std::vector<RealType> molFractions;
224 std::vector<RealType> shellRadii;
225 std::vector<RealType> molecularMasses;
226 std::vector<int> nMol;
227 std::map<int, int> componentFromSite;
228 nComponents = components.size();
229 // cerr << "nComponents = " << nComponents << "\n";
230
231 if (args_info.molFraction_given && args_info.shellRadius_given) {
232 sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
233 "arguments, but not both!");
234 painCave.isFatal = 1;
235 simError();
236 }
237
238 if (nComponents == 1) {
239 molFractions.push_back(1.0);
240 shellRadii.push_back(rodRadius);
241 } else if (args_info.molFraction_given) {
242 if ((int)args_info.molFraction_given == nComponents) {
243 for (int i = 0; i < nComponents; i++) {
244 molFractions.push_back(args_info.molFraction_arg[i]);
245 }
246 } else if ((int)args_info.molFraction_given == nComponents-1) {
247 RealType remainingFraction = 1.0;
248 for (int i = 0; i < nComponents-1; i++) {
249 molFractions.push_back(args_info.molFraction_arg[i]);
250 remainingFraction -= molFractions[i];
251 }
252 molFractions.push_back(remainingFraction);
253 } else {
254 sprintf(painCave.errMsg, "nanorodBuilder can't figure out molFractions "
255 "for all of the components in the <MetaData> block.");
256 painCave.isFatal = 1;
257 simError();
258 }
259 } else if ((int)args_info.shellRadius_given) {
260 if ((int)args_info.shellRadius_given == nComponents) {
261 for (int i = 0; i < nComponents; i++) {
262 shellRadii.push_back(args_info.shellRadius_arg[i]);
263 }
264 } else if ((int)args_info.shellRadius_given == nComponents-1) {
265 for (int i = 0; i < nComponents-1; i++) {
266 shellRadii.push_back(args_info.shellRadius_arg[i]);
267 }
268 shellRadii.push_back(rodRadius);
269 } else {
270 sprintf(painCave.errMsg, "nanorodBuilder can't figure out the\n"
271 "\tshell radii for all of the components in the <MetaData> block.");
272 painCave.isFatal = 1;
273 simError();
274 }
275 } else {
276 sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
277 "\tbut have not specified either molFraction or shellRadius arguments.");
278 painCave.isFatal = 1;
279 simError();
280 }
281
282 if (args_info.molFraction_given) {
283 RealType totalFraction = 0.0;
284
285 /* Do some simple sanity checking*/
286
287 for (int i = 0; i < nComponents; i++) {
288 if (molFractions.at(i) < 0.0) {
289 sprintf(painCave.errMsg, "One of the requested molFractions was"
290 " less than zero!");
291 painCave.isFatal = 1;
292 simError();
293 }
294 if (molFractions.at(i) > 1.0) {
295 sprintf(painCave.errMsg, "One of the requested molFractions was"
296 " greater than one!");
297 painCave.isFatal = 1;
298 simError();
299 }
300 totalFraction += molFractions.at(i);
301 }
302 if (abs(totalFraction - 1.0) > 1e-6) {
303 sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
304 painCave.isFatal = 1;
305 simError();
306 }
307
308 int remaining = nSites;
309 for (int i=0; i < nComponents-1; i++) {
310 nMol.push_back(int((RealType)nSites * molFractions.at(i)));
311 remaining -= nMol.at(i);
312 }
313 nMol.push_back(remaining);
314
315 // recompute actual mol fractions and perform final sanity check:
316
317 int totalMolecules = 0;
318 for (int i=0; i < nComponents; i++) {
319 molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
320 totalMolecules += nMol.at(i);
321 }
322 if (totalMolecules != nSites) {
323 sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
324 "to the number of lattice sites!");
325 painCave.isFatal = 1;
326 simError();
327 }
328 } else {
329
330 for (int i = 0; i < shellRadii.size(); i++) {
331 if (shellRadii.at(i) > rodRadius + 1e-6 ) {
332 sprintf(painCave.errMsg, "One of the shellRadius values exceeds the rod Radius.");
333 painCave.isFatal = 1;
334 simError();
335 }
336 if (shellRadii.at(i) <= 0.0 ) {
337 sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!");
338 painCave.isFatal = 1;
339 simError();
340 }
341 }
342 }
343
344 vector<int> ids;
345 if ((int)args_info.molFraction_given){
346 // cerr << "molFraction given 2" << "\n";
347 sprintf(painCave.errMsg, "Creating a randomized spherically-capped nanorod.");
348 painCave.isFatal = 0;
349 simError();
350 /* Random rod is the default case*/
351
352 for (int 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 simError();
361
362 // RealType smallestSoFar;
363 int myComponent = -1;
364 nMol.clear();
365 nMol.resize(nComponents);
366
367 // cerr << "shellRadii[0] " << shellRadii[0] << "\n";
368 // cerr << "rodRadius " << rodRadius << "\n";
369
370 for (int i = 0; i < sites.size(); i++) {
371 myLoc = sites[i];
372 myR = myLoc.length();
373 // smallestSoFar = rodRadius;
374 // cerr << "vac = " << isVacancy[i]<< "\n";
375
376 if (!isVacancy[i]) {
377
378
379 // for (int j = 0; j < nComponents; j++) {
380 // if (myR <= shellRadii[j]) {
381 // if (shellRadii[j] <= smallestSoFar) {
382 // smallestSoFar = shellRadii[j];
383 // myComponent = j;
384 // }
385 // }
386 // }
387 myComponent = 0;
388 componentFromSite[i] = myComponent;
389 nMol[myComponent]++;
390 // cerr << "nMol for myComp(" << myComponent<<") = " << nMol[myComponent] << "\n";
391 }
392 }
393 }
394 // cerr << "nMol = " << nMol.at(0) << "\n";
395
396 outputFileName = args_info.output_arg;
397
398 //creat new .md file on fly which corrects the number of molecule
399
400 createMdFile(inputFileName, outputFileName, nMol);
401
402 if (oldInfo != NULL)
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 (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 simError();
476 return 0;
477 }
478
479 void createMdFile(const std::string&oldMdFileName,
480 const std::string&newMdFileName,
481 std::vector<int> nMol) {
482 ifstream oldMdFile;
483 ofstream newMdFile;
484 const int MAXLEN = 65535;
485 char buffer[MAXLEN];
486
487 //create new .md file based on old .md file
488 oldMdFile.open(oldMdFileName.c_str());
489 newMdFile.open(newMdFileName.c_str());
490 oldMdFile.getline(buffer, MAXLEN);
491
492 int i = 0;
493 while (!oldMdFile.eof()) {
494
495 //correct molecule number
496 if (strstr(buffer, "nMol") != NULL) {
497 if(i<nMol.size()){
498 sprintf(buffer, "\tnMol = %i;", nMol.at(i));
499 newMdFile << buffer << std::endl;
500 i++;
501 }
502 } else
503 newMdFile << buffer << std::endl;
504
505 oldMdFile.getline(buffer, MAXLEN);
506 }
507
508 oldMdFile.close();
509 newMdFile.close();
510
511 if (i != nMol.size()) {
512 sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
513 "\tstatements in component blocks. Make sure that all\n"
514 "\tcomponents in the template file have nMol=1");
515 painCave.isFatal = 1;
516 simError();
517 }
518
519 }
520