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.
#	User	Rev	Content
1	chuckv	653	/*
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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	chuckv	653	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	chuckv	653	* 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	gezelter	1390	*
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	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	chuckv	653	*/
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	chuckv	911	#include <algorithm>
52	chuckv	653
53			#include "config.h"
54	chuckv	911	#include "shapedLatticeSpherical.hpp"
55	chuckv	653	#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	gezelter	1390	using namespace OpenMD;
69	chuckv	653	void createMdFile(const std::string&oldMdFileName,
70			const std::string&newMdFileName,
71	chuckv	1069	std::vector<int> numMol);
72	chuckv	653
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	chuckv	1069	std::string outputFileName;
81	chuckv	653
82	chuckv	911	MoLocator* locator;
83			int nComponents;
84	chuckv	653	double latticeConstant;
85			std::vector<double> lc;
86	gezelter	1075
87	gezelter	1077	RealType particleRadius;
88	chuckv	653
89			Mat3x3d hmat;
90			std::vector<Vector3d> latticePos;
91			std::vector<Vector3d> latticeOrt;
92	chuckv	1069
93	chuckv	653	DumpWriter *writer;
94
95	chuckv	875	// Parse Command Line Arguments
96	chuckv	653	if (cmdline_parser(argc, argv, &args_info) != 0)
97			exit(1);
98	gezelter	1070
99	chuckv	653	/* get lattice type */
100	chuckv	1069	latticeType = "FCC";
101
102	chuckv	653	/* get input file name */
103			if (args_info.inputs_num)
104			inputFileName = args_info.inputs[0];
105			else {
106	gezelter	1077	sprintf(painCave.errMsg, "No input .md file name was specified "
107	gezelter	1075	"on the command line");
108	chuckv	1069	painCave.isFatal = 1;
109	chuckv	653	cmdline_parser_print_help();
110	chuckv	1069	simError();
111	chuckv	653	}
112
113			/* parse md file and set up the system */
114			SimCreator oldCreator;
115			SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
116
117	gezelter	1077	latticeConstant = args_info.latticeConstant_arg;
118	chuckv	653	particleRadius = args_info.radius_arg;
119	chuckv	911	Globals* simParams = oldInfo->getSimParams();
120	chuckv	653
121	chuckv	911	/* Create nanoparticle */
122	gezelter	1070	shapedLatticeSpherical nanoParticle(latticeConstant, latticeType,
123			particleRadius);
124	chuckv	918
125	chuckv	911	/* Build a lattice and get lattice points for this lattice constant */
126	gezelter	1070	vector<Vector3d> sites = nanoParticle.getSites();
127			vector<Vector3d> orientations = nanoParticle.getOrientations();
128	gezelter	1704
129
130	gezelter	1077	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	chuckv	1069
139	gezelter	1077	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	gezelter	1390	"\tsites between %lf and %lf.", (int) vacancyTargets.size(),
175	gezelter	1077	vIR, vOR);
176			painCave.isFatal = 0;
177			simError();
178	chuckv	1069
179	gezelter	1077	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	chuckv	911	/* Get number of lattice sites */
198	gezelter	1077	int nSites = sites.size() - vacancyTargets.size();
199	chuckv	1069
200			std::vector<Component*> components = simParams->getComponents();
201			std::vector<RealType> molFractions;
202	gezelter	1075	std::vector<RealType> shellRadii;
203	chuckv	1069	std::vector<RealType> molecularMasses;
204			std::vector<int> nMol;
205	gezelter	1075	std::map<int, int> componentFromSite;
206	chuckv	1069	nComponents = components.size();
207
208	gezelter	1077	if (args_info.molFraction_given && args_info.shellRadius_given) {
209			sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
210	gezelter	1075	"arguments, but not both!");
211			painCave.isFatal = 1;
212			simError();
213			}
214
215			if (nComponents == 1) {
216	chuckv	1069	molFractions.push_back(1.0);
217	gezelter	1075	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	chuckv	1069	painCave.isFatal = 1;
234			simError();
235			}
236	gezelter	1077	} else if ((int)args_info.shellRadius_given) {
237			if ((int)args_info.shellRadius_given == nComponents) {
238	gezelter	1075	for (int i = 0; i < nComponents; i++) {
239	gezelter	1077	shellRadii.push_back(args_info.shellRadius_arg[i]);
240	gezelter	1075	}
241	gezelter	1077	} else if ((int)args_info.shellRadius_given == nComponents-1) {
242	gezelter	1075	for (int i = 0; i < nComponents-1; i++) {
243	gezelter	1077	shellRadii.push_back(args_info.shellRadius_arg[i]);
244	gezelter	1075	}
245			shellRadii.push_back(particleRadius);
246			} else {
247	gezelter	1077	sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out the\n"
248			"\tshell radii for all of the components in the <MetaData> block.");
249	chuckv	1069	painCave.isFatal = 1;
250			simError();
251			}
252	gezelter	1075	} else {
253	gezelter	1077	sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
254			"\tbut have not specified either molFraction or shellRadius arguments.");
255	chuckv	1069	painCave.isFatal = 1;
256			simError();
257			}
258	gezelter	1075
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	chuckv	1069
308	gezelter	1075	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	chuckv	1069	}
321	gezelter	1077
322			vector<int> ids;
323	gezelter	1075	if ((int)args_info.molFraction_given){
324			sprintf(painCave.errMsg, "Creating a randomized spherical nanoparticle.");
325			painCave.isFatal = 0;
326			simError();
327	gezelter	1077	/* 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	chuckv	1069	std::random_shuffle(ids.begin(), ids.end());
333	gezelter	1077
334	gezelter	1075	} else{
335			sprintf(painCave.errMsg, "Creating a core-shell spherical nanoparticle.");
336			painCave.isFatal = 0;
337			simError();
338	chuckv	653
339	gezelter	1075	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	gezelter	1077	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	gezelter	1075	}
356			}
357	gezelter	1077	componentFromSite[i] = myComponent;
358			nMol[myComponent]++;
359	gezelter	1075	}
360	gezelter	1077	}
361			}
362	chuckv	949
363	chuckv	1069	outputFileName = args_info.output_arg;
364	gezelter	1077
365	gezelter	1075	//creat new .md file on fly which corrects the number of molecule
366	chuckv	1069	createMdFile(inputFileName, outputFileName, nMol);
367	chuckv	653
368			if (oldInfo != NULL)
369			delete oldInfo;
370
371	chuckv	949	SimCreator newCreator;
372	chuckv	1069	SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
373	gezelter	1075
374	chuckv	911	// Place molecules
375	chuckv	653	Molecule* mol;
376			SimInfo::MoleculeIterator mi;
377			mol = NewInfo->beginMolecule(mi);
378	gezelter	1077
379	chuckv	949	int l = 0;
380	gezelter	1077	int whichSite = 0;
381	chuckv	1069
382			for (int i = 0; i < nComponents; i++){
383			locator = new MoLocator(NewInfo->getMoleculeStamp(i),
384			NewInfo->getForceField());
385	gezelter	1077
386			if (!args_info.molFraction_given) {
387	gezelter	1075	for (int n = 0; n < sites.size(); n++) {
388	gezelter	1077	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	gezelter	1075	}
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	chuckv	1069	}
403	gezelter	1077	}
404	chuckv	653
405			//fill Hmat
406	gezelter	1075	hmat(0, 0)= 10.0*particleRadius;
407	chuckv	653	hmat(0, 1) = 0.0;
408			hmat(0, 2) = 0.0;
409
410			hmat(1, 0) = 0.0;
411	gezelter	1075	hmat(1, 1) = 10.0*particleRadius;
412	chuckv	653	hmat(1, 2) = 0.0;
413
414			hmat(2, 0) = 0.0;
415			hmat(2, 1) = 0.0;
416	gezelter	1075	hmat(2, 2) = 10.0*particleRadius;
417	chuckv	653
418			//set Hmat
419			NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
420
421
422			//create dumpwriter and write out the coordinates
423	gezelter	1070	writer = new DumpWriter(NewInfo, outputFileName);
424	chuckv	653
425			if (writer == NULL) {
426	gezelter	1077	sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
427	gezelter	1075	painCave.isFatal = 1;
428			simError();
429	chuckv	653	}
430
431			writer->writeDump();
432	chuckv	1069
433			// deleting the writer will put the closing at the end of the dump file
434	gezelter	1070
435	chuckv	1069	delete writer;
436
437			// cleanup a by calling sim error.....
438	gezelter	1390	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
439	chuckv	1069	"generated.\n", outputFileName.c_str());
440			painCave.isFatal = 0;
441			simError();
442	chuckv	653	return 0;
443			}
444
445	gezelter	1075	void createMdFile(const std::string&oldMdFileName,
446			const std::string&newMdFileName,
447	chuckv	1069	std::vector<int> nMol) {
448	chuckv	653	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	gezelter	1077
458	chuckv	911	int i = 0;
459	chuckv	653	while (!oldMdFile.eof()) {
460	gezelter	1077
461	chuckv	653	//correct molecule number
462			if (strstr(buffer, "nMol") != NULL) {
463	chuckv	1069	if(i<nMol.size()){
464	gezelter	1077	sprintf(buffer, "\tnMol = %i;", nMol.at(i));
465	chuckv	949	newMdFile << buffer << std::endl;
466			i++;
467			}
468	chuckv	653	} else
469			newMdFile << buffer << std::endl;
470
471			oldMdFile.getline(buffer, MAXLEN);
472			}
473
474			oldMdFile.close();
475			newMdFile.close();
476	gezelter	1077
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	chuckv	653	}
486