applications/simpleBuilder/simpleBuilder.cpp

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

#include "io/Globals.hpp"
#include "brains/SimInfo.hpp"
#include "brains/SimSetup.hpp"
#include "applications/simpleBuilder/simpleBuilderCmd.h"
#include "utils/StringUtils.hpp"
#include "applications/simpleBuilder/LatticeFactory.hpp"
#include "math/Vector3.hpp"
#include "applications/simpleBuilder/MoLocator.hpp"
#include "applications/simpleBuilder/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;
  string outInitFileName;
  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) ;
  
  density = args_info.density_arg;

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

  
  //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){

    outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
    outMdFileName = outPrefix + ".md";

    //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);
  }
  
  //determine the output file names  
  if (args_info.output_given)
    outInitFileName = args_info.output_arg;
  else
    outInitFileName = getPrefix(inputFileName.c_str())  + ".in";
  
  
  //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
  oldInfo->finalName = outInitFileName;
  writer = new DumpWriter( oldInfo );
  if(writer == NULL){
    cerr << "error in creating DumpWriter" << endl;
    exit(1);    
  }
  writer->writeFinal(0);
  cout << "new initial configuration file: " << outInitFileName <<" is generated." << endl;
  //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:	1695
Committed:	Mon Nov 1 22:52:57 2004 UTC (21 years, 8 months ago) by tim
File size:	7398 byte(s)
Log Message:	Molecule, Atom, DirectionalAtom, RigidBody and StuntDouble classes get compiled
#	User	Rev	Content
1	tim	1695	#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 "io/Globals.hpp"
11			#include "brains/SimInfo.hpp"
12			#include "brains/SimSetup.hpp"
13			#include "applications/simpleBuilder/simpleBuilderCmd.h"
14			#include "utils/StringUtils.hpp"
15			#include "applications/simpleBuilder/LatticeFactory.hpp"
16			#include "math/Vector3.hpp"
17			#include "applications/simpleBuilder/MoLocator.hpp"
18			#include "applications/simpleBuilder/Lattice.hpp"
19
20			using namespace std;
21
22			void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol);
23			double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF);
24
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			string outInitFileName;
33			SimInfo* oldInfo;
34			SimSetup* oldSimSetup;
35			BaseLattice* simpleLat;
36			int numMol;
37			double latticeConstant;
38			vector<double> lc;
39			double mass;
40			const double rhoConvertConst = 1.661;
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			if (cmdline_parser (argc, argv, &args_info) != 0)
53			exit(1) ;
54
55			density = args_info.density_arg;
56
57			//get lattice type
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()->toString() << endl;
62			exit(1);
63			}
64
65			//get the number of unit cell
66			nx = args_info.nx_arg;
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			ny = args_info.ny_arg;
73			if(ny <= 0){
74			cerr << "The number of unit cell in l direction must be greater than 0" << endl;
75			exit(1);
76			}
77
78			nz = args_info.nz_arg;
79			if(nz <= 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			else {
88			cerr <<"You must specify a input file name.\n" << endl;
89			cmdline_parser_print_help();
90			exit(1);
91			}
92
93
94			//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			}
100
101			oldSimSetup = new SimSetup();
102			if(oldSimSetup == NULL){
103			cerr << "error in creating SimSetup" << endl;
104			exit(1);
105			}
106
107			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			}
127
128			numMolPerCell = simpleLat->getNumSitesPerCell();
129
130			//calculate lattice constant (in Angstrom)
131			latticeConstant = pow(rhoConvertConst * numMolPerCell * mass /density, 1.0/3.0);
132
133			//set lattice constant
134			lc.push_back(latticeConstant);
135			simpleLat->setLatticeConstant(lc);
136
137			//calculate the total number of molecules
138			numMol = nx * ny * nz * numMolPerCell;
139
140			if (oldInfo->n_mol != numMol){
141
142			outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
143			outMdFileName = outPrefix + ".md";
144
145			//creat new .md file on fly which corrects the number of molecule
146			createMdFile(inputFileName, outMdFileName, numMol);
147			cerr << "SimpleBuilder Error: the number of molecule and the density are not matched" <<endl;
148			cerr << "A new .md file: " << outMdFileName << " is generated, use it to rerun the simpleBuilder" << endl;
149			exit(1);
150			}
151
152			//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";
157
158
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();
163
164			//creat Molocator
165			locator = new MoLocator(oldInfo->compStamps[0], oldSimSetup->getForceField());
166
167			//fill Hmat
168			Hmat[0][0] = nx * latticeConstant;
169			Hmat[0][1] = 0.0;
170			Hmat[0][2] = 0.0;
171
172			Hmat[1][0] = 0.0;
173			Hmat[1][1] = ny * latticeConstant;
174			Hmat[1][2] = 0.0;
175
176			Hmat[2][0] = 0.0;
177			Hmat[2][1] = 0.0;
178			Hmat[2][2] = nz * latticeConstant ;
179
180			//set Hmat
181			oldInfo->setBoxM(Hmat);
182
183			//place the molecules
184
185			curMolIndex = 0;
186
187			//get the orientation of the cell sites
188			//for the same type of molecule in same lattice, it will not change
189			latticeOrt = simpleLat->getLatticePointsOrt();
190
191			for(int i =0; i < nx; i++){
192			for(int j=0; j < ny; j++){
193			for(int k = 0; k < nz; k++){
194
195			//get the position of the cell sites
196			simpleLat->getLatticePointsPos(latticePos, i, j, k);
197
198			for(int l = 0; l < numMolPerCell; l++)
199			locator->placeMol(latticePos[l], latticeOrt[l], &(oldInfo->molecules[curMolIndex++]));
200			}
201			}
202			}
203
204			//create dumpwriter and write out the coordinates
205			oldInfo->finalName = outInitFileName;
206			writer = new DumpWriter( oldInfo );
207			if(writer == NULL){
208			cerr << "error in creating DumpWriter" << endl;
209			exit(1);
210			}
211			writer->writeFinal(0);
212			cout << "new initial configuration file: " << outInitFileName <<" is generated." << endl;
213			//delete objects
214
215			//delete oldInfo and oldSimSetup
216			if(oldInfo != NULL)
217			delete oldInfo;
218
219			if(oldSimSetup != NULL)
220			delete oldSimSetup;
221
222			if (writer != NULL)
223			delete writer;
224			return 0;
225			}
226
227			void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol){
228			ifstream oldMdFile;
229			ofstream newMdFile;
230			const int MAXLEN = 65535;
231			char buffer[MAXLEN];
232
233			//create new .md file based on old .md file
234			oldMdFile.open(oldMdFileName.c_str());
235			newMdFile.open(newMdFileName.c_str());
236
237			oldMdFile.getline(buffer, MAXLEN);
238			while(!oldMdFile.eof()){
239
240			//correct molecule number
241			if(strstr(buffer, "nMol") !=NULL){
242			sprintf(buffer, "\t\tnMol = %d;", numMol);
243			newMdFile << buffer << endl;
244			}
245			else
246			newMdFile << buffer << endl;
247
248			oldMdFile.getline(buffer, MAXLEN);
249			}
250
251			oldMdFile.close();
252			newMdFile.close();
253
254			}
255
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			}