--- trunk/makeLipid/makeRandomLipid.cpp 2002/07/15 20:28:33 28 +++ trunk/makeLipid/makeRandomLipid.cpp 2002/07/17 20:33:56 32 @@ -2,6 +2,7 @@ #include #include #include +#include #include "SimSetup.hpp" #include "SimInfo.hpp" @@ -11,6 +12,12 @@ #include "ReadWrite.hpp" +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; @@ -33,15 +40,15 @@ int main(int argc,char* argv[]){ const double water_vol = 4.0 / water_rho; // volume occupied by 4 waters const double water_cell = 4.929; // fcc unit cell length - int n_lipidsX = 5; - int n_lipidsY = 10; - int n_lipids = n_lipidsX * n_lipidsY; + int n_lipids = 50; + double water_ratio = 25.0; // water to lipid ratio + int n_h2o_target = (int)( n_lipids * water_ratio + 0.5 ); std::cerr << "n_lipids = " << n_lipids << "\n"; 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. @@ -67,83 +74,22 @@ int main(int argc,char* argv[]){ lipidAtoms = entry_plug->atoms; lipidNAtoms = entry_plug->n_atoms; - int group_nAtoms = n_lipids * lipidNAtoms; - Atom** group_atoms = new Atom*[group_nAtoms]; - DirectionalAtom* dAtom; - DirectionalAtom* dAtomNew; - - double rotMat[3][3]; - - rotMat[0][0] = 1.0; - rotMat[0][1] = 0.0; - rotMat[0][2] = 0.0; - - rotMat[1][0] = 0.0; - rotMat[1][1] = 1.0; - rotMat[1][2] = 0.0; - - rotMat[2][0] = 0.0; - rotMat[2][1] = 0.0; - rotMat[2][2] = 1.0; - int index =0; - for(i=0; iisDirectional() ){ - dAtom = (DirectionalAtom *)lipidAtoms[j]; - - dAtomNew = new DirectionalAtom(); - dAtomNew->setSUx( dAtom->getSUx() ); - dAtomNew->setSUx( dAtom->getSUx() ); - dAtomNew->setSUx( dAtom->getSUx() ); - - dAtomNew->setA( rotMat ); - - group_atoms[index] = dAtomNew; - } - else{ - - group_atoms[index] = new GeneralAtom(); - } - - group_atoms[index]->setType( lipidAtoms[j]->getType() ); - - index++; - } - } + // find the width, height, and length of the molecule - index = 0; - for(i=0; isetX( lipidAtoms[l]->getX() + - i*lipid_spaceing ); - - group_atoms[index]->setY( lipidAtoms[l]->getY() + - j*lipid_spaceing ); - - group_atoms[index]->setZ( lipidAtoms[l]->getZ() ); - - index++; - } - } - } - double min_x, min_y, min_z; double max_x, max_y, max_z; double test_x, test_y, test_z; - max_x = min_x = group_atoms[0]->getX(); - max_y = min_y = group_atoms[0]->getY(); - max_z = min_z = group_atoms[0]->getZ(); + max_x = min_x = lipidAtoms[0]->getX(); + max_y = min_y = lipidAtoms[0]->getY(); + max_z = min_z = lipidAtoms[0]->getZ(); - for(i=0; igetX(); - test_y = group_atoms[i]->getY(); - test_z = group_atoms[i]->getZ(); + test_x = lipidAtoms[i]->getX(); + test_y = lipidAtoms[i]->getY(); + test_z = lipidAtoms[i]->getZ(); if( test_x < min_x ) min_x = test_x; if( test_y < min_y ) min_y = test_y; @@ -154,14 +100,24 @@ int main(int argc,char* argv[]){ if( test_z > max_z ) max_z = test_z; } - double box_x = max_x - min_x + 2*water_shell; - double box_y = max_y - min_y + 2*water_shell; - double box_z = max_z - min_z + 2*water_shell; + double ml2 = pow((max_x - min_x), 2 ) + pow((max_y - min_y), 2 ) + + pow((max_x - min_x), 2 ); + double max_length = sqrt( ml2 ); - int n_cellX = (int)(box_x / water_cell + 0.5 ); - int n_cellY = (int)(box_y / water_cell + 0.5 ); - int n_cellZ = (int)(box_z / water_cell + 0.5 ); + // from this information, create the test box + + + double box_x; + double box_y; + double box_z; + + box_x = box_y = box_z = max_length + water_cell * 4.0; // pad with 4 cells + + int n_cellX = (int)(box_x / water_cell + 1.0 ); + int n_cellY = (int)(box_y / water_cell + 1.0 ); + int n_cellZ = (int)(box_z / water_cell + 1.0 ); + box_x = water_cell * n_cellX; box_y = water_cell * n_cellY; box_z = water_cell * n_cellZ; @@ -172,67 +128,81 @@ int main(int argc,char* argv[]){ double *waterY = new double[n_water]; double *waterZ = new double[n_water]; + + // find the center of the test lipid, and make it the center of our + // soon to be created water box. + + double cx, cy, cz; cx = 0.0; cy = 0.0; cz = 0.0; - for(i=0; igetX(); - cy += group_atoms[i]->getY(); - cz += group_atoms[i]->getZ(); + for(i=0; igetX(); + cy += lipidAtoms[i]->getY(); + cz += lipidAtoms[i]->getZ(); } - cx /= group_nAtoms; - cy /= group_nAtoms; - cz /= group_nAtoms; + cx /= lipidNAtoms; + cy /= lipidNAtoms; + cz /= lipidNAtoms; double x0 = cx - ( box_x * 0.5 ); double y0 = cy - ( box_y * 0.5 ); double z0 = cz - ( box_z * 0.5 ); + + + // 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[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 + 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 + 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 + x0; - waterY[index] = j * water_cell + 0.5 * water_cell + y0; - waterZ[index] = k * water_cell + 0.5 * water_cell + z0; - index++; - - 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++; } } } + + // calculate the number of water's displaced by our molecule. + int *isActive = new int[n_water]; for(i=0; igetX(); - dy = waterY[i] - group_atoms[j]->getY(); - dz = waterZ[i] - group_atoms[j]->getZ(); + dx = waterX[i] - lipidAtoms[j]->getX(); + dy = waterY[i] - lipidAtoms[j]->getY(); + dz = waterZ[i] - lipidAtoms[j]->getZ(); + map( dx, dy, dz, box_x, box_y, box_z ); + dx2 = dx * dx; dy2 = dy * dy; dz2 = dz * dz; @@ -240,6 +210,280 @@ int main(int argc,char* argv[]){ dSqr = dx2 + dy2 + dz2; if( dSqr < rCutSqr ){ isActive[i] = 0; + n_deleted++; + } + } + } + + 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_water < n_h2o_target ){ + + int n_generated = n_cellX; + int n_test, nx, ny, nz; + nx = n_cellX; + ny = n_cellY; + nz = n_cellZ; + + n_test = 4 * nx * ny * nz; + + while( n_test < n_h2o_target ){ + + nz++; + n_test = 4 * nx * ny * nz; + } + + int n_diff, goodX, goodY, goodZ; + + n_diff = n_test - n_h2o_target; + goodX = nx; + goodY = ny; + goodZ = nz; + + int test_diff; + int n_limit = nz; + nz = n_cellZ; + + for( i=n_generated; i<=n_limit; i++ ){ + for( j=i; j<=n_limit; j++ ){ + for( k=j; k<=n_limit; k++ ){ + + n_test = 4 * i * j * k; + + if( n_test > n_h2o_target ){ + + test_diff = n_test - n_h2o_target; + + if( test_diff < n_diff ){ + + n_diff = test_diff; + goodX = nx; + goodY = ny; + goodZ = nz; + } + } + } + } + } + + n_cellX = goodX; + n_cellY = goodY; + n_cellZ = goodZ; + } + + // 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--; } } @@ -247,68 +491,68 @@ int main(int argc,char* argv[]){ 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++; } } @@ -354,3 +598,46 @@ int main(int argc,char* argv[]){ return 0; } + + +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)); + + 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)); +} + + +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); + + 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; +}