ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/simpleBuilder/simpleBuilder.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/applications/simpleBuilder/simpleBuilder.cpp (file contents):
Revision 2181 by tim, Tue Apr 12 21:58:09 2005 UTC vs.
Revision 3041 by gezelter, Tue Oct 10 18:34:12 2006 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 62 | Line 62 | void createMdFile(const std::string&oldMdFileName, con
62  
63   using namespace std;
64   using namespace oopse;
65 – void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
66 –                  int numMol);
65  
66 + void createMdFile(const std::string&oldMdFileName,
67 +                  const std::string&newMdFileName,
68 +                  int nMol);
69 +
70   int main(int argc, char *argv []) {
71  
72 <    //register force fields
73 <    registerForceFields();
74 <    registerLattice();
72 >  // register force fields
73 >  registerForceFields();
74 >  registerLattice();
75      
76 <    gengetopt_args_info args_info;
77 <    std::string latticeType;
78 <    std::string inputFileName;
79 <    std::string outPrefix;
80 <    std::string outMdFileName;
81 <    std::string outInitFileName;
82 <    BaseLattice *simpleLat;
83 <    int numMol;
84 <    double latticeConstant;
85 <    std::vector<double> lc;
86 <    double mass;
87 <    const double rhoConvertConst = 1.661;
88 <    double density;
89 <    int nx,
90 <    ny,
91 <    nz;
90 <    Mat3x3d hmat;
91 <    MoLocator *locator;
92 <    std::vector<Vector3d> latticePos;
93 <    std::vector<Vector3d> latticeOrt;
94 <    int numMolPerCell;
95 <    int curMolIndex;
96 <    DumpWriter *writer;
76 >  gengetopt_args_info args_info;
77 >  std::string latticeType;
78 >  std::string inputFileName;
79 >  std::string outputFileName;
80 >  Lattice *simpleLat;
81 >  RealType latticeConstant;
82 >  std::vector<RealType> lc;
83 >  const RealType rhoConvertConst = 1.661;
84 >  RealType density;
85 >  int nx, ny, nz;
86 >  Mat3x3d hmat;
87 >  MoLocator *locator;
88 >  std::vector<Vector3d> latticePos;
89 >  std::vector<Vector3d> latticeOrt;
90 >  int nMolPerCell;
91 >  DumpWriter *writer;
92  
93 <    // parse command line arguments
94 <    if (cmdline_parser(argc, argv, &args_info) != 0)
95 <        exit(1);
93 >  // parse command line arguments
94 >  if (cmdline_parser(argc, argv, &args_info) != 0)
95 >    exit(1);
96  
97 <    density = args_info.density_arg;
97 >  density = args_info.density_arg;
98  
99 <    //get lattice type
100 <    latticeType = UpperCase(args_info.latticetype_arg);
99 >  //get lattice type
100 >  latticeType = "FCC";
101  
102 <    if (!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)) {
103 <        std::cerr << latticeType << " is an invalid lattice type" << std::endl;
104 <        std::cerr << LatticeFactory::getInstance()->toString() << std::endl;
105 <        exit(1);
106 <    }
102 >  simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
103 >    
104 >  if (simpleLat == NULL) {
105 >    sprintf(painCave.errMsg, "Lattice Factory can not create %s lattice\n",
106 >            latticeType.c_str());
107 >    painCave.isFatal = 1;
108 >    simError();
109 >  }
110 >  nMolPerCell = simpleLat->getNumSitesPerCell();
111  
112 <    //get the number of unit cell
114 <    nx = args_info.nx_arg;
112 >  //get the number of unit cells in each direction:
113  
114 <    if (nx <= 0) {
117 <        std::cerr << "The number of unit cell in h direction must be greater than 0" << std::endl;
118 <        exit(1);
119 <    }
114 >  nx = args_info.nx_arg;
115  
116 <    ny = args_info.ny_arg;
116 >  if (nx <= 0) {
117 >    sprintf(painCave.errMsg, "The number of unit cells in the x direction "
118 >            "must be greater than 0.");
119 >    painCave.isFatal = 1;
120 >    simError();
121 >  }
122  
123 <    if (ny <= 0) {
124 <        std::cerr << "The number of unit cell in l direction must be greater than 0" << std::endl;
125 <        exit(1);
126 <    }
123 >  ny = args_info.ny_arg;
124  
125 <    nz = args_info.nz_arg;
125 >  if (ny <= 0) {
126 >    sprintf(painCave.errMsg, "The number of unit cells in the y direction "
127 >            "must be greater than 0.");
128 >    painCave.isFatal = 1;
129 >    simError();
130 >  }
131  
132 <    if (nz <= 0) {
131 <        std::cerr << "The number of unit cell in k direction must be greater than 0" << std::endl;
132 <        exit(1);
133 <    }
132 >  nz = args_info.nz_arg;
133  
134 <    //get input file name
135 <    if (args_info.inputs_num)
136 <        inputFileName = args_info.inputs[0];
137 <    else {
138 <        std::cerr << "You must specify a input file name.\n" << std::endl;
139 <        cmdline_parser_print_help();
141 <        exit(1);
142 <    }
134 >  if (nz <= 0) {
135 >    sprintf(painCave.errMsg, "The number of unit cells in the z direction "
136 >            "must be greater than 0.");
137 >    painCave.isFatal = 1;
138 >    simError();
139 >  }
140  
141 <    //parse md file and set up the system
145 <    SimCreator oldCreator;
146 <    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
141 >  int nSites = nMolPerCell * nx * ny * nz;
142  
143 <    if (oldInfo->getNMoleculeStamp()>= 2) {
144 <        std::cerr << "can not build the system with more than two components"
145 <            << std::endl;
146 <        exit(1);
147 <    }
143 >  //get input file name
144 >  if (args_info.inputs_num)
145 >    inputFileName = args_info.inputs[0];
146 >  else {
147 >    sprintf(painCave.errMsg, "No input .md file name was specified "
148 >            "on the command line");
149 >    painCave.isFatal = 1;
150 >    simError();
151 >  }
152  
153 <    //get mass of molecule.
153 >  //parse md file and set up the system
154  
155 <    mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
155 >  SimCreator oldCreator;
156 >  SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
157 >  Globals* simParams = oldInfo->getSimParams();
158  
159 <    //creat lattice
159 <    simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
159 >  // Calculate lattice constant (in Angstroms)
160  
161 <    if (simpleLat == NULL) {
162 <        std::cerr << "Error in creating lattice" << std::endl;
163 <        exit(1);
164 <    }
161 >  RealType avgMass = getMolMass(oldInfo->getMoleculeStamp(0),
162 >                                  oldInfo->getForceField());
163  
164 <    numMolPerCell = simpleLat->getNumSitesPerCell();
164 >  latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
165 >                        (RealType)(1.0 / 3.0));
166 >  
167 >  // Set the lattice constant
168 >  
169 >  lc.push_back(latticeConstant);
170 >  simpleLat->setLatticeConstant(lc);
171  
172 <    //calculate lattice constant (in Angstrom)
169 <    latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
170 <                          1.0 / 3.0);
172 >  // Calculate the lattice sites and fill the lattice vector.
173  
174 <    //set lattice constant
173 <    lc.push_back(latticeConstant);
174 <    simpleLat->setLatticeConstant(lc);
174 >  // Get the standard orientations of the cell sites
175  
176 <    //calculate the total number of molecules
177 <    numMol = nx * ny * nz * numMolPerCell;
176 >  latticeOrt = simpleLat->getLatticePointsOrt();
177  
178 <    if (oldInfo->getNGlobalMolecules() != numMol) {
179 <        outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
180 <        outMdFileName = outPrefix + ".md";
178 >  vector<Vector3d> sites;
179 >  vector<Vector3d> orientations;
180 >  
181 >  for(int i = 0; i < nx; i++) {
182 >    for(int j = 0; j < ny; j++) {
183 >      for(int k = 0; k < nz; k++) {
184  
185 <        //creat new .md file on fly which corrects the number of molecule    
186 <        createMdFile(inputFileName, outMdFileName, numMol);
187 <        std::cerr
188 <            << "SimpleBuilder Error: the number of molecule and the density are not matched"
189 <            << std::endl;
190 <        std::cerr << "A new .md file: " << outMdFileName
191 <            << " is generated, use it to rerun the simpleBuilder" << std::endl;
192 <        exit(1);
185 >        // Get the position of the cell sites
186 >        
187 >        simpleLat->getLatticePointsPos(latticePos, i, j, k);
188 >        
189 >        for(int l = 0; l < nMolPerCell; l++) {
190 >          sites.push_back(latticePos[l]);
191 >          orientations.push_back(latticeOrt[l]);
192 >        }
193 >      }
194      }
195 +  }
196 +  
197 +  outputFileName = args_info.output_arg;
198 +  
199 +  // create a new .md file on the fly which corrects the number of molecules
200  
201 <    //determine the output file names  
194 <    if (args_info.output_given)
195 <        outInitFileName = args_info.output_arg;
196 <    else
197 <        outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
198 <    
199 <    //creat Molocator
200 <    locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
201 >  createMdFile(inputFileName, outputFileName, nSites);
202  
203 <    //fill Hmat
204 <    hmat(0, 0)= nx * latticeConstant;
204 <    hmat(0, 1) = 0.0;
205 <    hmat(0, 2) = 0.0;
203 >  if (oldInfo != NULL)
204 >    delete oldInfo;
205  
206 <    hmat(1, 0) = 0.0;
207 <    hmat(1, 1) = ny * latticeConstant;
209 <    hmat(1, 2) = 0.0;
206 >  // We need to read in the new SimInfo object, then Parse the
207 >  // md file and set up the system
208  
209 <    hmat(2, 0) = 0.0;
210 <    hmat(2, 1) = 0.0;
213 <    hmat(2, 2) = nz * latticeConstant;
209 >  SimCreator newCreator;
210 >  SimInfo* newInfo = newCreator.createSim(outputFileName, false);
211  
212 <    //set Hmat
216 <    oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
212 >  // fill Hmat
213  
214 <    //place the molecules
214 >  hmat(0, 0) = nx * latticeConstant;
215 >  hmat(0, 1) = 0.0;
216 >  hmat(0, 2) = 0.0;
217  
218 <    curMolIndex = 0;
218 >  hmat(1, 0) = 0.0;
219 >  hmat(1, 1) = ny * latticeConstant;
220 >  hmat(1, 2) = 0.0;
221  
222 <    //get the orientation of the cell sites
223 <    //for the same type of molecule in same lattice, it will not change
224 <    latticeOrt = simpleLat->getLatticePointsOrt();
222 >  hmat(2, 0) = 0.0;
223 >  hmat(2, 1) = 0.0;
224 >  hmat(2, 2) = nz * latticeConstant;
225  
226 <    Molecule* mol;
227 <    SimInfo::MoleculeIterator mi;
228 <    mol = oldInfo->beginMolecule(mi);
229 <    for(int i = 0; i < nx; i++) {
230 <        for(int j = 0; j < ny; j++) {
231 <            for(int k = 0; k < nz; k++) {
226 >  // Set Hmat
227  
228 <                //get the position of the cell sites
234 <                simpleLat->getLatticePointsPos(latticePos, i, j, k);
228 >  newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
229  
230 <                for(int l = 0; l < numMolPerCell; l++) {
231 <                    if (mol != NULL) {
232 <                        locator->placeMol(latticePos[l], latticeOrt[l], mol);
233 <                    } else {
234 <                        std::cerr << std::endl;                    
235 <                    }
236 <                    mol = oldInfo->nextMolecule(mi);
237 <                }
238 <            }
239 <        }
240 <    }
230 >  // place the molecules
231 >  
232 >  Molecule* mol;
233 >  locator = new MoLocator(newInfo->getMoleculeStamp(0),
234 >                          newInfo->getForceField());
235 >  for (int n = 0; n < nSites; n++) {
236 >    mol = newInfo->getMoleculeByGlobalIndex(n);
237 >    locator->placeMol(sites[n], orientations[n], mol);
238 >  }
239 >  
240 >  // Create DumpWriter and write out the coordinates
241  
242 <    //create dumpwriter and write out the coordinates
243 <    oldInfo->setFinalConfigFileName(outInitFileName);
244 <    writer = new DumpWriter(oldInfo);
242 >  writer = new DumpWriter(newInfo, outputFileName);
243 >  
244 >  if (writer == NULL) {
245 >    sprintf(painCave.errMsg, "error in creating DumpWriter");
246 >    painCave.isFatal = 1;
247 >    simError();
248 >  }
249  
250 <    if (writer == NULL) {
253 <        std::cerr << "error in creating DumpWriter" << std::endl;
254 <        exit(1);
255 <    }
250 >  writer->writeDump();
251  
252 <    writer->writeDump();
258 <    std::cout << "new initial configuration file: " << outInitFileName
259 <        << " is generated." << std::endl;
252 >  // deleting the writer will put the closing at the end of the dump file.
253  
254 <    //delete objects
254 >  delete writer;
255  
256 <    //delete oldInfo and oldSimSetup
257 <    if (oldInfo != NULL)
258 <        delete oldInfo;
259 <
260 <    if (writer != NULL)
268 <        delete writer;
269 <
270 <    return 0;
256 >  sprintf(painCave.errMsg, "A new OOPSE MD file called \"%s\" has been "
257 >          "generated.\n", outputFileName.c_str());
258 >  painCave.isFatal = 0;
259 >  simError();
260 >  return 0;
261   }
262  
263 < void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
264 <                  int numMol) {
265 <    ifstream oldMdFile;
266 <    ofstream newMdFile;
267 <    const int MAXLEN = 65535;
268 <    char buffer[MAXLEN];
263 > void createMdFile(const std::string&oldMdFileName,
264 >                  const std::string&newMdFileName,
265 >                  int nMol) {
266 >  ifstream oldMdFile;
267 >  ofstream newMdFile;
268 >  const int MAXLEN = 65535;
269 >  char buffer[MAXLEN];
270  
271 <    //create new .md file based on old .md file
272 <    oldMdFile.open(oldMdFileName.c_str());
273 <    newMdFile.open(newMdFileName.c_str());
271 >  //create new .md file based on old .md file
272 >  oldMdFile.open(oldMdFileName.c_str());
273 >  newMdFile.open(newMdFileName.c_str());
274  
275 <    oldMdFile.getline(buffer, MAXLEN);
275 >  oldMdFile.getline(buffer, MAXLEN);
276  
277 <    while (!oldMdFile.eof()) {
277 >  while (!oldMdFile.eof()) {
278  
279 <        //correct molecule number
280 <        if (strstr(buffer, "nMol") != NULL) {
281 <            sprintf(buffer, "\t\tnMol = %d;", numMol);
282 <            newMdFile << buffer << std::endl;
283 <        } else
284 <            newMdFile << buffer << std::endl;
279 >    //correct molecule number
280 >    if (strstr(buffer, "nMol") != NULL) {
281 >      sprintf(buffer, "\t\tnMol = %d;", nMol);
282 >      newMdFile << buffer << std::endl;
283 >    } else
284 >      newMdFile << buffer << std::endl;
285  
286 <        oldMdFile.getline(buffer, MAXLEN);
287 <    }
286 >    oldMdFile.getline(buffer, MAXLEN);
287 >  }
288  
289 <    oldMdFile.close();
290 <    newMdFile.close();
289 >  oldMdFile.close();
290 >  newMdFile.close();
291   }
292  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines