ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-1.0/libmdtools/SimSetup.cpp
(Generate patch)

Comparing trunk/OOPSE-1.0/libmdtools/SimSetup.cpp (file contents):
Revision 1418 by gezelter, Tue Jul 27 16:13:29 2004 UTC vs.
Revision 1419 by gezelter, Tue Jul 27 18:14:16 2004 UTC

# Line 56 | Line 56 | bool isDivisible(double dividend, double divisor){
56    else
57      return false;  
58   }
59 +
60 + string getPrefix(const string& str ){
61 +  string prefix;
62 +  string suffix;
63 +  int pos;
64 +
65 +  pos = str.rfind(".");
66 +
67 +  if (pos >= 0) {
68 +     prefix = str.substr(0, pos);
69 +     suffix = str.substr(pos, str.size());
70 +
71 +     // leave .bass there in case we've reverted to old habits
72 +     if (LowerCase(suffix) == ".md" || LowerCase(suffix) == ".bass")
73 +       return prefix;
74 +     else
75 +       return str;
76 +    
77 +  } else
78 +    return str;
79 + };
80  
81 +
82   SimSetup::SimSetup(){
83    
84    initSuspend = false;
# Line 645 | Line 667 | void SimSetup::makeMolecules(void){
667   #endif // is_mpi
668  
669   }
648 –
649 – void SimSetup::initFromMetaDataFile(void){
650 –  int i, j, k;
651 –  int n_cells;
652 –  double cellx, celly, cellz;
653 –  double temp1, temp2, temp3;
654 –  int n_per_extra;
655 –  int n_extra;
656 –  int have_extra, done;
670  
658 –  double vel[3];
659 –  vel[0] = 0.0;
660 –  vel[1] = 0.0;
661 –  vel[2] = 0.0;
662 –
663 –  temp1 = (double) tot_nmol / 4.0;
664 –  temp2 = pow(temp1, (1.0 / 3.0));
665 –  temp3 = ceil(temp2);
666 –
667 –  have_extra = 0;
668 –  if (temp2 < temp3){
669 –    // we have a non-complete lattice
670 –    have_extra = 1;
671 –
672 –    n_cells = (int) temp3 - 1;
673 –    cellx = info[0].boxL[0] / temp3;
674 –    celly = info[0].boxL[1] / temp3;
675 –    cellz = info[0].boxL[2] / temp3;
676 –    n_extra = tot_nmol - (4 * n_cells * n_cells * n_cells);
677 –    temp1 = ((double) n_extra) / (pow(temp3, 3.0) - pow(n_cells, 3.0));
678 –    n_per_extra = (int) ceil(temp1);
679 –
680 –    if (n_per_extra > 4){
681 –      sprintf(painCave.errMsg,
682 –              "SimSetup error. There has been an error in constructing"
683 –              " the non-complete lattice.\n");
684 –      painCave.isFatal = 1;
685 –      simError();
686 –    }
687 –  }
688 –  else{
689 –    n_cells = (int) temp3;
690 –    cellx = info[0].boxL[0] / temp3;
691 –    celly = info[0].boxL[1] / temp3;
692 –    cellz = info[0].boxL[2] / temp3;
693 –  }
694 –
695 –  current_mol = 0;
696 –  current_comp_mol = 0;
697 –  current_comp = 0;
698 –  current_atom_ndx = 0;
699 –
700 –  for (i = 0; i < n_cells ; i++){
701 –    for (j = 0; j < n_cells; j++){
702 –      for (k = 0; k < n_cells; k++){
703 –        makeElement(i * cellx, j * celly, k * cellz);
704 –
705 –        makeElement(i * cellx + 0.5 * cellx, j * celly + 0.5 * celly, k * cellz);
706 –
707 –        makeElement(i * cellx, j * celly + 0.5 * celly, k * cellz + 0.5 * cellz);
708 –
709 –        makeElement(i * cellx + 0.5 * cellx, j * celly, k * cellz + 0.5 * cellz);
710 –      }
711 –    }
712 –  }
713 –
714 –  if (have_extra){
715 –    done = 0;
716 –
717 –    int start_ndx;
718 –    for (i = 0; i < (n_cells + 1) && !done; i++){
719 –      for (j = 0; j < (n_cells + 1) && !done; j++){
720 –        if (i < n_cells){
721 –          if (j < n_cells){
722 –            start_ndx = n_cells;
723 –          }
724 –          else
725 –            start_ndx = 0;
726 –        }
727 –        else
728 –          start_ndx = 0;
729 –
730 –        for (k = start_ndx; k < (n_cells + 1) && !done; k++){
731 –          makeElement(i * cellx, j * celly, k * cellz);
732 –          done = (current_mol >= tot_nmol);
733 –
734 –          if (!done && n_per_extra > 1){
735 –            makeElement(i * cellx + 0.5 * cellx, j * celly + 0.5 * celly,
736 –                        k * cellz);
737 –            done = (current_mol >= tot_nmol);
738 –          }
739 –
740 –          if (!done && n_per_extra > 2){
741 –            makeElement(i * cellx, j * celly + 0.5 * celly,
742 –                        k * cellz + 0.5 * cellz);
743 –            done = (current_mol >= tot_nmol);
744 –          }
745 –
746 –          if (!done && n_per_extra > 3){
747 –            makeElement(i * cellx + 0.5 * cellx, j * celly,
748 –                        k * cellz + 0.5 * cellz);
749 –            done = (current_mol >= tot_nmol);
750 –          }
751 –        }
752 –      }
753 –    }
754 –  }
755 –
756 –  for (i = 0; i < info[0].n_atoms; i++){
757 –    info[0].atoms[i]->setVel(vel);
758 –  }
759 – }
760 –
761 – void SimSetup::makeElement(double x, double y, double z){
762 –  int k;
763 –  AtomStamp* current_atom;
764 –  DirectionalAtom* dAtom;
765 –  double rotMat[3][3];
766 –  double pos[3];
767 –
768 –  for (k = 0; k < comp_stamps[current_comp]->getNAtoms(); k++){
769 –    current_atom = comp_stamps[current_comp]->getAtom(k);
770 –    if (!current_atom->havePosition()){
771 –      sprintf(painCave.errMsg,
772 –              "SimSetup:initFromMetaData error.\n"
773 –              "\tComponent %s, atom %s does not have a position specified.\n"
774 –              "\tThe initialization routine is unable to give a start"
775 –              " position.\n",
776 –              comp_stamps[current_comp]->getID(), current_atom->getType());
777 –      painCave.isFatal = 1;
778 –      simError();
779 –    }
780 –
781 –    pos[0] = x + current_atom->getPosX();
782 –    pos[1] = y + current_atom->getPosY();
783 –    pos[2] = z + current_atom->getPosZ();
784 –
785 –    info[0].atoms[current_atom_ndx]->setPos(pos);
786 –
787 –    if (info[0].atoms[current_atom_ndx]->isDirectional()){
788 –      dAtom = (DirectionalAtom *) info[0].atoms[current_atom_ndx];
789 –
790 –      rotMat[0][0] = 1.0;
791 –      rotMat[0][1] = 0.0;
792 –      rotMat[0][2] = 0.0;
793 –
794 –      rotMat[1][0] = 0.0;
795 –      rotMat[1][1] = 1.0;
796 –      rotMat[1][2] = 0.0;
797 –
798 –      rotMat[2][0] = 0.0;
799 –      rotMat[2][1] = 0.0;
800 –      rotMat[2][2] = 1.0;
801 –
802 –      dAtom->setA(rotMat);
803 –    }
804 –
805 –    current_atom_ndx++;
806 –  }
807 –
808 –  current_mol++;
809 –  current_comp_mol++;
810 –
811 –  if (current_comp_mol >= components_nmol[current_comp]){
812 –    current_comp_mol = 0;
813 –    current_comp++;
814 –  }
815 – }
816 –
817 –
671   void SimSetup::gatherInfo(void){
672    int i;
673  
# Line 1330 | Line 1183 | void SimSetup::makeOutNames(void){
1183      if (worldRank == 0){
1184   #endif // is_mpi
1185        
1186 <      if(globals->haveFinalConfig())
1186 >      if(globals->haveFinalConfig())
1187          prefix = getPrefix(globals->getFinalConfig());  
1188        else
1189 <        prefix = getPrefix(info[k].finalName);
1189 >        prefix = getPrefix(inFileName);
1190  
1191        info[k].finalName = prefix + ".eor";      
1192        info[k].sampleName = prefix + ".dump";

Diff Legend

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