1 |
+ |
#include <iostream> |
2 |
+ |
|
3 |
|
#include <cstdlib> |
4 |
|
#include <cstring> |
5 |
|
#include <cmath> |
13 |
|
#include "bilayerSys.hpp" |
14 |
|
|
15 |
|
|
16 |
+ |
|
17 |
+ |
void map( double &x, double &y, double &z, |
18 |
+ |
double boxX, double boxY, double boxZ ); |
19 |
+ |
|
20 |
|
int buildRandomBilayer( void ); |
21 |
|
|
22 |
|
void getRandomRot( double rot[3][3] ); |
36 |
|
} |
37 |
|
|
38 |
|
|
33 |
– |
|
39 |
|
int buildRandomBilayer( void ){ |
40 |
|
|
41 |
+ |
typedef struct{ |
42 |
+ |
double rot[3][3]; |
43 |
+ |
double pos[3]; |
44 |
+ |
} coord; |
45 |
+ |
|
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; |
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 |
> |
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 |
> |
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 |
> |
//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 |
> |
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 |
> |
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 |
> |
|
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 |
> |
// calculate the number of water's displaced by our lipid. |
171 |
> |
|
172 |
> |
testSite.rot[0][0] = 1.0; |
173 |
> |
testSite.rot[0][1] = 0.0; |
174 |
> |
testSite.rot[0][2] = 0.0; |
175 |
> |
|
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 |
> |
} |
218 |
> |
|
219 |
> |
int targetWaters = nWaters + n_deleted * nLipids; |
220 |
> |
|
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 |
> |
} |
280 |
> |
|
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 |
> |
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 |
> |
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 |
> |
if( n_active < nWaters ){ |
376 |
> |
|
377 |
> |
sprintf( painCave.errMsg, |
378 |
> |
"Too many waters were removed, edit code and try again.\n" ); |
379 |
> |
|
380 |
> |
painCave.isFatal = 1; |
381 |
> |
simError(); |
382 |
> |
} |
383 |
> |
|
384 |
> |
int quickKill; |
385 |
> |
while( n_active > nWaters ){ |
386 |
> |
|
387 |
> |
quickKill = (int)(drand48()*newWaters); |
388 |
> |
|
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; |
522 |
|
|
523 |
|
srand48( RAND_SEED ); |
524 |
|
molSeq = NULL; |
525 |
+ |
molStart = NULL; |
526 |
|
molMap = NULL; |
527 |
|
cardDeck = NULL; |
528 |
|
atoms = NULL; |
536 |
|
nSites = 0; |
537 |
|
while( nSites < bsInfo.totNmol ){ |
538 |
|
nCells++; |
539 |
< |
nSites = 4 * pow( nCells, 3 ); |
539 |
> |
nSites = 4.0 * pow( (double)nCells, 3.0 ); |
540 |
|
} |
541 |
|
|
542 |
|
|
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.compStamp[i]->getNAtoms(); |
581 |
> |
nAtoms += bsInfo.compStamps[i]->getNAtoms(); |
582 |
|
} |
583 |
|
} |
584 |
|
|
592 |
|
for(i=0; i<bsInfo.nComponents; i++){ |
593 |
|
if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength(); |
594 |
|
} |
595 |
< |
cell *= M_SQRT_2; |
595 |
> |
cell *= 1.2; // add a little buffer |
596 |
|
|
597 |
+ |
cell *= M_SQRT2; |
598 |
+ |
|
599 |
|
siteIndex = 0; |
127 |
– |
atomIndex = 0; |
600 |
|
for(i=0; i<nCells; i++){ |
601 |
|
for(j=0; j<nCells; j++){ |
602 |
|
for(k=0; k<nCells; k++){ |
608 |
|
|
609 |
|
getRandomRot( rot ); |
610 |
|
molID = molSeq[molMap[siteIndex]]; |
611 |
< |
locate[molID]->placeMol( pos, rot, &atoms[atomIndex] ); |
612 |
< |
|
141 |
< |
atomIndex += bsInfo.compStamps[molID]->getNatoms(); |
611 |
> |
atomIndex = molStart[ molMap[siteIndex] ]; |
612 |
> |
locate[molID]->placeMol( pos, rot, atoms, atomIndex ); |
613 |
|
} |
614 |
|
siteIndex++; |
615 |
|
|
617 |
|
pos[0] = i * cell + (0.5 * cell); |
618 |
|
pos[1] = j * cell; |
619 |
|
pos[2] = k * cell + (0.5 * cell); |
620 |
< |
|
620 |
> |
|
621 |
|
getRandomRot( rot ); |
622 |
|
molID = molSeq[molMap[siteIndex]]; |
623 |
< |
locate[molID]->placeMol( pos, rot, &atoms[atomIndex] ); |
624 |
< |
|
154 |
< |
atomIndex += bsInfo.compStamps[molID]->getNatoms(); |
623 |
> |
atomIndex = molStart[ molMap[siteIndex] ]; |
624 |
> |
locate[molID]->placeMol( pos, rot, atoms, atomIndex ); |
625 |
|
} |
626 |
|
siteIndex++; |
627 |
|
|
632 |
|
|
633 |
|
getRandomRot( rot ); |
634 |
|
molID = molSeq[molMap[siteIndex]]; |
635 |
< |
locate[molID]->placeMol( pos, rot, &atoms[atomIndex] ); |
636 |
< |
|
167 |
< |
atomIndex += bsInfo.compStamps[molID]->getNatoms(); |
635 |
> |
atomIndex = molStart[ molMap[siteIndex] ]; |
636 |
> |
locate[molID]->placeMol( pos, rot, atoms, atomIndex ); |
637 |
|
} |
638 |
|
siteIndex++; |
639 |
|
|
641 |
|
pos[0] = i * cell; |
642 |
|
pos[1] = j * cell + (0.5 * cell); |
643 |
|
pos[2] = k * cell + (0.5 * cell); |
644 |
< |
|
644 |
> |
|
645 |
|
getRandomRot( rot ); |
646 |
|
molID = molSeq[molMap[siteIndex]]; |
647 |
< |
locate[molID]->placeMol( pos, rot, &atoms[atomIndex] ); |
648 |
< |
|
180 |
< |
atomIndex += bsInfo.compStamps[molID]->getNatoms(); |
647 |
> |
atomIndex = molStart[ molMap[siteIndex] ]; |
648 |
> |
locate[molID]->placeMol( pos, rot, atoms, atomIndex ); |
649 |
|
} |
650 |
|
siteIndex++; |
651 |
|
} |
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; |
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.statusName, "%s.dump", bsInfo.outPrefix ); |
668 |
< |
sprintf( simnfo.finalName, "%s.init", bsInfo.outPrefix ); |
667 |
> |
sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix ); |
668 |
> |
sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix ); |
669 |
> |
|
670 |
> |
simnfo->atoms = atoms; |
671 |
|
|
202 |
– |
simnfo.atoms = atoms; |
203 |
– |
|
672 |
|
// set up the writer and write out |
673 |
|
|
674 |
< |
writer = new DumpWriter( &simnfo ); |
674 |
> |
writer = new DumpWriter( simnfo ); |
675 |
|
writer->writeFinal(); |
676 |
|
|
677 |
|
// clean up the memory |
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); |
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)); |
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 |
+ |
} |