| 1 | 
  | 
#ifdef IS_MPI | 
| 2 | 
< | 
 | 
| 2 | 
> | 
#include <iostream> | 
| 3 | 
  | 
#include <cstdlib> | 
| 4 | 
  | 
#include <cstring> | 
| 5 | 
  | 
#include <cmath> | 
| 85 | 
  | 
 | 
| 86 | 
  | 
  myRandom = new randomSPRNG( baseSeed ); | 
| 87 | 
  | 
 | 
| 88 | 
< | 
  a = (double)mpiPlug->nMolGlobal / (double)mpiPlug->nAtomsGlobal; | 
| 88 | 
> | 
  a = 3.0 * (double)mpiPlug->nMolGlobal / (double)mpiPlug->nAtomsGlobal; | 
| 89 | 
  | 
 | 
| 90 | 
  | 
  // Initialize things that we'll send out later: | 
| 91 | 
  | 
  for (i = 0; i < mpiPlug->numberProcessors; i++ ) { | 
| 133 | 
  | 
        // How many atoms does this processor have? | 
| 134 | 
  | 
         | 
| 135 | 
  | 
        old_atoms = AtomsPerProc[which_proc]; | 
| 136 | 
– | 
 | 
| 137 | 
– | 
        // If the processor already had too many atoms, just skip this | 
| 138 | 
– | 
        // processor and try again. | 
| 139 | 
– | 
 | 
| 140 | 
– | 
        if (old_atoms >= nTarget) continue; | 
| 141 | 
– | 
 | 
| 136 | 
  | 
        add_atoms = compStamps[MolComponentType[i]]->getNAtoms(); | 
| 137 | 
  | 
        new_atoms = old_atoms + add_atoms; | 
| 144 | 
– | 
     | 
| 145 | 
– | 
        // If we can add this molecule to this processor without sending | 
| 146 | 
– | 
        // it above nTarget, then go ahead and do it: | 
| 147 | 
– | 
     | 
| 148 | 
– | 
        if (new_atoms <= nTarget) { | 
| 149 | 
– | 
          MolToProcMap[i] = which_proc; | 
| 150 | 
– | 
          AtomsPerProc[which_proc] += add_atoms; | 
| 151 | 
– | 
          for (j = 0 ; j < add_atoms; j++ ) { | 
| 152 | 
– | 
            atomIndex++; | 
| 153 | 
– | 
            AtomToProcMap[atomIndex] = which_proc; | 
| 154 | 
– | 
          } | 
| 155 | 
– | 
          done = 1; | 
| 156 | 
– | 
          continue; | 
| 157 | 
– | 
        } | 
| 138 | 
  | 
 | 
| 139 | 
  | 
        // If we've been through this loop too many times, we need | 
| 140 | 
  | 
        // to just give up and assign the molecule to this processor | 
| 152 | 
  | 
          MolToProcMap[i] = which_proc; | 
| 153 | 
  | 
          AtomsPerProc[which_proc] += add_atoms; | 
| 154 | 
  | 
          for (j = 0 ; j < add_atoms; j++ ) { | 
| 155 | 
< | 
            atomIndex++; | 
| 156 | 
< | 
            AtomToProcMap[atomIndex] = which_proc; | 
| 155 | 
> | 
            AtomToProcMap[atomIndex] = which_proc; | 
| 156 | 
> | 
            atomIndex++; | 
| 157 | 
> | 
          } | 
| 158 | 
> | 
          done = 1; | 
| 159 | 
> | 
          continue; | 
| 160 | 
> | 
        } | 
| 161 | 
> | 
     | 
| 162 | 
> | 
        // If we can add this molecule to this processor without sending | 
| 163 | 
> | 
        // it above nTarget, then go ahead and do it: | 
| 164 | 
> | 
     | 
| 165 | 
> | 
        if (new_atoms <= nTarget) { | 
| 166 | 
> | 
          MolToProcMap[i] = which_proc; | 
| 167 | 
> | 
          AtomsPerProc[which_proc] += add_atoms; | 
| 168 | 
> | 
          for (j = 0 ; j < add_atoms; j++ ) { | 
| 169 | 
> | 
            AtomToProcMap[atomIndex] = which_proc; | 
| 170 | 
> | 
            atomIndex++; | 
| 171 | 
  | 
          } | 
| 172 | 
  | 
          done = 1; | 
| 173 | 
  | 
          continue; | 
| 174 | 
  | 
        } | 
| 175 | 
  | 
 | 
| 182 | 
– | 
        // The only situation left is where old_atoms < nTarget, but | 
| 183 | 
– | 
        // new_atoms > nTarget.   We want to accept this with some | 
| 184 | 
– | 
        // probability that dies off the farther we are from nTarget | 
| 176 | 
  | 
 | 
| 177 | 
+ | 
        // The only situation left is when new_atoms > nTarget.  We | 
| 178 | 
+ | 
        // want to accept this with some probability that dies off the | 
| 179 | 
+ | 
        // farther we are from nTarget | 
| 180 | 
+ | 
 | 
| 181 | 
  | 
        // roughly:  x = new_atoms - nTarget | 
| 182 | 
  | 
        //           Pacc(x) = exp(- a * x) | 
| 183 | 
< | 
        // where a = 1 / (average atoms per molecule) | 
| 183 | 
> | 
        // where a = penalty / (average atoms per molecule) | 
| 184 | 
  | 
 | 
| 185 | 
  | 
        x = (double) (new_atoms - nTarget); | 
| 186 | 
  | 
        y = myRandom->getRandom(); | 
| 187 | 
< | 
         | 
| 188 | 
< | 
        if (exp(- a * x) > y) { | 
| 187 | 
> | 
       | 
| 188 | 
> | 
        if (y < exp(- a * x)) { | 
| 189 | 
  | 
          MolToProcMap[i] = which_proc; | 
| 190 | 
  | 
          AtomsPerProc[which_proc] += add_atoms; | 
| 191 | 
  | 
          for (j = 0 ; j < add_atoms; j++ ) { | 
| 192 | 
< | 
            atomIndex++; | 
| 193 | 
< | 
            AtomToProcMap[atomIndex] = which_proc; | 
| 194 | 
< | 
          } | 
| 192 | 
> | 
            AtomToProcMap[atomIndex] = which_proc; | 
| 193 | 
> | 
            atomIndex++; | 
| 194 | 
> | 
           } | 
| 195 | 
  | 
          done = 1; | 
| 196 | 
  | 
          continue; | 
| 197 | 
  | 
        } else { | 
| 229 | 
  | 
     | 
| 230 | 
  | 
    MPI::COMM_WORLD.Bcast(AtomsPerProc, mpiPlug->numberProcessors, | 
| 231 | 
  | 
                          MPI_INT, 0); | 
| 232 | 
+ | 
 | 
| 233 | 
+ | 
 | 
| 234 | 
  | 
  } | 
| 235 | 
  | 
 | 
| 236 | 
  | 
 | 
| 282 | 
  | 
  local_index = 0; | 
| 283 | 
  | 
  for (i = 0; i < mpiPlug->nAtomsGlobal; i++) { | 
| 284 | 
  | 
    if (AtomToProcMap[i] == mpiPlug->myNode) { | 
| 288 | 
– | 
      local_index++; | 
| 285 | 
  | 
      globalIndex[local_index] = i; | 
| 286 | 
+ | 
      local_index++; | 
| 287 | 
  | 
    } | 
| 288 | 
  | 
  } | 
| 289 | 
< | 
  | 
| 289 | 
> | 
   | 
| 290 | 
  | 
  return globalIndex; | 
| 291 | 
  | 
} | 
| 292 | 
  | 
 | 
| 296 | 
  | 
  int isError, i; | 
| 297 | 
  | 
  int *globalIndex = new int[mpiPlug->myNlocal]; | 
| 298 | 
  | 
 | 
| 299 | 
< | 
  for(i=0; i<mpiPlug->myNlocal; i++) globalIndex[i] = entryPlug->atoms[i]->getGlobalIndex(); | 
| 299 | 
> | 
  // Fortran indexing needs to be increased by 1 in order to get the 2 languages to | 
| 300 | 
> | 
  // not barf | 
| 301 | 
  | 
 | 
| 302 | 
+ | 
  for(i=0; i<mpiPlug->myNlocal; i++) globalIndex[i] = entryPlug->atoms[i]->getGlobalIndex()+1; | 
| 303 | 
+ | 
 | 
| 304 | 
  | 
   | 
| 305 | 
  | 
  isError = 0; | 
| 306 | 
  | 
  setFsimParallel( mpiPlug, &(entryPlug->n_atoms), globalIndex, &isError ); |