ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/utils/bilayerSys.cpp
(Generate patch)

Comparing trunk/OOPSE/utils/bilayerSys.cpp (file contents):
Revision 498 by mmeineke, Mon Apr 14 21:22:54 2003 UTC vs.
Revision 537 by mmeineke, Fri May 16 21:37:56 2003 UTC

# Line 1 | Line 1
1 + #include <iostream>
2 +
3   #include <cstdlib>
4   #include <cstring>
5   #include <cmath>
6  
7   #include "simError.h"
6 – #include "parse_me.h"
7 – #include "MakeStamps.hpp"
8 – #include "Globals.hpp"
8   #include "SimInfo.hpp"
9   #include "ReadWrite.hpp"
10  
11 <
11 > #include "MoLocator.hpp"
12   #include "sysBuild.hpp"
13   #include "bilayerSys.hpp"
14  
15  
17 – // this routine is defined in BASS_interface.cpp
18 – extern void set_interface_stamps( MakeStamps* ms, Globals* g );
16  
17 < void buildRandomBilayer( sysBuildInfo info );
17 > void map( double &x, double &y, double &z,
18 >          double boxX, double boxY, double boxZ );
19  
20 + int buildRandomBilayer( void );
21  
22 < void buildBilayer( sysBuildInfo info ){
22 > void getRandomRot( double rot[3][3] );
23 >
24 > int buildBilayer( int isRandom ){
25    
26 <  if( info.isRandom ){
27 <    buildRandomBilayer( info );
26 >  if( isRandom ){
27 >    return buildRandomBilayer();
28    }
29    else{
29 –
30      sprintf( painCave.errMsg,
31               "Cannot currently create a non-random bilayer.\n" );
32      painCave.isFatal = 1;
33      simError();
34 +    return 0;
35    }
36   }
37  
38  
39 + int buildRandomBilayer( void ){
40  
41 < void buildRandomBilayer( sysBuildInfo info ){
41 >  typedef struct{
42 >    double rot[3][3];
43 >    double pos[3];
44 >  } coord;
45  
46 <  MakeStamps* the_stamps;
47 <  Globals* the_globals;
48 <  SimInfo* simnfo;
49 <  bwMolLinked bwInfo;
46 >
47 >  const double waterRho = 0.0334; // number density per cubic angstrom
48 >  const double waterVol = 4.0 / water_rho; // volume occupied by 4 waters
49 >  const double waterCell = 4.929; // fcc unit cell length
50 >
51 >  const double water_padding = 2.5;
52 >  const double lipid_spaceing = 5.0;
53 >
54 >
55 >  int i,j,k;
56 >  int nAtoms, atomIndex, molIndex, molID;
57 >  int* molSeq;
58 >  int* molMap;
59 >  int* molStart;
60 >  int* cardDeck;
61 >  int deckSize;
62 >  int rSite, rCard;
63 >  double cell;
64 >  int nCells, nSites, siteIndex;
65    
66 <  // init the bwInfo
66 >  coord testSite;
67 >
68 >  MoleculeStamp* lipidStamp;
69 >  MoleculeStamp* waterStamp;  
70 >  MoLocator *lipidLocate;
71 >  MoLocator *waterLocate
72 >  int foundLipid, foundWater;
73 >  int nLipids, lipiNatoms, nWaters, waterNatoms;
74 >  double testBox, maxLength;
75    
76 <  bwinfo.components = NULL;
76 >  srand48( RAND_SEED );
77 >
78 >
79 >  // set the the lipidStamp
80 >
81 >  foundLipid = 0;
82 >  for(i=0; i<bsInfo.nComponents; i++){
83 >    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
84 >      
85 >      foundlipid = 1;
86 >      lipidStamp = bsInfo.compStamps[i];
87 >      nLipids = bsInfo.componentsNmol[i];
88 >    }
89 >    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
90 >      
91 >      foundWater = 1;
92 >      waterStamp = bsInfo.compStamps[i];
93 >      nWaters = bsInfo.componentsNmol[i];
94 >    }
95 >  }
96 >  if( !foundLipid ){
97 >    sprintf(painCave.errMsg,
98 >            "Could not find lipid \"%s\" in the bass file.\n",
99 >            bsInfo.lipidName );
100 >    painCave.isFatal = 1;
101 >    simError();
102 >  }
103 >  if( !foundWater ){
104 >    sprintf(painCave.errMsg,
105 >            "Could not find solvent \"%s\" in the bass file.\n",
106 >            bsInfo.waterName );
107 >    painCave.isFatal = 1;
108 >    simError();
109 >  }
110    
111 <  bwInfo.havePressure = 0;
51 <  bwInfo.haveTauBarrostat = 0;
52 <  bwInfo.haveTauTemp = 0;
53 <  bwInfo.haveQmass = 0;
111 >  //create the temp Molocator and atom Arrays
112    
113 +  lipidLocate = new MoLocator( lipidStamp );
114 +  lipidNatoms = lipidStamp->getNAtoms();
115 +  maxLength = lipidLocate->getMaxLength();
116 +
117 +  waterLocate = new MoLocator( waterStamp );
118 +  waterNatoms = waterStamp->getNatoms();
119    
120 <  // create parser and read the Bas file
120 >  nAtoms = nLipids * lipidNatoms;
121 >
122 >  Atom::createArrays( nAtoms );
123 >  atoms = new Atom*[nAtoms];
124 >
125 >  // create the test box for initial water displacement
126 >
127 >  testBox = maxLength + waterCell * 4.0; // pad with 4 cells
128 >  int nCells = (int)( testBox / waterCell + 1.0 );
129 >  int testWaters = 4 * nCells * nCells * nCells;
130    
131 <  simnfo = new SimInfo();
132 <  the_stamps = new MakeStamps();
133 <  the_globals = new Globals();
134 <  set_interface_stamps( the_stamps, the_globals );
131 >  double* waterX = new double[testWaters];
132 >  double* waterX = new double[testWaters];
133 >  double* waterX = new double[testWaters];
134 >  
135 >  double x0 = 0.0 - ( testBox * 0.5 );
136 >  double y0 = 0.0 - ( testBox * 0.5 );
137 >  double z0 = 0.0 - ( testBox * 0.5 );
138  
63 –  yacc_BASS( info.in_name );  
139  
140 +  // create an fcc lattice in the water box.
141  
142 +  int ndx = 0;
143 +  for( i=0; i < nCells; i++ ){
144 +    for( j=0; j < nCells; j++ ){
145 +      for( k=0; k < nCells; k++ ){
146 +        
147 +        waterX[ndx] = i * waterCell + x0;
148 +        waterY[ndx] = j * waterCell + y0;
149 +        waterZ[ndx] = k * waterCell + z0;
150 +        ndx++;
151 +        
152 +        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
153 +        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
154 +        waterZ[ndx] = k * waterCell + z0;
155 +        ndx++;
156 +        
157 +        waterX[ndx] = i * waterCell + x0;
158 +        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
159 +        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
160 +        ndx++;
161 +        
162 +        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
163 +        waterY[ndx] = j * waterCell + y0;
164 +        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
165 +        ndx++;
166 +      }
167 +    }
168 +  }
169  
170 <  // set the easy ones first
68 <  bwInfo.targetTemp = the_globals->getTargetTemp();
69 <  bwInfo.dt = the_globals->getDt();
70 <  bwInfo.runTime = the_globals->getRunTime();
170 >  // calculate the number of water's displaced by our lipid.
171  
172 <  // get the ones we know are there, yet still may need some work.
173 <  bwInfo.nComponents = the_globals->getNComponents();
174 <  strcpy( bwInfo.forceField, the_globals->getForceField() );
172 >  testSite.rot[0][0] = 1.0;
173 >  testSite.rot[0][1] = 0.0;
174 >  testSite.rot[0][2] = 0.0;
175  
176 <  // get the ensemble:
177 <  strcpy( bwInfo.ensemble, the_globals->getEnsemble() );
178 <  if( !strcasecmp( bwInfo.ensemble, "NPT" ) ) {
179 <    
180 <    if (the_globals->haveTargetPressure()){
181 <      bwInfo.targetPressure = the_globals->getTargetPressure();
182 <      bwInfo.havePressure = 1;
176 >  testSite.rot[1][0] = 0.0;
177 >  testSite.rot[1][1] = 1.0;
178 >  testSite.rot[1][2] = 0.0;
179 >
180 >  testSite.rot[2][0] = 0.0;
181 >  testSite.rot[2][1] = 0.0;
182 >  testSite.rot[2][2] = 1.0;
183 >  
184 >  testSite.pos[0] = 0.0;
185 >  testSite.pos[1] = 0.0;
186 >  testSite.pos[2] = 0.0;
187 >
188 >  lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0 );
189 >
190 >  int *isActive = new int[testWaters];
191 >  for(i=0; i<testWaters; i++) isActive[i] = 1;
192 >  
193 >  int n_deleted = 0;
194 >  double dx, dy, dz;
195 >  double dx2, dy2, dz2, dSqr;
196 >  double rCutSqr = water_padding * water_padding;
197 >  
198 >  for(i=0; ( (i<testWaters) && isActive[i] ); i++){
199 >    for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
200 >
201 >      dx = waterX[i] - atoms[j]->getX();
202 >      dy = waterY[i] - atoms[j]->getY();
203 >      dz = waterZ[i] - atoms[j]->getZ();
204 >
205 >      map( dx, dy, dz, testBox, testBox, testBox );
206 >
207 >      dx2 = dx * dx;
208 >      dy2 = dy * dy;
209 >      dz2 = dz * dz;
210 >      
211 >      dSqr = dx2 + dy2 + dz2;
212 >      if( dSqr < rCutSqr ){
213 >        isActive[i] = 0;
214 >        n_deleted++;
215 >      }
216      }
217 <    else {
218 <      sprintf( painCave.errMsg,
219 <               "buildBilayer error: If you use the constant pressure\n"
87 <               "    ensemble, you must set targetPressure.\n"
88 <               "    This was found in the BASS file.\n");
89 <      painCave.isFatal = 1;
90 <      simError();
91 <    }
217 >  }
218 >  
219 >  int targetWaters = nWaters + n_deleted * nLipids;
220  
221 <    if (the_globals->haveTauThermostat()){
222 <      bwInfo.tauThermostat = the_globals->getTauThermostat();
223 <      bwInfo.haveTauThermostat = 1;;
221 >  // find the best box size for the sim
222 >
223 >  int testTot;
224 >  int done = 0;
225 >  ndx = 0;
226 >  while( !done ){
227 >    
228 >    ndx++;
229 >    testTot = 4 * ndx * ndx * ndx;
230 >    
231 >    if( testTot >= targetWaters ) done = 1;
232 >  }
233 >
234 >  nCells = ndx;
235 >    
236 >    
237 >  // create the new water box to the new specifications
238 >  
239 >  int newWaters = nCells * nCells * nCells * 4;
240 >  
241 >  delete[] waterX;
242 >  delete[] waterY;
243 >  delete[] waterZ;
244 >  
245 >  coord* waterSites = new coord[newWaters];
246 >  
247 >  double box_x = waterCell * nCells;
248 >  double box_y = waterCell * nCells;
249 >  double box_z = waterCell * nCells;
250 >    
251 >  // create an fcc lattice in the water box.
252 >  
253 >  ndx = 0;
254 >  for( i=0; i < nCells; i++ ){
255 >    for( j=0; j < nCells; j++ ){
256 >      for( k=0; k < nCells; k++ ){
257 >        
258 >        waterSites[ndx].pos[0] = i * waterCell;
259 >        waterSites[ndx].pos[1] = j * waterCell;
260 >        waterSites[ndx].pos[2] = k * waterCell;
261 >        ndx++;
262 >        
263 >        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
264 >        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
265 >        waterSites[ndx].pos[2] = k * waterCell;
266 >        ndx++;
267 >        
268 >        waterSites[ndx].pos[0] = i * waterCell;
269 >        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
270 >        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
271 >        ndx++;
272 >        
273 >        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
274 >        waterSites[ndx].pos[1] = j * waterCell;
275 >        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
276 >        ndx++;
277 >      }
278      }
279 <    else if (the_globals->haveQmass()){
98 <      bwinfo.Qmass = the_globals->getQmass();
99 <      bwInfo.haveQmass = 1;
100 <    }
101 <    else {
102 <      sprintf( painCave.errMsg,
103 <               "buildBilayer error: If you use one of the constant temperature\n"
104 <               "    ensembles, you must set either tauThermostat or qMass.\n"
105 <               "    Neither of these was found in the BASS file.\n");
106 <      painCave.isFatal = 1;
107 <      simError();
108 <    }
279 >  }  
280  
281 <    if (the_globals->haveTauBarostat()){
282 <      bwInfo.tauBarostat = the_globals->getTauBarostat();
283 <      bwInfo.haveTauBarostat = 1;
284 <    }                                    
285 <    else {
286 <      sprintf( painCave.errMsg,
287 <               "SimSetup error: If you use the constant pressure\n"
288 <               "    ensemble, you must set tauBarostat.\n"
289 <               "    This was found in the BASS file.\n");
290 <      painCave.isFatal = 1;
291 <      simError();
281 >
282 >  // clear up memory from the test box
283 >
284 >  for(i=0; i<lipidNatoms; i++ ) delete atoms[i];
285 >
286 >  coord* lipidSites = new coord[nLipids];
287 >
288 >  // start a 3D RSA for the for the lipid placements
289 >  
290 >  
291 >  int reject;
292 >  int testDX, acceptedDX;
293 >
294 >  rCutSqr = lipid_spaceing * lipid_spaceing;
295 >
296 >  for(i=0; i<nLipids; i++ ){
297 >    done = 0;
298 >    while( !done ){
299 >      
300 >      lipidSite[i].pos[0] = drand48() * box_x;
301 >      lipidSite[i].pos[1] = drand48() * box_y;
302 >      lipidSite[i].pos[2] = drand48() * box_z;
303 >      
304 >      getRandomRot( lipidSite[i].rot );
305 >      
306 >      ndx = i * lipidNatoms;
307 >
308 >      lipidLocate->placeMol( lipidSite[i].pos, lipidSite[i].rot, atoms, ndx );
309 >      
310 >      reject = 0;
311 >      for( j=0; !reject && j<i; j++){
312 >        for(k=0; !reject && k<lipidNatoms; k++){
313 >          
314 >          acceptedDX = j*lipidNatoms + k;
315 >          for(l=0; !reject && l<lipidNatoms; l++){
316 >            
317 >            testDX = ndx + l;
318 >
319 >            dx = atoms[testDX]->getX() - atoms[acceptedDX]->getX();
320 >            dy = atoms[testDX]->getY() - atoms[acceptedDX]->getY();
321 >            dz = atoms[testDX]->getZ() - atoms[acceptedDX]->getZ();
322 >            
323 >            map( dx, dy, dz, box_x, box_y, box_z );
324 >            
325 >            dx2 = dx * dx;
326 >            dy2 = dy * dy;
327 >            dz2 = dz * dz;
328 >            
329 >            dSqr = dx2 + dy2 + dz2;
330 >            if( dSqr < rCutSqr ) reject = 1;
331 >          }
332 >        }
333 >      }
334 >
335 >      if( reject ){
336 >
337 >        for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
338 >      }
339 >      else{
340 >        done = 1;
341 >        std::cout << i << " has been accepted\n";
342 >      }
343      }
344 +  }
345 +        
346 +  // cut out the waters that overlap with the lipids.
347 +  
348 +  delete[] isActive;
349 +  isActive = new int[newWaters];
350 +  for(i=0; i<newWaters; i++) isActive[i] = 1;
351 +  int n_active = newWaters;
352 +  rCutSqr = water_padding * water_padding;
353 +  
354 +  for(i=0; ( (i<newWaters) && isActive[i] ); i++){
355 +    for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
356  
357 <  } else if ( !strcasecmp( ensemble, "NVT") ) {
358 <    the_extendedsystem = new ExtendedSystem( simnfo );
359 <    the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
357 >      dx = waterSite[i].pos[0] - rsaAtoms[j]->getX();
358 >      dy = waterSite[i].pos[1] - rsaAtoms[j]->getY();
359 >      dz = waterSite[i].pos[2] - rsaAtoms[j]->getZ();
360  
361 <    if (the_globals->haveTauThermostat())
362 <      the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
363 <    else if (the_globals->haveQmass())
364 <      the_extendedsystem->setQmass(the_globals->getQmass());
365 <    else {
366 <      sprintf( painCave.errMsg,
367 <               "SimSetup error: If you use one of the constant temperature\n"
368 <               "    ensembles, you must set either tauThermostat or qMass.\n"
369 <               "    Neither of these was found in the BASS file.\n");
370 <      painCave.isFatal = 1;
371 <      simError();
361 >      map( dx, dy, dz, box_x, box_y, box_z );
362 >
363 >      dx2 = dx * dx;
364 >      dy2 = dy * dy;
365 >      dz2 = dz * dz;
366 >      
367 >      dSqr = dx2 + dy2 + dz2;
368 >      if( dSqr < rCutSqr ){
369 >        isActive[i] = 0;
370 >        n_active--;
371 >      }
372      }
373 +  }
374  
375 <  } else if ( !strcasecmp( ensemble, "NVE") ) {
376 <  } else {
375 >  if( n_active < nWaters ){
376 >    
377      sprintf( painCave.errMsg,
378 <             "SimSetup Warning. Unrecognized Ensemble -> %s, "
379 <             "reverting to NVE for this simulation.\n",
380 <             ensemble );
146 <    painCave.isFatal = 0;
378 >             "Too many waters were removed, edit code and try again.\n" );
379 >    
380 >    painCave.isFatal = 1;
381      simError();
382 <    strcpy( ensemble, "NVE" );
149 <  }  
150 <  strcpy( simnfo->ensemble, ensemble );
382 >  }
383  
384 +  int quickKill;
385 +  while( n_active > nWaters ){
386  
387 +    quickKill = (int)(drand48()*newWaters);
388  
389 <  delete stamps;
390 <  delete globals;
389 >    if( isActive[quickKill] ){
390 >      isActive[quickKill] = 0;
391 >      n_active--;
392 >    }
393 >  }
394 >
395 >  if( n_active != nWaters ){
396 >    
397 >    sprintf( painCave.errMsg,
398 >             "QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
399 >             n_active, nWaters );
400 >    painCave.isFatal = 1;
401 >    simError();
402 >  }
403 >
404 >  // clean up our messes before building the final system.
405 >
406 >  for(i=0; i<nAtoms; i++){
407 >    
408 >    delete atoms[i];
409 >  }
410 >  Atom::destroyArrays();
411 >  
412 >
413 >  // create the real Atom arrays
414 >  
415 >  nAtoms = 0;
416 >  molIndex = 0;
417 >  locate = new MoLocator*[2];
418 >  molSeq = new int[nLipids + nWaters];
419 >  molStart = new int[nLipids + nWaters];  
420 >  
421 >  locate[0] = lipidLocate;
422 >  for(j=0; j<nLipids; j++){
423 >    molSeq[molIndex] = 0;
424 >    molStart[molIndex] = nAtoms;
425 >    molIndex++;
426 >    nAtoms += lipidNatoms;
427 >  }
428 >
429 >  locate[1] = waterLocate;
430 >  for(j=0; j<nLipids; j++){
431 >    molSeq[molIndex] = 1;
432 >    molStart[molIndex] = nAtoms;
433 >    molIndex++;
434 >    nAtoms += waterNatoms;
435 >  }
436 >  
437 >  
438 >  Atom::createArrays( nAtoms );
439 >  atoms = new Atom*[nAtoms];
440 >
441 >  
442 >  // initialize lipid positions
443 >
444 >  
445 >  
446 >  
447 >
448 >  // set up the SimInfo object
449 >  
450 >  bsInfo.boxX = box_x;
451 >  bsInfo.boxY = box_y;
452 >  bsInfo.boxZ = box_z;
453 >  
454 >  simnfo = new SimInfo();
455 >  simnfo->n_atoms = nAtoms;
456 >  simnfo->box_x = bsInfo.boxX;
457 >  simnfo->box_y = bsInfo.boxY;
458 >  simnfo->box_z = bsInfo.boxZ;
459 >  
460 >  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
461 >  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
462 >
463 >  simnfo->atoms = atoms;
464 >  
465 >  // set up the writer and write out
466 >  
467 >  writer = new DumpWriter( simnfo );
468 >  writer->writeFinal();
469 >    
470 >  // clean up the memory
471 >  
472 >  if( molMap != NULL )   delete[] molMap;
473 >  if( cardDeck != NULL ) delete[] cardDeck;
474 >  if( locate != NULL ){
475 >    for(i=0; i<bsInfo.nComponents; i++){
476 >      delete locate[i];
477 >    }
478 >    delete[] locate;
479 >  }
480 >  if( atoms != NULL ){
481 >    for(i=0; i<nAtoms; i++){
482 >      delete atoms[i];
483 >    }
484 >    Atom::destroyArrays();
485 >    delete[] atoms;
486 >  }
487 >  if( molSeq != NULL ) delete[] molSeq;
488 >  if( simnfo != NULL ) delete simnfo;
489 >  if( writer != NULL ) delete writer;
490 >
491 >  return 1;
492   }
493 +
494 +  
495 +
496 +
497 + }
498 +
499 +
500 +
501 + int Old_buildRandomBilayer( void ){
502 +
503 +  int i,j,k;
504 +  int nAtoms, atomIndex, molIndex, molID;
505 +  int* molSeq;
506 +  int* molMap;
507 +  int* molStart;
508 +  int* cardDeck;
509 +  int deckSize;
510 +  int rSite, rCard;
511 +  double cell;
512 +  int nCells, nSites, siteIndex;
513 +  double rot[3][3];
514 +  double pos[3];
515 +
516 +  Atom** atoms;
517 +  SimInfo* simnfo;
518 +  DumpWriter* writer;
519 +  MoLocator** locate;
520 +
521 +  // initialize functions and variables
522 +  
523 +  srand48( RAND_SEED );
524 +  molSeq = NULL;
525 +  molStart = NULL;
526 +  molMap = NULL;
527 +  cardDeck = NULL;
528 +  atoms = NULL;
529 +  locate = NULL;
530 +  simnfo = NULL;
531 +  writer = NULL;
532 +
533 +  // calculate the number of cells in the fcc box
534 +
535 +  nCells = 0;
536 +  nSites = 0;
537 +  while( nSites < bsInfo.totNmol ){
538 +    nCells++;
539 +    nSites = 4.0 * pow( (double)nCells, 3.0 );
540 +  }
541 +
542 +
543 +  // create the molMap and cardDeck arrays
544 +  
545 +  molMap = new int[nSites];
546 +  cardDeck = new int[nSites];
547 +  
548 +  for(i=0; i<nSites; i++){
549 +    molMap[i] = -1;
550 +    cardDeck[i] = i;
551 +  }
552 +
553 +  // randomly place the molecules on the sites
554 +  
555 +  deckSize = nSites;
556 +  for(i=0; i<bsInfo.totNmol; i++){
557 +    rCard = (int)( deckSize * drand48() );
558 +    rSite = cardDeck[rCard];
559 +    molMap[rSite] = i;
560 +
561 +    // book keep the card deck;
562 +    
563 +    deckSize--;
564 +    cardDeck[rCard] = cardDeck[deckSize];
565 +  }
566 +  
567 +  
568 +  // create the MoLocator and Atom arrays
569 +  
570 +  nAtoms = 0;
571 +  molIndex = 0;
572 +  locate = new MoLocator*[bsInfo.nComponents];
573 +  molSeq = new int[bsInfo.totNmol];
574 +  molStart = new int[bsInfo.totNmol];  
575 +  for(i=0; i<bsInfo.nComponents; i++){
576 +    locate[i] = new MoLocator( bsInfo.compStamps[i] );
577 +    for(j=0; j<bsInfo.componentsNmol[i]; j++){
578 +      molSeq[molIndex] = i;
579 +      molStart[molIndex] = nAtoms;
580 +      molIndex++;
581 +      nAtoms += bsInfo.compStamps[i]->getNAtoms();
582 +    }
583 +  }
584 +  
585 +  Atom::createArrays( nAtoms );
586 +  atoms = new Atom*[nAtoms];
587 +  
588 +  
589 +  // place the molecules at each FCC site
590 +  
591 +  cell = 5.0;
592 +  for(i=0; i<bsInfo.nComponents; i++){
593 +    if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
594 +  }
595 +  cell *= 1.2; // add a little buffer
596 +
597 +  cell *= M_SQRT2;
598 +
599 +  siteIndex = 0;
600 +  for(i=0; i<nCells; i++){
601 +    for(j=0; j<nCells; j++){
602 +      for(k=0; k<nCells; k++){
603 +        
604 +        if( molMap[siteIndex] >= 0 ){
605 +          pos[0] = i * cell;
606 +          pos[1] = j * cell;
607 +          pos[2] = k * cell;
608 +          
609 +          getRandomRot( rot );
610 +          molID = molSeq[molMap[siteIndex]];
611 +          atomIndex = molStart[ molMap[siteIndex] ];
612 +          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
613 +        }
614 +        siteIndex++;
615 +
616 +        if( molMap[siteIndex] >= 0 ){
617 +          pos[0] = i * cell + (0.5 * cell);
618 +          pos[1] = j * cell;
619 +          pos[2] = k * cell + (0.5 * cell);
620 +
621 +          getRandomRot( rot );
622 +          molID = molSeq[molMap[siteIndex]];
623 +          atomIndex = molStart[ molMap[siteIndex] ];
624 +          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
625 +        }
626 +        siteIndex++;
627 +
628 +        if( molMap[siteIndex] >= 0 ){
629 +          pos[0] = i * cell + (0.5 * cell);
630 +          pos[1] = j * cell + (0.5 * cell);
631 +          pos[2] = k * cell;
632 +          
633 +          getRandomRot( rot );
634 +          molID = molSeq[molMap[siteIndex]];
635 +          atomIndex = molStart[ molMap[siteIndex] ];
636 +          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
637 +        }
638 +        siteIndex++;
639 +
640 +        if( molMap[siteIndex] >= 0 ){
641 +          pos[0] = i * cell;
642 +          pos[1] = j * cell + (0.5 * cell);
643 +          pos[2] = k * cell + (0.5 * cell);
644 +
645 +          getRandomRot( rot );
646 +          molID = molSeq[molMap[siteIndex]];
647 +          atomIndex = molStart[ molMap[siteIndex] ];
648 +          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
649 +        }
650 +        siteIndex++;
651 +      }
652 +    }
653 +  }
654 +
655 +  // set up the SimInfo object
656 +  
657 +  bsInfo.boxX = nCells * cell;
658 +  bsInfo.boxY = nCells * cell;
659 +  bsInfo.boxZ = nCells * cell;
660 +  
661 +  simnfo = new SimInfo();
662 +  simnfo->n_atoms = nAtoms;
663 +  simnfo->box_x = bsInfo.boxX;
664 +  simnfo->box_y = bsInfo.boxY;
665 +  simnfo->box_z = bsInfo.boxZ;
666 +  
667 +  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
668 +  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
669 +
670 +  simnfo->atoms = atoms;
671 +  
672 +  // set up the writer and write out
673 +  
674 +  writer = new DumpWriter( simnfo );
675 +  writer->writeFinal();
676 +    
677 +  // clean up the memory
678 +  
679 +  if( molMap != NULL )   delete[] molMap;
680 +  if( cardDeck != NULL ) delete[] cardDeck;
681 +  if( locate != NULL ){
682 +    for(i=0; i<bsInfo.nComponents; i++){
683 +      delete locate[i];
684 +    }
685 +    delete[] locate;
686 +  }
687 +  if( atoms != NULL ){
688 +    for(i=0; i<nAtoms; i++){
689 +      delete atoms[i];
690 +    }
691 +    Atom::destroyArrays();
692 +    delete[] atoms;
693 +  }
694 +  if( molSeq != NULL ) delete[] molSeq;
695 +  if( simnfo != NULL ) delete simnfo;
696 +  if( writer != NULL ) delete writer;
697 +
698 +  return 1;
699 + }
700 +
701 +
702 + void getRandomRot( double rot[3][3] ){
703 +
704 +  double theta, phi, psi;
705 +  double cosTheta;
706 +
707 +  // select random phi, psi, and cosTheta
708 +
709 +  phi = 2.0 * M_PI * drand48();
710 +  psi = 2.0 * M_PI * drand48();
711 +  cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
712 +
713 +  theta = acos( cosTheta );
714 +
715 +  rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
716 +  rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
717 +  rot[0][2] = sin(theta) * sin(psi);
718 +  
719 +  rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
720 +  rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
721 +  rot[1][2] = sin(theta) * cos(psi);
722 +
723 +  rot[2][0] = sin(phi) * sin(theta);
724 +  rot[2][1] = -cos(phi) * sin(theta);
725 +  rot[2][2] = cos(theta);  
726 + }
727 +                        
728 +
729 +
730 + void map( double &x, double &y, double &z,
731 +          double boxX, double boxY, double boxZ ){
732 +  
733 +  if(x < 0) x -= boxX * (double)( (int)( (x / boxX)  - 0.5 ) );
734 +  else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
735 +
736 +  if(y < 0) y -= boxY * (double)( (int)( (y / boxY)  - 0.5 ) );
737 +  else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
738 +  
739 +  if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ)  - 0.5 ) );
740 +  else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
741 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines