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 845 by gezelter, Thu Oct 30 18:59:20 2003 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  
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 103 | Line 131 | void SimSetup::createSim(void){
131  
132    sysObjectsCreation();
133  
134 +  // check on the post processing info
135 +
136 +  finalInfoCheck();
137 +
138    // initialize the system coordinates
139  
140    if ( !initSuspend ){
# Line 112 | Line 144 | void SimSetup::createSim(void){
144        info[0].currentTime = 0.0;
145    }  
146  
115 –  // check on the post processing info
116 –
117 –  finalInfoCheck();
118 –
147    // make the output filenames
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;
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  
158 –
198    // init the atoms
199  
200 <  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 172 | 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]);
217 <      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
218 <      molInfo.myBonds = new Bond * [molInfo.nBonds];
219 <      molInfo.myBends = new Bend * [molInfo.nBends];
220 <      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());
245            info[k].n_oriented++;
246            molInfo.myAtoms[j] = dAtom;
247  
248 <          ux = currentAtom->getOrntX();
249 <          uy = currentAtom->getOrntY();
250 <          uz = currentAtom->getOrntZ();
248 >          // Directional Atoms have standard unit vectors which are oriented
249 >          // in space using the three Euler angles.  We assume the standard
250 >          // unit vector was originally along the z axis below.
251  
252 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
252 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
253 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
254 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
255  
256 <          u = sqrt(uSqr);
257 <          ux = ux / u;
205 <          uy = uy / u;
206 <          uz = uz / u;
207 <
208 <          dAtom->setSUx(ux);
209 <          dAtom->setSUy(uy);
210 <          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 228 | 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  
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 –        }
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);
247 < #else  // isn't MPI
248 <
249 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
250 < #endif  //is_mpi
290 >        info[k].excludes->addPair(exI, exJ);
291        }
252 –      excludeOffset += molInfo.nBonds;
292  
293        //make the bends
294        for (j = 0; j < molInfo.nBends; j++){
# Line 299 | 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{
346 <          exI = theBends[j].a;
347 <          exJ = theBends[j].b;
348 <        }
341 >        if (theBends[j].isGhost) {
342 >          
343 >          tempI = theBends[j].a;
344 >          tempJ = theBends[j].b;
345 >          
346 > #ifdef IS_MPI
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 <        // exclude_I must always be the smaller of the pair
312 <        if (exI > exJ){
313 <          tempEx = exI;
314 <          exI = exJ;
315 <          exJ = tempEx;
316 <        }
317 < #ifdef IS_MPI
318 <        tempEx = exI;
319 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
320 <        tempEx = exJ;
321 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
355 >        } else {
356  
357 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
358 < #else  // isn't MPI
359 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
360 < #endif  //is_mpi
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        }
328 –      excludeOffset += molInfo.nBends;
376  
377        for (j = 0; j < molInfo.nTorsions; j++){
378          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 334 | 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  
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 –        }
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        }
357 –      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  
368 –
478        info[k].molecules[i].initialize(molInfo);
479  
480  
# Line 374 | 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 383 | 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 585 | Line 723 | void SimSetup::gatherInfo(void){
723    else if (!strcasecmp(force_field, "EAM")){
724      ffCase = FF_EAM;
725    }
726 +  else if (!strcasecmp(force_field, "WATER")){
727 +    ffCase = FF_H2O;
728 +  }
729    else{
730      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
731              force_field);
# Line 613 | Line 754 | void SimSetup::gatherInfo(void){
754    }
755    else{
756      sprintf(painCave.errMsg,
757 <            "SimSetup Warning. Unrecognized Ensemble -> %s, "
758 <            "reverting to NVE for this simulation.\n",
757 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
758 >            "\treverting to NVE for this simulation.\n",
759              ensemble);
760           painCave.isFatal = 0;
761           simError();
# Line 646 | Line 787 | void SimSetup::gatherInfo(void){
787        if (!the_components[i]->haveNMol()){
788          // we have a problem
789          sprintf(painCave.errMsg,
790 <                "SimSetup Error. No global NMol or component NMol"
791 <                " given. Cannot calculate the number of atoms.\n");
790 >                "SimSetup Error. No global NMol or component NMol given.\n"
791 >                "\tCannot calculate the number of atoms.\n");
792          painCave.isFatal = 1;
793          simError();
794        }
# Line 665 | Line 806 | void SimSetup::gatherInfo(void){
806              " Please give nMol in the components.\n");
807      painCave.isFatal = 1;
808      simError();
809 +  }
810 +
811 +  //check whether sample time, status time, thermal time and reset time are divisble by dt
812 +  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
813 +    sprintf(painCave.errMsg,
814 +            "Sample time is not divisible by dt.\n"
815 +            "\tThis will result in samples that are not uniformly\n"
816 +            "\tdistributed in time.  If this is a problem, change\n"
817 +            "\tyour sampleTime variable.\n");
818 +    painCave.isFatal = 0;
819 +    simError();    
820    }
821  
822 +  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
823 +    sprintf(painCave.errMsg,
824 +            "Status time is not divisible by dt.\n"
825 +            "\tThis will result in status reports that are not uniformly\n"
826 +            "\tdistributed in time.  If this is a problem, change \n"
827 +            "\tyour statusTime variable.\n");
828 +    painCave.isFatal = 0;
829 +    simError();    
830 +  }
831 +
832 +  if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
833 +    sprintf(painCave.errMsg,
834 +            "Thermal time is not divisible by dt.\n"
835 +            "\tThis will result in thermalizations that are not uniformly\n"
836 +            "\tdistributed in time.  If this is a problem, change \n"
837 +            "\tyour thermalTime variable.\n");
838 +    painCave.isFatal = 0;
839 +    simError();    
840 +  }  
841 +
842 +  if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
843 +    sprintf(painCave.errMsg,
844 +            "Reset time is not divisible by dt.\n"
845 +            "\tThis will result in integrator resets that are not uniformly\n"
846 +            "\tdistributed in time.  If this is a problem, change\n"
847 +            "\tyour resetTime variable.\n");
848 +    painCave.isFatal = 0;
849 +    simError();    
850 +  }
851 +
852    // set the status, sample, and thermal kick times
853  
854    for (i = 0; i < nInfo; i++){
# Line 699 | Line 881 | void SimSetup::gatherInfo(void){
881      
882      if (globals->haveTempSet())
883        info[i].setTemp = globals->getTempSet();
884 +
885 +    // check for the extended State init
886 +
887 +    info[i].useInitXSstate = globals->getUseInitXSstate();
888 +    info[i].orthoTolerance = globals->getOrthoBoxTolerance();
889      
890    }
891    
# Line 741 | Line 928 | void SimSetup::gatherInfo(void){
928    for (int i = 0; i < nInfo; i++){
929      info[i].setSeed(seedValue);
930    }
931 <
931 >  
932   #ifdef IS_MPI
933 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
933 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
934    MPIcheckPoint();
935   #endif // is_mpi
936   }
# Line 776 | Line 963 | void SimSetup::finalInfoCheck(void){
963  
964        if (!globals->haveECR()){
965          sprintf(painCave.errMsg,
966 <                "SimSetup Warning: using default value of 1/2 the smallest "
967 <                "box length for the electrostaticCutoffRadius.\n"
968 <                "I hope you have a very fast processor!\n");
966 >                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
967 >                "\tOOPSE will use a default value of 15.0 angstroms"
968 >                "\tfor the electrostaticCutoffRadius.\n");
969          painCave.isFatal = 0;
970          simError();
971 <        double smallest;
785 <        smallest = info[i].boxL[0];
786 <        if (info[i].boxL[1] <= smallest)
787 <          smallest = info[i].boxL[1];
788 <        if (info[i].boxL[2] <= smallest)
789 <          smallest = info[i].boxL[2];
790 <        theEcr = 0.5 * smallest;
971 >        theEcr = 15.0;
972        }
973        else{
974          theEcr = globals->getECR();
# Line 795 | Line 976 | void SimSetup::finalInfoCheck(void){
976  
977        if (!globals->haveEST()){
978          sprintf(painCave.errMsg,
979 <                "SimSetup Warning: using default value of 0.05 * the "
980 <                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
979 >                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
980 >                "\tOOPSE will use a default value of\n"
981 >                "\t0.05 * electrostaticCutoffRadius\n"
982 >                "\tfor the electrostaticSkinThickness\n");
983          painCave.isFatal = 0;
984          simError();
985          theEst = 0.05 * theEcr;
# Line 809 | Line 992 | void SimSetup::finalInfoCheck(void){
992  
993        if (!globals->haveDielectric()){
994          sprintf(painCave.errMsg,
995 <                "SimSetup Error: You are trying to use Reaction Field without"
996 <                "setting a dielectric constant!\n");
995 >                "SimSetup Error: No Dielectric constant was set.\n"
996 >                "\tYou are trying to use Reaction Field without"
997 >                "\tsetting a dielectric constant!\n");
998          painCave.isFatal = 1;
999          simError();
1000        }
# Line 820 | Line 1004 | void SimSetup::finalInfoCheck(void){
1004        if (usesDipoles){
1005          if (!globals->haveECR()){
1006            sprintf(painCave.errMsg,
1007 <                  "SimSetup Warning: using default value of 1/2 the smallest "
1008 <                  "box length for the electrostaticCutoffRadius.\n"
1009 <                  "I hope you have a very fast processor!\n");
1010 <          painCave.isFatal = 0;
1011 <          simError();
1012 <          double smallest;
829 <          smallest = info[i].boxL[0];
830 <          if (info[i].boxL[1] <= smallest)
831 <            smallest = info[i].boxL[1];
832 <          if (info[i].boxL[2] <= smallest)
833 <            smallest = info[i].boxL[2];
834 <          theEcr = 0.5 * smallest;
1007 >                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1008 >                  "\tOOPSE will use a default value of 15.0 angstroms"
1009 >                  "\tfor the electrostaticCutoffRadius.\n");
1010 >          painCave.isFatal = 0;
1011 >          simError();
1012 >          theEcr = 15.0;
1013          }
1014          else{
1015            theEcr = globals->getECR();
1016          }
1017 <
1017 >        
1018          if (!globals->haveEST()){
1019            sprintf(painCave.errMsg,
1020 <                  "SimSetup Warning: using default value of 0.05 * the "
1021 <                  "electrostaticCutoffRadius for the "
1022 <                  "electrostaticSkinThickness\n");
1020 >                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1021 >                  "\tOOPSE will use a default value of\n"
1022 >                  "\t0.05 * electrostaticCutoffRadius\n"
1023 >                  "\tfor the electrostaticSkinThickness\n");
1024            painCave.isFatal = 0;
1025            simError();
1026            theEst = 0.05 * theEcr;
# Line 849 | Line 1028 | void SimSetup::finalInfoCheck(void){
1028          else{
1029            theEst = globals->getEST();
1030          }
1031 <
1031 >        
1032          info[i].setDefaultEcr(theEcr, theEst);
1033        }
1034      }
856 –    info[i].checkCutOffs();
1035    }
858 –
1036   #ifdef IS_MPI
1037    strcpy(checkPointMsg, "post processing checks out");
1038    MPIcheckPoint();
1039   #endif // is_mpi
1040   }
1041 <
1041 >  
1042   void SimSetup::initSystemCoords(void){
1043    int i;
1044  
# Line 889 | Line 1066 | void SimSetup::initSystemCoords(void){
1066      delete fileInit;
1067    }
1068    else{
1069 < #ifdef IS_MPI
893 <
1069 >    
1070      // no init from bass
1071 <
1071 >    
1072      sprintf(painCave.errMsg,
1073 <            "Cannot intialize a parallel simulation without an initial configuration file.\n");
1073 >            "Cannot intialize a simulation without an initial configuration file.\n");
1074      painCave.isFatal = 1;;
1075      simError();
1076 <
901 < #else
902 <
903 <    initFromBass();
904 <
905 <
906 < #endif
1076 >    
1077    }
1078  
1079   #ifdef IS_MPI
# Line 1057 | Line 1227 | void SimSetup::createFF(void){
1227        the_ff = new EAM_FF();
1228        break;
1229  
1230 +    case FF_H2O:
1231 +      the_ff = new WATER();
1232 +      break;
1233 +
1234      default:
1235        sprintf(painCave.errMsg,
1236                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1129 | 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 1157 | 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 1173 | 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  
1178 –
1356    for (i = 0; i < n_components; i++){
1357      for (j = 0; j < components_nmol[i]; j++){
1358        if (mol2proc[allMol] == worldRank){
# Line 1183 | 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 1196 | 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,
1381 <            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1382 <            " localAtom (%d) are not equal.\n",
1381 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1382 >            "\tlocalAtom (%d) are not equal.\n",
1383              info[0].n_atoms, local_atoms);
1384      painCave.isFatal = 1;
1385      simError();
# Line 1228 | Line 1407 | void SimSetup::makeSysArrays(void){
1407  
1408    Atom** the_atoms;
1409    Molecule* the_molecules;
1231 –  Exclude** the_excludes;
1410  
1233 –
1411    for (l = 0; l < nInfo; l++){
1412      // create the atom and short range interaction arrays
1413  
# Line 1273 | 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);
1279 <      the_excludes = new Exclude * [info[l].n_SRI];
1280 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1281 <        the_excludes[ex] = new Exclude(ex);
1282 <      }
1283 <      info[l].globalExcludes = new int;
1284 <      info[l].n_exclude = info[l].n_SRI;
1285 <    }
1286 <    else{
1287 <      Exclude::createArray(1);
1288 <      the_excludes = new Exclude * ;
1289 <      the_excludes[0] = new Exclude(0);
1290 <      the_excludes[0]->setPair(0, 0);
1291 <      info[l].globalExcludes = new int;
1292 <      info[l].globalExcludes[0] = 0;
1293 <      info[l].n_exclude = 0;
1294 <    }
1295 <
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;
1301 –    info[l].excludes = the_excludes;
1461  
1462      the_ff->setSimInfo(info);
1463    }
# Line 1342 | Line 1501 | void SimSetup::makeIntegrator(void){
1501          else{
1502            sprintf(painCave.errMsg,
1503                    "SimSetup error: If you use the NVT\n"
1504 <                  "    ensemble, you must set tauThermostat.\n");
1504 >                  "\tensemble, you must set tauThermostat.\n");
1505            painCave.isFatal = 1;
1506            simError();
1507          }
# Line 1365 | Line 1524 | void SimSetup::makeIntegrator(void){
1524          else{
1525            sprintf(painCave.errMsg,
1526                    "SimSetup error: If you use a constant pressure\n"
1527 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1527 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1528            painCave.isFatal = 1;
1529            simError();
1530          }
# Line 1375 | Line 1534 | void SimSetup::makeIntegrator(void){
1534          else{
1535            sprintf(painCave.errMsg,
1536                    "SimSetup error: If you use an NPT\n"
1537 <                  "    ensemble, you must set tauThermostat.\n");
1537 >                  "\tensemble, you must set tauThermostat.\n");
1538            painCave.isFatal = 1;
1539            simError();
1540          }
# Line 1385 | Line 1544 | void SimSetup::makeIntegrator(void){
1544          else{
1545            sprintf(painCave.errMsg,
1546                    "SimSetup error: If you use an NPT\n"
1547 <                  "    ensemble, you must set tauBarostat.\n");
1547 >                  "\tensemble, you must set tauBarostat.\n");
1548            painCave.isFatal = 1;
1549            simError();
1550          }
# Line 1408 | Line 1567 | void SimSetup::makeIntegrator(void){
1567          else{
1568            sprintf(painCave.errMsg,
1569                    "SimSetup error: If you use a constant pressure\n"
1570 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1570 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1571            painCave.isFatal = 1;
1572            simError();
1573          }    
# Line 1419 | Line 1578 | void SimSetup::makeIntegrator(void){
1578          else{
1579            sprintf(painCave.errMsg,
1580                    "SimSetup error: If you use an NPT\n"
1581 <                  "    ensemble, you must set tauThermostat.\n");
1581 >                  "\tensemble, you must set tauThermostat.\n");
1582            painCave.isFatal = 1;
1583            simError();
1584          }
# Line 1430 | Line 1589 | void SimSetup::makeIntegrator(void){
1589          else{
1590            sprintf(painCave.errMsg,
1591                    "SimSetup error: If you use an NPT\n"
1592 <                  "    ensemble, you must set tauBarostat.\n");
1592 >                  "\tensemble, you must set tauBarostat.\n");
1593            painCave.isFatal = 1;
1594            simError();
1595          }
# Line 1453 | Line 1612 | void SimSetup::makeIntegrator(void){
1612          else{
1613            sprintf(painCave.errMsg,
1614                    "SimSetup error: If you use a constant pressure\n"
1615 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1615 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1616            painCave.isFatal = 1;
1617            simError();
1618          }    
# Line 1463 | Line 1622 | void SimSetup::makeIntegrator(void){
1622          else{
1623            sprintf(painCave.errMsg,
1624                    "SimSetup error: If you use an NPT\n"
1625 <                  "    ensemble, you must set tauThermostat.\n");
1625 >                  "\tensemble, you must set tauThermostat.\n");
1626            painCave.isFatal = 1;
1627            simError();
1628          }
# Line 1473 | Line 1632 | void SimSetup::makeIntegrator(void){
1632          else{
1633            sprintf(painCave.errMsg,
1634                    "SimSetup error: If you use an NPT\n"
1635 <                  "    ensemble, you must set tauBarostat.\n");
1635 >                  "\tensemble, you must set tauBarostat.\n");
1636            painCave.isFatal = 1;
1637            simError();
1638          }
# Line 1526 | Line 1685 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1685    }
1686    else{
1687      sprintf(painCave.errMsg,
1688 <            "ZConstraint error: If you use an ZConstraint\n"
1689 <            " , you must set sample time.\n");
1688 >            "ZConstraint error: If you use a ZConstraint,\n"
1689 >            "\tyou must set zconsTime.\n");
1690      painCave.isFatal = 1;
1691      simError();
1692    }
# Line 1542 | Line 1701 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1701    else{
1702      double defaultZConsTol = 0.01;
1703      sprintf(painCave.errMsg,
1704 <            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1705 <            " , default value %f is used.\n",
1704 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
1705 >            "\tOOPSE will use a default value of %f.\n"
1706 >            "\tTo set the tolerance, use the zconsTol variable.\n",
1707              defaultZConsTol);
1708      painCave.isFatal = 0;
1709      simError();      
# Line 1561 | Line 1721 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1721    }
1722    else{
1723      sprintf(painCave.errMsg,
1724 <            "ZConstraint Warning: User does not set force Subtraction policy, "
1725 <            "PolicyByMass is used\n");
1724 >            "ZConstraint Warning: No force subtraction policy was set.\n"
1725 >            "\tOOPSE will use PolicyByMass.\n"
1726 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
1727      painCave.isFatal = 0;
1728      simError();
1729      zconsForcePolicy->setData("BYMASS");
1730    }
1731  
1732    theInfo.addProperty(zconsForcePolicy);
1733 +
1734 +  //set zcons gap
1735 +  DoubleData* zconsGap = new DoubleData();
1736 +  zconsGap->setID(ZCONSGAP_ID);
1737 +
1738 +  if (globals->haveZConsGap()){
1739 +    zconsGap->setData(globals->getZconsGap());
1740 +    theInfo.addProperty(zconsGap);  
1741 +  }
1742 +
1743 +  //set zcons fixtime
1744 +  DoubleData* zconsFixtime = new DoubleData();
1745 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1746 +
1747 +  if (globals->haveZConsFixTime()){
1748 +    zconsFixtime->setData(globals->getZconsFixtime());
1749 +    theInfo.addProperty(zconsFixtime);  
1750 +  }
1751 +
1752 +  //set zconsUsingSMD
1753 +  IntData* zconsUsingSMD = new IntData();
1754 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1755 +
1756 +  if (globals->haveZConsUsingSMD()){
1757 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1758 +    theInfo.addProperty(zconsUsingSMD);  
1759 +  }
1760  
1761    //Determine the name of ouput file and add it into SimInfo's property list
1762    //Be careful, do not use inFileName, since it is a pointer which
# Line 1599 | Line 1787 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1787      tempParaItem.zPos = zconStamp[i]->getZpos();
1788      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1789      tempParaItem.kRatio = zconStamp[i]->getKratio();
1790 <
1790 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1791 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1792      zconsParaData->addItem(tempParaItem);
1793    }
1794  
1795    //check the uniqueness of index  
1796    if(!zconsParaData->isIndexUnique()){
1797      sprintf(painCave.errMsg,
1798 <            "ZConstraint Error: molIndex is not unique\n");
1798 >            "ZConstraint Error: molIndex is not unique!\n");
1799      painCave.isFatal = 1;
1800      simError();
1801    }
# Line 1617 | Line 1806 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1806    //push data into siminfo, therefore, we can retrieve later
1807    theInfo.addProperty(zconsParaData);
1808   }
1809 +
1810 + void SimSetup::makeMinimizer(){
1811 +
1812 +  OOPSEMinimizer* myOOPSEMinimizer;
1813 +  MinimizerParameterSet* param;
1814 +  char minimizerName[100];
1815 +  
1816 +  for (int i = 0; i < nInfo; i++){
1817 +    
1818 +    //prepare parameter set for minimizer
1819 +    param = new MinimizerParameterSet();
1820 +    param->setDefaultParameter();
1821 +
1822 +    if (globals->haveMinimizer()){
1823 +      param->setFTol(globals->getMinFTol());
1824 +    }
1825 +
1826 +    if (globals->haveMinGTol()){
1827 +      param->setGTol(globals->getMinGTol());
1828 +    }
1829 +
1830 +    if (globals->haveMinMaxIter()){
1831 +      param->setMaxIteration(globals->getMinMaxIter());
1832 +    }
1833 +
1834 +    if (globals->haveMinWriteFrq()){
1835 +      param->setMaxIteration(globals->getMinMaxIter());
1836 +    }
1837 +
1838 +    if (globals->haveMinWriteFrq()){
1839 +      param->setWriteFrq(globals->getMinWriteFrq());
1840 +    }
1841 +    
1842 +    if (globals->haveMinStepSize()){
1843 +      param->setStepSize(globals->getMinStepSize());
1844 +    }
1845 +
1846 +    if (globals->haveMinLSMaxIter()){
1847 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
1848 +    }    
1849 +
1850 +    if (globals->haveMinLSTol()){
1851 +      param->setLineSearchTol(globals->getMinLSTol());
1852 +    }    
1853 +
1854 +    strcpy(minimizerName, globals->getMinimizer());
1855 +
1856 +    if (!strcasecmp(minimizerName, "CG")){
1857 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1858 +    }
1859 +    else if (!strcasecmp(minimizerName, "SD")){
1860 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1861 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1862 +    }
1863 +    else{
1864 +          sprintf(painCave.errMsg,
1865 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1866 +          painCave.isFatal = 0;
1867 +          simError();
1868 +
1869 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1870 +    }
1871 +     info[i].the_integrator = myOOPSEMinimizer;
1872 +
1873 +     //store the minimizer into simInfo
1874 +     info[i].the_minimizer = myOOPSEMinimizer;
1875 +     info[i].has_minimizer = true;
1876 +  }
1877 +
1878 + }

Diff Legend

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