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 962 by tim, Mon Jan 19 18:36:21 2004 UTC vs.
Revision 1129 by tim, Thu Apr 22 03:29:30 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 +  RigidBodyStamp* currentRigidBody;
186  
187    bond_pair* theBonds;
188    bend_set* theBends;
189    torsion_set* theTorsions;
190  
191 +  set<int> skipList;
192  
193 +  double phi, theta, psi;
194 +  char* molName;
195 +  char rbName[100];
196 +
197    //init the forceField paramters
198  
199    the_ff->readParams();
200  
182 –
201    // init the atoms
202  
203 <  double ux, uy, uz, u, uSqr;
203 >  int nMembers, nNew, rb1, rb2;
204  
205    for (k = 0; k < nInfo; k++){
206      the_ff->setSimInfo(&(info[k]));
207  
208      atomOffset = 0;
209 <    excludeOffset = 0;
209 >
210      for (i = 0; i < info[k].n_mol; i++){
211        stampID = info[k].molecules[i].getStampID();
212 +      molName = comp_stamps[stampID]->getID();
213  
214        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
215        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
216        molInfo.nBends = comp_stamps[stampID]->getNBends();
217        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
218 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
219 <
218 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
219 >      
220        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
202 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
203 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
204 –      molInfo.myBends = new Bend * [molInfo.nBends];
205 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
221  
222 +      if (molInfo.nBonds > 0)
223 +        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
224 +      else
225 +        molInfo.myBonds = NULL;
226 +
227 +      if (molInfo.nBends > 0)
228 +        molInfo.myBends = new (Bend *) [molInfo.nBends];
229 +      else
230 +        molInfo.myBends = NULL;
231 +
232 +      if (molInfo.nTorsions > 0)
233 +        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
234 +      else
235 +        molInfo.myTorsions = NULL;
236 +
237        theBonds = new bond_pair[molInfo.nBonds];
238        theBends = new bend_set[molInfo.nBends];
239        theTorsions = new torsion_set[molInfo.nTorsions];
240 <
240 >      
241        // make the Atoms
242  
243        for (j = 0; j < molInfo.nAtoms; j++){
244          currentAtom = comp_stamps[stampID]->getAtom(j);
245 +
246          if (currentAtom->haveOrientation()){
247            dAtom = new DirectionalAtom((j + atomOffset),
248                                        info[k].getConfiguration());
249            info[k].n_oriented++;
250            molInfo.myAtoms[j] = dAtom;
251  
252 <          ux = currentAtom->getOrntX();
253 <          uy = currentAtom->getOrntY();
254 <          uz = currentAtom->getOrntZ();
252 >          // Directional Atoms have standard unit vectors which are oriented
253 >          // in space using the three Euler angles.  We assume the standard
254 >          // unit vector was originally along the z axis below.
255  
256 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
257 <
258 <          u = sqrt(uSqr);
259 <          ux = ux / u;
260 <          uy = uy / u;
261 <          uz = uz / u;
231 <
232 <          dAtom->setSUx(ux);
233 <          dAtom->setSUy(uy);
234 <          dAtom->setSUz(uz);
256 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
257 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
258 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
259 >
260 >          dAtom->setUnitFrameFromEuler(phi, theta, psi);
261 >            
262          }
263          else{
264 <          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
265 <                                               info[k].getConfiguration());
264 >
265 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
266 >
267          }
240 –        molInfo.myAtoms[j]->setType(currentAtom->getType());
268  
269 +        molInfo.myAtoms[j]->setType(currentAtom->getType());
270   #ifdef IS_MPI
271  
272 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
272 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
273  
274   #endif // is_mpi
275        }
# Line 252 | Line 280 | void SimSetup::makeMolecules(void){
280          theBonds[j].a = currentBond->getA() + atomOffset;
281          theBonds[j].b = currentBond->getB() + atomOffset;
282  
283 <        exI = theBonds[j].a;
284 <        exJ = theBonds[j].b;
283 >        tempI = theBonds[j].a;
284 >        tempJ = theBonds[j].b;
285  
258 –        // exclude_I must always be the smaller of the pair
259 –        if (exI > exJ){
260 –          tempEx = exI;
261 –          exI = exJ;
262 –          exJ = tempEx;
263 –        }
286   #ifdef IS_MPI
287 <        tempEx = exI;
288 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
289 <        tempEx = exJ;
290 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
287 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
288 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
289 > #else
290 >        exI = tempI + 1;
291 >        exJ = tempJ + 1;
292 > #endif
293  
294 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
271 < #else  // isn't MPI
272 <
273 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
274 < #endif  //is_mpi
294 >        info[k].excludes->addPair(exI, exJ);
295        }
276 –      excludeOffset += molInfo.nBonds;
296  
297        //make the bends
298        for (j = 0; j < molInfo.nBends; j++){
# Line 323 | Line 342 | void SimSetup::makeMolecules(void){
342            }
343          }
344  
345 <        if (!theBends[j].isGhost){
346 <          exI = theBends[j].a;
347 <          exJ = theBends[j].c;
348 <        }
349 <        else{
331 <          exI = theBends[j].a;
332 <          exJ = theBends[j].b;
333 <        }
334 <
335 <        // exclude_I must always be the smaller of the pair
336 <        if (exI > exJ){
337 <          tempEx = exI;
338 <          exI = exJ;
339 <          exJ = tempEx;
340 <        }
345 >        if (theBends[j].isGhost) {
346 >          
347 >          tempI = theBends[j].a;
348 >          tempJ = theBends[j].b;
349 >          
350   #ifdef IS_MPI
351 <        tempEx = exI;
352 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
353 <        tempEx = exJ;
354 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
351 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
352 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
353 > #else
354 >          exI = tempI + 1;
355 >          exJ = tempJ + 1;
356 > #endif          
357 >          info[k].excludes->addPair(exI, exJ);
358  
359 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
360 < #else  // isn't MPI
361 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
362 < #endif  //is_mpi
359 >        } else {
360 >
361 >          tempI = theBends[j].a;
362 >          tempJ = theBends[j].b;
363 >          tempK = theBends[j].c;
364 >          
365 > #ifdef IS_MPI
366 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
367 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
368 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
369 > #else
370 >          exI = tempI + 1;
371 >          exJ = tempJ + 1;
372 >          exK = tempK + 1;
373 > #endif
374 >          
375 >          info[k].excludes->addPair(exI, exK);
376 >          info[k].excludes->addPair(exI, exJ);
377 >          info[k].excludes->addPair(exJ, exK);
378 >        }
379        }
352 –      excludeOffset += molInfo.nBends;
380  
381        for (j = 0; j < molInfo.nTorsions; j++){
382          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 358 | Line 385 | void SimSetup::makeMolecules(void){
385          theTorsions[j].c = currentTorsion->getC() + atomOffset;
386          theTorsions[j].d = currentTorsion->getD() + atomOffset;
387  
388 <        exI = theTorsions[j].a;
389 <        exJ = theTorsions[j].d;
388 >        tempI = theTorsions[j].a;      
389 >        tempJ = theTorsions[j].b;
390 >        tempK = theTorsions[j].c;
391 >        tempL = theTorsions[j].d;
392  
393 <        // exclude_I must always be the smaller of the pair
394 <        if (exI > exJ){
395 <          tempEx = exI;
396 <          exI = exJ;
397 <          exJ = tempEx;
393 > #ifdef IS_MPI
394 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
395 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
396 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
397 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
398 > #else
399 >        exI = tempI + 1;
400 >        exJ = tempJ + 1;
401 >        exK = tempK + 1;
402 >        exL = tempL + 1;
403 > #endif
404 >
405 >        info[k].excludes->addPair(exI, exJ);
406 >        info[k].excludes->addPair(exI, exK);
407 >        info[k].excludes->addPair(exI, exL);        
408 >        info[k].excludes->addPair(exJ, exK);
409 >        info[k].excludes->addPair(exJ, exL);
410 >        info[k].excludes->addPair(exK, exL);
411 >      }
412 >
413 >      
414 >      molInfo.myRigidBodies.clear();
415 >      
416 >      for (j = 0; j < molInfo.nRigidBodies; j++){
417 >
418 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
419 >        nMembers = currentRigidBody->getNMembers();
420 >
421 >        // Create the Rigid Body:
422 >
423 >        myRB = new RigidBody();
424 >
425 >        sprintf(rbName,"%s_RB_%d", molName, j);
426 >        myRB->setType(rbName);
427 >        
428 >        for (rb1 = 0; rb1 < nMembers; rb1++) {
429 >
430 >          // molI is atom numbering inside this molecule
431 >          molI = currentRigidBody->getMember(rb1);    
432 >
433 >          // tempI is atom numbering on local processor
434 >          tempI = molI + atomOffset;
435 >
436 >          // currentAtom is the AtomStamp (which we need for
437 >          // rigid body reference positions)
438 >          currentAtom = comp_stamps[stampID]->getAtom(molI);
439 >
440 >          // When we add to the rigid body, add the atom itself and
441 >          // the stamp info:
442 >
443 >          myRB->addAtom(info[k].atoms[tempI], currentAtom);
444 >          
445 >          // Add this atom to the Skip List for the integrators
446 > #ifdef IS_MPI
447 >          slI = info[k].atoms[tempI]->getGlobalIndex();
448 > #else
449 >          slI = tempI;
450 > #endif
451 >          skipList.insert(slI);
452 >          
453          }
454 +        
455 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
456 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
457 +            
458 +            tempI = currentRigidBody->getMember(rb1);
459 +            tempJ = currentRigidBody->getMember(rb2);
460 +            
461 +            // Some explanation is required here.
462 +            // Fortran indexing starts at 1, while c indexing starts at 0
463 +            // Also, in parallel computations, the GlobalIndex is
464 +            // used for the exclude list:
465 +            
466   #ifdef IS_MPI
467 <        tempEx = exI;
468 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
469 <        tempEx = exJ;
470 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
467 >            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
468 >            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
469 > #else
470 >            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
471 >            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
472 > #endif
473 >            
474 >            info[k].excludes->addPair(exI, exJ);
475 >            
476 >          }
477 >        }
478  
479 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
480 < #else  // isn't MPI
378 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
379 < #endif  //is_mpi
479 >        molInfo.myRigidBodies.push_back(myRB);
480 >        info[k].rigidBodies.push_back(myRB);
481        }
482 <      excludeOffset += molInfo.nTorsions;
482 >      
483  
484 +      // After this is all set up, scan through the atoms to
485 +      // see if they can be added to the integrableObjects:
486  
487 <      // send the arrays off to the forceField for init.
487 >      molInfo.myIntegrableObjects.clear();
488 >      
489 >
490 >      for (j = 0; j < molInfo.nAtoms; j++){
491 >
492 > #ifdef IS_MPI
493 >        slJ = molInfo.myAtoms[j]->getGlobalIndex();
494 > #else
495 >        slJ = j+atomOffset;
496 > #endif
497 >
498 >        // if they aren't on the skip list, then they can be integrated
499 >
500 >        if (skipList.find(slJ) == skipList.end()) {
501 >          mySD = (StuntDouble *) molInfo.myAtoms[j];
502 >          info[k].integrableObjects.push_back(mySD);
503 >          molInfo.myIntegrableObjects.push_back(mySD);
504 >        }
505 >      }
506  
507 +      // all rigid bodies are integrated:
508 +
509 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
510 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
511 +        info[k].integrableObjects.push_back(mySD);      
512 +        molInfo.myIntegrableObjects.push_back(mySD);
513 +      }
514 +    
515 +      
516 +      // send the arrays off to the forceField for init.
517 +      
518        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
519        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
520        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
521        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
522                                   theTorsions);
523  
392 –
524        info[k].molecules[i].initialize(molInfo);
525  
526  
# Line 397 | Line 528 | void SimSetup::makeMolecules(void){
528        delete[] theBonds;
529        delete[] theBends;
530        delete[] theTorsions;
531 <    }
531 >    }    
532    }
533  
534   #ifdef IS_MPI
# Line 407 | Line 538 | void SimSetup::makeMolecules(void){
538  
539    // clean up the forcefield
540  
541 <  the_ff->calcRcut();
541 >  if (!globals->haveLJrcut()){
542 >
543 >    the_ff->calcRcut();
544 >
545 >  } else {
546 >    
547 >    the_ff->setRcut( globals->getLJrcut() );
548 >  }
549 >
550    the_ff->cleanMe();
551   }
552  
# Line 609 | Line 748 | void SimSetup::gatherInfo(void){
748    else if (!strcasecmp(force_field, "EAM")){
749      ffCase = FF_EAM;
750    }
751 +  else if (!strcasecmp(force_field, "WATER")){
752 +    ffCase = FF_H2O;
753 +  }
754    else{
755      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
756              force_field);
# Line 637 | Line 779 | void SimSetup::gatherInfo(void){
779    }
780    else{
781      sprintf(painCave.errMsg,
782 <            "SimSetup Warning. Unrecognized Ensemble -> %s, "
783 <            "reverting to NVE for this simulation.\n",
782 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
783 >            "\treverting to NVE for this simulation.\n",
784              ensemble);
785           painCave.isFatal = 0;
786           simError();
# Line 670 | Line 812 | void SimSetup::gatherInfo(void){
812        if (!the_components[i]->haveNMol()){
813          // we have a problem
814          sprintf(painCave.errMsg,
815 <                "SimSetup Error. No global NMol or component NMol"
816 <                " given. Cannot calculate the number of atoms.\n");
815 >                "SimSetup Error. No global NMol or component NMol given.\n"
816 >                "\tCannot calculate the number of atoms.\n");
817          painCave.isFatal = 1;
818          simError();
819        }
# Line 692 | Line 834 | void SimSetup::gatherInfo(void){
834    }
835  
836    //check whether sample time, status time, thermal time and reset time are divisble by dt
837 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
837 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
838      sprintf(painCave.errMsg,
839 <              "Sample time is not divisible by dt \n");
839 >            "Sample time is not divisible by dt.\n"
840 >            "\tThis will result in samples that are not uniformly\n"
841 >            "\tdistributed in time.  If this is a problem, change\n"
842 >            "\tyour sampleTime variable.\n");
843      painCave.isFatal = 0;
844      simError();    
845    }
846  
847 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
847 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
848      sprintf(painCave.errMsg,
849 <              "Status time is not divisible by dt\n");
849 >            "Status time is not divisible by dt.\n"
850 >            "\tThis will result in status reports that are not uniformly\n"
851 >            "\tdistributed in time.  If this is a problem, change \n"
852 >            "\tyour statusTime variable.\n");
853      painCave.isFatal = 0;
854      simError();    
855    }
856  
857    if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
858      sprintf(painCave.errMsg,
859 <              "Thermal time is not divisible by dt\n");
859 >            "Thermal time is not divisible by dt.\n"
860 >            "\tThis will result in thermalizations that are not uniformly\n"
861 >            "\tdistributed in time.  If this is a problem, change \n"
862 >            "\tyour thermalTime variable.\n");
863      painCave.isFatal = 0;
864      simError();    
865    }  
866  
867    if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
868      sprintf(painCave.errMsg,
869 <              "Reset time is not divisible by dt\n");
869 >            "Reset time is not divisible by dt.\n"
870 >            "\tThis will result in integrator resets that are not uniformly\n"
871 >            "\tdistributed in time.  If this is a problem, change\n"
872 >            "\tyour resetTime variable.\n");
873      painCave.isFatal = 0;
874      simError();    
875    }
# Line 726 | Line 880 | void SimSetup::gatherInfo(void){
880      if (globals->haveSampleTime()){
881        info[i].sampleTime = globals->getSampleTime();
882        info[i].statusTime = info[i].sampleTime;
729 –      info[i].thermalTime = info[i].sampleTime;
883      }
884      else{
885        info[i].sampleTime = globals->getRunTime();
886        info[i].statusTime = info[i].sampleTime;
734 –      info[i].thermalTime = info[i].sampleTime;
887      }
888  
889      if (globals->haveStatusTime()){
# Line 740 | Line 892 | void SimSetup::gatherInfo(void){
892  
893      if (globals->haveThermalTime()){
894        info[i].thermalTime = globals->getThermalTime();
895 +    } else {
896 +      info[i].thermalTime = globals->getRunTime();
897      }
898  
899      info[i].resetIntegrator = 0;
# Line 799 | Line 953 | void SimSetup::gatherInfo(void){
953    for (int i = 0; i < nInfo; i++){
954      info[i].setSeed(seedValue);
955    }
956 <
956 >  
957   #ifdef IS_MPI
958 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
958 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
959    MPIcheckPoint();
960   #endif // is_mpi
961   }
# Line 810 | Line 964 | void SimSetup::finalInfoCheck(void){
964   void SimSetup::finalInfoCheck(void){
965    int index;
966    int usesDipoles;
967 +  int usesCharges;
968    int i;
969  
970    for (i = 0; i < nInfo; i++){
# Line 821 | Line 976 | void SimSetup::finalInfoCheck(void){
976        usesDipoles = (info[i].atoms[index])->hasDipole();
977        index++;
978      }
979 <
979 >    index = 0;
980 >    usesCharges = 0;
981 >    while ((index < info[i].n_atoms) && !usesCharges){
982 >      usesCharges= (info[i].atoms[index])->hasCharge();
983 >      index++;
984 >    }
985   #ifdef IS_MPI
986      int myUse = usesDipoles;
987      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
# Line 834 | Line 994 | void SimSetup::finalInfoCheck(void){
994  
995        if (!globals->haveECR()){
996          sprintf(painCave.errMsg,
997 <                "SimSetup Warning: using default value of 15.0 angstroms"
998 <                "box length for the electrostaticCutoffRadius.\n");
997 >                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
998 >                "\tOOPSE will use a default value of 15.0 angstroms"
999 >                "\tfor the electrostaticCutoffRadius.\n");
1000          painCave.isFatal = 0;
1001          simError();
1002          theEcr = 15.0;
# Line 846 | Line 1007 | void SimSetup::finalInfoCheck(void){
1007  
1008        if (!globals->haveEST()){
1009          sprintf(painCave.errMsg,
1010 <                "SimSetup Warning: using default value of 0.05 * the "
1011 <                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
1010 >                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1011 >                "\tOOPSE will use a default value of\n"
1012 >                "\t0.05 * electrostaticCutoffRadius\n"
1013 >                "\tfor the electrostaticSkinThickness\n");
1014          painCave.isFatal = 0;
1015          simError();
1016          theEst = 0.05 * theEcr;
# Line 860 | Line 1023 | void SimSetup::finalInfoCheck(void){
1023  
1024        if (!globals->haveDielectric()){
1025          sprintf(painCave.errMsg,
1026 <                "SimSetup Error: You are trying to use Reaction Field without"
1027 <                "setting a dielectric constant!\n");
1026 >                "SimSetup Error: No Dielectric constant was set.\n"
1027 >                "\tYou are trying to use Reaction Field without"
1028 >                "\tsetting a dielectric constant!\n");
1029          painCave.isFatal = 1;
1030          simError();
1031        }
1032        info[i].dielectric = globals->getDielectric();
1033      }
1034      else{
1035 <      if (usesDipoles){
1035 >      if (usesDipoles || usesCharges){
1036          if (!globals->haveECR()){
1037            sprintf(painCave.errMsg,
1038 <                  "SimSetup Warning: using default value of 15.0 angstroms"
1039 <                  "box length for the electrostaticCutoffRadius.\n");
1038 >                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1039 >                  "\tOOPSE will use a default value of 15.0 angstroms"
1040 >                  "\tfor the electrostaticCutoffRadius.\n");
1041            painCave.isFatal = 0;
1042            simError();
1043            theEcr = 15.0;
# Line 883 | Line 1048 | void SimSetup::finalInfoCheck(void){
1048          
1049          if (!globals->haveEST()){
1050            sprintf(painCave.errMsg,
1051 <                  "SimSetup Warning: using default value of 0.05 * the "
1052 <                  "electrostaticCutoffRadius for the "
1053 <                  "electrostaticSkinThickness\n");
1051 >                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1052 >                  "\tOOPSE will use a default value of\n"
1053 >                  "\t0.05 * electrostaticCutoffRadius\n"
1054 >                  "\tfor the electrostaticSkinThickness\n");
1055            painCave.isFatal = 0;
1056            simError();
1057            theEst = 0.05 * theEcr;
# Line 1090 | Line 1256 | void SimSetup::createFF(void){
1256  
1257      case FF_EAM:
1258        the_ff = new EAM_FF();
1259 +      break;
1260 +
1261 +    case FF_H2O:
1262 +      the_ff = new WATER();
1263        break;
1264  
1265      default:
# Line 1164 | Line 1334 | void SimSetup::calcSysValues(void){
1334    tot_bonds = 0;
1335    tot_bends = 0;
1336    tot_torsions = 0;
1337 +  tot_rigid = 0;
1338    for (i = 0; i < n_components; i++){
1339      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1340      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1341      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1342      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1343 +    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1344    }
1345 <
1345 >  
1346    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1347    molMembershipArray = new int[tot_atoms];
1348  
# Line 1192 | Line 1364 | void SimSetup::mpiMolDivide(void){
1364    int i, j, k;
1365    int localMol, allMol;
1366    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1367 +  int local_rigid;
1368 +  vector<int> globalMolIndex;
1369  
1370    mpiSim = new mpiSimulation(info);
1371  
1372 <  globalIndex = mpiSim->divideLabor();
1372 >  mpiSim->divideLabor();
1373 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1374 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1375  
1376    // set up the local variables
1377  
# Line 1208 | Line 1384 | void SimSetup::mpiMolDivide(void){
1384    local_bonds = 0;
1385    local_bends = 0;
1386    local_torsions = 0;
1387 <  globalAtomIndex = 0;
1387 >  local_rigid = 0;
1388 >  globalAtomCounter = 0;
1389  
1213 –
1390    for (i = 0; i < n_components; i++){
1391      for (j = 0; j < components_nmol[i]; j++){
1392        if (mol2proc[allMol] == worldRank){
# Line 1218 | Line 1394 | void SimSetup::mpiMolDivide(void){
1394          local_bonds += comp_stamps[i]->getNBonds();
1395          local_bends += comp_stamps[i]->getNBends();
1396          local_torsions += comp_stamps[i]->getNTorsions();
1397 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1398          localMol++;
1399        }      
1400        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1401 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1402 <        globalAtomIndex++;
1401 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1402 >        globalAtomCounter++;
1403        }
1404  
1405        allMol++;
# Line 1231 | Line 1408 | void SimSetup::mpiMolDivide(void){
1408    local_SRI = local_bonds + local_bends + local_torsions;
1409  
1410    info[0].n_atoms = mpiSim->getMyNlocal();  
1411 +  
1412  
1413    if (local_atoms != info[0].n_atoms){
1414      sprintf(painCave.errMsg,
1415 <            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1416 <            " localAtom (%d) are not equal.\n",
1415 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1416 >            "\tlocalAtom (%d) are not equal.\n",
1417              info[0].n_atoms, local_atoms);
1418      painCave.isFatal = 1;
1419      simError();
# Line 1263 | Line 1441 | void SimSetup::makeSysArrays(void){
1441  
1442    Atom** the_atoms;
1443    Molecule* the_molecules;
1266 –  Exclude** the_excludes;
1444  
1268 –
1445    for (l = 0; l < nInfo; l++){
1446      // create the atom and short range interaction arrays
1447  
# Line 1291 | Line 1467 | void SimSetup::makeSysArrays(void){
1467   #else // is_mpi
1468  
1469      molIndex = 0;
1470 <    globalAtomIndex = 0;
1470 >    globalAtomCounter = 0;
1471      for (i = 0; i < n_components; i++){
1472        for (j = 0; j < components_nmol[i]; j++){
1473          the_molecules[molIndex].setStampID(i);
1474          the_molecules[molIndex].setMyIndex(molIndex);
1475          the_molecules[molIndex].setGlobalIndex(molIndex);
1476          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1477 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1478 <          globalAtomIndex++;
1477 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1478 >          globalAtomCounter++;
1479          }
1480          molIndex++;
1481        }
# Line 1308 | Line 1484 | void SimSetup::makeSysArrays(void){
1484  
1485   #endif // is_mpi
1486  
1487 <
1488 <    if (info[l].n_SRI){
1489 <      Exclude::createArray(info[l].n_SRI);
1314 <      the_excludes = new Exclude * [info[l].n_SRI];
1315 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1316 <        the_excludes[ex] = new Exclude(ex);
1317 <      }
1318 <      info[l].globalExcludes = new int;
1319 <      info[l].n_exclude = info[l].n_SRI;
1320 <    }
1321 <    else{
1322 <      Exclude::createArray(1);
1323 <      the_excludes = new Exclude * ;
1324 <      the_excludes[0] = new Exclude(0);
1325 <      the_excludes[0]->setPair(0, 0);
1326 <      info[l].globalExcludes = new int;
1327 <      info[l].globalExcludes[0] = 0;
1328 <      info[l].n_exclude = 0;
1329 <    }
1330 <
1487 >    info[l].globalExcludes = new int;
1488 >    info[l].globalExcludes[0] = 0;
1489 >    
1490      // set the arrays into the SimInfo object
1491  
1492      info[l].atoms = the_atoms;
1493      info[l].molecules = the_molecules;
1494      info[l].nGlobalExcludes = 0;
1336 –    info[l].excludes = the_excludes;
1495  
1496      the_ff->setSimInfo(info);
1497    }
# Line 1377 | Line 1535 | void SimSetup::makeIntegrator(void){
1535          else{
1536            sprintf(painCave.errMsg,
1537                    "SimSetup error: If you use the NVT\n"
1538 <                  "    ensemble, you must set tauThermostat.\n");
1538 >                  "\tensemble, you must set tauThermostat.\n");
1539            painCave.isFatal = 1;
1540            simError();
1541          }
# Line 1400 | Line 1558 | void SimSetup::makeIntegrator(void){
1558          else{
1559            sprintf(painCave.errMsg,
1560                    "SimSetup error: If you use a constant pressure\n"
1561 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1561 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1562            painCave.isFatal = 1;
1563            simError();
1564          }
# Line 1410 | Line 1568 | void SimSetup::makeIntegrator(void){
1568          else{
1569            sprintf(painCave.errMsg,
1570                    "SimSetup error: If you use an NPT\n"
1571 <                  "    ensemble, you must set tauThermostat.\n");
1571 >                  "\tensemble, you must set tauThermostat.\n");
1572            painCave.isFatal = 1;
1573            simError();
1574          }
# Line 1420 | 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 tauBarostat.\n");
1581 >                  "\tensemble, you must set tauBarostat.\n");
1582            painCave.isFatal = 1;
1583            simError();
1584          }
# Line 1443 | Line 1601 | void SimSetup::makeIntegrator(void){
1601          else{
1602            sprintf(painCave.errMsg,
1603                    "SimSetup error: If you use a constant pressure\n"
1604 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1604 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1605            painCave.isFatal = 1;
1606            simError();
1607          }    
# Line 1454 | Line 1612 | void SimSetup::makeIntegrator(void){
1612          else{
1613            sprintf(painCave.errMsg,
1614                    "SimSetup error: If you use an NPT\n"
1615 <                  "    ensemble, you must set tauThermostat.\n");
1615 >                  "\tensemble, you must set tauThermostat.\n");
1616            painCave.isFatal = 1;
1617            simError();
1618          }
# Line 1465 | Line 1623 | void SimSetup::makeIntegrator(void){
1623          else{
1624            sprintf(painCave.errMsg,
1625                    "SimSetup error: If you use an NPT\n"
1626 <                  "    ensemble, you must set tauBarostat.\n");
1626 >                  "\tensemble, you must set tauBarostat.\n");
1627            painCave.isFatal = 1;
1628            simError();
1629          }
# Line 1488 | Line 1646 | void SimSetup::makeIntegrator(void){
1646          else{
1647            sprintf(painCave.errMsg,
1648                    "SimSetup error: If you use a constant pressure\n"
1649 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1649 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1650            painCave.isFatal = 1;
1651            simError();
1652          }    
# Line 1498 | Line 1656 | void SimSetup::makeIntegrator(void){
1656          else{
1657            sprintf(painCave.errMsg,
1658                    "SimSetup error: If you use an NPT\n"
1659 <                  "    ensemble, you must set tauThermostat.\n");
1659 >                  "\tensemble, you must set tauThermostat.\n");
1660            painCave.isFatal = 1;
1661            simError();
1662          }
# Line 1508 | Line 1666 | void SimSetup::makeIntegrator(void){
1666          else{
1667            sprintf(painCave.errMsg,
1668                    "SimSetup error: If you use an NPT\n"
1669 <                  "    ensemble, you must set tauBarostat.\n");
1669 >                  "\tensemble, you must set tauBarostat.\n");
1670            painCave.isFatal = 1;
1671            simError();
1672          }
# Line 1561 | Line 1719 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1719    }
1720    else{
1721      sprintf(painCave.errMsg,
1722 <            "ZConstraint error: If you use an ZConstraint\n"
1723 <            " , you must set sample time.\n");
1722 >            "ZConstraint error: If you use a ZConstraint,\n"
1723 >            "\tyou must set zconsTime.\n");
1724      painCave.isFatal = 1;
1725      simError();
1726    }
# Line 1577 | Line 1735 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1735    else{
1736      double defaultZConsTol = 0.01;
1737      sprintf(painCave.errMsg,
1738 <            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1739 <            " , default value %f is used.\n",
1738 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
1739 >            "\tOOPSE will use a default value of %f.\n"
1740 >            "\tTo set the tolerance, use the zconsTol variable.\n",
1741              defaultZConsTol);
1742      painCave.isFatal = 0;
1743      simError();      
# Line 1596 | Line 1755 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1755    }
1756    else{
1757      sprintf(painCave.errMsg,
1758 <            "ZConstraint Warning: User does not set force Subtraction policy, "
1759 <            "PolicyByMass is used\n");
1758 >            "ZConstraint Warning: No force subtraction policy was set.\n"
1759 >            "\tOOPSE will use PolicyByMass.\n"
1760 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
1761      painCave.isFatal = 0;
1762      simError();
1763      zconsForcePolicy->setData("BYMASS");
1764    }
1765  
1766    theInfo.addProperty(zconsForcePolicy);
1767 +
1768 +  //set zcons gap
1769 +  DoubleData* zconsGap = new DoubleData();
1770 +  zconsGap->setID(ZCONSGAP_ID);
1771 +
1772 +  if (globals->haveZConsGap()){
1773 +    zconsGap->setData(globals->getZconsGap());
1774 +    theInfo.addProperty(zconsGap);  
1775 +  }
1776 +
1777 +  //set zcons fixtime
1778 +  DoubleData* zconsFixtime = new DoubleData();
1779 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1780 +
1781 +  if (globals->haveZConsFixTime()){
1782 +    zconsFixtime->setData(globals->getZconsFixtime());
1783 +    theInfo.addProperty(zconsFixtime);  
1784 +  }
1785 +
1786 +  //set zconsUsingSMD
1787 +  IntData* zconsUsingSMD = new IntData();
1788 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1789 +
1790 +  if (globals->haveZConsUsingSMD()){
1791 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1792 +    theInfo.addProperty(zconsUsingSMD);  
1793 +  }
1794  
1795    //Determine the name of ouput file and add it into SimInfo's property list
1796    //Be careful, do not use inFileName, since it is a pointer which
# Line 1634 | Line 1821 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1821      tempParaItem.zPos = zconStamp[i]->getZpos();
1822      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1823      tempParaItem.kRatio = zconStamp[i]->getKratio();
1824 <
1824 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1825 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1826      zconsParaData->addItem(tempParaItem);
1827    }
1828  
1829    //check the uniqueness of index  
1830    if(!zconsParaData->isIndexUnique()){
1831      sprintf(painCave.errMsg,
1832 <            "ZConstraint Error: molIndex is not unique\n");
1832 >            "ZConstraint Error: molIndex is not unique!\n");
1833      painCave.isFatal = 1;
1834      simError();
1835    }
# Line 1652 | Line 1840 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1840    //push data into siminfo, therefore, we can retrieve later
1841    theInfo.addProperty(zconsParaData);
1842   }
1843 +
1844 + void SimSetup::makeMinimizer(){
1845 +
1846 +  OOPSEMinimizer* myOOPSEMinimizer;
1847 +  MinimizerParameterSet* param;
1848 +  char minimizerName[100];
1849 +  
1850 +  for (int i = 0; i < nInfo; i++){
1851 +    
1852 +    //prepare parameter set for minimizer
1853 +    param = new MinimizerParameterSet();
1854 +    param->setDefaultParameter();
1855 +
1856 +    if (globals->haveMinimizer()){
1857 +      param->setFTol(globals->getMinFTol());
1858 +    }
1859 +
1860 +    if (globals->haveMinGTol()){
1861 +      param->setGTol(globals->getMinGTol());
1862 +    }
1863 +
1864 +    if (globals->haveMinMaxIter()){
1865 +      param->setMaxIteration(globals->getMinMaxIter());
1866 +    }
1867 +
1868 +    if (globals->haveMinWriteFrq()){
1869 +      param->setMaxIteration(globals->getMinMaxIter());
1870 +    }
1871 +
1872 +    if (globals->haveMinWriteFrq()){
1873 +      param->setWriteFrq(globals->getMinWriteFrq());
1874 +    }
1875 +    
1876 +    if (globals->haveMinStepSize()){
1877 +      param->setStepSize(globals->getMinStepSize());
1878 +    }
1879 +
1880 +    if (globals->haveMinLSMaxIter()){
1881 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
1882 +    }    
1883 +
1884 +    if (globals->haveMinLSTol()){
1885 +      param->setLineSearchTol(globals->getMinLSTol());
1886 +    }    
1887 +
1888 +    strcpy(minimizerName, globals->getMinimizer());
1889 +
1890 +    if (!strcasecmp(minimizerName, "CG")){
1891 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1892 +    }
1893 +    else if (!strcasecmp(minimizerName, "SD")){
1894 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1895 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1896 +    }
1897 +    else{
1898 +          sprintf(painCave.errMsg,
1899 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1900 +          painCave.isFatal = 0;
1901 +          simError();
1902 +
1903 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1904 +    }
1905 +     info[i].the_integrator = myOOPSEMinimizer;
1906 +
1907 +     //store the minimizer into simInfo
1908 +     info[i].the_minimizer = myOOPSEMinimizer;
1909 +     info[i].has_minimizer = true;
1910 +  }
1911 +
1912 + }

Diff Legend

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