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

Comparing trunk/OOPSE/utils/sysbuilder/bilayerSys.cpp (file contents):
Revision 820 by gezelter, Fri Oct 24 17:36:18 2003 UTC vs.
Revision 821 by mmeineke, Fri Oct 24 22:17:45 2003 UTC

# Line 33 | Line 33 | int buildRandomBilayer( void );
33                 double boxX, double boxY, double boxZ );
34  
35   int buildRandomBilayer( void );
36 – int buildLatticeBilayer( int isHexLattice,
37 –                         double hexSpacing,
38 –                         double aLat,
39 –                         double bLat,
40 –                         int targetNlipid,
41 –                         double targetWaterLipidRatio,
42 –                         double leafSpacing);
36  
44 – void getRandomRot( double rot[3][3] );
45 – void getEulerRot( double theta, double phi, double psi, double rot[3][3] );
46 – void getUnitRot( double unit[3], double rot[3][3] );
47 –
37   int buildBilayer( int isRandom ){
38    
39    if( isRandom ){
# Line 52 | Line 41 | int buildBilayer( int isRandom ){
41    }
42    else{
43  
44 <    
45 <
57 <    return buildLatticeBilayer();
44 >    std::cerr << "unsupported feature\n";
45 >    return 0;
46    }
47   }
48  
# Line 579 | Line 567 | int buildRandomBilayer( void ){
567    //     if( molSeq != NULL ) delete[] molSeq;
568    //     if( simnfo != NULL ) delete simnfo;
569    //     if( writer != NULL ) delete writer;
582 –
583 –  return 1;
584 – }
585 –
586 – int buildLatticeBilayer(int isHexLattice,
587 –                        double hexSpacing,
588 –                        double aLat,
589 –                        double bLat,
590 –                        int targetNlipid,
591 –                        double targetWaterLipidRatio,
592 –                        double leafSpacing){
593 –
594 –  typedef struct{
595 –    double rot[3][3];
596 –    double pos[3];
597 –  } coord;
598 –
599 –  const double waterRho = 0.0334; // number density per cubic angstrom
600 –  const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters
601 –  
602 –  double waterCell[3];
603 –
604 –  double *posX, *posY, *posZ;
605 –  double pos[3], posA[3], posB[3];
606 –
607 –  const double waterFudge = 5.0;
608 –
609 –  int i,j,k,l;
610 –  int nAtoms, atomIndex, molIndex, molID;
611 –  int* molSeq;
612 –  int* molMap;
613 –  int* molStart;
614 –  int testTot, done;
615 –  int nCells, nCellsX, nCellsY, nCellsZ;
616 –  int nx, ny;
617 –  double boxX, boxY, boxZ;
618 –  double unitVector[3];
619 –  int which;
620 –  int targetWaters;
621 –  
622 –  
623 –
624 –  coord testSite;
625 –
626 –  Atom** atoms;
627 –  SimInfo* simnfo;
628 –  SimState* theConfig;
629 –  DumpWriter* writer;
630 –
631 –  MoleculeStamp* lipidStamp;
632 –  MoleculeStamp* waterStamp;  
633 –  MoLocator *lipidLocate;
634 –  MoLocator *waterLocate;
635 –  int foundLipid, foundWater;
636 –  int nLipids, lipidNatoms, nWaters, waterNatoms;
637 –  
638 –  srand48( RAND_SEED );
639 –
640 –  // create the simInfo objects
641 –
642 –  simnfo = new SimInfo;
643 –
644 –  // set the the lipidStamp
645 –
646 –  foundLipid = 0;
647 –  foundWater = 0;
648 –  for(i=0; i<bsInfo.nComponents; i++){
649 –    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
650 –          
651 –      foundLipid = 1;
652 –      lipidStamp = bsInfo.compStamps[i];
653 –      nLipids = bsInfo.componentsNmol[i];
654 –      lipidNatoms = lipidStamp->getNAtoms();
655 –    }
656 –    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
657 –          
658 –      foundWater = 1;
659 –          
660 –      waterStamp = bsInfo.compStamps[i];
661 –      nWaters = bsInfo.componentsNmol[i];
662 –      waterNatoms = waterStamp->getNAtoms();
663 –    }
664 –  }
665 –  if( !foundLipid ){
666 –    sprintf(painCave.errMsg,
667 –            "Could not find lipid \"%s\" in the bass file.\n",
668 –            bsInfo.lipidName );
669 –    painCave.isFatal = 1;
670 –    simError();
671 –  }
672 –  if( !foundWater ){
673 –    sprintf(painCave.errMsg,
674 –            "Could not find solvent \"%s\" in the bass file.\n",
675 –            bsInfo.waterName );
676 –    painCave.isFatal = 1;
677 –    simError();
678 –  }
679 –  
680 –  //create the Molocator arrays
681 –  
682 –  lipidLocate = new MoLocator( lipidStamp );
683 –  waterLocate = new MoLocator( waterStamp );
684 –
685 –
686 –  // set up the bilayer leaves
687 –  
688 –  if (isHexLattice) {
689 –    aLat = sqrt(3.0)*hexSpacing;
690 –    bLat = hexSpacing;
691 –  }
692 –
693 –  nCells = (int) sqrt( (double)targetNlipid * bLat / (4.0 * aLat) );
694 –
695 –  nx = nCells;
696 –  ny = (int) ((double)nCells  * aLat / bLat);
697 –
698 –  boxX = nx * aLat;
699 –  boxY = ny * bLat;  
700 –
701 –  nLipids = 4 * nx * ny;
702 –  coord* lipidSites = new coord[nLipids];
703 –
704 –  unitVector[0] = 0.0;
705 –  unitVector[1] = 0.0;
706 –
707 –  which = 0;
708 –
709 –  for (i = 0; i < nx; i++) {
710 –    for (j = 0; j < ny; j++ ) {
711 –      for (k = 0; k < 2; k++) {
712 –        
713 –        lipidSites[which].pos[0] = (double)i * aLat;
714 –        lipidSites[which].pos[1] = (double)j * bLat;
715 –        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing / 2.0);
716 –        
717 –        unitVector[2] = 2.0 * (double)k  - 1.0;
718 –        
719 –        getUnitRot( unitVector, lipidSites[which].rot );
720 –
721 –        which++;
722 –
723 –        lipidSites[which].pos[0] = aLat * ((double)i + 0.5);
724 –        lipidSites[which].pos[1] = bLat * ((double)j + 0.5);
725 –        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing / 2.0);
726 –        
727 –        unitVector[2] = 2.0 * (double)k  - 1.0;
728 –        
729 –        getUnitRot( unitVector, lipidSites[which].rot );
730 –        
731 –        which++;
732 –      }
733 –    }
734 –  }
735 –
736 –  targetWaters = targetWaterLipidRatio * nLipids;
570  
738 –  // guess the size of the water box
739 –  
740 –  
741 –
742 –  nCellsX = (int)ceil(boxX / pow(waterVol, ( 1.0 / 3.0 )) );
743 –  nCellsY = (int)ceil(boxY / pow(waterVol, ( 1.0 / 3.0 )) );
744 –
745 –  done = 0;
746 –  nCellsZ = 0;
747 –  while( !done ){
748 –        
749 –    nCellsZ++;
750 –    testTot = 4 * nCellsX * nCellsY * nCellsZ;
751 –        
752 –    if( testTot >= targetWaters ) done = 1;
753 –  }
754 –  
755 –  nWaters = nCellsX * nCellsY * nCellsZ * 4;
756 –  
757 –  coord* waterSites = new coord[nWaters];
758 –
759 –  waterCell[0] = boxX / nCellsX;
760 –  waterCell[1] = boxY / nCellsY;
761 –  waterCell[2] = 4.0 / (waterRho * waterCell[0] * waterCell[1]);
762 –
763 –  Lattice *myORTHO;
764 –  myORTHO = new Lattice( ORTHORHOMBIC_LATTICE_TYPE, waterCell);
765 –  myORTHO->setStartZ( leafSpacing / 2.0 + waterFudge);
766 –
767 –  boxZ = waterCell[2] * nCellsZ;
768 –        
769 –  // create an fcc lattice in the water box.
770 –  
771 –  which = 0;
772 –  for( i=0; i < nCellsX; i++ ){
773 –    for( j=0; j < nCellsY; j++ ){
774 –      for( k=0; k < nCellsZ; k++ ){
775 –
776 –        myORTHO->getLatticePoints(&posX, &posY, &posZ, i, j, k);
777 –        for(l=0; l<4; l++){
778 –          waterSites[which].pos[0] = posX[l];
779 –          waterSites[which].pos[1] = posY[l];
780 –          waterSites[which].pos[2] = posZ[l];
781 –          which++;
782 –        }
783 –      }
784 –    }
785 –  }  
786 –
787 –  // create the real Atom arrays
788 –  
789 –  nAtoms = 0;
790 –  molIndex = 0;
791 –  molStart = new int[nLipids + nWaters];  
792 –  
793 –  for(j=0; j<nLipids; j++){
794 –    molStart[molIndex] = nAtoms;
795 –    molIndex++;
796 –    nAtoms += lipidNatoms;
797 –  }
798 –
799 –  for(j=0; j<nWaters; j++){
800 –    molStart[molIndex] = nAtoms;
801 –    molIndex++;
802 –    nAtoms += waterNatoms;
803 –  }
804 –  
805 –  theConfig = simnfo->getConfiguration();
806 –  theConfig->createArrays( nAtoms );
807 –  simnfo->atoms = new Atom*[nAtoms];
808 –  atoms = simnfo->atoms;
809 –
810 –  // initialize lipid positions
811 –
812 –  molIndex = 0;
813 –  for(i=0; i<nLipids; i++ ){
814 –    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
815 –                           molStart[molIndex], theConfig );
816 –    molIndex++;
817 –  }
818 –
819 –  // initialize the water positions
820 –
821 –  for(i=0; i<nWaters; i++){
822 –    
823 –    getRandomRot( waterSites[i].rot );
824 –    waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
825 –                           molStart[molIndex], theConfig );
826 –    molIndex++;
827 –  }
828 –
829 –  // set up the SimInfo object
830 –  
831 –  double Hmat[3][3];
832 –
833 –  Hmat[0][0] = boxX;
834 –  Hmat[0][1] = 0.0;
835 –  Hmat[0][2] = 0.0;
836 –
837 –  Hmat[1][0] = 0.0;
838 –  Hmat[1][1] = boxY;
839 –  Hmat[1][2] = 0.0;
840 –
841 –  Hmat[2][0] = 0.0;
842 –  Hmat[2][1] = 0.0;
843 –  Hmat[2][2] = boxZ;
844 –  
845 –
846 –  bsInfo.boxX = boxX;
847 –  bsInfo.boxY = boxY;
848 –  bsInfo.boxZ = boxZ;
849 –  
850 –  simnfo->setBoxM( Hmat );
851 –
852 –  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
853 –  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
854 –
855 –  // set up the writer and write out
856 –  
857 –  writer = new DumpWriter( simnfo );
858 –  writer->writeFinal( 0.0 );
859 –
571    return 1;
572   }
573  
863 –
864 – void getRandomRot( double rot[3][3] ){
865 –
866 –  double theta, phi, psi;
867 –  double cosTheta;
868 –
869 –  // select random phi, psi, and cosTheta
870 –
871 –  phi = 2.0 * M_PI * drand48();
872 –  psi = 2.0 * M_PI * drand48();
873 –  cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
874 –
875 –  theta = acos( cosTheta );
876 –
877 –  getEulerRot( theta, phi, psi, rot );
878 – }
879 –
880 –
881 – void getEulerRot( double theta, double phi, double psi, double rot[3][3] ){
882 –
883 –  rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
884 –  rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
885 –  rot[0][2] = sin(theta) * sin(psi);
886 –  
887 –  rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
888 –  rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
889 –  rot[1][2] = sin(theta) * cos(psi);
890 –
891 –  rot[2][0] = sin(phi) * sin(theta);
892 –  rot[2][1] = -cos(phi) * sin(theta);
893 –  rot[2][2] = cos(theta);  
894 – }
895 –
896 –
897 – void getUnitRot( double u[3], double rot[3][3] ){
898 –
899 –  double theta, phi, psi;
900 –
901 –  theta = acos(u[2]);
902 –  phi = atan(u[1] / u[0]);
903 –  psi = 0.0;
904 –
905 –  getEulerRot( theta, phi, psi, rot );
906 – }
907 –                        
908 –
909 –
574   void buildMap( double &x, double &y, double &z,
575                 double boxX, double boxY, double boxZ ){
576    

Diff Legend

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