OOPSE/utils/nanoBuilder.cpp

#include <iostream>

#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>

#include "simError.h"
#include "SimInfo.hpp"
#include "ReadWrite.hpp"

#include "latticeBuilder.hpp"
#include "MoLocator.hpp"
#include "sysBuild.hpp"
#include "nanoBuilder.hpp"


nanoBuilder::nanoBuilder(int thisIsRandom, int thisHasVacancies,
                         int thisLatticeType, double thisParticleRadius, 
                         double thisCoreRadius, double thisVacancyFraction,
                         double thisVacancyRadius,
                         double thisLatticeSpacing, 
                         double solute_X,
                         int &hasError){
  int Errors;
  int foundCore,foundShell;
  int i;

  //Zero variables
  particleRadius = 0.0;
  coreRadius = 0.0;
  vacancyFraction = 0.0;
  vacancyRadius = 0.0;
  shellRadius = 0.0;
  latticeSpacing = 0.0;

  buildNmol = 0;

  nCoreMolecules = 0;
  nShellMolecules = 0;

  atomCount = 0;
  coreAtomCount = 0;
  shellAtomCount = 0;

  
  moleculeCount = 0; 
  foundCore  = 0;
  foundShell = 0;
  totalMolecules = 0;
  coreHasOrientation = 0;
  shellHasOrientation = 0;
  nInterface = 0;
  nMol = 0;

  hasError = 0;
  Errors = 0;


  isRandom        = thisIsRandom;
  hasVacancies    = thisHasVacancies;
  latticeType     = thisLatticeType;
  particleRadius  = thisParticleRadius;
  coreRadius      = thisCoreRadius;
  vacancyFraction = thisVacancyFraction;
  latticeSpacing  = thisLatticeSpacing;
  soluteX = solute_X; //Mole fraction for random particle.


  for (i=0;bsInfo.nComponents;i++){
    if( !strcmp( bsInfo.compStamps[i]->getID(),bsInfo.coreName )){
      foundCore = 1;
      coreStamp = bsInfo.compStamps[i];
      nCoreMolecules = bsInfo.componentsNmol[i];
    }
    if( !strcmp( bsInfo.compStamps[i]->getID(),bsInfo.shellName)){
      foundShell = 1;
      shellStamp = bsInfo.compStamps[i];
      nShellMolecules = bsInfo.componentsNmol[i];

    }    

  }


  if( !foundCore ){
    hasError = 1;
    return;
  }
  if( !foundShell ){
    hasError = 1;
    return;
  }


  Errors = sanityCheck();

  if (Errors){
    hasError = 1;
    return;
  }


  nCoreModelAtoms  = coreStamp->getNAtoms();
  nShellModelAtoms = shellStamp->getNAtoms();
  

  // We assume that if the core or shell model has more then one atom
  // the model has an orientational component...
  if (nCoreModelAtoms > 1)    coreHasOrientation = 1;
  if (nShellModelAtoms > 1)   shellHasOrientation = 1;
  
  maxModelNatoms = std::max(nCoreModelAtoms,nShellModelAtoms);
  
  /* If we specify a number of atoms in bass, we will try to build a nanopartice
     with that number.
  */
 

  if ((nShellMolecules != 0) && (nCoreMolecules != 0)){
    totalMolecules = nShellMolecules + nCoreMolecules;
    nCells = ceil(pow((double)totalMolecules/4.0, 1/3));
    buildNmol = 1;
  }
  else {
    nCells = 2.0 * particleRadius/latticeSpacing;
    shellRadius = particleRadius - coreRadius;
  }


  // Initialize random seed
  srand48( RAND_SEED );

  
}


nanoBuilder::~nanoBuilder(){
}


// Checks to make sure we aren't doing something the builder can't do.
int nanoBuilder::sanityCheck(void){

  // Right now we only do bimetallic nanoparticles  
  if (bsInfo.nComponents > 2) return 1;

  //Check for vacancies and random
  if (hasVacancies && isRandom) return 1;

  // make sure we aren't trying to build a core larger then the total particle size
  if ((coreRadius >= particleRadius) && (particleRadius != 0)) return 1;

  // we initialize the lattice spacing to be 0.0, if the lattice spacing is still 0.0
  // we have a problem
  if (latticeSpacing == 0.0) return 1;

  // Check to see if we are specifing the number of atoms in the particle correctly.
  if ((nShellMolecules == 0) && (nCoreMolecules != 0)){
    cerr << "nShellParticles is zero and nCoreParticles != 0" << "\n";
    return 1;
  } 
  // Make sure there are more then two components if we are building a randomly mixed particle.
  if ((bsInfo.nComponents < 2) && (isRandom)){
    cerr << "Two Components are needed to build a random particle." << "\n";
  }
  // Make sure both the core and shell models specify a target nmol.
  if ((nShellMolecules != 0) && (nCoreMolecules == 0)){
    cerr << "nCoreParticles is zero and nShellParticles != 0" << "\n";
    return 1;
  } 
  
  return 0;

}

    
int nanoBuilder::buildNanoParticle( void ){

  int ix;
  int iy;
  int iz;
  double *rx;
  double *ry;
  double *rz;
  double pos[3];
  double A[3][3];
  double HmatI[9];
  
  int nCellSites;
  int iref;
  int appNMols;
  int latticeCount = 0;
 
  int nAtoms;
  int nCoreAtomCounter = 0;
  int nShellAtomCounter = 0;
  int hasError;

  int i;

  int interfaceIndex = 0;
  double dist;
  double distsq;
  int latticeNpoints;
  int shesActualSizetoMe = 0;

  DumpWriter* writer;
  SimInfo* simnfo;

  Lattice::Lattice *myLattice;
  MoLocator::MoLocator *coreLocate;
  MoLocator::MoLocator *shellLocate;


  Atom** atoms;

  hasError = 0;

  myLattice = new Lattice(FCC_LATTICE_TYPE,latticeSpacing);
  /*
  latticeNpoints = myLattice.getNpoints();

  // Initializd atom vector to approximate size. 
  switch (buildType){

  case BUILD_NMOL_PARTICLE:

    break;
  case BUILD_CORE_SHELL_VACANCY:
   // Make space in the vector for all atoms except the last full cells
    // We will have to add at most (latticeNpoints-1)^3 to vector 
    appNMols = latticeNPoints * pow((double)(nCells - 1),3);
    moleculeVector.pushBack();

  default:
    // Make space in the vector for all atoms except the last full cells
    // We will have to add at most (latticeNpoints-1)^3 to vector 
    appNMols = latticeNPoints * pow((double)(nCells - 1),3);

  }
  */
  

  // Create molocator and atom arrays.
  coreLocate  = new MoLocator(coreStamp);
  shellLocate = new MoLocator(shellStamp);
 
  
  for(iz=-nCells;iz < nCells;iz++){ 
    for(iy=-nCells;iy<nCells;iy++){      
      for(ix=-nCells;ix<nCells;ix++){
        nCellSites = myLattice->getLatticePoints(&rx,&ry,&rz,
                                                 ix,iy,iz);
        for (iref=1;iref<nCellSites;iref++){
          latticeCount++;
          
          pos[0] = rx[iref];
          pos[1] = ry[iref];
          pos[2] = rz[iref];
          
          distsq = rx[iref]*rx[iref] + ry[iref]*ry[iref] +rz[iref]*rz[iref];
          dist = sqrt(distsq);
          
          switch(buildType){

          case BUILD_CORE_SHELL:
            nanoBuilder::buildWithCoreShell(dist,pos);
            break;
          case BUILD_CORE_SHELL_VACANCY:
            nanoBuilder::buildWithVacancies(dist,pos);
            break;

          case BUILD_RANDOM_PARTICLE:
            nanoBuilder::buildRandomlyMixed(dist,pos);
            break;
          case BUILD_NMOL_PARTICLE:
            nanoBuilder::buildNmolParticle(dist,pos);
          }
        }
      }
    }
  }


  // Create vacancies
  if (hasVacancies) buildVacancies();

  // Find the size of the atom vector not including Null atoms
  for (i=0;i<moleculeVector.size();i++){
    if (! moleculeVector[i].isVacancy){
      shesActualSizetoMe++;
      nAtoms = moleculeVector[i].myStamp->getNAtoms();
    }
  }

// Make a random particle.
  if (isRandom){
    placeRandom(shesActualSizetoMe); 

  // Loop back thru and count natoms since they may have changed
    for (i=0;i<moleculeVector.size();i++){
      if (! moleculeVector[i].isVacancy){
        shesActualSizetoMe++;
        nAtoms = moleculeVector[i].myStamp->getNAtoms();
      }
    }
  }


  Atom::createArrays( nAtoms );
  atoms = new Atom*[nAtoms];

 
  shesActualSizetoMe = 0;
  /*  Use the information from the molecule vector to place the atoms.
   */
  for (i= 0;i<moleculeVector.size();i++){
    if (! moleculeVector[i].isVacancy) {
      orientationMunger( A );
      if( moleculeVector[i].isCore){
        nCoreAtomCounter =+ nCoreModelAtoms;
        coreLocate->placeMol(moleculeVector[i].pos,A,atoms,nShellAtomCounter);
      }
      else {
        nShellAtomCounter =+ nShellModelAtoms;
        shellLocate->placeMol(moleculeVector[i].pos,A,atoms,nCoreAtomCounter);
      }
      shesActualSizetoMe++;
    }
  }


  //      shellLocate.placeMol(pos, A, moleculeVector,shellAtomCount);

  for (i=0;i<9;i++) simnfo->Hmat[i] = 0;
  simnfo->Hmat[0] = 1;
  simnfo->Hmat[4] = 1;
  simnfo->Hmat[8] = 1;

  // set up the SimInfo object
  simnfo = new SimInfo();
  simnfo->n_atoms = nAtoms;

  
  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );

  simnfo->atoms = atoms;
  
  // set up the writer and write out
  
  writer = new DumpWriter( simnfo );
  writer->writeFinal(0.0);

    // clean up  

  delete[] myLattice;

  return hasError;
} 

// Begin Builder routines------------------------------->

/* Builds a standard core-shell nanoparticle.
*/
void nanoBuilder::buildWithCoreShell(double dist, double pos[3]){

  
  if ( dist <= particleRadius ){
    moleculeVector.push_back(myMol);
    
    if (dist <= coreRadius){
      coreAtomCount =+ nCoreModelAtoms;
      moleculeVector[moleculeCount].pos[0] = pos[0]; 
      moleculeVector[moleculeCount].pos[1] = pos[1]; 
      moleculeVector[moleculeCount].pos[2] = pos[2]; 
      moleculeVector[moleculeCount].myStamp = coreStamp;
      moleculeVector[moleculeCount].isCore = 1;
      moleculeVector[moleculeCount].isShell = 0;
      
    }
    // Place shell
    else{
      shellAtomCount =+ nShellModelAtoms;
      moleculeVector[moleculeCount].pos[0] = pos[0]; 
      moleculeVector[moleculeCount].pos[1] = pos[1]; 
      moleculeVector[moleculeCount].pos[2] = pos[2]; 
      moleculeVector[moleculeCount].myStamp = shellStamp;
      moleculeVector[moleculeCount].isCore = 0;
      moleculeVector[moleculeCount].isShell = 1;
      
    }
    moleculeCount++;
  }      
  
}
/*
Builds a core-shell nanoparticle and tracks the number of molecules at the
interface between the core-shell. These are recorded in vacancyInterface which is just
an integer vector. 
*/
void nanoBuilder::buildWithVacancies(double dist, double pos[3]){
  if ( dist <= particleRadius ){

    moleculeVector.push_back(myMol);
    if (dist <= coreRadius){
        
      coreAtomCount =+ nCoreModelAtoms;
      moleculeVector[moleculeCount].pos[0] = pos[0]; 
      moleculeVector[moleculeCount].pos[1] = pos[1]; 
      moleculeVector[moleculeCount].pos[2] = pos[2]; 
      moleculeVector[moleculeCount].myStamp = coreStamp;
      moleculeVector[moleculeCount].isCore = 1;
      moleculeVector[moleculeCount].isShell = 0;

      if ((dist >= coreRadius - vacancyRadius/2.0) && 
          (dist <= coreRadius + vacancyRadius/2.0)){
        
        vacancyInterface.push_back(moleculeCount);
        nInterface++;
      }
    } else {
      // Place shell
      shellAtomCount =+ nShellModelAtoms;
      moleculeVector[moleculeCount].pos[0] = pos[0]; 
      moleculeVector[moleculeCount].pos[1] = pos[1]; 
      moleculeVector[moleculeCount].pos[2] = pos[2]; 
      moleculeVector[moleculeCount].myStamp = shellStamp;
      moleculeVector[moleculeCount].isCore = 0;
      moleculeVector[moleculeCount].isShell = 1;

    }
    moleculeCount++;
  }


}

/* Builds a core-shell nanoparticle where the number of core and shell
 molecules is known.
*/
void nanoBuilder::buildNmolParticle(double dist, double pos[3]){
  static int nMolCounter = 0;
  static int nCoreMolCounter = 0;
  
  
  if (nMolCounter < totalMolecules){
    moleculeVector.push_back(myMol);
    if (nCoreMolCounter < nCoreMolecules){
      
      coreAtomCount =+ nCoreModelAtoms;
      moleculeVector[moleculeCount].pos[0] = pos[0]; 
      moleculeVector[moleculeCount].pos[1] = pos[1]; 
      moleculeVector[moleculeCount].pos[2] = pos[2]; 
      moleculeVector[moleculeCount].myStamp = coreStamp;
      moleculeVector[moleculeCount].isCore = 1;
      moleculeVector[moleculeCount].isShell = 0;
      
      
    } else {
      shellAtomCount =+ nShellModelAtoms;
      moleculeVector[moleculeCount].pos[0] = pos[0]; 
      moleculeVector[moleculeCount].pos[1] = pos[1]; 
      moleculeVector[moleculeCount].pos[2] = pos[2]; 
      moleculeVector[moleculeCount].myStamp = shellStamp;
      moleculeVector[moleculeCount].isCore = 0;
      moleculeVector[moleculeCount].isShell = 1;
      
        
    }

  }
}


/* Builds a randomly mixed nanoparticle. We build the particle to be 
 entirely the core model, then randomly switch identities after the particle is built.
*/ 
void nanoBuilder::buildRandomlyMixed(double dist, double pos[3]){


  if ( dist <= particleRadius ){
    moleculeCount++;


    moleculeVector[moleculeCount].pos[0] = pos[0]; 
    moleculeVector[moleculeCount].pos[1] = pos[1]; 
    moleculeVector[moleculeCount].pos[2] = pos[2]; 
    moleculeVector[moleculeCount].myStamp = coreStamp;
    moleculeVector[moleculeCount].isCore = 1;
    moleculeVector[moleculeCount].isShell = 0;
    
  }     


}


// -----------------------END Builder routines.


//------------------------Begin Helper routines.
void nanoBuilder::placeRandom(int totalMol){
  int nSolute;
  int nSolvent;
  int i;
  int notfound;
  double solute_x;
  double solvent_x;
  
  int tester;

  nSolute = floor(soluteX * (double)totalMolecules); //CHECK ME
  nSolvent = totalMolecules - nSolute;
  
  solute_x = (double)nSolute/(double)totalMolecules;
  solvent_x = 1.0 - solute_x;
  
  
  for(i=0;nSolute-1;i++){
    notfound = 1;
    
    while(notfound){
      
      tester = floor((double)totalMolecules * drand48()); //Pick a molecule
      
      if (moleculeVector[tester].isCore){ // Make sure we select a core atom to change
          
        moleculeVector[tester].isCore  = 0;
        moleculeVector[tester].isShell = 1;
        moleculeVector[tester].myStamp = shellStamp;
        notfound = 0; //set notfound = false.
      }
        
    }
      
  }
}


void nanoBuilder::buildVacancies(void){
  int i;  
  int* VacancyList; //logical nInterface long.
  int notfound;
  int index = 0;
  int nVacancies;
  int tester;

  if (nInterface != 0){
    nVacancies = floor((double)nInterface * vacancyFraction);

    VacancyList = new int[nInterface];
    
    // make vacancy list all false
    for(i=0;i<nInterface-1;i++){
      VacancyList[i] = 0;
    }
    
    // Build a vacancy list....    
    for(i=0;nVacancies-1;i++){
      notfound = 1;
      while(notfound){
        
        tester = floor((double)nInterface * drand48());
      
        if(! VacancyList[tester]){
          VacancyList[tester] = 1;
          notfound = 0;
        }
        
      }
    }
  }
  // Loop through and kill the vacancies from atom vector.

  for (i=0;i<nInterface;i++){
    if (VacancyList[i]){
      moleculeVector[vacancyInterface[i]].isVacancy = 1;   
    } // End Vacancy List
  } // for nInterface
  

    delete[] VacancyList;
}


void nanoBuilder::orientationMunger(double rot[3][3]){

  double theta, phi, psi;
  double cosTheta;

  // select random phi, psi, and cosTheta

  phi = 2.0 * M_PI * drand48();
  psi = 2.0 * M_PI * drand48();
  cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1

  theta = acos( cosTheta );

  rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
  rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
  rot[0][2] = sin(theta) * sin(psi);

  rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
  rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
  rot[1][2] = sin(theta) * cos(psi);

  rot[2][0] = sin(phi) * sin(theta);
  rot[2][1] = -cos(phi) * sin(theta);
  rot[2][2] = cos(theta);

}


Revision:	589
Committed:	Thu Jul 10 19:53:50 2003 UTC (21 years ago) by chuckv
File size:	15076 byte(s)
Log Message:	Added nanoBuilder and a general Lattice builder.
#	User	Rev	Content
1	chuckv	589	#include <iostream>
2
3			#include <cstdlib>
4			#include <cstring>
5			#include <cmath>
6			#include <vector>
7
8			#include "simError.h"
9			#include "SimInfo.hpp"
10			#include "ReadWrite.hpp"
11
12			#include "latticeBuilder.hpp"
13			#include "MoLocator.hpp"
14			#include "sysBuild.hpp"
15			#include "nanoBuilder.hpp"
16
17
18
19			nanoBuilder::nanoBuilder(int thisIsRandom, int thisHasVacancies,
20			int thisLatticeType, double thisParticleRadius,
21			double thisCoreRadius, double thisVacancyFraction,
22			double thisVacancyRadius,
23			double thisLatticeSpacing,
24			double solute_X,
25			int &hasError){
26			int Errors;
27			int foundCore,foundShell;
28			int i;
29
30			//Zero variables
31			particleRadius = 0.0;
32			coreRadius = 0.0;
33			vacancyFraction = 0.0;
34			vacancyRadius = 0.0;
35			shellRadius = 0.0;
36			latticeSpacing = 0.0;
37
38			buildNmol = 0;
39
40			nCoreMolecules = 0;
41			nShellMolecules = 0;
42
43			atomCount = 0;
44			coreAtomCount = 0;
45			shellAtomCount = 0;
46
47
48
49			moleculeCount = 0;
50			foundCore = 0;
51			foundShell = 0;
52			totalMolecules = 0;
53			coreHasOrientation = 0;
54			shellHasOrientation = 0;
55			nInterface = 0;
56			nMol = 0;
57
58			hasError = 0;
59			Errors = 0;
60
61
62			isRandom = thisIsRandom;
63			hasVacancies = thisHasVacancies;
64			latticeType = thisLatticeType;
65			particleRadius = thisParticleRadius;
66			coreRadius = thisCoreRadius;
67			vacancyFraction = thisVacancyFraction;
68			latticeSpacing = thisLatticeSpacing;
69			soluteX = solute_X; //Mole fraction for random particle.
70
71
72
73
74
75			for (i=0;bsInfo.nComponents;i++){
76			if( !strcmp( bsInfo.compStamps[i]->getID(),bsInfo.coreName )){
77			foundCore = 1;
78			coreStamp = bsInfo.compStamps[i];
79			nCoreMolecules = bsInfo.componentsNmol[i];
80			}
81			if( !strcmp( bsInfo.compStamps[i]->getID(),bsInfo.shellName)){
82			foundShell = 1;
83			shellStamp = bsInfo.compStamps[i];
84			nShellMolecules = bsInfo.componentsNmol[i];
85
86			}
87
88			}
89
90
91
92			if( !foundCore ){
93			hasError = 1;
94			return;
95			}
96			if( !foundShell ){
97			hasError = 1;
98			return;
99			}
100
101
102
103			Errors = sanityCheck();
104
105			if (Errors){
106			hasError = 1;
107			return;
108			}
109
110
111
112
113
114
115
116			nCoreModelAtoms = coreStamp->getNAtoms();
117			nShellModelAtoms = shellStamp->getNAtoms();
118
119
120			// We assume that if the core or shell model has more then one atom
121			// the model has an orientational component...
122			if (nCoreModelAtoms > 1) coreHasOrientation = 1;
123			if (nShellModelAtoms > 1) shellHasOrientation = 1;
124
125			maxModelNatoms = std::max(nCoreModelAtoms,nShellModelAtoms);
126
127			/* If we specify a number of atoms in bass, we will try to build a nanopartice
128			with that number.
129			*/
130
131
132			if ((nShellMolecules != 0) && (nCoreMolecules != 0)){
133			totalMolecules = nShellMolecules + nCoreMolecules;
134			nCells = ceil(pow((double)totalMolecules/4.0, 1/3));
135			buildNmol = 1;
136			}
137			else {
138			nCells = 2.0 * particleRadius/latticeSpacing;
139			shellRadius = particleRadius - coreRadius;
140			}
141
142
143
144
145			// Initialize random seed
146			srand48( RAND_SEED );
147
148
149			}
150
151
152			nanoBuilder::~nanoBuilder(){
153			}
154
155
156			// Checks to make sure we aren't doing something the builder can't do.
157			int nanoBuilder::sanityCheck(void){
158
159			// Right now we only do bimetallic nanoparticles
160			if (bsInfo.nComponents > 2) return 1;
161
162			//Check for vacancies and random
163			if (hasVacancies && isRandom) return 1;
164
165			// make sure we aren't trying to build a core larger then the total particle size
166			if ((coreRadius >= particleRadius) && (particleRadius != 0)) return 1;
167
168			// we initialize the lattice spacing to be 0.0, if the lattice spacing is still 0.0
169			// we have a problem
170			if (latticeSpacing == 0.0) return 1;
171
172			// Check to see if we are specifing the number of atoms in the particle correctly.
173			if ((nShellMolecules == 0) && (nCoreMolecules != 0)){
174			cerr << "nShellParticles is zero and nCoreParticles != 0" << "\n";
175			return 1;
176			}
177			// Make sure there are more then two components if we are building a randomly mixed particle.
178			if ((bsInfo.nComponents < 2) && (isRandom)){
179			cerr << "Two Components are needed to build a random particle." << "\n";
180			}
181			// Make sure both the core and shell models specify a target nmol.
182			if ((nShellMolecules != 0) && (nCoreMolecules == 0)){
183			cerr << "nCoreParticles is zero and nShellParticles != 0" << "\n";
184			return 1;
185			}
186
187			return 0;
188
189			}
190
191
192
193			int nanoBuilder::buildNanoParticle( void ){
194
195			int ix;
196			int iy;
197			int iz;
198			double *rx;
199			double *ry;
200			double *rz;
201			double pos[3];
202			double A[3][3];
203			double HmatI[9];
204
205			int nCellSites;
206			int iref;
207			int appNMols;
208			int latticeCount = 0;
209
210			int nAtoms;
211			int nCoreAtomCounter = 0;
212			int nShellAtomCounter = 0;
213			int hasError;
214
215			int i;
216
217			int interfaceIndex = 0;
218			double dist;
219			double distsq;
220			int latticeNpoints;
221			int shesActualSizetoMe = 0;
222
223			DumpWriter* writer;
224			SimInfo* simnfo;
225
226			Lattice::Lattice *myLattice;
227			MoLocator::MoLocator *coreLocate;
228			MoLocator::MoLocator *shellLocate;
229
230
231			Atom** atoms;
232
233			hasError = 0;
234
235			myLattice = new Lattice(FCC_LATTICE_TYPE,latticeSpacing);
236			/*
237			latticeNpoints = myLattice.getNpoints();
238
239			// Initializd atom vector to approximate size.
240			switch (buildType){
241
242			case BUILD_NMOL_PARTICLE:
243
244			break;
245			case BUILD_CORE_SHELL_VACANCY:
246			// Make space in the vector for all atoms except the last full cells
247			// We will have to add at most (latticeNpoints-1)^3 to vector
248			appNMols = latticeNPoints * pow((double)(nCells - 1),3);
249			moleculeVector.pushBack();
250
251			default:
252			// Make space in the vector for all atoms except the last full cells
253			// We will have to add at most (latticeNpoints-1)^3 to vector
254			appNMols = latticeNPoints * pow((double)(nCells - 1),3);
255
256			}
257			*/
258
259
260
261
262			// Create molocator and atom arrays.
263			coreLocate = new MoLocator(coreStamp);
264			shellLocate = new MoLocator(shellStamp);
265
266
267
268
269
270
271			for(iz=-nCells;iz < nCells;iz++){
272			for(iy=-nCells;iy<nCells;iy++){
273			for(ix=-nCells;ix<nCells;ix++){
274			nCellSites = myLattice->getLatticePoints(&rx,&ry,&rz,
275			ix,iy,iz);
276			for (iref=1;iref<nCellSites;iref++){
277			latticeCount++;
278
279			pos[0] = rx[iref];
280			pos[1] = ry[iref];
281			pos[2] = rz[iref];
282
283			distsq = rx[iref]rx[iref] + ry[iref]ry[iref] +rz[iref]*rz[iref];
284			dist = sqrt(distsq);
285
286			switch(buildType){
287
288			case BUILD_CORE_SHELL:
289			nanoBuilder::buildWithCoreShell(dist,pos);
290			break;
291			case BUILD_CORE_SHELL_VACANCY:
292			nanoBuilder::buildWithVacancies(dist,pos);
293			break;
294
295			case BUILD_RANDOM_PARTICLE:
296			nanoBuilder::buildRandomlyMixed(dist,pos);
297			break;
298			case BUILD_NMOL_PARTICLE:
299			nanoBuilder::buildNmolParticle(dist,pos);
300			}
301			}
302			}
303			}
304			}
305
306
307
308			// Create vacancies
309			if (hasVacancies) buildVacancies();
310
311			// Find the size of the atom vector not including Null atoms
312			for (i=0;i<moleculeVector.size();i++){
313			if (! moleculeVector[i].isVacancy){
314			shesActualSizetoMe++;
315			nAtoms = moleculeVector[i].myStamp->getNAtoms();
316			}
317			}
318
319			// Make a random particle.
320			if (isRandom){
321			placeRandom(shesActualSizetoMe);
322
323			// Loop back thru and count natoms since they may have changed
324			for (i=0;i<moleculeVector.size();i++){
325			if (! moleculeVector[i].isVacancy){
326			shesActualSizetoMe++;
327			nAtoms = moleculeVector[i].myStamp->getNAtoms();
328			}
329			}
330			}
331
332
333			Atom::createArrays( nAtoms );
334			atoms = new Atom*[nAtoms];
335
336
337
338			shesActualSizetoMe = 0;
339			/* Use the information from the molecule vector to place the atoms.
340			*/
341			for (i= 0;i<moleculeVector.size();i++){
342			if (! moleculeVector[i].isVacancy) {
343			orientationMunger( A );
344			if( moleculeVector[i].isCore){
345			nCoreAtomCounter =+ nCoreModelAtoms;
346			coreLocate->placeMol(moleculeVector[i].pos,A,atoms,nShellAtomCounter);
347			}
348			else {
349			nShellAtomCounter =+ nShellModelAtoms;
350			shellLocate->placeMol(moleculeVector[i].pos,A,atoms,nCoreAtomCounter);
351			}
352			shesActualSizetoMe++;
353			}
354			}
355
356
357			// shellLocate.placeMol(pos, A, moleculeVector,shellAtomCount);
358
359			for (i=0;i<9;i++) simnfo->Hmat[i] = 0;
360			simnfo->Hmat[0] = 1;
361			simnfo->Hmat[4] = 1;
362			simnfo->Hmat[8] = 1;
363
364			// set up the SimInfo object
365			simnfo = new SimInfo();
366			simnfo->n_atoms = nAtoms;
367
368
369			sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
370			sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
371
372			simnfo->atoms = atoms;
373
374			// set up the writer and write out
375
376			writer = new DumpWriter( simnfo );
377			writer->writeFinal(0.0);
378
379			// clean up
380
381			delete[] myLattice;
382
383			return hasError;
384			}
385
386			// Begin Builder routines------------------------------->
387
388			/* Builds a standard core-shell nanoparticle.
389			*/
390			void nanoBuilder::buildWithCoreShell(double dist, double pos[3]){
391
392
393			if ( dist <= particleRadius ){
394			moleculeVector.push_back(myMol);
395
396			if (dist <= coreRadius){
397			coreAtomCount =+ nCoreModelAtoms;
398			moleculeVector[moleculeCount].pos[0] = pos[0];
399			moleculeVector[moleculeCount].pos[1] = pos[1];
400			moleculeVector[moleculeCount].pos[2] = pos[2];
401			moleculeVector[moleculeCount].myStamp = coreStamp;
402			moleculeVector[moleculeCount].isCore = 1;
403			moleculeVector[moleculeCount].isShell = 0;
404
405			}
406			// Place shell
407			else{
408			shellAtomCount =+ nShellModelAtoms;
409			moleculeVector[moleculeCount].pos[0] = pos[0];
410			moleculeVector[moleculeCount].pos[1] = pos[1];
411			moleculeVector[moleculeCount].pos[2] = pos[2];
412			moleculeVector[moleculeCount].myStamp = shellStamp;
413			moleculeVector[moleculeCount].isCore = 0;
414			moleculeVector[moleculeCount].isShell = 1;
415
416			}
417			moleculeCount++;
418			}
419
420			}
421			/*
422			Builds a core-shell nanoparticle and tracks the number of molecules at the
423			interface between the core-shell. These are recorded in vacancyInterface which is just
424			an integer vector.
425			*/
426			void nanoBuilder::buildWithVacancies(double dist, double pos[3]){
427			if ( dist <= particleRadius ){
428
429			moleculeVector.push_back(myMol);
430			if (dist <= coreRadius){
431
432			coreAtomCount =+ nCoreModelAtoms;
433			moleculeVector[moleculeCount].pos[0] = pos[0];
434			moleculeVector[moleculeCount].pos[1] = pos[1];
435			moleculeVector[moleculeCount].pos[2] = pos[2];
436			moleculeVector[moleculeCount].myStamp = coreStamp;
437			moleculeVector[moleculeCount].isCore = 1;
438			moleculeVector[moleculeCount].isShell = 0;
439
440			if ((dist >= coreRadius - vacancyRadius/2.0) &&
441			(dist <= coreRadius + vacancyRadius/2.0)){
442
443			vacancyInterface.push_back(moleculeCount);
444			nInterface++;
445			}
446			} else {
447			// Place shell
448			shellAtomCount =+ nShellModelAtoms;
449			moleculeVector[moleculeCount].pos[0] = pos[0];
450			moleculeVector[moleculeCount].pos[1] = pos[1];
451			moleculeVector[moleculeCount].pos[2] = pos[2];
452			moleculeVector[moleculeCount].myStamp = shellStamp;
453			moleculeVector[moleculeCount].isCore = 0;
454			moleculeVector[moleculeCount].isShell = 1;
455
456			}
457			moleculeCount++;
458			}
459
460
461
462			}
463
464			/* Builds a core-shell nanoparticle where the number of core and shell
465			molecules is known.
466			*/
467			void nanoBuilder::buildNmolParticle(double dist, double pos[3]){
468			static int nMolCounter = 0;
469			static int nCoreMolCounter = 0;
470
471
472			if (nMolCounter < totalMolecules){
473			moleculeVector.push_back(myMol);
474			if (nCoreMolCounter < nCoreMolecules){
475
476			coreAtomCount =+ nCoreModelAtoms;
477			moleculeVector[moleculeCount].pos[0] = pos[0];
478			moleculeVector[moleculeCount].pos[1] = pos[1];
479			moleculeVector[moleculeCount].pos[2] = pos[2];
480			moleculeVector[moleculeCount].myStamp = coreStamp;
481			moleculeVector[moleculeCount].isCore = 1;
482			moleculeVector[moleculeCount].isShell = 0;
483
484
485			} else {
486			shellAtomCount =+ nShellModelAtoms;
487			moleculeVector[moleculeCount].pos[0] = pos[0];
488			moleculeVector[moleculeCount].pos[1] = pos[1];
489			moleculeVector[moleculeCount].pos[2] = pos[2];
490			moleculeVector[moleculeCount].myStamp = shellStamp;
491			moleculeVector[moleculeCount].isCore = 0;
492			moleculeVector[moleculeCount].isShell = 1;
493
494
495			}
496
497			}
498			}
499
500
501			/* Builds a randomly mixed nanoparticle. We build the particle to be
502			entirely the core model, then randomly switch identities after the particle is built.
503			*/
504			void nanoBuilder::buildRandomlyMixed(double dist, double pos[3]){
505
506
507			if ( dist <= particleRadius ){
508			moleculeCount++;
509
510
511			moleculeVector[moleculeCount].pos[0] = pos[0];
512			moleculeVector[moleculeCount].pos[1] = pos[1];
513			moleculeVector[moleculeCount].pos[2] = pos[2];
514			moleculeVector[moleculeCount].myStamp = coreStamp;
515			moleculeVector[moleculeCount].isCore = 1;
516			moleculeVector[moleculeCount].isShell = 0;
517
518			}
519
520
521
522			}
523
524
525			// -----------------------END Builder routines.
526
527
528
529			//------------------------Begin Helper routines.
530			void nanoBuilder::placeRandom(int totalMol){
531			int nSolute;
532			int nSolvent;
533			int i;
534			int notfound;
535			double solute_x;
536			double solvent_x;
537
538			int tester;
539
540			nSolute = floor(soluteX * (double)totalMolecules); //CHECK ME
541			nSolvent = totalMolecules - nSolute;
542
543			solute_x = (double)nSolute/(double)totalMolecules;
544			solvent_x = 1.0 - solute_x;
545
546
547
548
549			for(i=0;nSolute-1;i++){
550			notfound = 1;
551
552			while(notfound){
553
554			tester = floor((double)totalMolecules * drand48()); //Pick a molecule
555
556			if (moleculeVector[tester].isCore){ // Make sure we select a core atom to change
557
558			moleculeVector[tester].isCore = 0;
559			moleculeVector[tester].isShell = 1;
560			moleculeVector[tester].myStamp = shellStamp;
561			notfound = 0; //set notfound = false.
562			}
563
564			}
565
566			}
567			}
568
569
570			void nanoBuilder::buildVacancies(void){
571			int i;
572			int* VacancyList; //logical nInterface long.
573			int notfound;
574			int index = 0;
575			int nVacancies;
576			int tester;
577
578			if (nInterface != 0){
579			nVacancies = floor((double)nInterface * vacancyFraction);
580
581			VacancyList = new int[nInterface];
582
583			// make vacancy list all false
584			for(i=0;i<nInterface-1;i++){
585			VacancyList[i] = 0;
586			}
587
588			// Build a vacancy list....
589			for(i=0;nVacancies-1;i++){
590			notfound = 1;
591			while(notfound){
592
593			tester = floor((double)nInterface * drand48());
594
595			if(! VacancyList[tester]){
596			VacancyList[tester] = 1;
597			notfound = 0;
598			}
599
600			}
601			}
602			}
603			// Loop through and kill the vacancies from atom vector.
604
605			for (i=0;i<nInterface;i++){
606			if (VacancyList[i]){
607			moleculeVector[vacancyInterface[i]].isVacancy = 1;
608			} // End Vacancy List
609			} // for nInterface
610
611
612			delete[] VacancyList;
613			}
614
615
616
617
618			void nanoBuilder::orientationMunger(double rot[3][3]){
619
620			double theta, phi, psi;
621			double cosTheta;
622
623			// select random phi, psi, and cosTheta
624
625			phi = 2.0 * M_PI * drand48();
626			psi = 2.0 * M_PI * drand48();
627			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
628
629			theta = acos( cosTheta );
630
631			rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
632			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
633			rot[0][2] = sin(theta) * sin(psi);
634
635			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
636			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
637			rot[1][2] = sin(theta) * cos(psi);
638
639			rot[2][0] = sin(phi) * sin(theta);
640			rot[2][1] = -cos(phi) * sin(theta);
641			rot[2][2] = cos(theta);
642
643			}
644
645
646
647
648
649
650