| 40 |
|
*/ |
| 41 |
|
|
| 42 |
|
#include "nonbonded/InteractionManager.hpp" |
| 43 |
– |
#include "UseTheForce/doForces_interface.h" |
| 43 |
|
|
| 44 |
|
namespace OpenMD { |
| 45 |
|
|
| 53 |
|
RealType InteractionManager::listRadius_ = 0.0; |
| 54 |
|
CutoffMethod InteractionManager::cutoffMethod_ = SHIFTED_FORCE; |
| 55 |
|
SwitchingFunctionType InteractionManager::sft_ = cubic; |
| 56 |
< |
RealType InteractionManager::vdwScale_[4] = {0.0, 0.0, 0.0, 0.0}; |
| 57 |
< |
RealType InteractionManager::electrostaticScale_[4] = {0.0, 0.0, 0.0, 0.0}; |
| 56 |
> |
RealType InteractionManager::vdwScale_[4] = {1.0, 0.0, 0.0, 0.0}; |
| 57 |
> |
RealType InteractionManager::electrostaticScale_[4] = {1.0, 0.0, 0.0, 0.0}; |
| 58 |
|
|
| 59 |
|
map<int, AtomType*> InteractionManager::typeMap_; |
| 60 |
|
map<pair<AtomType*, AtomType*>, set<NonBondedInteraction*> > InteractionManager::interactions_; |
| 94 |
|
// Force fields can set options on how to scale van der Waals and electrostatic |
| 95 |
|
// interactions for atoms connected via bonds, bends and torsions |
| 96 |
|
// in this case the topological distance between atoms is: |
| 97 |
< |
// 0 = the atom itself |
| 97 |
> |
// 0 = topologically unconnected |
| 98 |
|
// 1 = bonded together |
| 99 |
|
// 2 = connected via a bend |
| 100 |
|
// 3 = connected via a torsion |
| 101 |
|
|
| 102 |
< |
vdwScale_[0] = 0.0; |
| 102 |
> |
vdwScale_[0] = 1.0; |
| 103 |
|
vdwScale_[1] = fopts.getvdw12scale(); |
| 104 |
|
vdwScale_[2] = fopts.getvdw13scale(); |
| 105 |
|
vdwScale_[3] = fopts.getvdw14scale(); |
| 106 |
|
|
| 107 |
< |
electrostaticScale_[0] = 0.0; |
| 107 |
> |
electrostaticScale_[0] = 1.0; |
| 108 |
|
electrostaticScale_[1] = fopts.getelectrostatic12scale(); |
| 109 |
|
electrostaticScale_[2] = fopts.getelectrostatic13scale(); |
| 110 |
|
electrostaticScale_[3] = fopts.getelectrostatic14scale(); |
| 306 |
|
setupCutoffs(); |
| 307 |
|
setupSwitching(); |
| 308 |
|
setupNeighborlists(); |
| 310 |
– |
notifyFortranSkinThickness(&skinThickness_); |
| 309 |
|
|
| 310 |
< |
int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
| 311 |
< |
int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
| 312 |
< |
notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf); |
| 313 |
< |
|
| 316 |
< |
int isError; |
| 317 |
< |
initFortranFF(&isError); |
| 310 |
> |
//int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
| 311 |
> |
//int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
| 312 |
> |
//notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf); |
| 313 |
> |
//notifyFortranSkinThickness(&skinThickness_); |
| 314 |
|
|
| 315 |
|
initialized_ = true; |
| 316 |
|
} |
| 453 |
|
} |
| 454 |
|
} |
| 455 |
|
} |
| 456 |
+ |
|
| 457 |
+ |
switcher_->setSwitchType(sft_); |
| 458 |
+ |
switcher_->setSwitch(rSwitch_, rCut_); |
| 459 |
|
} |
| 460 |
|
|
| 461 |
|
/** |
| 482 |
|
} |
| 483 |
|
|
| 484 |
|
listRadius_ = rCut_ + skinThickness_; |
| 486 |
– |
|
| 485 |
|
} |
| 486 |
|
|
| 487 |
|
|
| 488 |
< |
void InteractionManager::doPrePair(int *atid1, int *atid2, RealType *rij, RealType *rho_i_at_j, RealType *rho_j_at_i){ |
| 488 |
> |
void InteractionManager::doPrePair(InteractionData idat){ |
| 489 |
|
|
| 490 |
|
if (!initialized_) initialize(); |
| 491 |
< |
|
| 494 |
< |
DensityData ddat; |
| 495 |
< |
|
| 496 |
< |
ddat.atype1 = typeMap_[*atid1]; |
| 497 |
< |
ddat.atype2 = typeMap_[*atid2]; |
| 498 |
< |
ddat.rij = *rij; |
| 499 |
< |
ddat.rho_i_at_j = *rho_i_at_j; |
| 500 |
< |
ddat.rho_j_at_i = *rho_j_at_i; |
| 501 |
< |
|
| 502 |
< |
pair<AtomType*, AtomType*> key = make_pair(ddat.atype1, ddat.atype2); |
| 491 |
> |
|
| 492 |
|
set<NonBondedInteraction*>::iterator it; |
| 493 |
|
|
| 494 |
< |
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ |
| 494 |
> |
for (it = interactions_[ *(idat.atypes) ].begin(); |
| 495 |
> |
it != interactions_[ *(idat.atypes) ].end(); ++it){ |
| 496 |
|
if ((*it)->getFamily() == METALLIC_FAMILY) { |
| 497 |
< |
dynamic_cast<MetallicInteraction*>(*it)->calcDensity(ddat); |
| 497 |
> |
dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat); |
| 498 |
|
} |
| 499 |
|
} |
| 500 |
|
|
| 501 |
|
return; |
| 502 |
|
} |
| 503 |
|
|
| 504 |
< |
void InteractionManager::doPreForce(int *atid, RealType *rho, RealType *frho, RealType *dfrhodrho){ |
| 504 |
> |
void InteractionManager::doPreForce(SelfData sdat){ |
| 505 |
|
|
| 506 |
|
if (!initialized_) initialize(); |
| 507 |
< |
|
| 508 |
< |
FunctionalData fdat; |
| 519 |
< |
|
| 520 |
< |
fdat.atype = typeMap_[*atid]; |
| 521 |
< |
fdat.rho = *rho; |
| 522 |
< |
fdat.frho = *frho; |
| 523 |
< |
fdat.dfrhodrho = *dfrhodrho; |
| 524 |
< |
|
| 525 |
< |
pair<AtomType*, AtomType*> key = make_pair(fdat.atype, fdat.atype); |
| 507 |
> |
|
| 508 |
> |
pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype); |
| 509 |
|
set<NonBondedInteraction*>::iterator it; |
| 510 |
|
|
| 511 |
|
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ |
| 512 |
|
if ((*it)->getFamily() == METALLIC_FAMILY) { |
| 513 |
< |
dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(fdat); |
| 513 |
> |
dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(sdat); |
| 514 |
|
} |
| 515 |
|
} |
| 516 |
|
|
| 517 |
|
return; |
| 518 |
|
} |
| 519 |
|
|
| 520 |
< |
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){ |
| 520 |
> |
void InteractionManager::doPair(InteractionData idat){ |
| 521 |
|
|
| 522 |
|
if (!initialized_) initialize(); |
| 523 |
< |
|
| 541 |
< |
InteractionData idat; |
| 542 |
< |
|
| 543 |
< |
idat.atype1 = typeMap_[*atid1]; |
| 544 |
< |
idat.atype2 = typeMap_[*atid2]; |
| 545 |
< |
idat.d = Vector3d(d); |
| 546 |
< |
idat.rij = *r; |
| 547 |
< |
idat.r2 = *r2; |
| 548 |
< |
idat.rcut = *rcut; |
| 549 |
< |
idat.sw = *sw; |
| 550 |
< |
idat.vdwMult = vdwScale_[*topoDist]; |
| 551 |
< |
idat.electroMult = electrostaticScale_[*topoDist]; |
| 552 |
< |
idat.pot = *pot; |
| 553 |
< |
idat.vpair = *vpair; |
| 554 |
< |
idat.f1 = Vector3d(f1); |
| 555 |
< |
idat.eFrame1 = Mat3x3d(eFrame1); |
| 556 |
< |
idat.eFrame2 = Mat3x3d(eFrame2); |
| 557 |
< |
idat.A1 = RotMat3x3d(A1); |
| 558 |
< |
idat.A2 = RotMat3x3d(A2); |
| 559 |
< |
idat.t1 = Vector3d(t1); |
| 560 |
< |
idat.t2 = Vector3d(t2); |
| 561 |
< |
idat.rho1 = *rho1; |
| 562 |
< |
idat.rho2 = *rho2; |
| 563 |
< |
idat.dfrho1 = *dfrho1; |
| 564 |
< |
idat.dfrho2 = *dfrho2; |
| 565 |
< |
idat.fshift1 = *fshift1; |
| 566 |
< |
idat.fshift2 = *fshift2; |
| 567 |
< |
|
| 568 |
< |
pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2); |
| 523 |
> |
|
| 524 |
|
set<NonBondedInteraction*>::iterator it; |
| 525 |
|
|
| 526 |
< |
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) |
| 526 |
> |
for (it = interactions_[ *(idat.atypes) ].begin(); |
| 527 |
> |
it != interactions_[ *(idat.atypes) ].end(); ++it) |
| 528 |
|
(*it)->calcForce(idat); |
| 529 |
|
|
| 574 |
– |
f1[0] = idat.f1.x(); |
| 575 |
– |
f1[1] = idat.f1.y(); |
| 576 |
– |
f1[2] = idat.f1.z(); |
| 577 |
– |
|
| 578 |
– |
t1[0] = idat.t1.x(); |
| 579 |
– |
t1[1] = idat.t1.y(); |
| 580 |
– |
t1[2] = idat.t1.z(); |
| 581 |
– |
|
| 582 |
– |
t2[0] = idat.t2.x(); |
| 583 |
– |
t2[1] = idat.t2.y(); |
| 584 |
– |
t2[2] = idat.t2.z(); |
| 585 |
– |
|
| 530 |
|
return; |
| 531 |
|
} |
| 532 |
|
|
| 533 |
< |
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){ |
| 533 |
> |
void InteractionManager::doSkipCorrection(InteractionData idat){ |
| 534 |
|
|
| 535 |
< |
if (!initialized_) initialize(); |
| 535 |
> |
if (!initialized_) initialize(); |
| 536 |
|
|
| 593 |
– |
SkipCorrectionData skdat; |
| 594 |
– |
|
| 595 |
– |
skdat.atype1 = typeMap_[*atid1]; |
| 596 |
– |
skdat.atype2 = typeMap_[*atid2]; |
| 597 |
– |
skdat.d = Vector3d(d); |
| 598 |
– |
skdat.rij = *r; |
| 599 |
– |
skdat.skippedCharge1 = *skippedCharge1; |
| 600 |
– |
skdat.skippedCharge2 = *skippedCharge2; |
| 601 |
– |
skdat.sw = *sw; |
| 602 |
– |
skdat.electroMult = *electroMult; |
| 603 |
– |
skdat.pot = *pot; |
| 604 |
– |
skdat.vpair = *vpair; |
| 605 |
– |
skdat.f1 = Vector3d(f1); |
| 606 |
– |
skdat.eFrame1 = Mat3x3d(eFrame1); |
| 607 |
– |
skdat.eFrame2 = Mat3x3d(eFrame2); |
| 608 |
– |
skdat.t1 = Vector3d(t1); |
| 609 |
– |
skdat.t2 = Vector3d(t2); |
| 610 |
– |
|
| 611 |
– |
pair<AtomType*, AtomType*> key = make_pair(skdat.atype1, skdat.atype2); |
| 537 |
|
set<NonBondedInteraction*>::iterator it; |
| 538 |
|
|
| 539 |
< |
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ |
| 539 |
> |
for (it = interactions_[ *(idat.atypes) ].begin(); |
| 540 |
> |
it != interactions_[ *(idat.atypes) ].end(); ++it){ |
| 541 |
|
if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
| 542 |
< |
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(skdat); |
| 542 |
> |
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat); |
| 543 |
|
} |
| 544 |
|
} |
| 545 |
|
|
| 620 |
– |
f1[0] = skdat.f1.x(); |
| 621 |
– |
f1[1] = skdat.f1.y(); |
| 622 |
– |
f1[2] = skdat.f1.z(); |
| 623 |
– |
|
| 624 |
– |
t1[0] = skdat.t1.x(); |
| 625 |
– |
t1[1] = skdat.t1.y(); |
| 626 |
– |
t1[2] = skdat.t1.z(); |
| 627 |
– |
|
| 628 |
– |
t2[0] = skdat.t2.x(); |
| 629 |
– |
t2[1] = skdat.t2.y(); |
| 630 |
– |
t2[2] = skdat.t2.z(); |
| 631 |
– |
|
| 546 |
|
return; |
| 547 |
|
} |
| 548 |
|
|
| 549 |
< |
void InteractionManager::doSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, RealType *pot, RealType *t){ |
| 549 |
> |
void InteractionManager::doSelfCorrection(SelfData sdat){ |
| 550 |
|
|
| 551 |
|
if (!initialized_) initialize(); |
| 552 |
|
|
| 553 |
< |
SelfCorrectionData scdat; |
| 640 |
< |
|
| 641 |
< |
scdat.atype = typeMap_[*atid]; |
| 642 |
< |
scdat.eFrame = Mat3x3d(eFrame); |
| 643 |
< |
scdat.skippedCharge = *skippedCharge; |
| 644 |
< |
scdat.pot = *pot; |
| 645 |
< |
scdat.t = Vector3d(t); |
| 646 |
< |
|
| 647 |
< |
pair<AtomType*, AtomType*> key = make_pair(scdat.atype, scdat.atype); |
| 553 |
> |
pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype); |
| 554 |
|
set<NonBondedInteraction*>::iterator it; |
| 555 |
|
|
| 556 |
|
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ |
| 557 |
|
if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
| 558 |
< |
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(scdat); |
| 558 |
> |
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(sdat); |
| 559 |
|
} |
| 560 |
|
} |
| 561 |
< |
|
| 656 |
< |
t[0] = scdat.t.x(); |
| 657 |
< |
t[1] = scdat.t.y(); |
| 658 |
< |
t[2] = scdat.t.z(); |
| 659 |
< |
|
| 561 |
> |
|
| 562 |
|
return; |
| 563 |
|
} |
| 564 |
|
|
| 663 |
– |
|
| 565 |
|
RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { |
| 566 |
|
if (!initialized_) initialize(); |
| 567 |
|
|
| 572 |
|
RealType cutoff = 0.0; |
| 573 |
|
|
| 574 |
|
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) |
| 575 |
< |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype)); |
| 575 |
> |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); |
| 576 |
|
return cutoff; |
| 577 |
|
} |
| 578 |
|
|
| 584 |
|
RealType cutoff = 0.0; |
| 585 |
|
|
| 586 |
|
for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) |
| 587 |
< |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype)); |
| 587 |
> |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); |
| 588 |
|
return cutoff; |
| 589 |
|
} |
| 689 |
– |
|
| 690 |
– |
|
| 691 |
– |
void InteractionManager::setSwitch(RealType *rIn, RealType *rOut) { |
| 692 |
– |
switcher_->setSwitch(*rIn, *rOut); |
| 693 |
– |
} |
| 694 |
– |
|
| 695 |
– |
void InteractionManager::getSwitch(RealType *r2, RealType *sw, RealType *dswdr, RealType *r, |
| 696 |
– |
int *in_switching_region) { |
| 697 |
– |
bool isr = switcher_->getSwitch(*r2, *sw, *dswdr, *r); |
| 698 |
– |
*in_switching_region = (int)isr; |
| 699 |
– |
} |
| 700 |
– |
|
| 590 |
|
} //end namespace OpenMD |
| 702 |
– |
|
| 703 |
– |
extern "C" { |
| 704 |
– |
|
| 705 |
– |
#define fortranDoPrePair FC_FUNC(do_prepair, DO_PREPAIR) |
| 706 |
– |
#define fortranDoPreForce FC_FUNC(do_preforce, DO_PREFORCE) |
| 707 |
– |
#define fortranDoPair FC_FUNC(do_pair, DO_PAIR) |
| 708 |
– |
#define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION) |
| 709 |
– |
#define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION) |
| 710 |
– |
#define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF) |
| 711 |
– |
#define fortranSetSwitch FC_FUNC(set_switch, SET_SWITCH) |
| 712 |
– |
#define fortranGetSwitch FC_FUNC(get_switch, GET_SWITCH) |
| 713 |
– |
|
| 714 |
– |
void fortranDoPrePair(int *atid1, int *atid2, RealType *rij, |
| 715 |
– |
RealType *rho_i_at_j, RealType *rho_j_at_i) { |
| 716 |
– |
|
| 717 |
– |
return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij, |
| 718 |
– |
rho_i_at_j, |
| 719 |
– |
rho_j_at_i); |
| 720 |
– |
} |
| 721 |
– |
void fortranDoPreForce(int *atid, RealType *rho, RealType *frho, |
| 722 |
– |
RealType *dfrhodrho) { |
| 723 |
– |
|
| 724 |
– |
return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho, |
| 725 |
– |
dfrhodrho); |
| 726 |
– |
} |
| 727 |
– |
|
| 728 |
– |
void fortranDoPair(int *atid1, int *atid2, RealType *d, RealType *r, |
| 729 |
– |
RealType *r2, RealType *rcut, RealType *sw, int *topoDist, |
| 730 |
– |
RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, |
| 731 |
– |
RealType *eFrame2, RealType *A1, RealType *A2, |
| 732 |
– |
RealType *t1, RealType *t2, RealType *rho1, RealType *rho2, |
| 733 |
– |
RealType *dfrho1, RealType *dfrho2, RealType *fshift1, |
| 734 |
– |
RealType *fshift2){ |
| 735 |
– |
|
| 736 |
– |
return OpenMD::InteractionManager::Instance()->doPair(atid1, atid2, d, r, |
| 737 |
– |
r2, rcut, sw, topoDist, |
| 738 |
– |
pot, vpair, f1, |
| 739 |
– |
eFrame1, eFrame2, |
| 740 |
– |
A1, A2, t1, t2, rho1, |
| 741 |
– |
rho2, dfrho1, dfrho2, |
| 742 |
– |
fshift1, fshift2); |
| 743 |
– |
} |
| 744 |
– |
|
| 745 |
– |
void fortranDoSkipCorrection(int *atid1, int *atid2, RealType *d, |
| 746 |
– |
RealType *r, RealType *skippedCharge1, |
| 747 |
– |
RealType *skippedCharge2, RealType *sw, |
| 748 |
– |
RealType *electroMult, RealType *pot, |
| 749 |
– |
RealType *vpair, RealType *f1, |
| 750 |
– |
RealType *eFrame1, RealType *eFrame2, |
| 751 |
– |
RealType *t1, RealType *t2){ |
| 752 |
– |
|
| 753 |
– |
return OpenMD::InteractionManager::Instance()->doSkipCorrection(atid1, |
| 754 |
– |
atid2, d, |
| 755 |
– |
r, |
| 756 |
– |
skippedCharge1, |
| 757 |
– |
skippedCharge2, |
| 758 |
– |
sw, electroMult, pot, |
| 759 |
– |
vpair, f1, eFrame1, |
| 760 |
– |
eFrame2, t1, t2); |
| 761 |
– |
} |
| 762 |
– |
|
| 763 |
– |
void fortranDoSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, |
| 764 |
– |
RealType *pot, RealType *t) { |
| 765 |
– |
|
| 766 |
– |
return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid, |
| 767 |
– |
eFrame, |
| 768 |
– |
skippedCharge, |
| 769 |
– |
pot, t); |
| 770 |
– |
} |
| 771 |
– |
RealType fortranGetCutoff(int *atid) { |
| 772 |
– |
return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid); |
| 773 |
– |
} |
| 774 |
– |
|
| 775 |
– |
void fortranGetSwitch(RealType *r2, RealType *sw, RealType *dswdr, RealType *r, |
| 776 |
– |
int *in_switching_region) { |
| 777 |
– |
|
| 778 |
– |
return OpenMD::InteractionManager::Instance()->getSwitch(r2, sw, dswdr, r, |
| 779 |
– |
in_switching_region); |
| 780 |
– |
} |
| 781 |
– |
|
| 782 |
– |
void fortranSetSwitch(RealType *rIn, RealType *rOut) { |
| 783 |
– |
return OpenMD::InteractionManager::Instance()->setSwitch(rIn, rOut); |
| 784 |
– |
} |
| 785 |
– |
|
| 786 |
– |
} |
