55#include "brains/Register.hpp"
58#include "io/DumpWriter.hpp"
59#include "lattice/Lattice.hpp"
64#include "utils/MoLocator.hpp"
70void createMdFile(
const std::string& oldMdFileName,
71 const std::string& newMdFileName, std::vector<int> nMol);
73int main(
int argc,
char* argv[]) {
77 std::string latticeType;
78 std::string inputFileName;
79 std::string outputFileName;
81 RealType latticeConstant;
82 std::vector<RealType> lc;
83 const RealType rhoConvertConst = 1.661;
88 std::vector<Vector3d> latticePos;
89 std::vector<Vector3d> latticeOrt;
94 if (cmdline_parser(argc, argv, &args_info) != 0) exit(1);
96 density = args_info.density_arg;
100 if (args_info.lattice_given) { latticeType = args_info.lattice_arg; }
104 if (simpleLat == NULL) {
105 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
106 "Lattice Factory can not create %s lattice\n",
107 latticeType.c_str());
108 painCave.isFatal = 1;
111 nMolPerCell = simpleLat->getNumSitesPerCell();
115 nx = args_info.nx_arg;
118 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
119 "The number of unit cells in the x direction "
120 "must be greater than 0.");
121 painCave.isFatal = 1;
125 ny = args_info.ny_arg;
128 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
129 "The number of unit cells in the y direction "
130 "must be greater than 0.");
131 painCave.isFatal = 1;
135 nz = args_info.nz_arg;
138 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
139 "The number of unit cells in the z direction "
140 "must be greater than 0.");
141 painCave.isFatal = 1;
145 int nSites = nMolPerCell * nx * ny * nz;
148 if (args_info.inputs_num)
149 inputFileName = args_info.inputs[0];
151 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
152 "No input .omd file name was specified "
153 "on the command line");
154 painCave.isFatal = 1;
161 SimInfo* oldInfo = oldCreator.createSim(inputFileName,
false);
162 Globals* simParams = oldInfo->getSimParams();
166 std::vector<Component*> components = simParams->getComponents();
167 std::vector<RealType> molFractions;
168 std::vector<RealType> molecularMasses;
169 std::vector<int> nMol;
170 std::size_t nComponents = components.size();
172 if (nComponents == 1) {
173 molFractions.push_back(1.0);
175 if (args_info.molFraction_given == nComponents) {
176 for (std::size_t i = 0; i < nComponents; i++) {
177 molFractions.push_back(args_info.molFraction_arg[i]);
179 }
else if (args_info.molFraction_given == nComponents - 1) {
180 RealType remainingFraction = 1.0;
181 for (std::size_t i = 0; i < nComponents - 1; i++) {
182 molFractions.push_back(args_info.molFraction_arg[i]);
183 remainingFraction -= molFractions[i];
185 molFractions.push_back(remainingFraction);
187 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
188 "randomBuilder can't figure out molFractions "
189 "for all of the components in the <MetaData> block.");
190 painCave.isFatal = 1;
197 RealType totalFraction = 0.0;
199 for (std::size_t i = 0; i < nComponents; i++) {
200 if (molFractions.at(i) < 0.0) {
201 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
202 "One of the requested molFractions was"
204 painCave.isFatal = 1;
207 if (molFractions.at(i) > 1.0) {
208 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
209 "One of the requested molFractions was"
210 " greater than one!");
211 painCave.isFatal = 1;
214 totalFraction += molFractions.at(i);
216 if (abs(totalFraction - 1.0) > 1e-6) {
217 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
218 "The sum of molFractions was not close enough to 1.0");
219 painCave.isFatal = 1;
223 int remaining = nSites;
224 for (std::size_t i = 0; i < nComponents - 1; i++) {
225 nMol.push_back(
int((RealType)nSites * molFractions.at(i)));
226 remaining -= nMol.at(i);
228 nMol.push_back(remaining);
232 int totalMolecules = 0;
233 RealType totalMass = 0.0;
234 for (std::size_t i = 0; i < nComponents; i++) {
235 molFractions[i] = (RealType)(nMol.at(i)) / (RealType)nSites;
236 totalMolecules += nMol.at(i);
237 molecularMasses.push_back(MoLocator::getMolMass(
238 oldInfo->getMoleculeStamp(i), oldInfo->getForceField()));
239 totalMass += (RealType)(nMol.at(i)) * molecularMasses.at(i);
241 RealType avgMass = totalMass / (RealType)totalMolecules;
243 if (totalMolecules != nSites) {
244 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
245 "Computed total number of molecules is not equal "
246 "to the number of lattice sites!");
247 painCave.isFatal = 1;
251 latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
252 (RealType)(1.0 / 3.0));
256 lc.push_back(latticeConstant);
257 simpleLat->setLatticeConstant(lc);
263 latticeOrt = simpleLat->getLatticePointsOrt();
265 vector<Vector3d> sites;
266 vector<Vector3d> orientations;
268 for (
int i = 0; i < nx; i++) {
269 for (
int j = 0; j < ny; j++) {
270 for (
int k = 0; k < nz; k++) {
273 simpleLat->getLatticePointsPos(latticePos, i, j, k);
275 for (
int l = 0; l < nMolPerCell; l++) {
276 sites.push_back(latticePos[l]);
277 orientations.push_back(latticeOrt[l]);
283 outputFileName = args_info.output_arg;
287 createMdFile(inputFileName, outputFileName, nMol);
295 SimInfo* newInfo = newCreator.createSim(outputFileName,
false);
299 hmat(0, 0) = nx * latticeConstant;
304 hmat(1, 1) = ny * latticeConstant;
309 hmat(2, 2) = nz * latticeConstant;
313 newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
320 for (std::size_t i = 0; i < sites.size(); i++)
324 std::random_device rd;
326 std::mt19937 gen(rd());
328 std::shuffle(ids.begin(), ids.end(), gen);
332 for (std::size_t i = 0; i < nComponents; i++) {
334 new MoLocator(newInfo->getMoleculeStamp(i), newInfo->getForceField());
335 for (
int n = 0; n < nMol.at(i); n++) {
336 mol = newInfo->getMoleculeByGlobalIndex(l);
337 locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
344 writer =
new DumpWriter(newInfo, outputFileName);
346 if (writer == NULL) {
347 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
348 "error in creating DumpWriter");
349 painCave.isFatal = 1;
359 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
360 "A new OpenMD file called \"%s\" has been "
362 outputFileName.c_str());
363 painCave.isFatal = 0;
364 painCave.severity = OPENMD_INFO;
369void createMdFile(
const std::string& oldMdFileName,
370 const std::string& newMdFileName, std::vector<int> nMol) {
373 const int MAXLEN = 65535;
378 oldMdFile.open(oldMdFileName.c_str());
379 newMdFile.open(newMdFileName.c_str());
381 oldMdFile.getline(buffer, MAXLEN);
384 while (!oldMdFile.eof()) {
386 if (strstr(buffer,
"nMol") != NULL) {
387 if (i < nMol.size()) {
388 snprintf(buffer, MAXLEN,
"\tnMol = %i;", nMol.at(i));
389 newMdFile << buffer << std::endl;
393 newMdFile << buffer << std::endl;
395 oldMdFile.getline(buffer, MAXLEN);
401 if (i != nMol.size()) {
402 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
403 "Couldn't replace the correct number of nMol\n"
404 "\tstatements in component blocks. Make sure that all\n"
405 "\tcomponents in the template file have nMol=1");
406 painCave.isFatal = 1;
static LatticeFactory & getInstance()
Returns an instance of Lattice factory.
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...
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.