applications/nanoparticleBuilder/nanorod_pentBuilder.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 4/5/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 "shapedLatticePentRod.hpp"
#include "nanorod_pentBuilderCmd.h"
#include "shapedLatticeRod.hpp"
#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/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 */
  shapedLatticePentRod nanoRod(latticeConstant, latticeType,
                           rodRadius, rodLength);
  
  /* Build a lattice and get lattice points for this lattice constant */

  //Rotation angles for lattice
  RealType phi, theta, psi;

  RealType cphi, sphi, ctheta, stheta, cpsi, spsi;

    
  /*cphi = cos(phi);
  sphi = sin(phi);
  ctheta = cos(theta);
  stheta = sin(theta);
  cpsi = cos(psi);
  spsi = sin(psi);*/

  //Rotates 45 degrees about z-axis
  RotMat3x3d rotation45( 45.0 * M_PI / 180.0, 0.0, 0.0);
 
  /*rotation45[0][0] = sqrt(2)/2;
  rotation45[0][1] = -sqrt(2)/2;
  rotation45[0][2] = 0;
  rotation45[1][0] = sqrt(2)/2;
  rotation45[1][1] = sqrt(2)/2;
  rotation45[1][2] = 0;
  rotation45[2][0] = 0;
  rotation45[2][1] = 0;
  rotation45[2][2] = 1;*/

  phi =  0.0;
  theta = 72 * M_PI / 180.0;
  psi = 0.0;

  //Rotates 72 degrees about y-axis
  RotMat3x3d rotation72(phi, theta, psi);

  /*rotation72[0][0] = sqrt(5)/4 - 0.25;
  rotation72[0][1] = 0;
  rotation72[0][2] = sqrt(2*(sqrt(5) + 5))/4;
  rotation72[1][0] = 0;
  rotation72[1][1] = 1;
  rotation72[1][2] = 0;
  rotation72[2][0] = -sqrt(2*(sqrt(5) + 5))/4;
  rotation72[2][1] = 0;
  rotation72[2][2] = sqrt(5)/4 - 0.25;*/

  vector<Vector3d> getsites = nanoRod.getSites();
  vector<Vector3d> getorientations = nanoRod.getOrientations();
  vector<Vector3d> sites;
  vector<Vector3d> orientations;
 
  int index;

  for (index = 0; index < getsites.size(); index++) {
    Vector3d mySite = getsites[index];
    Vector3d myOrient = getorientations[index];
    Vector3d mySite2 = rotation45 * mySite;
    Vector3d o2 = rotation45 * myOrient;
    sites.push_back(  mySite2 );
    orientations.push_back( o2 );

    mySite2 = rotation72 * mySite2;
    o2 = rotation72 * o2;
    sites.push_back( mySite2 );
    orientations.push_back( o2 );

    mySite2 = rotation72 * mySite2;
    o2 = rotation72 * o2;
    sites.push_back( mySite2 );
    orientations.push_back( o2 );

    mySite2 = rotation72 * mySite2;
    o2 = rotation72 * o2;
    sites.push_back( mySite2 );
    orientations.push_back( o2 );

    mySite2 = rotation72 * mySite2;
    o2 = rotation72 * o2;
    sites.push_back( mySite2 );
    orientations.push_back( o2 );
  }

  int nCenter = (rodLength + 1.154700538*rodRadius)/2.88;

  for (index = 0; index <= 0.5*nCenter; index++) {
    Vector3d myLoc_top(2.88*index, 0.0, 0.0);
    sites.push_back(myLoc_top);
    orientations.push_back(Vector3d(0.0));
  }

  for (index = 1; index <= 0.5*nCenter; index++) {
    Vector3d myLoc_bottom(-2.88*index, 0.0, 0.0);
    sites.push_back(myLoc_bottom);
    orientations.push_back(Vector3d(0.0));
  }

  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
Original Path:	*branches/development/src/applications/nanoparticleBuilder/nanorod_pentBuilder.cpp*
File size:	18335 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 4/5/12.
42			* @author Kelsey M. Stocker
43			*
44			*/
45
46
47
48			#include <cstdlib>
49			#include <cstdio>
50			#include <cstring>
51			#include <cmath>
52			#include <iostream>
53			#include <string>
54			#include <map>
55			#include <fstream>
56			#include <algorithm>
57
58			#include "config.h"
59			#include "shapedLatticePentRod.hpp"
60			#include "nanorod_pentBuilderCmd.h"
61			#include "shapedLatticeRod.hpp"
62			#include "lattice/LatticeFactory.hpp"
63			#include "utils/MoLocator.hpp"
64			#include "lattice/Lattice.hpp"
65			#include "brains/Register.hpp"
66			#include "brains/SimInfo.hpp"
67			#include "brains/SimCreator.hpp"
68			#include "io/DumpWriter.hpp"
69			#include "math/SquareMatrix3.hpp"
70			#include "utils/StringUtils.hpp"
71
72			using namespace std;
73			using namespace OpenMD;
74			void createMdFile(const std::string&oldMdFileName,
75			const std::string&newMdFileName,
76			std::vector<int> numMol);
77
78			int main(int argc, char *argv []) {
79
80			registerLattice();
81
82			gengetopt_args_info args_info;
83			std::string latticeType;
84			std::string inputFileName;
85			std::string outputFileName;
86
87			MoLocator* locator;
88			int nComponents;
89			double latticeConstant;
90			std::vector<double> lc;
91
92			RealType rodRadius;
93			RealType rodLength;
94
95			Mat3x3d hmat;
96			std::vector<Vector3d> latticePos;
97			std::vector<Vector3d> latticeOrt;
98
99			DumpWriter *writer;
100
101			// Parse Command Line Arguments
102			if (cmdline_parser(argc, argv, &args_info) != 0)
103			exit(1);
104
105			/* get lattice type */
106			latticeType = "FCC";
107
108			/* get input file name */
109			if (args_info.inputs_num)
110			inputFileName = args_info.inputs[0];
111			else {
112			sprintf(painCave.errMsg, "No input .md file name was specified "
113			"on the command line");
114			painCave.isFatal = 1;
115			cmdline_parser_print_help();
116			simError();
117			}
118
119			/* parse md file and set up the system */
120			SimCreator oldCreator;
121			SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
122
123			latticeConstant = args_info.latticeConstant_arg;
124			rodRadius = args_info.radius_arg;
125			rodLength = args_info.length_arg;
126			Globals* simParams = oldInfo->getSimParams();
127
128			/* Create nanorod */
129			shapedLatticePentRod nanoRod(latticeConstant, latticeType,
130			rodRadius, rodLength);
131
132			/* Build a lattice and get lattice points for this lattice constant */
133
134			//Rotation angles for lattice
135			RealType phi, theta, psi;
136
137			RealType cphi, sphi, ctheta, stheta, cpsi, spsi;
138
139
140			/*cphi = cos(phi);
141			sphi = sin(phi);
142			ctheta = cos(theta);
143			stheta = sin(theta);
144			cpsi = cos(psi);
145			spsi = sin(psi);*/
146
147			//Rotates 45 degrees about z-axis
148			RotMat3x3d rotation45( 45.0 * M_PI / 180.0, 0.0, 0.0);
149
150			/*rotation45[0][0] = sqrt(2)/2;
151			rotation45[0][1] = -sqrt(2)/2;
152			rotation45[0][2] = 0;
153			rotation45[1][0] = sqrt(2)/2;
154			rotation45[1][1] = sqrt(2)/2;
155			rotation45[1][2] = 0;
156			rotation45[2][0] = 0;
157			rotation45[2][1] = 0;
158			rotation45[2][2] = 1;*/
159
160			phi = 0.0;
161			theta = 72 * M_PI / 180.0;
162			psi = 0.0;
163
164			//Rotates 72 degrees about y-axis
165			RotMat3x3d rotation72(phi, theta, psi);
166
167			/*rotation72[0][0] = sqrt(5)/4 - 0.25;
168			rotation72[0][1] = 0;
169			rotation72[0][2] = sqrt(2*(sqrt(5) + 5))/4;
170			rotation72[1][0] = 0;
171			rotation72[1][1] = 1;
172			rotation72[1][2] = 0;
173			rotation72[2][0] = -sqrt(2*(sqrt(5) + 5))/4;
174			rotation72[2][1] = 0;
175			rotation72[2][2] = sqrt(5)/4 - 0.25;*/
176
177			vector<Vector3d> getsites = nanoRod.getSites();
178			vector<Vector3d> getorientations = nanoRod.getOrientations();
179			vector<Vector3d> sites;
180			vector<Vector3d> orientations;
181
182			int index;
183
184			for (index = 0; index < getsites.size(); index++) {
185			Vector3d mySite = getsites[index];
186			Vector3d myOrient = getorientations[index];
187			Vector3d mySite2 = rotation45 * mySite;
188			Vector3d o2 = rotation45 * myOrient;
189			sites.push_back( mySite2 );
190			orientations.push_back( o2 );
191
192			mySite2 = rotation72 * mySite2;
193			o2 = rotation72 * o2;
194			sites.push_back( mySite2 );
195			orientations.push_back( o2 );
196
197			mySite2 = rotation72 * mySite2;
198			o2 = rotation72 * o2;
199			sites.push_back( mySite2 );
200			orientations.push_back( o2 );
201
202			mySite2 = rotation72 * mySite2;
203			o2 = rotation72 * o2;
204			sites.push_back( mySite2 );
205			orientations.push_back( o2 );
206
207			mySite2 = rotation72 * mySite2;
208			o2 = rotation72 * o2;
209			sites.push_back( mySite2 );
210			orientations.push_back( o2 );
211			}
212
213			int nCenter = (rodLength + 1.154700538*rodRadius)/2.88;
214
215			for (index = 0; index <= 0.5*nCenter; index++) {
216			Vector3d myLoc_top(2.88*index, 0.0, 0.0);
217			sites.push_back(myLoc_top);
218			orientations.push_back(Vector3d(0.0));
219			}
220
221			for (index = 1; index <= 0.5*nCenter; index++) {
222			Vector3d myLoc_bottom(-2.88*index, 0.0, 0.0);
223			sites.push_back(myLoc_bottom);
224			orientations.push_back(Vector3d(0.0));
225			}
226
227			std::vector<int> vacancyTargets;
228			vector<bool> isVacancy;
229
230			Vector3d myLoc;
231			RealType myR;
232
233			for (int i = 0; i < sites.size(); i++)
234			isVacancy.push_back(false);
235
236			// cerr << "checking vacancyPercent" << "\n";
237			if (args_info.vacancyPercent_given) {
238			// cerr << "vacancyPercent given" << "\n";
239			if (args_info.vacancyPercent_arg < 0.0 \|\| args_info.vacancyPercent_arg > 100.0) {
240			sprintf(painCave.errMsg, "vacancyPercent was set to a non-sensical value.");
241			painCave.isFatal = 1;
242			simError();
243			} else {
244			RealType vF = args_info.vacancyPercent_arg / 100.0;
245			// cerr << "vacancyPercent = " << vF << "\n";
246			RealType vIR;
247			RealType vOR;
248			if (args_info.vacancyInnerRadius_given) {
249			vIR = args_info.vacancyInnerRadius_arg;
250			} else {
251			vIR = 0.0;
252			}
253			if (args_info.vacancyOuterRadius_given) {
254			vOR = args_info.vacancyOuterRadius_arg;
255			} else {
256			vOR = rodRadius;
257			}
258			if (vIR >= 0.0 && vOR <= rodRadius && vOR >= vIR) {
259
260			for (int i = 0; i < sites.size(); i++) {
261			myLoc = sites[i];
262			myR = myLoc.length();
263			if (myR >= vIR && myR <= vOR) {
264			vacancyTargets.push_back(i);
265			}
266			}
267			std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end());
268
269			int nTargets = vacancyTargets.size();
270			vacancyTargets.resize((int)(vF * nTargets));
271
272
273			sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n"
274			"\tsites between %lf and %lf.", (int) vacancyTargets.size(),
275			vIR, vOR);
276			painCave.isFatal = 0;
277			simError();
278
279			isVacancy.clear();
280			for (int i = 0; i < sites.size(); i++) {
281			bool vac = false;
282			for (int j = 0; j < vacancyTargets.size(); j++) {
283			if (i == vacancyTargets[j]) vac = true;
284			}
285			isVacancy.push_back(vac);
286			}
287
288			} else {
289			sprintf(painCave.errMsg, "Something is strange about the vacancy\n"
290			"\tinner or outer radii. Check their values.");
291			painCave.isFatal = 1;
292			simError();
293			}
294			}
295			}
296
297			/* Get number of lattice sites */
298			int nSites = sites.size() - vacancyTargets.size();
299
300			// cerr << "sites.size() = " << sites.size() << "\n";
301			// cerr << "nSites = " << nSites << "\n";
302			// cerr << "vacancyTargets = " << vacancyTargets.size() << "\n";
303
304			std::vector<Component*> components = simParams->getComponents();
305			std::vector<RealType> molFractions;
306			std::vector<RealType> shellRadii;
307			std::vector<RealType> molecularMasses;
308			std::vector<int> nMol;
309			std::map<int, int> componentFromSite;
310			nComponents = components.size();
311			// cerr << "nComponents = " << nComponents << "\n";
312
313			if (args_info.molFraction_given && args_info.shellRadius_given) {
314			sprintf(painCave.errMsg, "Specify either molFraction or shellRadius "
315			"arguments, but not both!");
316			painCave.isFatal = 1;
317			simError();
318			}
319
320			if (nComponents == 1) {
321			molFractions.push_back(1.0);
322			shellRadii.push_back(rodRadius);
323			} else if (args_info.molFraction_given) {
324			if ((int)args_info.molFraction_given == nComponents) {
325			for (int i = 0; i < nComponents; i++) {
326			molFractions.push_back(args_info.molFraction_arg[i]);
327			}
328			} else if ((int)args_info.molFraction_given == nComponents-1) {
329			RealType remainingFraction = 1.0;
330			for (int i = 0; i < nComponents-1; i++) {
331			molFractions.push_back(args_info.molFraction_arg[i]);
332			remainingFraction -= molFractions[i];
333			}
334			molFractions.push_back(remainingFraction);
335			} else {
336			sprintf(painCave.errMsg, "nanorodBuilder can't figure out molFractions "
337			"for all of the components in the <MetaData> block.");
338			painCave.isFatal = 1;
339			simError();
340			}
341			} else if ((int)args_info.shellRadius_given) {
342			if ((int)args_info.shellRadius_given == nComponents) {
343			for (int i = 0; i < nComponents; i++) {
344			shellRadii.push_back(args_info.shellRadius_arg[i]);
345			}
346			} else if ((int)args_info.shellRadius_given == nComponents-1) {
347			for (int i = 0; i < nComponents-1; i++) {
348			shellRadii.push_back(args_info.shellRadius_arg[i]);
349			}
350			shellRadii.push_back(rodRadius);
351			} else {
352			sprintf(painCave.errMsg, "nanorodBuilder can't figure out the\n"
353			"\tshell radii for all of the components in the <MetaData> block.");
354			painCave.isFatal = 1;
355			simError();
356			}
357			} else {
358			sprintf(painCave.errMsg, "You have a multi-component <MetaData> block,\n"
359			"\tbut have not specified either molFraction or shellRadius arguments.");
360			painCave.isFatal = 1;
361			simError();
362			}
363
364			if (args_info.molFraction_given) {
365			RealType totalFraction = 0.0;
366
367			/* Do some simple sanity checking*/
368
369			for (int i = 0; i < nComponents; i++) {
370			if (molFractions.at(i) < 0.0) {
371			sprintf(painCave.errMsg, "One of the requested molFractions was"
372			" less than zero!");
373			painCave.isFatal = 1;
374			simError();
375			}
376			if (molFractions.at(i) > 1.0) {
377			sprintf(painCave.errMsg, "One of the requested molFractions was"
378			" greater than one!");
379			painCave.isFatal = 1;
380			simError();
381			}
382			totalFraction += molFractions.at(i);
383			}
384			if (abs(totalFraction - 1.0) > 1e-6) {
385			sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
386			painCave.isFatal = 1;
387			simError();
388			}
389
390			int remaining = nSites;
391			for (int i=0; i < nComponents-1; i++) {
392			nMol.push_back(int((RealType)nSites * molFractions.at(i)));
393			remaining -= nMol.at(i);
394			}
395			nMol.push_back(remaining);
396
397			// recompute actual mol fractions and perform final sanity check:
398
399			int totalMolecules = 0;
400			for (int i=0; i < nComponents; i++) {
401			molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
402			totalMolecules += nMol.at(i);
403			}
404			if (totalMolecules != nSites) {
405			sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
406			"to the number of lattice sites!");
407			painCave.isFatal = 1;
408			simError();
409			}
410			} else {
411
412			for (int i = 0; i < shellRadii.size(); i++) {
413			if (shellRadii.at(i) > rodRadius + 1e-6 ) {
414			sprintf(painCave.errMsg, "One of the shellRadius values exceeds the rod Radius.");
415			painCave.isFatal = 1;
416			simError();
417			}
418			if (shellRadii.at(i) <= 0.0 ) {
419			sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!");
420			painCave.isFatal = 1;
421			simError();
422			}
423			}
424			}
425
426			vector<int> ids;
427			if ((int)args_info.molFraction_given){
428			// cerr << "molFraction given 2" << "\n";
429			sprintf(painCave.errMsg, "Creating a randomized spherically-capped nanorod.");
430			painCave.isFatal = 0;
431			simError();
432			/* Random rod is the default case*/
433
434			for (int i = 0; i < sites.size(); i++)
435			if (!isVacancy[i]) ids.push_back(i);
436
437			std::random_shuffle(ids.begin(), ids.end());
438
439			} else{
440			sprintf(painCave.errMsg, "Creating an fcc nanorod.");
441			painCave.isFatal = 0;
442			simError();
443
444			RealType smallestSoFar;
445			int myComponent = -1;
446			nMol.clear();
447			nMol.resize(nComponents);
448
449			// cerr << "shellRadii[0] " << shellRadii[0] << "\n";
450			// cerr << "rodRadius " << rodRadius << "\n";
451
452			for (int i = 0; i < sites.size(); i++) {
453			myLoc = sites[i];
454			myR = myLoc.length();
455			smallestSoFar = rodRadius;
456			//cerr << "vac = " << isVacancy[i]<< "\n";
457
458			if (!isVacancy[i]) {
459
460
461			// for (int j = 0; j < nComponents; j++) {
462			// if (myR <= shellRadii[j]) {
463			// if (shellRadii[j] <= smallestSoFar) {
464			// smallestSoFar = shellRadii[j];
465			// myComponent = j;
466			// }
467			// }
468			// }
469			myComponent = 0;
470			componentFromSite[i] = myComponent;
471			nMol[myComponent]++;
472			// cerr << "nMol for myComp(" << myComponent<<") = " << nMol[myComponent] << "\n";
473			}
474			}
475			}
476			// cerr << "nMol = " << nMol.at(0) << "\n";
477
478			outputFileName = args_info.output_arg;
479
480			//creat new .md file on fly which corrects the number of molecule
481
482			createMdFile(inputFileName, outputFileName, nMol);
483
484			if (oldInfo != NULL)
485			delete oldInfo;
486
487			SimCreator newCreator;
488			SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
489
490			// Place molecules
491			Molecule* mol;
492			SimInfo::MoleculeIterator mi;
493			mol = NewInfo->beginMolecule(mi);
494
495			int l = 0;
496			int whichSite = 0;
497
498			for (int i = 0; i < nComponents; i++){
499			locator = new MoLocator(NewInfo->getMoleculeStamp(i),
500			NewInfo->getForceField());
501
502			// cerr << "nMol = " << nMol.at(i) << "\n";
503			if (!args_info.molFraction_given) {
504			for (int n = 0; n < sites.size(); n++) {
505			if (!isVacancy[n]) {
506			if (componentFromSite[n] == i) {
507			mol = NewInfo->getMoleculeByGlobalIndex(l);
508			locator->placeMol(sites[n], orientations[n], mol);
509			l++;
510			}
511			}
512			}
513			} else {
514			for (int n = 0; n < nMol.at(i); n++) {
515			mol = NewInfo->getMoleculeByGlobalIndex(l);
516			locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
517			l++;
518			}
519			}
520			}
521
522			//fill Hmat
523			hmat(0, 0)= 10.0*rodRadius;
524			hmat(0, 1) = 0.0;
525			hmat(0, 2) = 0.0;
526
527			hmat(1, 0) = 0.0;
528			hmat(1, 1) = 10.0*rodRadius;
529			hmat(1, 2) = 0.0;
530
531			hmat(2, 0) = 0.0;
532			hmat(2, 1) = 0.0;
533			hmat(2, 2) = 5.0rodLength + 2.0rodRadius;
534
535			//set Hmat
536			NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
537
538
539			//create dumpwriter and write out the coordinates
540			writer = new DumpWriter(NewInfo, outputFileName);
541
542			if (writer == NULL) {
543			sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
544			painCave.isFatal = 1;
545			simError();
546			}
547
548			writer->writeDump();
549
550			// deleting the writer will put the closing at the end of the dump file
551
552			delete writer;
553
554			// cleanup a by calling sim error.....
555			sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
556			"generated.\n", outputFileName.c_str());
557			painCave.isFatal = 0;
558			simError();
559			return 0;
560			}
561
562			void createMdFile(const std::string&oldMdFileName,
563			const std::string&newMdFileName,
564			std::vector<int> nMol) {
565			ifstream oldMdFile;
566			ofstream newMdFile;
567			const int MAXLEN = 65535;
568			char buffer[MAXLEN];
569
570			//create new .md file based on old .md file
571			oldMdFile.open(oldMdFileName.c_str());
572			newMdFile.open(newMdFileName.c_str());
573			oldMdFile.getline(buffer, MAXLEN);
574
575			int i = 0;
576			while (!oldMdFile.eof()) {
577
578			//correct molecule number
579			if (strstr(buffer, "nMol") != NULL) {
580			if(i<nMol.size()){
581			sprintf(buffer, "\tnMol = %i;", nMol.at(i));
582			newMdFile << buffer << std::endl;
583			i++;
584			}
585			} else
586			newMdFile << buffer << std::endl;
587
588			oldMdFile.getline(buffer, MAXLEN);
589			}
590
591			oldMdFile.close();
592			newMdFile.close();
593
594			if (i != nMol.size()) {
595			sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
596			"\tstatements in component blocks. Make sure that all\n"
597			"\tcomponents in the template file have nMol=1");
598			painCave.isFatal = 1;
599			simError();
600			}
601
602			}
603