ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/InteractionManager.cpp
(Generate patch)

Comparing branches/development/src/nonbonded/InteractionManager.cpp (file contents):
Revision 1502 by gezelter, Sat Oct 2 19:53:32 2010 UTC vs.
Revision 1530 by gezelter, Tue Dec 28 21:47:55 2010 UTC

# Line 44 | Line 44 | namespace OpenMD {
44   namespace OpenMD {
45  
46    bool InteractionManager::initialized_ = false;
47 +  RealType InteractionManager::rCut_ = -1.0;
48 +  RealType InteractionManager::rSwitch_ = -1.0;
49    ForceField* InteractionManager::forceField_ = NULL;  
50    InteractionManager* InteractionManager::_instance = NULL;
51    map<int, AtomType*> InteractionManager::typeMap_;
52    map<pair<AtomType*, AtomType*>, set<NonBondedInteraction*> > InteractionManager::interactions_;
53 +
54 +  LJ* InteractionManager::lj_ = new LJ();
55 +  GB* InteractionManager::gb_ = new GB();
56 +  Sticky* InteractionManager::sticky_ = new Sticky();
57 +  Morse* InteractionManager::morse_ = new Morse();
58 +  EAM* InteractionManager::eam_ = new EAM();
59 +  SC* InteractionManager::sc_ = new SC();
60 +  Electrostatic* InteractionManager::electrostatic_ = new Electrostatic();
61 +  SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction();
62  
63    InteractionManager* InteractionManager::Instance() {
64      if (!_instance) {
# Line 58 | Line 69 | namespace OpenMD {
69  
70    void InteractionManager::initialize() {
71      
61 –    lj_ = new LJ();
62 –    gb_ = new GB();
63 –    sticky_ = new Sticky();
64 –    eam_ = new EAM();
65 –    sc_ = new SC();
66 –    morse_ = new Morse();
67 –    electrostatic_ = new Electrostatic();
68 –
72      lj_->setForceField(forceField_);
73      gb_->setForceField(forceField_);
74      sticky_->setForceField(forceField_);
# Line 74 | Line 77 | namespace OpenMD {
77      morse_->setForceField(forceField_);
78      electrostatic_->setForceField(forceField_);
79  
80 +    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
81 +    // Force fields can set options on how to scale van der Waals and electrostatic
82 +    // interactions for atoms connected via bonds, bends and torsions
83 +    // in this case the topological distance between atoms is:
84 +    // 0 = the atom itself
85 +    // 1 = bonded together
86 +    // 2 = connected via a bend
87 +    // 3 = connected via a torsion
88 +
89 +    vdwScale_[0] = 0.0;
90 +    vdwScale_[1] = fopts.getvdw12scale();
91 +    vdwScale_[2] = fopts.getvdw13scale();
92 +    vdwScale_[3] = fopts.getvdw14scale();
93 +
94 +    electrostaticScale_[0] = 0.0;
95 +    electrostaticScale_[1] = fopts.getelectrostatic12scale();
96 +    electrostaticScale_[2] = fopts.getelectrostatic13scale();
97 +    electrostaticScale_[3] = fopts.getelectrostatic14scale();
98 +
99      ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
100      ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
101      AtomType* atype1;
# Line 231 | Line 253 | namespace OpenMD {
253        }
254      }
255    }
234 –
235 –
236 –  void InteractionManager::doPrePair(AtomType* atype1,
237 –                                     AtomType* atype2,
238 –                                     RealType rij,
239 –                                     RealType &rho_i_at_j,
240 –                                     RealType &rho_j_at_i) {
241 –    
242 –  }
256    
257 <  void InteractionManager::doPreForce(AtomType* atype,
258 <                                      RealType rho,      
259 <                                      RealType &frho,
260 <                                      RealType &dfrhodrho) {
261 <  }
257 >  void InteractionManager::doPrePair(int *atid1, int *atid2, RealType *rij, RealType *rho_i_at_j, RealType *rho_j_at_i){
258 >    
259 >    if (!initialized_) initialize();
260 >          
261 >    DensityData ddat;
262  
263 <  void InteractionManager::doSkipCorrection(AtomType* atype1,      
264 <                                            AtomType* atype2,
265 <                                            Vector3d d,
266 <                                            RealType rij,
267 <                                            RealType &skippedCharge1,
255 <                                            RealType &skippedCharge2,
256 <                                            RealType sw,
257 <                                            RealType electroMult,
258 <                                            RealType &pot,
259 <                                            RealType &vpair,
260 <                                            Vector3d &f1,
261 <                                            Mat3x3d eFrame1,
262 <                                            Mat3x3d eFrame2,
263 <                                            Vector3d &t1,
264 <                                            Vector3d &t2) {
265 <  }
266 <  
267 <  void InteractionManager::doSelfCorrection(AtomType* atype,
268 <                                            Mat3x3d eFrame,
269 <                                            RealType skippedCharge,
270 <                                            RealType &pot,
271 <                                            Vector3d &t) {
272 <  }
263 >    ddat.atype1 = typeMap_[*atid1];
264 >    ddat.atype2 = typeMap_[*atid2];
265 >    ddat.rij = *rij;
266 >    ddat.rho_i_at_j = *rho_i_at_j;
267 >    ddat.rho_j_at_i = *rho_j_at_i;
268  
269 <  void InteractionManager::do_prepair(int *atid1, int *atid2, RealType *rij, RealType *rho_i_at_j, RealType *rho_j_at_i){
270 <    
276 <    if (!initialized_) initialize();
277 <    AtomType* atype1 = typeMap_[*atid1];
278 <    AtomType* atype2 = typeMap_[*atid2];
279 <    
280 <    doPrePair(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
281 <    
282 <    return;    
283 <  }
269 >    pair<AtomType*, AtomType*> key = make_pair(ddat.atype1, ddat.atype2);
270 >    set<NonBondedInteraction*>::iterator it;
271  
272 <  void InteractionManager::do_preforce(int *atid, RealType *rho, RealType *frho, RealType *dfrhodrho){
273 <
274 <    if (!initialized_) initialize();    
275 <    AtomType* atype = typeMap_[*atid];
276 <
290 <    doPreForce(atype, *rho, *frho, *dfrhodrho);
272 >    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
273 >      if ((*it)->getFamily() == METALLIC_FAMILY) {
274 >        dynamic_cast<MetallicInteraction*>(*it)->calcDensity(ddat);
275 >      }
276 >    }
277      
278      return;    
279    }
280 +  
281 +  void InteractionManager::doPreForce(int *atid, RealType *rho, RealType *frho, RealType *dfrhodrho){
282  
295 –  void InteractionManager::do_pair(int *atid1, int *atid2, RealType *d, RealType *r, RealType *r2, RealType *rcut, RealType *sw, RealType *vdwMult,RealType *electroMult, RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, RealType *eFrame2, RealType *A1, RealType *A2, RealType *t1, RealType *t2, RealType *rho1, RealType *rho2, RealType *dfrho1, RealType *dfrho2, RealType *fshift1, RealType *fshift2){
296 –
283      if (!initialized_) initialize();
284            
285 <    InteractionData idat;
285 >    FunctionalData fdat;
286  
287 +    fdat.atype = typeMap_[*atid];
288 +    fdat.rho = *rho;
289 +    fdat.frho = *frho;
290 +    fdat.dfrhodrho = *dfrhodrho;
291 +
292 +    pair<AtomType*, AtomType*> key = make_pair(fdat.atype, fdat.atype);
293 +    set<NonBondedInteraction*>::iterator it;
294 +    
295 +    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
296 +      if ((*it)->getFamily() == METALLIC_FAMILY) {
297 +        dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(fdat);
298 +      }
299 +    }
300 +    
301 +    return;    
302 +  }
303 +
304 +  void InteractionManager::doPair(int *atid1, int *atid2, RealType *d, RealType *r, RealType *r2, RealType *rcut, RealType *sw, int *topoDist, RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, RealType *eFrame2, RealType *A1, RealType *A2, RealType *t1, RealType *t2, RealType *rho1, RealType *rho2, RealType *dfrho1, RealType *dfrho2, RealType *fshift1, RealType *fshift2){
305 +    
306 +    if (!initialized_) initialize();
307 +    
308 +    InteractionData idat;
309 +    
310      idat.atype1 = typeMap_[*atid1];
311      idat.atype2 = typeMap_[*atid2];
312      idat.d = Vector3d(d);
# Line 305 | Line 314 | namespace OpenMD {
314      idat.r2 = *r2;
315      idat.rcut = *rcut;
316      idat.sw = *sw;
317 <    idat.vdwMult = *vdwMult;
318 <    idat.electroMult = *electroMult;
317 >    idat.vdwMult = vdwScale_[*topoDist];
318 >    idat.electroMult = electrostaticScale_[*topoDist];
319      idat.pot = *pot;
320      idat.vpair = *vpair;
321      idat.f1 = Vector3d(f1);
# Line 344 | Line 353 | namespace OpenMD {
353      return;    
354    }
355  
356 <  void InteractionManager::do_skip_correction(int *atid1, int *atid2, RealType *d, RealType *r, RealType *skippedCharge1, RealType *skippedCharge2, RealType *sw, RealType *electroMult, RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, RealType *eFrame2, RealType *t1, RealType *t2){
356 >  void InteractionManager::doSkipCorrection(int *atid1, int *atid2, RealType *d, RealType *r, RealType *skippedCharge1, RealType *skippedCharge2, RealType *sw, RealType *electroMult, RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, RealType *eFrame2, RealType *t1, RealType *t2){
357  
358      if (!initialized_) initialize();
350 –          
351 –    AtomType* atype1 = typeMap_[*atid1];
352 –    AtomType* atype2 = typeMap_[*atid2];
353 –    Vector3d disp(d);
354 –    Vector3d frc(f1);
355 –    Vector3d trq1(t1);
356 –    Vector3d trq2(t2);
357 –    Mat3x3d eFi(eFrame1);
358 –    Mat3x3d eFj(eFrame2);
359 –      
360 –    doSkipCorrection(atype1, atype2, disp, *r, *skippedCharge1, *skippedCharge2, *sw,
361 –                     *electroMult, *pot,  *vpair, frc, eFi, eFj, trq1, trq2);
362 –
363 –    f1[0] = frc.x();
364 –    f1[1] = frc.y();
365 –    f1[2] = frc.z();
359      
360 <    t1[0] = trq1.x();
368 <    t1[1] = trq1.y();
369 <    t1[2] = trq1.z();
360 >    SkipCorrectionData skdat;
361      
362 <    t2[0] = trq2.x();
363 <    t2[1] = trq2.y();
364 <    t2[2] = trq2.z();
362 >    skdat.atype1 = typeMap_[*atid1];
363 >    skdat.atype2 = typeMap_[*atid2];
364 >    skdat.d = Vector3d(d);
365 >    skdat.rij = *r;
366 >    skdat.skippedCharge1 = *skippedCharge1;
367 >    skdat.skippedCharge2 = *skippedCharge2;
368 >    skdat.sw = *sw;
369 >    skdat.electroMult = *electroMult;
370 >    skdat.pot = *pot;
371 >    skdat.vpair = *vpair;
372 >    skdat.f1 = Vector3d(f1);
373 >    skdat.eFrame1 = Mat3x3d(eFrame1);
374 >    skdat.eFrame2 = Mat3x3d(eFrame2);
375 >    skdat.t1 = Vector3d(t1);
376 >    skdat.t2 = Vector3d(t2);
377  
378 <    return;
379 <  }
378 >    pair<AtomType*, AtomType*> key = make_pair(skdat.atype1, skdat.atype2);
379 >    set<NonBondedInteraction*>::iterator it;
380  
381 <  void InteractionManager::do_self_correction(int *atid, RealType *eFrame, RealType *skippedCharge, RealType *pot, RealType *t){
381 >    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
382 >      if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
383 >        dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(skdat);
384 >      }
385 >    }
386 >    
387 >    f1[0] = skdat.f1.x();
388 >    f1[1] = skdat.f1.y();
389 >    f1[2] = skdat.f1.z();
390 >    
391 >    t1[0] = skdat.t1.x();
392 >    t1[1] = skdat.t1.y();
393 >    t1[2] = skdat.t1.z();
394 >    
395 >    t2[0] = skdat.t2.x();
396 >    t2[1] = skdat.t2.y();
397 >    t2[2] = skdat.t2.z();
398  
399 +    return;    
400 +  }
401 +
402 +  void InteractionManager::doSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, RealType *pot, RealType *t){
403 +
404      if (!initialized_) initialize();
405 <          
405 >    
406 >    SelfCorrectionData scdat;
407 >    
408 >    scdat.atype = typeMap_[*atid];
409 >    scdat.eFrame = Mat3x3d(eFrame);
410 >    scdat.skippedCharge = *skippedCharge;
411 >    scdat.pot = *pot;
412 >    scdat.t = Vector3d(t);
413 >
414 >    pair<AtomType*, AtomType*> key = make_pair(scdat.atype, scdat.atype);
415 >    set<NonBondedInteraction*>::iterator it;
416 >
417 >    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
418 >      if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
419 >        dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(scdat);
420 >      }
421 >    }
422 >        
423 >    t[0] = scdat.t.x();
424 >    t[1] = scdat.t.y();
425 >    t[2] = scdat.t.z();
426 >
427 >    return;    
428 >  }
429 >
430 >
431 >  RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
432 >    if (!initialized_) initialize();
433 >    
434      AtomType* atype = typeMap_[*atid];
383 –    Mat3x3d eFi(eFrame);
384 –    Vector3d trq1(t);
385 –      
386 –    doSelfCorrection(atype, eFi, *skippedCharge, *pot, trq1);
435  
436 <    t[0] = trq1.x();
437 <    t[1] = trq1.y();
438 <    t[2] = trq1.z();
436 >    pair<AtomType*, AtomType*> key = make_pair(atype, atype);
437 >    set<NonBondedInteraction*>::iterator it;
438 >    RealType cutoff = 0.0;
439 >    
440 >    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
441 >      cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype));  
442 >    return cutoff;    
443 >  }
444  
445 <    return;
445 >  RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) {
446 >    if (!initialized_) initialize();
447 >    
448 >    pair<AtomType*, AtomType*> key = make_pair(atype, atype);
449 >    set<NonBondedInteraction*>::iterator it;
450 >    RealType cutoff = 0.0;
451 >    
452 >    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
453 >      cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype));  
454 >    return cutoff;    
455    }
456  
457 +
458 +  void InteractionManager::setSwitch(RealType *rIn, RealType *rOut) {
459 +    switcher_->setSwitch(*rIn, *rOut);    
460 +  }
461 +
462 +  void InteractionManager::getSwitch(RealType *r2, RealType *sw, RealType *dswdr, RealType *r,
463 +                                     int *in_switching_region) {
464 +    bool isr = switcher_->getSwitch(*r2, *sw, *dswdr, *r);    
465 +    *in_switching_region = (int)isr;
466 +  }
467 +
468   } //end namespace OpenMD
469  
470   extern "C" {
# Line 402 | Line 475 | extern "C" {
475   #define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION)
476   #define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION)
477   #define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF)
478 + #define fortranSetSwitch FC_FUNC(set_switch, SET_SWITCH)
479 + #define fortranGetSwitch FC_FUNC(get_switch, GET_SWITCH)
480  
481    void fortranDoPrePair(int *atid1, int *atid2, RealType *rij,
482                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
483              
484 <    return OpenMD::InteractionManager::Instance()->do_prepair(atid1, atid2, rij,
485 <                                                              rho_i_at_j,  
486 <                                                              rho_j_at_i);
484 >    return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij,
485 >                                                             rho_i_at_j,  
486 >                                                             rho_j_at_i);
487    }
488 <  void fortranDoPreforce(int *atid, RealType *rho, RealType *frho,
488 >  void fortranDoPreForce(int *atid, RealType *rho, RealType *frho,
489                           RealType *dfrhodrho) {  
490      
491 <    return OpenMD::InteractionManager::Instance()->do_preforce(atid, rho, frho,
492 <                                                               dfrhodrho);    
491 >    return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho,
492 >                                                              dfrhodrho);    
493    }
494    
495    void fortranDoPair(int *atid1, int *atid2, RealType *d, RealType *r,
496 <                     RealType *r2, RealType *rcut, RealType *sw, RealType *vdwMult,
497 <                     RealType *electroMult, RealType *pot, RealType *vpair, RealType *f1,
498 <                     RealType *eFrame1, RealType *eFrame2, RealType *A1, RealType *A2,
499 <                     RealType *t1, RealType *t2,
500 <                     RealType *rho1, RealType *rho2, RealType *dfrho1, RealType *dfrho2,
501 <                     RealType *fshift1, RealType *fshift2){
496 >                     RealType *r2, RealType *rcut, RealType *sw, int *topoDist,
497 >                     RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1,
498 >                     RealType *eFrame2, RealType *A1, RealType *A2,
499 >                     RealType *t1, RealType *t2, RealType *rho1, RealType *rho2,
500 >                     RealType *dfrho1, RealType *dfrho2, RealType *fshift1,
501 >                     RealType *fshift2){
502      
503 <    return OpenMD::InteractionManager::Instance()->do_pair(atid1, atid2, d, r, r2, rcut,
504 <                                                           sw, vdwMult, electroMult, pot,
505 <                                                           vpair, f1, eFrame1, eFrame2,
506 <                                                           A1, A2, t1, t2, rho1, rho2,
507 <                                                           dfrho1, dfrho2, fshift1, fshift2);
503 >    return OpenMD::InteractionManager::Instance()->doPair(atid1, atid2, d, r,
504 >                                                          r2, rcut, sw, topoDist,
505 >                                                          pot, vpair, f1,
506 >                                                          eFrame1, eFrame2,
507 >                                                          A1, A2, t1, t2, rho1,
508 >                                                          rho2, dfrho1, dfrho2,
509 >                                                          fshift1, fshift2);
510    }
511 <
512 <  void fortranDoSkipCorrection(int *atid1, int *atid2, RealType *d, RealType *r,
513 <                               RealType *skippedCharge1, RealType *skippedCharge2,
514 <                               RealType *sw, RealType *electroMult, RealType *pot,
511 >  
512 >  void fortranDoSkipCorrection(int *atid1, int *atid2, RealType *d,
513 >                               RealType *r, RealType *skippedCharge1,
514 >                               RealType *skippedCharge2, RealType *sw,
515 >                               RealType *electroMult, RealType *pot,
516                                 RealType *vpair, RealType *f1,
517                                 RealType *eFrame1, RealType *eFrame2,
518                                 RealType *t1, RealType *t2){
519      
520 <    return OpenMD::InteractionManager::Instance()->do_skip_correction(atid1, atid2, d, r,
521 <                                                                      skippedCharge1,
522 <                                                                      skippedCharge2,
523 <                                                                      sw, electroMult, pot,
524 <                                                                      vpair, f1, eFrame1,
525 <                                                                      eFrame2, t1, t2);
520 >    return OpenMD::InteractionManager::Instance()->doSkipCorrection(atid1,
521 >                                                                    atid2, d,
522 >                                                                    r,
523 >                                                                    skippedCharge1,
524 >                                                                    skippedCharge2,
525 >                                                                    sw, electroMult, pot,
526 >                                                                    vpair, f1, eFrame1,
527 >                                                                    eFrame2, t1, t2);
528    }
529 <
529 >  
530    void fortranDoSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge,
531                                 RealType *pot, RealType *t) {
532      
533 <    return OpenMD::InteractionManager::Instance()->do_self_correction(atid, eFrame,
534 <                                                                      skippedCharge,
535 <                                                                      pot, t);
533 >    return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid,
534 >                                                                    eFrame,
535 >                                                                    skippedCharge,
536 >                                                                    pot, t);
537    }
538 <  RealType fortranGetCutoff() {
539 <    return OpenMD::InteractionManager::Instance()->getCutoff();
538 >  RealType fortranGetCutoff(int *atid) {    
539 >    return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid);
540    }
541 +
542 +  void fortranGetSwitch(RealType *r2, RealType *sw, RealType *dswdr, RealType *r,
543 +                        int *in_switching_region) {
544 +    
545 +    return OpenMD::InteractionManager::Instance()->getSwitch(r2, sw, dswdr, r,
546 +                                                             in_switching_region);
547 +  }
548 +
549 +  void fortranSetSwitch(RealType *rIn, RealType *rOut) {    
550 +    return OpenMD::InteractionManager::Instance()->setSwitch(rIn, rOut);
551 +  }
552 +  
553   }

Diff Legend

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