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"
#include "Vector3d.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;
  vector lc;
  double mass;
  double density;
  int nx, ny, nz;
  double Hmat[3][3];
  MoLocator *locator;
  vector<Vector3d> latticePos;
  vector<Vector3d> latticeOrt;
  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
  lc.push_back(latticeConstant);
  simpleLat->setLatticeConstant(lc);
  
  //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 Hmat
  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], newSimSetup->the_ff);

  //place the molecules

  curMolIndex = 0;

  //get the orientation of the cell sites
  //for the same type of molecule in same lattice, it will not change
  latticeOrt = simpleLat->getLatticePointsOrt();

  for(int i =0; i < nx; i++){
    for(int j=0; j < ny; j++){
       for(int k = 0; k < nz; k++){

          //get the position of the cell sites
          simpleLat->getLatticePointsPos(latticePos, i, j, k);

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