ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/utils/sysbuilder/latticeBilayer.cpp
(Generate patch)

Comparing trunk/OOPSE/utils/sysbuilder/latticeBilayer.cpp (file contents):
Revision 821 by mmeineke, Fri Oct 24 22:17:45 2003 UTC vs.
Revision 874 by mmeineke, Fri Nov 21 20:10:02 2003 UTC

# Line 7 | Line 7
7   #include "simError.h"
8   #include "SimInfo.hpp"
9   #include "ReadWrite.hpp"
10 + #include "SimSetup.hpp"
11  
12   #include "MoLocator.hpp"
13   #include "latticeBuilder.hpp"
13 – #include "QuickBass.hpp"
14  
15 +
16   #define VERSION_MAJOR 0
17   #define VERSION_MINOR 1
18  
# Line 19 | Line 20 | int buildLatticeBilayer( double aLat,
20   void usage(void);
21   int buildLatticeBilayer( double aLat,
22                           double bLat,
23 <                         double leafSpacing);
23 >                         double leafSpacing,
24 >                         char* waterName,
25 >                         char* lipidName);
26   using namespace std;
27 + SimInfo* mainInfo;
28  
29   int main(int argC,char* argV[]){
30    
# Line 45 | Line 49 | int main(int argC,char* argV[]){
49    
50    char* inName;
51    
52 +  SimSetup* simInit;
53 +  
54    // first things first, all of the initializations
55  
56    fflush(stdout);
# Line 280 | Line 286 | int main(int argC,char* argV[]){
286                simError();
287              }    
288  
289 <            bLat = strtod( argV[i], &conversionCheck);
290 <            if( conversionCheck == argV[i] ) conversionError = true;
291 <            if( *conversionCheck != '\0' ) conversionError = true;
289 >
290 >            bLat = atof( argV[i] );
291 > //          bLat = strtod( argV[i], &conversionCheck);
292 > //          if( conversionCheck == argV[i] ) conversionError = true;
293 > //          if( *conversionCheck != '\0' ) conversionError = true;
294              
295              if( conversionError ){
296                sprintf( painCave.errMsg,
# Line 328 | Line 336 | int main(int argC,char* argV[]){
336                simError();
337              }    
338  
339 <            aLat = strtod( argV[i], &conversionCheck);
340 <            if( conversionCheck == argV[i] ) conversionError = true;
341 <            if( *conversionCheck != '\0' ) conversionError = true;
339 >            aLat = atof( argV[i] );
340 > //          aLat = strtod( argV[i], &conversionCheck);
341 > //          if( conversionCheck == argV[i] ) conversionError = true;
342 > //          if( *conversionCheck != '\0' ) conversionError = true;
343              
344              if( conversionError ){
345                sprintf( painCave.errMsg,
# Line 374 | Line 383 | int main(int argC,char* argV[]){
383                simError();
384              }    
385  
386 <            bLat = strtod( argV[i], &conversionCheck);
387 <            if( conversionCheck == argV[i] ) conversionError = true;
388 <            if( *conversionCheck != '\0' ) conversionError = true;
386 >            bLat = atof( argV[i] );
387 >
388 > //          bLat = strtod( argV[i], &conversionCheck);
389 > //          if( conversionCheck == argV[i] ) conversionError = true;
390 > //          if( *conversionCheck != '\0' ) conversionError = true;
391              
392              if( conversionError ){
393                sprintf( painCave.errMsg,
# Line 420 | Line 431 | int main(int argC,char* argV[]){
431                simError();
432              }    
433  
434 <            leafSpacing = strtod( argV[i], &conversionCheck);
435 <            if( conversionCheck == argV[i] ) conversionError = true;
436 <            if( *conversionCheck != '\0' ) conversionError = true;
434 >            leafSpacing = atof( argV[i] );
435 >
436 > //          leafSpacing = strtod( argV[i], &conversionCheck);
437 > //          if( conversionCheck == argV[i] ) conversionError = true;
438 > //          if( *conversionCheck != '\0' ) conversionError = true;
439              
440              if( conversionError ){
441                sprintf( painCave.errMsg,
# Line 531 | Line 544 | int main(int argC,char* argV[]){
544      simError();
545    }
546  
547 <  bsInfo.outPrefix = outPrefix;
548 <  strcpy(bsInfo.waterName, waterName);
549 <  strcpy(bsInfo.lipidName, lipidName);
547 >  mainInfo = new SimInfo();
548 >  simInit = new SimSetup();
549 >  simInit->setSimInfo( mainInfo );
550 >  simInit->suspendInit();
551 >  simInit->parseFile( inName );
552 >  simInit->createSim();
553  
554 <  parseBuildBass( inName );
554 >  delete simInit;
555  
556 <  buildLatticeBilayer( aLat, bLat, leafSpacing );
556 >  sprintf( mainInfo->statusName, "%s.stat", outPrefix );
557 >  sprintf( mainInfo->sampleName, "%s.dump", outPrefix );
558 >  sprintf( mainInfo->finalName, "%s.init", outPrefix );
559  
560 +  buildLatticeBilayer( aLat, bLat, leafSpacing, waterName, lipidName );
561 +
562    return 0;  
563   }
564  
565   int buildLatticeBilayer(double aLat,
566                          double bLat,
567 <                        double leafSpacing){
567 >                        double leafSpacing,
568 >                        char* waterName,
569 >                        char* lipidName){
570  
571    typedef struct{
572      double rot[3][3];
# Line 559 | Line 581 | int buildLatticeBilayer(double aLat,
581    double *posX, *posY, *posZ;
582    double pos[3], posA[3], posB[3];
583  
584 <  const double waterFudge = 5.0;
584 >  const double waterFudge = 6.0;
585  
586    int i,j,k,l;
587    int nAtoms, atomIndex, molIndex, molID;
# Line 570 | Line 592 | int buildLatticeBilayer(double aLat,
592    int nCells, nCellsX, nCellsY, nCellsZ;
593    int nx, ny;
594    double boxX, boxY, boxZ;
595 <  double unitVector[3];
595 >  double theta, phi, psi;
596    int which;
597    int targetWaters;
598    
599    Atom** atoms;
578 –  SimInfo* simnfo;
600    SimInfo* testInfo;
601    coord testSite;
602    SimState* theConfig;
# Line 590 | Line 611 | int buildLatticeBilayer(double aLat,
611    int targetNlipids, targetNwaters;
612    double targetWaterLipidRatio;
613    double maxZ, minZ, zHeight;
614 +  double maxY, minY;
615 +  double maxX, minX;
616    
617 <  // create the simInfo objects
617 >  molStart = NULL;
618  
596 –  simnfo = new SimInfo;
597 –
619    // set the the lipidStamp
620  
621    foundLipid = 0;
622    foundWater = 0;
623 <  for(i=0; i<bsInfo.nComponents; i++){
624 <    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
623 >
624 >  for(i=0; i<mainInfo->nComponents; i++){
625 >    
626 >    if( !strcmp( mainInfo->compStamps[i]->getID(), lipidName ) ){
627            
628        foundLipid = 1;
629 <      lipidStamp = bsInfo.compStamps[i];
630 <      targetNlipids = bsInfo.componentsNmol[i];
629 >      lipidStamp = mainInfo->compStamps[i];
630 >      targetNlipids = mainInfo->componentsNmol[i];
631        lipidNatoms = lipidStamp->getNAtoms();
632      }
633 <    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
633 >    if( !strcmp( mainInfo->compStamps[i]->getID(), waterName ) ){
634            
635        foundWater = 1;
636            
637 <      waterStamp = bsInfo.compStamps[i];
638 <      targetNwaters = bsInfo.componentsNmol[i];
637 >      waterStamp = mainInfo->compStamps[i];
638 >      targetNwaters = mainInfo->componentsNmol[i];
639        waterNatoms = waterStamp->getNAtoms();
640      }
641    }
642    if( !foundLipid ){
643      sprintf(painCave.errMsg,
644 <            "Could not find lipid \"%s\" in the bass file.\n",
645 <            bsInfo.lipidName );
644 >            "latticeBilayer error: Could not find lipid \"%s\" in the bass file.\n",
645 >            lipidName );
646      painCave.isFatal = 1;
647      simError();
648    }
649    if( !foundWater ){
650      sprintf(painCave.errMsg,
651 <            "Could not find solvent \"%s\" in the bass file.\n",
652 <            bsInfo.waterName );
651 >            "latticeBilayer error: Could not find solvent \"%s\" in the bass file.\n",
652 >            waterName );
653      painCave.isFatal = 1;
654      simError();
655    }
# Line 649 | Line 672 | int buildLatticeBilayer(double aLat,
672    testSite.pos[1] = 0.0;
673    testSite.pos[2] = 0.0;
674    
675 <  unitVector[0] = 0.0;
676 <  unitVector[1] = 0.0;
677 <  unitVector[2] = 1.0;
678 <  
679 <  getUnitRot(unitVector, testSite.rot );
675 >  theta = 0.0;
676 >  phi = 0.0;
677 >  psi = 0.0;
678 >
679 >  getEulerRot(theta, phi, psi, testSite.rot );
680    lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig );
681  
682    minZ = 0.0;
# Line 666 | Line 689 | int buildLatticeBilayer(double aLat,
689    }
690    zHeight = maxZ - minZ;
691      
692 <  nCells = (int) sqrt( (double)targetNlipids * bLat / (4.0 * aLat) );
692 >  std::cerr << "aLat = " << aLat << "; bLat = " << bLat << "\n";
693  
694 +  nCells = (int)ceil( sqrt( (double)targetNlipids * bLat / (4.0 * aLat) ));
695 +
696    nx = nCells;
697    ny = (int) ((double)nCells  * aLat / bLat);
698  
699 +  std::cerr << "nx = " << nx << "; ny = " << ny << "\n";
700 +
701    boxX = nx * aLat;
702    boxY = ny * bLat;  
703  
704    nLipids = 4 * nx * ny;
705    coord* lipidSites = new coord[nLipids];
706  
707 <  unitVector[0] = 0.0;
708 <  unitVector[1] = 0.0;
707 >  phi = 0.0;
708 >  psi = 0.0;
709  
710    which = 0;
711  
# Line 688 | Line 715 | int buildLatticeBilayer(double aLat,
715          
716          lipidSites[which].pos[0] = (double)i * aLat;
717          lipidSites[which].pos[1] = (double)j * bLat;
718 <        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing - maxZ);
718 >        lipidSites[which].pos[2] = (2.0* (double)k - 1.0) *
719 >          ((leafSpacing / 2.0) - maxZ);
720          
721 <        unitVector[2] = 2.0 * (double)k  - 1.0;
721 >        theta = (1.0 - (double)k) * M_PI;
722          
723 <        getUnitRot( unitVector, lipidSites[which].rot );
723 >        getEulerRot( theta, phi, psi, lipidSites[which].rot );
724  
725          which++;
726  
727          lipidSites[which].pos[0] = aLat * ((double)i + 0.5);
728          lipidSites[which].pos[1] = bLat * ((double)j + 0.5);
729 <        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing - maxZ);
729 >        lipidSites[which].pos[2] = (2.0* (double)k - 1.0) *
730 >          ((leafSpacing / 2.0) - maxZ);
731          
732 <        unitVector[2] = 2.0 * (double)k  - 1.0;
732 >        theta = (1.0 - (double)k) * M_PI;
733          
734 <        getUnitRot( unitVector, lipidSites[which].rot );
734 >        getEulerRot( theta, phi, psi, lipidSites[which].rot );
735          
736          which++;
737        }
# Line 730 | Line 759 | int buildLatticeBilayer(double aLat,
759    
760    nWaters = nCellsX * nCellsY * nCellsZ * 4;
761    
762 +  // calc current system size;
763 +
764 +  nAtoms = 0;
765 +  molIndex = 0;
766 +  if(molStart != NULL ) delete[] molStart;
767 +  molStart = new int[nLipids];  
768 +  
769 +  for(j=0; j<nLipids; j++){
770 +    molStart[molIndex] = nAtoms;
771 +    molIndex++;
772 +    nAtoms += lipidNatoms;
773 +  }
774 +
775 +  testInfo->n_atoms = nAtoms;
776 +  theConfig = testInfo->getConfiguration();
777 +  theConfig->destroyArrays();
778 +  theConfig->createArrays( nAtoms );
779 +  testInfo->atoms = new Atom*[nAtoms];
780 +  atoms = testInfo->atoms;
781 +
782 +  molIndex = 0;
783 +  for(i=0; i<nLipids; i++ ){
784 +    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
785 +                           molStart[molIndex], theConfig );
786 +    molIndex++;
787 +  }
788 +
789 +  atoms[0]->getPos( myPos );
790 +
791 +  maxX = myPos[0];
792 +  minX = myPos[0];
793 +
794 +  maxY = myPos[1];
795 +  minY = myPos[1];
796 +
797 +  maxZ = myPos[2];
798 +  minZ = myPos[2];
799 +
800 +  for(i=0;i<nAtoms;i++){
801 +    atoms[i]->getPos( myPos );
802 +    minX = (minX > myPos[0]) ? myPos[0] : minX;
803 +    maxX = (maxX < myPos[0]) ? myPos[0] : maxX;
804 +
805 +    minY = (minY > myPos[1]) ? myPos[1] : minY;
806 +    maxY = (maxY < myPos[1]) ? myPos[1] : maxY;
807 +
808 +    minZ = (minZ > myPos[2]) ? myPos[2] : minZ;
809 +    maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ;
810 +  }
811 +
812 +  boxX = (maxX - minX)+2.0;
813 +  boxY = (maxY - minY)+2.0;
814 +  boxZ = (maxZ - minZ)+2.0;
815 +
816 +  double centerX, centerY, centerZ;
817 +  
818 +  centerX = ((maxX - minX) / 2.0) + minX;
819 +  centerY = ((maxY - minY) / 2.0) + minY;
820 +  centerZ = ((maxZ - minZ) / 2.0) + minZ;
821 +
822 +  // set up water coordinates
823 +
824    coord* waterSites = new coord[nWaters];
825  
826    waterCell[0] = boxX / nCellsX;
# Line 738 | Line 829 | int buildLatticeBilayer(double aLat,
829  
830    Lattice *myORTHO;
831    myORTHO = new Lattice( ORTHORHOMBIC_LATTICE_TYPE, waterCell);
832 <  myORTHO->setStartZ( leafSpacing / 2.0 + waterFudge);
832 >  myORTHO->setStartZ( maxZ + waterFudge);
833  
743 –  boxZ = waterCell[2] * nCellsZ + leafSpacing;
744 –        
834    // create an fcc lattice in the water box.
835    
836    which = 0;
# Line 754 | Line 843 | int buildLatticeBilayer(double aLat,
843            waterSites[which].pos[0] = posX[l];
844            waterSites[which].pos[1] = posY[l];
845            waterSites[which].pos[2] = posZ[l];
846 +          getRandomRot( waterSites[which].rot );
847            which++;
848          }
849        }
850      }
851    }  
852  
853 <  // create the real Atom arrays
853 >  // calc the system sizes
854    
855    nAtoms = 0;
856    molIndex = 0;
857 +  if(molStart != NULL ) delete[] molStart;
858    molStart = new int[nLipids + nWaters];  
859    
860    for(j=0; j<nLipids; j++){
# Line 778 | Line 869 | int buildLatticeBilayer(double aLat,
869      nAtoms += waterNatoms;
870    }
871    
872 <  theConfig = simnfo->getConfiguration();
872 >  testInfo->n_atoms = nAtoms;
873 >  theConfig = testInfo->getConfiguration();
874 >  theConfig->destroyArrays();
875    theConfig->createArrays( nAtoms );
876 <  simnfo->atoms = new Atom*[nAtoms];
877 <  atoms = simnfo->atoms;
876 >  testInfo->atoms = new Atom*[nAtoms];
877 >  atoms = testInfo->atoms;
878  
879 +  molIndex = 0;
880 +  for(i=0; i<nLipids; i++ ){
881 +    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
882 +                           molStart[molIndex], theConfig );
883 +    molIndex++;
884 +  }
885  
886 +  for(i=0; i<nWaters; i++){
887 +
888 +    getRandomRot( waterSites[i].rot );
889 +    waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
890 +                           molStart[molIndex], theConfig );
891 +    molIndex++;
892 +  }
893 +
894 +  atoms[0]->getPos( myPos );
895 +  maxZ = myPos[2];
896 +  minZ = myPos[2];
897 +  for(i=0;i<nAtoms;i++){
898 +    atoms[i]->getPos( myPos );
899 +    minZ = (minZ > myPos[2]) ? myPos[2] : minZ;
900 +    maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ;
901 +  }
902 +  boxZ = (maxZ - (minZ - waterFudge / 2.0));
903 +  
904 +  // create the real Atom arrays
905 +
906 +  delete[] (mainInfo->atoms);
907 +  theConfig = mainInfo->getConfiguration();
908 +  theConfig->createArrays( nAtoms );
909 +  mainInfo->atoms = new Atom*[nAtoms];
910 +  mainInfo->n_atoms = nAtoms;
911 +  atoms = mainInfo->atoms;
912 +
913    // wrap back to <0,0,0> as center
914  
915    double Hmat[3][3];
# Line 800 | Line 926 | int buildLatticeBilayer(double aLat,
926    Hmat[2][1] = 0.0;
927    Hmat[2][2] = boxZ;
928    
929 <  bsInfo.boxX = boxX;
804 <  bsInfo.boxY = boxY;
805 <  bsInfo.boxZ = boxZ;
806 <  
807 <  simnfo->setBoxM( Hmat );
929 >  mainInfo->setBoxM( Hmat );
930  
931 <  for(j=0;j<nLipids;j++)
810 <    simnfo->wrapVector( lipidSites[j].pos );
931 >  for(j=0;j<nLipids;j++){
932  
933 <  for(j=0;j<nWaters;j++)
934 <    simnfo->wrapVector( waterSites[j].pos );
933 >    lipidSites[j].pos[0] -= centerX;
934 >    lipidSites[j].pos[1] -= centerY;
935 >    lipidSites[j].pos[2] -= centerZ;
936  
937 +    mainInfo->wrapVector( lipidSites[j].pos );
938 +  }
939 +
940 +  for(j=0;j<nWaters;j++){
941 +
942 +    waterSites[j].pos[0] -= centerX;
943 +    waterSites[j].pos[1] -= centerY;
944 +    waterSites[j].pos[2] -= centerZ;
945 +
946 +    mainInfo->wrapVector( waterSites[j].pos );
947 +  }
948 +
949    // initialize lipid positions
950  
951    molIndex = 0;
# Line 824 | Line 958 | int buildLatticeBilayer(double aLat,
958    // initialize the water positions
959  
960    for(i=0; i<nWaters; i++){
961 <    
961 >
962      getRandomRot( waterSites[i].rot );
963      waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
964                             molStart[molIndex], theConfig );
965      molIndex++;
966    }
967  
834 –  // set up the SimInfo object
835 –  
836 –  Hmat[0][0] = boxX;
837 –  Hmat[0][1] = 0.0;
838 –  Hmat[0][2] = 0.0;
839 –
840 –  Hmat[1][0] = 0.0;
841 –  Hmat[1][1] = boxY;
842 –  Hmat[1][2] = 0.0;
843 –
844 –  Hmat[2][0] = 0.0;
845 –  Hmat[2][1] = 0.0;
846 –  Hmat[2][2] = boxZ;
847 –  
848 –
849 –  bsInfo.boxX = boxX;
850 –  bsInfo.boxY = boxY;
851 –  bsInfo.boxZ = boxZ;
852 –  
853 –  simnfo->setBoxM( Hmat );
854 –
855 –  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
856 –  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
857 –
968    // set up the writer and write out
969    
970 <  writer = new DumpWriter( simnfo );
970 >  writer = new DumpWriter( mainInfo );
971    writer->writeFinal( 0.0 );
972  
973 +  std::cout << "\n"
974 +            << "----------------------------------------\n"
975 +            << "\n"
976 +            << " final nLipids = " << nLipids << "\n"
977 +            << " final nWaters = " << nWaters << "\n"
978 +            << "\n";
979 +      
980    return 1;
981   }
982  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines