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 1674 by tim, Thu Oct 28 19:06:59 2004 UTC vs.
Revision 3041 by gezelter, Tue Oct 10 18:34:12 2006 UTC

# Line 1 | Line 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
5 + * non-exclusive, royalty free, license to use, modify and
6 + * redistribute this software in source and binary code form, provided
7 + * that the following conditions are met:
8 + *
9 + * 1. Acknowledgement of the program authors must be made in any
10 + *    publication of scientific results based in part on use of the
11 + *    program.  An acceptable form of acknowledgement is citation of
12 + *    the article in which the program was described (Matthew
13 + *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 + *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 + *    Parallel Simulation Engine for Molecular Dynamics,"
16 + *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 + *
18 + * 2. Redistributions of source code must retain the above copyright
19 + *    notice, this list of conditions and the following disclaimer.
20 + *
21 + * 3. Redistributions in binary form must reproduce the above copyright
22 + *    notice, this list of conditions and the following disclaimer in the
23 + *    documentation and/or other materials provided with the
24 + *    distribution.
25 + *
26 + * This software is provided "AS IS," without a warranty of any
27 + * kind. All express or implied conditions, representations and
28 + * warranties, including any implied warranty of merchantability,
29 + * fitness for a particular purpose or non-infringement, are hereby
30 + * excluded.  The University of Notre Dame and its licensors shall not
31 + * be liable for any damages suffered by licensee as a result of
32 + * using, modifying or distributing the software or its
33 + * derivatives. In no event will the University of Notre Dame or its
34 + * licensors be liable for any lost revenue, profit or data, or for
35 + * direct, indirect, special, consequential, incidental or punitive
36 + * damages, however caused and regardless of the theory of liability,
37 + * arising out of the use of or inability to use software, even if the
38 + * University of Notre Dame has been advised of the possibility of
39 + * such damages.
40 + */
41 +
42   #include <cstdlib>
43   #include <cstdio>
44   #include <cstring>
# Line 7 | Line 48
48   #include <map>
49   #include <fstream>
50  
10 #include "io/Globals.hpp"
11 #include "brains/SimInfo.hpp"
12 #include "brains/SimSetup.hpp"
51   #include "applications/simpleBuilder/simpleBuilderCmd.h"
52 < #include "utils/StringUtils.hpp"
53 < #include "applications/simpleBuilder/LatticeFactory.hpp"
52 > #include "lattice/LatticeFactory.hpp"
53 > #include "utils/MoLocator.hpp"
54 > #include "lattice/Lattice.hpp"
55 > #include "brains/Register.hpp"
56 > #include "brains/SimInfo.hpp"
57 > #include "brains/SimCreator.hpp"
58 > #include "io/DumpWriter.hpp"
59   #include "math/Vector3.hpp"
60 < #include "applications/simpleBuilder/MoLocator.hpp"
61 < #include "applications/simpleBuilder/Lattice.hpp"
60 > #include "math/SquareMatrix3.hpp"
61 > #include "utils/StringUtils.hpp"
62  
63   using namespace std;
64 + using namespace oopse;
65  
66 < void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol);
67 < double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF);
66 > void createMdFile(const std::string&oldMdFileName,
67 >                  const std::string&newMdFileName,
68 >                  int nMol);
69  
70 < int main( int argc, char* argv[]){
70 > int main(int argc, char *argv []) {
71  
72 +  // register force fields
73 +  registerForceFields();
74 +  registerLattice();
75 +    
76    gengetopt_args_info args_info;
77 <  string latticeType;
78 <  string inputFileName;
79 <  string outPrefix;
80 <  string outMdFileName;
81 <  string outInitFileName;
82 <  SimInfo* oldInfo;
83 <  SimSetup* oldSimSetup;
84 <  BaseLattice* simpleLat;
36 <  int numMol;
37 <  double latticeConstant;
38 <  vector<double> lc;
39 <  double mass;
40 <  const double rhoConvertConst = 1.661;
41 <  double density;
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 <  double Hmat[3][3];
86 >  Mat3x3d hmat;
87    MoLocator *locator;
88 <  vector<Vector3d> latticePos;
89 <  vector<Vector3d> latticeOrt;
90 <  int numMolPerCell;
91 <  int curMolIndex;
92 <  DumpWriter* writer;
50 <  
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) ;
96 <  
94 >  if (cmdline_parser(argc, argv, &args_info) != 0)
95 >    exit(1);
96 >
97    density = args_info.density_arg;
98  
99    //get lattice type
100 <  latticeType = UpperCase(args_info.latticetype_arg);
101 <  if(!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)){
102 <    cerr << latticeType << " is an invalid lattice type" << endl;
103 <    cerr << LatticeFactory::getInstance()->toString() << endl;
104 <    exit(1);
100 >  latticeType = "FCC";
101 >
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
112 >  //get the number of unit cells in each direction:
113 >
114    nx = args_info.nx_arg;
115 <  if(nx <= 0){
116 <    cerr << "The number of unit cell in h direction must be greater than 0" << endl;
117 <    exit(1);
115 >
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    ny = args_info.ny_arg;
124 <  if(ny <= 0){
125 <    cerr << "The number of unit cell in l direction must be greater than 0" << endl;
126 <    exit(1);
124 >
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    nz = args_info.nz_arg;
133 <  if(nz <= 0){
134 <    cerr << "The number of unit cell in k direction must be greater than 0" << endl;
135 <    exit(1);
133 >
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 <        
140 >
141 >  int nSites = nMolPerCell * nx * ny * nz;
142 >
143    //get input file name
144 <  if (args_info.inputs_num)
144 >  if (args_info.inputs_num)
145      inputFileName = args_info.inputs[0];
146 <  else {                
147 <    cerr <<"You must specify a input file name.\n" << endl;
148 <    cmdline_parser_print_help();
149 <    exit(1);
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  
93  
153    //parse md file and set up the system
95  oldInfo = new SimInfo;
96  if(oldInfo == NULL){
97     cerr << "error in creating SimInfo" << endl;
98     exit(1);
99  }
154  
155 <  oldSimSetup = new SimSetup();  
156 <  if(oldSimSetup == NULL){
157 <     cerr << "error in creating SimSetup" << endl;
104 <     exit(1);
105 <  }
155 >  SimCreator oldCreator;
156 >  SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
157 >  Globals* simParams = oldInfo->getSimParams();
158  
159 <  oldSimSetup->suspendInit();
108 <  oldSimSetup->setSimInfo(oldInfo );
109 <  oldSimSetup->parseFile(&inputFileName[0] );
110 <  oldSimSetup->createSim();
111 <  
112 <  if(oldInfo->nComponents >=2){
113 <      cerr << "can not build the system with more than two components" << endl;
114 <      exit(1);
115 <  }
116 <  
117 <  //get mass of molecule.
118 <  //Due to the design of OOPSE, given atom type, we have to query forcefiled to get the mass
119 <  mass = getMolMass(oldInfo->compStamps[0], oldSimSetup->getForceField());
120 <  
121 <  //creat lattice
122 <        simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
123 <        if(simpleLat == NULL){
124 <                cerr << "Error in creating lattice" << endl;
125 <                exit(1);
126 <        }
159 >  // Calculate lattice constant (in Angstroms)
160  
161 <  numMolPerCell = simpleLat->getNumSitesPerCell();
161 >  RealType avgMass = getMolMass(oldInfo->getMoleculeStamp(0),
162 >                                  oldInfo->getForceField());
163 >
164 >  latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
165 >                        (RealType)(1.0 / 3.0));
166    
167 <  //calculate lattice constant (in Angstrom)
131 <  latticeConstant = pow(rhoConvertConst * numMolPerCell * mass /density, 1.0/3.0);
167 >  // Set the lattice constant
168    
133  //set lattice constant
169    lc.push_back(latticeConstant);
170    simpleLat->setLatticeConstant(lc);
136  
137  //calculate the total number of molecules
138  numMol = nx * ny * nz * numMolPerCell;
171  
172 <  if (oldInfo->n_mol != numMol){
172 >  // Calculate the lattice sites and fill the lattice vector.
173  
174 <    outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
143 <    outMdFileName = outPrefix + ".md";
174 >  // Get the standard orientations of the cell sites
175  
176 <    //creat new .md file on fly which corrects the number of molecule    
177 <    createMdFile(inputFileName, outMdFileName, numMol);
178 <    cerr << "SimpleBuilder Error: the number of molecule and the density are not matched" <<endl;
179 <    cerr << "A new .md file: " << outMdFileName << " is generated, use it to rerun the simpleBuilder" << endl;
180 <    exit(1);
176 >  latticeOrt = simpleLat->getLatticePointsOrt();
177 >
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 >        // 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 <  //determine the output file names  
153 <  if (args_info.output_given)
154 <    outInitFileName = args_info.output_arg;
155 <  else
156 <    outInitFileName = getPrefix(inputFileName.c_str())  + ".in";
197 >  outputFileName = args_info.output_arg;
198    
199 <  
159 <  //allocat memory for storing pos, vel and etc
160 <  oldInfo->getConfiguration()->createArrays(oldInfo->n_atoms);
161 <  for (int i = 0; i < oldInfo->n_atoms; i++)
162 <    oldInfo->atoms[i]->setCoords();  
199 >  // create a new .md file on the fly which corrects the number of molecules
200  
201 <  //creat Molocator
165 <  locator = new MoLocator(oldInfo->compStamps[0], oldSimSetup->getForceField());
201 >  createMdFile(inputFileName, outputFileName, nSites);
202  
203 <  //fill Hmat
204 <  Hmat[0][0] = nx * latticeConstant;
169 <  Hmat[0][1] = 0.0;
170 <  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;
174 <  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;
178 <  Hmat[2][2] = nz * latticeConstant ;
209 >  SimCreator newCreator;
210 >  SimInfo* newInfo = newCreator.createSim(outputFileName, false);
211  
212 <  //set Hmat
181 <  oldInfo->setBoxM(Hmat);
182 <  
183 <  //place the molecules
212 >  // fill Hmat
213  
214 <  curMolIndex = 0;
214 >  hmat(0, 0) = nx * latticeConstant;
215 >  hmat(0, 1) = 0.0;
216 >  hmat(0, 2) = 0.0;
217  
218 <  //get the orientation of the cell sites
219 <  //for the same type of molecule in same lattice, it will not change
220 <  latticeOrt = simpleLat->getLatticePointsOrt();
218 >  hmat(1, 0) = 0.0;
219 >  hmat(1, 1) = ny * latticeConstant;
220 >  hmat(1, 2) = 0.0;
221  
222 <  for(int i =0; i < nx; i++){
223 <    for(int j=0; j < ny; j++){
224 <       for(int k = 0; k < nz; k++){
222 >  hmat(2, 0) = 0.0;
223 >  hmat(2, 1) = 0.0;
224 >  hmat(2, 2) = nz * latticeConstant;
225  
226 <          //get the position of the cell sites
196 <          simpleLat->getLatticePointsPos(latticePos, i, j, k);
226 >  // Set Hmat
227  
228 <          for(int l = 0; l < numMolPerCell; l++)
199 <            locator->placeMol(latticePos[l], latticeOrt[l], &(oldInfo->molecules[curMolIndex++]));
200 <       }
201 <    }
202 <  }
228 >  newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
229  
230 <  //create dumpwriter and write out the coordinates
231 <  oldInfo->finalName = outInitFileName;
232 <  writer = new DumpWriter( oldInfo );
233 <  if(writer == NULL){
234 <    cerr << "error in creating DumpWriter" << endl;
235 <    exit(1);    
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 <  writer->writeFinal(0);
240 <  cout << "new initial configuration file: " << outInitFileName <<" is generated." << endl;
213 <  //delete objects
239 >  
240 >  // Create DumpWriter and write out the coordinates
241  
242 <  //delete oldInfo and oldSimSetup
216 <  if(oldInfo != NULL)
217 <     delete oldInfo;
242 >  writer = new DumpWriter(newInfo, outputFileName);
243    
244 <  if(oldSimSetup != NULL)
245 <     delete oldSimSetup;
246 <  
247 <  if (writer != NULL)
248 <    delete writer;
244 >  if (writer == NULL) {
245 >    sprintf(painCave.errMsg, "error in creating DumpWriter");
246 >    painCave.isFatal = 1;
247 >    simError();
248 >  }
249 >
250 >  writer->writeDump();
251 >
252 >  // deleting the writer will put the closing at the end of the dump file.
253 >
254 >  delete writer;
255 >
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 string& oldMdFileName, const string& newMdFileName, int numMol){
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;
# Line 235 | Line 273 | void createMdFile(const string& oldMdFileName, const s
273    newMdFile.open(newMdFileName.c_str());
274  
275    oldMdFile.getline(buffer, MAXLEN);
238  while(!oldMdFile.eof()){
276  
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 << endl;
283 <    }
284 <    else
246 <      newMdFile << buffer << endl;
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    }
288  
289    oldMdFile.close();
290    newMdFile.close();
253
291   }
292  
256 double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF){
257  int nAtoms;
258  AtomStamp* currAtomStamp;
259  double totMass;
260  
261  totMass = 0;
262  nAtoms = molStamp->getNAtoms();
263
264  for(size_t i=0; i<nAtoms; i++){
265    currAtomStamp = molStamp->getAtom(i);
266    totMass += myFF->getAtomTypeMass(currAtomStamp->getType());
267  }
268
269  return totMass;
270 }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines