utils/sysbuilder/simpleBuilder.cpp

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

#include "Globals.hpp"
#include "SimInfo.hpp"
#include "SimSetup.hpp"
#include "simpleBuilderCmd.h"
#include "StringUtils.hpp"
#include "LatticeFactory.hpp"
#include "Vector3d.hpp"
#include "MoLocator.hpp"
#include "Lattice.hpp"

using namespace std;

void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol);

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<double> 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()->toString() << endl;
    exit(1);
  }
        
  nx = args_info.nx_arg;
  if(nx <= 0){
    cerr << "The number of unit cell in h direction must be greater than 0" << endl;
    exit(1);
  }

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

  nz = args_info.nz_arg;
  if(nz <= 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];
  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();  
  if(oldSimSetup == NULL){
     cerr << "error in creating SimSetup" << endl;
     exit(1);
  }
    
  oldSimSetup->setSimInfo(oldInfo );
  oldSimSetup->parseFile(&inputFileName[0] );


  if(oldInfo->nComponents >=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->getNumSitesPerCell();
  //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();  
  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->getForceField());

  //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(0);

  
  //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];
  string line;

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

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

    if(line.find("nMol") < line.size()){
      
      sprintf(buffer, "\t\tnMol = %d;", numMol);
      newMdFile << buffer << endl;
    }
    else
      newMdFile << buffer << endl;

    oldMdFile.getline(buffer, MAXLEN);
  }

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

}
Revision:	1429
Committed:	Wed Jul 28 20:29:49 2004 UTC (19 years, 11 months ago) by tim
File size:	6340 byte(s)
Log Message:	simpleBuilder in progress
#	Content
1	#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 "Globals.hpp"
11	#include "SimInfo.hpp"
12	#include "SimSetup.hpp"
13	#include "simpleBuilderCmd.h"
14	#include "StringUtils.hpp"
15	#include "LatticeFactory.hpp"
16	#include "Vector3d.hpp"
17	#include "MoLocator.hpp"
18	#include "Lattice.hpp"
19
20	using namespace std;
21
22	void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol);
23
24	int main( int argc, char* argv[]){
25
26	gengetopt_args_info args_info;
27	string latticeType;
28	string inputFileName;
29	string outPrefix;
30	string outMdFileName;
31	string outInitFileName;
32	SimInfo* oldInfo;
33	SimInfo* newInfo;
34	SimSetup* oldSimSetup;
35	SimSetup* newSimSetup;
36	BaseLattice* simpleLat;
37	int numMol;
38	double latticeConstant;
39	vector<double> lc;
40	double mass;
41	double density;
42	int nx, ny, nz;
43	double Hmat[3][3];
44	MoLocator *locator;
45	vector<Vector3d> latticePos;
46	vector<Vector3d> latticeOrt;
47	int numMolPerCell;
48	int curMolIndex;
49	DumpWriter* writer;
50
51	// parse command line arguments
52
53	if (cmdline_parser (argc, argv, &args_info) != 0)
54	exit(1) ;
55
56	if(!args_info.density_given && !args_info.ndensity_given){
57	cerr << "density or number density must be given" << endl;
58	return false;
59	}
60
61	latticeType = UpperCase(args_info.latticetype_arg);
62	if(!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)){
63	cerr << latticeType << " is an invalid lattice type" << endl;
64	cerr << LatticeFactory::getInstance()->toString() << endl;
65	exit(1);
66	}
67
68	nx = args_info.nx_arg;
69	if(nx <= 0){
70	cerr << "The number of unit cell in h direction must be greater than 0" << endl;
71	exit(1);
72	}
73
74	ny = args_info.ny_arg;
75	if(ny <= 0){
76	cerr << "The number of unit cell in l direction must be greater than 0" << endl;
77	exit(1);
78	}
79
80	nz = args_info.nz_arg;
81	if(nz <= 0){
82	cerr << "The number of unit cell in k direction must be greater than 0" << endl;
83	exit(1);
84	}
85
86	//get input file name
87	if (args_info.inputs_num)
88	inputFileName = args_info.inputs[0];
89	else {
90	cerr <<"You must specify a input file name.\n" << endl;
91	cmdline_parser_print_help();
92	exit(1);
93	}
94
95	//determine the output file names
96
97	if (args_info.output_given){
98	outPrefix = args_info.output_arg;
99	}
100	else
101	outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
102
103	outMdFileName = outPrefix + ".md";
104	outInitFileName = outPrefix + ".in";
105
106	//parse md file and set up the system
107	oldInfo = new SimInfo;
108	if(oldInfo == NULL){
109	cerr << "error in creating SimInfo" << endl;
110	exit(1);
111	}
112
113	oldSimSetup = new SimSetup();
114	if(oldSimSetup == NULL){
115	cerr << "error in creating SimSetup" << endl;
116	exit(1);
117	}
118
119	oldSimSetup->setSimInfo(oldInfo );
120	oldSimSetup->parseFile(&inputFileName[0] );
121
122
123	if(oldInfo->nComponents >=2){
124	cerr << "can not build the system with more than two components" << endl;
125	exit(1);
126	}
127
128	//creat lattice
129	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
130	if(simpleLat == NULL){
131	cerr << "Error in creating lattice" << endl;
132	exit(1);
133	}
134
135
136	numMolPerCell = simpleLat->getNumSitesPerCell();
137	//calculate lattice constant
138	latticeConstant = pow(numMolPerCell * mass /density, 1.0/3.0);
139
140	//set lattice constant
141	lc.push_back(latticeConstant);
142	simpleLat->setLatticeConstant(lc);
143
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();
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	//set Hmat
181	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	locator = new MoLocator(newInfo->compStamps[0], newSimSetup->getForceField());
190
191	//place the molecules
192
193	curMolIndex = 0;
194
195	//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
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	//get the position of the cell sites
204	simpleLat->getLatticePointsPos(latticePos, i, j, k);
205
206	for(int l = 0; l < numMolPerCell; l++)
207	locator->placeMol(latticePos[l], latticeOrt[l], &(newInfo->molecules[curMolIndex++]));
208	}
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(0);
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	string line;
245
246	//create new .md file based on old .md file
247	oldMdFile.open(oldMdFileName.c_str());
248	newMdFile.open(newMdFileName.c_str());
249
250	oldMdFile.getline(buffer, MAXLEN);
251	while(!oldMdFile.eof()){
252
253	if(line.find("nMol") < line.size()){
254
255	sprintf(buffer, "\t\tnMol = %d;", numMol);
256	newMdFile << buffer << endl;
257	}
258	else
259	newMdFile << buffer << endl;
260
261	oldMdFile.getline(buffer, MAXLEN);
262	}
263
264	oldMdFile.close();
265	newMdFile.close();
266
267	}