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

Diff Legend

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