ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/nanorodBuilder/nanoBuilder.cpp
Revision: 2148
Committed: Wed Apr 6 16:34:09 2005 UTC (19 years, 2 months ago) by chuckv
File size: 15071 byte(s)
Log Message:
Initial import of nanorod builder.

File Contents

# Content
1 #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