utils/sysbuilder/simpleBuilder.cpp

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <string>
#include <map>
#include <fstream>

#include "simError.h"
#include "Globals.hpp"
#include "SimInfo.hpp"
#include "SimSetup.hpp"
#include "simpleBuilderCmd.hpp"
#include "StringUtils.hpp"
#include "LatticeFactory.hpp"

using namespace std;


int main( int argc, char* argv[]){

  gengetopt_args_info args_info;
  string latticeType;
  string inputFileName;
  string outPrefix;
  string outMdFileName;
  string outInitFileName;
  SimInfo* oldInfo;
  SimInfo* newInfo;
  SimSetup* oldSimSetup;
  SimSetup* newSimSetup;
  BaseLattice* simpleLat;
  int numMol;
  double latticeConstant;
  double mass;
  double density;
  int nx, ny, nz;
  double Hmat[3][3];
  MoLocator *locator;
  double* posX;
  double* posY;
  double* posZ;
  int numMolPerCell;
  int curMolIndex;
  DumpWriter* writer;

  // parse command line arguments

  if (cmdline_parser (argc, argv, &args_info) != 0)
    exit(1) ;

  if(!args_info.density_given && !args_info.ndensity_given){
    cerr << "density or number density must be given" << endl;
    return false;
  }

  latticeType = UpperCase(args_info.latticetype_arg);
  if(!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)){
    cerr << latticeType << " is an invalid lattice type" << endl;
    cerr << LatticeFactory::getInstance()->oString << endl;
    exit(1);
  }
        
  nx = args_info.nx_args;
  if(nx <= 0){
    cerr << "The number of unit cell in h direction must be greater than 0" << endl;
    exit(1);
  }

  nx = args_info.nx_args;
  if(nx <= 0){
    cerr << "The number of unit cell in l direction must be greater than 0" << endl;
    exit(1);
  }

  nx = args_info.nx_args;
  if(nx <= 0){
    cerr << "The number of unit cell in k direction must be greater than 0" << endl;
    exit(1);
  }
        
  //get input file name
  if (args_info.inputs_num) {
    inputFileName = args_info.inputs[0];
    cerr << in_name << "\n";      }
  else {                 
    cerr <<"You must specify a input file name.\n" << endl;
    cmdline_parser_print_help();
    exit(1);
  }

  //determine the output file names
  
  if (args_info.output_given){
    outPrefix = args_info.output_arg;
  }
  else
    outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
  
  outMdFileName = outPrefix + ".md";
  outInitFileName = outPrefix + ".in";

  //parse md file and set up the system
  oldInfo = new SimInfo;
  if(oldInfo == NULL){
     cerr << "error in creating SimInfo" << endl;
     exit(1);
  }

  oldSimSetup = new SimSetup(oldInfo);  
  if(oldSimSetup == NULL){
     cerr << "error in creating SimSetup" << endl;
     exit(1);
  }
    
  oldSimSetup->setSimInfo(oldInfo );
  oldSimSetup->parseFile(&inputFileName[0] );
  oldSimSetup->createSim();


  if(oldInfo->nComponets >=2){
      cerr << "can not build the system with more than two components" << endl;
      exit(1);
  }

  //creat lattice
        simpleLat = LatticeFactory::getInstance()       ->createLattice(latticeType);
        if(simpleLat == NULL){
                cerr << "Error in creating lattice" << endl;
                exit(1);
        }


  numMolPerCell = simpleLat->getNpoints();
  //calculate lattice constant 
  latticeConstant = pow(numMolPerCell * mass /density, 1.0/3.0);
  
  //set lattice constant
  simpleLat->setLatticeConstant(latticeConstant);
  
  //calculate the total
  numMol = args_info.nx_arg * args_info.ny_arg * args_info.nz_arg * numMolPerCell;

  //fill Hmat
  Hmat[0][0] = nx * latticeConstant;
  Hmat[0][1] = 0.0;
  Hmat[0][2] = 0.0;

  Hmat[1][0] = ny * latticeConstant;
  Hmat[1][1] = 0.0;
  Hmat[1][2] = 0.0;

  Hmat[2][0] = nz * latticeConstant;
  Hmat[2][1] = 0.0;
  Hmat[2][2] = 0.0;
  
  //creat new .md file on fly 
  createMdFile(inputFileName, outMdFileName, numMol);

  //parse new .md file and set up the system
  newInfo = new SimInfo;
  if(newInfo == NULL){
     cerr << "error in creating SimInfo" << endl;
     exit(1);
  }

  newSimSetup = new SimSetup(newInfo);  
  if(newSimSetup == NULL){
     cerr << "error in creating SimSetup" << endl;
     exit(1);
  }
    
  newSimSetup->setSimInfo(newInfo );
  newSimSetup->parseFile(&outMdFileName[0] );
  newSimSetup->createSim(); 

  //set Hamt
  newInfo->setBoxM(Hmat);

  //allocat memory for storing pos, vel and etc
  newInfo.getConfiguration()->createArrays(newInfo.n_atoms);
  for (int i = 0; i < newInfo.n_atoms; i++)
    newInfo.atoms[i]->setCoords();  

  //creat Molocator
  locator = new MoLocator(newInfo->compStamps[0]);

  //allocate memory for posX, posY, posZ

  posX = new doule[numMolPerCell];
  if(posX == NULL){
    cerr << "memory allocation error" << endl;
    exit(1);
  }

  posY = new doule[numMolPerCell];
  if(posX == NULL){
    cerr << "memory allocation error" << endl;
    exit(1);
  }

  posZ = new doule[numMolPerCell];
  if(posX == NULL){
    cerr << "memory allocation error" << endl;
    exit(1);
  }

  //place the molecules

  curMolIndex = 0;
  for(int i =0; i < nx; i++){
    for(int j=0; j < ny; j++){
       for(int k = 0; k < nz; k++){
        
          simpleLat->getLatticePoints(&posX, &posY, &posZ, i, j, k);

          for(int l = 0; l < numMolPerCell; l++)
            locator->placeMol(posX[l], posY[l], posZ[l], newInfo->molecules[curMolIndex++]);
       }
    }
  }

  //create dumpwriter and write out the coordinates
  writer = new DumpWriter( newInfo );
  if(writer == NULL){
    cerr << "error in creating DumpWriter" << endl;
    exit(1);    
  }
  writer->writeFinal();

  
  //delete objects
  if(!oldInfo)
     delete oldInfo;
  
  if(!oldSimSetup)
     delete oldSimSetup;

  if(!newInfo)
     delete newInfo;
  
  if(!newSimSetup)
     delete newSimSetup;

  if (posX)
     delete[] posX;

  if (posY)
     delete[] posY;

  if (posZ)
   delete[] posZ;

  if (writer != NULL)
    delete writer;
  return 0;
}

void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol){
  ifstream oldMdFile;
  ofstream newMdFile;
  const int MAXLEN = 65535;
  char buffer[MAXLEN];


  //create new .md file based on old .md file
  oldMdFile.open(oldMdFileName.c_str())
  newMdFile.open(newMdFileName);

  oldMdFile.getline(buffer, MAXLEN);
  while(!oldMdFile.eof()){

    if(line.find("nMol") < line.size())
      sprintf(buffer, "nMol = %s;", numMol);
    else
      newMdFile << buffer << endl;

    oldMdFile.getline(buffer, MAXLEN);
  }

  oldMdFile.close();
  newMdFile.close();

}
Revision:	1424
Committed:	Wed Jul 28 04:59:35 2004 UTC (19 years, 11 months ago) by tim
File size:	6405 byte(s)
Log Message:	add simpleBuilder
#	User	Rev	Content
1	tim	1424	#include <cstdlib>
2			#include <cstdio>
3			#include <cstring>
4			#include <cmath>
5			#include <iostream>
6			#include <string>
7			#include <map>
8			#include <fstream>
9
10			#include "simError.h"
11			#include "Globals.hpp"
12			#include "SimInfo.hpp"
13			#include "SimSetup.hpp"
14			#include "simpleBuilderCmd.hpp"
15			#include "StringUtils.hpp"
16			#include "LatticeFactory.hpp"
17
18			using namespace std;
19
20
21			int main( int argc, char* argv[]){
22
23			gengetopt_args_info args_info;
24			string latticeType;
25			string inputFileName;
26			string outPrefix;
27			string outMdFileName;
28			string outInitFileName;
29			SimInfo* oldInfo;
30			SimInfo* newInfo;
31			SimSetup* oldSimSetup;
32			SimSetup* newSimSetup;
33			BaseLattice* simpleLat;
34			int numMol;
35			double latticeConstant;
36			double mass;
37			double density;
38			int nx, ny, nz;
39			double Hmat[3][3];
40			MoLocator *locator;
41			double* posX;
42			double* posY;
43			double* posZ;
44			int numMolPerCell;
45			int curMolIndex;
46			DumpWriter* writer;
47
48			// parse command line arguments
49
50			if (cmdline_parser (argc, argv, &args_info) != 0)
51			exit(1) ;
52
53			if(!args_info.density_given && !args_info.ndensity_given){
54			cerr << "density or number density must be given" << endl;
55			return false;
56			}
57
58			latticeType = UpperCase(args_info.latticetype_arg);
59			if(!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)){
60			cerr << latticeType << " is an invalid lattice type" << endl;
61			cerr << LatticeFactory::getInstance()->oString << endl;
62			exit(1);
63			}
64
65			nx = args_info.nx_args;
66			if(nx <= 0){
67			cerr << "The number of unit cell in h direction must be greater than 0" << endl;
68			exit(1);
69			}
70
71			nx = args_info.nx_args;
72			if(nx <= 0){
73			cerr << "The number of unit cell in l direction must be greater than 0" << endl;
74			exit(1);
75			}
76
77			nx = args_info.nx_args;
78			if(nx <= 0){
79			cerr << "The number of unit cell in k direction must be greater than 0" << endl;
80			exit(1);
81			}
82
83			//get input file name
84			if (args_info.inputs_num) {
85			inputFileName = args_info.inputs[0];
86			cerr << in_name << "\n"; }
87			else {
88			cerr <<"You must specify a input file name.\n" << endl;
89			cmdline_parser_print_help();
90			exit(1);
91			}
92
93			//determine the output file names
94
95			if (args_info.output_given){
96			outPrefix = args_info.output_arg;
97			}
98			else
99			outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
100
101			outMdFileName = outPrefix + ".md";
102			outInitFileName = outPrefix + ".in";
103
104			//parse md file and set up the system
105			oldInfo = new SimInfo;
106			if(oldInfo == NULL){
107			cerr << "error in creating SimInfo" << endl;
108			exit(1);
109			}
110
111			oldSimSetup = new SimSetup(oldInfo);
112			if(oldSimSetup == NULL){
113			cerr << "error in creating SimSetup" << endl;
114			exit(1);
115			}
116
117			oldSimSetup->setSimInfo(oldInfo );
118			oldSimSetup->parseFile(&inputFileName[0] );
119			oldSimSetup->createSim();
120
121
122			if(oldInfo->nComponets >=2){
123			cerr << "can not build the system with more than two components" << endl;
124			exit(1);
125			}
126
127			//creat lattice
128			simpleLat = LatticeFactory::getInstance() ->createLattice(latticeType);
129			if(simpleLat == NULL){
130			cerr << "Error in creating lattice" << endl;
131			exit(1);
132			}
133
134
135			numMolPerCell = simpleLat->getNpoints();
136			//calculate lattice constant
137			latticeConstant = pow(numMolPerCell * mass /density, 1.0/3.0);
138
139			//set lattice constant
140			simpleLat->setLatticeConstant(latticeConstant);
141
142			//calculate the total
143			numMol = args_info.nx_arg * args_info.ny_arg * args_info.nz_arg * numMolPerCell;
144
145			//fill Hmat
146			Hmat[0][0] = nx * latticeConstant;
147			Hmat[0][1] = 0.0;
148			Hmat[0][2] = 0.0;
149
150			Hmat[1][0] = ny * latticeConstant;
151			Hmat[1][1] = 0.0;
152			Hmat[1][2] = 0.0;
153
154			Hmat[2][0] = nz * latticeConstant;
155			Hmat[2][1] = 0.0;
156			Hmat[2][2] = 0.0;
157
158			//creat new .md file on fly
159			createMdFile(inputFileName, outMdFileName, numMol);
160
161			//parse new .md file and set up the system
162			newInfo = new SimInfo;
163			if(newInfo == NULL){
164			cerr << "error in creating SimInfo" << endl;
165			exit(1);
166			}
167
168			newSimSetup = new SimSetup(newInfo);
169			if(newSimSetup == NULL){
170			cerr << "error in creating SimSetup" << endl;
171			exit(1);
172			}
173
174			newSimSetup->setSimInfo(newInfo );
175			newSimSetup->parseFile(&outMdFileName[0] );
176			newSimSetup->createSim();
177
178			//set Hamt
179			newInfo->setBoxM(Hmat);
180
181			//allocat memory for storing pos, vel and etc
182			newInfo.getConfiguration()->createArrays(newInfo.n_atoms);
183			for (int i = 0; i < newInfo.n_atoms; i++)
184			newInfo.atoms[i]->setCoords();
185
186			//creat Molocator
187			locator = new MoLocator(newInfo->compStamps[0]);
188
189			//allocate memory for posX, posY, posZ
190
191			posX = new doule[numMolPerCell];
192			if(posX == NULL){
193			cerr << "memory allocation error" << endl;
194			exit(1);
195			}
196
197			posY = new doule[numMolPerCell];
198			if(posX == NULL){
199			cerr << "memory allocation error" << endl;
200			exit(1);
201			}
202
203			posZ = new doule[numMolPerCell];
204			if(posX == NULL){
205			cerr << "memory allocation error" << endl;
206			exit(1);
207			}
208
209			//place the molecules
210
211			curMolIndex = 0;
212			for(int i =0; i < nx; i++){
213			for(int j=0; j < ny; j++){
214			for(int k = 0; k < nz; k++){
215
216			simpleLat->getLatticePoints(&posX, &posY, &posZ, i, j, k);
217
218			for(int l = 0; l < numMolPerCell; l++)
219			locator->placeMol(posX[l], posY[l], posZ[l], newInfo->molecules[curMolIndex++]);
220			}
221			}
222			}
223
224			//create dumpwriter and write out the coordinates
225			writer = new DumpWriter( newInfo );
226			if(writer == NULL){
227			cerr << "error in creating DumpWriter" << endl;
228			exit(1);
229			}
230			writer->writeFinal();
231
232
233			//delete objects
234			if(!oldInfo)
235			delete oldInfo;
236
237			if(!oldSimSetup)
238			delete oldSimSetup;
239
240			if(!newInfo)
241			delete newInfo;
242
243			if(!newSimSetup)
244			delete newSimSetup;
245
246			if (posX)
247			delete[] posX;
248
249			if (posY)
250			delete[] posY;
251
252			if (posZ)
253			delete[] posZ;
254
255			if (writer != NULL)
256			delete writer;
257			return 0;
258			}
259
260			void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol){
261			ifstream oldMdFile;
262			ofstream newMdFile;
263			const int MAXLEN = 65535;
264			char buffer[MAXLEN];
265
266
267			//create new .md file based on old .md file
268			oldMdFile.open(oldMdFileName.c_str())
269			newMdFile.open(newMdFileName);
270
271			oldMdFile.getline(buffer, MAXLEN);
272			while(!oldMdFile.eof()){
273
274			if(line.find("nMol") < line.size())
275			sprintf(buffer, "nMol = %s;", numMol);
276			else
277			newMdFile << buffer << endl;
278
279			oldMdFile.getline(buffer, MAXLEN);
280			}
281
282			oldMdFile.close();
283			newMdFile.close();
284
285			}