src/brains/mpiSimulation.cpp

#ifdef IS_MPI
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>

#include "brains/mpiSimulation.hpp"
#include "utils/simError.h"
#include "UseTheForce/DarkSide/simParallel_interface.h"
#include "math/randomSPRNG.hpp"

mpiSimulation* mpiSim;

mpiSimulation::mpiSimulation(SimInfo* the_entryPlug)
{
  entryPlug = the_entryPlug;
  parallelData = new mpiSimData;
  
  MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors) );
  parallelData->myNode = worldRank;

  MolToProcMap = new int[entryPlug->n_mol];
  MolComponentType = new int[entryPlug->n_mol];
  AtomToProcMap = new int[entryPlug->n_atoms];
  GroupToProcMap = new int[entryPlug->ngroup];

  mpiSim = this;
}


mpiSimulation::~mpiSimulation(){
  
  delete[] MolToProcMap;
  delete[] MolComponentType;
  delete[] AtomToProcMap;
  delete[] GroupToProcMap;

  delete parallelData;
  // perhaps we should let fortran know the party is over.
  
}

void mpiSimulation::divideLabor( ){

  int nComponents;
  MoleculeStamp** compStamps;
  randomSPRNG *myRandom;
  int* componentsNmol;
  int* AtomsPerProc;
  int* GroupsPerProc;

  double numerator;
  double denominator;
  double precast;
  double x, y, a;
  int old_atoms, add_atoms, new_atoms;
  int old_groups, add_groups, new_groups;

  int nTarget;
  int molIndex, atomIndex, groupIndex;
  int done;
  int i, j, loops, which_proc;
  int nmol_global, nmol_local;
  int ngroups_global, ngroups_local;
  int natoms_global, natoms_local;
  int ncutoff_groups, nAtomsInGroups;
  int local_index;
  int baseSeed = entryPlug->getSeed();
  CutoffGroupStamp* cg;

  nComponents = entryPlug->nComponents;
  compStamps = entryPlug->compStamps;
  componentsNmol = entryPlug->componentsNmol;
  AtomsPerProc = new int[parallelData->nProcessors];
  GroupsPerProc = new int[parallelData->nProcessors];
  
  parallelData->nAtomsGlobal = entryPlug->n_atoms;
  parallelData->nBondsGlobal = entryPlug->n_bonds;
  parallelData->nBendsGlobal = entryPlug->n_bends;
  parallelData->nTorsionsGlobal = entryPlug->n_torsions;
  parallelData->nSRIGlobal = entryPlug->n_SRI;
  parallelData->nGroupsGlobal = entryPlug->ngroup;
  parallelData->nMolGlobal = entryPlug->n_mol;

  if (parallelData->nProcessors > parallelData->nMolGlobal) {
    sprintf( painCave.errMsg,
             "nProcessors (%d) > nMol (%d)\n"
             "\tThe number of processors is larger than\n"
             "\tthe number of molecules.  This will not result in a \n"
             "\tusable division of atoms for force decomposition.\n"
             "\tEither try a smaller number of processors, or run the\n"
             "\tsingle-processor version of OOPSE.\n",
             parallelData->nProcessors, parallelData->nMolGlobal );
    painCave.isFatal = 1;
    simError();
  }

  myRandom = new randomSPRNG( baseSeed );  


  a = 3.0 * (double)parallelData->nMolGlobal / (double)parallelData->nAtomsGlobal;

  // Initialize things that we'll send out later:
  for (i = 0; i < parallelData->nProcessors; i++ ) {
    AtomsPerProc[i] = 0;
    GroupsPerProc[i] = 0;
  }
  for (i = 0; i < parallelData->nMolGlobal; i++ ) {
    // default to an error condition:
    MolToProcMap[i] = -1;
    MolComponentType[i] = -1;
  }
  for (i = 0; i < parallelData->nAtomsGlobal; i++ ) {
    // default to an error condition:
    AtomToProcMap[i] = -1;
  }
  for (i = 0; i < parallelData->nGroupsGlobal; i++ ) {
    // default to an error condition:
    GroupToProcMap[i] = -1;
  }
    
  if (parallelData->myNode == 0) {
    numerator = (double) entryPlug->n_atoms;
    denominator = (double) parallelData->nProcessors;
    precast = numerator / denominator;
    nTarget = (int)( precast + 0.5 );

    // Build the array of molecule component types first
    molIndex = 0;
    for (i=0; i < nComponents; i++) {
      for (j=0; j < componentsNmol[i]; j++) {        
        MolComponentType[molIndex] = i;
        molIndex++;
      }
    }

    atomIndex = 0;
    groupIndex = 0;

    for (i = 0; i < molIndex; i++ ) {

      done = 0;
      loops = 0;

      while( !done ){
        loops++;
        
        // Pick a processor at random

        which_proc = (int) (myRandom->getRandom() * parallelData->nProcessors);

        // How many atoms does this processor have?
        
        old_atoms = AtomsPerProc[which_proc];
        add_atoms = compStamps[MolComponentType[i]]->getNAtoms();
        new_atoms = old_atoms + add_atoms;

        old_groups = GroupsPerProc[which_proc];
        ncutoff_groups = compStamps[MolComponentType[i]]->getNCutoffGroups(); 
        nAtomsInGroups = 0;
        for (j = 0; j < ncutoff_groups; j++) {
          cg = compStamps[MolComponentType[i]]->getCutoffGroup(j);
          nAtomsInGroups += cg->getNMembers();
        }
        add_groups = add_atoms - nAtomsInGroups + ncutoff_groups;        
        new_groups = old_groups + add_groups;

        // If we've been through this loop too many times, we need
        // to just give up and assign the molecule to this processor
        // and be done with it. 
        
        if (loops > 100) {          
          sprintf( painCave.errMsg,
                   "I've tried 100 times to assign molecule %d to a "
                   " processor, but can't find a good spot.\n"  
                   "I'm assigning it at random to processor %d.\n",
                   i, which_proc);
          painCave.isFatal = 0;
          simError();
          
          MolToProcMap[i] = which_proc;
          AtomsPerProc[which_proc] += add_atoms;
          for (j = 0 ; j < add_atoms; j++ ) {
            AtomToProcMap[atomIndex] = which_proc;
            atomIndex++;
          }
          GroupsPerProc[which_proc] += add_groups;
          for (j=0; j < add_groups; j++) {
            GroupToProcMap[groupIndex] = which_proc;
            groupIndex++;
          }
          done = 1;
          continue;
        }
    
        // If we can add this molecule to this processor without sending
        // it above nTarget, then go ahead and do it:
    
        if (new_atoms <= nTarget) {
          MolToProcMap[i] = which_proc;
          AtomsPerProc[which_proc] += add_atoms;
          for (j = 0 ; j < add_atoms; j++ ) {
            AtomToProcMap[atomIndex] = which_proc;
            atomIndex++;
          }
          GroupsPerProc[which_proc] += add_groups;
          for (j=0; j < add_groups; j++) {
            GroupToProcMap[groupIndex] = which_proc;
            groupIndex++;
          }
          done = 1;
          continue;
        }


        // The only situation left is when new_atoms > nTarget.  We
        // want to accept this with some probability that dies off the
        // farther we are from nTarget

        // roughly:  x = new_atoms - nTarget
        //           Pacc(x) = exp(- a * x)
        // where a = penalty / (average atoms per molecule)

        x = (double) (new_atoms - nTarget);
        y = myRandom->getRandom();
      
        if (y < exp(- a * x)) {
          MolToProcMap[i] = which_proc;
          AtomsPerProc[which_proc] += add_atoms;
          for (j = 0 ; j < add_atoms; j++ ) {
            AtomToProcMap[atomIndex] = which_proc;
            atomIndex++;
           }
          GroupsPerProc[which_proc] += add_groups;
          for (j=0; j < add_groups; j++) {
            GroupToProcMap[groupIndex] = which_proc;
            groupIndex++;
          }
          done = 1;
          continue;
        } else {
          continue;
        }       
        
      }
    }


    // Spray out this nonsense to all other processors:

    //std::cerr << "node 0 mol2proc = \n";
    //for (i = 0; i < parallelData->nMolGlobal; i++) 
    //  std::cerr << i << "\t" << MolToProcMap[i] << "\n";

    MPI_Bcast(MolToProcMap, parallelData->nMolGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(AtomToProcMap, parallelData->nAtomsGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(GroupToProcMap, parallelData->nGroupsGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(MolComponentType, parallelData->nMolGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(AtomsPerProc, parallelData->nProcessors,
              MPI_INT, 0, MPI_COMM_WORLD);    

    MPI_Bcast(GroupsPerProc, parallelData->nProcessors,
              MPI_INT, 0, MPI_COMM_WORLD);    
  } else {

    // Listen to your marching orders from processor 0:
    
    MPI_Bcast(MolToProcMap, parallelData->nMolGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);
    
    MPI_Bcast(AtomToProcMap, parallelData->nAtomsGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(GroupToProcMap, parallelData->nGroupsGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(MolComponentType, parallelData->nMolGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);
    
    MPI_Bcast(AtomsPerProc, parallelData->nProcessors,
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(GroupsPerProc, parallelData->nProcessors,
              MPI_INT, 0, MPI_COMM_WORLD);


  }

  // Let's all check for sanity:

  nmol_local = 0;
  for (i = 0 ; i < parallelData->nMolGlobal; i++ ) {
    if (MolToProcMap[i] == parallelData->myNode) {
      nmol_local++;
    }
  }

  natoms_local = 0;
  for (i = 0; i < parallelData->nAtomsGlobal; i++) {
    if (AtomToProcMap[i] == parallelData->myNode) {
      natoms_local++;      
    }
  }

  ngroups_local = 0;
  for (i = 0; i < parallelData->nGroupsGlobal; i++) {
    if (GroupToProcMap[i] == parallelData->myNode) {
      ngroups_local++;      
    }
  }

  MPI_Allreduce(&nmol_local,&nmol_global,1,MPI_INT,MPI_SUM, 
                MPI_COMM_WORLD);

  MPI_Allreduce(&natoms_local,&natoms_global,1,MPI_INT,
                MPI_SUM, MPI_COMM_WORLD);

  MPI_Allreduce(&ngroups_local,&ngroups_global,1,MPI_INT,
                MPI_SUM, MPI_COMM_WORLD);
  
  if( nmol_global != entryPlug->n_mol ){
    sprintf( painCave.errMsg,
             "The sum of all nmol_local, %d, did not equal the "
             "total number of molecules, %d.\n",
             nmol_global, entryPlug->n_mol );
    painCave.isFatal = 1;
    simError();
  }
  
  if( natoms_global != entryPlug->n_atoms ){
    sprintf( painCave.errMsg,
             "The sum of all natoms_local, %d, did not equal the "
             "total number of atoms, %d.\n",
             natoms_global, entryPlug->n_atoms );
    painCave.isFatal = 1;
    simError();
  }

  if( ngroups_global != entryPlug->ngroup ){
    sprintf( painCave.errMsg,
             "The sum of all ngroups_local, %d, did not equal the "
             "total number of cutoffGroups, %d.\n",
             ngroups_global, entryPlug->ngroup );
    painCave.isFatal = 1;
    simError();
  }

  sprintf( checkPointMsg,
           "Successfully divided the molecules among the processors.\n" );
  MPIcheckPoint();

  parallelData->nMolLocal = nmol_local;
  parallelData->nAtomsLocal = natoms_local;
  parallelData->nGroupsLocal = ngroups_local;

  globalAtomIndex.resize(parallelData->nAtomsLocal);
  globalToLocalAtom.resize(parallelData->nAtomsGlobal);
  local_index = 0;
  for (i = 0; i < parallelData->nAtomsGlobal; i++) {
    if (AtomToProcMap[i] == parallelData->myNode) {
      globalAtomIndex[local_index] = i;

      globalToLocalAtom[i] = local_index;
      local_index++;
      
    }
    else
       globalToLocalAtom[i] = -1;
  }

  globalGroupIndex.resize(parallelData->nGroupsLocal);
  globalToLocalGroup.resize(parallelData->nGroupsGlobal);
  local_index = 0;
  for (i = 0; i < parallelData->nGroupsGlobal; i++) {
    if (GroupToProcMap[i] == parallelData->myNode) {
      globalGroupIndex[local_index] = i;

      globalToLocalGroup[i] = local_index;
      local_index++;
      
    }
    else
       globalToLocalGroup[i] = -1;
  }

  globalMolIndex.resize(parallelData->nMolLocal);
  globalToLocalMol.resize(parallelData->nMolGlobal);  
  local_index = 0;
  for (i = 0; i < parallelData->nMolGlobal; i++) {
    if (MolToProcMap[i] == parallelData->myNode) {
      globalMolIndex[local_index] = i;
      globalToLocalMol[i] = local_index;
      local_index++;
    }
    else
      globalToLocalMol[i] = -1;
  }
  
}


void mpiSimulation::mpiRefresh( void ){

  int isError, i;
  int *localToGlobalAtomIndex = new int[parallelData->nAtomsLocal];
  int *localToGlobalGroupIndex = new int[parallelData->nGroupsLocal];

  // Fortran indexing needs to be increased by 1 in order to get the 2
  // languages to not barf

  for(i = 0; i < parallelData->nAtomsLocal; i++) 
    localToGlobalAtomIndex[i] = globalAtomIndex[i] + 1;

  for(i = 0; i < parallelData->nGroupsLocal; i++) 
    localToGlobalGroupIndex[i] = globalGroupIndex[i] + 1;
  
  isError = 0;

  setFsimParallel( parallelData, 
                   &(parallelData->nAtomsLocal), localToGlobalAtomIndex, 
                   &(parallelData->nGroupsLocal),  localToGlobalGroupIndex, 
                   &isError );

  if( isError ){

    sprintf( painCave.errMsg,
             "mpiRefresh errror: fortran didn't like something we gave it.\n" );
    painCave.isFatal = 1;
    simError();
  }

  delete[] localToGlobalGroupIndex;
  delete[] localToGlobalAtomIndex;


  sprintf( checkPointMsg,
           " mpiRefresh successful.\n" );
  MPIcheckPoint();
}
 

#endif // is_mpi
Revision:	1619
Committed:	Wed Oct 20 21:16:01 2004 UTC (19 years, 8 months ago) by chuckv
File size:	13042 byte(s)
Log Message:	Fortran/C++ interface de-obfuscation project appears to be complete! Woo hoo!
#	User	Rev	Content
1	gezelter	1490	#ifdef IS_MPI
2			#include <iostream>
3			#include <stdlib.h>
4			#include <string.h>
5			#include <math.h>
6			#include <mpi.h>
7
8	tim	1492	#include "brains/mpiSimulation.hpp"
9			#include "utils/simError.h"
10	chuckv	1619	#include "UseTheForce/DarkSide/simParallel_interface.h"
11	tim	1492	#include "math/randomSPRNG.hpp"
12	gezelter	1490
13			mpiSimulation* mpiSim;
14
15			mpiSimulation::mpiSimulation(SimInfo* the_entryPlug)
16			{
17			entryPlug = the_entryPlug;
18			parallelData = new mpiSimData;
19
20			MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors) );
21			parallelData->myNode = worldRank;
22
23			MolToProcMap = new int[entryPlug->n_mol];
24			MolComponentType = new int[entryPlug->n_mol];
25			AtomToProcMap = new int[entryPlug->n_atoms];
26			GroupToProcMap = new int[entryPlug->ngroup];
27
28			mpiSim = this;
29			}
30
31
32			mpiSimulation::~mpiSimulation(){
33
34			delete[] MolToProcMap;
35			delete[] MolComponentType;
36			delete[] AtomToProcMap;
37			delete[] GroupToProcMap;
38
39			delete parallelData;
40			// perhaps we should let fortran know the party is over.
41
42			}
43
44			void mpiSimulation::divideLabor( ){
45
46			int nComponents;
47			MoleculeStamp** compStamps;
48			randomSPRNG *myRandom;
49			int* componentsNmol;
50			int* AtomsPerProc;
51			int* GroupsPerProc;
52
53			double numerator;
54			double denominator;
55			double precast;
56			double x, y, a;
57			int old_atoms, add_atoms, new_atoms;
58			int old_groups, add_groups, new_groups;
59
60			int nTarget;
61			int molIndex, atomIndex, groupIndex;
62			int done;
63			int i, j, loops, which_proc;
64			int nmol_global, nmol_local;
65			int ngroups_global, ngroups_local;
66			int natoms_global, natoms_local;
67			int ncutoff_groups, nAtomsInGroups;
68			int local_index;
69			int baseSeed = entryPlug->getSeed();
70			CutoffGroupStamp* cg;
71
72			nComponents = entryPlug->nComponents;
73			compStamps = entryPlug->compStamps;
74			componentsNmol = entryPlug->componentsNmol;
75			AtomsPerProc = new int[parallelData->nProcessors];
76			GroupsPerProc = new int[parallelData->nProcessors];
77
78			parallelData->nAtomsGlobal = entryPlug->n_atoms;
79			parallelData->nBondsGlobal = entryPlug->n_bonds;
80			parallelData->nBendsGlobal = entryPlug->n_bends;
81			parallelData->nTorsionsGlobal = entryPlug->n_torsions;
82			parallelData->nSRIGlobal = entryPlug->n_SRI;
83			parallelData->nGroupsGlobal = entryPlug->ngroup;
84			parallelData->nMolGlobal = entryPlug->n_mol;
85
86			if (parallelData->nProcessors > parallelData->nMolGlobal) {
87			sprintf( painCave.errMsg,
88			"nProcessors (%d) > nMol (%d)\n"
89			"\tThe number of processors is larger than\n"
90			"\tthe number of molecules. This will not result in a \n"
91			"\tusable division of atoms for force decomposition.\n"
92			"\tEither try a smaller number of processors, or run the\n"
93			"\tsingle-processor version of OOPSE.\n",
94			parallelData->nProcessors, parallelData->nMolGlobal );
95			painCave.isFatal = 1;
96			simError();
97			}
98
99			myRandom = new randomSPRNG( baseSeed );
100
101
102			a = 3.0 * (double)parallelData->nMolGlobal / (double)parallelData->nAtomsGlobal;
103
104			// Initialize things that we'll send out later:
105			for (i = 0; i < parallelData->nProcessors; i++ ) {
106			AtomsPerProc[i] = 0;
107			GroupsPerProc[i] = 0;
108			}
109			for (i = 0; i < parallelData->nMolGlobal; i++ ) {
110			// default to an error condition:
111			MolToProcMap[i] = -1;
112			MolComponentType[i] = -1;
113			}
114			for (i = 0; i < parallelData->nAtomsGlobal; i++ ) {
115			// default to an error condition:
116			AtomToProcMap[i] = -1;
117			}
118			for (i = 0; i < parallelData->nGroupsGlobal; i++ ) {
119			// default to an error condition:
120			GroupToProcMap[i] = -1;
121			}
122
123			if (parallelData->myNode == 0) {
124			numerator = (double) entryPlug->n_atoms;
125			denominator = (double) parallelData->nProcessors;
126			precast = numerator / denominator;
127			nTarget = (int)( precast + 0.5 );
128
129			// Build the array of molecule component types first
130			molIndex = 0;
131			for (i=0; i < nComponents; i++) {
132			for (j=0; j < componentsNmol[i]; j++) {
133			MolComponentType[molIndex] = i;
134			molIndex++;
135			}
136			}
137
138			atomIndex = 0;
139			groupIndex = 0;
140
141			for (i = 0; i < molIndex; i++ ) {
142
143			done = 0;
144			loops = 0;
145
146			while( !done ){
147			loops++;
148
149			// Pick a processor at random
150
151			which_proc = (int) (myRandom->getRandom() * parallelData->nProcessors);
152
153			// How many atoms does this processor have?
154
155			old_atoms = AtomsPerProc[which_proc];
156			add_atoms = compStamps[MolComponentType[i]]->getNAtoms();
157			new_atoms = old_atoms + add_atoms;
158
159			old_groups = GroupsPerProc[which_proc];
160			ncutoff_groups = compStamps[MolComponentType[i]]->getNCutoffGroups();
161			nAtomsInGroups = 0;
162			for (j = 0; j < ncutoff_groups; j++) {
163			cg = compStamps[MolComponentType[i]]->getCutoffGroup(j);
164			nAtomsInGroups += cg->getNMembers();
165			}
166			add_groups = add_atoms - nAtomsInGroups + ncutoff_groups;
167			new_groups = old_groups + add_groups;
168
169			// If we've been through this loop too many times, we need
170			// to just give up and assign the molecule to this processor
171			// and be done with it.
172
173			if (loops > 100) {
174			sprintf( painCave.errMsg,
175			"I've tried 100 times to assign molecule %d to a "
176			" processor, but can't find a good spot.\n"
177			"I'm assigning it at random to processor %d.\n",
178			i, which_proc);
179			painCave.isFatal = 0;
180			simError();
181
182			MolToProcMap[i] = which_proc;
183			AtomsPerProc[which_proc] += add_atoms;
184			for (j = 0 ; j < add_atoms; j++ ) {
185			AtomToProcMap[atomIndex] = which_proc;
186			atomIndex++;
187			}
188			GroupsPerProc[which_proc] += add_groups;
189			for (j=0; j < add_groups; j++) {
190			GroupToProcMap[groupIndex] = which_proc;
191			groupIndex++;
192			}
193			done = 1;
194			continue;
195			}
196
197			// If we can add this molecule to this processor without sending
198			// it above nTarget, then go ahead and do it:
199
200			if (new_atoms <= nTarget) {
201			MolToProcMap[i] = which_proc;
202			AtomsPerProc[which_proc] += add_atoms;
203			for (j = 0 ; j < add_atoms; j++ ) {
204			AtomToProcMap[atomIndex] = which_proc;
205			atomIndex++;
206			}
207			GroupsPerProc[which_proc] += add_groups;
208			for (j=0; j < add_groups; j++) {
209			GroupToProcMap[groupIndex] = which_proc;
210			groupIndex++;
211			}
212			done = 1;
213			continue;
214			}
215
216
217			// The only situation left is when new_atoms > nTarget. We
218			// want to accept this with some probability that dies off the
219			// farther we are from nTarget
220
221			// roughly: x = new_atoms - nTarget
222			// Pacc(x) = exp(- a * x)
223			// where a = penalty / (average atoms per molecule)
224
225			x = (double) (new_atoms - nTarget);
226			y = myRandom->getRandom();
227
228			if (y < exp(- a * x)) {
229			MolToProcMap[i] = which_proc;
230			AtomsPerProc[which_proc] += add_atoms;
231			for (j = 0 ; j < add_atoms; j++ ) {
232			AtomToProcMap[atomIndex] = which_proc;
233			atomIndex++;
234			}
235			GroupsPerProc[which_proc] += add_groups;
236			for (j=0; j < add_groups; j++) {
237			GroupToProcMap[groupIndex] = which_proc;
238			groupIndex++;
239			}
240			done = 1;
241			continue;
242			} else {
243			continue;
244			}
245
246			}
247			}
248
249
250			// Spray out this nonsense to all other processors:
251
252			//std::cerr << "node 0 mol2proc = \n";
253			//for (i = 0; i < parallelData->nMolGlobal; i++)
254			// std::cerr << i << "\t" << MolToProcMap[i] << "\n";
255
256			MPI_Bcast(MolToProcMap, parallelData->nMolGlobal,
257			MPI_INT, 0, MPI_COMM_WORLD);
258
259			MPI_Bcast(AtomToProcMap, parallelData->nAtomsGlobal,
260			MPI_INT, 0, MPI_COMM_WORLD);
261
262			MPI_Bcast(GroupToProcMap, parallelData->nGroupsGlobal,
263			MPI_INT, 0, MPI_COMM_WORLD);
264
265			MPI_Bcast(MolComponentType, parallelData->nMolGlobal,
266			MPI_INT, 0, MPI_COMM_WORLD);
267
268			MPI_Bcast(AtomsPerProc, parallelData->nProcessors,
269			MPI_INT, 0, MPI_COMM_WORLD);
270
271			MPI_Bcast(GroupsPerProc, parallelData->nProcessors,
272			MPI_INT, 0, MPI_COMM_WORLD);
273			} else {
274
275			// Listen to your marching orders from processor 0:
276
277			MPI_Bcast(MolToProcMap, parallelData->nMolGlobal,
278			MPI_INT, 0, MPI_COMM_WORLD);
279
280			MPI_Bcast(AtomToProcMap, parallelData->nAtomsGlobal,
281			MPI_INT, 0, MPI_COMM_WORLD);
282
283			MPI_Bcast(GroupToProcMap, parallelData->nGroupsGlobal,
284			MPI_INT, 0, MPI_COMM_WORLD);
285
286			MPI_Bcast(MolComponentType, parallelData->nMolGlobal,
287			MPI_INT, 0, MPI_COMM_WORLD);
288
289			MPI_Bcast(AtomsPerProc, parallelData->nProcessors,
290			MPI_INT, 0, MPI_COMM_WORLD);
291
292			MPI_Bcast(GroupsPerProc, parallelData->nProcessors,
293			MPI_INT, 0, MPI_COMM_WORLD);
294
295
296			}
297
298			// Let's all check for sanity:
299
300			nmol_local = 0;
301			for (i = 0 ; i < parallelData->nMolGlobal; i++ ) {
302			if (MolToProcMap[i] == parallelData->myNode) {
303			nmol_local++;
304			}
305			}
306
307			natoms_local = 0;
308			for (i = 0; i < parallelData->nAtomsGlobal; i++) {
309			if (AtomToProcMap[i] == parallelData->myNode) {
310			natoms_local++;
311			}
312			}
313
314			ngroups_local = 0;
315			for (i = 0; i < parallelData->nGroupsGlobal; i++) {
316			if (GroupToProcMap[i] == parallelData->myNode) {
317			ngroups_local++;
318			}
319			}
320
321			MPI_Allreduce(&nmol_local,&nmol_global,1,MPI_INT,MPI_SUM,
322			MPI_COMM_WORLD);
323
324			MPI_Allreduce(&natoms_local,&natoms_global,1,MPI_INT,
325			MPI_SUM, MPI_COMM_WORLD);
326
327			MPI_Allreduce(&ngroups_local,&ngroups_global,1,MPI_INT,
328			MPI_SUM, MPI_COMM_WORLD);
329
330			if( nmol_global != entryPlug->n_mol ){
331			sprintf( painCave.errMsg,
332			"The sum of all nmol_local, %d, did not equal the "
333			"total number of molecules, %d.\n",
334			nmol_global, entryPlug->n_mol );
335			painCave.isFatal = 1;
336			simError();
337			}
338
339			if( natoms_global != entryPlug->n_atoms ){
340			sprintf( painCave.errMsg,
341			"The sum of all natoms_local, %d, did not equal the "
342			"total number of atoms, %d.\n",
343			natoms_global, entryPlug->n_atoms );
344			painCave.isFatal = 1;
345			simError();
346			}
347
348			if( ngroups_global != entryPlug->ngroup ){
349			sprintf( painCave.errMsg,
350			"The sum of all ngroups_local, %d, did not equal the "
351			"total number of cutoffGroups, %d.\n",
352			ngroups_global, entryPlug->ngroup );
353			painCave.isFatal = 1;
354			simError();
355			}
356
357			sprintf( checkPointMsg,
358			"Successfully divided the molecules among the processors.\n" );
359			MPIcheckPoint();
360
361			parallelData->nMolLocal = nmol_local;
362			parallelData->nAtomsLocal = natoms_local;
363			parallelData->nGroupsLocal = ngroups_local;
364
365			globalAtomIndex.resize(parallelData->nAtomsLocal);
366			globalToLocalAtom.resize(parallelData->nAtomsGlobal);
367			local_index = 0;
368			for (i = 0; i < parallelData->nAtomsGlobal; i++) {
369			if (AtomToProcMap[i] == parallelData->myNode) {
370			globalAtomIndex[local_index] = i;
371
372			globalToLocalAtom[i] = local_index;
373			local_index++;
374
375			}
376			else
377			globalToLocalAtom[i] = -1;
378			}
379
380			globalGroupIndex.resize(parallelData->nGroupsLocal);
381			globalToLocalGroup.resize(parallelData->nGroupsGlobal);
382			local_index = 0;
383			for (i = 0; i < parallelData->nGroupsGlobal; i++) {
384			if (GroupToProcMap[i] == parallelData->myNode) {
385			globalGroupIndex[local_index] = i;
386
387			globalToLocalGroup[i] = local_index;
388			local_index++;
389
390			}
391			else
392			globalToLocalGroup[i] = -1;
393			}
394
395			globalMolIndex.resize(parallelData->nMolLocal);
396			globalToLocalMol.resize(parallelData->nMolGlobal);
397			local_index = 0;
398			for (i = 0; i < parallelData->nMolGlobal; i++) {
399			if (MolToProcMap[i] == parallelData->myNode) {
400			globalMolIndex[local_index] = i;
401			globalToLocalMol[i] = local_index;
402			local_index++;
403			}
404			else
405			globalToLocalMol[i] = -1;
406			}
407
408			}
409
410
411			void mpiSimulation::mpiRefresh( void ){
412
413			int isError, i;
414			int *localToGlobalAtomIndex = new int[parallelData->nAtomsLocal];
415			int *localToGlobalGroupIndex = new int[parallelData->nGroupsLocal];
416
417			// Fortran indexing needs to be increased by 1 in order to get the 2
418			// languages to not barf
419
420			for(i = 0; i < parallelData->nAtomsLocal; i++)
421			localToGlobalAtomIndex[i] = globalAtomIndex[i] + 1;
422
423			for(i = 0; i < parallelData->nGroupsLocal; i++)
424			localToGlobalGroupIndex[i] = globalGroupIndex[i] + 1;
425
426			isError = 0;
427
428			setFsimParallel( parallelData,
429			&(parallelData->nAtomsLocal), localToGlobalAtomIndex,
430			&(parallelData->nGroupsLocal), localToGlobalGroupIndex,
431			&isError );
432
433			if( isError ){
434
435			sprintf( painCave.errMsg,
436			"mpiRefresh errror: fortran didn't like something we gave it.\n" );
437			painCave.isFatal = 1;
438			simError();
439			}
440
441			delete[] localToGlobalGroupIndex;
442			delete[] localToGlobalAtomIndex;
443
444
445			sprintf( checkPointMsg,
446			" mpiRefresh successful.\n" );
447			MPIcheckPoint();
448			}
449
450
451			#endif // is_mpi