58#include "brains/Register.hpp"
61#include "io/DumpWriter.hpp"
62#include "lattice/Lattice.hpp"
67#include "shapedLatticeEllipsoid.hpp"
68#include "shapedLatticeRod.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();
119 vector<Vector3d> sites;
120 vector<Vector3d> orientations;
122 if (args_info.ellipsoid_flag) {
124 rodLength, rodRadius);
125 sites = nanoEllipsoid.getSites();
126 orientations = nanoEllipsoid.getOrientations();
132 sites = nanoRod.getSites();
133 orientations = nanoRod.getOrientations();
137 std::random_device rd;
139 std::mt19937 gen(rd());
141 std::vector<std::size_t> vacancyTargets;
142 vector<bool> isVacancy;
147 for (std::size_t i = 0; i < sites.size(); i++)
148 isVacancy.push_back(
false);
151 if (args_info.vacancyPercent_given) {
153 if (args_info.vacancyPercent_arg < 0.0 ||
154 args_info.vacancyPercent_arg > 100.0) {
155 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
156 "vacancyPercent was set to a non-sensical value.");
157 painCave.isFatal = 1;
160 RealType vF = args_info.vacancyPercent_arg / 100.0;
164 if (args_info.vacancyInnerRadius_given) {
165 vIR = args_info.vacancyInnerRadius_arg;
169 if (args_info.vacancyOuterRadius_given) {
170 vOR = args_info.vacancyOuterRadius_arg;
174 if (vIR >= 0.0 && vOR <= rodRadius && vOR >= vIR) {
175 for (std::size_t i = 0; i < sites.size(); i++) {
177 myR = myLoc.length();
178 if (myR >= vIR && myR <= vOR) { vacancyTargets.push_back(i); }
180 std::shuffle(vacancyTargets.begin(), vacancyTargets.end(), gen);
182 std::size_t nTargets = vacancyTargets.size();
183 vacancyTargets.resize((
int)(vF * nTargets));
185 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
186 "Removing %d atoms from randomly-selected\n"
187 "\tsites between %lf and %lf.",
188 (
int)vacancyTargets.size(), vIR, vOR);
189 painCave.isFatal = 0;
190 painCave.severity = OPENMD_INFO;
194 for (std::size_t i = 0; i < sites.size(); i++) {
196 for (std::size_t j = 0; j < vacancyTargets.size(); j++) {
197 if (i == vacancyTargets[j]) vac =
true;
199 isVacancy.push_back(vac);
203 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
204 "Something is strange about the vacancy\n"
205 "\tinner or outer radii. Check their values.");
206 painCave.isFatal = 1;
213 std::size_t nSites = sites.size() - vacancyTargets.size();
219 std::vector<Component*> components = simParams->getComponents();
220 std::vector<RealType> molFractions;
221 std::vector<RealType> shellRadii;
222 std::vector<int> nMol;
223 std::map<int, int> componentFromSite;
224 nComponents = components.size();
227 if (args_info.molFraction_given && args_info.shellRadius_given) {
228 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
229 "Specify either molFraction or shellRadius "
230 "arguments, but not both!");
231 painCave.isFatal = 1;
235 if (nComponents == 1) {
236 molFractions.push_back(1.0);
237 shellRadii.push_back(rodRadius);
238 }
else if (args_info.molFraction_given) {
239 if ((
int)args_info.molFraction_given == nComponents) {
240 for (
int i = 0; i < nComponents; i++) {
241 molFractions.push_back(args_info.molFraction_arg[i]);
243 }
else if ((
int)args_info.molFraction_given == nComponents - 1) {
244 RealType remainingFraction = 1.0;
245 for (
int i = 0; i < nComponents - 1; i++) {
246 molFractions.push_back(args_info.molFraction_arg[i]);
247 remainingFraction -= molFractions[i];
249 molFractions.push_back(remainingFraction);
251 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
252 "nanorodBuilder can't figure out molFractions "
253 "for all of the components in the <MetaData> block.");
254 painCave.isFatal = 1;
257 }
else if ((
int)args_info.shellRadius_given) {
258 if ((
int)args_info.shellRadius_given == nComponents) {
259 for (
int i = 0; i < nComponents; i++) {
260 shellRadii.push_back(args_info.shellRadius_arg[i]);
262 }
else if ((
int)args_info.shellRadius_given == nComponents - 1) {
263 for (
int i = 0; i < nComponents - 1; i++) {
264 shellRadii.push_back(args_info.shellRadius_arg[i]);
266 shellRadii.push_back(rodRadius);
269 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
270 "nanorodBuilder can't figure out the\n"
271 "\tshell radii for all of the components in the <MetaData> block.");
272 painCave.isFatal = 1;
276 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
277 "You have a multi-component <MetaData> block,\n"
278 "\tbut have not specified either molFraction or shellRadius "
280 painCave.isFatal = 1;
284 if (args_info.molFraction_given) {
285 RealType totalFraction = 0.0;
289 for (
int i = 0; i < nComponents; i++) {
290 if (molFractions.at(i) < 0.0) {
291 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
292 "One of the requested molFractions was"
294 painCave.isFatal = 1;
297 if (molFractions.at(i) > 1.0) {
298 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
299 "One of the requested molFractions was"
300 " greater than one!");
301 painCave.isFatal = 1;
304 totalFraction += molFractions.at(i);
306 if (abs(totalFraction - 1.0) > 1e-6) {
307 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
308 "The sum of molFractions was not close enough to 1.0");
309 painCave.isFatal = 1;
313 int remaining = nSites;
314 for (
int i = 0; i < nComponents - 1; i++) {
315 nMol.push_back(
int((RealType)nSites * molFractions.at(i)));
316 remaining -= nMol.at(i);
318 nMol.push_back(remaining);
322 std::size_t totalMolecules = 0;
323 for (
int i = 0; i < nComponents; i++) {
324 molFractions[i] = (RealType)(nMol.at(i)) / (RealType)nSites;
325 totalMolecules += nMol.at(i);
327 if (totalMolecules != nSites) {
328 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
329 "Computed total number of molecules is not equal "
330 "to the number of lattice sites!");
331 painCave.isFatal = 1;
335 for (
unsigned int i = 0; i < shellRadii.size(); i++) {
336 if (shellRadii.at(i) > rodRadius + 1e-6) {
337 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
338 "One of the shellRadius values exceeds the rod Radius.");
339 painCave.isFatal = 1;
342 if (shellRadii.at(i) <= 0.0) {
343 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
344 "One of the shellRadius values is smaller than zero!");
345 painCave.isFatal = 1;
352 if ((
int)args_info.molFraction_given) {
354 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
355 "Creating a randomized spherically-capped nanorod.");
356 painCave.isFatal = 0;
357 painCave.severity = OPENMD_INFO;
361 for (std::size_t i = 0; i < sites.size(); i++)
362 if (!isVacancy[i]) ids.push_back(i);
364 std::shuffle(ids.begin(), ids.end(), gen);
367 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
368 "Creating an fcc nanorod.");
369 painCave.isFatal = 0;
370 painCave.severity = OPENMD_INFO;
374 int myComponent = -1;
376 nMol.resize(nComponents);
381 for (
unsigned int i = 0; i < sites.size(); i++) {
383 myR = myLoc.length();
397 componentFromSite[i] = myComponent;
407 outputFileName = args_info.output_arg;
411 createMdFile(inputFileName, outputFileName, nMol);
416 SimInfo* NewInfo = newCreator.createSim(outputFileName,
false);
420 SimInfo::MoleculeIterator mi;
421 mol = NewInfo->beginMolecule(mi);
425 for (
int i = 0; i < nComponents; i++) {
427 new MoLocator(NewInfo->getMoleculeStamp(i), NewInfo->getForceField());
430 if (!args_info.molFraction_given) {
431 for (
unsigned int n = 0; n < sites.size(); n++) {
433 if (componentFromSite[n] == i) {
434 mol = NewInfo->getMoleculeByGlobalIndex(l);
435 locator->placeMol(sites[n], orientations[n], mol);
441 for (
int n = 0; n < nMol.at(i); n++) {
442 mol = NewInfo->getMoleculeByGlobalIndex(l);
443 locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
450 hmat(0, 0) = 10.0 * rodRadius;
455 hmat(1, 1) = 10.0 * rodRadius;
460 hmat(2, 2) = 5.0 * rodLength + 2.0 * rodRadius;
463 NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
466 writer =
new DumpWriter(NewInfo, outputFileName);
468 if (writer == NULL) {
469 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
470 "Error in creating dumpwriter object ");
471 painCave.isFatal = 1;
482 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
483 "A new OpenMD file called \"%s\" has been "
485 outputFileName.c_str());
486 painCave.isFatal = 0;
487 painCave.severity = OPENMD_INFO;
492void createMdFile(
const std::string& oldMdFileName,
493 const std::string& newMdFileName, std::vector<int> nMol) {
496 const int MAXLEN = 65535;
500 oldMdFile.open(oldMdFileName.c_str());
501 newMdFile.open(newMdFileName.c_str());
502 oldMdFile.getline(buffer, MAXLEN);
505 while (!oldMdFile.eof()) {
507 if (strstr(buffer,
"nMol") != NULL) {
508 if (i < nMol.size()) {
509 snprintf(buffer, MAXLEN,
"\tnMol = %i;", nMol.at(i));
510 newMdFile << buffer << std::endl;
514 newMdFile << buffer << std::endl;
516 oldMdFile.getline(buffer, MAXLEN);
522 if (i != nMol.size()) {
523 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
524 "Couldn't replace the correct number of nMol\n"
525 "\tstatements in component blocks. Make sure that all\n"
526 "\tcomponents in the template file have nMol=1");
527 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 an ellipsoid-shaped lattice.
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.