utils/sysbuilder/nanoBuilder.cpp

#include <iostream>

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


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

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

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