50#include "brains/Register.hpp"
53#include "brains/Thermo.hpp"
55#include "io/DumpWriter.hpp"
58#include "types/FluctuatingChargeAdapter.hpp"
59#include "utils/Constants.hpp"
60#include "utils/simError.h"
66void createMdFile(
const std::string& oldMdFileName,
67 const std::string& newMdFileName, std::vector<int> nMol);
69int main(
int argc,
char* argv[]) {
75 if (cmdline_parser(argc, argv, &args_info) != 0) { exit(1); }
78 if (args_info.input_given) {
79 dumpFileName = args_info.input_arg;
81 strcpy(painCave.errMsg,
"No input file name was specified.\n");
86 if (args_info.output_given) {
87 outFileName = args_info.output_arg;
89 strcpy(painCave.errMsg,
"No output file name was specified.\n");
95 double phi = args_info.rotatePhi_arg * (Constants::PI / 180.0);
96 double theta = args_info.rotateTheta_arg * (Constants::PI / 180.0);
97 double psi = args_info.rotatePsi_arg * (Constants::PI / 180.0);
101 rotMatrix.setupRotMat(phi, theta, psi);
104 args_info.repeatZ_arg);
107 repeatD(0, 0) = repeat.x();
108 repeatD(1, 1) = repeat.y();
109 repeatD(2, 2) = repeat.z();
112 Vector3d(args_info.translateX_arg, args_info.translateY_arg,
113 args_info.translateZ_arg);
118 SimInfo* oldInfo = oldCreator.createSim(dumpFileName,
false);
119 Globals* simParams = oldInfo->getSimParams();
120 std::vector<Component*> components = simParams->getComponents();
121 std::vector<int> nMol;
122 for (vector<Component*>::iterator i = components.begin();
123 i != components.end(); ++i) {
124 int nMolOld = (*i)->getNMol();
125 int nMolNew = nMolOld * repeat.x() * repeat.y() * repeat.z();
126 nMol.push_back(nMolNew);
129 createMdFile(dumpFileName, outFileName, nMol);
132 SimInfo* newInfo = newCreator.createSim(outFileName,
false);
135 int nframes = dumpReader->getNFrames();
138 if (writer == NULL) {
139 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
140 "error in creating DumpWriter");
141 painCave.isFatal = 1;
145 SimInfo::MoleculeIterator miter;
146 Molecule::IntegrableObjectIterator iiter;
147 Molecule::RigidBodyIterator rbIter;
165 for (
int i = 0; i < nframes; i++) {
166 cerr <<
"frame = " << i <<
"\n";
167 dumpReader->readFrame(i);
168 oldSnap = oldInfo->getSnapshotManager()->getCurrentSnapshot();
169 newSnap = newInfo->getSnapshotManager()->getCurrentSnapshot();
171 newSnap->setID(oldSnap->getID());
172 newSnap->setTime(oldSnap->getTime());
174 oldHmat = oldSnap->getHmat();
175 rotHmat = rotMatrix * oldHmat;
176 newHmat = repeatD * rotHmat;
177 newSnap->setHmat(newHmat);
179 newSnap->setThermostat(oldSnap->getThermostat());
180 newSnap->setBarostat(oldSnap->getBarostat());
182 COM = thermo.getCom();
185 for (mol = oldInfo->beginMolecule(miter); mol != NULL;
186 mol = oldInfo->nextMolecule(miter)) {
187 if (args_info.repairMolecules_arg == 1) { molCOM = mol->
getCom(); }
189 for (
int ii = 0; ii < repeat.x(); ii++) {
190 for (
int jj = 0; jj < repeat.y(); jj++) {
191 for (
int kk = 0; kk < repeat.z(); kk++) {
193 for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
194 sd = mol->nextIntegrableObject(iiter)) {
195 if (args_info.repairMolecules_arg == 1) {
196 relPos = sd->
getPos() - molCOM;
197 oldPos = molCOM - COM + translate;
198 oldSnap->wrapVector(relPos);
199 oldSnap->wrapVector(oldPos);
202 oldPos = sd->
getPos() - COM + translate;
203 oldSnap->wrapVector(oldPos);
206 newPos = rotMatrix * oldPos + trans * oldHmat;
207 sdNew = newInfo->getIOIndexToIntegrableObject(newIndex);
208 sdNew->setPos(newPos);
209 sdNew->setVel(rotMatrix * sd->
getVel());
213 bodyRotMat = bodyRotMat * rotMatrix.
inverse();
214 sdNew->setA(bodyRotMat);
216 sdNew->setJ(rotMatrix * sd->
getJ());
220 atype =
static_cast<Atom*
>(sd)->getAtomType();
222 if (fqa.isFluctuatingCharge()) {
224 sdNew->setFlucQPos(charge);
226 sdNew->setFlucQVel(cv);
238 for (mol = newInfo->beginMolecule(miter); mol != NULL;
239 mol = newInfo->nextMolecule(miter)) {
241 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
242 rb = mol->nextRigidBody(rbIter)) {
255void createMdFile(
const std::string& oldMdFileName,
256 const std::string& newMdFileName, std::vector<int> nMol) {
259 const int MAXLEN = 65535;
264 oldMdFile.open(oldMdFileName.c_str());
265 newMdFile.open(newMdFileName.c_str());
267 oldMdFile.getline(buffer, MAXLEN);
270 while (!oldMdFile.eof()) {
272 if (strstr(buffer,
"nMol") != NULL) {
273 if (i < nMol.size()) {
274 snprintf(buffer, MAXLEN,
"\tnMol = %i;", nMol.at(i));
275 newMdFile << buffer << std::endl;
279 newMdFile << buffer << std::endl;
281 oldMdFile.getline(buffer, MAXLEN);
287 if (i != nMol.size()) {
288 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
289 "Couldn't replace the correct number of nMol\n"
290 "\tstatements in component blocks.");
291 painCave.isFatal = 1;
AtomType is what OpenMD looks to for unchanging data about an atom.
Vector3d getCom()
Returns the current center of mass position of this molecule.
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...
The Snapshot class is a repository storing dynamic data during a Simulation.
SquareMatrix3< Real > inverse() const
Sets the value of this matrix to the inverse of itself.
"Don't move, or you're dead! Stand up! Captain, we've got them!"
RotMat3x3d getA()
Returns the current rotation matrix of this stuntDouble.
Vector3d getVel()
Returns the current velocity of this stuntDouble.
Vector3d getPos()
Returns the current position of this stuntDouble.
RealType getFlucQPos()
Returns the current fluctuating charge of this stuntDouble.
Vector3d getJ()
Returns the current angular momentum of this stuntDouble (body -fixed).
bool isAtom()
Tests if this stuntDouble is an atom.
bool isDirectional()
Tests if this stuntDouble is a directional one.
RealType getFlucQVel()
Returns the current charge velocity of this stuntDouble.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
The header file for the command line option parser generated by GNU Gengetopt version 2....
Where the command line options are stored.