applications/nanoparticleBuilder/nanoparticleBuilder.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).
 */

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

#include "config.h"
#include "shapedLatticeSpherical.hpp"
#include "nanoparticleBuilderCmd.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 particleRadius;

  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;
  particleRadius = args_info.radius_arg;
  Globals* simParams = oldInfo->getSimParams();
  
  /* Create nanoparticle */
  shapedLatticeSpherical nanoParticle(latticeConstant, latticeType,
                                      particleRadius);
  
  /* Build a lattice and get lattice points for this lattice constant */
  vector<Vector3d> sites = nanoParticle.getSites();
  vector<Vector3d> orientations = nanoParticle.getOrientations();


  std::vector<int> vacancyTargets;
  vector<bool> isVacancy;
  
  Vector3d myLoc;
  RealType myR;
 
  for (int i = 0; i < sites.size(); i++) 
    isVacancy.push_back(false);

  if (args_info.vacancyPercent_given) {
    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;
      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 = particleRadius;
      }
      if (vIR >= 0.0 && vOR <= particleRadius && 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();

  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();

  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(particleRadius);
  } 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, "nanoparticleBuilder 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(particleRadius);
    } else {    
      sprintf(painCave.errMsg, "nanoparticleBuilder 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) > particleRadius + 1e-6 ) {
        sprintf(painCave.errMsg, "One of the shellRadius values exceeds the particle 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){
    sprintf(painCave.errMsg, "Creating a randomized spherical nanoparticle.");
    painCave.isFatal = 0;
    simError();
    /* Random particle 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 a core-shell spherical nanoparticle.");
    painCave.isFatal = 0;
    simError();

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

    for (int i = 0; i < sites.size(); i++) {
      myLoc = sites[i];
      myR = myLoc.length();
      smallestSoFar = particleRadius;      
      if (!isVacancy[i]) {
        for (int j = 0; j < nComponents; j++) {
          if (myR <= shellRadii[j]) {
            if (shellRadii[j] <= smallestSoFar) {
              smallestSoFar = shellRadii[j];
              myComponent = j;
            }
          }
        }
        componentFromSite[i] = myComponent;
        nMol[myComponent]++;
      }
    }       
  }
  
  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());
    
    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*particleRadius;
  hmat(0, 1) = 0.0;
  hmat(0, 2) = 0.0;
  
  hmat(1, 0) = 0.0;
  hmat(1, 1) =  10.0*particleRadius;
  hmat(1, 2) = 0.0;
  
  hmat(2, 0) = 0.0;
  hmat(2, 1) = 0.0;
  hmat(2, 2) =  10.0*particleRadius;
  
  //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:	15185 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	*/
42
43	#include <cstdlib>
44	#include <cstdio>
45	#include <cstring>
46	#include <cmath>
47	#include <iostream>
48	#include <string>
49	#include <map>
50	#include <fstream>
51	#include <algorithm>
52
53	#include "config.h"
54	#include "shapedLatticeSpherical.hpp"
55	#include "nanoparticleBuilderCmd.h"
56	#include "lattice/LatticeFactory.hpp"
57	#include "utils/MoLocator.hpp"
58	#include "lattice/Lattice.hpp"
59	#include "brains/Register.hpp"
60	#include "brains/SimInfo.hpp"
61	#include "brains/SimCreator.hpp"
62	#include "io/DumpWriter.hpp"
63	#include "math/Vector3.hpp"
64	#include "math/SquareMatrix3.hpp"
65	#include "utils/StringUtils.hpp"
66
67	using namespace std;
68	using namespace OpenMD;
69	void createMdFile(const std::string&oldMdFileName,
70	const std::string&newMdFileName,
71	std::vector<int> numMol);
72
73	int main(int argc, char *argv []) {
74
75	registerLattice();
76
77	gengetopt_args_info args_info;
78	std::string latticeType;
79	std::string inputFileName;
80	std::string outputFileName;
81
82	MoLocator* locator;
83	int nComponents;
84	double latticeConstant;
85	std::vector<double> lc;
86
87	RealType particleRadius;
88
89	Mat3x3d hmat;
90	std::vector<Vector3d> latticePos;
91	std::vector<Vector3d> latticeOrt;
92
93	DumpWriter *writer;
94
95	// Parse Command Line Arguments
96	if (cmdline_parser(argc, argv, &args_info) != 0)
97	exit(1);
98
99	/* get lattice type */
100	latticeType = "FCC";
101
102	/* get input file name */
103	if (args_info.inputs_num)
104	inputFileName = args_info.inputs[0];
105	else {
106	sprintf(painCave.errMsg, "No input .md file name was specified "
107	"on the command line");
108	painCave.isFatal = 1;
109	cmdline_parser_print_help();
110	simError();
111	}
112
113	/* parse md file and set up the system */
114	SimCreator oldCreator;
115	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
116
117	latticeConstant = args_info.latticeConstant_arg;
118	particleRadius = args_info.radius_arg;
119	Globals* simParams = oldInfo->getSimParams();
120
121	/* Create nanoparticle */
122	shapedLatticeSpherical nanoParticle(latticeConstant, latticeType,
123	particleRadius);
124
125	/* Build a lattice and get lattice points for this lattice constant */
126	vector<Vector3d> sites = nanoParticle.getSites();
127	vector<Vector3d> orientations = nanoParticle.getOrientations();
128
129
130	std::vector<int> vacancyTargets;
131	vector<bool> isVacancy;
132
133	Vector3d myLoc;
134	RealType myR;
135
136	for (int i = 0; i < sites.size(); i++)
137	isVacancy.push_back(false);
138
139	if (args_info.vacancyPercent_given) {
140	if (args_info.vacancyPercent_arg < 0.0 \|\| args_info.vacancyPercent_arg > 100.0) {
141	sprintf(painCave.errMsg, "vacancyPercent was set to a non-sensical value.");
142	painCave.isFatal = 1;
143	simError();
144	} else {
145	RealType vF = args_info.vacancyPercent_arg / 100.0;
146	RealType vIR;
147	RealType vOR;
148	if (args_info.vacancyInnerRadius_given) {
149	vIR = args_info.vacancyInnerRadius_arg;
150	} else {
151	vIR = 0.0;
152	}
153	if (args_info.vacancyOuterRadius_given) {
154	vOR = args_info.vacancyOuterRadius_arg;
155	} else {
156	vOR = particleRadius;
157	}
158	if (vIR >= 0.0 && vOR <= particleRadius && vOR >= vIR) {
159
160	for (int i = 0; i < sites.size(); i++) {
161	myLoc = sites[i];
162	myR = myLoc.length();
163	if (myR >= vIR && myR <= vOR) {
164	vacancyTargets.push_back(i);
165	}
166	}
167	std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
168
169	int nTargets = vacancyTargets.size();
170	vacancyTargets.resize((int)(vF * nTargets));
171
172
173	sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
174	"\tsites between %lf and %lf.", (int) vacancyTargets.size(),
175	vIR, vOR);
176	painCave.isFatal = 0;
177	simError();
178
179	isVacancy.clear();
180	for (int i = 0; i < sites.size(); i++) {
181	bool vac = false;
182	for (int j = 0; j < vacancyTargets.size(); j++) {
183	if (i == vacancyTargets[j]) vac = true;
184	}
185	isVacancy.push_back(vac);
186	}
187
188	} else {
189	sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
190	"\tinner or outer radii. Check their values.");
191	painCave.isFatal = 1;
192	simError();
193	}
194	}
195	}
196
197	/* Get number of lattice sites */
198	int nSites = sites.size() - vacancyTargets.size();
199
200	std::vector<Component*> components = simParams->getComponents();
201	std::vector<RealType> molFractions;
202	std::vector<RealType> shellRadii;
203	std::vector<RealType> molecularMasses;
204	std::vector<int> nMol;
205	std::map<int, int> componentFromSite;
206	nComponents = components.size();
207
208	if (args_info.molFraction_given && args_info.shellRadius_given) {
209	sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
210	"arguments, but not both!");
211	painCave.isFatal = 1;
212	simError();
213	}
214
215	if (nComponents == 1) {
216	molFractions.push_back(1.0);
217	shellRadii.push_back(particleRadius);
218	} else if (args_info.molFraction_given) {
219	if ((int)args_info.molFraction_given == nComponents) {
220	for (int i = 0; i < nComponents; i++) {
221	molFractions.push_back(args_info.molFraction_arg[i]);
222	}
223	} else if ((int)args_info.molFraction_given == nComponents-1) {
224	RealType remainingFraction = 1.0;
225	for (int i = 0; i < nComponents-1; i++) {
226	molFractions.push_back(args_info.molFraction_arg[i]);
227	remainingFraction -= molFractions[i];
228	}
229	molFractions.push_back(remainingFraction);
230	} else {
231	sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out molFractions "
232	"for all of the components in the <MetaData> block.");
233	painCave.isFatal = 1;
234	simError();
235	}
236	} else if ((int)args_info.shellRadius_given) {
237	if ((int)args_info.shellRadius_given == nComponents) {
238	for (int i = 0; i < nComponents; i++) {
239	shellRadii.push_back(args_info.shellRadius_arg[i]);
240	}
241	} else if ((int)args_info.shellRadius_given == nComponents-1) {
242	for (int i = 0; i < nComponents-1; i++) {
243	shellRadii.push_back(args_info.shellRadius_arg[i]);
244	}
245	shellRadii.push_back(particleRadius);
246	} else {
247	sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out the\n"
248	"\tshell radii for all of the components in the <MetaData> block.");
249	painCave.isFatal = 1;
250	simError();
251	}
252	} else {
253	sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
254	"\tbut have not specified either molFraction or shellRadius arguments.");
255	painCave.isFatal = 1;
256	simError();
257	}
258
259	if (args_info.molFraction_given) {
260	RealType totalFraction = 0.0;
261
262	/* Do some simple sanity checking*/
263
264	for (int i = 0; i < nComponents; i++) {
265	if (molFractions.at(i) < 0.0) {
266	sprintf(painCave.errMsg, "One of the requested molFractions was"
267	" less than zero!");
268	painCave.isFatal = 1;
269	simError();
270	}
271	if (molFractions.at(i) > 1.0) {
272	sprintf(painCave.errMsg, "One of the requested molFractions was"
273	" greater than one!");
274	painCave.isFatal = 1;
275	simError();
276	}
277	totalFraction += molFractions.at(i);
278	}
279	if (abs(totalFraction - 1.0) > 1e-6) {
280	sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
281	painCave.isFatal = 1;
282	simError();
283	}
284
285	int remaining = nSites;
286	for (int i=0; i < nComponents-1; i++) {
287	nMol.push_back(int((RealType)nSites * molFractions.at(i)));
288	remaining -= nMol.at(i);
289	}
290	nMol.push_back(remaining);
291
292	// recompute actual mol fractions and perform final sanity check:
293
294	int totalMolecules = 0;
295	for (int i=0; i < nComponents; i++) {
296	molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
297	totalMolecules += nMol.at(i);
298	}
299
300	if (totalMolecules != nSites) {
301	sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
302	"to the number of lattice sites!");
303	painCave.isFatal = 1;
304	simError();
305	}
306	} else {
307
308	for (int i = 0; i < shellRadii.size(); i++) {
309	if (shellRadii.at(i) > particleRadius + 1e-6 ) {
310	sprintf(painCave.errMsg, "One of the shellRadius values exceeds the particle Radius.");
311	painCave.isFatal = 1;
312	simError();
313	}
314	if (shellRadii.at(i) <= 0.0 ) {
315	sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!");
316	painCave.isFatal = 1;
317	simError();
318	}
319	}
320	}
321
322	vector<int> ids;
323	if ((int)args_info.molFraction_given){
324	sprintf(painCave.errMsg, "Creating a randomized spherical nanoparticle.");
325	painCave.isFatal = 0;
326	simError();
327	/* Random particle is the default case*/
328
329	for (int i = 0; i < sites.size(); i++)
330	if (!isVacancy[i]) ids.push_back(i);
331
332	std::random_shuffle(ids.begin(), ids.end());
333
334	} else{
335	sprintf(painCave.errMsg, "Creating a core-shell spherical nanoparticle.");
336	painCave.isFatal = 0;
337	simError();
338
339	RealType smallestSoFar;
340	int myComponent = -1;
341	nMol.clear();
342	nMol.resize(nComponents);
343
344	for (int i = 0; i < sites.size(); i++) {
345	myLoc = sites[i];
346	myR = myLoc.length();
347	smallestSoFar = particleRadius;
348	if (!isVacancy[i]) {
349	for (int j = 0; j < nComponents; j++) {
350	if (myR <= shellRadii[j]) {
351	if (shellRadii[j] <= smallestSoFar) {
352	smallestSoFar = shellRadii[j];
353	myComponent = j;
354	}
355	}
356	}
357	componentFromSite[i] = myComponent;
358	nMol[myComponent]++;
359	}
360	}
361	}
362
363	outputFileName = args_info.output_arg;
364
365	//creat new .md file on fly which corrects the number of molecule
366	createMdFile(inputFileName, outputFileName, nMol);
367
368	if (oldInfo != NULL)
369	delete oldInfo;
370
371	SimCreator newCreator;
372	SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
373
374	// Place molecules
375	Molecule* mol;
376	SimInfo::MoleculeIterator mi;
377	mol = NewInfo->beginMolecule(mi);
378
379	int l = 0;
380	int whichSite = 0;
381
382	for (int i = 0; i < nComponents; i++){
383	locator = new MoLocator(NewInfo->getMoleculeStamp(i),
384	NewInfo->getForceField());
385
386	if (!args_info.molFraction_given) {
387	for (int n = 0; n < sites.size(); n++) {
388	if (!isVacancy[n]) {
389	if (componentFromSite[n] == i) {
390	mol = NewInfo->getMoleculeByGlobalIndex(l);
391	locator->placeMol(sites[n], orientations[n], mol);
392	l++;
393	}
394	}
395	}
396	} else {
397	for (int n = 0; n < nMol.at(i); n++) {
398	mol = NewInfo->getMoleculeByGlobalIndex(l);
399	locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
400	l++;
401	}
402	}
403	}
404
405	//fill Hmat
406	hmat(0, 0)= 10.0*particleRadius;
407	hmat(0, 1) = 0.0;
408	hmat(0, 2) = 0.0;
409
410	hmat(1, 0) = 0.0;
411	hmat(1, 1) = 10.0*particleRadius;
412	hmat(1, 2) = 0.0;
413
414	hmat(2, 0) = 0.0;
415	hmat(2, 1) = 0.0;
416	hmat(2, 2) = 10.0*particleRadius;
417
418	//set Hmat
419	NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
420
421
422	//create dumpwriter and write out the coordinates
423	writer = new DumpWriter(NewInfo, outputFileName);
424
425	if (writer == NULL) {
426	sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
427	painCave.isFatal = 1;
428	simError();
429	}
430
431	writer->writeDump();
432
433	// deleting the writer will put the closing at the end of the dump file
434
435	delete writer;
436
437	// cleanup a by calling sim error.....
438	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
439	"generated.\n", outputFileName.c_str());
440	painCave.isFatal = 0;
441	simError();
442	return 0;
443	}
444
445	void createMdFile(const std::string&oldMdFileName,
446	const std::string&newMdFileName,
447	std::vector<int> nMol) {
448	ifstream oldMdFile;
449	ofstream newMdFile;
450	const int MAXLEN = 65535;
451	char buffer[MAXLEN];
452
453	//create new .md file based on old .md file
454	oldMdFile.open(oldMdFileName.c_str());
455	newMdFile.open(newMdFileName.c_str());
456	oldMdFile.getline(buffer, MAXLEN);
457
458	int i = 0;
459	while (!oldMdFile.eof()) {
460
461	//correct molecule number
462	if (strstr(buffer, "nMol") != NULL) {
463	if(i<nMol.size()){
464	sprintf(buffer, "\tnMol = %i;", nMol.at(i));
465	newMdFile << buffer << std::endl;
466	i++;
467	}
468	} else
469	newMdFile << buffer << std::endl;
470
471	oldMdFile.getline(buffer, MAXLEN);
472	}
473
474	oldMdFile.close();
475	newMdFile.close();
476
477	if (i != nMol.size()) {
478	sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
479	"\tstatements in component blocks. Make sure that all\n"
480	"\tcomponents in the template file have nMol=1");
481	painCave.isFatal = 1;
482	simError();
483	}
484
485	}
486