| 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 | 
+ | 
} |