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);
double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF);

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

  gengetopt_args_info args_info;
  string latticeType;
  string inputFileName;
  string outPrefix;
  string outMdFileName;
  SimInfo* oldInfo;
  SimSetup* oldSimSetup;
  BaseLattice* simpleLat;
  int numMol;
  double latticeConstant;
  vector<double> lc;
  double mass;
  const double rhoConvertConst = 1.661;
  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) ;

  //process density
  if(!args_info.density_given && !args_info.ndensity_given){
    cerr << "SimpleBuilder Error: density or number density must be given" << endl;
    cmdline_parser_print_help();
    exit(1);
  }
  else if(args_info.density_given){
    if( args_info.ndensity_given)
      cerr << "SimpleBuilder Warning: both of density and number density are given, use density" << endl;
    
    density = args_info.density_arg;
  }
  else if(args_info.ndensity_given){
    //convert number density to density 
  }


  //get lattice type
  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);
  }

  //get the number of unit cell
  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";

  //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->suspendInit();
  oldSimSetup->setSimInfo(oldInfo );
  oldSimSetup->parseFile(&inputFileName[0] );
  oldSimSetup->createSim(); 
  
  if(oldInfo->nComponents >=2){
      cerr << "can not build the system with more than two components" << endl;
      exit(1);
  }
  
  //get mass of molecule. 
  //Due to the design of OOPSE, given atom type, we have to query forcefiled to get the mass
  mass = getMolMass(oldInfo->compStamps[0], oldSimSetup->getForceField());
  
  //creat lattice
        simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
        if(simpleLat == NULL){
                cerr << "Error in creating lattice" << endl;
                exit(1);
        }

  numMolPerCell = simpleLat->getNumSitesPerCell();
  
  //calculate lattice constant (in Angstrom)
  latticeConstant = pow(rhoConvertConst * numMolPerCell * mass /density, 1.0/3.0);
  
  //set lattice constant
  lc.push_back(latticeConstant);
  simpleLat->setLatticeConstant(lc);
  
  //calculate the total number of molecules
  numMol = nx * ny * nz * numMolPerCell;

  if (oldInfo->n_mol != numMol){
    //creat new .md file on fly which corrects the number of molecule 
    createMdFile(inputFileName, outMdFileName, numMol);
    cerr << "SimpleBuilder Error: the number of molecule and the density are not matched" <<endl;
    cerr << "A new .md file: " << outMdFileName << " is generated, use it to rerun the simpleBuilder" << endl;
    exit(1);
  }

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

  //creat Molocator
  locator = new MoLocator(oldInfo->compStamps[0], oldSimSetup->getForceField());

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

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

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

  //set Hmat
  oldInfo->setBoxM(Hmat);
  
  //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], &(oldInfo->molecules[curMolIndex++]));
       }
    }
  }

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

  //delete oldInfo and oldSimSetup
  if(oldInfo != NULL)
     delete oldInfo;
  
  if(oldSimSetup != NULL)
     delete oldSimSetup; 
  
  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.c_str());

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

    //correct molecule number
    if(strstr(buffer, "nMol") !=NULL){      
      sprintf(buffer, "\t\tnMol = %d;", numMol);
      newMdFile << buffer << endl;
    }
    else
      newMdFile << buffer << endl;

    oldMdFile.getline(buffer, MAXLEN);
  }

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

}

double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF){
  int nAtoms;
  AtomStamp* currAtomStamp;
  double totMass;
  
  totMass = 0;
  nAtoms = molStamp->getNAtoms();

  for(size_t i=0; i<nAtoms; i++){
    currAtomStamp = molStamp->getAtom(i);
    totMass += myFF->getAtomTypeMass(currAtomStamp->getType());
  }

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