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 983 by gezelter, Mon Jan 26 21:45:03 2004 UTC vs.
Revision 1157 by tim, Tue May 11 20:33:41 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  
16   #ifdef IS_MPI
17   #include "mpiBASS.h"
# Line 24 | Line 27
27   #define NPTxyz_ENS     4
28  
29  
30 < #define FF_DUFF 0
31 < #define FF_LJ   1
32 < #define FF_EAM  2
30 > #define FF_DUFF  0
31 > #define FF_LJ    1
32 > #define FF_EAM   2
33 > #define FF_H2O   3
34  
35   using namespace std;
36  
# Line 143 | Line 147 | void SimSetup::createSim(void){
147    // make the output filenames
148  
149    makeOutNames();
150 <
147 <  // make the integrator
148 <
149 <  makeIntegrator();
150 <
150 >  
151   #ifdef IS_MPI
152    mpiSim->mpiRefresh();
153   #endif
# Line 155 | Line 155 | void SimSetup::createSim(void){
155    // initialize the Fortran
156  
157    initFortran();
158 +
159 +  if (globals->haveMinimizer())
160 +    // make minimizer
161 +    makeMinimizer();
162 +  else
163 +    // make the integrator
164 +    makeIntegrator();
165 +
166   }
167  
168  
169   void SimSetup::makeMolecules(void){
170 <  int k;
171 <  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
170 >  int i, j, k;
171 >  int exI, exJ, exK, exL, slI, slJ;
172 >  int tempI, tempJ, tempK, tempL;
173 >  int molI;
174 >  int stampID, atomOffset, rbOffset;
175    molInit molInfo;
176    DirectionalAtom* dAtom;
177 +  RigidBody* myRB;
178 +  StuntDouble* mySD;
179    LinkedAssign* extras;
180    LinkedAssign* current_extra;
181    AtomStamp* currentAtom;
182    BondStamp* currentBond;
183    BendStamp* currentBend;
184    TorsionStamp* currentTorsion;
185 <
185 >  RigidBodyStamp* currentRigidBody;
186 >  CutoffGroupStamp* currentCutoffGroup;
187 >  CutoffGroup* myCutoffGroup;
188 >  
189    bond_pair* theBonds;
190    bend_set* theBends;
191    torsion_set* theTorsions;
192  
193 +  set<int> skipList;
194 +
195 +  double phi, theta, psi;
196 +  char* molName;
197 +  char rbName[100];
198 +
199    //init the forceField paramters
200  
201    the_ff->readParams();
202  
181 –
203    // init the atoms
204  
205 <  double phi, theta, psi;
185 <  double sux, suy, suz;
186 <  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
187 <  double ux, uy, uz, u, uSqr;
205 >  int nMembers, nNew, rb1, rb2;
206  
207    for (k = 0; k < nInfo; k++){
208      the_ff->setSimInfo(&(info[k]));
209  
210      atomOffset = 0;
211 <    excludeOffset = 0;
211 >
212      for (i = 0; i < info[k].n_mol; i++){
213        stampID = info[k].molecules[i].getStampID();
214 +      molName = comp_stamps[stampID]->getID();
215  
216        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
217        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
218        molInfo.nBends = comp_stamps[stampID]->getNBends();
219        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
220 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
221 <
220 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
221 >      molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
222 >      
223        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
204 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
205 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
206 –      molInfo.myBends = new Bend * [molInfo.nBends];
207 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
224  
225 +      if (molInfo.nBonds > 0)
226 +        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
227 +      else
228 +        molInfo.myBonds = NULL;
229 +
230 +      if (molInfo.nBends > 0)
231 +        molInfo.myBends = new (Bend *) [molInfo.nBends];
232 +      else
233 +        molInfo.myBends = NULL;
234 +
235 +      if (molInfo.nTorsions > 0)
236 +        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
237 +      else
238 +        molInfo.myTorsions = NULL;
239 +
240        theBonds = new bond_pair[molInfo.nBonds];
241        theBends = new bend_set[molInfo.nBends];
242        theTorsions = new torsion_set[molInfo.nTorsions];
243 <
243 >      
244        // make the Atoms
245  
246        for (j = 0; j < molInfo.nAtoms; j++){
247          currentAtom = comp_stamps[stampID]->getAtom(j);
248 +
249          if (currentAtom->haveOrientation()){
250            dAtom = new DirectionalAtom((j + atomOffset),
251                                        info[k].getConfiguration());
# Line 224 | Line 256 | void SimSetup::makeMolecules(void){
256            // in space using the three Euler angles.  We assume the standard
257            // unit vector was originally along the z axis below.
258  
259 <          phi = currentAtom->getEulerPhi();
260 <          theta = currentAtom->getEulerTheta();
261 <          psi = currentAtom->getEulerPsi();
259 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
260 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
261 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
262 >
263 >          dAtom->setUnitFrameFromEuler(phi, theta, psi);
264              
265 <          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
266 <          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
233 <          Axz = sin(theta) * sin(psi);
234 <          
235 <          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
236 <          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
237 <          Ayz = sin(theta) * cos(psi);
238 <          
239 <          Azx = sin(phi) * sin(theta);
240 <          Azy = -cos(phi) * sin(theta);
241 <          Azz = cos(theta);
265 >        }
266 >        else{
267  
268 <          sux = 0.0;
244 <          suy = 0.0;
245 <          suz = 1.0;
268 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
269  
247 –          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
248 –          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
249 –          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
250 –
251 –          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
252 –
253 –          u = sqrt(uSqr);
254 –          ux = ux / u;
255 –          uy = uy / u;
256 –          uz = uz / u;
257 –
258 –          dAtom->setSUx(ux);
259 –          dAtom->setSUy(uy);
260 –          dAtom->setSUz(uz);
270          }
262 –        else{
263 –          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
264 –                                               info[k].getConfiguration());
265 –        }
266 –        molInfo.myAtoms[j]->setType(currentAtom->getType());
271  
272 +        molInfo.myAtoms[j]->setType(currentAtom->getType());
273   #ifdef IS_MPI
274  
275 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
275 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
276  
277   #endif // is_mpi
278        }
# Line 278 | Line 283 | void SimSetup::makeMolecules(void){
283          theBonds[j].a = currentBond->getA() + atomOffset;
284          theBonds[j].b = currentBond->getB() + atomOffset;
285  
286 <        exI = theBonds[j].a;
287 <        exJ = theBonds[j].b;
286 >        tempI = theBonds[j].a;
287 >        tempJ = theBonds[j].b;
288  
284 –        // exclude_I must always be the smaller of the pair
285 –        if (exI > exJ){
286 –          tempEx = exI;
287 –          exI = exJ;
288 –          exJ = tempEx;
289 –        }
289   #ifdef IS_MPI
290 <        tempEx = exI;
291 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
292 <        tempEx = exJ;
293 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
290 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
291 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
292 > #else
293 >        exI = tempI + 1;
294 >        exJ = tempJ + 1;
295 > #endif
296  
297 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
297 < #else  // isn't MPI
298 <
299 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
300 < #endif  //is_mpi
297 >        info[k].excludes->addPair(exI, exJ);
298        }
302 –      excludeOffset += molInfo.nBonds;
299  
300        //make the bends
301        for (j = 0; j < molInfo.nBends; j++){
# Line 349 | Line 345 | void SimSetup::makeMolecules(void){
345            }
346          }
347  
348 <        if (!theBends[j].isGhost){
349 <          exI = theBends[j].a;
350 <          exJ = theBends[j].c;
351 <        }
352 <        else{
353 <          exI = theBends[j].a;
354 <          exJ = theBends[j].b;
355 <        }
348 >        if (theBends[j].isGhost) {
349 >          
350 >          tempI = theBends[j].a;
351 >          tempJ = theBends[j].b;
352 >          
353 > #ifdef IS_MPI
354 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
355 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
356 > #else
357 >          exI = tempI + 1;
358 >          exJ = tempJ + 1;
359 > #endif          
360 >          info[k].excludes->addPair(exI, exJ);
361  
362 <        // exclude_I must always be the smaller of the pair
362 <        if (exI > exJ){
363 <          tempEx = exI;
364 <          exI = exJ;
365 <          exJ = tempEx;
366 <        }
367 < #ifdef IS_MPI
368 <        tempEx = exI;
369 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
370 <        tempEx = exJ;
371 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
362 >        } else {
363  
364 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
365 < #else  // isn't MPI
366 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
367 < #endif  //is_mpi
364 >          tempI = theBends[j].a;
365 >          tempJ = theBends[j].b;
366 >          tempK = theBends[j].c;
367 >          
368 > #ifdef IS_MPI
369 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
370 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
371 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
372 > #else
373 >          exI = tempI + 1;
374 >          exJ = tempJ + 1;
375 >          exK = tempK + 1;
376 > #endif
377 >          
378 >          info[k].excludes->addPair(exI, exK);
379 >          info[k].excludes->addPair(exI, exJ);
380 >          info[k].excludes->addPair(exJ, exK);
381 >        }
382        }
378 –      excludeOffset += molInfo.nBends;
383  
384        for (j = 0; j < molInfo.nTorsions; j++){
385          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 384 | Line 388 | void SimSetup::makeMolecules(void){
388          theTorsions[j].c = currentTorsion->getC() + atomOffset;
389          theTorsions[j].d = currentTorsion->getD() + atomOffset;
390  
391 <        exI = theTorsions[j].a;
392 <        exJ = theTorsions[j].d;
391 >        tempI = theTorsions[j].a;      
392 >        tempJ = theTorsions[j].b;
393 >        tempK = theTorsions[j].c;
394 >        tempL = theTorsions[j].d;
395  
396 <        // exclude_I must always be the smaller of the pair
397 <        if (exI > exJ){
398 <          tempEx = exI;
399 <          exI = exJ;
400 <          exJ = tempEx;
396 > #ifdef IS_MPI
397 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
398 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
399 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
400 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
401 > #else
402 >        exI = tempI + 1;
403 >        exJ = tempJ + 1;
404 >        exK = tempK + 1;
405 >        exL = tempL + 1;
406 > #endif
407 >
408 >        info[k].excludes->addPair(exI, exJ);
409 >        info[k].excludes->addPair(exI, exK);
410 >        info[k].excludes->addPair(exI, exL);        
411 >        info[k].excludes->addPair(exJ, exK);
412 >        info[k].excludes->addPair(exJ, exL);
413 >        info[k].excludes->addPair(exK, exL);
414 >      }
415 >
416 >      
417 >      molInfo.myRigidBodies.clear();
418 >      
419 >      for (j = 0; j < molInfo.nRigidBodies; j++){
420 >
421 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
422 >        nMembers = currentRigidBody->getNMembers();
423 >
424 >        // Create the Rigid Body:
425 >
426 >        myRB = new RigidBody();
427 >
428 >        sprintf(rbName,"%s_RB_%d", molName, j);
429 >        myRB->setType(rbName);
430 >        
431 >        for (rb1 = 0; rb1 < nMembers; rb1++) {
432 >
433 >          // molI is atom numbering inside this molecule
434 >          molI = currentRigidBody->getMember(rb1);    
435 >
436 >          // tempI is atom numbering on local processor
437 >          tempI = molI + atomOffset;
438 >
439 >          // currentAtom is the AtomStamp (which we need for
440 >          // rigid body reference positions)
441 >          currentAtom = comp_stamps[stampID]->getAtom(molI);
442 >
443 >          // When we add to the rigid body, add the atom itself and
444 >          // the stamp info:
445 >
446 >          myRB->addAtom(info[k].atoms[tempI], currentAtom);
447 >          
448 >          // Add this atom to the Skip List for the integrators
449 > #ifdef IS_MPI
450 >          slI = info[k].atoms[tempI]->getGlobalIndex();
451 > #else
452 >          slI = tempI;
453 > #endif
454 >          skipList.insert(slI);
455 >          
456          }
457 +        
458 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
459 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
460 +            
461 +            tempI = currentRigidBody->getMember(rb1);
462 +            tempJ = currentRigidBody->getMember(rb2);
463 +            
464 +            // Some explanation is required here.
465 +            // Fortran indexing starts at 1, while c indexing starts at 0
466 +            // Also, in parallel computations, the GlobalIndex is
467 +            // used for the exclude list:
468 +            
469   #ifdef IS_MPI
470 <        tempEx = exI;
471 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
472 <        tempEx = exJ;
473 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
470 >            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
471 >            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
472 > #else
473 >            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
474 >            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
475 > #endif
476 >            
477 >            info[k].excludes->addPair(exI, exJ);
478 >            
479 >          }
480 >        }
481  
482 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
483 < #else  // isn't MPI
404 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
405 < #endif  //is_mpi
482 >        molInfo.myRigidBodies.push_back(myRB);
483 >        info[k].rigidBodies.push_back(myRB);
484        }
485 <      excludeOffset += molInfo.nTorsions;
485 >      
486  
487 +      //creat cutoff group for molecule
488 +      molInfo.myCutoffGroups.clear();
489 +      for (j = 0; j < molInfo.nCutoffGroups; j++){
490  
491 <      // send the arrays off to the forceField for init.
491 >        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
492 >        nMembers = currentCutoffGroup->getNMembers();
493  
494 +        myCutoffGroup = new CutoffGroup();
495 +        
496 +        for (int cg = 0; cg < nMembers; cg++) {
497 +
498 +          // molI is atom numbering inside this molecule
499 +          molI = currentCutoffGroup->getMember(cg);    
500 +
501 +          // tempI is atom numbering on local processor
502 +          tempI = molI + atomOffset;
503 +
504 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
505 +        }
506 +
507 +        molInfo.myCutoffGroups.push_back(myCutoffGroup);
508 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
509 +      
510 +
511 +
512 +      // After this is all set up, scan through the atoms to
513 +      // see if they can be added to the integrableObjects:
514 +
515 +      molInfo.myIntegrableObjects.clear();
516 +      
517 +
518 +      for (j = 0; j < molInfo.nAtoms; j++){
519 +
520 + #ifdef IS_MPI
521 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
522 + #else
523 +        slJ = j+atomOffset;
524 + #endif
525 +
526 +        // if they aren't on the skip list, then they can be integrated
527 +
528 +        if (skipList.find(slJ) == skipList.end()) {
529 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
530 +          info[k].integrableObjects.push_back(mySD);
531 +          molInfo.myIntegrableObjects.push_back(mySD);
532 +        }
533 +      }
534 +
535 +      // all rigid bodies are integrated:
536 +
537 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
538 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
539 +        info[k].integrableObjects.push_back(mySD);      
540 +        molInfo.myIntegrableObjects.push_back(mySD);
541 +      }
542 +    
543 +      
544 +      // send the arrays off to the forceField for init.
545 +      
546        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
547        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
548        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
549        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
550                                   theTorsions);
551  
418 –
552        info[k].molecules[i].initialize(molInfo);
553  
554  
# Line 423 | Line 556 | void SimSetup::makeMolecules(void){
556        delete[] theBonds;
557        delete[] theBends;
558        delete[] theTorsions;
559 <    }
559 >    }    
560    }
561  
562   #ifdef IS_MPI
# Line 433 | Line 566 | void SimSetup::makeMolecules(void){
566  
567    // clean up the forcefield
568  
569 <  the_ff->calcRcut();
569 >  if (!globals->haveRcut()){
570 >
571 >    the_ff->calcRcut();
572 >
573 >  } else {
574 >    
575 >    the_ff->setRcut( globals->getRcut() );
576 >  }
577 >
578    the_ff->cleanMe();
579   }
580  
# Line 634 | Line 775 | void SimSetup::gatherInfo(void){
775    }
776    else if (!strcasecmp(force_field, "EAM")){
777      ffCase = FF_EAM;
778 +  }
779 +  else if (!strcasecmp(force_field, "WATER")){
780 +    ffCase = FF_H2O;
781    }
782    else{
783      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
# Line 718 | Line 862 | void SimSetup::gatherInfo(void){
862    }
863  
864    //check whether sample time, status time, thermal time and reset time are divisble by dt
865 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
865 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
866      sprintf(painCave.errMsg,
867              "Sample time is not divisible by dt.\n"
868              "\tThis will result in samples that are not uniformly\n"
# Line 728 | Line 872 | void SimSetup::gatherInfo(void){
872      simError();    
873    }
874  
875 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
875 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
876      sprintf(painCave.errMsg,
877              "Status time is not divisible by dt.\n"
878              "\tThis will result in status reports that are not uniformly\n"
# Line 764 | Line 908 | void SimSetup::gatherInfo(void){
908      if (globals->haveSampleTime()){
909        info[i].sampleTime = globals->getSampleTime();
910        info[i].statusTime = info[i].sampleTime;
767 –      info[i].thermalTime = info[i].sampleTime;
911      }
912      else{
913        info[i].sampleTime = globals->getRunTime();
914        info[i].statusTime = info[i].sampleTime;
772 –      info[i].thermalTime = info[i].sampleTime;
915      }
916  
917      if (globals->haveStatusTime()){
# Line 778 | Line 920 | void SimSetup::gatherInfo(void){
920  
921      if (globals->haveThermalTime()){
922        info[i].thermalTime = globals->getThermalTime();
923 +    } else {
924 +      info[i].thermalTime = globals->getRunTime();
925      }
926  
927      info[i].resetIntegrator = 0;
# Line 837 | Line 981 | void SimSetup::gatherInfo(void){
981    for (int i = 0; i < nInfo; i++){
982      info[i].setSeed(seedValue);
983    }
984 <
984 >  
985   #ifdef IS_MPI
986 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
986 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
987    MPIcheckPoint();
988   #endif // is_mpi
989   }
# Line 848 | Line 992 | void SimSetup::finalInfoCheck(void){
992   void SimSetup::finalInfoCheck(void){
993    int index;
994    int usesDipoles;
995 +  int usesCharges;
996    int i;
997  
998    for (i = 0; i < nInfo; i++){
# Line 859 | Line 1004 | void SimSetup::finalInfoCheck(void){
1004        usesDipoles = (info[i].atoms[index])->hasDipole();
1005        index++;
1006      }
1007 <
1007 >    index = 0;
1008 >    usesCharges = 0;
1009 >    while ((index < info[i].n_atoms) && !usesCharges){
1010 >      usesCharges= (info[i].atoms[index])->hasCharge();
1011 >      index++;
1012 >    }
1013   #ifdef IS_MPI
1014      int myUse = usesDipoles;
1015      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1016   #endif //is_mpi
1017  
1018 <    double theEcr, theEst;
1018 >    double theRcut, theRsw;
1019  
1020      if (globals->getUseRF()){
1021        info[i].useReactionField = 1;
1022  
1023 <      if (!globals->haveECR()){
1023 >      if (!globals->haveRcut()){
1024          sprintf(painCave.errMsg,
1025 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1025 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1026                  "\tOOPSE will use a default value of 15.0 angstroms"
1027 <                "\tfor the electrostaticCutoffRadius.\n");
1027 >                "\tfor the cutoffRadius.\n");
1028          painCave.isFatal = 0;
1029          simError();
1030 <        theEcr = 15.0;
1030 >        theRcut = 15.0;
1031        }
1032        else{
1033 <        theEcr = globals->getECR();
1033 >        theRcut = globals->getRcut();
1034        }
1035  
1036 <      if (!globals->haveEST()){
1036 >      if (!globals->haveRsw()){
1037          sprintf(painCave.errMsg,
1038 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1038 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1039                  "\tOOPSE will use a default value of\n"
1040 <                "\t0.05 * electrostaticCutoffRadius\n"
891 <                "\tfor the electrostaticSkinThickness\n");
1040 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1041          painCave.isFatal = 0;
1042          simError();
1043 <        theEst = 0.05 * theEcr;
1043 >        theRsw = 0.95 * theRcut;
1044        }
1045        else{
1046 <        theEst = globals->getEST();
1046 >        theRsw = globals->getRsw();
1047        }
1048  
1049 <      info[i].setDefaultEcr(theEcr, theEst);
1049 >      info[i].setDefaultRcut(theRcut, theRsw);
1050  
1051        if (!globals->haveDielectric()){
1052          sprintf(painCave.errMsg,
# Line 910 | Line 1059 | void SimSetup::finalInfoCheck(void){
1059        info[i].dielectric = globals->getDielectric();
1060      }
1061      else{
1062 <      if (usesDipoles){
1063 <        if (!globals->haveECR()){
1062 >      if (usesDipoles || usesCharges){
1063 >
1064 >        if (!globals->haveRcut()){
1065            sprintf(painCave.errMsg,
1066 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1066 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1067                    "\tOOPSE will use a default value of 15.0 angstroms"
1068 <                  "\tfor the electrostaticCutoffRadius.\n");
1069 <          painCave.isFatal = 0;
1070 <          simError();
1071 <          theEcr = 15.0;
1072 <        }
1068 >                  "\tfor the cutoffRadius.\n");
1069 >          painCave.isFatal = 0;
1070 >          simError();
1071 >          theRcut = 15.0;
1072 >      }
1073          else{
1074 <          theEcr = globals->getECR();
1074 >          theRcut = globals->getRcut();
1075          }
1076 <        
1077 <        if (!globals->haveEST()){
1078 <          sprintf(painCave.errMsg,
1079 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1076 >        
1077 >        if (!globals->haveRsw()){
1078 >          sprintf(painCave.errMsg,
1079 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1080                    "\tOOPSE will use a default value of\n"
1081 <                  "\t0.05 * electrostaticCutoffRadius\n"
932 <                  "\tfor the electrostaticSkinThickness\n");
1081 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1082            painCave.isFatal = 0;
1083            simError();
1084 <          theEst = 0.05 * theEcr;
1084 >          theRsw = 0.95 * theRcut;
1085          }
1086          else{
1087 <          theEst = globals->getEST();
1087 >          theRsw = globals->getRsw();
1088          }
1089 +        
1090 +        info[i].setDefaultRcut(theRcut, theRsw);
1091          
941 –        info[i].setDefaultEcr(theEcr, theEst);
1092        }
1093      }
1094    }
# Line 1136 | Line 1286 | void SimSetup::createFF(void){
1286        the_ff = new EAM_FF();
1287        break;
1288  
1289 +    case FF_H2O:
1290 +      the_ff = new WATER();
1291 +      break;
1292 +
1293      default:
1294        sprintf(painCave.errMsg,
1295                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1156 | Line 1310 | void SimSetup::compList(void){
1310    LinkedMolStamp* headStamp = new LinkedMolStamp();
1311    LinkedMolStamp* currentStamp = NULL;
1312    comp_stamps = new MoleculeStamp * [n_components];
1313 +  bool haveCutoffGroups;
1314  
1315 +  haveCutoffGroups = false;
1316 +  
1317    // make an array of molecule stamps that match the components used.
1318    // also extract the used stamps out into a separate linked list
1319  
# Line 1191 | Line 1348 | void SimSetup::compList(void){
1348        headStamp->add(currentStamp);
1349        comp_stamps[i] = headStamp->match(id);
1350      }
1351 +
1352 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1353 +      haveCutoffGroups = true;    
1354    }
1355 +    
1356 +  for (i = 0; i < nInfo; i++)
1357 +    info[i].haveCutoffGroups = haveCutoffGroups;
1358  
1359   #ifdef IS_MPI
1360    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1208 | Line 1371 | void SimSetup::calcSysValues(void){
1371    tot_bonds = 0;
1372    tot_bends = 0;
1373    tot_torsions = 0;
1374 +  tot_rigid = 0;
1375    for (i = 0; i < n_components; i++){
1376      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1377      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1378      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1379      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1380 +    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1381    }
1382 <
1382 >  
1383    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1384    molMembershipArray = new int[tot_atoms];
1385  
# Line 1236 | Line 1401 | void SimSetup::mpiMolDivide(void){
1401    int i, j, k;
1402    int localMol, allMol;
1403    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1404 +  int local_rigid;
1405 +  vector<int> globalMolIndex;
1406  
1407    mpiSim = new mpiSimulation(info);
1408  
1409 <  globalIndex = mpiSim->divideLabor();
1409 >  mpiSim->divideLabor();
1410 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1411 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1412  
1413    // set up the local variables
1414  
# Line 1252 | Line 1421 | void SimSetup::mpiMolDivide(void){
1421    local_bonds = 0;
1422    local_bends = 0;
1423    local_torsions = 0;
1424 <  globalAtomIndex = 0;
1424 >  local_rigid = 0;
1425 >  globalAtomCounter = 0;
1426  
1257 –
1427    for (i = 0; i < n_components; i++){
1428      for (j = 0; j < components_nmol[i]; j++){
1429        if (mol2proc[allMol] == worldRank){
# Line 1262 | Line 1431 | void SimSetup::mpiMolDivide(void){
1431          local_bonds += comp_stamps[i]->getNBonds();
1432          local_bends += comp_stamps[i]->getNBends();
1433          local_torsions += comp_stamps[i]->getNTorsions();
1434 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1435          localMol++;
1436        }      
1437        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1438 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1439 <        globalAtomIndex++;
1438 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1439 >        globalAtomCounter++;
1440        }
1441  
1442        allMol++;
# Line 1275 | Line 1445 | void SimSetup::mpiMolDivide(void){
1445    local_SRI = local_bonds + local_bends + local_torsions;
1446  
1447    info[0].n_atoms = mpiSim->getMyNlocal();  
1448 +  
1449  
1450    if (local_atoms != info[0].n_atoms){
1451      sprintf(painCave.errMsg,
# Line 1307 | Line 1478 | void SimSetup::makeSysArrays(void){
1478  
1479    Atom** the_atoms;
1480    Molecule* the_molecules;
1310 –  Exclude** the_excludes;
1481  
1312 –
1482    for (l = 0; l < nInfo; l++){
1483      // create the atom and short range interaction arrays
1484  
# Line 1335 | Line 1504 | void SimSetup::makeSysArrays(void){
1504   #else // is_mpi
1505  
1506      molIndex = 0;
1507 <    globalAtomIndex = 0;
1507 >    globalAtomCounter = 0;
1508      for (i = 0; i < n_components; i++){
1509        for (j = 0; j < components_nmol[i]; j++){
1510          the_molecules[molIndex].setStampID(i);
1511          the_molecules[molIndex].setMyIndex(molIndex);
1512          the_molecules[molIndex].setGlobalIndex(molIndex);
1513          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1514 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1515 <          globalAtomIndex++;
1514 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1515 >          globalAtomCounter++;
1516          }
1517          molIndex++;
1518        }
# Line 1352 | Line 1521 | void SimSetup::makeSysArrays(void){
1521  
1522   #endif // is_mpi
1523  
1524 <
1525 <    if (info[l].n_SRI){
1526 <      Exclude::createArray(info[l].n_SRI);
1358 <      the_excludes = new Exclude * [info[l].n_SRI];
1359 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1360 <        the_excludes[ex] = new Exclude(ex);
1361 <      }
1362 <      info[l].globalExcludes = new int;
1363 <      info[l].n_exclude = info[l].n_SRI;
1364 <    }
1365 <    else{
1366 <      Exclude::createArray(1);
1367 <      the_excludes = new Exclude * ;
1368 <      the_excludes[0] = new Exclude(0);
1369 <      the_excludes[0]->setPair(0, 0);
1370 <      info[l].globalExcludes = new int;
1371 <      info[l].globalExcludes[0] = 0;
1372 <      info[l].n_exclude = 0;
1373 <    }
1374 <
1524 >    info[l].globalExcludes = new int;
1525 >    info[l].globalExcludes[0] = 0;
1526 >    
1527      // set the arrays into the SimInfo object
1528  
1529      info[l].atoms = the_atoms;
1530      info[l].molecules = the_molecules;
1531      info[l].nGlobalExcludes = 0;
1532 <    info[l].excludes = the_excludes;
1381 <
1532 >    
1533      the_ff->setSimInfo(info);
1534    }
1535   }
# Line 1650 | Line 1801 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1801    }
1802  
1803    theInfo.addProperty(zconsForcePolicy);
1804 +
1805 +  //set zcons gap
1806 +  DoubleData* zconsGap = new DoubleData();
1807 +  zconsGap->setID(ZCONSGAP_ID);
1808 +
1809 +  if (globals->haveZConsGap()){
1810 +    zconsGap->setData(globals->getZconsGap());
1811 +    theInfo.addProperty(zconsGap);  
1812 +  }
1813 +
1814 +  //set zcons fixtime
1815 +  DoubleData* zconsFixtime = new DoubleData();
1816 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1817 +
1818 +  if (globals->haveZConsFixTime()){
1819 +    zconsFixtime->setData(globals->getZconsFixtime());
1820 +    theInfo.addProperty(zconsFixtime);  
1821 +  }
1822 +
1823 +  //set zconsUsingSMD
1824 +  IntData* zconsUsingSMD = new IntData();
1825 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1826 +
1827 +  if (globals->haveZConsUsingSMD()){
1828 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1829 +    theInfo.addProperty(zconsUsingSMD);  
1830 +  }
1831  
1832    //Determine the name of ouput file and add it into SimInfo's property list
1833    //Be careful, do not use inFileName, since it is a pointer which
# Line 1680 | Line 1858 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1858      tempParaItem.zPos = zconStamp[i]->getZpos();
1859      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1860      tempParaItem.kRatio = zconStamp[i]->getKratio();
1861 <
1861 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1862 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1863      zconsParaData->addItem(tempParaItem);
1864    }
1865  
# Line 1698 | Line 1877 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1877    //push data into siminfo, therefore, we can retrieve later
1878    theInfo.addProperty(zconsParaData);
1879   }
1880 +
1881 + void SimSetup::makeMinimizer(){
1882 +
1883 +  OOPSEMinimizer* myOOPSEMinimizer;
1884 +  MinimizerParameterSet* param;
1885 +  char minimizerName[100];
1886 +  
1887 +  for (int i = 0; i < nInfo; i++){
1888 +    
1889 +    //prepare parameter set for minimizer
1890 +    param = new MinimizerParameterSet();
1891 +    param->setDefaultParameter();
1892 +
1893 +    if (globals->haveMinimizer()){
1894 +      param->setFTol(globals->getMinFTol());
1895 +    }
1896 +
1897 +    if (globals->haveMinGTol()){
1898 +      param->setGTol(globals->getMinGTol());
1899 +    }
1900 +
1901 +    if (globals->haveMinMaxIter()){
1902 +      param->setMaxIteration(globals->getMinMaxIter());
1903 +    }
1904 +
1905 +    if (globals->haveMinWriteFrq()){
1906 +      param->setMaxIteration(globals->getMinMaxIter());
1907 +    }
1908 +
1909 +    if (globals->haveMinWriteFrq()){
1910 +      param->setWriteFrq(globals->getMinWriteFrq());
1911 +    }
1912 +    
1913 +    if (globals->haveMinStepSize()){
1914 +      param->setStepSize(globals->getMinStepSize());
1915 +    }
1916 +
1917 +    if (globals->haveMinLSMaxIter()){
1918 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
1919 +    }    
1920 +
1921 +    if (globals->haveMinLSTol()){
1922 +      param->setLineSearchTol(globals->getMinLSTol());
1923 +    }    
1924 +
1925 +    strcpy(minimizerName, globals->getMinimizer());
1926 +
1927 +    if (!strcasecmp(minimizerName, "CG")){
1928 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1929 +    }
1930 +    else if (!strcasecmp(minimizerName, "SD")){
1931 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1932 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1933 +    }
1934 +    else{
1935 +          sprintf(painCave.errMsg,
1936 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1937 +          painCave.isFatal = 0;
1938 +          simError();
1939 +
1940 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1941 +    }
1942 +     info[i].the_integrator = myOOPSEMinimizer;
1943 +
1944 +     //store the minimizer into simInfo
1945 +     info[i].the_minimizer = myOOPSEMinimizer;
1946 +     info[i].has_minimizer = true;
1947 +  }
1948 +
1949 + }

Diff Legend

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