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[3][3];
  
  int nCellSites;
  int iref;
  int appNMols;
  int latticeCount = 0;
 
  int nAtoms;
  int nCoreAtomCounter = 0;
  int nShellAtomCounter = 0;
  int hasError;

  int i, j;

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

  DumpWriter* writer;
  SimInfo* simnfo;

  Lattice *myLattice;
  MoLocator *coreLocate;
  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<3;i++) 
    for (j=0; j<3; j++) 
      simnfo->Hmat[i][j] = 0.0;

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