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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 983 by gezelter, Mon Jan 26 21:45:03 2004 UTC vs.
Revision 1104 by gezelter, Tue Apr 13 16:26:03 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 144 | 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 159 | 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 +
194    //init the forceField paramters
195  
196    the_ff->readParams();
197  
181 –
198    // init the atoms
199  
200 <  double phi, theta, psi;
185 <  double sux, suy, suz;
186 <  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
187 <  double ux, uy, uz, u, uSqr;
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 198 | 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]);
204 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
205 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
206 –      molInfo.myBends = new Bend * [molInfo.nBends];
207 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
217  
218 +      if (molInfo.nBonds > 0)
219 +        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
220 +      else
221 +        molInfo.myBonds = NULL;
222 +
223 +      if (molInfo.nBends > 0)
224 +        molInfo.myBends = new (Bend *) [molInfo.nBends];
225 +      else
226 +        molInfo.myBends = NULL;
227 +
228 +      if (molInfo.nTorsions > 0)
229 +        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
230 +      else
231 +        molInfo.myTorsions = NULL;
232 +
233        theBonds = new bond_pair[molInfo.nBonds];
234        theBends = new bend_set[molInfo.nBends];
235        theTorsions = new torsion_set[molInfo.nTorsions];
236 <
236 >      
237        // make the Atoms
238  
239        for (j = 0; j < molInfo.nAtoms; j++){
240          currentAtom = comp_stamps[stampID]->getAtom(j);
241 +
242          if (currentAtom->haveOrientation()){
243            dAtom = new DirectionalAtom((j + atomOffset),
244                                        info[k].getConfiguration());
# Line 224 | Line 249 | void SimSetup::makeMolecules(void){
249            // in space using the three Euler angles.  We assume the standard
250            // unit vector was originally along the z axis below.
251  
252 <          phi = currentAtom->getEulerPhi();
253 <          theta = currentAtom->getEulerTheta();
254 <          psi = currentAtom->getEulerPsi();
230 <            
231 <          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
232 <          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
233 <          Axz = sin(theta) * sin(psi);
234 <          
235 <          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
236 <          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
237 <          Ayz = sin(theta) * cos(psi);
238 <          
239 <          Azx = sin(phi) * sin(theta);
240 <          Azy = -cos(phi) * sin(theta);
241 <          Azz = cos(theta);
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 <          sux = 0.0;
257 <          suy = 0.0;
245 <          suz = 1.0;
246 <
247 <          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
248 <          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
249 <          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
250 <
251 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
252 <
253 <          u = sqrt(uSqr);
254 <          ux = ux / u;
255 <          uy = uy / u;
256 <          uz = uz / u;
257 <
258 <          dAtom->setSUx(ux);
259 <          dAtom->setSUy(uy);
260 <          dAtom->setSUz(uz);
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 278 | 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  
284 –        // exclude_I must always be the smaller of the pair
285 –        if (exI > exJ){
286 –          tempEx = exI;
287 –          exI = exJ;
288 –          exJ = tempEx;
289 –        }
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;
287 <
288 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
297 < #else  // isn't MPI
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 + 1), (exJ + 1));
300 < #endif  //is_mpi
290 >        info[k].excludes->addPair(exI, exJ);
291        }
302 –      excludeOffset += molInfo.nBonds;
292  
293        //make the bends
294        for (j = 0; j < molInfo.nBends; j++){
# Line 349 | 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{
357 <          exI = theBends[j].a;
358 <          exJ = theBends[j].b;
359 <        }
360 <
361 <        // exclude_I must always be the smaller of the pair
362 <        if (exI > exJ){
363 <          tempEx = exI;
364 <          exI = exJ;
365 <          exJ = tempEx;
366 <        }
341 >        if (theBends[j].isGhost) {
342 >          
343 >          tempI = theBends[j].a;
344 >          tempJ = theBends[j].b;
345 >          
346   #ifdef IS_MPI
347 <        tempEx = exI;
348 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
349 <        tempEx = exJ;
350 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
347 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
348 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
349 > #else
350 >          exI = tempI + 1;
351 >          exJ = tempJ + 1;
352 > #endif          
353 >          info[k].excludes->addPair(exI, exJ);
354  
355 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
356 < #else  // isn't MPI
357 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
358 < #endif  //is_mpi
355 >        } else {
356 >
357 >          tempI = theBends[j].a;
358 >          tempJ = theBends[j].b;
359 >          tempK = theBends[j].c;
360 >          
361 > #ifdef IS_MPI
362 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
363 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
364 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
365 > #else
366 >          exI = tempI + 1;
367 >          exJ = tempJ + 1;
368 >          exK = tempK + 1;
369 > #endif
370 >          
371 >          info[k].excludes->addPair(exI, exK);
372 >          info[k].excludes->addPair(exI, exJ);
373 >          info[k].excludes->addPair(exJ, exK);
374 >        }
375        }
378 –      excludeOffset += molInfo.nBends;
376  
377        for (j = 0; j < molInfo.nTorsions; j++){
378          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 384 | 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  
390 –        // exclude_I must always be the smaller of the pair
391 –        if (exI > exJ){
392 –          tempEx = exI;
393 –          exI = exJ;
394 –          exJ = tempEx;
395 –        }
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        }
407 –      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 >        molInfo.myRigidBodies.push_back(myRB);
470 >        info[k].rigidBodies.push_back(myRB);
471 >      }
472 >      
473 >
474 >      // After this is all set up, scan through the atoms to
475 >      // see if they can be added to the integrableObjects:
476 >
477 >      for (j = 0; j < molInfo.nAtoms; j++){
478 >
479 > #ifdef IS_MPI
480 >        slJ = molInfo.myAtoms[j]->getGlobalIndex();
481 > #else
482 >        slJ = j+atomOffset;
483 > #endif
484  
485 +        // if they aren't on the skip list, then they can be integrated
486 +
487 +        if (skipList.find(slJ) == skipList.end()) {
488 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
489 +          info[k].integrableObjects.push_back(mySD);
490 +          molInfo.myIntegrableObjects.push_back(mySD);
491 +        }
492 +      }
493 +
494 +      // all rigid bodies are integrated:
495 +
496 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
497 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
498 +        info[k].integrableObjects.push_back(mySD);      
499 +        molInfo.myIntegrableObjects.push_back(mySD);
500 +      }
501 +    
502 +      
503 +      // send the arrays off to the forceField for init.
504 +      
505        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
506        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
507        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
508        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
509                                   theTorsions);
510  
418 –
511        info[k].molecules[i].initialize(molInfo);
512  
513  
# Line 423 | Line 515 | void SimSetup::makeMolecules(void){
515        delete[] theBonds;
516        delete[] theBends;
517        delete[] theTorsions;
518 <    }
518 >    }    
519    }
520  
521   #ifdef IS_MPI
# Line 433 | Line 525 | void SimSetup::makeMolecules(void){
525  
526    // clean up the forcefield
527  
528 <  the_ff->calcRcut();
528 >  if (!globals->haveLJrcut()){
529 >
530 >    the_ff->calcRcut();
531 >
532 >  } else {
533 >    
534 >    the_ff->setRcut( globals->getLJrcut() );
535 >  }
536 >
537    the_ff->cleanMe();
538   }
539  
# Line 635 | Line 735 | void SimSetup::gatherInfo(void){
735    else if (!strcasecmp(force_field, "EAM")){
736      ffCase = FF_EAM;
737    }
738 +  else if (!strcasecmp(force_field, "WATER")){
739 +    ffCase = FF_H2O;
740 +  }
741    else{
742      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
743              force_field);
# Line 837 | Line 940 | void SimSetup::gatherInfo(void){
940    for (int i = 0; i < nInfo; i++){
941      info[i].setSeed(seedValue);
942    }
943 <
943 >  
944   #ifdef IS_MPI
945 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
945 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
946    MPIcheckPoint();
947   #endif // is_mpi
948   }
# Line 1136 | Line 1239 | void SimSetup::createFF(void){
1239        the_ff = new EAM_FF();
1240        break;
1241  
1242 +    case FF_H2O:
1243 +      the_ff = new WATER();
1244 +      break;
1245 +
1246      default:
1247        sprintf(painCave.errMsg,
1248                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1208 | Line 1315 | void SimSetup::calcSysValues(void){
1315    tot_bonds = 0;
1316    tot_bends = 0;
1317    tot_torsions = 0;
1318 +  tot_rigid = 0;
1319    for (i = 0; i < n_components; i++){
1320      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1321      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1322      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1323      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1324 +    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1325    }
1326 <
1326 >  
1327    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1328    molMembershipArray = new int[tot_atoms];
1329  
# Line 1236 | Line 1345 | void SimSetup::mpiMolDivide(void){
1345    int i, j, k;
1346    int localMol, allMol;
1347    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1348 +  int local_rigid;
1349  
1350    mpiSim = new mpiSimulation(info);
1351  
# Line 1252 | Line 1362 | void SimSetup::mpiMolDivide(void){
1362    local_bonds = 0;
1363    local_bends = 0;
1364    local_torsions = 0;
1365 +  local_rigid = 0;
1366    globalAtomIndex = 0;
1367  
1257 –
1368    for (i = 0; i < n_components; i++){
1369      for (j = 0; j < components_nmol[i]; j++){
1370        if (mol2proc[allMol] == worldRank){
# Line 1262 | Line 1372 | void SimSetup::mpiMolDivide(void){
1372          local_bonds += comp_stamps[i]->getNBonds();
1373          local_bends += comp_stamps[i]->getNBends();
1374          local_torsions += comp_stamps[i]->getNTorsions();
1375 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1376          localMol++;
1377        }      
1378        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
# Line 1275 | Line 1386 | void SimSetup::mpiMolDivide(void){
1386    local_SRI = local_bonds + local_bends + local_torsions;
1387  
1388    info[0].n_atoms = mpiSim->getMyNlocal();  
1389 +  
1390  
1391    if (local_atoms != info[0].n_atoms){
1392      sprintf(painCave.errMsg,
# Line 1307 | Line 1419 | void SimSetup::makeSysArrays(void){
1419  
1420    Atom** the_atoms;
1421    Molecule* the_molecules;
1310 –  Exclude** the_excludes;
1422  
1312 –
1423    for (l = 0; l < nInfo; l++){
1424      // create the atom and short range interaction arrays
1425  
# Line 1352 | Line 1462 | void SimSetup::makeSysArrays(void){
1462  
1463   #endif // is_mpi
1464  
1465 <
1466 <    if (info[l].n_SRI){
1467 <      Exclude::createArray(info[l].n_SRI);
1358 <      the_excludes = new Exclude * [info[l].n_SRI];
1359 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1360 <        the_excludes[ex] = new Exclude(ex);
1361 <      }
1362 <      info[l].globalExcludes = new int;
1363 <      info[l].n_exclude = info[l].n_SRI;
1364 <    }
1365 <    else{
1366 <      Exclude::createArray(1);
1367 <      the_excludes = new Exclude * ;
1368 <      the_excludes[0] = new Exclude(0);
1369 <      the_excludes[0]->setPair(0, 0);
1370 <      info[l].globalExcludes = new int;
1371 <      info[l].globalExcludes[0] = 0;
1372 <      info[l].n_exclude = 0;
1373 <    }
1374 <
1465 >    info[l].globalExcludes = new int;
1466 >    info[l].globalExcludes[0] = 0;
1467 >    
1468      // set the arrays into the SimInfo object
1469  
1470      info[l].atoms = the_atoms;
1471      info[l].molecules = the_molecules;
1472      info[l].nGlobalExcludes = 0;
1380 –    info[l].excludes = the_excludes;
1473  
1474      the_ff->setSimInfo(info);
1475    }
# Line 1650 | Line 1742 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1742    }
1743  
1744    theInfo.addProperty(zconsForcePolicy);
1745 +
1746 +  //set zcons gap
1747 +  DoubleData* zconsGap = new DoubleData();
1748 +  zconsGap->setID(ZCONSGAP_ID);
1749 +
1750 +  if (globals->haveZConsGap()){
1751 +    zconsGap->setData(globals->getZconsGap());
1752 +    theInfo.addProperty(zconsGap);  
1753 +  }
1754 +
1755 +  //set zcons fixtime
1756 +  DoubleData* zconsFixtime = new DoubleData();
1757 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1758  
1759 +  if (globals->haveZConsFixTime()){
1760 +    zconsFixtime->setData(globals->getZconsFixtime());
1761 +    theInfo.addProperty(zconsFixtime);  
1762 +  }
1763 +
1764 +  //set zconsUsingSMD
1765 +  IntData* zconsUsingSMD = new IntData();
1766 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1767 +
1768 +  if (globals->haveZConsUsingSMD()){
1769 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1770 +    theInfo.addProperty(zconsUsingSMD);  
1771 +  }
1772 +
1773    //Determine the name of ouput file and add it into SimInfo's property list
1774    //Be careful, do not use inFileName, since it is a pointer which
1775    //point to a string at master node, and slave nodes do not contain that string
# Line 1680 | Line 1799 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1799      tempParaItem.zPos = zconStamp[i]->getZpos();
1800      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1801      tempParaItem.kRatio = zconStamp[i]->getKratio();
1802 <
1802 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1803 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1804      zconsParaData->addItem(tempParaItem);
1805    }
1806  
# Line 1697 | Line 1817 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1817    
1818    //push data into siminfo, therefore, we can retrieve later
1819    theInfo.addProperty(zconsParaData);
1820 + }
1821 +
1822 + void SimSetup::makeMinimizer(){
1823 +
1824 +  OOPSEMinimizer* myOOPSEMinimizer;
1825 +  MinimizerParameterSet* param;
1826 +  char minimizerName[100];
1827 +  
1828 +  for (int i = 0; i < nInfo; i++){
1829 +    
1830 +    //prepare parameter set for minimizer
1831 +    param = new MinimizerParameterSet();
1832 +    param->setDefaultParameter();
1833 +
1834 +    if (globals->haveMinimizer()){
1835 +      param->setFTol(globals->getMinFTol());
1836 +    }
1837 +
1838 +    if (globals->haveMinGTol()){
1839 +      param->setGTol(globals->getMinGTol());
1840 +    }
1841 +
1842 +    if (globals->haveMinMaxIter()){
1843 +      param->setMaxIteration(globals->getMinMaxIter());
1844 +    }
1845 +
1846 +    if (globals->haveMinWriteFrq()){
1847 +      param->setMaxIteration(globals->getMinMaxIter());
1848 +    }
1849 +
1850 +    if (globals->haveMinWriteFrq()){
1851 +      param->setWriteFrq(globals->getMinWriteFrq());
1852 +    }
1853 +    
1854 +    if (globals->haveMinStepSize()){
1855 +      param->setStepSize(globals->getMinStepSize());
1856 +    }
1857 +
1858 +    if (globals->haveMinLSMaxIter()){
1859 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
1860 +    }    
1861 +
1862 +    if (globals->haveMinLSTol()){
1863 +      param->setLineSearchTol(globals->getMinLSTol());
1864 +    }    
1865 +
1866 +    strcpy(minimizerName, globals->getMinimizer());
1867 +
1868 +    if (!strcasecmp(minimizerName, "CG")){
1869 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1870 +    }
1871 +    else if (!strcasecmp(minimizerName, "SD")){
1872 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1873 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1874 +    }
1875 +    else{
1876 +          sprintf(painCave.errMsg,
1877 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1878 +          painCave.isFatal = 0;
1879 +          simError();
1880 +
1881 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1882 +    }
1883 +     info[i].the_integrator = myOOPSEMinimizer;
1884 +
1885 +     //store the minimizer into simInfo
1886 +     info[i].the_minimizer = myOOPSEMinimizer;
1887 +     info[i].has_minimizer = true;
1888 +  }
1889 +
1890   }

Diff Legend

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