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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 1093 by tim, Wed Mar 17 14:22:59 2004 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

# Line 9 | Line 9
9   #include "parse_me.h"
10   #include "Integrator.hpp"
11   #include "simError.h"
12 + #include "RigidBody.hpp"
13   //#include "ConjugateMinimizer.hpp"
14   #include "OOPSEMinimizer.hpp"
15  
# Line 165 | Line 166 | void SimSetup::makeMolecules(void){
166  
167  
168   void SimSetup::makeMolecules(void){
169 <  int k;
170 <  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
169 >  int i, j, k;
170 >  int exI, exJ, exK, exL, slI;
171 >  int tempI, tempJ, tempK, tempL;
172 >  int molI;
173 >  int stampID, atomOffset, rbOffset;
174    molInit molInfo;
175    DirectionalAtom* dAtom;
176 +  RigidBody* myRB;
177 +  StuntDouble* mySD;
178    LinkedAssign* extras;
179    LinkedAssign* current_extra;
180    AtomStamp* currentAtom;
181    BondStamp* currentBond;
182    BendStamp* currentBend;
183    TorsionStamp* currentTorsion;
184 +  RigidBodyStamp* currentRigidBody;
185  
186    bond_pair* theBonds;
187    bend_set* theBends;
188    torsion_set* theTorsions;
189  
190 +  set<int> skipList;
191 +
192 +  double phi, theta, psi;
193 +
194    //init the forceField paramters
195  
196    the_ff->readParams();
197  
187 –
198    // init the atoms
199  
200 <  double phi, theta, psi;
191 <  double sux, suy, suz;
192 <  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
193 <  double ux, uy, uz, u, uSqr;
200 >  int nMembers, nNew, rb1, rb2;
201  
202    for (k = 0; k < nInfo; k++){
203      the_ff->setSimInfo(&(info[k]));
204  
205      atomOffset = 0;
206 <    excludeOffset = 0;
206 >
207      for (i = 0; i < info[k].n_mol; i++){
208        stampID = info[k].molecules[i].getStampID();
209  
# Line 204 | Line 211 | void SimSetup::makeMolecules(void){
211        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
212        molInfo.nBends = comp_stamps[stampID]->getNBends();
213        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
214 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
215 <
214 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
215 >      
216        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
210 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
211 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
212 –      molInfo.myBends = new Bend * [molInfo.nBends];
213 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
217  
218 +      if (molInfo.nBonds > 0)
219 +        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
220 +      else
221 +        molInfo.myBonds = NULL;
222 +
223 +      if (molInfo.nBends > 0)
224 +        molInfo.myBends = new (Bend *) [molInfo.nBends];
225 +      else
226 +        molInfo.myBends = NULL;
227 +
228 +      if (molInfo.nTorsions > 0)
229 +        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
230 +      else
231 +        molInfo.myTorsions = NULL;
232 +
233        theBonds = new bond_pair[molInfo.nBonds];
234        theBends = new bend_set[molInfo.nBends];
235        theTorsions = new torsion_set[molInfo.nTorsions];
236 <
236 >      
237        // make the Atoms
238  
239        for (j = 0; j < molInfo.nAtoms; j++){
240          currentAtom = comp_stamps[stampID]->getAtom(j);
241 +
242          if (currentAtom->haveOrientation()){
243            dAtom = new DirectionalAtom((j + atomOffset),
244                                        info[k].getConfiguration());
# Line 233 | Line 252 | void SimSetup::makeMolecules(void){
252            phi = currentAtom->getEulerPhi() * M_PI / 180.0;
253            theta = currentAtom->getEulerTheta() * M_PI / 180.0;
254            psi = currentAtom->getEulerPsi()* M_PI / 180.0;
236 –            
237 –          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
238 –          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
239 –          Axz = sin(theta) * sin(psi);
240 –          
241 –          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
242 –          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
243 –          Ayz = sin(theta) * cos(psi);
244 –          
245 –          Azx = sin(phi) * sin(theta);
246 –          Azy = -cos(phi) * sin(theta);
247 –          Azz = cos(theta);
255  
256 <          sux = 0.0;
257 <          suy = 0.0;
251 <          suz = 1.0;
252 <
253 <          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
254 <          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
255 <          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
256 <
257 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
258 <
259 <          u = sqrt(uSqr);
260 <          ux = ux / u;
261 <          uy = uy / u;
262 <          uz = uz / u;
263 <
264 <          dAtom->setSUx(ux);
265 <          dAtom->setSUy(uy);
266 <          dAtom->setSUz(uz);
256 >          dAtom->setUnitFrameFromEuler(phi, theta, psi);
257 >            
258          }
259          else{
260 <          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
261 <                                               info[k].getConfiguration());
260 >
261 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
262          }
263 +
264          molInfo.myAtoms[j]->setType(currentAtom->getType());
265  
266   #ifdef IS_MPI
# Line 284 | Line 276 | void SimSetup::makeMolecules(void){
276          theBonds[j].a = currentBond->getA() + atomOffset;
277          theBonds[j].b = currentBond->getB() + atomOffset;
278  
279 <        exI = theBonds[j].a;
280 <        exJ = theBonds[j].b;
279 >        tempI = theBonds[j].a;
280 >        tempJ = theBonds[j].b;
281  
290 –        // exclude_I must always be the smaller of the pair
291 –        if (exI > exJ){
292 –          tempEx = exI;
293 –          exI = exJ;
294 –          exJ = tempEx;
295 –        }
282   #ifdef IS_MPI
283 <        tempEx = exI;
284 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
285 <        tempEx = exJ;
286 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
283 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
284 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
285 > #else
286 >        exI = tempI + 1;
287 >        exJ = tempJ + 1;
288 > #endif
289  
290 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
303 < #else  // isn't MPI
304 <
305 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
306 < #endif  //is_mpi
290 >        info[k].excludes->addPair(exI, exJ);
291        }
308 –      excludeOffset += molInfo.nBonds;
292  
293        //make the bends
294        for (j = 0; j < molInfo.nBends; j++){
# Line 355 | Line 338 | void SimSetup::makeMolecules(void){
338            }
339          }
340  
341 <        if (!theBends[j].isGhost){
342 <          exI = theBends[j].a;
343 <          exJ = theBends[j].c;
344 <        }
345 <        else{
363 <          exI = theBends[j].a;
364 <          exJ = theBends[j].b;
365 <        }
366 <
367 <        // exclude_I must always be the smaller of the pair
368 <        if (exI > exJ){
369 <          tempEx = exI;
370 <          exI = exJ;
371 <          exJ = tempEx;
372 <        }
341 >        if (theBends[j].isGhost) {
342 >          
343 >          tempI = theBends[j].a;
344 >          tempJ = theBends[j].b;
345 >          
346   #ifdef IS_MPI
347 <        tempEx = exI;
348 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
349 <        tempEx = exJ;
350 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
347 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
348 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
349 > #else
350 >          exI = tempI + 1;
351 >          exJ = tempJ + 1;
352 > #endif          
353 >          info[k].excludes->addPair(exI, exJ);
354  
355 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
356 < #else  // isn't MPI
357 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
358 < #endif  //is_mpi
355 >        } else {
356 >
357 >          tempI = theBends[j].a;
358 >          tempJ = theBends[j].b;
359 >          tempK = theBends[j].c;
360 >          
361 > #ifdef IS_MPI
362 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
363 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
364 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
365 > #else
366 >          exI = tempI + 1;
367 >          exJ = tempJ + 1;
368 >          exK = tempK + 1;
369 > #endif
370 >          
371 >          info[k].excludes->addPair(exI, exK);
372 >          info[k].excludes->addPair(exI, exJ);
373 >          info[k].excludes->addPair(exJ, exK);
374 >        }
375        }
384 –      excludeOffset += molInfo.nBends;
376  
377        for (j = 0; j < molInfo.nTorsions; j++){
378          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 390 | Line 381 | void SimSetup::makeMolecules(void){
381          theTorsions[j].c = currentTorsion->getC() + atomOffset;
382          theTorsions[j].d = currentTorsion->getD() + atomOffset;
383  
384 <        exI = theTorsions[j].a;
385 <        exJ = theTorsions[j].d;
384 >        tempI = theTorsions[j].a;      
385 >        tempJ = theTorsions[j].b;
386 >        tempK = theTorsions[j].c;
387 >        tempL = theTorsions[j].d;
388  
396 –        // exclude_I must always be the smaller of the pair
397 –        if (exI > exJ){
398 –          tempEx = exI;
399 –          exI = exJ;
400 –          exJ = tempEx;
401 –        }
389   #ifdef IS_MPI
390 <        tempEx = exI;
391 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
392 <        tempEx = exJ;
393 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
390 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
391 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
392 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
393 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
394 > #else
395 >        exI = tempI + 1;
396 >        exJ = tempJ + 1;
397 >        exK = tempK + 1;
398 >        exL = tempL + 1;
399 > #endif
400  
401 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
402 < #else  // isn't MPI
403 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
404 < #endif  //is_mpi
401 >        info[k].excludes->addPair(exI, exJ);
402 >        info[k].excludes->addPair(exI, exK);
403 >        info[k].excludes->addPair(exI, exL);        
404 >        info[k].excludes->addPair(exJ, exK);
405 >        info[k].excludes->addPair(exJ, exL);
406 >        info[k].excludes->addPair(exK, exL);
407        }
413 –      excludeOffset += molInfo.nTorsions;
408  
409 +      for (j = 0; j < molInfo.nRigidBodies; j++){
410  
411 <      // send the arrays off to the forceField for init.
411 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
412 >        nMembers = currentRigidBody->getNMembers();
413  
414 +        // Create the Rigid Body:
415 +
416 +        myRB = new RigidBody();
417 +        
418 +        for (rb1 = 0; rb1 < nMembers; rb1++) {
419 +
420 +          // molI is atom numbering inside this molecule
421 +          molI = currentRigidBody->getMember(rb1);    
422 +
423 +          // tempI is atom numbering on local processor
424 +          tempI = molI + atomOffset;
425 +
426 +          // currentAtom is the AtomStamp (which we need for
427 +          // rigid body reference positions)
428 +          currentAtom = comp_stamps[stampID]->getAtom(molI);
429 +
430 +          // When we add to the rigid body, add the atom itself and
431 +          // the stamp info:
432 +
433 +          myRB->addAtom(info[k].atoms[tempI], currentAtom);
434 +          
435 +          // Add this atom to the Skip List for the integrators
436 + #ifdef IS_MPI
437 +          slI = info[k].atoms[tempI]->getGlobalIndex();
438 + #else
439 +          slI = tempI;
440 + #endif
441 +          skipList.insert(slI);
442 +          
443 +        }
444 +        
445 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
446 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
447 +            
448 +            tempI = currentRigidBody->getMember(rb1);
449 +            tempJ = currentRigidBody->getMember(rb2);
450 +            
451 +            // Some explanation is required here.
452 +            // Fortran indexing starts at 1, while c indexing starts at 0
453 +            // Also, in parallel computations, the GlobalIndex is
454 +            // used for the exclude list:
455 +            
456 + #ifdef IS_MPI
457 +            exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
458 +            exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
459 + #else
460 +            exI = tempI + 1;
461 +            exJ = tempJ + 1;
462 + #endif
463 +            
464 +            info[k].excludes->addPair(exI, exJ);
465 +            
466 +          }
467 +        }
468 +      }
469 +      
470 +      // send the arrays off to the forceField for init.
471 +      
472        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
473        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
474        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
475        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
476                                   theTorsions);
477  
424 –
478        info[k].molecules[i].initialize(molInfo);
479  
480  
# Line 430 | Line 483 | void SimSetup::makeMolecules(void){
483        delete[] theBends;
484        delete[] theTorsions;
485      }
486 +
487 +    // build up the integrableObjects vector:
488 +
489 +    for (i = 0; i < info[k].n_atoms; i++) {
490 +      
491 + #ifdef IS_MPI
492 +      slI = info[k].atoms[i]->getGlobalIndex();
493 + #else
494 +      slI = i;
495 + #endif
496 +
497 +      if (skipList.find(slI) == skipList.end()) {
498 +        mySD = (StuntDouble *) info[k].atoms[i];
499 +        info[k].integrableObjects.push_back(mySD);
500 +      }
501 +    }
502 +    for (i = 0; i < info[k].rigidBodies.size(); i++) {
503 +      mySD = (StuntDouble *) info[k].rigidBodies[i];
504 +      info[k].integrableObjects.push_back(mySD);      
505 +    }
506 +    
507    }
508  
509   #ifdef IS_MPI
# Line 439 | Line 513 | void SimSetup::makeMolecules(void){
513  
514    // clean up the forcefield
515  
516 <  the_ff->calcRcut();
516 >  if (!globals->haveLJrcut()){
517 >
518 >    the_ff->calcRcut();
519 >
520 >  } else {
521 >    
522 >    the_ff->setRcut( globals->getLJrcut() );
523 >  }
524 >
525    the_ff->cleanMe();
526   }
527  
# Line 1221 | Line 1303 | void SimSetup::calcSysValues(void){
1303    tot_bonds = 0;
1304    tot_bends = 0;
1305    tot_torsions = 0;
1306 +  tot_rigid = 0;
1307    for (i = 0; i < n_components; i++){
1308      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1309      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1310      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1311      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1312 +    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1313    }
1314 <
1314 >  
1315    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1316    molMembershipArray = new int[tot_atoms];
1317  
# Line 1249 | Line 1333 | void SimSetup::mpiMolDivide(void){
1333    int i, j, k;
1334    int localMol, allMol;
1335    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1336 +  int local_rigid;
1337  
1338    mpiSim = new mpiSimulation(info);
1339  
# Line 1265 | Line 1350 | void SimSetup::mpiMolDivide(void){
1350    local_bonds = 0;
1351    local_bends = 0;
1352    local_torsions = 0;
1353 +  local_rigid = 0;
1354    globalAtomIndex = 0;
1355  
1270 –
1356    for (i = 0; i < n_components; i++){
1357      for (j = 0; j < components_nmol[i]; j++){
1358        if (mol2proc[allMol] == worldRank){
# Line 1275 | Line 1360 | void SimSetup::mpiMolDivide(void){
1360          local_bonds += comp_stamps[i]->getNBonds();
1361          local_bends += comp_stamps[i]->getNBends();
1362          local_torsions += comp_stamps[i]->getNTorsions();
1363 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1364          localMol++;
1365        }      
1366        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
# Line 1288 | Line 1374 | void SimSetup::mpiMolDivide(void){
1374    local_SRI = local_bonds + local_bends + local_torsions;
1375  
1376    info[0].n_atoms = mpiSim->getMyNlocal();  
1377 +  
1378  
1379    if (local_atoms != info[0].n_atoms){
1380      sprintf(painCave.errMsg,
# Line 1320 | Line 1407 | void SimSetup::makeSysArrays(void){
1407  
1408    Atom** the_atoms;
1409    Molecule* the_molecules;
1323 –  Exclude** the_excludes;
1410  
1325 –
1411    for (l = 0; l < nInfo; l++){
1412      // create the atom and short range interaction arrays
1413  
# Line 1365 | Line 1450 | void SimSetup::makeSysArrays(void){
1450  
1451   #endif // is_mpi
1452  
1453 <
1454 <    if (info[l].n_SRI){
1455 <      Exclude::createArray(info[l].n_SRI);
1371 <      the_excludes = new Exclude * [info[l].n_SRI];
1372 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1373 <        the_excludes[ex] = new Exclude(ex);
1374 <      }
1375 <      info[l].globalExcludes = new int;
1376 <      info[l].n_exclude = info[l].n_SRI;
1377 <    }
1378 <    else{
1379 <      Exclude::createArray(1);
1380 <      the_excludes = new Exclude * ;
1381 <      the_excludes[0] = new Exclude(0);
1382 <      the_excludes[0]->setPair(0, 0);
1383 <      info[l].globalExcludes = new int;
1384 <      info[l].globalExcludes[0] = 0;
1385 <      info[l].n_exclude = 0;
1386 <    }
1387 <
1453 >    info[l].globalExcludes = new int;
1454 >    info[l].globalExcludes[0] = 0;
1455 >    
1456      // set the arrays into the SimInfo object
1457  
1458      info[l].atoms = the_atoms;
1459      info[l].molecules = the_molecules;
1460      info[l].nGlobalExcludes = 0;
1393 –    info[l].excludes = the_excludes;
1461  
1462      the_ff->setSimInfo(info);
1463    }

Diff Legend

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