| 1 |
+ |
|
| 2 |
|
#include <iostream> |
| 3 |
+ |
#include <vector> |
| 4 |
+ |
#include <algorithm> |
| 5 |
|
|
| 6 |
|
#include <cstdlib> |
| 7 |
|
#include <cstring> |
| 8 |
|
#include <cmath> |
| 9 |
|
|
| 10 |
+ |
|
| 11 |
|
#include "simError.h" |
| 12 |
|
#include "SimInfo.hpp" |
| 13 |
|
#include "ReadWrite.hpp" |
| 18 |
|
|
| 19 |
|
#include "latticeBuilder.hpp" |
| 20 |
|
|
| 21 |
+ |
class SortCond{ |
| 22 |
+ |
|
| 23 |
+ |
public: |
| 24 |
+ |
bool operator()(const pair<int, double>& p1, const pair<int, double>& p2){ |
| 25 |
+ |
return p1.second < p2.second; |
| 26 |
+ |
} |
| 27 |
+ |
|
| 28 |
+ |
|
| 29 |
+ |
}; |
| 30 |
+ |
|
| 31 |
+ |
|
| 32 |
|
void buildMap( double &x, double &y, double &z, |
| 33 |
|
double boxX, double boxY, double boxZ ); |
| 34 |
|
|
| 35 |
|
int buildRandomBilayer( void ); |
| 36 |
|
|
| 22 |
– |
void getRandomRot( double rot[3][3] ); |
| 23 |
– |
|
| 37 |
|
int buildBilayer( int isRandom ){ |
| 38 |
|
|
| 39 |
|
if( isRandom ){ |
| 40 |
|
return buildRandomBilayer(); |
| 41 |
|
} |
| 42 |
|
else{ |
| 43 |
< |
sprintf( painCave.errMsg, |
| 44 |
< |
"Cannot currently create a non-random bilayer.\n" ); |
| 32 |
< |
painCave.isFatal = 1; |
| 33 |
< |
simError(); |
| 43 |
> |
|
| 44 |
> |
std::cerr << "unsupported feature\n"; |
| 45 |
|
return 0; |
| 46 |
|
} |
| 47 |
|
} |
| 245 |
|
targetWaters = (int) ( targetWaters * 1.2 ); |
| 246 |
|
|
| 247 |
|
// find the best box size for the sim |
| 248 |
+ |
|
| 249 |
+ |
int nCellsX, nCellsY, nCellsZ; |
| 250 |
+ |
|
| 251 |
+ |
const double boxTargetX = 66.22752; |
| 252 |
+ |
const double boxTargetY = 60.53088; |
| 253 |
|
|
| 254 |
+ |
nCellsX = (int)ceil(boxTargetX / waterCell); |
| 255 |
+ |
nCellsY = (int)ceil(boxTargetY / waterCell); |
| 256 |
+ |
|
| 257 |
|
int testTot; |
| 258 |
|
int done = 0; |
| 259 |
< |
ndx = 0; |
| 259 |
> |
nCellsZ = 0; |
| 260 |
|
while( !done ){ |
| 261 |
|
|
| 262 |
< |
ndx++; |
| 263 |
< |
testTot = 4 * ndx * ndx * ndx; |
| 262 |
> |
nCellsZ++; |
| 263 |
> |
testTot = 4 * nCellsX * nCellsY * nCellsZ; |
| 264 |
|
|
| 265 |
|
if( testTot >= targetWaters ) done = 1; |
| 266 |
|
} |
| 267 |
|
|
| 249 |
– |
nCells = ndx; |
| 250 |
– |
|
| 251 |
– |
|
| 268 |
|
// create the new water box to the new specifications |
| 269 |
|
|
| 270 |
< |
int newWaters = nCells * nCells * nCells * 4; |
| 270 |
> |
int newWaters = nCellsX * nCellsY * nCellsZ * 4; |
| 271 |
|
|
| 272 |
|
delete[] waterX; |
| 273 |
|
delete[] waterY; |
| 275 |
|
|
| 276 |
|
coord* waterSites = new coord[newWaters]; |
| 277 |
|
|
| 278 |
< |
double box_x = waterCell * nCells; |
| 279 |
< |
double box_y = waterCell * nCells; |
| 280 |
< |
double box_z = waterCell * nCells; |
| 278 |
> |
double box_x = waterCell * nCellsX; |
| 279 |
> |
double box_y = waterCell * nCellsY; |
| 280 |
> |
double box_z = waterCell * nCellsZ; |
| 281 |
|
|
| 282 |
|
// create an fcc lattice in the water box. |
| 283 |
|
|
| 284 |
|
ndx = 0; |
| 285 |
< |
for( i=0; i < nCells; i++ ){ |
| 286 |
< |
for( j=0; j < nCells; j++ ){ |
| 287 |
< |
for( k=0; k < nCells; k++ ){ |
| 285 |
> |
for( i=0; i < nCellsX; i++ ){ |
| 286 |
> |
for( j=0; j < nCellsY; j++ ){ |
| 287 |
> |
for( k=0; k < nCellsZ; k++ ){ |
| 288 |
|
|
| 289 |
|
myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k); |
| 290 |
|
for(l=0; l<4; l++){ |
| 371 |
|
} |
| 372 |
|
} |
| 373 |
|
|
| 374 |
+ |
|
| 375 |
+ |
// zSort of the lipid positions |
| 376 |
+ |
|
| 377 |
+ |
|
| 378 |
+ |
vector< pair<int,double> >zSortArray; |
| 379 |
+ |
for(i=0;i<nLipids;i++) |
| 380 |
+ |
zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) ); |
| 381 |
+ |
|
| 382 |
+ |
sort(zSortArray.begin(),zSortArray.end(),SortCond()); |
| 383 |
+ |
|
| 384 |
+ |
ofstream outFile( "./zipper.bass", ios::app); |
| 385 |
+ |
|
| 386 |
+ |
for(i=0; i<nLipids; i++){ |
| 387 |
+ |
outFile << "zConstraint[" << i << "]{\n" |
| 388 |
+ |
<< " molIndex = " << zSortArray[i].first << ";\n" |
| 389 |
+ |
<< " zPos = "; |
| 390 |
+ |
|
| 391 |
+ |
if(i<32) outFile << "60.0;\n"; |
| 392 |
+ |
else outFile << "100.0;\n"; |
| 393 |
+ |
|
| 394 |
+ |
outFile << " kRatio = 0.5;\n" |
| 395 |
+ |
<< "}\n"; |
| 396 |
+ |
} |
| 397 |
+ |
|
| 398 |
+ |
outFile.close(); |
| 399 |
+ |
|
| 400 |
+ |
|
| 401 |
|
// cut out the waters that overlap with the lipids. |
| 402 |
|
|
| 403 |
|
|
| 571 |
|
return 1; |
| 572 |
|
} |
| 573 |
|
|
| 531 |
– |
void getRandomRot( double rot[3][3] ){ |
| 532 |
– |
|
| 533 |
– |
double theta, phi, psi; |
| 534 |
– |
double cosTheta; |
| 535 |
– |
|
| 536 |
– |
// select random phi, psi, and cosTheta |
| 537 |
– |
|
| 538 |
– |
phi = 2.0 * M_PI * drand48(); |
| 539 |
– |
psi = 2.0 * M_PI * drand48(); |
| 540 |
– |
cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1 |
| 541 |
– |
|
| 542 |
– |
theta = acos( cosTheta ); |
| 543 |
– |
|
| 544 |
– |
rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); |
| 545 |
– |
rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); |
| 546 |
– |
rot[0][2] = sin(theta) * sin(psi); |
| 547 |
– |
|
| 548 |
– |
rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); |
| 549 |
– |
rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); |
| 550 |
– |
rot[1][2] = sin(theta) * cos(psi); |
| 551 |
– |
|
| 552 |
– |
rot[2][0] = sin(phi) * sin(theta); |
| 553 |
– |
rot[2][1] = -cos(phi) * sin(theta); |
| 554 |
– |
rot[2][2] = cos(theta); |
| 555 |
– |
} |
| 556 |
– |
|
| 557 |
– |
|
| 558 |
– |
|
| 574 |
|
void buildMap( double &x, double &y, double &z, |
| 575 |
|
double boxX, double boxY, double boxZ ){ |
| 576 |
|
|
