ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/nanoparticleBuilder/nanorodBuilder.cpp
Revision: 1977
Committed: Wed Mar 12 21:35:23 2014 UTC (11 years, 1 month ago) by gezelter
File size: 16416 byte(s)
Log Message:
Removed some scary errors that are really supposed to be OPENMD_INFO messages

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<int> vacancyTargets;
139 vector<bool> isVacancy;
140
141 Vector3d myLoc;
142 RealType myR;
143
144 for (unsigned int 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, "vacancyPercent was set to a non-sensical value.");
152 painCave.isFatal = 1;
153 simError();
154 } else {
155 RealType vF = args_info.vacancyPercent_arg / 100.0;
156 // cerr << "vacancyPercent = " << vF << "\n";
157 RealType vIR;
158 RealType vOR;
159 if (args_info.vacancyInnerRadius_given) {
160 vIR = args_info.vacancyInnerRadius_arg;
161 } else {
162 vIR = 0.0;
163 }
164 if (args_info.vacancyOuterRadius_given) {
165 vOR = args_info.vacancyOuterRadius_arg;
166 } else {
167 vOR = rodRadius;
168 }
169 if (vIR >= 0.0 && vOR <= rodRadius && vOR >= vIR) {
170
171 for (unsigned int i = 0; i < sites.size(); i++) {
172 myLoc = sites[i];
173 myR = myLoc.length();
174 if (myR >= vIR && myR <= vOR) {
175 vacancyTargets.push_back(i);
176 }
177 }
178 std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
179
180 int nTargets = vacancyTargets.size();
181 vacancyTargets.resize((int)(vF * nTargets));
182
183
184 sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
185 "\tsites between %lf and %lf.", (int) vacancyTargets.size(),
186 vIR, vOR);
187 painCave.isFatal = 0;
188 painCave.severity = OPENMD_INFO;
189 simError();
190
191 isVacancy.clear();
192 for (unsigned int i = 0; i < sites.size(); i++) {
193 bool vac = false;
194 for (unsigned int j = 0; j < vacancyTargets.size(); j++) {
195 if (i == vacancyTargets[j]) vac = true;
196 }
197 isVacancy.push_back(vac);
198 }
199
200 } else {
201 sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
202 "\tinner or outer radii. Check their values.");
203 painCave.isFatal = 1;
204 simError();
205 }
206 }
207 }
208
209 /* Get number of lattice sites */
210 int nSites = sites.size() - vacancyTargets.size();
211
212 // cerr << "sites.size() = " << sites.size() << "\n";
213 // cerr << "nSites = " << nSites << "\n";
214 // cerr << "vacancyTargets = " << vacancyTargets.size() << "\n";
215
216 std::vector<Component*> components = simParams->getComponents();
217 std::vector<RealType> molFractions;
218 std::vector<RealType> shellRadii;
219 std::vector<int> nMol;
220 std::map<int, int> componentFromSite;
221 nComponents = components.size();
222 // cerr << "nComponents = " << nComponents << "\n";
223
224 if (args_info.molFraction_given && args_info.shellRadius_given) {
225 sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
226 "arguments, but not both!");
227 painCave.isFatal = 1;
228 simError();
229 }
230
231 if (nComponents == 1) {
232 molFractions.push_back(1.0);
233 shellRadii.push_back(rodRadius);
234 } else if (args_info.molFraction_given) {
235 if ((int)args_info.molFraction_given == nComponents) {
236 for (int i = 0; i < nComponents; i++) {
237 molFractions.push_back(args_info.molFraction_arg[i]);
238 }
239 } else if ((int)args_info.molFraction_given == nComponents-1) {
240 RealType remainingFraction = 1.0;
241 for (int i = 0; i < nComponents-1; i++) {
242 molFractions.push_back(args_info.molFraction_arg[i]);
243 remainingFraction -= molFractions[i];
244 }
245 molFractions.push_back(remainingFraction);
246 } else {
247 sprintf(painCave.errMsg, "nanorodBuilder can't figure out molFractions "
248 "for all of the components in the <MetaData> block.");
249 painCave.isFatal = 1;
250 simError();
251 }
252 } else if ((int)args_info.shellRadius_given) {
253 if ((int)args_info.shellRadius_given == nComponents) {
254 for (int i = 0; i < nComponents; i++) {
255 shellRadii.push_back(args_info.shellRadius_arg[i]);
256 }
257 } else if ((int)args_info.shellRadius_given == nComponents-1) {
258 for (int i = 0; i < nComponents-1; i++) {
259 shellRadii.push_back(args_info.shellRadius_arg[i]);
260 }
261 shellRadii.push_back(rodRadius);
262 } else {
263 sprintf(painCave.errMsg, "nanorodBuilder can't figure out the\n"
264 "\tshell radii for all of the components in the <MetaData> block.");
265 painCave.isFatal = 1;
266 simError();
267 }
268 } else {
269 sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
270 "\tbut have not specified either molFraction or shellRadius arguments.");
271 painCave.isFatal = 1;
272 simError();
273 }
274
275 if (args_info.molFraction_given) {
276 RealType totalFraction = 0.0;
277
278 /* Do some simple sanity checking*/
279
280 for (int i = 0; i < nComponents; i++) {
281 if (molFractions.at(i) < 0.0) {
282 sprintf(painCave.errMsg, "One of the requested molFractions was"
283 " less than zero!");
284 painCave.isFatal = 1;
285 simError();
286 }
287 if (molFractions.at(i) > 1.0) {
288 sprintf(painCave.errMsg, "One of the requested molFractions was"
289 " greater than one!");
290 painCave.isFatal = 1;
291 simError();
292 }
293 totalFraction += molFractions.at(i);
294 }
295 if (abs(totalFraction - 1.0) > 1e-6) {
296 sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
297 painCave.isFatal = 1;
298 simError();
299 }
300
301 int remaining = nSites;
302 for (int i=0; i < nComponents-1; i++) {
303 nMol.push_back(int((RealType)nSites * molFractions.at(i)));
304 remaining -= nMol.at(i);
305 }
306 nMol.push_back(remaining);
307
308 // recompute actual mol fractions and perform final sanity check:
309
310 int totalMolecules = 0;
311 for (int i=0; i < nComponents; i++) {
312 molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
313 totalMolecules += nMol.at(i);
314 }
315 if (totalMolecules != nSites) {
316 sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
317 "to the number of lattice sites!");
318 painCave.isFatal = 1;
319 simError();
320 }
321 } else {
322
323 for (unsigned int i = 0; i < shellRadii.size(); i++) {
324 if (shellRadii.at(i) > rodRadius + 1e-6 ) {
325 sprintf(painCave.errMsg, "One of the shellRadius values exceeds the rod Radius.");
326 painCave.isFatal = 1;
327 simError();
328 }
329 if (shellRadii.at(i) <= 0.0 ) {
330 sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!");
331 painCave.isFatal = 1;
332 simError();
333 }
334 }
335 }
336
337 vector<int> ids;
338 if ((int)args_info.molFraction_given){
339 // cerr << "molFraction given 2" << "\n";
340 sprintf(painCave.errMsg, "Creating a randomized spherically-capped nanorod.");
341 painCave.isFatal = 0;
342 painCave.severity = OPENMD_INFO;
343 simError();
344 /* Random rod is the default case*/
345
346 for (unsigned int i = 0; i < sites.size(); i++)
347 if (!isVacancy[i]) ids.push_back(i);
348
349 std::random_shuffle(ids.begin(), ids.end());
350
351 } else{
352 sprintf(painCave.errMsg, "Creating an fcc nanorod.");
353 painCave.isFatal = 0;
354 painCave.severity = OPENMD_INFO;
355 simError();
356
357 // RealType smallestSoFar;
358 int myComponent = -1;
359 nMol.clear();
360 nMol.resize(nComponents);
361
362 // cerr << "shellRadii[0] " << shellRadii[0] << "\n";
363 // cerr << "rodRadius " << rodRadius << "\n";
364
365 for (unsigned int i = 0; i < sites.size(); i++) {
366 myLoc = sites[i];
367 myR = myLoc.length();
368 // smallestSoFar = rodRadius;
369 // cerr << "vac = " << isVacancy[i]<< "\n";
370
371 if (!isVacancy[i]) {
372
373
374 // for (int j = 0; j < nComponents; j++) {
375 // if (myR <= shellRadii[j]) {
376 // if (shellRadii[j] <= smallestSoFar) {
377 // smallestSoFar = shellRadii[j];
378 // myComponent = j;
379 // }
380 // }
381 // }
382 myComponent = 0;
383 componentFromSite[i] = myComponent;
384 nMol[myComponent]++;
385 // cerr << "nMol for myComp(" << myComponent<<") = " << nMol[myComponent] << "\n";
386 }
387 }
388 }
389 // cerr << "nMol = " << nMol.at(0) << "\n";
390
391 outputFileName = args_info.output_arg;
392
393 //creat new .md file on fly which corrects the number of molecule
394
395 createMdFile(inputFileName, outputFileName, nMol);
396
397 delete oldInfo;
398
399 SimCreator newCreator;
400 SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
401
402 // Place molecules
403 Molecule* mol;
404 SimInfo::MoleculeIterator mi;
405 mol = NewInfo->beginMolecule(mi);
406
407 int l = 0;
408
409 for (unsigned int i = 0; i < nComponents; i++){
410 locator = new MoLocator(NewInfo->getMoleculeStamp(i),
411 NewInfo->getForceField());
412
413 // cerr << "nMol = " << nMol.at(i) << "\n";
414 if (!args_info.molFraction_given) {
415 for (unsigned int n = 0; n < sites.size(); n++) {
416 if (!isVacancy[n]) {
417 if (componentFromSite[n] == i) {
418 mol = NewInfo->getMoleculeByGlobalIndex(l);
419 locator->placeMol(sites[n], orientations[n], mol);
420 l++;
421 }
422 }
423 }
424 } else {
425 for (unsigned int n = 0; n < nMol.at(i); n++) {
426 mol = NewInfo->getMoleculeByGlobalIndex(l);
427 locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
428 l++;
429 }
430 }
431 }
432
433 //fill Hmat
434 hmat(0, 0)= 10.0*rodRadius;
435 hmat(0, 1) = 0.0;
436 hmat(0, 2) = 0.0;
437
438 hmat(1, 0) = 0.0;
439 hmat(1, 1) = 10.0*rodRadius;
440 hmat(1, 2) = 0.0;
441
442 hmat(2, 0) = 0.0;
443 hmat(2, 1) = 0.0;
444 hmat(2, 2) = 5.0*rodLength + 2.0*rodRadius;
445
446 //set Hmat
447 NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
448
449
450 //create dumpwriter and write out the coordinates
451 writer = new DumpWriter(NewInfo, outputFileName);
452
453 if (writer == NULL) {
454 sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
455 painCave.isFatal = 1;
456 simError();
457 }
458
459 writer->writeDump();
460
461 // deleting the writer will put the closing at the end of the dump file
462
463 delete writer;
464
465 // cleanup a by calling sim error.....
466 sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
467 "generated.\n", outputFileName.c_str());
468 painCave.isFatal = 0;
469 painCave.severity = OPENMD_INFO;
470 simError();
471 return 0;
472 }
473
474 void createMdFile(const std::string&oldMdFileName,
475 const std::string&newMdFileName,
476 std::vector<int> nMol) {
477 ifstream oldMdFile;
478 ofstream newMdFile;
479 const int MAXLEN = 65535;
480 char buffer[MAXLEN];
481
482 //create new .md file based on old .md file
483 oldMdFile.open(oldMdFileName.c_str());
484 newMdFile.open(newMdFileName.c_str());
485 oldMdFile.getline(buffer, MAXLEN);
486
487 unsigned int i = 0;
488 while (!oldMdFile.eof()) {
489
490 //correct molecule number
491 if (strstr(buffer, "nMol") != NULL) {
492 if(i<nMol.size()){
493 sprintf(buffer, "\tnMol = %i;", nMol.at(i));
494 newMdFile << buffer << std::endl;
495 i++;
496 }
497 } else
498 newMdFile << buffer << std::endl;
499
500 oldMdFile.getline(buffer, MAXLEN);
501 }
502
503 oldMdFile.close();
504 newMdFile.close();
505
506 if (i != nMol.size()) {
507 sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
508 "\tstatements in component blocks. Make sure that all\n"
509 "\tcomponents in the template file have nMol=1");
510 painCave.isFatal = 1;
511 simError();
512 }
513
514 }
515