applications/nanoparticleBuilder/nanorodBuilder.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 *  Created by Kelsey M. Stocker on 2/9/12.
 *  @author  Kelsey M. Stocker 
 *
 */

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

#include "config.h"
#include "shapedLatticeRod.hpp"
#include "nanorodBuilderCmd.h"
#include "lattice/LatticeFactory.hpp"
#include "utils/MoLocator.hpp"
#include "lattice/Lattice.hpp"
#include "brains/Register.hpp"
#include "brains/SimInfo.hpp"
#include "brains/SimCreator.hpp"
#include "io/DumpWriter.hpp"
#include "math/Vector3.hpp"
#include "math/SquareMatrix3.hpp"
#include "utils/StringUtils.hpp"

using namespace std;
using namespace OpenMD;
void createMdFile(const std::string&oldMdFileName, 
                  const std::string&newMdFileName,
                  std::vector<int> numMol);

int main(int argc, char *argv []) {
  
  registerLattice();
  
  gengetopt_args_info args_info;
  std::string latticeType;
  std::string inputFileName;
  std::string outputFileName;

  MoLocator* locator;
  int nComponents;
  double latticeConstant;
  std::vector<double> lc;

  RealType rodRadius;
  RealType rodLength;

  Mat3x3d hmat;
  std::vector<Vector3d> latticePos;
  std::vector<Vector3d> latticeOrt;
 
  DumpWriter *writer;
  
  // Parse Command Line Arguments
  if (cmdline_parser(argc, argv, &args_info) != 0)
    exit(1);
         
  /* get lattice type */
  latticeType = "FCC";

  /* get input file name */
  if (args_info.inputs_num)
    inputFileName = args_info.inputs[0];
  else {
    sprintf(painCave.errMsg, "No input .md file name was specified "
            "on the command line");
    painCave.isFatal = 1;
    cmdline_parser_print_help();
    simError();
  }
  
  /* parse md file and set up the system */
  SimCreator oldCreator;
  SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
  
  latticeConstant = args_info.latticeConstant_arg;
  rodRadius = args_info.radius_arg;
  rodLength = args_info.length_arg;
  Globals* simParams = oldInfo->getSimParams();
  
  /* Create nanorod */
  shapedLatticeRod nanoRod(latticeConstant, latticeType,
                           rodRadius, rodLength);
  
  /* Build a lattice and get lattice points for this lattice constant */
  vector<Vector3d> sites = nanoRod.getSites();
  vector<Vector3d> orientations = nanoRod.getOrientations();
  std::vector<int> vacancyTargets;
  vector<bool> isVacancy;
  
  Vector3d myLoc;
  RealType myR;
 
  for (int i = 0; i < sites.size(); i++) 
    isVacancy.push_back(false);

  // cerr << "checking vacancyPercent" << "\n";
  if (args_info.vacancyPercent_given) {
    // cerr << "vacancyPercent given" << "\n";
    if (args_info.vacancyPercent_arg < 0.0 || args_info.vacancyPercent_arg > 100.0) {
      sprintf(painCave.errMsg, "vacancyPercent was set to a non-sensical value.");
      painCave.isFatal = 1;
      simError();
    } else {
      RealType vF = args_info.vacancyPercent_arg / 100.0;
      //  cerr << "vacancyPercent = " << vF << "\n";
      RealType vIR;
      RealType vOR;
      if (args_info.vacancyInnerRadius_given) {
        vIR = args_info.vacancyInnerRadius_arg;
      } else {
        vIR = 0.0;
      }
      if (args_info.vacancyOuterRadius_given) {
        vOR = args_info.vacancyOuterRadius_arg;
      } else {
        vOR = rodRadius;
      }
      if (vIR >= 0.0 && vOR <= rodRadius && vOR >= vIR) {
        
        for (int i = 0; i < sites.size(); i++) {
          myLoc = sites[i];
          myR = myLoc.length();
          if (myR >= vIR && myR <= vOR) {
            vacancyTargets.push_back(i);
          }          
        }
        std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
        
        int nTargets = vacancyTargets.size();
        vacancyTargets.resize((int)(vF * nTargets));
        
                  
        sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
                "\tsites between %lf and %lf.", (int) vacancyTargets.size(), 
                vIR, vOR); 
        painCave.isFatal = 0;
        simError();

        isVacancy.clear();
        for (int i = 0; i < sites.size(); i++) {
          bool vac = false;
          for (int j = 0; j < vacancyTargets.size(); j++) {
            if (i == vacancyTargets[j]) vac = true;
          }
          isVacancy.push_back(vac);
        }
               
      } else {
        sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
                "\tinner or outer radii.  Check their values.");
        painCave.isFatal = 1;
        simError();
      }
    }
  }

  /* Get number of lattice sites */
  int nSites = sites.size() - vacancyTargets.size();

  // cerr << "sites.size() = " << sites.size() << "\n";
  // cerr << "nSites = " << nSites << "\n";
  // cerr << "vacancyTargets = " << vacancyTargets.size() << "\n";

  std::vector<Component*> components = simParams->getComponents();
  std::vector<RealType> molFractions;
  std::vector<RealType> shellRadii;
  std::vector<RealType> molecularMasses;
  std::vector<int> nMol;
  std::map<int, int> componentFromSite;
  nComponents = components.size();
  // cerr << "nComponents = " << nComponents << "\n";

  if (args_info.molFraction_given && args_info.shellRadius_given) {
    sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
            "arguments, but not both!");
    painCave.isFatal = 1;
    simError();
  }
  
  if (nComponents == 1) {
    molFractions.push_back(1.0);    
    shellRadii.push_back(rodRadius);
  } else if (args_info.molFraction_given) {
    if ((int)args_info.molFraction_given == nComponents) {
      for (int i = 0; i < nComponents; i++) {
        molFractions.push_back(args_info.molFraction_arg[i]);
      }
    } else if ((int)args_info.molFraction_given == nComponents-1) {
      RealType remainingFraction = 1.0;
      for (int i = 0; i < nComponents-1; i++) {
        molFractions.push_back(args_info.molFraction_arg[i]);
        remainingFraction -= molFractions[i];
      }
      molFractions.push_back(remainingFraction);
    } else {    
      sprintf(painCave.errMsg, "nanorodBuilder can't figure out molFractions "
              "for all of the components in the <MetaData> block.");
      painCave.isFatal = 1;
      simError();
    }
  } else if ((int)args_info.shellRadius_given) {
    if ((int)args_info.shellRadius_given == nComponents) {
      for (int i = 0; i < nComponents; i++) {
        shellRadii.push_back(args_info.shellRadius_arg[i]);
      }
    } else if ((int)args_info.shellRadius_given == nComponents-1) {
      for (int i = 0; i < nComponents-1; i++) {
        shellRadii.push_back(args_info.shellRadius_arg[i]);
      }
      shellRadii.push_back(rodRadius);
    } else {    
      sprintf(painCave.errMsg, "nanorodBuilder can't figure out the\n"
              "\tshell radii for all of the components in the <MetaData> block.");
      painCave.isFatal = 1;
      simError();
    }
  } else {
    sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
            "\tbut have not specified either molFraction or shellRadius arguments.");
    painCave.isFatal = 1;
    simError();
  }
  
  if (args_info.molFraction_given) {
    RealType totalFraction = 0.0;
    
    /* Do some simple sanity checking*/
    
    for (int i = 0; i < nComponents; i++) {
      if (molFractions.at(i) < 0.0) {
        sprintf(painCave.errMsg, "One of the requested molFractions was"
                " less than zero!");
        painCave.isFatal = 1;
        simError();
      }
      if (molFractions.at(i) > 1.0) {
        sprintf(painCave.errMsg, "One of the requested molFractions was"
                " greater than one!");
        painCave.isFatal = 1;
        simError();
      }
      totalFraction += molFractions.at(i);
    }
    if (abs(totalFraction - 1.0) > 1e-6) {
      sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
      painCave.isFatal = 1;
      simError();
    }
    
    int remaining = nSites;
    for (int i=0; i < nComponents-1; i++) {    
      nMol.push_back(int((RealType)nSites * molFractions.at(i)));
      remaining -= nMol.at(i);
    }
    nMol.push_back(remaining);
    
    // recompute actual mol fractions and perform final sanity check:
    
    int totalMolecules = 0;
    for (int i=0; i < nComponents; i++) {
      molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
      totalMolecules += nMol.at(i);
    }
    if (totalMolecules != nSites) {
      sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
              "to the number of lattice sites!");
      painCave.isFatal = 1;
      simError();
    }
  } else {

    for (int i = 0; i < shellRadii.size(); i++) {
      if (shellRadii.at(i) > rodRadius + 1e-6 ) {
        sprintf(painCave.errMsg, "One of the shellRadius values exceeds the rod Radius.");
        painCave.isFatal = 1;
        simError();
      } 
      if (shellRadii.at(i) <= 0.0 ) {
        sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!");
        painCave.isFatal = 1;
        simError();
      }
    }
  }

  vector<int> ids;           
  if ((int)args_info.molFraction_given){
    //  cerr << "molFraction given 2" << "\n";
    sprintf(painCave.errMsg, "Creating a randomized spherically-capped nanorod.");
    painCave.isFatal = 0;
    simError();
    /* Random rod is the default case*/

    for (int i = 0; i < sites.size(); i++) 
      if (!isVacancy[i]) ids.push_back(i);
    
    std::random_shuffle(ids.begin(), ids.end());
    
  } else{ 
    sprintf(painCave.errMsg, "Creating an fcc nanorod.");
    painCave.isFatal = 0;
    simError();

    RealType smallestSoFar;
    int myComponent = -1;
    nMol.clear();
    nMol.resize(nComponents);

    // cerr << "shellRadii[0] " << shellRadii[0] << "\n";
    //  cerr << "rodRadius " << rodRadius << "\n";

    for (int i = 0; i < sites.size(); i++) {
      myLoc = sites[i];
      myR = myLoc.length();
      smallestSoFar = rodRadius;  
      //cerr << "vac = " << isVacancy[i]<< "\n";
    
      if (!isVacancy[i]) {


        // for (int j = 0; j < nComponents; j++) {
        //   if (myR <= shellRadii[j]) {
        //     if (shellRadii[j] <= smallestSoFar) {
        //       smallestSoFar = shellRadii[j];
        //       myComponent = j;
        //     }
        //   }
        // }
        myComponent = 0;
        componentFromSite[i] = myComponent;
        nMol[myComponent]++;
        //      cerr << "nMol for myComp(" << myComponent<<") = " << nMol[myComponent] << "\n";
      }
    }       
  }
  //     cerr << "nMol = " << nMol.at(0) << "\n";

  outputFileName = args_info.output_arg;
   
  //creat new .md file on fly which corrects the number of molecule    

  createMdFile(inputFileName, outputFileName, nMol);
  
  if (oldInfo != NULL)
    delete oldInfo;
  
  SimCreator newCreator;
  SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
    
  // Place molecules
  Molecule* mol;
  SimInfo::MoleculeIterator mi;
  mol = NewInfo->beginMolecule(mi);

  int l = 0;
  int whichSite = 0;

  for (int i = 0; i < nComponents; i++){
    locator = new MoLocator(NewInfo->getMoleculeStamp(i), 
                            NewInfo->getForceField());
    
    //   cerr << "nMol = " << nMol.at(i) << "\n";
    if (!args_info.molFraction_given) {
      for (int n = 0; n < sites.size(); n++) {
        if (!isVacancy[n]) {
          if (componentFromSite[n] == i) {
            mol = NewInfo->getMoleculeByGlobalIndex(l);
            locator->placeMol(sites[n], orientations[n], mol);
            l++;
          }
        }
      }
    } else {
      for (int n = 0; n < nMol.at(i); n++) {
        mol = NewInfo->getMoleculeByGlobalIndex(l);
        locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
        l++;
      }
    }
  }
  
  //fill Hmat
  hmat(0, 0)=  10.0*rodRadius;
  hmat(0, 1) = 0.0;
  hmat(0, 2) = 0.0;
  
  hmat(1, 0) = 0.0;
  hmat(1, 1) = 10.0*rodRadius;
  hmat(1, 2) = 0.0;
  
  hmat(2, 0) = 0.0;
  hmat(2, 1) = 0.0;
  hmat(2, 2) = 5.0*rodLength + 2.0*rodRadius;
  
  //set Hmat
  NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
  
  
  //create dumpwriter and write out the coordinates
  writer = new DumpWriter(NewInfo, outputFileName);
  
  if (writer == NULL) {
    sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
    painCave.isFatal = 1;
    simError();
  }
  
  writer->writeDump();

  // deleting the writer will put the closing at the end of the dump file

  delete writer;

  // cleanup a by calling sim error.....
  sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
          "generated.\n", outputFileName.c_str());
  painCave.isFatal = 0;
  simError();
  return 0;
}

void createMdFile(const std::string&oldMdFileName, 
                  const std::string&newMdFileName,
                  std::vector<int> nMol) {
  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);

  int i = 0;
  while (!oldMdFile.eof()) {

    //correct molecule number
    if (strstr(buffer, "nMol") != NULL) {
      if(i<nMol.size()){
        sprintf(buffer, "\tnMol = %i;", nMol.at(i));
        newMdFile << buffer << std::endl;
        i++;
      }
    } else
      newMdFile << buffer << std::endl;
    
    oldMdFile.getline(buffer, MAXLEN);
  }
  
  oldMdFile.close();
  newMdFile.close();

  if (i != nMol.size()) {
    sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
            "\tstatements in component blocks.  Make sure that all\n"
            "\tcomponents in the template file have nMol=1");
    painCave.isFatal = 1;
    simError();
  }
    
}

Revision:	1725
Committed:	Sat May 26 18:13:43 2012 UTC (12 years, 11 months ago) by gezelter
File size:	15974 byte(s)
Log Message:	Individual ForceField classes have been removed (they were essentially all duplicates anyway). ForceField has moved to brains, and since only one force field is in play at any time, the ForceFieldFactory and Register methods have been removed.
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
12	* 2. Redistributions in binary form must reproduce the above copyright
13	* notice, this list of conditions and the following disclaimer in the
14	* documentation and/or other materials provided with the
15	* distribution.
16	*
17	* This software is provided "AS IS," without a warranty of any
18	* kind. All express or implied conditions, representations and
19	* warranties, including any implied warranty of merchantability,
20	* fitness for a particular purpose or non-infringement, are hereby
21	* excluded. The University of Notre Dame and its licensors shall not
22	* be liable for any damages suffered by licensee as a result of
23	* using, modifying or distributing the software or its
24	* derivatives. In no event will the University of Notre Dame or its
25	* licensors be liable for any lost revenue, profit or data, or for
26	* direct, indirect, special, consequential, incidental or punitive
27	* damages, however caused and regardless of the theory of liability,
28	* arising out of the use of or inability to use software, even if the
29	* University of Notre Dame has been advised of the possibility of
30	* such damages.
31	*
32	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33	* research, please cite the appropriate papers when you publish your
34	* work. Good starting points are:
35	*
36	* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	* Created by Kelsey M. Stocker on 2/9/12.
42	* @author Kelsey M. Stocker
43	*
44	*/
45
46	#include <cstdlib>
47	#include <cstdio>
48	#include <cstring>
49	#include <cmath>
50	#include <iostream>
51	#include <string>
52	#include <map>
53	#include <fstream>
54	#include <algorithm>
55
56	#include "config.h"
57	#include "shapedLatticeRod.hpp"
58	#include "nanorodBuilderCmd.h"
59	#include "lattice/LatticeFactory.hpp"
60	#include "utils/MoLocator.hpp"
61	#include "lattice/Lattice.hpp"
62	#include "brains/Register.hpp"
63	#include "brains/SimInfo.hpp"
64	#include "brains/SimCreator.hpp"
65	#include "io/DumpWriter.hpp"
66	#include "math/Vector3.hpp"
67	#include "math/SquareMatrix3.hpp"
68	#include "utils/StringUtils.hpp"
69
70	using namespace std;
71	using namespace OpenMD;
72	void createMdFile(const std::string&oldMdFileName,
73	const std::string&newMdFileName,
74	std::vector<int> numMol);
75
76	int main(int argc, char *argv []) {
77
78	registerLattice();
79
80	gengetopt_args_info args_info;
81	std::string latticeType;
82	std::string inputFileName;
83	std::string outputFileName;
84
85	MoLocator* locator;
86	int nComponents;
87	double latticeConstant;
88	std::vector<double> lc;
89
90	RealType rodRadius;
91	RealType rodLength;
92
93	Mat3x3d hmat;
94	std::vector<Vector3d> latticePos;
95	std::vector<Vector3d> latticeOrt;
96
97	DumpWriter *writer;
98
99	// Parse Command Line Arguments
100	if (cmdline_parser(argc, argv, &args_info) != 0)
101	exit(1);
102
103	/* get lattice type */
104	latticeType = "FCC";
105
106	/* get input file name */
107	if (args_info.inputs_num)
108	inputFileName = args_info.inputs[0];
109	else {
110	sprintf(painCave.errMsg, "No input .md file name was specified "
111	"on the command line");
112	painCave.isFatal = 1;
113	cmdline_parser_print_help();
114	simError();
115	}
116
117	/* parse md file and set up the system */
118	SimCreator oldCreator;
119	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
120
121	latticeConstant = args_info.latticeConstant_arg;
122	rodRadius = args_info.radius_arg;
123	rodLength = args_info.length_arg;
124	Globals* simParams = oldInfo->getSimParams();
125
126	/* Create nanorod */
127	shapedLatticeRod nanoRod(latticeConstant, latticeType,
128	rodRadius, rodLength);
129
130	/* Build a lattice and get lattice points for this lattice constant */
131	vector<Vector3d> sites = nanoRod.getSites();
132	vector<Vector3d> orientations = nanoRod.getOrientations();
133	std::vector<int> vacancyTargets;
134	vector<bool> isVacancy;
135
136	Vector3d myLoc;
137	RealType myR;
138
139	for (int i = 0; i < sites.size(); i++)
140	isVacancy.push_back(false);
141
142	// cerr << "checking vacancyPercent" << "\n";
143	if (args_info.vacancyPercent_given) {
144	// cerr << "vacancyPercent given" << "\n";
145	if (args_info.vacancyPercent_arg < 0.0 \|\| args_info.vacancyPercent_arg > 100.0) {
146	sprintf(painCave.errMsg, "vacancyPercent was set to a non-sensical value.");
147	painCave.isFatal = 1;
148	simError();
149	} else {
150	RealType vF = args_info.vacancyPercent_arg / 100.0;
151	// cerr << "vacancyPercent = " << vF << "\n";
152	RealType vIR;
153	RealType vOR;
154	if (args_info.vacancyInnerRadius_given) {
155	vIR = args_info.vacancyInnerRadius_arg;
156	} else {
157	vIR = 0.0;
158	}
159	if (args_info.vacancyOuterRadius_given) {
160	vOR = args_info.vacancyOuterRadius_arg;
161	} else {
162	vOR = rodRadius;
163	}
164	if (vIR >= 0.0 && vOR <= rodRadius && vOR >= vIR) {
165
166	for (int i = 0; i < sites.size(); i++) {
167	myLoc = sites[i];
168	myR = myLoc.length();
169	if (myR >= vIR && myR <= vOR) {
170	vacancyTargets.push_back(i);
171	}
172	}
173	std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
174
175	int nTargets = vacancyTargets.size();
176	vacancyTargets.resize((int)(vF * nTargets));
177
178
179	sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
180	"\tsites between %lf and %lf.", (int) vacancyTargets.size(),
181	vIR, vOR);
182	painCave.isFatal = 0;
183	simError();
184
185	isVacancy.clear();
186	for (int i = 0; i < sites.size(); i++) {
187	bool vac = false;
188	for (int j = 0; j < vacancyTargets.size(); j++) {
189	if (i == vacancyTargets[j]) vac = true;
190	}
191	isVacancy.push_back(vac);
192	}
193
194	} else {
195	sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
196	"\tinner or outer radii. Check their values.");
197	painCave.isFatal = 1;
198	simError();
199	}
200	}
201	}
202
203	/* Get number of lattice sites */
204	int nSites = sites.size() - vacancyTargets.size();
205
206	// cerr << "sites.size() = " << sites.size() << "\n";
207	// cerr << "nSites = " << nSites << "\n";
208	// cerr << "vacancyTargets = " << vacancyTargets.size() << "\n";
209
210	std::vector<Component*> components = simParams->getComponents();
211	std::vector<RealType> molFractions;
212	std::vector<RealType> shellRadii;
213	std::vector<RealType> molecularMasses;
214	std::vector<int> nMol;
215	std::map<int, int> componentFromSite;
216	nComponents = components.size();
217	// cerr << "nComponents = " << nComponents << "\n";
218
219	if (args_info.molFraction_given && args_info.shellRadius_given) {
220	sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
221	"arguments, but not both!");
222	painCave.isFatal = 1;
223	simError();
224	}
225
226	if (nComponents == 1) {
227	molFractions.push_back(1.0);
228	shellRadii.push_back(rodRadius);
229	} else if (args_info.molFraction_given) {
230	if ((int)args_info.molFraction_given == nComponents) {
231	for (int i = 0; i < nComponents; i++) {
232	molFractions.push_back(args_info.molFraction_arg[i]);
233	}
234	} else if ((int)args_info.molFraction_given == nComponents-1) {
235	RealType remainingFraction = 1.0;
236	for (int i = 0; i < nComponents-1; i++) {
237	molFractions.push_back(args_info.molFraction_arg[i]);
238	remainingFraction -= molFractions[i];
239	}
240	molFractions.push_back(remainingFraction);
241	} else {
242	sprintf(painCave.errMsg, "nanorodBuilder can't figure out molFractions "
243	"for all of the components in the <MetaData> block.");
244	painCave.isFatal = 1;
245	simError();
246	}
247	} else if ((int)args_info.shellRadius_given) {
248	if ((int)args_info.shellRadius_given == nComponents) {
249	for (int i = 0; i < nComponents; i++) {
250	shellRadii.push_back(args_info.shellRadius_arg[i]);
251	}
252	} else if ((int)args_info.shellRadius_given == nComponents-1) {
253	for (int i = 0; i < nComponents-1; i++) {
254	shellRadii.push_back(args_info.shellRadius_arg[i]);
255	}
256	shellRadii.push_back(rodRadius);
257	} else {
258	sprintf(painCave.errMsg, "nanorodBuilder can't figure out the\n"
259	"\tshell radii for all of the components in the <MetaData> block.");
260	painCave.isFatal = 1;
261	simError();
262	}
263	} else {
264	sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
265	"\tbut have not specified either molFraction or shellRadius arguments.");
266	painCave.isFatal = 1;
267	simError();
268	}
269
270	if (args_info.molFraction_given) {
271	RealType totalFraction = 0.0;
272
273	/* Do some simple sanity checking*/
274
275	for (int i = 0; i < nComponents; i++) {
276	if (molFractions.at(i) < 0.0) {
277	sprintf(painCave.errMsg, "One of the requested molFractions was"
278	" less than zero!");
279	painCave.isFatal = 1;
280	simError();
281	}
282	if (molFractions.at(i) > 1.0) {
283	sprintf(painCave.errMsg, "One of the requested molFractions was"
284	" greater than one!");
285	painCave.isFatal = 1;
286	simError();
287	}
288	totalFraction += molFractions.at(i);
289	}
290	if (abs(totalFraction - 1.0) > 1e-6) {
291	sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
292	painCave.isFatal = 1;
293	simError();
294	}
295
296	int remaining = nSites;
297	for (int i=0; i < nComponents-1; i++) {
298	nMol.push_back(int((RealType)nSites * molFractions.at(i)));
299	remaining -= nMol.at(i);
300	}
301	nMol.push_back(remaining);
302
303	// recompute actual mol fractions and perform final sanity check:
304
305	int totalMolecules = 0;
306	for (int i=0; i < nComponents; i++) {
307	molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
308	totalMolecules += nMol.at(i);
309	}
310	if (totalMolecules != nSites) {
311	sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
312	"to the number of lattice sites!");
313	painCave.isFatal = 1;
314	simError();
315	}
316	} else {
317
318	for (int i = 0; i < shellRadii.size(); i++) {
319	if (shellRadii.at(i) > rodRadius + 1e-6 ) {
320	sprintf(painCave.errMsg, "One of the shellRadius values exceeds the rod Radius.");
321	painCave.isFatal = 1;
322	simError();
323	}
324	if (shellRadii.at(i) <= 0.0 ) {
325	sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!");
326	painCave.isFatal = 1;
327	simError();
328	}
329	}
330	}
331
332	vector<int> ids;
333	if ((int)args_info.molFraction_given){
334	// cerr << "molFraction given 2" << "\n";
335	sprintf(painCave.errMsg, "Creating a randomized spherically-capped nanorod.");
336	painCave.isFatal = 0;
337	simError();
338	/* Random rod is the default case*/
339
340	for (int i = 0; i < sites.size(); i++)
341	if (!isVacancy[i]) ids.push_back(i);
342
343	std::random_shuffle(ids.begin(), ids.end());
344
345	} else{
346	sprintf(painCave.errMsg, "Creating an fcc nanorod.");
347	painCave.isFatal = 0;
348	simError();
349
350	RealType smallestSoFar;
351	int myComponent = -1;
352	nMol.clear();
353	nMol.resize(nComponents);
354
355	// cerr << "shellRadii[0] " << shellRadii[0] << "\n";
356	// cerr << "rodRadius " << rodRadius << "\n";
357
358	for (int i = 0; i < sites.size(); i++) {
359	myLoc = sites[i];
360	myR = myLoc.length();
361	smallestSoFar = rodRadius;
362	//cerr << "vac = " << isVacancy[i]<< "\n";
363
364	if (!isVacancy[i]) {
365
366
367	// for (int j = 0; j < nComponents; j++) {
368	// if (myR <= shellRadii[j]) {
369	// if (shellRadii[j] <= smallestSoFar) {
370	// smallestSoFar = shellRadii[j];
371	// myComponent = j;
372	// }
373	// }
374	// }
375	myComponent = 0;
376	componentFromSite[i] = myComponent;
377	nMol[myComponent]++;
378	// cerr << "nMol for myComp(" << myComponent<<") = " << nMol[myComponent] << "\n";
379	}
380	}
381	}
382	// cerr << "nMol = " << nMol.at(0) << "\n";
383
384	outputFileName = args_info.output_arg;
385
386	//creat new .md file on fly which corrects the number of molecule
387
388	createMdFile(inputFileName, outputFileName, nMol);
389
390	if (oldInfo != NULL)
391	delete oldInfo;
392
393	SimCreator newCreator;
394	SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
395
396	// Place molecules
397	Molecule* mol;
398	SimInfo::MoleculeIterator mi;
399	mol = NewInfo->beginMolecule(mi);
400
401	int l = 0;
402	int whichSite = 0;
403
404	for (int i = 0; i < nComponents; i++){
405	locator = new MoLocator(NewInfo->getMoleculeStamp(i),
406	NewInfo->getForceField());
407
408	// cerr << "nMol = " << nMol.at(i) << "\n";
409	if (!args_info.molFraction_given) {
410	for (int n = 0; n < sites.size(); n++) {
411	if (!isVacancy[n]) {
412	if (componentFromSite[n] == i) {
413	mol = NewInfo->getMoleculeByGlobalIndex(l);
414	locator->placeMol(sites[n], orientations[n], mol);
415	l++;
416	}
417	}
418	}
419	} else {
420	for (int n = 0; n < nMol.at(i); n++) {
421	mol = NewInfo->getMoleculeByGlobalIndex(l);
422	locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
423	l++;
424	}
425	}
426	}
427
428	//fill Hmat
429	hmat(0, 0)= 10.0*rodRadius;
430	hmat(0, 1) = 0.0;
431	hmat(0, 2) = 0.0;
432
433	hmat(1, 0) = 0.0;
434	hmat(1, 1) = 10.0*rodRadius;
435	hmat(1, 2) = 0.0;
436
437	hmat(2, 0) = 0.0;
438	hmat(2, 1) = 0.0;
439	hmat(2, 2) = 5.0rodLength + 2.0rodRadius;
440
441	//set Hmat
442	NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
443
444
445	//create dumpwriter and write out the coordinates
446	writer = new DumpWriter(NewInfo, outputFileName);
447
448	if (writer == NULL) {
449	sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
450	painCave.isFatal = 1;
451	simError();
452	}
453
454	writer->writeDump();
455
456	// deleting the writer will put the closing at the end of the dump file
457
458	delete writer;
459
460	// cleanup a by calling sim error.....
461	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
462	"generated.\n", outputFileName.c_str());
463	painCave.isFatal = 0;
464	simError();
465	return 0;
466	}
467
468	void createMdFile(const std::string&oldMdFileName,
469	const std::string&newMdFileName,
470	std::vector<int> nMol) {
471	ifstream oldMdFile;
472	ofstream newMdFile;
473	const int MAXLEN = 65535;
474	char buffer[MAXLEN];
475
476	//create new .md file based on old .md file
477	oldMdFile.open(oldMdFileName.c_str());
478	newMdFile.open(newMdFileName.c_str());
479	oldMdFile.getline(buffer, MAXLEN);
480
481	int i = 0;
482	while (!oldMdFile.eof()) {
483
484	//correct molecule number
485	if (strstr(buffer, "nMol") != NULL) {
486	if(i<nMol.size()){
487	sprintf(buffer, "\tnMol = %i;", nMol.at(i));
488	newMdFile << buffer << std::endl;
489	i++;
490	}
491	} else
492	newMdFile << buffer << std::endl;
493
494	oldMdFile.getline(buffer, MAXLEN);
495	}
496
497	oldMdFile.close();
498	newMdFile.close();
499
500	if (i != nMol.size()) {
501	sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
502	"\tstatements in component blocks. Make sure that all\n"
503	"\tcomponents in the template file have nMol=1");
504	painCave.isFatal = 1;
505	simError();
506	}
507
508	}
509