| 307 | 
  | 
 | 
| 308 | 
  | 
 | 
| 309 | 
  | 
  if( simnfo->n_SRI ){ | 
| 310 | 
< | 
    the_sris = new SRI*[simnfo->n_SRI]; | 
| 311 | 
< | 
    the_excludes = new int[2 * simnfo->n_SRI]; | 
| 310 | 
> | 
    Exclude::createArray(simnfo->n_SRI); | 
| 311 | 
> | 
    the_excludes = new Exclude*[simnfo->n_SRI]; | 
| 312 | 
  | 
    simnfo->globalExcludes = new int; | 
| 313 | 
  | 
    simnfo->n_exclude = tot_SRI; | 
| 314 | 
  | 
  } | 
| 315 | 
  | 
  else{ | 
| 316 | 
  | 
     | 
| 317 | 
< | 
    the_excludes = new int[2]; | 
| 318 | 
< | 
    the_excludes[0] = 0; | 
| 319 | 
< | 
    the_excludes[1] = 0; | 
| 317 | 
> | 
    Exclude::createArray( 1 ); | 
| 318 | 
> | 
    the_excludes = new Exclude*; | 
| 319 | 
> | 
    the_excludes[0] = new Exclude(0); | 
| 320 | 
> | 
    the_excludes[0]->setPair( 0,0 ); | 
| 321 | 
  | 
    simnfo->globalExcludes = new int; | 
| 322 | 
  | 
    simnfo->globalExcludes[0] = 0; | 
| 323 | 
< | 
 | 
| 323 | 
< | 
    simnfo->n_exclude = 1; | 
| 323 | 
> | 
    simnfo->n_exclude = 0; | 
| 324 | 
  | 
  } | 
| 325 | 
  | 
 | 
| 326 | 
  | 
  // set the arrays into the SimInfo object | 
| 670 | 
  | 
  MPIcheckPoint(); | 
| 671 | 
  | 
#endif // is_mpi | 
| 672 | 
  | 
} | 
| 673 | 
+ | 
 | 
| 674 | 
+ | 
 | 
| 675 | 
+ | 
void SimSetup::makeMolecules( void ){ | 
| 676 | 
+ | 
 | 
| 677 | 
+ | 
  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; | 
| 678 | 
+ | 
  molInit info; | 
| 679 | 
+ | 
  DirectionalAtom* dAtom; | 
| 680 | 
+ | 
  LinkedAssign* extras; | 
| 681 | 
+ | 
  LinkedAssign* current_extra; | 
| 682 | 
+ | 
  AtomStamp* currentAtom; | 
| 683 | 
+ | 
  BondStamp* currentBond; | 
| 684 | 
+ | 
  BendStamp* currentBend; | 
| 685 | 
+ | 
  TorsionStamp* currentTorsion; | 
| 686 | 
+ | 
   | 
| 687 | 
+ | 
  //init the forceField paramters | 
| 688 | 
+ | 
 | 
| 689 | 
+ | 
  the_ff->readParams(); | 
| 690 | 
+ | 
 | 
| 691 | 
+ | 
   | 
| 692 | 
+ | 
  // init the molecules | 
| 693 | 
+ | 
 | 
| 694 | 
+ | 
  atomOffset = 0; | 
| 695 | 
+ | 
  excludeOffset = 0; | 
| 696 | 
+ | 
  for(i=0; i<simnfo->n_mol; i++){ | 
| 697 | 
+ | 
     | 
| 698 | 
+ | 
    stampID = the_molecules[i].getStampID(); | 
| 699 | 
+ | 
 | 
| 700 | 
+ | 
    info.nAtoms    = comp_stamps[stampID]->getNAtoms(); | 
| 701 | 
+ | 
    info.nBonds    = comp_stamps[stampID]->getNBonds(); | 
| 702 | 
+ | 
    info.nBends    = comp_stamps[stampID]->getNBends(); | 
| 703 | 
+ | 
    info.nTorsions = comp_stamps[stampID]->getNTorsions(); | 
| 704 | 
+ | 
    info.nExcludes = info.nBonds + info.nBends + info.nTorsions; | 
| 705 | 
+ | 
 | 
| 706 | 
+ | 
    info.myAtoms = &the_atoms[atomOffset]; | 
| 707 | 
+ | 
    info.myExcludes = &the_excludes[excludeOffset]; | 
| 708 | 
+ | 
    info.myBonds = new Bond*[info.nBonds]; | 
| 709 | 
+ | 
    info.myBends = new Bend*[info.nBends]; | 
| 710 | 
+ | 
    info.myTorsions = new Torsions*[info.nTorsions]; | 
| 711 | 
+ | 
 | 
| 712 | 
+ | 
    theBonds = new bond_pair[info.nBonds]; | 
| 713 | 
+ | 
    theBends = new bend_set[info.nBends]; | 
| 714 | 
+ | 
    theTorsions = new torsion_set[info.nTorsions]; | 
| 715 | 
+ | 
     | 
| 716 | 
+ | 
    // make the Atoms | 
| 717 | 
+ | 
     | 
| 718 | 
+ | 
    for(j=0; j<info.nAtoms; j++){ | 
| 719 | 
+ | 
       | 
| 720 | 
+ | 
      currentAtom = theComponents[stampID]->getAtom( j ); | 
| 721 | 
+ | 
      if( currentAtom->haveOrientation() ){ | 
| 722 | 
+ | 
         | 
| 723 | 
+ | 
        dAtom = new DirectionalAtom(j + atomOffset); | 
| 724 | 
+ | 
        simnfo->n_oriented++; | 
| 725 | 
+ | 
        info.myAtoms[j] = dAtom; | 
| 726 | 
+ | 
         | 
| 727 | 
+ | 
        ux = currentAtom->getOrntX(); | 
| 728 | 
+ | 
        uy = currentAtom->getOrntY(); | 
| 729 | 
+ | 
        uz = currentAtom->getOrntZ(); | 
| 730 | 
+ | 
         | 
| 731 | 
+ | 
        uSqr = (ux * ux) + (uy * uy) + (uz * uz); | 
| 732 | 
+ | 
         | 
| 733 | 
+ | 
        u = sqrt( uSqr ); | 
| 734 | 
+ | 
        ux = ux / u; | 
| 735 | 
+ | 
        uy = uy / u; | 
| 736 | 
+ | 
        uz = uz / u; | 
| 737 | 
+ | 
         | 
| 738 | 
+ | 
        dAtom->setSUx( ux ); | 
| 739 | 
+ | 
        dAtom->setSUy( uy ); | 
| 740 | 
+ | 
        dAtom->setSUz( uz ); | 
| 741 | 
+ | 
      } | 
| 742 | 
+ | 
      else{ | 
| 743 | 
+ | 
        info.myAtoms[j] = new GeneralAtom(j + atomOffset); | 
| 744 | 
+ | 
      } | 
| 745 | 
+ | 
      info.myAtoms[j]->setType( currentAtom->getType() ); | 
| 746 | 
+ | 
     | 
| 747 | 
+ | 
#ifdef IS_MPI | 
| 748 | 
+ | 
       | 
| 749 | 
+ | 
      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] ); | 
| 750 | 
+ | 
       | 
| 751 | 
+ | 
#endif // is_mpi | 
| 752 | 
+ | 
    }  | 
| 753 | 
+ | 
     | 
| 754 | 
+ | 
    // make the bonds | 
| 755 | 
+ | 
    for(j=0; j<info.nBonds; j++){ | 
| 756 | 
+ | 
       | 
| 757 | 
+ | 
      currentBond = comp_stamps[stampID]->getBond( j ); | 
| 758 | 
+ | 
      theBonds[j].a = currentBond->getA() + atomOffset; | 
| 759 | 
+ | 
      theBonds[j].b = currentBond->getB() + atomOffset; | 
| 760 | 
+ | 
 | 
| 761 | 
+ | 
      exI = theBonds[i].a; | 
| 762 | 
+ | 
      exJ = theBonds[i].b; | 
| 763 | 
+ | 
 | 
| 764 | 
+ | 
      // exclude_I must always be the smaller of the pair | 
| 765 | 
+ | 
      if( exI > exJ ){ | 
| 766 | 
+ | 
        tempEx = exI; | 
| 767 | 
+ | 
        exI = exJ; | 
| 768 | 
+ | 
        exJ = tempEx; | 
| 769 | 
+ | 
      } | 
| 770 | 
+ | 
#ifdef IS_MPI | 
| 771 | 
+ | 
      tempEx = exI; | 
| 772 | 
+ | 
      exI = the_atoms[tempEx]->getGlobalIndex() + 1; | 
| 773 | 
+ | 
      tempEx = exJ; | 
| 774 | 
+ | 
      exJ = the_atoms[tempEx]->getGlobalIndex() + 1; | 
| 775 | 
+ | 
       | 
| 776 | 
+ | 
      the_excludes[j+excludeOffset]->setPair( exI, exJ ); | 
| 777 | 
+ | 
#else  // isn't MPI | 
| 778 | 
+ | 
      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); | 
| 779 | 
+ | 
#endif  //is_mpi | 
| 780 | 
+ | 
    } | 
| 781 | 
+ | 
    excludeOffset += info.nBonds; | 
| 782 | 
+ | 
 | 
| 783 | 
+ | 
    //make the bends | 
| 784 | 
+ | 
    for(j=0; j<info.nBends; j++){ | 
| 785 | 
+ | 
       | 
| 786 | 
+ | 
      currentBend = comp_stamps[stampID]->getBend( j ); | 
| 787 | 
+ | 
      theBends[j].a = currentBend->getA() + atomOffset; | 
| 788 | 
+ | 
      theBends[j].b = currentBend->getB() + atomOffset; | 
| 789 | 
+ | 
      theBends[j].c = currentBend->getC() + atomOffset; | 
| 790 | 
+ | 
           | 
| 791 | 
+ | 
      if( currentBend->haveExtras() ){ | 
| 792 | 
+ | 
             | 
| 793 | 
+ | 
        extras = current_bend->getExtras(); | 
| 794 | 
+ | 
        current_extra = extras; | 
| 795 | 
+ | 
             | 
| 796 | 
+ | 
        while( current_extra != NULL ){ | 
| 797 | 
+ | 
          if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){ | 
| 798 | 
+ | 
                 | 
| 799 | 
+ | 
            switch( current_extra->getType() ){ | 
| 800 | 
+ | 
               | 
| 801 | 
+ | 
            case 0: | 
| 802 | 
+ | 
              theBends[j].ghost = | 
| 803 | 
+ | 
                current_extra->getInt() + atomOffset; | 
| 804 | 
+ | 
              theBends[j].isGhost = 1; | 
| 805 | 
+ | 
              break; | 
| 806 | 
+ | 
                   | 
| 807 | 
+ | 
            case 1: | 
| 808 | 
+ | 
              theBends[j].ghost =  | 
| 809 | 
+ | 
                (int)current_extra->getDouble() + atomOffset; | 
| 810 | 
+ | 
              theBends[j].isGhost = 1; | 
| 811 | 
+ | 
              break; | 
| 812 | 
+ | 
               | 
| 813 | 
+ | 
            default: | 
| 814 | 
+ | 
              sprintf( painCave.errMsg, | 
| 815 | 
+ | 
                       "SimSetup Error: ghostVectorSource was neiter a " | 
| 816 | 
+ | 
                       "double nor an int.\n" | 
| 817 | 
+ | 
                       "-->Bend[%d] in %s\n", | 
| 818 | 
+ | 
                       j, comp_stamps[stampID]->getID() ); | 
| 819 | 
+ | 
              painCave.isFatal = 1; | 
| 820 | 
+ | 
              simError(); | 
| 821 | 
+ | 
            } | 
| 822 | 
+ | 
          } | 
| 823 | 
+ | 
           | 
| 824 | 
+ | 
          else{ | 
| 825 | 
+ | 
             | 
| 826 | 
+ | 
            sprintf( painCave.errMsg, | 
| 827 | 
+ | 
                     "SimSetup Error: unhandled bend assignment:\n" | 
| 828 | 
+ | 
                     "    -->%s in Bend[%d] in %s\n", | 
| 829 | 
+ | 
                     current_extra->getlhs(), | 
| 830 | 
+ | 
                     j, comp_stamps[stampID]->getID() ); | 
| 831 | 
+ | 
            painCave.isFatal = 1; | 
| 832 | 
+ | 
            simError(); | 
| 833 | 
+ | 
          } | 
| 834 | 
+ | 
           | 
| 835 | 
+ | 
          current_extra = current_extra->getNext(); | 
| 836 | 
+ | 
        } | 
| 837 | 
+ | 
      } | 
| 838 | 
+ | 
           | 
| 839 | 
+ | 
      if( !theBends[j].isGhost ){ | 
| 840 | 
+ | 
             | 
| 841 | 
+ | 
        exI = theBends[j].a; | 
| 842 | 
+ | 
        exJ = theBends[j].c; | 
| 843 | 
+ | 
      } | 
| 844 | 
+ | 
      else{ | 
| 845 | 
+ | 
         | 
| 846 | 
+ | 
        exI = theBends[j].a; | 
| 847 | 
+ | 
        exJ = theBends[j].b; | 
| 848 | 
+ | 
      } | 
| 849 | 
+ | 
       | 
| 850 | 
+ | 
      // exclude_I must always be the smaller of the pair | 
| 851 | 
+ | 
      if( exI > exJ ){ | 
| 852 | 
+ | 
        tempEx = exI; | 
| 853 | 
+ | 
        exI = exJ; | 
| 854 | 
+ | 
        exJ = tempEx; | 
| 855 | 
+ | 
      } | 
| 856 | 
+ | 
#ifdef IS_MPI | 
| 857 | 
+ | 
      tempEx = exI; | 
| 858 | 
+ | 
      exI = the_atoms[tempEx]->getGlobalIndex() + 1; | 
| 859 | 
+ | 
      tempEx = exJ; | 
| 860 | 
+ | 
      exJ = the_atoms[tempEx]->getGlobalIndex() + 1; | 
| 861 | 
+ | 
       | 
| 862 | 
+ | 
      the_excludes[j+excludeOffset]->setPair( exI, exJ ); | 
| 863 | 
+ | 
#else  // isn't MPI | 
| 864 | 
+ | 
      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); | 
| 865 | 
+ | 
#endif  //is_mpi | 
| 866 | 
+ | 
    } | 
| 867 | 
+ | 
    excludeOffset += info.nBends; | 
| 868 | 
+ | 
 | 
| 869 | 
+ | 
    for(j=0; j<info.nTorsions; j++){ | 
| 870 | 
+ | 
       | 
| 871 | 
+ | 
      currentTorsion = comp_stamps[stampID]->getTorsion( j ); | 
| 872 | 
+ | 
      theTorsions[j].a = currentTorsion->getA() + atomOffset; | 
| 873 | 
+ | 
      theTorsions[j].b = currentTorsion->getB() + atomOffset; | 
| 874 | 
+ | 
      theTorsions[j].c = currentTorsion->getC() + atomOffset; | 
| 875 | 
+ | 
      theTorsions[j].d = currentTorsion->getD() + atomOffset; | 
| 876 | 
+ | 
       | 
| 877 | 
+ | 
      exI = theTorsions[j].a; | 
| 878 | 
+ | 
      exJ = theTorsions[j].d; | 
| 879 | 
  | 
 | 
| 880 | 
+ | 
      // exclude_I must always be the smaller of the pair | 
| 881 | 
+ | 
      if( exI > exJ ){ | 
| 882 | 
+ | 
        tempEx = exI; | 
| 883 | 
+ | 
        exI = exJ; | 
| 884 | 
+ | 
        exJ = tempEx; | 
| 885 | 
+ | 
      } | 
| 886 | 
+ | 
#ifdef IS_MPI | 
| 887 | 
+ | 
      tempEx = exI; | 
| 888 | 
+ | 
      exI = the_atoms[tempEx]->getGlobalIndex() + 1; | 
| 889 | 
+ | 
      tempEx = exJ; | 
| 890 | 
+ | 
      exJ = the_atoms[tempEx]->getGlobalIndex() + 1; | 
| 891 | 
+ | 
       | 
| 892 | 
+ | 
      the_excludes[j+excludeOffset]->setPair( exI, exJ ); | 
| 893 | 
+ | 
#else  // isn't MPI | 
| 894 | 
+ | 
      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); | 
| 895 | 
+ | 
#endif  //is_mpi | 
| 896 | 
+ | 
    } | 
| 897 | 
+ | 
    excludeOffset += info.nTorsions; | 
| 898 | 
+ | 
 | 
| 899 | 
+ | 
     | 
| 900 | 
+ | 
 | 
| 901 | 
+ | 
 | 
| 902 | 
+ | 
 | 
| 903 | 
+ | 
 | 
| 904 | 
+ | 
 | 
| 905 | 
+ | 
 | 
| 906 | 
  | 
void SimSetup::makeAtoms( void ){ | 
| 907 | 
  | 
 | 
| 908 | 
  | 
  int i, j, k, index; |