--- trunk/makeLipid/makeRandomLipid.cpp 2002/07/15 20:28:33 28 +++ trunk/makeLipid/makeRandomLipid.cpp 2002/07/16 01:22:04 29 @@ -2,6 +2,7 @@ #include #include #include +#include #include "SimSetup.hpp" #include "SimInfo.hpp" @@ -11,6 +12,9 @@ char* program_name; #include "ReadWrite.hpp" +void map( double *x, double *y, double *z, double centerX, double centerY, + double centerZ, double boxX, double boxY, double boxZ ); + char* program_name; using namespace std; @@ -33,9 +37,9 @@ 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"; @@ -67,83 +71,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 +97,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,23 +125,32 @@ 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; for( i=0; i < n_cellX; i++ ){ @@ -218,21 +180,26 @@ int main(int argc,char* argv[]){ } } + + // 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, cx, cy, cz, box_x, box_y, box_z ); + dx2 = dx * dx; dy2 = dy * dy; dz2 = dz * dz; @@ -240,13 +207,80 @@ int main(int argc,char* argv[]){ dSqr = dx2 + dy2 + dz2; if( dSqr < rCutSqr ){ isActive[i] = 0; - n_active--; + n_deleted++; } } } - std::cerr << "final n_waters = " << n_active << "\n"; + n_h2o_target += n_deleted * n_lipids; + // find a box size that best suits the number of waters we need. + + int done = 0; + + if( n_waters < n_h2o_target ){ + + int n_generated = n_cellX; + int n_test, nx, ny, nz; + nx = n_cellX; + ny = n_cellY; + nz = n_cellZ; + + while( n_test < n_h2o_target ){ + + nz++; + n_test = 4 * nx * ny * nz; + } + + int n_diff, goodX, goodY, goodZ; + + n_diff = ntest - n_h2o_target; + goodX = nx; + goodY = ny; + goodZ = nz; + + int test_diff; + int n_limit = n_z; + n_z = 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. + + + + + + + + + + int new_nAtoms = group_nAtoms + n_active; Atom** new_atoms = new Atom*[new_nAtoms]; @@ -354,3 +388,24 @@ int main(int argc,char* argv[]){ return 0; } + + +void map( x, y, z, centerX, centerY, centerZ, boxX, boxY, boxZ ) + double *x, *y, *z; + double centerX, centerY, centerZ; + double boxX, boxY, boxZ; +{ + + *x -= centerX; + *y -= centerY; + *z -= centerZ; + + 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)); +}