src/applications/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;
  SimState* theConfig;

  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();
      }
    }
  }


  // set up the SimInfo object

  simnfo = new SimInfo();
  simnfo->n_atoms = nAtoms;

  theConfig = simnfo->getConfiguration();
  theConfig->createArrays( nAtoms );
  simnfo->atoms = new Atom*[nAtoms];
  atoms = simnfo->atoms;
 

  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, theConfig);
      }
      else {
        nShellAtomCounter += nShellModelAtoms;
        shellLocate->placeMol(moleculeVector[i].pos,A,atoms,nCoreAtomCounter, theConfig);
      }
      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;
  

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

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