58#include "brains/Register.hpp"
61#include "io/DumpWriter.hpp"
62#include "lattice/Lattice.hpp"
66#include "shapedLatticePentRod.hpp"
67#include "shapedLatticeRod.hpp"
68#include "utils/Constants.hpp"
69#include "utils/MoLocator.hpp"
74void createMdFile(
const std::string& oldMdFileName,
75 const std::string& newMdFileName, std::vector<int> numMol);
77int main(
int argc,
char* argv[]) {
81 std::string latticeType;
82 std::string inputFileName;
83 std::string outputFileName;
86 double latticeConstant;
93 if (cmdline_parser(argc, argv, &args_info) != 0) exit(1);
99 if (args_info.inputs_num)
100 inputFileName = args_info.inputs[0];
102 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
103 "No input .omd file name was specified "
104 "on the command line");
105 painCave.isFatal = 1;
106 cmdline_parser_print_help();
112 SimInfo* oldInfo = oldCreator.createSim(inputFileName,
false);
114 latticeConstant = args_info.latticeConstant_arg;
115 rodRadius = args_info.radius_arg;
116 rodLength = args_info.length_arg;
117 Globals* simParams = oldInfo->getSimParams();
124 std::random_device rd;
126 std::mt19937 gen(rd());
131 RealType phi, theta, psi;
145 RotMat3x3d rotation45(45.0 * Constants::PI / 180.0, 0.0, 0.0);
158 theta = 72.0 * Constants::PI / 180.0;
174 vector<Vector3d> getsites = nanoRod.getSites();
175 vector<Vector3d> getorientations = nanoRod.getOrientations();
176 vector<Vector3d> sites;
177 vector<Vector3d> orientations;
179 for (
unsigned int index = 0; index < getsites.size(); index++) {
181 Vector3d myOrient = getorientations[index];
182 Vector3d mySite2 = rotation45 * mySite;
183 Vector3d o2 = rotation45 * myOrient;
184 sites.push_back(mySite2);
185 orientations.push_back(o2);
187 mySite2 = rotation72 * mySite2;
188 o2 = rotation72 * o2;
189 sites.push_back(mySite2);
190 orientations.push_back(o2);
192 mySite2 = rotation72 * mySite2;
193 o2 = rotation72 * o2;
194 sites.push_back(mySite2);
195 orientations.push_back(o2);
197 mySite2 = rotation72 * mySite2;
198 o2 = rotation72 * o2;
199 sites.push_back(mySite2);
200 orientations.push_back(o2);
202 mySite2 = rotation72 * mySite2;
203 o2 = rotation72 * o2;
204 sites.push_back(mySite2);
205 orientations.push_back(o2);
208 int nCenter = int((rodLength + 1.154700538 * rodRadius) / 2.88);
210 for (
unsigned int index = 0; index <= 0.5 * nCenter; index++) {
211 Vector3d myLoc_top(2.88 * index, 0.0, 0.0);
212 sites.push_back(myLoc_top);
213 orientations.push_back(
Vector3d(0.0));
216 for (
unsigned int index = 1; index <= 0.5 * nCenter; index++) {
217 Vector3d myLoc_bottom(-2.88 * index, 0.0, 0.0);
218 sites.push_back(myLoc_bottom);
219 orientations.push_back(
Vector3d(0.0));
222 std::vector<std::size_t> vacancyTargets;
223 vector<bool> isVacancy;
228 for (
unsigned int i = 0; i < sites.size(); i++)
229 isVacancy.push_back(
false);
232 if (args_info.vacancyPercent_given) {
234 if (args_info.vacancyPercent_arg < 0.0 ||
235 args_info.vacancyPercent_arg > 100.0) {
236 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
237 "vacancyPercent was set to a non-sensical value.");
238 painCave.isFatal = 1;
241 RealType vF = args_info.vacancyPercent_arg / 100.0;
245 if (args_info.vacancyInnerRadius_given) {
246 vIR = args_info.vacancyInnerRadius_arg;
250 if (args_info.vacancyOuterRadius_given) {
251 vOR = args_info.vacancyOuterRadius_arg;
255 if (vIR >= 0.0 && vOR <= rodRadius && vOR >= vIR) {
256 for (
unsigned int i = 0; i < sites.size(); i++) {
258 myR = myLoc.length();
259 if (myR >= vIR && myR <= vOR) { vacancyTargets.push_back(i); }
261 std::shuffle(vacancyTargets.begin(), vacancyTargets.end(), gen);
263 int nTargets = vacancyTargets.size();
264 vacancyTargets.resize((
int)(vF * nTargets));
266 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
267 "Removing %d atoms from randomly-selected\n"
268 "\tsites between %lf and %lf.",
269 (
int)vacancyTargets.size(), vIR, vOR);
270 painCave.severity = OPENMD_INFO;
271 painCave.isFatal = 0;
275 for (std::size_t i = 0; i < sites.size(); i++) {
277 for (std::size_t j = 0; j < vacancyTargets.size(); j++) {
278 if (i == vacancyTargets[j]) vac =
true;
280 isVacancy.push_back(vac);
284 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
285 "Something is strange about the vacancy\n"
286 "\tinner or outer radii. Check their values.");
287 painCave.isFatal = 1;
294 int nSites = sites.size() - vacancyTargets.size();
300 std::vector<Component*> components = simParams->getComponents();
301 std::vector<RealType> molFractions;
302 std::vector<RealType> shellRadii;
303 std::vector<int> nMol;
304 std::map<int, int> componentFromSite;
305 nComponents = components.size();
308 if (args_info.molFraction_given && args_info.shellRadius_given) {
309 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
310 "Specify either molFraction or shellRadius "
311 "arguments, but not both!");
312 painCave.isFatal = 1;
316 if (nComponents == 1) {
317 molFractions.push_back(1.0);
318 shellRadii.push_back(rodRadius);
319 }
else if (args_info.molFraction_given) {
320 if ((
int)args_info.molFraction_given == nComponents) {
321 for (
int i = 0; i < nComponents; i++) {
322 molFractions.push_back(args_info.molFraction_arg[i]);
324 }
else if ((
int)args_info.molFraction_given == nComponents - 1) {
325 RealType remainingFraction = 1.0;
326 for (
int i = 0; i < nComponents - 1; i++) {
327 molFractions.push_back(args_info.molFraction_arg[i]);
328 remainingFraction -= molFractions[i];
330 molFractions.push_back(remainingFraction);
332 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
333 "nanorodBuilder can't figure out molFractions "
334 "for all of the components in the <MetaData> block.");
335 painCave.isFatal = 1;
338 }
else if ((
int)args_info.shellRadius_given) {
339 if ((
int)args_info.shellRadius_given == nComponents) {
340 for (
int i = 0; i < nComponents; i++) {
341 shellRadii.push_back(args_info.shellRadius_arg[i]);
343 }
else if ((
int)args_info.shellRadius_given == nComponents - 1) {
344 for (
int i = 0; i < nComponents - 1; i++) {
345 shellRadii.push_back(args_info.shellRadius_arg[i]);
347 shellRadii.push_back(rodRadius);
350 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
351 "nanorodBuilder can't figure out the\n"
352 "\tshell radii for all of the components in the <MetaData> block.");
353 painCave.isFatal = 1;
357 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
358 "You have a multi-component <MetaData> block,\n"
359 "\tbut have not specified either molFraction or shellRadius "
361 painCave.isFatal = 1;
365 if (args_info.molFraction_given) {
366 RealType totalFraction = 0.0;
370 for (
int i = 0; i < nComponents; i++) {
371 if (molFractions.at(i) < 0.0) {
372 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
373 "One of the requested molFractions was"
375 painCave.isFatal = 1;
378 if (molFractions.at(i) > 1.0) {
379 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
380 "One of the requested molFractions was"
381 " greater than one!");
382 painCave.isFatal = 1;
385 totalFraction += molFractions.at(i);
387 if (abs(totalFraction - 1.0) > 1e-6) {
388 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
389 "The sum of molFractions was not close enough to 1.0");
390 painCave.isFatal = 1;
394 int remaining = nSites;
395 for (
int i = 0; i < nComponents - 1; i++) {
396 nMol.push_back(
int((RealType)nSites * molFractions.at(i)));
397 remaining -= nMol.at(i);
399 nMol.push_back(remaining);
403 int totalMolecules = 0;
404 for (
int i = 0; i < nComponents; i++) {
405 molFractions[i] = (RealType)(nMol.at(i)) / (RealType)nSites;
406 totalMolecules += nMol.at(i);
408 if (totalMolecules != nSites) {
409 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
410 "Computed total number of molecules is not equal "
411 "to the number of lattice sites!");
412 painCave.isFatal = 1;
416 for (
unsigned int i = 0; i < shellRadii.size(); i++) {
417 if (shellRadii.at(i) > rodRadius + 1e-6) {
418 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
419 "One of the shellRadius values exceeds the rod Radius.");
420 painCave.isFatal = 1;
423 if (shellRadii.at(i) <= 0.0) {
424 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
425 "One of the shellRadius values is smaller than zero!");
426 painCave.isFatal = 1;
433 if ((
int)args_info.molFraction_given) {
435 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
436 "Creating a randomized spherically-capped nanorod.");
437 painCave.isFatal = 0;
438 painCave.severity = OPENMD_INFO;
442 for (
unsigned int i = 0; i < sites.size(); i++)
443 if (!isVacancy[i]) ids.push_back(i);
445 std::shuffle(ids.begin(), ids.end(), gen);
448 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
449 "Creating an fcc nanorod.");
450 painCave.isFatal = 0;
451 painCave.severity = OPENMD_INFO;
455 int myComponent = -1;
457 nMol.resize(nComponents);
462 for (
unsigned int i = 0; i < sites.size(); i++) {
464 myR = myLoc.length();
478 componentFromSite[i] = myComponent;
488 outputFileName = args_info.output_arg;
492 createMdFile(inputFileName, outputFileName, nMol);
497 SimInfo* NewInfo = newCreator.createSim(outputFileName,
false);
501 SimInfo::MoleculeIterator mi;
502 mol = NewInfo->beginMolecule(mi);
506 for (
int i = 0; i < nComponents; i++) {
508 new MoLocator(NewInfo->getMoleculeStamp(i), NewInfo->getForceField());
511 if (!args_info.molFraction_given) {
512 for (
unsigned int n = 0; n < sites.size(); n++) {
514 if (componentFromSite[n] == i) {
515 mol = NewInfo->getMoleculeByGlobalIndex(l);
516 locator->placeMol(sites[n], orientations[n], mol);
522 for (
int n = 0; n < nMol.at(i); n++) {
523 mol = NewInfo->getMoleculeByGlobalIndex(l);
524 locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
531 hmat(0, 0) = 10.0 * rodRadius;
536 hmat(1, 1) = 10.0 * rodRadius;
541 hmat(2, 2) = 5.0 * rodLength + 2.0 * rodRadius;
544 NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
547 writer =
new DumpWriter(NewInfo, outputFileName);
549 if (writer == NULL) {
550 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
551 "Error in creating dumpwriter object ");
552 painCave.isFatal = 1;
563 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
564 "A new OpenMD file called \"%s\" has been "
566 outputFileName.c_str());
567 painCave.isFatal = 0;
568 painCave.severity = OPENMD_INFO;
573void createMdFile(
const std::string& oldMdFileName,
574 const std::string& newMdFileName, std::vector<int> nMol) {
577 const int MAXLEN = 65535;
581 oldMdFile.open(oldMdFileName.c_str());
582 newMdFile.open(newMdFileName.c_str());
583 oldMdFile.getline(buffer, MAXLEN);
586 while (!oldMdFile.eof()) {
588 if (strstr(buffer,
"nMol") != NULL) {
589 if (i < nMol.size()) {
590 snprintf(buffer, MAXLEN,
"\tnMol = %i;", nMol.at(i));
591 newMdFile << buffer << std::endl;
595 newMdFile << buffer << std::endl;
597 oldMdFile.getline(buffer, MAXLEN);
603 if (i != nMol.size()) {
604 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
605 "Couldn't replace the correct number of nMol\n"
606 "\tstatements in component blocks. Make sure that all\n"
607 "\tcomponents in the template file have nMol=1");
608 painCave.isFatal = 1;
The only responsibility of SimCreator is to parse the meta-data file and create a SimInfo instance ba...
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Implements a spherically-capped rod-shaped lattice.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
void registerLattice()
Register all lattice.
The header file for the command line option parser generated by GNU Gengetopt version 2....
Where the command line options are stored.