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 |
|
|
39 |
|
double latticeConstant; |
40 |
|
vector<double> lc; |
41 |
|
double mass; |
42 |
+ |
const double rhoConvertConst = 1.661; |
43 |
|
double density; |
44 |
|
int nx, ny, nz; |
45 |
|
double Hmat[3][3]; |
49 |
|
int numMolPerCell; |
50 |
|
int curMolIndex; |
51 |
|
DumpWriter* writer; |
52 |
< |
|
52 |
> |
|
53 |
|
// parse command line arguments |
52 |
– |
|
54 |
|
if (cmdline_parser (argc, argv, &args_info) != 0) |
55 |
|
exit(1) ; |
56 |
|
|
57 |
+ |
//process density |
58 |
|
if(!args_info.density_given && !args_info.ndensity_given){ |
59 |
< |
cerr << "density or number density must be given" << endl; |
60 |
< |
return false; |
59 |
> |
cerr << "SimpleBuilder Error: density or number density must be given" << endl; |
60 |
> |
cmdline_parser_print_help(); |
61 |
> |
exit(1); |
62 |
|
} |
63 |
+ |
else if(args_info.density_given){ |
64 |
+ |
if( args_info.ndensity_given) |
65 |
+ |
cerr << "SimpleBuilder Warning: both of density and number density are given, use density" << endl; |
66 |
+ |
|
67 |
+ |
density = args_info.density_arg; |
68 |
+ |
} |
69 |
+ |
else if(args_info.ndensity_given){ |
70 |
+ |
//convert number density to density |
71 |
+ |
} |
72 |
|
|
73 |
+ |
|
74 |
+ |
//get lattice type |
75 |
|
latticeType = UpperCase(args_info.latticetype_arg); |
76 |
|
if(!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)){ |
77 |
|
cerr << latticeType << " is an invalid lattice type" << endl; |
78 |
|
cerr << LatticeFactory::getInstance()->toString() << endl; |
79 |
|
exit(1); |
80 |
|
} |
81 |
< |
|
81 |
> |
|
82 |
> |
//get the number of unit cell |
83 |
|
nx = args_info.nx_arg; |
84 |
|
if(nx <= 0){ |
85 |
|
cerr << "The number of unit cell in h direction must be greater than 0" << endl; |
133 |
|
|
134 |
|
oldSimSetup->setSimInfo(oldInfo ); |
135 |
|
oldSimSetup->parseFile(&inputFileName[0] ); |
136 |
< |
|
122 |
< |
|
136 |
> |
oldSimSetup->createSim(); |
137 |
|
if(oldInfo->nComponents >=2){ |
138 |
|
cerr << "can not build the system with more than two components" << endl; |
139 |
|
exit(1); |
140 |
|
} |
141 |
|
|
142 |
+ |
|
143 |
+ |
//get mass of molecule. |
144 |
+ |
//Due to the design of OOPSE, given atom type, we have to query forcefiled to get the mass |
145 |
+ |
mass = getMolMass(oldInfo->compStamps[0], oldSimSetup->getForceField()); |
146 |
+ |
|
147 |
|
//creat lattice |
148 |
|
simpleLat = LatticeFactory::getInstance()->createLattice(latticeType); |
149 |
|
if(simpleLat == NULL){ |
151 |
|
exit(1); |
152 |
|
} |
153 |
|
|
135 |
– |
|
154 |
|
numMolPerCell = simpleLat->getNumSitesPerCell(); |
137 |
– |
//calculate lattice constant |
138 |
– |
latticeConstant = pow(numMolPerCell * mass /density, 1.0/3.0); |
155 |
|
|
156 |
+ |
//calculate lattice constant (in Angstrom) |
157 |
+ |
latticeConstant = pow(rhoConvertConst * numMolPerCell * mass /density, 1.0/3.0); |
158 |
+ |
|
159 |
|
//set lattice constant |
160 |
|
lc.push_back(latticeConstant); |
161 |
|
simpleLat->setLatticeConstant(lc); |
162 |
|
|
163 |
< |
//calculate the total |
163 |
> |
//calculate the total number of molecules |
164 |
|
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; |
165 |
|
|
166 |
|
//creat new .md file on fly |
167 |
|
createMdFile(inputFileName, outMdFileName, numMol); |
172 |
|
cerr << "error in creating SimInfo" << endl; |
173 |
|
exit(1); |
174 |
|
} |
175 |
< |
|
175 |
> |
|
176 |
|
newSimSetup = new SimSetup(); |
177 |
|
if(newSimSetup == NULL){ |
178 |
|
cerr << "error in creating SimSetup" << endl; |
183 |
|
newSimSetup->parseFile(&outMdFileName[0] ); |
184 |
|
newSimSetup->createSim(); |
185 |
|
|
180 |
– |
//set Hmat |
181 |
– |
newInfo->setBoxM(Hmat); |
182 |
– |
|
186 |
|
//allocat memory for storing pos, vel and etc |
187 |
|
newInfo->getConfiguration()->createArrays(newInfo->n_atoms); |
188 |
|
for (int i = 0; i < newInfo->n_atoms; i++) |
191 |
|
//creat Molocator |
192 |
|
locator = new MoLocator(newInfo->compStamps[0], newSimSetup->getForceField()); |
193 |
|
|
194 |
+ |
//fill Hmat |
195 |
+ |
Hmat[0][0] = nx * latticeConstant; |
196 |
+ |
Hmat[0][1] = 0.0; |
197 |
+ |
Hmat[0][2] = 0.0; |
198 |
+ |
|
199 |
+ |
Hmat[1][0] = ny * latticeConstant; |
200 |
+ |
Hmat[1][1] = 0.0; |
201 |
+ |
Hmat[1][2] = 0.0; |
202 |
+ |
|
203 |
+ |
Hmat[2][0] = nz * latticeConstant; |
204 |
+ |
Hmat[2][1] = 0.0; |
205 |
+ |
Hmat[2][2] = 0.0; |
206 |
+ |
|
207 |
+ |
//set Hmat |
208 |
+ |
newInfo->setBoxM(Hmat); |
209 |
+ |
|
210 |
|
//place the molecules |
211 |
|
|
212 |
|
curMolIndex = 0; |
235 |
|
exit(1); |
236 |
|
} |
237 |
|
writer->writeFinal(0); |
219 |
– |
|
238 |
|
|
239 |
|
//delete objects |
240 |
|
if(!oldInfo) |
283 |
|
newMdFile.close(); |
284 |
|
|
285 |
|
} |
286 |
+ |
|
287 |
+ |
double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF){ |
288 |
+ |
int nAtoms; |
289 |
+ |
AtomStamp* currAtomStamp; |
290 |
+ |
double totMass; |
291 |
+ |
|
292 |
+ |
totMass = 0; |
293 |
+ |
nAtoms = molStamp->getNAtoms(); |
294 |
+ |
|
295 |
+ |
for(size_t i=0; i<nAtoms; i++){ |
296 |
+ |
currAtomStamp = molStamp->getAtom(i); |
297 |
+ |
totMass += myFF->getAtomTypeMass(currAtomStamp->getType()); |
298 |
+ |
} |
299 |
+ |
|
300 |
+ |
return totMass; |
301 |
+ |
} |