--- trunk/makeLipid/makeRandomLipid.cpp 2002/07/16 01:22:04 29 +++ trunk/makeLipid/makeRandomLipid.cpp 2002/07/17 20:33:56 32 @@ -12,9 +12,12 @@ void map( double *x, double *y, double *z, double cent #include "ReadWrite.hpp" -void map( double *x, double *y, double *z, double centerX, double centerY, - double centerZ, double boxX, double boxY, double boxZ ); +void map( double &x, double &y, double &z, + double boxX, double boxY, double boxZ ); +void rotate( double &x, double &y, double &z, + double theta, double phi, double psi ); + char* program_name; using namespace std; @@ -45,7 +48,7 @@ int main(int argc,char* argv[]){ double water_shell = 10.0; double water_padding = 2.5; - double lipid_spaceing = 4.0; + double lipid_spaceing = 2.5; srand48( 1337 ); // initialize the random number generator. @@ -152,30 +155,30 @@ int main(int argc,char* argv[]){ // create an fcc lattice in the water box. - index = 0; + int ndx = 0; for( i=0; i < n_cellX; i++ ){ for( j=0; j < n_cellY; j++ ){ for( k=0; k < n_cellZ; k++ ){ - waterX[index] = i * water_cell + x0; - waterY[index] = j * water_cell + y0; - waterZ[index] = k * water_cell + z0; - index++; + waterX[ndx] = i * water_cell + x0; + waterY[ndx] = j * water_cell + y0; + waterZ[ndx] = k * water_cell + z0; + ndx++; - waterX[index] = i * water_cell + 0.5 * water_cell + x0; - waterY[index] = j * water_cell + 0.5 * water_cell + y0; - waterZ[index] = k * water_cell + z0; - index++; + waterX[ndx] = i * water_cell + 0.5 * water_cell + x0; + waterY[ndx] = j * water_cell + 0.5 * water_cell + y0; + waterZ[ndx] = k * water_cell + z0; + ndx++; - waterX[index] = i * water_cell + x0; - waterY[index] = j * water_cell + 0.5 * water_cell + y0; - waterZ[index] = k * water_cell + 0.5 * water_cell + z0; - index++; + waterX[ndx] = i * water_cell + x0; + waterY[ndx] = j * water_cell + 0.5 * water_cell + y0; + waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0; + ndx++; - waterX[index] = i * water_cell + 0.5 * water_cell + x0; - waterY[index] = j * water_cell + y0; - waterZ[index] = k * water_cell + 0.5 * water_cell + z0; - index++; + waterX[ndx] = i * water_cell + 0.5 * water_cell + x0; + waterY[ndx] = j * water_cell + y0; + waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0; + ndx++; } } } @@ -198,7 +201,7 @@ int main(int argc,char* argv[]){ dy = waterY[i] - lipidAtoms[j]->getY(); dz = waterZ[i] - lipidAtoms[j]->getZ(); - map( &dx, &dy, &dz, cx, cy, cz, box_x, box_y, box_z ); + map( dx, dy, dz, box_x, box_y, box_z ); dx2 = dx * dx; dy2 = dy * dy; @@ -211,14 +214,19 @@ int main(int argc,char* argv[]){ } } } - + + std::cerr << "nTarget before: " << n_h2o_target; + n_h2o_target += n_deleted * n_lipids; + + std::cerr << ", after: " << n_h2o_target << ", n_deleted: " << n_deleted + << "\n"; // find a box size that best suits the number of waters we need. int done = 0; - if( n_waters < n_h2o_target ){ + if( n_water < n_h2o_target ){ int n_generated = n_cellX; int n_test, nx, ny, nz; @@ -226,6 +234,8 @@ int main(int argc,char* argv[]){ ny = n_cellY; nz = n_cellZ; + n_test = 4 * nx * ny * nz; + while( n_test < n_h2o_target ){ nz++; @@ -234,14 +244,14 @@ int main(int argc,char* argv[]){ int n_diff, goodX, goodY, goodZ; - n_diff = ntest - n_h2o_target; + n_diff = n_test - n_h2o_target; goodX = nx; goodY = ny; goodZ = nz; int test_diff; - int n_limit = n_z; - n_z = n_cellZ; + int n_limit = nz; + nz = n_cellZ; for( i=n_generated; i<=n_limit; i++ ){ for( j=i; j<=n_limit; j++ ){ @@ -270,79 +280,279 @@ int main(int argc,char* argv[]){ n_cellZ = goodZ; } - // we now have the best box size for the simulation. - - + // we now have the best box size for the simulation. Next we + // recreate the water box to the new specifications. + + n_water = n_cellX * n_cellY * n_cellZ * 4; + + std::cerr << "new waters = " << n_water << "\n"; + + delete[] waterX; + delete[] waterY; + delete[] waterZ; + + waterX = new double[n_water]; + waterY = new double[n_water]; + waterZ = new double[n_water]; + + box_x = water_cell * n_cellX; + box_y = water_cell * n_cellY; + box_z = water_cell * n_cellZ; + + x0 = 0.0; + y0 = 0.0; + z0 = 0.0; + + cx = ( box_x * 0.5 ); + cy = ( box_y * 0.5 ); + cz = ( box_z * 0.5 ); + + // create an fcc lattice in the water box. + + ndx = 0; + for( i=0; i < n_cellX; i++ ){ + for( j=0; j < n_cellY; j++ ){ + for( k=0; k < n_cellZ; k++ ){ + + waterX[ndx] = i * water_cell + x0; + waterY[ndx] = j * water_cell + y0; + waterZ[ndx] = k * water_cell + z0; + ndx++; + + waterX[ndx] = i * water_cell + 0.5 * water_cell + x0; + waterY[ndx] = j * water_cell + 0.5 * water_cell + y0; + waterZ[ndx] = k * water_cell + z0; + ndx++; + + waterX[ndx] = i * water_cell + x0; + waterY[ndx] = j * water_cell + 0.5 * water_cell + y0; + waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0; + ndx++; + + waterX[ndx] = i * water_cell + 0.5 * water_cell + x0; + waterY[ndx] = j * water_cell + y0; + waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0; + ndx++; + } + } + } + + // ************************************************************** + + + + // start a 3D RSA for the for the lipid placements + + srand48( 1337 ); + + int rsaNAtoms = n_lipids * lipidNAtoms; + Atom** rsaAtoms = new Atom*[rsaNAtoms]; + + DirectionalAtom* dAtom; + DirectionalAtom* dAtomNew; + + double rotMat[3][3]; + double unitRotMat[3][3]; + + unitRotMat[0][0] = 1.0; + unitRotMat[0][1] = 0.0; + unitRotMat[0][2] = 0.0; + + unitRotMat[1][0] = 0.0; + unitRotMat[1][1] = 1.0; + unitRotMat[1][2] = 0.0; + + unitRotMat[2][0] = 0.0; + unitRotMat[2][1] = 0.0; + unitRotMat[2][2] = 1.0; + + ndx = 0; + for(i=0; iisDirectional() ){ + dAtom = (DirectionalAtom *)lipidAtoms[j]; + + dAtomNew = new DirectionalAtom(); + dAtomNew->setSUx( dAtom->getSUx() ); + dAtomNew->setSUx( dAtom->getSUx() ); + dAtomNew->setSUx( dAtom->getSUx() ); + + dAtom->getA( rotMat ); + dAtomNew->setA( rotMat ); + + rsaAtoms[ndx] = dAtomNew; + } + else{ + + rsaAtoms[ndx] = new GeneralAtom(); + } + + rsaAtoms[ndx]->setType( lipidAtoms[j]->getType() ); + + ndx++; + } + } + + double testX, testY, testZ; + double theta, phi, psi; + double tempX, tempY, tempZ; + int reject; + int testDX, acceptedDX; + + rCutSqr = lipid_spaceing * lipid_spaceing; + + for(i=0; igetX(); + tempY = lipidAtoms[j]->getY(); + tempZ = lipidAtoms[j]->getZ(); + + rotate( tempX, tempY, tempZ, theta, phi, psi ); + + rsaAtoms[ndx + j]->setX( tempX + testX ); + rsaAtoms[ndx + j]->setY( tempY + testY ); + rsaAtoms[ndx + j]->setZ( tempZ + testZ ); + } + + reject = 0; + for( j=0; !reject && jgetX() - rsaAtoms[acceptedDX]->getX(); + dy = rsaAtoms[testDX]->getY() - rsaAtoms[acceptedDX]->getY(); + dz = rsaAtoms[testDX]->getZ() - rsaAtoms[acceptedDX]->getZ(); + + map( dx, dy, dz, box_x, box_y, box_z ); + + dx2 = dx * dx; + dy2 = dy * dy; + dz2 = dz * dz; + + dSqr = dx2 + dy2 + dz2; + if( dSqr < rCutSqr ) reject = 1; + } + } + } + if( !reject ){ + done = 1; + std::cerr << i << " has been accepted\n"; + } + } + } + + // cut out the waters that overlap with the lipids. + + delete[] isActive; + isActive = new int[n_water]; + for(i=0; igetX(); + dy = waterY[i] - rsaAtoms[j]->getY(); + dz = waterZ[i] - rsaAtoms[j]->getZ(); + map( dx, dy, dz, box_x, box_y, box_z ); + dx2 = dx * dx; + dy2 = dy * dy; + dz2 = dz * dz; + + dSqr = dx2 + dy2 + dz2; + if( dSqr < rCutSqr ){ + isActive[i] = 0; + n_active--; + } + } + } + std::cerr << "final n_waters = " << n_active << "\n"; - int new_nAtoms = group_nAtoms + n_active; + // place all of the waters and lipids into one new array + + int new_nAtoms = rsaNAtoms + n_active; Atom** new_atoms = new Atom*[new_nAtoms]; - index = 0; - for(i=0; iisDirectional() ){ - dAtom = (DirectionalAtom *)group_atoms[i]; + if( rsaAtoms[i]->isDirectional() ){ + dAtom = (DirectionalAtom *)rsaAtoms[i]; dAtomNew = new DirectionalAtom(); dAtomNew->setSUx( dAtom->getSUx() ); dAtomNew->setSUx( dAtom->getSUx() ); dAtomNew->setSUx( dAtom->getSUx() ); + dAtom->getA( rotMat ); dAtomNew->setA( rotMat ); - new_atoms[index] = dAtomNew; + new_atoms[ndx] = dAtomNew; } else{ - new_atoms[index] = new GeneralAtom(); + new_atoms[ndx] = new GeneralAtom(); } - new_atoms[index]->setType( group_atoms[i]->getType() ); + new_atoms[ndx]->setType( rsaAtoms[i]->getType() ); - new_atoms[index]->setX( group_atoms[i]->getX() ); - new_atoms[index]->setY( group_atoms[i]->getY() ); - new_atoms[index]->setZ( group_atoms[i]->getZ() ); + new_atoms[ndx]->setX( rsaAtoms[i]->getX() ); + new_atoms[ndx]->setY( rsaAtoms[i]->getY() ); + new_atoms[ndx]->setZ( rsaAtoms[i]->getZ() ); - new_atoms[index]->set_vx( 0.0 ); - new_atoms[index]->set_vy( 0.0 ); - new_atoms[index]->set_vz( 0.0 ); + new_atoms[ndx]->set_vx( 0.0 ); + new_atoms[ndx]->set_vy( 0.0 ); + new_atoms[ndx]->set_vz( 0.0 ); - index++; + ndx++; } - - - for(i=0; isetType( "SSD" ); + new_atoms[ndx] = new DirectionalAtom(); + new_atoms[ndx]->setType( "SSD" ); - new_atoms[index]->setX( waterX[i] ); - new_atoms[index]->setY( waterY[i] ); - new_atoms[index]->setZ( waterZ[i] ); + new_atoms[ndx]->setX( waterX[i] ); + new_atoms[ndx]->setY( waterY[i] ); + new_atoms[ndx]->setZ( waterZ[i] ); - new_atoms[index]->set_vx( 0.0 ); - new_atoms[index]->set_vy( 0.0 ); - new_atoms[index]->set_vz( 0.0 ); + new_atoms[ndx]->set_vx( 0.0 ); + new_atoms[ndx]->set_vy( 0.0 ); + new_atoms[ndx]->set_vz( 0.0 ); - dAtom = (DirectionalAtom *) new_atoms[index]; + dAtom = (DirectionalAtom *) new_atoms[ndx]; dAtom->setSUx( 0.0 ); dAtom->setSUy( 0.0 ); dAtom->setSUz( 1.0 ); - dAtom->setA( rotMat ); + dAtom->setA( unitRotMat ); - index++; + ndx++; } } @@ -390,22 +600,44 @@ void map( x, y, z, centerX, centerY, centerZ, boxX, bo } -void map( x, y, z, centerX, centerY, centerZ, boxX, boxY, boxZ ) - double *x, *y, *z; - double centerX, centerY, centerZ; - double boxX, boxY, boxZ; -{ +void map( double &x, double &y, double &z, + double boxX, double boxY, double boxZ ){ + + if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) ); + else x -= boxX * (double)( (int)( (x / boxX ) + 0.5)); - *x -= centerX; - *y -= centerY; - *z -= centerZ; + if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) ); + else y -= boxY * (double)( (int)( (y / boxY ) + 0.5)); + + if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) ); + else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5)); +} - if(*x < 0) *x -= boxX * (double)( (int)( (*x / boxX) - 0.5 ) ); - else *x -= boxX * (double)( (int)( (*x / boxX ) + 0.5)); - if(*y < 0) *y -= boxY * (double)( (int)( (*y / boxY) - 0.5 ) ); - else *y -= boxY * (double)( (int)( (*y / boxY ) + 0.5)); +void rotate( double &x, double &y, double &z, + double theta, double phi, double psi ){ + + double newX, newY, newZ; + + double A[3][3]; + + A[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); + A[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); + A[0][2] = sin(theta) * sin(psi); - if(*z < 0) *z -= boxZ * (double)( (int)( (*z / boxZ) - 0.5 ) ); - else *z -= boxZ * (double)( (int)( (*z / boxZ ) + 0.5)); + A[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); + A[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); + A[1][2] = sin(theta) * cos(psi); + + A[2][0] = sin(phi) * sin(theta); + A[2][1] = -cos(phi) * sin(theta); + A[2][2] = cos(theta); + + newX = (x * A[0][0]) + (y * A[0][1]) + (z * A[0][2]); + newY = (x * A[1][0]) + (y * A[1][1]) + (z * A[1][2]); + newZ = (x * A[2][0]) + (y * A[2][1]) + (z * A[2][2]); + + x = newX; + y = newY; + z = newZ; }