| 88 | 
  | 
           | 
| 89 | 
  | 
 | 
| 90 | 
  | 
 | 
| 91 | 
< | 
  if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF(); | 
| 92 | 
< | 
  else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF(); | 
| 93 | 
< | 
  else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF(); | 
| 91 | 
> | 
  if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF(); | 
| 92 | 
  | 
  else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF(); | 
| 93 | 
  | 
  else{ | 
| 94 | 
  | 
    sprintf( painCave.errMsg, | 
| 240 | 
  | 
 | 
| 241 | 
  | 
  globalIndex = mpiSim->divideLabor(); | 
| 242 | 
  | 
 | 
| 245 | 
– | 
 | 
| 246 | 
– | 
 | 
| 243 | 
  | 
  // set up the local variables  | 
| 244 | 
  | 
   | 
| 245 | 
  | 
  int localMol, allMol; | 
| 246 | 
  | 
  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | 
| 247 | 
+ | 
 | 
| 248 | 
+ | 
  int* mol2proc = mpiSim->getMolToProcMap(); | 
| 249 | 
+ | 
  int* molCompType = mpiSim->getMolComponentType(); | 
| 250 | 
  | 
   | 
| 251 | 
  | 
  allMol = 0; | 
| 252 | 
  | 
  localMol = 0; | 
| 258 | 
  | 
 | 
| 259 | 
  | 
    for( j=0; j<components_nmol[i]; j++ ){ | 
| 260 | 
  | 
       | 
| 261 | 
< | 
      if( mpiSim->getMyMolStart() <= allMol && | 
| 263 | 
< | 
          allMol <= mpiSim->getMyMolEnd() ){ | 
| 261 | 
> | 
      if( mol2proc[j] == worldRank ){ | 
| 262 | 
  | 
         | 
| 263 | 
  | 
        local_atoms +=    comp_stamps[i]->getNAtoms(); | 
| 264 | 
  | 
        local_bonds +=    comp_stamps[i]->getNBonds(); | 
| 277 | 
  | 
  if( local_atoms != simnfo->n_atoms ){ | 
| 278 | 
  | 
    sprintf( painCave.errMsg, | 
| 279 | 
  | 
             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" | 
| 280 | 
< | 
             " localAtom (%d) are note equal.\n", | 
| 280 | 
> | 
             " localAtom (%d) are not equal.\n", | 
| 281 | 
  | 
             simnfo->n_atoms, | 
| 282 | 
  | 
             local_atoms ); | 
| 283 | 
  | 
    painCave.isFatal = 1; | 
| 302 | 
  | 
  Atom::createArrays(simnfo->n_atoms); | 
| 303 | 
  | 
  the_atoms = new Atom*[simnfo->n_atoms]; | 
| 304 | 
  | 
  the_molecules = new Molecule[simnfo->n_mol]; | 
| 305 | 
+ | 
  int molIndex; | 
| 306 | 
+ | 
 | 
| 307 | 
+ | 
  // initialize the molecule's stampID's | 
| 308 | 
+ | 
 | 
| 309 | 
+ | 
#ifdef IS_MPI | 
| 310 | 
+ | 
   | 
| 311 | 
+ | 
 | 
| 312 | 
+ | 
  molIndex = 0; | 
| 313 | 
+ | 
  for(i=0; i<mpiSim->getTotNmol(); i++){ | 
| 314 | 
+ | 
     | 
| 315 | 
+ | 
    if(mol2proc[i] == worldRank ){ | 
| 316 | 
+ | 
      the_molecules[molIndex].setStampID( molCompType[i] ); | 
| 317 | 
+ | 
      the_molecules[molIndex].setMyIndex( molIndex ); | 
| 318 | 
+ | 
      molIndex++; | 
| 319 | 
+ | 
    } | 
| 320 | 
+ | 
  } | 
| 321 | 
  | 
 | 
| 322 | 
+ | 
#else // is_mpi | 
| 323 | 
+ | 
   | 
| 324 | 
+ | 
  molIndex = 0; | 
| 325 | 
+ | 
  for(i=0; i<n_components; i++){ | 
| 326 | 
+ | 
    for(j=0; j<components_nmol[i]; j++ ){ | 
| 327 | 
+ | 
      the_molecules[molIndex].setStampID( i ); | 
| 328 | 
+ | 
      the_molecules[molIndex].setMyIndex( molIndex ); | 
| 329 | 
+ | 
      molIndex++; | 
| 330 | 
+ | 
    } | 
| 331 | 
+ | 
  } | 
| 332 | 
+ | 
     | 
| 333 | 
  | 
 | 
| 334 | 
+ | 
#endif // is_mpi | 
| 335 | 
+ | 
 | 
| 336 | 
+ | 
 | 
| 337 | 
  | 
  if( simnfo->n_SRI ){ | 
| 338 | 
+ | 
     | 
| 339 | 
  | 
    Exclude::createArray(simnfo->n_SRI); | 
| 340 | 
  | 
    the_excludes = new Exclude*[simnfo->n_SRI]; | 
| 341 | 
+ | 
    for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex); | 
| 342 | 
  | 
    simnfo->globalExcludes = new int; | 
| 343 | 
< | 
    simnfo->n_exclude = tot_SRI; | 
| 343 | 
> | 
    simnfo->n_exclude = simnfo->n_SRI; | 
| 344 | 
  | 
  } | 
| 345 | 
  | 
  else{ | 
| 346 | 
  | 
     | 
| 356 | 
  | 
  // set the arrays into the SimInfo object | 
| 357 | 
  | 
 | 
| 358 | 
  | 
  simnfo->atoms = the_atoms; | 
| 359 | 
< | 
  simnfo->sr_interactions = the_sris; | 
| 359 | 
> | 
  simnfo->molecules = the_molecules; | 
| 360 | 
  | 
  simnfo->nGlobalExcludes = 0; | 
| 361 | 
  | 
  simnfo->excludes = the_excludes; | 
| 362 | 
  | 
 | 
| 413 | 
  | 
 | 
| 414 | 
  | 
  the_ff->setSimInfo( simnfo ); | 
| 415 | 
  | 
 | 
| 416 | 
< | 
  makeAtoms(); | 
| 416 | 
> | 
  makeMolecules(); | 
| 417 | 
  | 
  simnfo->identArray = new int[simnfo->n_atoms]; | 
| 418 | 
  | 
  for(i=0; i<simnfo->n_atoms; i++){ | 
| 419 | 
  | 
    simnfo->identArray[i] = the_atoms[i]->getIdent(); | 
| 420 | 
  | 
  } | 
| 421 | 
  | 
   | 
| 392 | 
– | 
  if( tot_bonds ){ | 
| 393 | 
– | 
    makeBonds(); | 
| 394 | 
– | 
  } | 
| 395 | 
– | 
 | 
| 396 | 
– | 
  if( tot_bends ){ | 
| 397 | 
– | 
    makeBends(); | 
| 398 | 
– | 
  } | 
| 399 | 
– | 
 | 
| 400 | 
– | 
  if( tot_torsions ){ | 
| 401 | 
– | 
    makeTorsions(); | 
| 402 | 
– | 
  } | 
| 403 | 
– | 
 | 
| 404 | 
– | 
 | 
| 422 | 
  | 
  if (the_globals->getUseRF() ) { | 
| 423 | 
  | 
    simnfo->useReactionField = 1; | 
| 424 | 
  | 
   | 
| 655 | 
  | 
 | 
| 656 | 
  | 
//   new AllLong( simnfo ); | 
| 657 | 
  | 
 | 
| 641 | 
– | 
  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff ); | 
| 642 | 
– | 
  if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff ); | 
| 658 | 
  | 
  if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff ); | 
| 659 | 
  | 
  if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff ); | 
| 645 | 
– | 
 | 
| 660 | 
  | 
 | 
| 661 | 
+ | 
#ifdef IS_MPI | 
| 662 | 
+ | 
  mpiSim->mpiRefresh(); | 
| 663 | 
+ | 
#endif | 
| 664 | 
  | 
 | 
| 665 | 
  | 
  // initialize the Fortran | 
| 666 | 
< | 
   | 
| 666 | 
> | 
 | 
| 667 | 
> | 
 | 
| 668 | 
  | 
  simnfo->refreshSim(); | 
| 669 | 
  | 
   | 
| 670 | 
  | 
  if( !strcmp( simnfo->mixingRule, "standard") ){ | 
| 701 | 
  | 
  BondStamp* currentBond; | 
| 702 | 
  | 
  BendStamp* currentBend; | 
| 703 | 
  | 
  TorsionStamp* currentTorsion; | 
| 704 | 
+ | 
 | 
| 705 | 
+ | 
  bond_pair* theBonds; | 
| 706 | 
+ | 
  bend_set* theBends; | 
| 707 | 
+ | 
  torsion_set* theTorsions; | 
| 708 | 
+ | 
 | 
| 709 | 
  | 
   | 
| 710 | 
  | 
  //init the forceField paramters | 
| 711 | 
  | 
 | 
| 712 | 
  | 
  the_ff->readParams(); | 
| 713 | 
  | 
 | 
| 714 | 
  | 
   | 
| 715 | 
< | 
  // init the molecules | 
| 715 | 
> | 
  // init the atoms | 
| 716 | 
  | 
 | 
| 717 | 
+ | 
  double ux, uy, uz, u, uSqr; | 
| 718 | 
+ | 
   | 
| 719 | 
  | 
  atomOffset = 0; | 
| 720 | 
  | 
  excludeOffset = 0; | 
| 721 | 
  | 
  for(i=0; i<simnfo->n_mol; i++){ | 
| 732 | 
  | 
    info.myExcludes = &the_excludes[excludeOffset]; | 
| 733 | 
  | 
    info.myBonds = new Bond*[info.nBonds]; | 
| 734 | 
  | 
    info.myBends = new Bend*[info.nBends]; | 
| 735 | 
< | 
    info.myTorsions = new Torsions*[info.nTorsions]; | 
| 735 | 
> | 
    info.myTorsions = new Torsion*[info.nTorsions]; | 
| 736 | 
  | 
 | 
| 737 | 
  | 
    theBonds = new bond_pair[info.nBonds]; | 
| 738 | 
  | 
    theBends = new bend_set[info.nBends]; | 
| 742 | 
  | 
     | 
| 743 | 
  | 
    for(j=0; j<info.nAtoms; j++){ | 
| 744 | 
  | 
       | 
| 745 | 
< | 
      currentAtom = theComponents[stampID]->getAtom( j ); | 
| 745 | 
> | 
      currentAtom = comp_stamps[stampID]->getAtom( j ); | 
| 746 | 
  | 
      if( currentAtom->haveOrientation() ){ | 
| 747 | 
  | 
         | 
| 748 | 
  | 
        dAtom = new DirectionalAtom(j + atomOffset); | 
| 783 | 
  | 
      theBonds[j].a = currentBond->getA() + atomOffset; | 
| 784 | 
  | 
      theBonds[j].b = currentBond->getB() + atomOffset; | 
| 785 | 
  | 
 | 
| 786 | 
< | 
      exI = theBonds[i].a; | 
| 787 | 
< | 
      exJ = theBonds[i].b; | 
| 786 | 
> | 
      exI = theBonds[j].a; | 
| 787 | 
> | 
      exJ = theBonds[j].b; | 
| 788 | 
  | 
 | 
| 789 | 
  | 
      // exclude_I must always be the smaller of the pair | 
| 790 | 
  | 
      if( exI > exJ ){ | 
| 800 | 
  | 
       | 
| 801 | 
  | 
      the_excludes[j+excludeOffset]->setPair( exI, exJ ); | 
| 802 | 
  | 
#else  // isn't MPI | 
| 803 | 
+ | 
 | 
| 804 | 
  | 
      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); | 
| 805 | 
  | 
#endif  //is_mpi | 
| 806 | 
  | 
    } | 
| 816 | 
  | 
           | 
| 817 | 
  | 
      if( currentBend->haveExtras() ){ | 
| 818 | 
  | 
             | 
| 819 | 
< | 
        extras = current_bend->getExtras(); | 
| 819 | 
> | 
        extras = currentBend->getExtras(); | 
| 820 | 
  | 
        current_extra = extras; | 
| 821 | 
  | 
             | 
| 822 | 
  | 
        while( current_extra != NULL ){ | 
| 838 | 
  | 
               | 
| 839 | 
  | 
            default: | 
| 840 | 
  | 
              sprintf( painCave.errMsg, | 
| 841 | 
< | 
                       "SimSetup Error: ghostVectorSource was neiter a " | 
| 841 | 
> | 
                       "SimSetup Error: ghostVectorSource was neither a " | 
| 842 | 
  | 
                       "double nor an int.\n" | 
| 843 | 
  | 
                       "-->Bend[%d] in %s\n", | 
| 844 | 
  | 
                       j, comp_stamps[stampID]->getID() ); | 
| 932 | 
  | 
 | 
| 933 | 
  | 
 | 
| 934 | 
  | 
    the_molecules[i].initialize( info ); | 
| 935 | 
+ | 
 | 
| 936 | 
+ | 
 | 
| 937 | 
  | 
    atomOffset += info.nAtoms; | 
| 938 | 
+ | 
    delete[] theBonds; | 
| 939 | 
+ | 
    delete[] theBends; | 
| 940 | 
+ | 
    delete[] theTorsions; | 
| 941 | 
  | 
  } | 
| 942 | 
  | 
 | 
| 943 | 
+ | 
#ifdef IS_MPI | 
| 944 | 
+ | 
  sprintf( checkPointMsg, "all molecules initialized succesfully" ); | 
| 945 | 
+ | 
  MPIcheckPoint(); | 
| 946 | 
+ | 
#endif // is_mpi | 
| 947 | 
+ | 
 | 
| 948 | 
  | 
  // clean up the forcefield | 
| 949 | 
  | 
  the_ff->calcRcut(); | 
| 950 | 
  | 
  the_ff->cleanMe(); | 
| 951 | 
+ | 
 | 
| 952 | 
  | 
} | 
| 953 | 
  | 
 | 
| 954 | 
  | 
void SimSetup::initFromBass( void ){ |