--- trunk/OOPSE/utils/sysbuilder/latticeBilayer.cpp 2003/10/27 22:07:49 824 +++ trunk/OOPSE/utils/sysbuilder/latticeBilayer.cpp 2003/10/27 22:08:36 825 @@ -7,11 +7,12 @@ #include "simError.h" #include "SimInfo.hpp" #include "ReadWrite.hpp" +#include "SimSetup.hpp" #include "MoLocator.hpp" #include "latticeBuilder.hpp" -#include "QuickBass.hpp" + #define VERSION_MAJOR 0 #define VERSION_MINOR 1 @@ -19,8 +20,11 @@ int buildLatticeBilayer( double aLat, void usage(void); int buildLatticeBilayer( double aLat, double bLat, - double leafSpacing); + double leafSpacing, + char* waterName, + char* lipidName); using namespace std; +SimInfo* mainInfo; int main(int argC,char* argV[]){ @@ -45,6 +49,8 @@ int main(int argC,char* argV[]){ char* inName; + SimSetup* simInit; + // first things first, all of the initializations fflush(stdout); @@ -280,9 +286,11 @@ int main(int argC,char* argV[]){ simError(); } - bLat = strtod( argV[i], &conversionCheck); - if( conversionCheck == argV[i] ) conversionError = true; - if( *conversionCheck != '\0' ) conversionError = true; + + bLat = atof( argV[i] ); +// bLat = strtod( argV[i], &conversionCheck); +// if( conversionCheck == argV[i] ) conversionError = true; +// if( *conversionCheck != '\0' ) conversionError = true; if( conversionError ){ sprintf( painCave.errMsg, @@ -328,9 +336,10 @@ int main(int argC,char* argV[]){ simError(); } - aLat = strtod( argV[i], &conversionCheck); - if( conversionCheck == argV[i] ) conversionError = true; - if( *conversionCheck != '\0' ) conversionError = true; + aLat = atof( argV[i] ); +// aLat = strtod( argV[i], &conversionCheck); +// if( conversionCheck == argV[i] ) conversionError = true; +// if( *conversionCheck != '\0' ) conversionError = true; if( conversionError ){ sprintf( painCave.errMsg, @@ -374,9 +383,11 @@ int main(int argC,char* argV[]){ simError(); } - bLat = strtod( argV[i], &conversionCheck); - if( conversionCheck == argV[i] ) conversionError = true; - if( *conversionCheck != '\0' ) conversionError = true; + bLat = atof( argV[i] ); + +// bLat = strtod( argV[i], &conversionCheck); +// if( conversionCheck == argV[i] ) conversionError = true; +// if( *conversionCheck != '\0' ) conversionError = true; if( conversionError ){ sprintf( painCave.errMsg, @@ -420,9 +431,11 @@ int main(int argC,char* argV[]){ simError(); } - leafSpacing = strtod( argV[i], &conversionCheck); - if( conversionCheck == argV[i] ) conversionError = true; - if( *conversionCheck != '\0' ) conversionError = true; + leafSpacing = atof( argV[i] ); + +// leafSpacing = strtod( argV[i], &conversionCheck); +// if( conversionCheck == argV[i] ) conversionError = true; +// if( *conversionCheck != '\0' ) conversionError = true; if( conversionError ){ sprintf( painCave.errMsg, @@ -531,20 +544,29 @@ int main(int argC,char* argV[]){ simError(); } - bsInfo.outPrefix = outPrefix; - strcpy(bsInfo.waterName, waterName); - strcpy(bsInfo.lipidName, lipidName); + mainInfo = new SimInfo(); + simInit = new SimSetup(); + simInit->setSimInfo( mainInfo ); + simInit->suspendInit(); + simInit->parseFile( inName ); + simInit->createSim(); - parseBuildBass( inName ); + delete simInit; - buildLatticeBilayer( aLat, bLat, leafSpacing ); + sprintf( mainInfo->statusName, "%s.stat", outPrefix ); + sprintf( mainInfo->sampleName, "%s.dump", outPrefix ); + sprintf( mainInfo->finalName, "%s.init", outPrefix ); + buildLatticeBilayer( aLat, bLat, leafSpacing, waterName, lipidName ); + return 0; } int buildLatticeBilayer(double aLat, double bLat, - double leafSpacing){ + double leafSpacing, + char* waterName, + char* lipidName){ typedef struct{ double rot[3][3]; @@ -570,7 +592,7 @@ int buildLatticeBilayer(double aLat, int nCells, nCellsX, nCellsY, nCellsZ; int nx, ny; double boxX, boxY, boxZ; - double unitVector[3]; + double theta, phi, psi; int which; int targetWaters; @@ -590,7 +612,11 @@ int buildLatticeBilayer(double aLat, int targetNlipids, targetNwaters; double targetWaterLipidRatio; double maxZ, minZ, zHeight; + double maxY, minY; + double maxX, minX; + molStart = NULL; + // create the simInfo objects simnfo = new SimInfo; @@ -599,34 +625,36 @@ int buildLatticeBilayer(double aLat, foundLipid = 0; foundWater = 0; - for(i=0; igetID(), bsInfo.lipidName ) ){ + + for(i=0; inComponents; i++){ + + if( !strcmp( mainInfo->compStamps[i]->getID(), lipidName ) ){ foundLipid = 1; - lipidStamp = bsInfo.compStamps[i]; - targetNlipids = bsInfo.componentsNmol[i]; + lipidStamp = mainInfo->compStamps[i]; + targetNlipids = mainInfo->componentsNmol[i]; lipidNatoms = lipidStamp->getNAtoms(); } - if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){ + if( !strcmp( mainInfo->compStamps[i]->getID(), waterName ) ){ foundWater = 1; - waterStamp = bsInfo.compStamps[i]; - targetNwaters = bsInfo.componentsNmol[i]; + waterStamp = mainInfo->compStamps[i]; + targetNwaters = mainInfo->componentsNmol[i]; waterNatoms = waterStamp->getNAtoms(); } } if( !foundLipid ){ sprintf(painCave.errMsg, - "Could not find lipid \"%s\" in the bass file.\n", - bsInfo.lipidName ); + "latticeBilayer error: Could not find lipid \"%s\" in the bass file.\n", + lipidName ); painCave.isFatal = 1; simError(); } if( !foundWater ){ sprintf(painCave.errMsg, - "Could not find solvent \"%s\" in the bass file.\n", - bsInfo.waterName ); + "latticeBilayer error: Could not find solvent \"%s\" in the bass file.\n", + waterName ); painCave.isFatal = 1; simError(); } @@ -649,11 +677,11 @@ int buildLatticeBilayer(double aLat, testSite.pos[1] = 0.0; testSite.pos[2] = 0.0; - unitVector[0] = 0.0; - unitVector[1] = 0.0; - unitVector[2] = 1.0; - - getUnitRot(unitVector, testSite.rot ); + theta = 0.0; + phi = 0.0; + psi = 0.0; + + getEulerRot(theta, phi, psi, testSite.rot ); lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig ); minZ = 0.0; @@ -677,8 +705,8 @@ int buildLatticeBilayer(double aLat, nLipids = 4 * nx * ny; coord* lipidSites = new coord[nLipids]; - unitVector[0] = 0.0; - unitVector[1] = 0.0; + phi = 0.0; + psi = 0.0; which = 0; @@ -688,21 +716,23 @@ int buildLatticeBilayer(double aLat, lipidSites[which].pos[0] = (double)i * aLat; lipidSites[which].pos[1] = (double)j * bLat; - lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing - maxZ); + lipidSites[which].pos[2] = (2.0* (double)k - 1.0) * + ((leafSpacing / 2.0) - maxZ); - unitVector[2] = 2.0 * (double)k - 1.0; + theta = (1.0 - (double)k) * M_PI; - getUnitRot( unitVector, lipidSites[which].rot ); + getEulerRot( theta, phi, psi, lipidSites[which].rot ); which++; lipidSites[which].pos[0] = aLat * ((double)i + 0.5); lipidSites[which].pos[1] = bLat * ((double)j + 0.5); - lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing - maxZ); + lipidSites[which].pos[2] = (2.0* (double)k - 1.0) * + ((leafSpacing / 2.0) - maxZ); - unitVector[2] = 2.0 * (double)k - 1.0; + theta = (1.0 - (double)k) * M_PI; - getUnitRot( unitVector, lipidSites[which].rot ); + getEulerRot( theta, phi, psi, lipidSites[which].rot ); which++; } @@ -729,7 +759,69 @@ int buildLatticeBilayer(double aLat, } nWaters = nCellsX * nCellsY * nCellsZ * 4; + + // calc current system size; + + nAtoms = 0; + molIndex = 0; + if(molStart != NULL ) delete[] molStart; + molStart = new int[nLipids]; + + for(j=0; jn_atoms = nAtoms; + theConfig = testInfo->getConfiguration(); + theConfig->destroyArrays(); + theConfig->createArrays( nAtoms ); + testInfo->atoms = new Atom*[nAtoms]; + atoms = testInfo->atoms; + + molIndex = 0; + for(i=0; iplaceMol( lipidSites[i].pos, lipidSites[i].rot, atoms, + molStart[molIndex], theConfig ); + molIndex++; + } + + atoms[0]->getPos( myPos ); + + maxX = myPos[0]; + minX = myPos[0]; + + maxY = myPos[1]; + minY = myPos[1]; + + maxZ = myPos[2]; + minZ = myPos[2]; + + for(i=0;igetPos( myPos ); + minX = (minX > myPos[0]) ? myPos[0] : minX; + maxX = (maxX < myPos[0]) ? myPos[0] : maxX; + + minY = (minY > myPos[1]) ? myPos[1] : minY; + maxY = (maxY < myPos[1]) ? myPos[1] : maxY; + + minZ = (minZ > myPos[2]) ? myPos[2] : minZ; + maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ; + } + + boxX = (maxX - minX)+2.0; + boxY = (maxY - minY)+2.0; + boxZ = (maxZ - minZ)+2.0; + + double centerX, centerY, centerZ; + centerX = ((maxX - minX) / 2.0) + minX; + centerY = ((maxY - minY) / 2.0) + minY; + centerZ = ((maxZ - minZ) / 2.0) + minZ; + + // set up water coordinates + coord* waterSites = new coord[nWaters]; waterCell[0] = boxX / nCellsX; @@ -738,10 +830,8 @@ int buildLatticeBilayer(double aLat, Lattice *myORTHO; myORTHO = new Lattice( ORTHORHOMBIC_LATTICE_TYPE, waterCell); - myORTHO->setStartZ( leafSpacing / 2.0 + waterFudge); + myORTHO->setStartZ( maxZ + waterFudge); - boxZ = waterCell[2] * nCellsZ + leafSpacing; - // create an fcc lattice in the water box. which = 0; @@ -754,16 +844,18 @@ int buildLatticeBilayer(double aLat, waterSites[which].pos[0] = posX[l]; waterSites[which].pos[1] = posY[l]; waterSites[which].pos[2] = posZ[l]; + getRandomRot( waterSites[which].rot ); which++; } } } } - // create the real Atom arrays + // calc the system sizes nAtoms = 0; molIndex = 0; + if(molStart != NULL ) delete[] molStart; molStart = new int[nLipids + nWaters]; for(j=0; jn_atoms = nAtoms; + theConfig = testInfo->getConfiguration(); + theConfig->destroyArrays(); + theConfig->createArrays( nAtoms ); + testInfo->atoms = new Atom*[nAtoms]; + atoms = testInfo->atoms; + + molIndex = 0; + for(i=0; iplaceMol( lipidSites[i].pos, lipidSites[i].rot, atoms, + molStart[molIndex], theConfig ); + molIndex++; + } + + for(i=0; iplaceMol( waterSites[i].pos, waterSites[i].rot, atoms, + molStart[molIndex], theConfig ); + molIndex++; + } + + atoms[0]->getPos( myPos ); + maxZ = myPos[2]; + minZ = myPos[2]; + for(i=0;igetPos( myPos ); + minZ = (minZ > myPos[2]) ? myPos[2] : minZ; + maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ; + } + boxZ = (maxZ - (minZ - waterFudge / 2.0)); + + // create the real Atom arrays + theConfig = simnfo->getConfiguration(); theConfig->createArrays( nAtoms ); simnfo->atoms = new Atom*[nAtoms]; + simnfo->n_atoms = nAtoms; atoms = simnfo->atoms; - // wrap back to <0,0,0> as center double Hmat[3][3]; @@ -800,17 +926,26 @@ int buildLatticeBilayer(double aLat, Hmat[2][1] = 0.0; Hmat[2][2] = boxZ; - bsInfo.boxX = boxX; - bsInfo.boxY = boxY; - bsInfo.boxZ = boxZ; - + mainInfo->setBoxM( Hmat ); simnfo->setBoxM( Hmat ); - for(j=0;jwrapVector( lipidSites[j].pos ); + } - for(j=0;jwrapVector( waterSites[j].pos ); + } // initialize lipid positions @@ -824,37 +959,16 @@ int buildLatticeBilayer(double aLat, // initialize the water positions for(i=0; iplaceMol( waterSites[i].pos, waterSites[i].rot, atoms, molStart[molIndex], theConfig ); molIndex++; } - // set up the SimInfo object - - Hmat[0][0] = boxX; - Hmat[0][1] = 0.0; - Hmat[0][2] = 0.0; + strcpy( simnfo->sampleName, mainInfo->sampleName ); + strcpy( simnfo->finalName, mainInfo->finalName ); - Hmat[1][0] = 0.0; - Hmat[1][1] = boxY; - Hmat[1][2] = 0.0; - - Hmat[2][0] = 0.0; - Hmat[2][1] = 0.0; - Hmat[2][2] = boxZ; - - - bsInfo.boxX = boxX; - bsInfo.boxY = boxY; - bsInfo.boxZ = boxZ; - - simnfo->setBoxM( Hmat ); - - sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix ); - sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix ); - // set up the writer and write out writer = new DumpWriter( simnfo );