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 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC vs.
Revision 1108 by tim, Wed Apr 14 15:37: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  
37 + /**
38 + * Check whether dividend is divisble by divisor or not
39 + */
40 + bool isDivisible(double dividend, double divisor){
41 +  double tolerance = 0.000001;
42 +  double quotient;
43 +  double diff;
44 +  int intQuotient;
45 +  
46 +  quotient = dividend / divisor;
47 +
48 +  if (quotient < 0)
49 +    quotient = -quotient;
50 +
51 +  intQuotient = int (quotient + tolerance);
52 +
53 +  diff = fabs(fabs(dividend) - intQuotient  * fabs(divisor));
54 +
55 +  if (diff <= tolerance)
56 +    return true;
57 +  else
58 +    return false;  
59 + }
60 +
61   SimSetup::SimSetup(){
62    
63    initSuspend = false;
# Line 120 | Line 148 | void SimSetup::createSim(void){
148  
149    makeOutNames();
150  
151 <  // make the integrator
152 <
153 <  makeIntegrator();
154 <
151 >  if (globals->haveMinimizer())
152 >    // make minimizer
153 >    makeMinimizer();
154 >  else
155 >    // make the integrator
156 >    makeIntegrator();
157 >  
158   #ifdef IS_MPI
159    mpiSim->mpiRefresh();
160   #endif
# Line 135 | 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, slJ;
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 +  char* molName;
194 +  char rbName[100];
195 +
196    //init the forceField paramters
197  
198    the_ff->readParams();
199  
158 –
200    // init the atoms
201  
202 <  double ux, uy, uz, u, uSqr;
202 >  int nMembers, nNew, rb1, rb2;
203  
204    for (k = 0; k < nInfo; k++){
205      the_ff->setSimInfo(&(info[k]));
206  
207      atomOffset = 0;
208 <    excludeOffset = 0;
208 >
209      for (i = 0; i < info[k].n_mol; i++){
210        stampID = info[k].molecules[i].getStampID();
211 +      molName = comp_stamps[stampID]->getID();
212  
213        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
214        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
215        molInfo.nBends = comp_stamps[stampID]->getNBends();
216        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
217 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
218 <
217 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
218 >      
219        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
178 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
179 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
180 –      molInfo.myBends = new Bend * [molInfo.nBends];
181 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
220  
221 +      if (molInfo.nBonds > 0)
222 +        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
223 +      else
224 +        molInfo.myBonds = NULL;
225 +
226 +      if (molInfo.nBends > 0)
227 +        molInfo.myBends = new (Bend *) [molInfo.nBends];
228 +      else
229 +        molInfo.myBends = NULL;
230 +
231 +      if (molInfo.nTorsions > 0)
232 +        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
233 +      else
234 +        molInfo.myTorsions = NULL;
235 +
236        theBonds = new bond_pair[molInfo.nBonds];
237        theBends = new bend_set[molInfo.nBends];
238        theTorsions = new torsion_set[molInfo.nTorsions];
239 <
239 >      
240        // make the Atoms
241  
242        for (j = 0; j < molInfo.nAtoms; j++){
243          currentAtom = comp_stamps[stampID]->getAtom(j);
244 +
245          if (currentAtom->haveOrientation()){
246            dAtom = new DirectionalAtom((j + atomOffset),
247                                        info[k].getConfiguration());
248            info[k].n_oriented++;
249            molInfo.myAtoms[j] = dAtom;
250  
251 <          ux = currentAtom->getOrntX();
252 <          uy = currentAtom->getOrntY();
253 <          uz = currentAtom->getOrntZ();
251 >          // Directional Atoms have standard unit vectors which are oriented
252 >          // in space using the three Euler angles.  We assume the standard
253 >          // unit vector was originally along the z axis below.
254  
255 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
255 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
256 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
257 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
258  
259 <          u = sqrt(uSqr);
260 <          ux = ux / u;
205 <          uy = uy / u;
206 <          uz = uz / u;
207 <
208 <          dAtom->setSUx(ux);
209 <          dAtom->setSUy(uy);
210 <          dAtom->setSUz(uz);
259 >          dAtom->setUnitFrameFromEuler(phi, theta, psi);
260 >            
261          }
262          else{
263 <          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
264 <                                               info[k].getConfiguration());
263 >
264 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
265          }
266 +
267          molInfo.myAtoms[j]->setType(currentAtom->getType());
268  
269   #ifdef IS_MPI
270  
271 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
271 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
272  
273   #endif // is_mpi
274        }
# Line 228 | Line 279 | void SimSetup::makeMolecules(void){
279          theBonds[j].a = currentBond->getA() + atomOffset;
280          theBonds[j].b = currentBond->getB() + atomOffset;
281  
282 <        exI = theBonds[j].a;
283 <        exJ = theBonds[j].b;
282 >        tempI = theBonds[j].a;
283 >        tempJ = theBonds[j].b;
284  
234 –        // exclude_I must always be the smaller of the pair
235 –        if (exI > exJ){
236 –          tempEx = exI;
237 –          exI = exJ;
238 –          exJ = tempEx;
239 –        }
285   #ifdef IS_MPI
286 <        tempEx = exI;
287 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
288 <        tempEx = exJ;
289 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
286 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
287 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
288 > #else
289 >        exI = tempI + 1;
290 >        exJ = tempJ + 1;
291 > #endif
292  
293 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
247 < #else  // isn't MPI
248 <
249 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
250 < #endif  //is_mpi
293 >        info[k].excludes->addPair(exI, exJ);
294        }
252 –      excludeOffset += molInfo.nBonds;
295  
296        //make the bends
297        for (j = 0; j < molInfo.nBends; j++){
# Line 299 | Line 341 | void SimSetup::makeMolecules(void){
341            }
342          }
343  
344 <        if (!theBends[j].isGhost){
345 <          exI = theBends[j].a;
346 <          exJ = theBends[j].c;
347 <        }
348 <        else{
307 <          exI = theBends[j].a;
308 <          exJ = theBends[j].b;
309 <        }
310 <
311 <        // exclude_I must always be the smaller of the pair
312 <        if (exI > exJ){
313 <          tempEx = exI;
314 <          exI = exJ;
315 <          exJ = tempEx;
316 <        }
344 >        if (theBends[j].isGhost) {
345 >          
346 >          tempI = theBends[j].a;
347 >          tempJ = theBends[j].b;
348 >          
349   #ifdef IS_MPI
350 <        tempEx = exI;
351 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
352 <        tempEx = exJ;
353 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
350 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
351 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
352 > #else
353 >          exI = tempI + 1;
354 >          exJ = tempJ + 1;
355 > #endif          
356 >          info[k].excludes->addPair(exI, exJ);
357  
358 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
359 < #else  // isn't MPI
360 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
361 < #endif  //is_mpi
358 >        } else {
359 >
360 >          tempI = theBends[j].a;
361 >          tempJ = theBends[j].b;
362 >          tempK = theBends[j].c;
363 >          
364 > #ifdef IS_MPI
365 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
366 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
367 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
368 > #else
369 >          exI = tempI + 1;
370 >          exJ = tempJ + 1;
371 >          exK = tempK + 1;
372 > #endif
373 >          
374 >          info[k].excludes->addPair(exI, exK);
375 >          info[k].excludes->addPair(exI, exJ);
376 >          info[k].excludes->addPair(exJ, exK);
377 >        }
378        }
328 –      excludeOffset += molInfo.nBends;
379  
380        for (j = 0; j < molInfo.nTorsions; j++){
381          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 334 | Line 384 | void SimSetup::makeMolecules(void){
384          theTorsions[j].c = currentTorsion->getC() + atomOffset;
385          theTorsions[j].d = currentTorsion->getD() + atomOffset;
386  
387 <        exI = theTorsions[j].a;
388 <        exJ = theTorsions[j].d;
387 >        tempI = theTorsions[j].a;      
388 >        tempJ = theTorsions[j].b;
389 >        tempK = theTorsions[j].c;
390 >        tempL = theTorsions[j].d;
391  
340 –        // exclude_I must always be the smaller of the pair
341 –        if (exI > exJ){
342 –          tempEx = exI;
343 –          exI = exJ;
344 –          exJ = tempEx;
345 –        }
392   #ifdef IS_MPI
393 <        tempEx = exI;
394 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
395 <        tempEx = exJ;
396 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
393 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
394 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
395 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
396 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
397 > #else
398 >        exI = tempI + 1;
399 >        exJ = tempJ + 1;
400 >        exK = tempK + 1;
401 >        exL = tempL + 1;
402 > #endif
403  
404 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
405 < #else  // isn't MPI
406 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
407 < #endif  //is_mpi
404 >        info[k].excludes->addPair(exI, exJ);
405 >        info[k].excludes->addPair(exI, exK);
406 >        info[k].excludes->addPair(exI, exL);        
407 >        info[k].excludes->addPair(exJ, exK);
408 >        info[k].excludes->addPair(exJ, exL);
409 >        info[k].excludes->addPair(exK, exL);
410        }
357 –      excludeOffset += molInfo.nTorsions;
411  
412 +      for (j = 0; j < molInfo.nRigidBodies; j++){
413  
414 <      // send the arrays off to the forceField for init.
414 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
415 >        nMembers = currentRigidBody->getNMembers();
416  
417 +        // Create the Rigid Body:
418 +
419 +        myRB = new RigidBody();
420 +
421 +        sprintf(rbName,"%s_RB_%s", molName, j);
422 +        myRB->setType(rbName);
423 +        
424 +        for (rb1 = 0; rb1 < nMembers; rb1++) {
425 +
426 +          // molI is atom numbering inside this molecule
427 +          molI = currentRigidBody->getMember(rb1);    
428 +
429 +          // tempI is atom numbering on local processor
430 +          tempI = molI + atomOffset;
431 +
432 +          // currentAtom is the AtomStamp (which we need for
433 +          // rigid body reference positions)
434 +          currentAtom = comp_stamps[stampID]->getAtom(molI);
435 +
436 +          // When we add to the rigid body, add the atom itself and
437 +          // the stamp info:
438 +
439 +          myRB->addAtom(info[k].atoms[tempI], currentAtom);
440 +          
441 +          // Add this atom to the Skip List for the integrators
442 + #ifdef IS_MPI
443 +          slI = info[k].atoms[tempI]->getGlobalIndex();
444 + #else
445 +          slI = tempI;
446 + #endif
447 +          skipList.insert(slI);
448 +          
449 +        }
450 +        
451 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
452 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
453 +            
454 +            tempI = currentRigidBody->getMember(rb1);
455 +            tempJ = currentRigidBody->getMember(rb2);
456 +            
457 +            // Some explanation is required here.
458 +            // Fortran indexing starts at 1, while c indexing starts at 0
459 +            // Also, in parallel computations, the GlobalIndex is
460 +            // used for the exclude list:
461 +            
462 + #ifdef IS_MPI
463 +            exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
464 +            exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
465 + #else
466 +            exI = tempI + 1;
467 +            exJ = tempJ + 1;
468 + #endif
469 +            
470 +            info[k].excludes->addPair(exI, exJ);
471 +            
472 +          }
473 +        }
474 +
475 +        molInfo.myRigidBodies.push_back(myRB);
476 +        info[k].rigidBodies.push_back(myRB);
477 +      }
478 +      
479 +
480 +      // After this is all set up, scan through the atoms to
481 +      // see if they can be added to the integrableObjects:
482 +
483 +      for (j = 0; j < molInfo.nAtoms; j++){
484 +
485 + #ifdef IS_MPI
486 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
487 + #else
488 +        slJ = j+atomOffset;
489 + #endif
490 +
491 +        // if they aren't on the skip list, then they can be integrated
492 +
493 +        if (skipList.find(slJ) == skipList.end()) {
494 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
495 +          info[k].integrableObjects.push_back(mySD);
496 +          molInfo.myIntegrableObjects.push_back(mySD);
497 +        }
498 +      }
499 +
500 +      // all rigid bodies are integrated:
501 +
502 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
503 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
504 +        info[k].integrableObjects.push_back(mySD);      
505 +        molInfo.myIntegrableObjects.push_back(mySD);
506 +      }
507 +    
508 +      
509 +      // send the arrays off to the forceField for init.
510 +      
511        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
512        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
513        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
514        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
515                                   theTorsions);
516  
368 –
517        info[k].molecules[i].initialize(molInfo);
518  
519  
# Line 373 | Line 521 | void SimSetup::makeMolecules(void){
521        delete[] theBonds;
522        delete[] theBends;
523        delete[] theTorsions;
524 <    }
524 >    }    
525    }
526  
527   #ifdef IS_MPI
# Line 383 | Line 531 | void SimSetup::makeMolecules(void){
531  
532    // clean up the forcefield
533  
534 <  the_ff->calcRcut();
534 >  if (!globals->haveLJrcut()){
535 >
536 >    the_ff->calcRcut();
537 >
538 >  } else {
539 >    
540 >    the_ff->setRcut( globals->getLJrcut() );
541 >  }
542 >
543    the_ff->cleanMe();
544   }
545  
# Line 584 | Line 740 | void SimSetup::gatherInfo(void){
740    }
741    else if (!strcasecmp(force_field, "EAM")){
742      ffCase = FF_EAM;
743 +  }
744 +  else if (!strcasecmp(force_field, "WATER")){
745 +    ffCase = FF_H2O;
746    }
747    else{
748      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
# Line 613 | Line 772 | void SimSetup::gatherInfo(void){
772    }
773    else{
774      sprintf(painCave.errMsg,
775 <            "SimSetup Warning. Unrecognized Ensemble -> %s, "
776 <            "reverting to NVE for this simulation.\n",
775 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
776 >            "\treverting to NVE for this simulation.\n",
777              ensemble);
778           painCave.isFatal = 0;
779           simError();
# Line 646 | Line 805 | void SimSetup::gatherInfo(void){
805        if (!the_components[i]->haveNMol()){
806          // we have a problem
807          sprintf(painCave.errMsg,
808 <                "SimSetup Error. No global NMol or component NMol"
809 <                " given. Cannot calculate the number of atoms.\n");
808 >                "SimSetup Error. No global NMol or component NMol given.\n"
809 >                "\tCannot calculate the number of atoms.\n");
810          painCave.isFatal = 1;
811          simError();
812        }
# Line 665 | Line 824 | void SimSetup::gatherInfo(void){
824              " Please give nMol in the components.\n");
825      painCave.isFatal = 1;
826      simError();
827 +  }
828 +
829 +  //check whether sample time, status time, thermal time and reset time are divisble by dt
830 +  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
831 +    sprintf(painCave.errMsg,
832 +            "Sample time is not divisible by dt.\n"
833 +            "\tThis will result in samples that are not uniformly\n"
834 +            "\tdistributed in time.  If this is a problem, change\n"
835 +            "\tyour sampleTime variable.\n");
836 +    painCave.isFatal = 0;
837 +    simError();    
838 +  }
839 +
840 +  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
841 +    sprintf(painCave.errMsg,
842 +            "Status time is not divisible by dt.\n"
843 +            "\tThis will result in status reports that are not uniformly\n"
844 +            "\tdistributed in time.  If this is a problem, change \n"
845 +            "\tyour statusTime variable.\n");
846 +    painCave.isFatal = 0;
847 +    simError();    
848    }
849  
850 +  if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
851 +    sprintf(painCave.errMsg,
852 +            "Thermal time is not divisible by dt.\n"
853 +            "\tThis will result in thermalizations that are not uniformly\n"
854 +            "\tdistributed in time.  If this is a problem, change \n"
855 +            "\tyour thermalTime variable.\n");
856 +    painCave.isFatal = 0;
857 +    simError();    
858 +  }  
859 +
860 +  if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
861 +    sprintf(painCave.errMsg,
862 +            "Reset time is not divisible by dt.\n"
863 +            "\tThis will result in integrator resets that are not uniformly\n"
864 +            "\tdistributed in time.  If this is a problem, change\n"
865 +            "\tyour resetTime variable.\n");
866 +    painCave.isFatal = 0;
867 +    simError();    
868 +  }
869 +
870    // set the status, sample, and thermal kick times
871  
872    for (i = 0; i < nInfo; i++){
# Line 696 | Line 896 | void SimSetup::gatherInfo(void){
896      }
897  
898      // check for the temperature set flag
899 <
899 >    
900      if (globals->haveTempSet())
901        info[i].setTemp = globals->getTempSet();
902  
903 <    // get some of the tricky things that may still be in the globals
903 >    // check for the extended State init
904  
905 <    double boxVector[3];
906 <    if (globals->haveBox()){
907 <      boxVector[0] = globals->getBox();
708 <      boxVector[1] = globals->getBox();
709 <      boxVector[2] = globals->getBox();
710 <
711 <      info[i].setBox(boxVector);
712 <    }
713 <    else if (globals->haveDensity()){
714 <      double vol;
715 <      vol = (double) tot_nmol / globals->getDensity();
716 <      boxVector[0] = pow(vol, (1.0 / 3.0));
717 <      boxVector[1] = boxVector[0];
718 <      boxVector[2] = boxVector[0];
719 <
720 <      info[i].setBox(boxVector);
721 <    }
722 <    else{
723 <      if (!globals->haveBoxX()){
724 <        sprintf(painCave.errMsg,
725 <                "SimSetup error, no periodic BoxX size given.\n");
726 <        painCave.isFatal = 1;
727 <        simError();
728 <      }
729 <      boxVector[0] = globals->getBoxX();
730 <
731 <      if (!globals->haveBoxY()){
732 <        sprintf(painCave.errMsg,
733 <                "SimSetup error, no periodic BoxY size given.\n");
734 <        painCave.isFatal = 1;
735 <        simError();
736 <      }
737 <      boxVector[1] = globals->getBoxY();
738 <
739 <      if (!globals->haveBoxZ()){
740 <        sprintf(painCave.errMsg,
741 <                "SimSetup error, no periodic BoxZ size given.\n");
742 <        painCave.isFatal = 1;
743 <        simError();
744 <      }
745 <      boxVector[2] = globals->getBoxZ();
746 <
747 <      info[i].setBox(boxVector);
748 <    }
905 >    info[i].useInitXSstate = globals->getUseInitXSstate();
906 >    info[i].orthoTolerance = globals->getOrthoBoxTolerance();
907 >    
908    }
909 <
909 >  
910    //setup seed for random number generator
911    int seedValue;
912  
# Line 787 | Line 946 | void SimSetup::gatherInfo(void){
946    for (int i = 0; i < nInfo; i++){
947      info[i].setSeed(seedValue);
948    }
949 <
949 >  
950   #ifdef IS_MPI
951 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
951 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
952    MPIcheckPoint();
953   #endif // is_mpi
954   }
# Line 822 | Line 981 | void SimSetup::finalInfoCheck(void){
981  
982        if (!globals->haveECR()){
983          sprintf(painCave.errMsg,
984 <                "SimSetup Warning: using default value of 1/2 the smallest "
985 <                "box length for the electrostaticCutoffRadius.\n"
986 <                "I hope you have a very fast processor!\n");
984 >                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
985 >                "\tOOPSE will use a default value of 15.0 angstroms"
986 >                "\tfor the electrostaticCutoffRadius.\n");
987          painCave.isFatal = 0;
988          simError();
989 <        double smallest;
831 <        smallest = info[i].boxL[0];
832 <        if (info[i].boxL[1] <= smallest)
833 <          smallest = info[i].boxL[1];
834 <        if (info[i].boxL[2] <= smallest)
835 <          smallest = info[i].boxL[2];
836 <        theEcr = 0.5 * smallest;
989 >        theEcr = 15.0;
990        }
991        else{
992          theEcr = globals->getECR();
# Line 841 | Line 994 | void SimSetup::finalInfoCheck(void){
994  
995        if (!globals->haveEST()){
996          sprintf(painCave.errMsg,
997 <                "SimSetup Warning: using default value of 0.05 * the "
998 <                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
997 >                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
998 >                "\tOOPSE will use a default value of\n"
999 >                "\t0.05 * electrostaticCutoffRadius\n"
1000 >                "\tfor the electrostaticSkinThickness\n");
1001          painCave.isFatal = 0;
1002          simError();
1003          theEst = 0.05 * theEcr;
# Line 851 | Line 1006 | void SimSetup::finalInfoCheck(void){
1006          theEst = globals->getEST();
1007        }
1008  
1009 <      info[i].setEcr(theEcr, theEst);
1009 >      info[i].setDefaultEcr(theEcr, theEst);
1010  
1011        if (!globals->haveDielectric()){
1012          sprintf(painCave.errMsg,
1013 <                "SimSetup Error: You are trying to use Reaction Field without"
1014 <                "setting a dielectric constant!\n");
1013 >                "SimSetup Error: No Dielectric constant was set.\n"
1014 >                "\tYou are trying to use Reaction Field without"
1015 >                "\tsetting a dielectric constant!\n");
1016          painCave.isFatal = 1;
1017          simError();
1018        }
# Line 866 | Line 1022 | void SimSetup::finalInfoCheck(void){
1022        if (usesDipoles){
1023          if (!globals->haveECR()){
1024            sprintf(painCave.errMsg,
1025 <                  "SimSetup Warning: using default value of 1/2 the smallest "
1026 <                  "box length for the electrostaticCutoffRadius.\n"
1027 <                  "I hope you have a very fast processor!\n");
1028 <          painCave.isFatal = 0;
1029 <          simError();
1030 <          double smallest;
875 <          smallest = info[i].boxL[0];
876 <          if (info[i].boxL[1] <= smallest)
877 <            smallest = info[i].boxL[1];
878 <          if (info[i].boxL[2] <= smallest)
879 <            smallest = info[i].boxL[2];
880 <          theEcr = 0.5 * smallest;
1025 >                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1026 >                  "\tOOPSE will use a default value of 15.0 angstroms"
1027 >                  "\tfor the electrostaticCutoffRadius.\n");
1028 >          painCave.isFatal = 0;
1029 >          simError();
1030 >          theEcr = 15.0;
1031          }
1032          else{
1033            theEcr = globals->getECR();
1034          }
1035 <
1035 >        
1036          if (!globals->haveEST()){
1037            sprintf(painCave.errMsg,
1038 <                  "SimSetup Warning: using default value of 0.05 * the "
1039 <                  "electrostaticCutoffRadius for the "
1040 <                  "electrostaticSkinThickness\n");
1038 >                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1039 >                  "\tOOPSE will use a default value of\n"
1040 >                  "\t0.05 * electrostaticCutoffRadius\n"
1041 >                  "\tfor the electrostaticSkinThickness\n");
1042            painCave.isFatal = 0;
1043            simError();
1044            theEst = 0.05 * theEcr;
# Line 895 | Line 1046 | void SimSetup::finalInfoCheck(void){
1046          else{
1047            theEst = globals->getEST();
1048          }
1049 <
1050 <        info[i].setEcr(theEcr, theEst);
1049 >        
1050 >        info[i].setDefaultEcr(theEcr, theEst);
1051        }
1052      }
1053    }
903 –
1054   #ifdef IS_MPI
1055    strcpy(checkPointMsg, "post processing checks out");
1056    MPIcheckPoint();
1057   #endif // is_mpi
1058   }
1059 <
1059 >  
1060   void SimSetup::initSystemCoords(void){
1061    int i;
1062  
# Line 934 | Line 1084 | void SimSetup::initSystemCoords(void){
1084      delete fileInit;
1085    }
1086    else{
1087 < #ifdef IS_MPI
938 <
1087 >    
1088      // no init from bass
1089 <
1089 >    
1090      sprintf(painCave.errMsg,
1091 <            "Cannot intialize a parallel simulation without an initial configuration file.\n");
1091 >            "Cannot intialize a simulation without an initial configuration file.\n");
1092      painCave.isFatal = 1;;
1093      simError();
1094 <
946 < #else
947 <
948 <    initFromBass();
949 <
950 <
951 < #endif
1094 >    
1095    }
1096  
1097   #ifdef IS_MPI
# Line 1100 | Line 1243 | void SimSetup::createFF(void){
1243  
1244      case FF_EAM:
1245        the_ff = new EAM_FF();
1246 +      break;
1247 +
1248 +    case FF_H2O:
1249 +      the_ff = new WATER();
1250        break;
1251  
1252      default:
# Line 1174 | Line 1321 | void SimSetup::calcSysValues(void){
1321    tot_bonds = 0;
1322    tot_bends = 0;
1323    tot_torsions = 0;
1324 +  tot_rigid = 0;
1325    for (i = 0; i < n_components; i++){
1326      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1327      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1328      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1329      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1330 +    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1331    }
1332 <
1332 >  
1333    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1334    molMembershipArray = new int[tot_atoms];
1335  
# Line 1202 | Line 1351 | void SimSetup::mpiMolDivide(void){
1351    int i, j, k;
1352    int localMol, allMol;
1353    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1354 +  int local_rigid;
1355 +  vector<int> globalAtomIndex;
1356 +  vector<int> globalMolIndex;
1357  
1358    mpiSim = new mpiSimulation(info);
1359  
1360 <  globalIndex = mpiSim->divideLabor();
1360 >  mpiSim->divideLabor();
1361 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1362 >  globalMolIndex = mpiSim->getGlobalMolIndex();
1363  
1364    // set up the local variables
1365  
# Line 1218 | Line 1372 | void SimSetup::mpiMolDivide(void){
1372    local_bonds = 0;
1373    local_bends = 0;
1374    local_torsions = 0;
1375 <  globalAtomIndex = 0;
1376 <
1375 >  local_rigid = 0;
1376 >  globalAtomCounter = 0;
1377  
1378    for (i = 0; i < n_components; i++){
1379      for (j = 0; j < components_nmol[i]; j++){
# Line 1228 | Line 1382 | void SimSetup::mpiMolDivide(void){
1382          local_bonds += comp_stamps[i]->getNBonds();
1383          local_bends += comp_stamps[i]->getNBends();
1384          local_torsions += comp_stamps[i]->getNTorsions();
1385 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1386          localMol++;
1387        }      
1388        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1389 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1390 <        globalAtomIndex++;
1389 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1390 >        globalAtomCounter++;
1391        }
1392  
1393        allMol++;
# Line 1241 | Line 1396 | void SimSetup::mpiMolDivide(void){
1396    local_SRI = local_bonds + local_bends + local_torsions;
1397  
1398    info[0].n_atoms = mpiSim->getMyNlocal();  
1399 +  
1400  
1401    if (local_atoms != info[0].n_atoms){
1402      sprintf(painCave.errMsg,
1403 <            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1404 <            " localAtom (%d) are not equal.\n",
1403 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1404 >            "\tlocalAtom (%d) are not equal.\n",
1405              info[0].n_atoms, local_atoms);
1406      painCave.isFatal = 1;
1407      simError();
# Line 1273 | Line 1429 | void SimSetup::makeSysArrays(void){
1429  
1430    Atom** the_atoms;
1431    Molecule* the_molecules;
1276 –  Exclude** the_excludes;
1432  
1278 –
1433    for (l = 0; l < nInfo; l++){
1434      // create the atom and short range interaction arrays
1435  
# Line 1301 | Line 1455 | void SimSetup::makeSysArrays(void){
1455   #else // is_mpi
1456  
1457      molIndex = 0;
1458 <    globalAtomIndex = 0;
1458 >    globalAtomCounter = 0;
1459      for (i = 0; i < n_components; i++){
1460        for (j = 0; j < components_nmol[i]; j++){
1461          the_molecules[molIndex].setStampID(i);
1462          the_molecules[molIndex].setMyIndex(molIndex);
1463          the_molecules[molIndex].setGlobalIndex(molIndex);
1464          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1465 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1466 <          globalAtomIndex++;
1465 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1466 >          globalAtomCounter++;
1467          }
1468          molIndex++;
1469        }
# Line 1318 | Line 1472 | void SimSetup::makeSysArrays(void){
1472  
1473   #endif // is_mpi
1474  
1475 <
1476 <    if (info[l].n_SRI){
1477 <      Exclude::createArray(info[l].n_SRI);
1324 <      the_excludes = new Exclude * [info[l].n_SRI];
1325 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1326 <        the_excludes[ex] = new Exclude(ex);
1327 <      }
1328 <      info[l].globalExcludes = new int;
1329 <      info[l].n_exclude = info[l].n_SRI;
1330 <    }
1331 <    else{
1332 <      Exclude::createArray(1);
1333 <      the_excludes = new Exclude * ;
1334 <      the_excludes[0] = new Exclude(0);
1335 <      the_excludes[0]->setPair(0, 0);
1336 <      info[l].globalExcludes = new int;
1337 <      info[l].globalExcludes[0] = 0;
1338 <      info[l].n_exclude = 0;
1339 <    }
1340 <
1475 >    info[l].globalExcludes = new int;
1476 >    info[l].globalExcludes[0] = 0;
1477 >    
1478      // set the arrays into the SimInfo object
1479  
1480      info[l].atoms = the_atoms;
1481      info[l].molecules = the_molecules;
1482      info[l].nGlobalExcludes = 0;
1346 –    info[l].excludes = the_excludes;
1483  
1484      the_ff->setSimInfo(info);
1485    }
# Line 1387 | Line 1523 | void SimSetup::makeIntegrator(void){
1523          else{
1524            sprintf(painCave.errMsg,
1525                    "SimSetup error: If you use the NVT\n"
1526 <                  "    ensemble, you must set tauThermostat.\n");
1526 >                  "\tensemble, you must set tauThermostat.\n");
1527            painCave.isFatal = 1;
1528            simError();
1529          }
# Line 1410 | Line 1546 | void SimSetup::makeIntegrator(void){
1546          else{
1547            sprintf(painCave.errMsg,
1548                    "SimSetup error: If you use a constant pressure\n"
1549 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1549 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1550            painCave.isFatal = 1;
1551            simError();
1552          }
# Line 1420 | Line 1556 | void SimSetup::makeIntegrator(void){
1556          else{
1557            sprintf(painCave.errMsg,
1558                    "SimSetup error: If you use an NPT\n"
1559 <                  "    ensemble, you must set tauThermostat.\n");
1559 >                  "\tensemble, you must set tauThermostat.\n");
1560            painCave.isFatal = 1;
1561            simError();
1562          }
# Line 1430 | Line 1566 | void SimSetup::makeIntegrator(void){
1566          else{
1567            sprintf(painCave.errMsg,
1568                    "SimSetup error: If you use an NPT\n"
1569 <                  "    ensemble, you must set tauBarostat.\n");
1569 >                  "\tensemble, you must set tauBarostat.\n");
1570            painCave.isFatal = 1;
1571            simError();
1572          }
# Line 1453 | Line 1589 | void SimSetup::makeIntegrator(void){
1589          else{
1590            sprintf(painCave.errMsg,
1591                    "SimSetup error: If you use a constant pressure\n"
1592 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1592 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1593            painCave.isFatal = 1;
1594            simError();
1595          }    
1596  
1597          if (globals->haveTauThermostat())
1598            myNPTf->setTauThermostat(globals->getTauThermostat());
1599 +
1600          else{
1601            sprintf(painCave.errMsg,
1602                    "SimSetup error: If you use an NPT\n"
1603 <                  "    ensemble, you must set tauThermostat.\n");
1603 >                  "\tensemble, you must set tauThermostat.\n");
1604            painCave.isFatal = 1;
1605            simError();
1606          }
1607  
1608          if (globals->haveTauBarostat())
1609            myNPTf->setTauBarostat(globals->getTauBarostat());
1610 +
1611          else{
1612            sprintf(painCave.errMsg,
1613                    "SimSetup error: If you use an NPT\n"
1614 <                  "    ensemble, you must set tauBarostat.\n");
1614 >                  "\tensemble, you must set tauBarostat.\n");
1615            painCave.isFatal = 1;
1616            simError();
1617          }
# Line 1496 | Line 1634 | void SimSetup::makeIntegrator(void){
1634          else{
1635            sprintf(painCave.errMsg,
1636                    "SimSetup error: If you use a constant pressure\n"
1637 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1637 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1638            painCave.isFatal = 1;
1639            simError();
1640          }    
# Line 1506 | Line 1644 | void SimSetup::makeIntegrator(void){
1644          else{
1645            sprintf(painCave.errMsg,
1646                    "SimSetup error: If you use an NPT\n"
1647 <                  "    ensemble, you must set tauThermostat.\n");
1647 >                  "\tensemble, you must set tauThermostat.\n");
1648            painCave.isFatal = 1;
1649            simError();
1650          }
# Line 1516 | Line 1654 | void SimSetup::makeIntegrator(void){
1654          else{
1655            sprintf(painCave.errMsg,
1656                    "SimSetup error: If you use an NPT\n"
1657 <                  "    ensemble, you must set tauBarostat.\n");
1657 >                  "\tensemble, you must set tauBarostat.\n");
1658            painCave.isFatal = 1;
1659            simError();
1660          }
# Line 1569 | Line 1707 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1707    }
1708    else{
1709      sprintf(painCave.errMsg,
1710 <            "ZConstraint error: If you use an ZConstraint\n"
1711 <            " , you must set sample time.\n");
1710 >            "ZConstraint error: If you use a ZConstraint,\n"
1711 >            "\tyou must set zconsTime.\n");
1712      painCave.isFatal = 1;
1713      simError();
1714    }
# Line 1585 | Line 1723 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1723    else{
1724      double defaultZConsTol = 0.01;
1725      sprintf(painCave.errMsg,
1726 <            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1727 <            " , default value %f is used.\n",
1726 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
1727 >            "\tOOPSE will use a default value of %f.\n"
1728 >            "\tTo set the tolerance, use the zconsTol variable.\n",
1729              defaultZConsTol);
1730      painCave.isFatal = 0;
1731      simError();      
# Line 1604 | Line 1743 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1743    }
1744    else{
1745      sprintf(painCave.errMsg,
1746 <            "ZConstraint Warning: User does not set force Subtraction policy, "
1747 <            "PolicyByMass is used\n");
1746 >            "ZConstraint Warning: No force subtraction policy was set.\n"
1747 >            "\tOOPSE will use PolicyByMass.\n"
1748 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
1749      painCave.isFatal = 0;
1750      simError();
1751      zconsForcePolicy->setData("BYMASS");
1752    }
1753  
1754    theInfo.addProperty(zconsForcePolicy);
1755 +
1756 +  //set zcons gap
1757 +  DoubleData* zconsGap = new DoubleData();
1758 +  zconsGap->setID(ZCONSGAP_ID);
1759 +
1760 +  if (globals->haveZConsGap()){
1761 +    zconsGap->setData(globals->getZconsGap());
1762 +    theInfo.addProperty(zconsGap);  
1763 +  }
1764 +
1765 +  //set zcons fixtime
1766 +  DoubleData* zconsFixtime = new DoubleData();
1767 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1768  
1769 +  if (globals->haveZConsFixTime()){
1770 +    zconsFixtime->setData(globals->getZconsFixtime());
1771 +    theInfo.addProperty(zconsFixtime);  
1772 +  }
1773 +
1774 +  //set zconsUsingSMD
1775 +  IntData* zconsUsingSMD = new IntData();
1776 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1777 +
1778 +  if (globals->haveZConsUsingSMD()){
1779 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1780 +    theInfo.addProperty(zconsUsingSMD);  
1781 +  }
1782 +
1783    //Determine the name of ouput file and add it into SimInfo's property list
1784    //Be careful, do not use inFileName, since it is a pointer which
1785    //point to a string at master node, and slave nodes do not contain that string
# Line 1642 | Line 1809 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1809      tempParaItem.zPos = zconStamp[i]->getZpos();
1810      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1811      tempParaItem.kRatio = zconStamp[i]->getKratio();
1812 <
1812 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1813 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1814      zconsParaData->addItem(tempParaItem);
1815    }
1816  
1817    //check the uniqueness of index  
1818    if(!zconsParaData->isIndexUnique()){
1819      sprintf(painCave.errMsg,
1820 <            "ZConstraint Error: molIndex is not unique\n");
1820 >            "ZConstraint Error: molIndex is not unique!\n");
1821      painCave.isFatal = 1;
1822      simError();
1823    }
# Line 1660 | Line 1828 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1828    //push data into siminfo, therefore, we can retrieve later
1829    theInfo.addProperty(zconsParaData);
1830   }
1831 +
1832 + void SimSetup::makeMinimizer(){
1833 +
1834 +  OOPSEMinimizer* myOOPSEMinimizer;
1835 +  MinimizerParameterSet* param;
1836 +  char minimizerName[100];
1837 +  
1838 +  for (int i = 0; i < nInfo; i++){
1839 +    
1840 +    //prepare parameter set for minimizer
1841 +    param = new MinimizerParameterSet();
1842 +    param->setDefaultParameter();
1843 +
1844 +    if (globals->haveMinimizer()){
1845 +      param->setFTol(globals->getMinFTol());
1846 +    }
1847 +
1848 +    if (globals->haveMinGTol()){
1849 +      param->setGTol(globals->getMinGTol());
1850 +    }
1851 +
1852 +    if (globals->haveMinMaxIter()){
1853 +      param->setMaxIteration(globals->getMinMaxIter());
1854 +    }
1855 +
1856 +    if (globals->haveMinWriteFrq()){
1857 +      param->setMaxIteration(globals->getMinMaxIter());
1858 +    }
1859 +
1860 +    if (globals->haveMinWriteFrq()){
1861 +      param->setWriteFrq(globals->getMinWriteFrq());
1862 +    }
1863 +    
1864 +    if (globals->haveMinStepSize()){
1865 +      param->setStepSize(globals->getMinStepSize());
1866 +    }
1867 +
1868 +    if (globals->haveMinLSMaxIter()){
1869 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
1870 +    }    
1871 +
1872 +    if (globals->haveMinLSTol()){
1873 +      param->setLineSearchTol(globals->getMinLSTol());
1874 +    }    
1875 +
1876 +    strcpy(minimizerName, globals->getMinimizer());
1877 +
1878 +    if (!strcasecmp(minimizerName, "CG")){
1879 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1880 +    }
1881 +    else if (!strcasecmp(minimizerName, "SD")){
1882 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1883 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1884 +    }
1885 +    else{
1886 +          sprintf(painCave.errMsg,
1887 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1888 +          painCave.isFatal = 0;
1889 +          simError();
1890 +
1891 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1892 +    }
1893 +     info[i].the_integrator = myOOPSEMinimizer;
1894 +
1895 +     //store the minimizer into simInfo
1896 +     info[i].the_minimizer = myOOPSEMinimizer;
1897 +     info[i].has_minimizer = true;
1898 +  }
1899 +
1900 + }

Diff Legend

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