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 1031 by tim, Fri Feb 6 18:58:06 2004 UTC vs.
Revision 1108 by tim, Wed Apr 14 15:37:41 2004 UTC

# Line 9 | Line 9
9   #include "parse_me.h"
10   #include "Integrator.hpp"
11   #include "simError.h"
12 < #include "ConjugateMinimizer.hpp"
12 > #include "RigidBody.hpp"
13 > //#include "ConjugateMinimizer.hpp"
14 > #include "OOPSEMinimizer.hpp"
15  
16   #ifdef IS_MPI
17   #include "mpiBASS.h"
# Line 28 | Line 30
30   #define FF_DUFF  0
31   #define FF_LJ    1
32   #define FF_EAM   2
33 < #define FF_H2O 3
33 > #define FF_H2O   3
34  
35   using namespace std;
36  
# Line 164 | Line 166 | void SimSetup::makeMolecules(void){
166  
167  
168   void SimSetup::makeMolecules(void){
169 <  int k;
170 <  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
169 >  int i, j, k;
170 >  int exI, exJ, exK, exL, slI, slJ;
171 >  int tempI, tempJ, tempK, tempL;
172 >  int molI;
173 >  int stampID, atomOffset, rbOffset;
174    molInit molInfo;
175    DirectionalAtom* dAtom;
176 +  RigidBody* myRB;
177 +  StuntDouble* mySD;
178    LinkedAssign* extras;
179    LinkedAssign* current_extra;
180    AtomStamp* currentAtom;
181    BondStamp* currentBond;
182    BendStamp* currentBend;
183    TorsionStamp* currentTorsion;
184 +  RigidBodyStamp* currentRigidBody;
185  
186    bond_pair* theBonds;
187    bend_set* theBends;
188    torsion_set* theTorsions;
189  
190 +  set<int> skipList;
191 +
192 +  double phi, theta, psi;
193 +  char* molName;
194 +  char rbName[100];
195 +
196    //init the forceField paramters
197  
198    the_ff->readParams();
199  
186 –
200    // init the atoms
201  
202 <  double phi, theta, psi;
190 <  double sux, suy, suz;
191 <  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
192 <  double ux, uy, uz, u, uSqr;
202 >  int nMembers, nNew, rb1, rb2;
203  
204    for (k = 0; k < nInfo; k++){
205      the_ff->setSimInfo(&(info[k]));
206  
207      atomOffset = 0;
208 <    excludeOffset = 0;
208 >
209      for (i = 0; i < info[k].n_mol; i++){
210        stampID = info[k].molecules[i].getStampID();
211 +      molName = comp_stamps[stampID]->getID();
212  
213        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
214        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
215        molInfo.nBends = comp_stamps[stampID]->getNBends();
216        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
217 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
218 <
217 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
218 >      
219        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
209 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
210 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
211 –      molInfo.myBends = new Bend * [molInfo.nBends];
212 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
220  
221 +      if (molInfo.nBonds > 0)
222 +        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
223 +      else
224 +        molInfo.myBonds = NULL;
225 +
226 +      if (molInfo.nBends > 0)
227 +        molInfo.myBends = new (Bend *) [molInfo.nBends];
228 +      else
229 +        molInfo.myBends = NULL;
230 +
231 +      if (molInfo.nTorsions > 0)
232 +        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
233 +      else
234 +        molInfo.myTorsions = NULL;
235 +
236        theBonds = new bond_pair[molInfo.nBonds];
237        theBends = new bend_set[molInfo.nBends];
238        theTorsions = new torsion_set[molInfo.nTorsions];
239 <
239 >      
240        // make the Atoms
241  
242        for (j = 0; j < molInfo.nAtoms; j++){
243          currentAtom = comp_stamps[stampID]->getAtom(j);
244 +
245          if (currentAtom->haveOrientation()){
246            dAtom = new DirectionalAtom((j + atomOffset),
247                                        info[k].getConfiguration());
# Line 232 | Line 255 | void SimSetup::makeMolecules(void){
255            phi = currentAtom->getEulerPhi() * M_PI / 180.0;
256            theta = currentAtom->getEulerTheta() * M_PI / 180.0;
257            psi = currentAtom->getEulerPsi()* M_PI / 180.0;
258 +
259 +          dAtom->setUnitFrameFromEuler(phi, theta, psi);
260              
261 <          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
262 <          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
238 <          Axz = sin(theta) * sin(psi);
239 <          
240 <          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
241 <          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
242 <          Ayz = sin(theta) * cos(psi);
243 <          
244 <          Azx = sin(phi) * sin(theta);
245 <          Azy = -cos(phi) * sin(theta);
246 <          Azz = cos(theta);
261 >        }
262 >        else{
263  
264 <          sux = 0.0;
249 <          suy = 0.0;
250 <          suz = 1.0;
251 <
252 <          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
253 <          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
254 <          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
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;
262 <
263 <          dAtom->setSUx(ux);
264 <          dAtom->setSUy(uy);
265 <          dAtom->setSUz(uz);
264 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
265          }
266 <        else{
268 <          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
269 <                                               info[k].getConfiguration());
270 <        }
266 >
267          molInfo.myAtoms[j]->setType(currentAtom->getType());
268  
269   #ifdef IS_MPI
270  
271 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
271 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
272  
273   #endif // is_mpi
274        }
# Line 283 | Line 279 | void SimSetup::makeMolecules(void){
279          theBonds[j].a = currentBond->getA() + atomOffset;
280          theBonds[j].b = currentBond->getB() + atomOffset;
281  
282 <        exI = theBonds[j].a;
283 <        exJ = theBonds[j].b;
282 >        tempI = theBonds[j].a;
283 >        tempJ = theBonds[j].b;
284  
289 –        // exclude_I must always be the smaller of the pair
290 –        if (exI > exJ){
291 –          tempEx = exI;
292 –          exI = exJ;
293 –          exJ = tempEx;
294 –        }
285   #ifdef IS_MPI
286 <        tempEx = exI;
287 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
288 <        tempEx = exJ;
289 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
286 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
287 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
288 > #else
289 >        exI = tempI + 1;
290 >        exJ = tempJ + 1;
291 > #endif
292  
293 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
302 < #else  // isn't MPI
303 <
304 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
305 < #endif  //is_mpi
293 >        info[k].excludes->addPair(exI, exJ);
294        }
307 –      excludeOffset += molInfo.nBonds;
295  
296        //make the bends
297        for (j = 0; j < molInfo.nBends; j++){
# Line 354 | Line 341 | void SimSetup::makeMolecules(void){
341            }
342          }
343  
344 <        if (!theBends[j].isGhost){
345 <          exI = theBends[j].a;
346 <          exJ = theBends[j].c;
347 <        }
348 <        else{
362 <          exI = theBends[j].a;
363 <          exJ = theBends[j].b;
364 <        }
365 <
366 <        // exclude_I must always be the smaller of the pair
367 <        if (exI > exJ){
368 <          tempEx = exI;
369 <          exI = exJ;
370 <          exJ = tempEx;
371 <        }
344 >        if (theBends[j].isGhost) {
345 >          
346 >          tempI = theBends[j].a;
347 >          tempJ = theBends[j].b;
348 >          
349   #ifdef IS_MPI
350 <        tempEx = exI;
351 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
352 <        tempEx = exJ;
353 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
350 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
351 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
352 > #else
353 >          exI = tempI + 1;
354 >          exJ = tempJ + 1;
355 > #endif          
356 >          info[k].excludes->addPair(exI, exJ);
357  
358 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
359 < #else  // isn't MPI
360 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
361 < #endif  //is_mpi
358 >        } else {
359 >
360 >          tempI = theBends[j].a;
361 >          tempJ = theBends[j].b;
362 >          tempK = theBends[j].c;
363 >          
364 > #ifdef IS_MPI
365 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
366 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
367 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
368 > #else
369 >          exI = tempI + 1;
370 >          exJ = tempJ + 1;
371 >          exK = tempK + 1;
372 > #endif
373 >          
374 >          info[k].excludes->addPair(exI, exK);
375 >          info[k].excludes->addPair(exI, exJ);
376 >          info[k].excludes->addPair(exJ, exK);
377 >        }
378        }
383 –      excludeOffset += molInfo.nBends;
379  
380        for (j = 0; j < molInfo.nTorsions; j++){
381          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 389 | Line 384 | void SimSetup::makeMolecules(void){
384          theTorsions[j].c = currentTorsion->getC() + atomOffset;
385          theTorsions[j].d = currentTorsion->getD() + atomOffset;
386  
387 <        exI = theTorsions[j].a;
388 <        exJ = theTorsions[j].d;
387 >        tempI = theTorsions[j].a;      
388 >        tempJ = theTorsions[j].b;
389 >        tempK = theTorsions[j].c;
390 >        tempL = theTorsions[j].d;
391  
395 –        // exclude_I must always be the smaller of the pair
396 –        if (exI > exJ){
397 –          tempEx = exI;
398 –          exI = exJ;
399 –          exJ = tempEx;
400 –        }
392   #ifdef IS_MPI
393 <        tempEx = exI;
394 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
395 <        tempEx = exJ;
396 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
393 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
394 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
395 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
396 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
397 > #else
398 >        exI = tempI + 1;
399 >        exJ = tempJ + 1;
400 >        exK = tempK + 1;
401 >        exL = tempL + 1;
402 > #endif
403  
404 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
405 < #else  // isn't MPI
406 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
407 < #endif  //is_mpi
404 >        info[k].excludes->addPair(exI, exJ);
405 >        info[k].excludes->addPair(exI, exK);
406 >        info[k].excludes->addPair(exI, exL);        
407 >        info[k].excludes->addPair(exJ, exK);
408 >        info[k].excludes->addPair(exJ, exL);
409 >        info[k].excludes->addPair(exK, exL);
410        }
412 –      excludeOffset += molInfo.nTorsions;
411  
412 +      for (j = 0; j < molInfo.nRigidBodies; j++){
413  
414 <      // send the arrays off to the forceField for init.
414 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
415 >        nMembers = currentRigidBody->getNMembers();
416 >
417 >        // Create the Rigid Body:
418 >
419 >        myRB = new RigidBody();
420 >
421 >        sprintf(rbName,"%s_RB_%s", molName, j);
422 >        myRB->setType(rbName);
423 >        
424 >        for (rb1 = 0; rb1 < nMembers; rb1++) {
425 >
426 >          // molI is atom numbering inside this molecule
427 >          molI = currentRigidBody->getMember(rb1);    
428 >
429 >          // tempI is atom numbering on local processor
430 >          tempI = molI + atomOffset;
431 >
432 >          // currentAtom is the AtomStamp (which we need for
433 >          // rigid body reference positions)
434 >          currentAtom = comp_stamps[stampID]->getAtom(molI);
435  
436 +          // When we add to the rigid body, add the atom itself and
437 +          // the stamp info:
438 +
439 +          myRB->addAtom(info[k].atoms[tempI], currentAtom);
440 +          
441 +          // Add this atom to the Skip List for the integrators
442 + #ifdef IS_MPI
443 +          slI = info[k].atoms[tempI]->getGlobalIndex();
444 + #else
445 +          slI = tempI;
446 + #endif
447 +          skipList.insert(slI);
448 +          
449 +        }
450 +        
451 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
452 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
453 +            
454 +            tempI = currentRigidBody->getMember(rb1);
455 +            tempJ = currentRigidBody->getMember(rb2);
456 +            
457 +            // Some explanation is required here.
458 +            // Fortran indexing starts at 1, while c indexing starts at 0
459 +            // Also, in parallel computations, the GlobalIndex is
460 +            // used for the exclude list:
461 +            
462 + #ifdef IS_MPI
463 +            exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
464 +            exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
465 + #else
466 +            exI = tempI + 1;
467 +            exJ = tempJ + 1;
468 + #endif
469 +            
470 +            info[k].excludes->addPair(exI, exJ);
471 +            
472 +          }
473 +        }
474 +
475 +        molInfo.myRigidBodies.push_back(myRB);
476 +        info[k].rigidBodies.push_back(myRB);
477 +      }
478 +      
479 +
480 +      // After this is all set up, scan through the atoms to
481 +      // see if they can be added to the integrableObjects:
482 +
483 +      for (j = 0; j < molInfo.nAtoms; j++){
484 +
485 + #ifdef IS_MPI
486 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
487 + #else
488 +        slJ = j+atomOffset;
489 + #endif
490 +
491 +        // if they aren't on the skip list, then they can be integrated
492 +
493 +        if (skipList.find(slJ) == skipList.end()) {
494 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
495 +          info[k].integrableObjects.push_back(mySD);
496 +          molInfo.myIntegrableObjects.push_back(mySD);
497 +        }
498 +      }
499 +
500 +      // all rigid bodies are integrated:
501 +
502 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
503 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
504 +        info[k].integrableObjects.push_back(mySD);      
505 +        molInfo.myIntegrableObjects.push_back(mySD);
506 +      }
507 +    
508 +      
509 +      // send the arrays off to the forceField for init.
510 +      
511        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
512        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
513        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
514        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
515                                   theTorsions);
516  
423 –
517        info[k].molecules[i].initialize(molInfo);
518  
519  
# Line 428 | Line 521 | void SimSetup::makeMolecules(void){
521        delete[] theBonds;
522        delete[] theBends;
523        delete[] theTorsions;
524 <    }
524 >    }    
525    }
526  
527   #ifdef IS_MPI
# Line 438 | Line 531 | void SimSetup::makeMolecules(void){
531  
532    // clean up the forcefield
533  
534 <  the_ff->calcRcut();
534 >  if (!globals->haveLJrcut()){
535 >
536 >    the_ff->calcRcut();
537 >
538 >  } else {
539 >    
540 >    the_ff->setRcut( globals->getLJrcut() );
541 >  }
542 >
543    the_ff->cleanMe();
544   }
545  
# Line 1220 | Line 1321 | void SimSetup::calcSysValues(void){
1321    tot_bonds = 0;
1322    tot_bends = 0;
1323    tot_torsions = 0;
1324 +  tot_rigid = 0;
1325    for (i = 0; i < n_components; i++){
1326      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1327      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1328      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1329      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1330 +    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1331    }
1332 <
1332 >  
1333    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1334    molMembershipArray = new int[tot_atoms];
1335  
# Line 1248 | Line 1351 | void SimSetup::mpiMolDivide(void){
1351    int i, j, k;
1352    int localMol, allMol;
1353    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1354 +  int local_rigid;
1355 +  vector<int> globalAtomIndex;
1356 +  vector<int> globalMolIndex;
1357  
1358    mpiSim = new mpiSimulation(info);
1359  
1360 <  globalIndex = mpiSim->divideLabor();
1360 >  mpiSim->divideLabor();
1361 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1362 >  globalMolIndex = mpiSim->getGlobalMolIndex();
1363  
1364    // set up the local variables
1365  
# Line 1264 | Line 1372 | void SimSetup::mpiMolDivide(void){
1372    local_bonds = 0;
1373    local_bends = 0;
1374    local_torsions = 0;
1375 <  globalAtomIndex = 0;
1375 >  local_rigid = 0;
1376 >  globalAtomCounter = 0;
1377  
1269 –
1378    for (i = 0; i < n_components; i++){
1379      for (j = 0; j < components_nmol[i]; j++){
1380        if (mol2proc[allMol] == worldRank){
# Line 1274 | Line 1382 | void SimSetup::mpiMolDivide(void){
1382          local_bonds += comp_stamps[i]->getNBonds();
1383          local_bends += comp_stamps[i]->getNBends();
1384          local_torsions += comp_stamps[i]->getNTorsions();
1385 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1386          localMol++;
1387        }      
1388        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1389 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1390 <        globalAtomIndex++;
1389 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1390 >        globalAtomCounter++;
1391        }
1392  
1393        allMol++;
# Line 1287 | Line 1396 | void SimSetup::mpiMolDivide(void){
1396    local_SRI = local_bonds + local_bends + local_torsions;
1397  
1398    info[0].n_atoms = mpiSim->getMyNlocal();  
1399 +  
1400  
1401    if (local_atoms != info[0].n_atoms){
1402      sprintf(painCave.errMsg,
# Line 1319 | Line 1429 | void SimSetup::makeSysArrays(void){
1429  
1430    Atom** the_atoms;
1431    Molecule* the_molecules;
1322 –  Exclude** the_excludes;
1432  
1324 –
1433    for (l = 0; l < nInfo; l++){
1434      // create the atom and short range interaction arrays
1435  
# Line 1347 | Line 1455 | void SimSetup::makeSysArrays(void){
1455   #else // is_mpi
1456  
1457      molIndex = 0;
1458 <    globalAtomIndex = 0;
1458 >    globalAtomCounter = 0;
1459      for (i = 0; i < n_components; i++){
1460        for (j = 0; j < components_nmol[i]; j++){
1461          the_molecules[molIndex].setStampID(i);
1462          the_molecules[molIndex].setMyIndex(molIndex);
1463          the_molecules[molIndex].setGlobalIndex(molIndex);
1464          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1465 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1466 <          globalAtomIndex++;
1465 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1466 >          globalAtomCounter++;
1467          }
1468          molIndex++;
1469        }
# Line 1364 | Line 1472 | void SimSetup::makeSysArrays(void){
1472  
1473   #endif // is_mpi
1474  
1475 <
1476 <    if (info[l].n_SRI){
1477 <      Exclude::createArray(info[l].n_SRI);
1370 <      the_excludes = new Exclude * [info[l].n_SRI];
1371 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1372 <        the_excludes[ex] = new Exclude(ex);
1373 <      }
1374 <      info[l].globalExcludes = new int;
1375 <      info[l].n_exclude = info[l].n_SRI;
1376 <    }
1377 <    else{
1378 <      Exclude::createArray(1);
1379 <      the_excludes = new Exclude * ;
1380 <      the_excludes[0] = new Exclude(0);
1381 <      the_excludes[0]->setPair(0, 0);
1382 <      info[l].globalExcludes = new int;
1383 <      info[l].globalExcludes[0] = 0;
1384 <      info[l].n_exclude = 0;
1385 <    }
1386 <
1475 >    info[l].globalExcludes = new int;
1476 >    info[l].globalExcludes[0] = 0;
1477 >    
1478      // set the arrays into the SimInfo object
1479  
1480      info[l].atoms = the_atoms;
1481      info[l].molecules = the_molecules;
1482      info[l].nGlobalExcludes = 0;
1392 –    info[l].excludes = the_excludes;
1483  
1484      the_ff->setSimInfo(info);
1485    }
# Line 1662 | Line 1752 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1752    }
1753  
1754    theInfo.addProperty(zconsForcePolicy);
1755 +
1756 +  //set zcons gap
1757 +  DoubleData* zconsGap = new DoubleData();
1758 +  zconsGap->setID(ZCONSGAP_ID);
1759  
1760 +  if (globals->haveZConsGap()){
1761 +    zconsGap->setData(globals->getZconsGap());
1762 +    theInfo.addProperty(zconsGap);  
1763 +  }
1764 +
1765 +  //set zcons fixtime
1766 +  DoubleData* zconsFixtime = new DoubleData();
1767 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1768 +
1769 +  if (globals->haveZConsFixTime()){
1770 +    zconsFixtime->setData(globals->getZconsFixtime());
1771 +    theInfo.addProperty(zconsFixtime);  
1772 +  }
1773 +
1774 +  //set zconsUsingSMD
1775 +  IntData* zconsUsingSMD = new IntData();
1776 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1777 +
1778 +  if (globals->haveZConsUsingSMD()){
1779 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1780 +    theInfo.addProperty(zconsUsingSMD);  
1781 +  }
1782 +
1783    //Determine the name of ouput file and add it into SimInfo's property list
1784    //Be careful, do not use inFileName, since it is a pointer which
1785    //point to a string at master node, and slave nodes do not contain that string
# Line 1692 | Line 1809 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1809      tempParaItem.zPos = zconStamp[i]->getZpos();
1810      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1811      tempParaItem.kRatio = zconStamp[i]->getKratio();
1812 <
1812 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1813 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1814      zconsParaData->addItem(tempParaItem);
1815    }
1816  
# Line 1712 | Line 1830 | void SimSetup::makeMinimizer(){
1830   }
1831  
1832   void SimSetup::makeMinimizer(){
1833 < /*  
1834 <  OOPSEMinimizerBase* myOOPSEMinimizerBase;
1717 <  ObjFunctor1 * objFunc;
1718 <  OutputFunctor* outputFunc;
1719 <  ConcreteNLModel1* nlp;
1833 >
1834 >  OOPSEMinimizer* myOOPSEMinimizer;
1835    MinimizerParameterSet* param;
1836 <  ConjugateMinimizerBase* minimizer;
1722 <  int dim;
1836 >  char minimizerName[100];
1837    
1838    for (int i = 0; i < nInfo; i++){
1839 <    //creat
1726 <    myOOPSEMinimizerBase = new OOPSEMinimizerBase(&(info[i]), the_ff);
1727 <
1728 <    //creat the object functor;
1729 <    objFunc = (ObjFunctor1*) new ClassMemObjFunctor1<OOPSEMinimizerBase>
1730 <                                              (myOOPSEMinimizerBase, &OOPSEMinimizerBase::calcGradient);
1731 <
1732 <    //creat output functor;
1733 <    outputFunc =  new ClassMemOutputFunctor<OOPSEMinimizerBase>
1734 <                               (myOOPSEMinimizerBase, &OOPSEMinimizerBase::output);
1735 <
1736 <    //creat nonlinear model
1737 <    dim = myOOPSEMinimizerBase->getDim();    
1738 <    nlp = new ConcreteNLModel1(dim, objFunc);
1739 <
1839 >    
1840      //prepare parameter set for minimizer
1841      param = new MinimizerParameterSet();
1842      param->setDefaultParameter();
# Line 1761 | Line 1861 | void SimSetup::makeMinimizer(){
1861        param->setWriteFrq(globals->getMinWriteFrq());
1862      }
1863      
1864 <    if (globals->haveMinResetFrq()){
1865 <      param->setResetFrq(globals->getMinResetFrq());
1864 >    if (globals->haveMinStepSize()){
1865 >      param->setStepSize(globals->getMinStepSize());
1866      }
1867  
1868      if (globals->haveMinLSMaxIter()){
# Line 1772 | Line 1872 | void SimSetup::makeMinimizer(){
1872      if (globals->haveMinLSTol()){
1873        param->setLineSearchTol(globals->getMinLSTol());
1874      }    
1775 –    
1776 –     //creat the minimizer
1777 –     minimizer = new PRCGMinimizer(nlp, param);
1778 –     minimizer->setLineSearchStrategy(nlp, GoldenSection);
1779 –     minimizer->setOutputFunctor(outputFunc);
1875  
1876 +    strcpy(minimizerName, globals->getMinimizer());
1877 +
1878 +    if (!strcasecmp(minimizerName, "CG")){
1879 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1880 +    }
1881 +    else if (!strcasecmp(minimizerName, "SD")){
1882 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1883 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1884 +    }
1885 +    else{
1886 +          sprintf(painCave.errMsg,
1887 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1888 +          painCave.isFatal = 0;
1889 +          simError();
1890 +
1891 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1892 +    }
1893 +     info[i].the_integrator = myOOPSEMinimizer;
1894 +
1895       //store the minimizer into simInfo
1896 <     info[i].the_minimizer = minimizer;
1896 >     info[i].the_minimizer = myOOPSEMinimizer;
1897       info[i].has_minimizer = true;
1898    }
1899 < */
1899 >
1900   }

Diff Legend

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