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 851 by mmeineke, Wed Nov 5 19:18:17 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 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    
# Line 590 | Line 612 | int buildLatticeBilayer(double aLat,
612    int targetNlipids, targetNwaters;
613    double targetWaterLipidRatio;
614    double maxZ, minZ, zHeight;
615 +  double maxY, minY;
616 +  double maxX, minX;
617    
618 +  molStart = NULL;
619 +
620    // create the simInfo objects
621  
622    simnfo = new SimInfo;
# Line 599 | Line 625 | int buildLatticeBilayer(double aLat,
625  
626    foundLipid = 0;
627    foundWater = 0;
628 <  for(i=0; i<bsInfo.nComponents; i++){
629 <    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
628 >
629 >  for(i=0; i<mainInfo->nComponents; i++){
630 >    
631 >    if( !strcmp( mainInfo->compStamps[i]->getID(), lipidName ) ){
632            
633        foundLipid = 1;
634 <      lipidStamp = bsInfo.compStamps[i];
635 <      targetNlipids = bsInfo.componentsNmol[i];
634 >      lipidStamp = mainInfo->compStamps[i];
635 >      targetNlipids = mainInfo->componentsNmol[i];
636        lipidNatoms = lipidStamp->getNAtoms();
637      }
638 <    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
638 >    if( !strcmp( mainInfo->compStamps[i]->getID(), waterName ) ){
639            
640        foundWater = 1;
641            
642 <      waterStamp = bsInfo.compStamps[i];
643 <      targetNwaters = bsInfo.componentsNmol[i];
642 >      waterStamp = mainInfo->compStamps[i];
643 >      targetNwaters = mainInfo->componentsNmol[i];
644        waterNatoms = waterStamp->getNAtoms();
645      }
646    }
647    if( !foundLipid ){
648      sprintf(painCave.errMsg,
649 <            "Could not find lipid \"%s\" in the bass file.\n",
650 <            bsInfo.lipidName );
649 >            "latticeBilayer error: Could not find lipid \"%s\" in the bass file.\n",
650 >            lipidName );
651      painCave.isFatal = 1;
652      simError();
653    }
654    if( !foundWater ){
655      sprintf(painCave.errMsg,
656 <            "Could not find solvent \"%s\" in the bass file.\n",
657 <            bsInfo.waterName );
656 >            "latticeBilayer error: Could not find solvent \"%s\" in the bass file.\n",
657 >            waterName );
658      painCave.isFatal = 1;
659      simError();
660    }
# Line 649 | Line 677 | int buildLatticeBilayer(double aLat,
677    testSite.pos[1] = 0.0;
678    testSite.pos[2] = 0.0;
679    
680 <  unitVector[0] = 0.0;
681 <  unitVector[1] = 0.0;
682 <  unitVector[2] = 1.0;
683 <  
684 <  getUnitRot(unitVector, testSite.rot );
680 >  theta = 0.0;
681 >  phi = 0.0;
682 >  psi = 0.0;
683 >
684 >  getEulerRot(theta, phi, psi, testSite.rot );
685    lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig );
686  
687    minZ = 0.0;
# Line 677 | Line 705 | int buildLatticeBilayer(double aLat,
705    nLipids = 4 * nx * ny;
706    coord* lipidSites = new coord[nLipids];
707  
708 <  unitVector[0] = 0.0;
709 <  unitVector[1] = 0.0;
708 >  phi = 0.0;
709 >  psi = 0.0;
710  
711    which = 0;
712  
# Line 688 | Line 716 | int buildLatticeBilayer(double aLat,
716          
717          lipidSites[which].pos[0] = (double)i * aLat;
718          lipidSites[which].pos[1] = (double)j * bLat;
719 <        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing - maxZ);
719 >        lipidSites[which].pos[2] = (2.0* (double)k - 1.0) *
720 >          ((leafSpacing / 2.0) - maxZ);
721          
722 <        unitVector[2] = 2.0 * (double)k  - 1.0;
722 >        theta = (1.0 - (double)k) * M_PI;
723          
724 <        getUnitRot( unitVector, lipidSites[which].rot );
724 >        getEulerRot( theta, phi, psi, lipidSites[which].rot );
725  
726          which++;
727  
728          lipidSites[which].pos[0] = aLat * ((double)i + 0.5);
729          lipidSites[which].pos[1] = bLat * ((double)j + 0.5);
730 <        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing - maxZ);
730 >        lipidSites[which].pos[2] = (2.0* (double)k - 1.0) *
731 >          ((leafSpacing / 2.0) - maxZ);
732          
733 <        unitVector[2] = 2.0 * (double)k  - 1.0;
733 >        theta = (1.0 - (double)k) * M_PI;
734          
735 <        getUnitRot( unitVector, lipidSites[which].rot );
735 >        getEulerRot( theta, phi, psi, lipidSites[which].rot );
736          
737          which++;
738        }
# Line 730 | Line 760 | int buildLatticeBilayer(double aLat,
760    
761    nWaters = nCellsX * nCellsY * nCellsZ * 4;
762    
763 +  // calc current system size;
764 +
765 +  nAtoms = 0;
766 +  molIndex = 0;
767 +  if(molStart != NULL ) delete[] molStart;
768 +  molStart = new int[nLipids];  
769 +  
770 +  for(j=0; j<nLipids; j++){
771 +    molStart[molIndex] = nAtoms;
772 +    molIndex++;
773 +    nAtoms += lipidNatoms;
774 +  }
775 +
776 +  testInfo->n_atoms = nAtoms;
777 +  theConfig = testInfo->getConfiguration();
778 +  theConfig->destroyArrays();
779 +  theConfig->createArrays( nAtoms );
780 +  testInfo->atoms = new Atom*[nAtoms];
781 +  atoms = testInfo->atoms;
782 +
783 +  molIndex = 0;
784 +  for(i=0; i<nLipids; i++ ){
785 +    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
786 +                           molStart[molIndex], theConfig );
787 +    molIndex++;
788 +  }
789 +
790 +  atoms[0]->getPos( myPos );
791 +
792 +  maxX = myPos[0];
793 +  minX = myPos[0];
794 +
795 +  maxY = myPos[1];
796 +  minY = myPos[1];
797 +
798 +  maxZ = myPos[2];
799 +  minZ = myPos[2];
800 +
801 +  for(i=0;i<nAtoms;i++){
802 +    atoms[i]->getPos( myPos );
803 +    minX = (minX > myPos[0]) ? myPos[0] : minX;
804 +    maxX = (maxX < myPos[0]) ? myPos[0] : maxX;
805 +
806 +    minY = (minY > myPos[1]) ? myPos[1] : minY;
807 +    maxY = (maxY < myPos[1]) ? myPos[1] : maxY;
808 +
809 +    minZ = (minZ > myPos[2]) ? myPos[2] : minZ;
810 +    maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ;
811 +  }
812 +
813 +  boxX = (maxX - minX)+2.0;
814 +  boxY = (maxY - minY)+2.0;
815 +  boxZ = (maxZ - minZ)+2.0;
816 +
817 +  double centerX, centerY, centerZ;
818 +  
819 +  centerX = ((maxX - minX) / 2.0) + minX;
820 +  centerY = ((maxY - minY) / 2.0) + minY;
821 +  centerZ = ((maxZ - minZ) / 2.0) + minZ;
822 +
823 +  // set up water coordinates
824 +
825    coord* waterSites = new coord[nWaters];
826  
827    waterCell[0] = boxX / nCellsX;
# Line 738 | Line 830 | int buildLatticeBilayer(double aLat,
830  
831    Lattice *myORTHO;
832    myORTHO = new Lattice( ORTHORHOMBIC_LATTICE_TYPE, waterCell);
833 <  myORTHO->setStartZ( leafSpacing / 2.0 + waterFudge);
833 >  myORTHO->setStartZ( maxZ + waterFudge);
834  
743 –  boxZ = waterCell[2] * nCellsZ + leafSpacing;
744 –        
835    // create an fcc lattice in the water box.
836    
837    which = 0;
# Line 754 | Line 844 | int buildLatticeBilayer(double aLat,
844            waterSites[which].pos[0] = posX[l];
845            waterSites[which].pos[1] = posY[l];
846            waterSites[which].pos[2] = posZ[l];
847 +          getRandomRot( waterSites[which].rot );
848            which++;
849          }
850        }
851      }
852    }  
853  
854 <  // create the real Atom arrays
854 >  // calc the system sizes
855    
856    nAtoms = 0;
857    molIndex = 0;
858 +  if(molStart != NULL ) delete[] molStart;
859    molStart = new int[nLipids + nWaters];  
860    
861    for(j=0; j<nLipids; j++){
# Line 778 | Line 870 | int buildLatticeBilayer(double aLat,
870      nAtoms += waterNatoms;
871    }
872    
873 +  testInfo->n_atoms = nAtoms;
874 +  theConfig = testInfo->getConfiguration();
875 +  theConfig->destroyArrays();
876 +  theConfig->createArrays( nAtoms );
877 +  testInfo->atoms = new Atom*[nAtoms];
878 +  atoms = testInfo->atoms;
879 +
880 +  molIndex = 0;
881 +  for(i=0; i<nLipids; i++ ){
882 +    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
883 +                           molStart[molIndex], theConfig );
884 +    molIndex++;
885 +  }
886 +
887 +  for(i=0; i<nWaters; i++){
888 +
889 +    getRandomRot( waterSites[i].rot );
890 +    waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
891 +                           molStart[molIndex], theConfig );
892 +    molIndex++;
893 +  }
894 +
895 +  atoms[0]->getPos( myPos );
896 +  maxZ = myPos[2];
897 +  minZ = myPos[2];
898 +  for(i=0;i<nAtoms;i++){
899 +    atoms[i]->getPos( myPos );
900 +    minZ = (minZ > myPos[2]) ? myPos[2] : minZ;
901 +    maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ;
902 +  }
903 +  boxZ = (maxZ - (minZ - waterFudge / 2.0));
904 +  
905 +  // create the real Atom arrays
906 +
907    theConfig = simnfo->getConfiguration();
908    theConfig->createArrays( nAtoms );
909    simnfo->atoms = new Atom*[nAtoms];
910 +  simnfo->n_atoms = nAtoms;
911    atoms = simnfo->atoms;
912  
786 –
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 <  
929 >  mainInfo->setBoxM( Hmat );
930    simnfo->setBoxM( Hmat );
931  
932 <  for(j=0;j<nLipids;j++)
932 >  for(j=0;j<nLipids;j++){
933 >
934 >    lipidSites[j].pos[0] -= centerX;
935 >    lipidSites[j].pos[1] -= centerY;
936 >    lipidSites[j].pos[2] -= centerZ;
937 >
938      simnfo->wrapVector( lipidSites[j].pos );
939 +  }
940  
941 <  for(j=0;j<nWaters;j++)
941 >  for(j=0;j<nWaters;j++){
942 >
943 >    waterSites[j].pos[0] -= centerX;
944 >    waterSites[j].pos[1] -= centerY;
945 >    waterSites[j].pos[2] -= centerZ;
946 >
947      simnfo->wrapVector( waterSites[j].pos );
948 +  }
949  
950    // initialize lipid positions
951  
# Line 824 | Line 959 | int buildLatticeBilayer(double aLat,
959    // initialize the water positions
960  
961    for(i=0; i<nWaters; i++){
962 <    
962 >
963      getRandomRot( waterSites[i].rot );
964      waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
965                             molStart[molIndex], theConfig );
966      molIndex++;
967    }
968  
969 <  // set up the SimInfo object
970 <  
836 <  Hmat[0][0] = boxX;
837 <  Hmat[0][1] = 0.0;
838 <  Hmat[0][2] = 0.0;
969 >  strcpy( simnfo->sampleName, mainInfo->sampleName );
970 >  strcpy( simnfo->finalName, mainInfo->finalName );
971  
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 –
972    // set up the writer and write out
973    
974    writer = new DumpWriter( simnfo );
975    writer->writeFinal( 0.0 );
976  
977 +  std::cout << "\n"
978 +            << "----------------------------------------\n"
979 +            << "\n"
980 +            << " final nLipids = " << nLipids << "\n"
981 +            << " final nWaters = " << nWaters << "\n"
982 +            << "\n";
983 +      
984    return 1;
985   }
986  

Diff Legend

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