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.
#	User	Rev	Content
1	kstocke1	1701	/*
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